atmdrvmod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

       end if

! Get variable names 

       status = nf_inq_nvars(ncid, nvar)
       if (status /= nf_noerr) then
          write (6,*) ' ATM_OPENFILE netcdf error = ',nf_strerror(status)
          call endrun
       end if
       do i = 1, nvar
          call wrap_inq_varname (ncid, i, varnam(i))
       end do

    endif     !end of if-masterproc block

! Calculate temporal resolution of the dataset in minutes and
! reset time index

#if (defined SPMD) 
    call mpi_bcast (ntim, 1, mpiint, 0, mpicom, ier)
#endif
    atmmin = ndaypm(kmo) * minpday / ntim
    itim = 0                 

    return
  end subroutine atm_openfile

!=======================================================================

  subroutine atm_readdata (fname, kmo, itim)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! read atmospheric forcing from netCDF file
!
! Method: 
! Read a netCDF data file.
!
! Every file that references a netCDF function must have the
! include statement: include 'netcdf.inc'
!
! netCDF datasets -
!    o referenced by a dataset id that is obtained when the
!      dataset is first created or opened
!    o nf_create (fname, nf_clobber, ncid): creates a new netCDF data file
!      with name [fname], overwritting if already exists, and returning
!      a dataset id [ncid] that is used to refer to dataset in other
!      netCDF function calls
!    o nf_open (fname, nf_nowrite, ncid): opens existing netCDF file
!      [fname] in read mode only, returning id [ncid]
!    o nf_enddef (ncid): takes open netCDF dataset, referenced by [ncid],
!      out of define mode
!    o nf_close (ncid): closes an open netCDF dataset, referenced by [ncid]
!
! netCDF dimensions -
!    o have both a name and a length
!    o one dimension in the dataset can have length unlimited (e.g., time
!      dimension can be unlimited to have multiple time slices in dataset)
!
! netCDF variables -
!    o A netCDF variable has a name, type [nf_char, nf_int, nf_float,
!      nf_double], shape (dimension), and attributes (e.g., long name, units).
!      These are specified when the variable is defined. A variable also
!      has values, which are specified in data mode.
!    o A netCDF variable is referenced by an integer variable (1,2,3,...),
!      which is in the order in which the variables are defined
!    o Character string is treated as an array of characters
!    o A coordinate variable is a special netCDF variable that has the
!      same name as a dimension (e.g., lat(lat), lon(lon)). It is used
!      by some appication packages to define the physical coordinate
!      corresponding to that dimension
!    o Variables are dimensioned in Fortran opposite of how listed
!      in netCDF file, e.g.:
!               Fortran             netCDF
!              --------------    --------------
!              x(lon,lat)     -> x(lat,lon)
!              x(lon,lat,lev) -> x(lev,lat,lon)
!              x(lon,lat,tim) -> x(tim,lat,lon)
! 
! Author: Sam Levis
! 
!-----------------------------------------------------------------------

    use precision
    use clm_varcon  , only : tfrz, sb
    use fileutils   , only : getfil
#if (defined SPMD)
    use mpishorthand, only : mpiint, mpir8, mpicom
#endif
    implicit none
    include 'netcdf.inc'

! ------------------------ arguments ---------------------------------
    character(len=*), intent(in) :: fname           !history file to open and read
    integer, intent(in)  :: kmo, itim               !current month and time index
! --------------------------------------------------------------------

! ------------------------ local variables ---------------------------
    integer :: i,j,k,n                  !do loop indices
    integer :: ier                      !error status
    integer :: varid                    !netCDF variable id
    integer :: status                   !netCDF error status
    integer :: beg3d(3)                 !netCDF 3-d start index (where to read first value)
    integer :: len3d(3)                 !netCDF 3-d count index (number of values to read)

! atm input field names 

    character(len=8) :: fldlst(14)      !name of possible atm fields in input file
    data fldlst( 1) /'TBOT    '/
    data fldlst( 2) /'WIND    '/
    data fldlst( 3) /'QBOT    '/ 
    data fldlst( 4) /'Tdew    '/ 
    data fldlst( 5) /'RH      '/ 
    data fldlst( 6) /'ZBOT    '/
    data fldlst( 7) /'PSRF    '/
    data fldlst( 8) /'FSDS    '/
    data fldlst( 9) /'FSDSdir '/
    data fldlst(10) /'FSDSdif '/
    data fldlst(11) /'FLDS    '/
    data fldlst(12) /'PRECTmms'/
    data fldlst(13) /'PRECCmms'/
    data fldlst(14) /'PRECLmms'/

    real(r8) ea                    !atmospheric emissivity

    logical atmread_err

! use polynomials to calculate saturation vapor pressure and derivative with
! respect to temperature: over water when t > 0 c and over ice when t <= 0 c
! required to convert relative humidity to specific humidity

    real(r8) esatw                 !saturation vapor pressure over water (Pa)
    real(r8) esati                 !saturation vapor pressure over ice (Pa)
    real(r8) e                     !vapor pressure (Pa)
    real(r8) qsat                  !saturation specific humidity (kg/kg)
    real(r8) a0,a1,a2,a3,a4,a5,a6  !coefficients for esat over water
    real(r8) b0,b1,b2,b3,b4,b5,b6  !coefficients for esat over ice
    real(r8) tdc, t                !Kelvins to Celcius function and its input

    parameter (a0=6.107799961    , a1=4.436518521e-01, &
               a2=1.428945805e-02, a3=2.650648471e-04, &
               a4=3.031240396e-06, a5=2.034080948e-08, &
               a6=6.136820929e-11)

    parameter (b0=6.109177956    , b1=5.034698970e-01, &
               b2=1.886013408e-02, b3=4.176223716e-04, &
               b4=5.824720280e-06, b5=4.838803174e-08, &
               b6=1.838826904e-10)


! function declarations

    tdc(t) = min( 50._r8, max(-50._r8,(t-SHR_CONST_TKFRZ)) )
    esatw(t) = 100.*(a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))))
    esati(t) = 100.*(b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))))

! --------------------------------------------------------------------

! --------------------------------------------------------------------
! Read single-level fields
!
! this is done at every subsequent call to atm_readdata
! this includes a second call at the 1st tstep of the run or of the month
!
!    Use wrap_get_vara_realx (a wrapper to nf_get_vara_real)
!    to read data for variable referenced by
!
!    variable id = [i]
!
!    nf_get_vara_real (ncid, i, beg, len, varval)
!    integer beg - a vector of integers specifying the index in the
!                  variable where the first data value is read.
!                  must be dimensioned same as the variable's dimension
!                  and a starting value must be given for each dimension
!    integer len - a vector of integers specifying the number of data
!                  values, for each of the variable's dimensions, to read.
!                  must be dimensioned same as the variable's dimension
!                  and a length must be given for each dimension
! --------------------------------------------------------------------

    if (masterproc) then

! initialize fields to the flag value

       x(:,:,:) = -1.         

! read input data single-level fields

       beg3d(1) = 1     ;  len3d(1) = atmlon
       beg3d(2) = 1     ;  len3d(2) = atmlat
       beg3d(3) = itim  ;  len3d(3) = 1
       do k = 1, 14
          do n = 1, nvar
             if (varnam(n) == fldlst(k)) then
                call wrap_get_vara_realx (ncid,n,beg3d,len3d,x(1,1,k))
             end if
          end do              !end loop of fields in input file
       end do                 !end loop of fields expected in input file

! Close file at the end of the month
! NOTE: as written will not close file if run ends mid-month

       if (itim == ntim) then
          call wrap_close (ncid)
          write (6,*) '---------------------------------------'
          write (6,*) 'ATM_READDATA: closing data for ',trim(fname)
          write (6,*) '---------------------------------------'
          write (6,*)
       end if

    endif	!end of if-masterproc block

#if (defined SPMD) 
    call mpi_bcast (x, size(x), mpir8, 0, mpicom, ier)
#endif

! ----------------------------------------------------------------------
! Determine 2d atmospheric fields
! Follow order in fldlst(14) to determine what was read and what was not 
! ----------------------------------------------------------------------

! Loop over atmospheric longitudes and latitudes

    atmread_err = .false.
!$OMP PARALLEL DO PRIVATE (i,j,e,ea,qsat)
    do j = 1, atmlat
       do i = 1, atmlon

! FORC_TXY
          if (nint(x(i,j,1)) == -1) then
             write(6,*)'ATM error: TBOT has not been read by atm_readdata'
             atmread_err = .true.
          else if (x(i,j,1) < 50.) then
             write(6,*)'ATM error: TBOT appears to be in deg C'
             write(6,*)'Converting to Kelvins now'
             forc_txy_a(i,j) = x(i,j,1) + tfrz
          else
             forc_txy_a(i,j) = x(i,j,1)
          end if

! FORC_UXY, FORC_VXY
          if (nint(x(i,j,2)) == -1) then
             write(6,*)'ATM error: WIND has not been read by atm_readdata'
             atmread_err = .true.
          else
             forc_uxy_a(i,j) = x(i,j,2) / sqrt(2.)
             forc_vxy_a(i,j) = x(i,j,2) / sqrt(2.)
          end if

! FORC_PSRFXY, FORC_PBOTXY
          if (nint(x(i,j,7)) == -1) then
             forc_psrfxy_a(i,j) = SHR_CONST_PSTD
          else
             forc_psrfxy_a(i,j) = x(i,j,7)
          end if
          forc_pbotxy_a(i,j)  = forc_psrfxy_a(i,j)

!FORC_QXY
          if (nint(x(i,j,3)) == -1) then
             if (nint(x(i,j,4)) == -1) then
                if (nint(x(i,j,5)) == -1) then
                   write(6,*)'ATM error: Humidity has not been'
                   write(6,*)'read by atm_readdata'
                   atmread_err = .true.
                else          !using RH as %
                   if (forc_txy_a(i,j) > tfrz) then
                      e = x(i,j,5)/100. * esatw(tdc(forc_txy_a(i,j)))
                   else
                      e = x(i,j,5)/100. * esati(tdc(forc_txy_a(i,j)))
                   end if
                end if
                forc_qxy_a(i,j) = 0.622*e / (forc_pbotxy_a(i,j) - 0.378*e)
             else             !using Tdew
                if (x(i,j,4) < 50.) then
                   write(6,*)'ATM warning: Tdew appears to be in'
                   write(6,*)'deg C, so converting to Kelvin'
                   x(i,j,4) = x(i,j,4) + tfrz
                end if
                if (x(i,j,4) > forc_txy_a(i,j)) then
                   write(6,*)'ATM warning: Dewpt temp > temp!'
                end if
                if (x(i,j,4) > tfrz) then
                   e = esatw(tdc(x(i,j,4)))
                else
                   e = esati(tdc(x(i,j,4)))
                end if
                forc_qxy_a(i,j) = 0.622*e / (forc_pbotxy_a(i,j) - 0.378*e)
             end if
          else                !using QBOT in kg/kg
             if (forc_txy_a(i,j) > tfrz) then
                e = esatw(tdc(forc_txy_a(i,j)))
             else
                e = esati(tdc(forc_txy_a(i,j)))
             end if
             qsat = 0.622*e / (forc_pbotxy_a(i,j) - 0.378*e)
             if (qsat < x(i,j,3)) then
                forc_qxy_a(i,j) = qsat
!                 write(6,*)'ATM warning: qsat < q!'
             else
                forc_qxy_a(i,j) = x(i,j,3)
             end if
          end if

! ZGCMXY
          if (nint(x(i,j,6)) == -1) then
             zgcmxy_a(i,j) = 30.
          else
             zgcmxy_a(i,j) = x(i,j,6)
          end if

! FORC_SOLSXY, FORC_SOLLXY, FORC_SOLSDXY, FORC_SOLLDXY