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  returnend subroutine sstintend module sst_data