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