sst_data.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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
  return
end 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 endrun
10   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)

  return
end 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, &