#include <misc.h>#include <params.h>!----------------------------------------------------------------------- !! BOP!! !MODULE: sst_data!! !DESCRIPTION:	Module to handle dealing with the Sea-Surface Temperature !			datasets.  This module also figures out the location of!			sea-ice from these datasets where it is assumed that ! 	seawater at freezing or below is a flag for the existence of sea-ice.!	SST datasets that are created for use with the stand-alone CCM should!	take this into account and set grid-points where sea-ice fraction is!	greater than 50% to -1.8C and ensure that other grid points where sea-ice!	is less than 50% have SST's greater than -1.8C.!! Public interfaces:!!	sstini -- Initialization and reading of dataset.!	sstint -- Interpolate dataset SST to current time.!	sstan --- Apply the interpolated SST to the model state.!! $Id: sst_data.F90,v 1.9.2.3 2002/05/02 21:11:29 rosinski Exp $!!----------------------------------------------------------------------- module sst_data!! USES:!  use precision, only: r8  use pmgrid,    only: plon, plat, masterproc  use ppgrid,    only: pcols, begchunk, endchunk  use phys_grid, only: scatter_field_to_chunk, get_ncols_p  use comsrf,    only: plevmx, icefrac  use physconst, only: tmelt  use commap,    only: clat, clon  implicit none  public :: sst ! needed in ice_data!----------------------------------------------------------------------- ! PUBLIC: Make default data and interfaces private!----------------------------------------------------------------------- !! ! PUBLIC MEMBER FUNCTIONS:!  public sstan    ! Set the surface temperature, oro, and sea-ice fraction  public sstini   ! Initialization  public sstint   ! Time interpolation of SST data!===============================================================================!EOP!===============================================================================!----------------------------------------------------------------------- ! PRIVATE: Everthing else is private to this module!-----------------------------------------------------------------------   private   ! By default all data is private to this module  integer, parameter :: totsstsz=2000  real(r8), parameter :: daysperyear = 365.0  ! Number of days in a year  real(r8), allocatable, dimension(:,:,:) :: &      sstbdy         ! SST values on boundary dataset (pcols,begchunk:endchunk,2)  real(r8), allocatable, dimension(:,:) :: &      sst            ! Interpolated model sst values (pcols,begchunk:endchunk)  real(r8) :: cdaysstm   ! Calendar day for prv. month SST values read in  real(r8) :: cdaysstp   ! Calendar day for nxt. month SST values read in        integer :: nm,np   ! Array indices for prv., nxt month sst data  integer :: sstid   ! netcdf id for sst variable  integer :: lonsiz  ! size of longitude dimension on sst dataset  integer :: levsiz  ! size of level dimension on sst dataset  integer :: latsiz  ! size of latitude dimension on sst dataset  integer :: timesiz ! size of time dimension on sst dataset  integer :: np1     ! current forward time index of sst dataset  integer :: date_sst(totsstsz)! Date on sst dataset (YYYYMMDD)  integer :: sec_sst(totsstsz) ! seconds of date on sst dataset (0-86399)  real(r8), parameter :: tsice = -1.7999 ! Freezing point of sea ice degrees C                                         ! Use this with global sst data !===============================================================================CONTAINS!===============================================================================!======================================================================! PUBLIC ROUTINES: Following routines are publically accessable!======================================================================!----------------------------------------------------------------------- ! ! BOP!! !IROUTINE: sstan!! !DESCRIPTION: ! Update sea surface temperatures (sst's) and sea ice distribution!! Method: ! Assume that the sst data exists in a two dimensional field! encoded as follows:!     Land               values where oro field says so ("valid sst's!                        are provided globally, the model's land mask!                        determines whether the sst is used or not)!     Ocean without      values degrees celcius (greater than tsice)!      sea ice!     Ocean with         values less than tsice!      sea ice! New sea ice has a constant 0.5 cm value for snow cover prescribed! ! Author: CCM1!!-----------------------------------------------------------------------!! !INTERFACE:!subroutine sstan(lchnk   ,ncol    ,ocnfrac ,ts      )!! !INPUT PARAMETERS:!  integer , intent(in)    :: lchnk                ! chunk identifier  integer , intent(in)    :: ncol                 ! number of atmospheric columns!! !INPUT/OUTPUT PARAMETERS:!  real(r8), intent(in) :: ocnfrac(pcols)           ! Surface type flag array  real(r8), intent(inout) :: ts(pcols)            ! Surface temperature!-----------------------------------------------------------------------! EOP!---------------------------Local variables-----------------------------  integer i                 ! Column index!-----------------------------------------------------------------------! Open ocean!  do i=1,ncol     if (ocnfrac(i) > 0.) then        ts(i) = sst(i,lchnk) + tmelt     end if  end do  returnend subroutine sstan!----------------------------------------------------------------------- ! ! BOP!! !IROUTINE: sstini!! !DESCRIPTION:!! Initialize the procedure for specifying sea surface temperatures! Do initial read of time-varying sst boundary dataset, reading two! consecutive months on either side of the current model date.!! Method: ! ! Author: L.Bath! !-----------------------------------------------------------------------!! !INTERFACE!subroutine sstini!! !USES:!  use rgrid, only: nlon  use error_messages, only: alloc_err, handle_ncerr  use time_manager, only: get_curr_date, get_curr_calday, &                          is_perpetual, get_perp_date#if ( defined SPMD )  use mpishorthand, only: mpicom, mpiint, mpir8#endif!! EOP!!---------------------------Common blocks-------------------------------#include <comctl.h>#include <comlun.h>!---------------------------Local variables-----------------------------  integer dateid                ! netcdf id for date variable  integer secid                 ! netcdf id for seconds variable  integer londimid              ! netcdf id for longitude variable  integer latdimid              ! netcdf id for latitude variable  integer lonid                 ! netcdf id for longitude variable  integer latid                 ! netcdf id for latitude variable  integer timeid                ! netcdf id for time variable  integer nlonid                ! netcdf id for nlon variable (rgrid)  integer cnt3(3)               ! array of counts for each dimension  integer strt3(3)              ! array of starting indices  integer n                     ! indices  integer nlon_sst(plat)        ! number of lons per lat on bdy dataset  integer j                     ! latitude index  integer istat                 ! error return  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               ! calendar day (includes yr if no cycling)  real(r8) caldayloc            ! calendar day (includes yr if no cycling)  real(r8) xvar(plon,plat,2)    ! work space !-----------------------------------------------------------------------!! Initialize time indices!  nm = 1  np = 2!! Allocate space for data.!  allocate( sst(pcols,begchunk:endchunk), stat=istat )  call alloc_err( istat, 'sstini', 'sst', &       pcols*(endchunk-begchunk+1) )  if(aqua_planet) return  allocate( sstbdy(pcols,begchunk:endchunk,2), stat=istat )  call alloc_err( istat, 'sstini', 'sstbdy', &       pcols*(endchunk-begchunk+1)*2 )!! SPMD: Master does all the work.!  if (masterproc) then!! Use year information only if not cycling sst 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!! Get and check dimension info!     call wrap_inq_dimid( ncid_sst, 'lon', londimid   )     call wrap_inq_dimid( ncid_sst, 'time', timeid  )     call wrap_inq_dimid( ncid_sst, 'lat', latdimid   )     call wrap_inq_dimlen( ncid_sst, londimid, lonsiz   )     if (lonsiz /= plon) then        write(6,*)'SSTINI: lonsiz=',lonsiz,' must = plon=',plon        call endrun     end if     call wrap_inq_dimlen( ncid_sst, latdimid, latsiz   )     if (latsiz /= plat) then        write(6,*)'SSTINI: latsiz=',latsiz,' must = plat=',plat        call endrun     end if     call wrap_inq_dimlen( ncid_sst, timeid, timesiz   )!! Check to make sure space allocated for time variables is sufficient!     if (timesiz>totsstsz) then        write(6,*)'SSTINI:  Allocated space for sst data is insufficient.'        write(6,*)'Please increase parameter totsstsz to',timesiz,' and recompile.'        call endrun     end if!! Check to ensure reduced or not grid of dataset matches that of model!     if (fullgrid) then        call wrap_inq_varid( ncid_sst, 'lon', lonid   )     else        call wrap_inq_varid (ncid_sst, 'nlon', nlonid)        call wrap_get_var_int (ncid_sst, nlonid, nlon_sst)        do j=1,plat           if (nlon_sst(j) /= nlon(j)) then              write(6,*)'SSTINI: model grid does not match dataset grid'              call endrun           end if        end do     end if     call wrap_inq_varid( ncid_sst, 'date', dateid   )     call wrap_inq_varid( ncid_sst, 'datesec', secid   )     call wrap_inq_varid( ncid_sst, 'SST_cpl', sstid   )     call wrap_inq_varid( ncid_sst, 'lat', latid   )!! Retrieve entire date and sec variables.!     call wrap_get_var_int (ncid_sst,dateid,date_sst)     call wrap_get_var_int (ncid_sst,secid,sec_sst)     if (sstcyc) then        if (timesiz<12) then            write(6,*)'SSTINI: ERROR'            write(6,*)'When cycling sst, sst data set must have 12'            write(6,*)'consecutive months of data starting with Jan'           write(6,*)'Current dataset has only ',timesiz,' months'           call endrun        end if        do n = 1,12           if (mod(date_sst(n),10000)/100/=n) then              write(6,*)'SSTINI: ERROR'               write(6,*)'When cycling sst, sst data set must have 12'               write(6,*)'consecutive months of data starting with Jan'              write(6,*)'Month ',n,' of sst data set is out of order'              call endrun           end if        end do     end if     strt3(1) = 1     strt3(2) = 1     strt3(3) = 1     cnt3(1)  = lonsiz     cnt3(2)  = latsiz     cnt3(3)  = 1!! Special code for interpolation between December and January!     if (sstcyc) then        n = 12        np1 = 1        call bnddyi(date_sst(n  ), sec_sst(n  ), cdaysstm)        call bnddyi(date_sst(np1), sec_sst(np1), cdaysstp)        if (caldayloc<=cdaysstp .or. caldayloc>cdaysstm) then           strt3(3) = n           call wrap_get_vara_realx (ncid_sst,sstid,strt3,cnt3,xvar(1,1,nm))           strt3(3) = np1                                                 call wrap_get_vara_realx (ncid_sst,sstid,strt3,cnt3,xvar(1,1,np))           goto 10        end if     end if!! Normal interpolation between consecutive time slices.!     do n=1,timesiz-1        np1 = n + 1        call bnddyi(date_sst(n  ), sec_sst(n  ), cdaysstm)        call bnddyi(date_sst(np1), sec_sst(np1), cdaysstp)        if (.not.sstcyc) then           yr = date_sst(n)/10000           cdaysstm = cdaysstm + yr*daysperyear           yr = date_sst(np1)/10000           cdaysstp = cdaysstp + yr*daysperyear        end if        if (caldayloc>cdaysstm .and. caldayloc<=cdaysstp) then           strt3(3) = n           call wrap_get_vara_realx (ncid_sst,sstid,strt3,cnt3,xvar(1,1,nm))           strt3(3) = np1                                                 call wrap_get_vara_realx (ncid_sst,sstid,strt3,cnt3,xvar(1,1,np))           goto 10        end if     end do     write(6,*)'SSTINI: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec     call endrun10   continue     write(6,*)'SSTINI: Read sst data for dates ',date_sst(n),sec_sst(n), &          ' and ',date_sst(np1),sec_sst(np1)#if (defined SPMD )     call mpibcast( timesiz, 1, mpiint, 0, mpicom )     call mpibcast( date_sst, totsstsz, mpiint, 0, mpicom )     call mpibcast( sec_sst, totsstsz, mpiint, 0, mpicom )     call mpibcast( cdaysstm, 1, mpir8, 0, mpicom )     call mpibcast( cdaysstp, 1, mpir8, 0, mpicom )     call mpibcast( np1, 1, mpiint, 0, mpicom )  else     call mpibcast( timesiz, 1, mpiint, 0, mpicom )     call mpibcast( date_sst, totsstsz, mpiint, 0, mpicom )     call mpibcast( sec_sst, totsstsz, mpiint, 0, mpicom )     call mpibcast( cdaysstm, 1, mpir8, 0, mpicom )     call mpibcast( cdaysstp, 1, mpir8, 0, mpicom )     call mpibcast( np1, 1, mpiint, 0, mpicom )#endif  end if  call scatter_field_to_chunk(1,1,2,plon,xvar,sstbdy)  returnend subroutine sstini!----------------------------------------------------------------------- ! ! BOP!! !IROUTINE: sstint!! !DESCRIPTION:!! if "aqua_planet", specify SST's analytically (Jerry Olson).! Otherwise, time interpolate SST's to current time, reading in new monthly data if! necessary.!! Method: ! ! Author: L.Bath! !-----------------------------------------------------------------------!! !INTERFACE:!subroutine sstint!! !USES:!  use rgrid, only: nlon  use comsrf, only: ocnfrac  use time_manager, only: get_curr_date, get_curr_calday, &