so4bnd.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>
module so4bnd
!----------------------------------------------------------------------- 
! 
! Purpose: 
!
! SO4 boundary module.  Deals with interpolating SO4 datasets.
! 
! Author: Brian Eaton
! 
!-----------------------------------------------------------------------
   use precision
!
!JR Stuck this "only" business in because Compaq compiler barfed on pcnst, pnats having dual
!JR declarations when radctl.F90 gets compiled.
!
   use pmgrid,    only: plon, plat, masterproc
   use ppgrid,    only: pcols, pver, begchunk, endchunk
   use phys_grid, only: scatter_field_to_chunk, get_ncols_p
   implicit none
   save
!
! Floating point data
!
   real(r8), private, allocatable, dimension(:,:,:,:) :: &
      sulfbioi    ! input sulfate bio mixing ratios (pcols,pver,begchunk:endchunk,2)
   real(r8), private, allocatable, dimension(:,:,:) :: &
      sulfbio     ! time interpolated sulfate bio mixing ratios (pcols,pver,begchunk:endchunk)
   real(r8), private, allocatable, dimension(:,:,:,:) :: &
      sulfanti    ! input sulfate ant mixing ratios (pcols,pver,begchunk:endchunk,2)
   real(r8), private, allocatable, dimension(:,:,:) :: &
      sulfant     ! time interpolated sulfate ant mixing ratios (pcols,pver,begchunk:endchunk)

   real(r8), private :: sulfscalef                  ! Sulfate scale factor (for 1870->1990 ramp) 
   real(r8), private :: cdaysulfm           ! calendar day for prv. month sulfate values read in
   real(r8), private :: cdaysulfp           ! calendar day for nxt. month sulfate values read in

   integer, private :: date_sulf(1000)              ! Date on sulfate dataset (YYYYMMDD)
   integer, private :: sec_sulf(1000)               ! seconds of date on sulfate dataset (0-86399)
!
! just check that hard-wired size is big enough
!
!
! Integer data
!
   integer, private :: nm,np      ! Array indices for prv., nxt month sulfate data
   integer, private :: np1        ! current forward time index of sulfate dataset
   integer, private :: ncid_sulf  ! sulfate dataset id
   integer, private :: sulfbio_id ! netcdf id for sulfate mmr bio variable
   integer, private :: sulfant_id ! netcdf id for sulfate mmr anth variable
   integer, private :: lonsiz     ! size of longitude dimension on sulfate dataset
   integer, private :: levsiz     ! size of level dimension on sulfate dataset
   integer, private :: latsiz     ! size of latitude dimension on sulfate dataset
   integer, private :: timsiz     ! size of time dimension on sulfate dataset
 
!
! Logical variables
!
   logical, private :: sulfcyc    ! If sulfur cycle code turned on or not

   character*80, private ::  sulfdata ! full pathname for sulfate dataset

contains

subroutine so4bndnl( xsulfdata )
!----------------------------------------------------------------------- 
! 
! Purpose: Set variables from namelist input.
! 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------

   character*80, intent(in):: xsulfdata ! full pathname for sulfate dataset

!-----------------------------------------------------------------------
   sulfdata = xsulfdata
   if (masterproc) &
      write(6,*)'Time-variant sulfate dataset is: ',trim(sulfdata)

   return
end subroutine so4bndnl

!###############################################################################

subroutine sulfini
!----------------------------------------------------------------------- 
! 
! Purpose: Do initial read of time-variant sulfate dataset, containing
!          sulfate mixing ratios as a function of time.  It is currently
!          required that the sulfate dataset have the *SAME* horizontal
!          and vertical resolution as the model. Therefore, ONLY a time
!          interpolation of the dataset is currently performed.
! 
!-----------------------------------------------------------------------
   use ioFileMod
   use error_messages, only: alloc_err, handle_ncerr
   use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
                           is_perpetual
#if ( defined SPMD )
   use mpishorthand
#endif
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comlun.h>
!-----------------------------------------------------------------------
   include 'netcdf.inc'
!
! Local workspace
!
   character(len=256) locfn      ! local filename        
!
   integer dateid                ! netcdf id for date variable
   integer secid                 ! netcdf id for seconds variable
   integer londimid              ! netcdf id for longitude dimension
   integer latdimid              ! netcdf id for latitude dimension
   integer levdimid              ! netcdf id for level dimension
   integer lonid                 ! netcdf id for longitude variable
   integer latid                 ! netcdf id for latitude variable
   integer levid                 ! netcdf id for level variable
   integer timeid                ! netcdf id for time variable
   integer cnt4(4)               ! array of counts for each dimension
   integer strt4(4)              ! array of starting indices
   integer i, k, lat, n          ! longitude, level, latitude, time indices
   integer istat                 ! error return
   integer dimids(nf_max_var_dims) ! netcdf variable shape
   integer  :: yr, mon, day, ncsec ! components of a date
   integer  :: ncdate              ! current date in integer format [yyyymmdd]
   real(r8) :: calday              ! current calendar day
   real(r8) caldayloc                ! calendar day (includes yr if no cycling)
   real(r8) xsulfbioi(plon,pver,plat,2)  ! input sulfate bio mixing ratios
   real(r8) xsulfanti(plon,pver,plat,2)  ! input sulfate ant mixing ratios
!
!-----------------------------------------------------------------------
!
! Initialize time counters
!
   nm = 1
   np = 2
!
! Allocate space for data.
!
  allocate( sulfbioi(pcols,pver,begchunk:endchunk,2), stat=istat )
  call alloc_err( istat, 'sulfini', 'sulfbioi', &
       pcols*pver*(endchunk-begchunk+1)*2 )
  allocate( sulfbio(pcols,pver,begchunk:endchunk), stat=istat )
  call alloc_err( istat, 'sulfini', 'sulfbio', &
       pcols*pver*(endchunk-begchunk+1) )
  allocate( sulfanti(pcols,pver,begchunk:endchunk,2), stat=istat )
  call alloc_err( istat, 'sulfini', 'sulfanti', &
       pcols*pver*(endchunk-begchunk+1)*2 )
  allocate( sulfant(pcols,pver,begchunk:endchunk), stat=istat )
  call alloc_err( istat, 'sulfini', 'sulfant', &
       pcols*pver*(endchunk-begchunk+1) )
!
! SPMD: Master does all the work.  Sends needed info to slaves
!
   if (masterproc) then
!
! Obtain dataset
!
      call getfil(sulfdata, locfn)
      call wrap_open(locfn, 0, ncid_sulf)
!
! Currently assume that cycle over 12 months of data
!
      sulfcyc = .true.
!
! Use year information only if not cycling sulfate 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 (sulfcyc) then
         caldayloc = calday
      else
         caldayloc = calday + yr*365.
      end if
!
! Obtain dimension id's
!
      call wrap_inq_dimid( ncid_sulf, 'lon', londimid)
      call wrap_inq_dimid( ncid_sulf, 'lat', latdimid)
      call wrap_inq_dimid( ncid_sulf, 'lev', levdimid)
      call wrap_inq_dimid( ncid_sulf, 'time',timeid  )
!
! Obtain size of dimensions.
! Check that horizontal and vertical dimensions are same as model's
!
      call wrap_inq_dimlen( ncid_sulf, londimid, lonsiz   )
      if (lonsiz /= plon) then
         write(6,*)'SULFINI: lonsiz=',lonsiz,' must = ',plon
         call endrun
      end if

      call wrap_inq_dimlen( ncid_sulf, latdimid, latsiz   )
      if (latsiz /= plat) then
         write(6,*)'SULFINI: latsiz=',latsiz,' must = ',plat
         call endrun
      end if

      call wrap_inq_dimlen( ncid_sulf, levdimid, levsiz   )
      if (levsiz /= pver) then
         write(6,*)'SULFINI: levsiz=',levsiz,' must = ',pver
         call endrun
      end if

      call wrap_inq_dimlen( ncid_sulf, timeid, timsiz   )
!
! Obtain date info id's
!
      call wrap_inq_varid( ncid_sulf, 'date'   , dateid  )
      call wrap_inq_varid( ncid_sulf, 'datesec', secid   )
! 
! Obtain sulfate mixing ratio id
!
      call wrap_inq_varid( ncid_sulf, 'sulfmmrbio' , sulfbio_id )
      call wrap_inq_varid( ncid_sulf, 'sulfmmranth', sulfant_id )

      call wrap_inq_vardimid (ncid_sulf, sulfbio_id, dimids)
      if (dimids(1) /= londimid .and. dimids(2) /= levdimid .and. dimids(3) /= latdimid) then
         write(6,*)'SULFINI: Data must be ordered lon, lev, lat, time'
         call endrun
      end if
!
! just check that hard-wired size is big enough
!
      if (timsiz > 1000) then
         write(6,*)'SO4BND: timsiz=',timsiz,' too small'
         call endrun
      end if
!
! Determine date ids
!
      call wrap_get_var_int (ncid_sulf, dateid, date_sulf)
      call wrap_get_var_int (ncid_sulf, secid, sec_sulf)
!
! If cycling data first do error checks
!
      if (sulfcyc) then
         if (timsiz.lt.12) then 
            write(6,*)'SULFINI: When cycling sulfate dataset must have 12 consecutive ', &
                      'months of data starting with Jan'
            write(6,*)'Current dataset has only ',timsiz,' months'
            call endrun
         end if
         do n = 1,12
            if (mod(date_sulf(n),10000)/100/=n) then
               write(6,*)'SULFINI: When cycling sulfate dataset must have 12 consecutive ', &
                         'months of data starting with Jan'
               write(6,*)'Month ',n,' of dataset says date= ', date_sulf(n)
               call endrun
            end if
         end do
      end if
!
! Set up hyperslab corners
!
      strt4(1) = 1
      strt4(2) = 1
      strt4(3) = 1
      cnt4(1)  = lonsiz
      cnt4(2)  = levsiz
      cnt4(3)  = latsiz
      cnt4(4)  = 1
!
! Special code for interpolation between December and January
!
      if (sulfcyc) then
         n = 12
         np1 = 1
         call bnddyi(date_sulf(n  ), sec_sulf(n  ), cdaysulfm)
         call bnddyi(date_sulf(np1), sec_sulf(np1), cdaysulfp)
         if (caldayloc.le.cdaysulfp .or. caldayloc.gt.cdaysulfm) then
            strt4(4) = n
            call wrap_get_vara_realx (ncid_sulf,sulfbio_id,strt4,cnt4,xsulfbioi(1,1,1,nm))
            call wrap_get_vara_realx (ncid_sulf,sulfant_id,strt4,cnt4,xsulfanti(1,1,1,nm))
            strt4(4) = np1
            call wrap_get_vara_realx (ncid_sulf,sulfbio_id,strt4,cnt4,xsulfbioi(1,1,1,np))
            call wrap_get_vara_realx (ncid_sulf,sulfant_id,strt4,cnt4,xsulfanti(1,1,1,np))
            goto 10
         end if
      end if