sst_data.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                          is_perpetual, get_perp_date
!
! EOP
!
!---------------------------Common blocks-------------------------------
#include <comctl.h>
#include <comlun.h>
!---------------------------Local variables-----------------------------
  integer cnt3(3)        ! array of counts for each dimension
  integer strt3(3)       ! array of starting indices
  integer i,j,lchnk      ! indices
  integer ncol           ! number of columns in current chunk
  integer ntmp           ! temporary
  real(r8) fact1, fact2  ! time interpolation factors
  integer :: yr, mon, day! components of a date
  integer :: ncdate      ! current date in integer format [yyyymmdd]
  integer :: ncsec       ! current time of day [seconds]
  real(r8) :: calday     ! current calendar day
  real(r8) caldayloc     ! calendar day (includes yr if no cycling)
  real(r8) deltat        ! time (days) between interpolating sst data
!
! Aqua planet variables
!
  real(r8) pi            ! 3.14159...
  real(r8) pio180        ! pi/180.
  real(r8) tmp           ! temporary
  real(r8) tmp1          ! temporary
  real(r8) t0_max        ! max reference temperature
  real(r8) t0_min        ! min reference temperature
  real(r8) t0_max6       ! max asymmetric reference temperature for option 6
  real(r8) t0_max7       ! max asymmetric reference temperature for option 7
  real(r8) maxlat        ! cutoff latitude poleward of which SST = 0 deg C
  real(r8) shift         ! number of degrees peak SST is shifted off equator
  real(r8) shift9        ! number of degrees peak SST is shifted off equator for opt. 9
  real(r8) shift10       ! number of degrees peak SST is shifted off equator for opt. 10
  real(r8) latcen        ! center of asymmetric SST forcing
  real(r8) latrad6       ! radius of asymmetric SST forcing for option 6
  real(r8) latrad8       ! radius of asymmetric SST forcing for option 8
  real(r8) loncen        ! center of asymmetric SST forcing
  real(r8) lonrad        ! radius of asymmetric SST forcing
  real(r8) xvar(plon,plat,2)    ! work space 
  integer  sst_option    ! option of analytical SST algorithm
!
!-----------------------------------------------------------------------
!
! SPMD: Master does all the work.  Sends needed info to slaves
!
  if (aqua_planet) then

     if (masterproc) then

        sst_option = 1
        pi         = 4.*atan(1.)
        pio180     = pi/180.
        if(sst_option .lt. 1 .or. sst_option .gt. 10) then
          write(6,*) 'ERROR SSTINT:  sst_option must be between 1 and 10'
          call endrun
        endif
!
! Parameters for zonally symmetric experiments
!
        t0_max     = 27.
        t0_min     = 0.
        maxlat     = 60.
        shift      = 5.
        shift9     = 10.
        shift10    = 15.
!
! Parameters for zonally asymmetric experiments
!
        t0_max6    = 1.
        t0_max7    = 3.
        latcen     = 0.
        loncen     = 90.
        latrad6    = 15.
        latrad8    = 30.
        lonrad     = 30.

        maxlat     = maxlat *pio180
        shift      = shift  *pio180
        shift9     = shift9 *pio180
        shift10    = shift10*pio180
        latcen     = latcen *pio180
        loncen     = loncen *pio180
        latrad6    = latrad6*pio180
        latrad8    = latrad8*pio180
        lonrad     = lonrad *pio180

        if(sst_option == 1 .or. sst_option == 6 .or. &
             sst_option == 7 .or. sst_option == 8     ) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              else
                 tmp = clat(j)*pi*0.5/maxlat
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        end if

        if(sst_option == 2) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              else
                 tmp = clat(j)*pi*0.5/maxlat
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp*tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        end if

        if(sst_option == 3) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              else
                 tmp = clat(j)*pi*0.5/maxlat
                 tmp = sin(tmp)
                 tmp = (2. - tmp*tmp*tmp*tmp - tmp*tmp)*0.5
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        end if

        if(sst_option == 4) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              else
                 tmp  = (maxlat - abs(clat(j)))/maxlat
                 tmp1 = 1. - tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_max*tmp + t0_min*tmp1
                 end do
              end if
           end do
        end if

        if(sst_option == 5) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              elseif(clat(j) > shift) then
                 tmp = (clat(j)-shift)*pi*0.5/(maxlat-shift)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              else
                 tmp = (clat(j)-shift)*pi*0.5/(maxlat+shift)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        end if

        if(sst_option == 6) then
           do j = 1,plat
              if(abs(clat(j)-latcen) <= latrad6) then
                 tmp1 = (clat(j)-latcen)*pi*0.5/latrad6
                 tmp1 = cos(tmp1)
                 tmp1 = tmp1*tmp1
                 do i=1,nlon(j)
                    if(abs(clon(i,j)-loncen) <= lonrad) then
                       tmp = (clon(i,j)-loncen)*pi*0.5/lonrad
                       tmp = cos(tmp)
                       tmp = tmp*tmp
                       xvar(i,j,1) = xvar(i,j,1) + t0_max6*tmp*tmp1
                    end if
                 end do
              end if
           end do
        end if

        if(sst_option == 7) then
           do j = 1,plat
              if(abs(clat(j)-latcen) <= latrad6) then
                 tmp1 = (clat(j)-latcen)*pi*0.5/latrad6
                 tmp1 = cos(tmp1)
                 tmp1 = tmp1*tmp1
                 do i=1,nlon(j)
                    if(abs(clon(i,j)-loncen) <= lonrad) then
                       tmp = (clon(i,j)-loncen)*pi*0.5/lonrad
                       tmp = cos(tmp)
                       tmp = tmp*tmp
                       xvar(i,j,1) = xvar(i,j,1) + t0_max7*tmp*tmp1
                    end if
                 end do
              end if
           end do
        end if

        if(sst_option == 8) then
           do j = 1,plat
              if(abs(clat(j)-latcen) <= latrad8) then
                 tmp1 = (clat(j)-latcen)*pi*0.5/latrad8
                 tmp1 = cos(tmp1)
                 tmp1 = tmp1*tmp1
                 do i=1,nlon(j)
                    tmp = cos(clon(i,j)-loncen)
                    xvar(i,j,1) = xvar(i,j,1) + t0_max7*tmp*tmp1
                 end do
              end if
           end do
        end if

        if(sst_option == 9) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              elseif(clat(j) > shift9) then
                 tmp = (clat(j)-shift9)*pi*0.5/(maxlat-shift9)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              else
                 tmp = (clat(j)-shift9)*pi*0.5/(maxlat+shift9)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        end if

        if(sst_option == 10) then
           do j = 1,plat
              if(abs(clat(j)) > maxlat) then
                 do i=1,nlon(j)
                    xvar(i,j,1) = t0_min
                 end do
              elseif(clat(j) > shift10) then
                 tmp = (clat(j)-shift10)*pi*0.5/(maxlat-shift10)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                 do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              else
                 tmp = (clat(j)-shift10)*pi*0.5/(maxlat+shift10)
                 tmp = sin(tmp)
                 tmp = 1. - tmp*tmp
                  do i=1,nlon(j)
                    xvar(i,j,1) = tmp*(t0_max - t0_min) + t0_min
                 end do
              end if
           end do
        endif

     endif

     call scatter_field_to_chunk(1,1,1,plon,xvar(1,1,1),sst(1,begchunk))

  else

!
! Use year information only if a multiyear dataset
!
     calday = get_curr_calday()
     if ( is_perpetual() ) then
        call get_perp_date(yr, mon, day, ncsec)
     else
        call get_curr_date(yr, mon, day, ncsec)
     end if
     ncdate = yr*10000 + mon*100 + day
     if (sstcyc) then
        caldayloc = calday
     else
        caldayloc = calday + yr*daysperyear
     end if

     if (masterproc) then

        strt3(1) = 1
        strt3(2) = 1
        strt3(3) = 1
        cnt3(1)  = lonsiz
        cnt3(2)  = latsiz
        cnt3(3)  = 1

     endif
!
! If model time is past current forward sst timeslice, read in the next
! timeslice for time interpolation.  Messy logic is for sstcyc = .true. 
! interpolation between December and January (np1==1).  Note that 
! np1 is never 1 when sstcyc is .false.
!
     if (caldayloc > cdaysstp .and. .not. (np1==1 .and. caldayloc>cdaysstm)) then
        if (sstcyc) then
           np1 = mod(np1,12) + 1
        else
           np1 = np1 + 1
        end if
        if (np1>timesiz) then
           if (masterproc) then
              write(6,*)'SSTINT: Attempt to read past end of SST dataset'
           endif
           call endrun
        end if
        cdaysstm = cdaysstp
        call bnddyi(date_sst(np1), sec_sst(np1), cdaysstp)

        if (.not.sstcyc) then
           yr = date_sst(np1)/10000
           cdaysstp = cdaysstp + yr*daysperyear
        end if

        if (np1==1 .or. caldayloc<=cdaysstp) then
           ntmp = nm
           nm = np
           np = ntmp
           if (masterproc) then
              strt3(3) = np1
              call wrap_get_vara_realx (ncid_sst,sstid,strt3,cnt3,xvar(1,1,np))
              write(6,*)'SSTINT: Read sst for date (yyyymmdd) ',date_sst(np1), &
                   ' sec ',sec_sst(np1)
           endif
           call scatter_field_to_chunk(1,1,1,plon,xvar(1,1,np),sstbdy(1,begchunk,np))
        else
           if (masterproc) then
              write(6,*)'SSTINT: Input sst for date',date_sst(np1), &
                   ' sec ',sec_sst(np1), 'does not exceed model date',ncdate,&
                   ' sec ',ncsec,' Stopping.'
           endif
           call endrun
        end if
     end if
!
! Time interpolation.  Account for December-January interpolation if
! cycling sst dataset.  Again note that np1 is never 1 when sstcyc is false
!
     if (np1==1) then                    ! Dec-Jan interpolation
        deltat = cdaysstp + daysperyear - cdaysstm
        if (caldayloc>cdaysstp) then      ! We're in December
           fact1 = (cdaysstp + daysperyear - caldayloc)/deltat
           fact2 = (caldayloc - cdaysstm)/deltat
        else                                ! We're in January
           fact1 = (cdaysstp - caldayloc)/deltat
           fact2 = (caldayloc + daysperyear - cdaysstm)/deltat
        end if
     else
        deltat = cdaysstp - cdaysstm
        fact1 = (cdaysstp - caldayloc)/deltat
        fact2 = (caldayloc - cdaysstm)/deltat
     end if
!
! Check sanity of time interpolation calculation to within 32-bit roundoff
!
     if (abs(fact1+fact2-1.)>1.e-6 .or. &
             fact1>1.000001 .or. fact1<-1.e-6 .or. &
             fact2>1.000001 .or. fact2<-1.e-6) then 
        if (masterproc) then
           write(6,*)'SSTINT: Bad fact1 and/or fact2=',fact1,fact2
        endif
        call endrun
     end if

     do lchnk=begchunk,endchunk
        ncol = get_ncols_p(lchnk)
        do i=1,ncol
           if (ocnfrac(i,lchnk)>=0.) then
              sst(i,lchnk) = sstbdy(i,lchnk,nm)*fact1 + sstbdy(i,lchnk,np)*fact2
!  Bound the sst temp by the freezing point of sea water.
              sst(i,lchnk)=max(sst(i,lchnk),tsice)
           end if
        end do
     end do
  endif

  return
end subroutine sstint


end module sst_data