#include <misc.h>#include <preproc.h>module atmdrvMod#if (defined OFFLINE)!----------------------------------------------------------------------- ! ! Purpose: ! Read atmospheric grid data and generate atmospheric forcing!! Method: ! ! Author: Gordon Bonan, Sam Levis and Mariana Vertenstein! !-----------------------------------------------------------------------! $Id: atmdrvMod.F90,v 1.6.6.7 2002/05/13 18:00:01 erik Exp $!-----------------------------------------------------------------------  use precision  use shr_const_mod, only : SHR_CONST_TKFRZ,SHR_CONST_PSTD  use clm_varpar   , only : lsmlon, lsmlat  use spmdMod      , only : masterproc  implicit none! logical variables for file manipuation  logical :: open_data=.true.             !true => open data file (first tstep of the run or month)  logical :: allocated_data=.false.       !true => allocate dynamic data! atmospheric grid data  integer  :: atmlon                      !number of atm longitudes  integer  :: atmlat                      !number of atm latitudes  real(r8) :: edge_a(4)                   !N,E,S,W edges of atm grid  integer , allocatable :: numlon_a(:)    !number of lon points at each lat  real(r8), allocatable :: latixy_a(:,:)  !latitude of grid cell (degrees)  real(r8), allocatable :: longxy_a(:,:)  !longitude of grid cell (degrees)! atmospheric forcing variables on atmospheric grid  real(r8), allocatable :: x(:,:,:)            !temp. array in which atm data is stored  real(r8), allocatable :: forc_txy_a (:,:)    !atm bottom level temperature (Kelvin)  real(r8), allocatable :: forc_uxy_a (:,:)    !atm bottom level zonal wind (m/s)  real(r8), allocatable :: forc_vxy_a (:,:)    !atm bottom level meridional wind (m/s)  real(r8), allocatable :: forc_qxy_a (:,:)    !atm bottom level specific humidity (kg/kg)  real(r8), allocatable :: zgcmxy_a (:,:)      !atm bottom level height above surface (m)  real(r8), allocatable :: prcxy_a  (:,:)      !convective precipitation rate (mm H2O/s)  real(r8), allocatable :: prlxy_a  (:,:)      !large-scale precipitation rate (mm H2O/s)  real(r8), allocatable :: flwdsxy_a(:,:)      !downward longwave rad onto surface (W/m**2)  real(r8), allocatable :: forc_solsxy_a (:,:) !vis direct beam solar rad onto srf (W/m**2)  real(r8), allocatable :: forc_sollxy_a (:,:) !nir direct beam solar rad onto srf (W/m**2)  real(r8), allocatable :: forc_solsdxy_a(:,:) !vis diffuse solar rad onto srf (W/m**2)  real(r8), allocatable :: forc_solldxy_a(:,:) !nir diffuse solar rad onto srf(W/m**2)  real(r8), allocatable :: forc_pbotxy_a (:,:) !atm bottom level pressure (Pa)  real(r8), allocatable :: forc_psrfxy_a (:,:) !atm surface pressure (Pa)! atmospheric forcing variables on land model grid  real(r8) :: forc_txy (lsmlon,lsmlat)         !atm bottom level temperature (Kelvin)  real(r8) :: forc_uxy (lsmlon,lsmlat)         !atm bottom level zonal wind (m/s)  real(r8) :: forc_vxy (lsmlon,lsmlat)         !atm bottom level meridional wind (m/s)  real(r8) :: forc_qxy (lsmlon,lsmlat)         !atm bottom level specific humidity (kg/kg)  real(r8) :: zgcmxy (lsmlon,lsmlat)           !atm bottom level height above surface (m)  real(r8) :: prcxy  (lsmlon,lsmlat)           !convective precipitation rate (mm H2O/s)  real(r8) :: prlxy  (lsmlon,lsmlat)           !large-scale precipitation rate (mm H2O/s)  real(r8) :: flwdsxy(lsmlon,lsmlat)           !downward longwave rad onto surface (W/m**2)  real(r8) :: forc_solsxy (lsmlon,lsmlat)      !vis direct beam solar rad onto srf (W/m**2)  real(r8) :: forc_sollxy (lsmlon,lsmlat)      !nir direct beam solar rad onto srf (W/m**2)  real(r8) :: forc_solsdxy(lsmlon,lsmlat)      !vis diffuse solar rad onto srf (W/m**2)  real(r8) :: forc_solldxy(lsmlon,lsmlat)      !nir diffuse solar rad onto srf(W/m**2)  real(r8) :: forc_pbotxy (lsmlon,lsmlat)      !atm bottom level pressure (Pa)  real(r8) :: forc_psrfxy (lsmlon,lsmlat)      !atm surface pressure (Pa)! atmosphere grid to land model surface grid mapping for each land grid cell:  integer, parameter :: mxovr =10          !maximum number of overlapping cells  integer :: novr_a2s(lsmlon,lsmlat)       !number    of overlapping atm cells  integer :: iovr_a2s(lsmlon,lsmlat,mxovr) !lon index of overlapping atm cells  integer :: jovr_a2s(lsmlon,lsmlat,mxovr) !lat index of overlapping atm cells  real(r8):: wovr_a2s(lsmlon,lsmlat,mxovr) !weight    of overlapping atm cells! 1d temporarys  real(r8), allocatable, private :: forc_sols(:)   !vis direct beam solar rad onto srf (W/m**2)  real(r8), allocatable, private :: forc_soll(:)   !nir direct beam solar rad onto srf (W/m**2)  real(r8), allocatable, private :: forc_solsd(:)  !vis diffuse solar rad onto srf (W/m**2)  real(r8), allocatable, private :: forc_solld(:)  !nir diffuse solar rad onto srf(W/m**2)! file netCDF id's  integer :: ncid                !netCDF dataset id  integer :: nvar                !number of variables in the data file  integer :: nlon                !number of atm longitude points  integer :: nlat                !number of atm latitude points  integer :: ntim                !number of atm time slices per data file  character(len=8) :: varnam(99) !variable names of atm. fields  SAVE!=======================================================================CONTAINS!=======================================================================  subroutine atmdrv(nstep)!----------------------------------------------------------------------- ! ! Purpose: ! read atmospheric data and generate required atmospheric forcing !! Method: ! This code reads in atmospheric fields from an input file and generates ! the required atmospheric forcing. These data files have [atmmin] minute ! average data for each month. Input data files are named in month-year ! format (e.g., 09-0001 contains 240 3-hour time slices of data, 30*8, for ! September of year one). The model will cycle through however many full ! years of data are available [pyr]. At least one full year of data is! necessary for cycling. The model may start on any base date, as long as! this date corresponds to an existing data file. A run need not be an! exact multiple of a year.!! ============================! Possible atmospheric fields:! ============================! Name     Description                              Required/Optional! -----------------------------------------------------------------------------! TBOT     temperature (K)                          Required! WIND     wind:sqrt(u**2+v**2) (m/s)               Required! QBOT     specific humidity (kg/kg)                Required! Tdew     dewpoint temperature (K)                 Alternative to Q! RH       relative humidity (percent)              Alternative to Q! ZBOT     reference height (m)                     optional ! PSRF     surface pressure (Pa)                    optional! FSDS     total incident solar radiation (W/m**2)  Required! FSDSdir  direct incident solar radiation (W/m**2) optional (replaces FSDS)! FSDSdif  diffuse incident solar rad (W/m**2)      optional (replaces FSDS)! FLDS     incident longwave radiation (W/m**2)     optional! PRECTmms total precipitation (mm H2O / sec)       Required! PRECCmms convective precipitation (mm H2O / sec)  optional (replaces PRECT)! PRECLmms large-scale precipitation (mm H2O / sec) optional (replaces PRECT)! ! ============! Data format:! ============! Data format is netCDF with dimensions longitude x latitude! for each time slice and field. Variable names can be as in above list! or can be reset to desired names using [fldlst] in code below.!! ===============! Namelist input:! ===============! character*256 offline_atmdir = directory for input atm data files !                                (can be Mass Store)! ! Author: Sam Levis! !-----------------------------------------------------------------------    use precision    use infnan    use clm_varder    use clm_varpar  , only : lsmlon, lsmlat    use clm_varctl  , only : offline_atmdir    use clm_varcon  , only : rair, cpair, po2, pco2, tcrit, tfrz    use clm_varmap  , only : begpatch, endpatch    use time_manager, only : get_step_size, get_curr_calday, get_curr_date    use fileutils   , only : getfil    implicit none    include 'netcdf.inc'! ------------------------ arguments ------------------------------    integer, intent(in) :: nstep    !current time step! -----------------------------------------------------------------! ------------------------ local variables ------------------------    integer :: i,j,k                  !loop indices    integer :: itimlast               !last time index used in atm_readdata    real(r8):: calday                 !calendar day at Greenwich (1.00 -> 365.99)    integer :: kda                    !day (1 -> 31)    integer :: kmo                    !month (1 -> 12)    integer :: kyr                    !year (0 -> ...)    integer :: ksec                   !current seconds of current date (0 -> 86400)    integer :: mcdate                 !current date in integer format [yyyymmdd]    integer :: mcsec                  !current time of day [seconds]    integer :: dtime                  !time step size      integer :: minpday = 1440         !minutes per day    integer :: secpmin = 60           !seconds per minute    integer :: ki,kf                  !temporaries	    integer, SAVE :: itim             !time index used in atm_readdata    integer, SAVE :: atmmin           !temporal resolution of atm data (in minutes)    character(len=256), SAVE :: locfn !full file name in case atmdir is in MSS    real(r8) :: temp(begpatch:endpatch)     real(r8) :: forc_precc(begpatch:endpatch)    real(r8) :: forc_precl(begpatch:endpatch)! -----------------------------------------------------------------! -----------------------------------------------------------------! allocate space for temporaries! -----------------------------------------------------------------    if (.not. allocated(forc_sols )) then       allocate (forc_sols(begpatch:endpatch )) ; forc_sols(:) = inf       endif    if (.not. allocated(forc_soll )) then       allocate (forc_soll(begpatch:endpatch )) ; forc_soll(:) = inf       endif    if (.not. allocated(forc_solsd)) then       allocate (forc_solsd(begpatch:endpatch)) ; forc_solsd(:) = inf       endif    if (.not. allocated(forc_solld)) then       allocate (forc_solld(begpatch:endpatch)) ; forc_solld(:) = inf       endif! -----------------------------------------------------------------! Open netcdf file and read data every [atmmin] minutes! -----------------------------------------------------------------! Calendar information for current [nstep]    dtime = get_step_size()    calday = get_curr_calday()    call get_curr_date(kyr, kmo, kda, mcsec)    mcdate = kyr*10000 + kmo*100 + kda! If 1st tstep of the run or of the month, then enter next if-block! Rest flag to open file and set flag to read data    if (open_data) then       call atm_openfile (kda, kmo, kyr, locfn, itim, atmmin)    endif! Calculate time index    itimlast = itim    itim = 1 + ((kda - 1)*minpday + (mcsec - dtime)/secpmin)/atmmin    if (dtime == int(atmmin*secpmin)) then          if (kda == 1 .and. mcsec == 0) itim = itimlast+1    else       if (kda == 1 .and. mcsec == 0) itim = itimlast    endif! Determine if new data is to be read    if (open_data .or. mod(nstep-1,atmmin*secpmin/dtime) ==0 ) then! Read data for current time slice       if (masterproc) then          write (6,*)          write (6,'(72a1)') ("-",i=1,60)          write (6,*)'nstep= ',nstep,' date= ',mcdate,' sec= ',mcsec          if ( len_trim(locfn) > 0 )then             write (6,*)'ATM: attempting to read data from ',trim(locfn)          end if          write (6,'(72a1)') ("-",i=1,60)          write (6,*)       endif                 call atm_readdata (locfn, kmo, itim)! Map 2d atmospheric fields from atmospheric grid to land model surface grid. ! Area-average absolute value of winds (i.e., regardless of! direction) since land model cares about magnitude not direction.! Then need to adjust resultant stresses for direction of wind.       call interpa2s (forc_txy_a    , forc_txy    , zgcmxy_a      , zgcmxy      , &                       forc_uxy_a    , forc_uxy    , forc_vxy_a    , forc_vxy    , &                       forc_qxy_a    , forc_qxy    , prcxy_a       , prcxy       , &                       prlxy_a       , prlxy       , flwdsxy_a     , flwdsxy     , &                       forc_solsxy_a , forc_solsxy , forc_sollxy_a , forc_sollxy , &                       forc_solsdxy_a, forc_solsdxy, forc_solldxy_a, forc_solldxy, &                       forc_pbotxy_a , forc_pbotxy , forc_psrfxy_a , forc_psrfxy )! Map atmospheric fields to force land model: ! [lsmlon] x [lsmlat] grid -> [numpatch] vector of subgrid patches       ki = begpatch; kf = endpatch       call xy2v (lsmlon,lsmlat,forc_txy    ,ki,kf,temp ) ; clm(ki:kf)%forc_t   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_uxy    ,ki,kf,temp ) ; clm(ki:kf)%forc_u   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_vxy    ,ki,kf,temp ) ; clm(ki:kf)%forc_v   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_qxy    ,ki,kf,temp ) ; clm(ki:kf)%forc_q   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,zgcmxy      ,ki,kf,temp ) ; clm(ki:kf)%forc_hgt = temp(ki:kf)       call xy2v (lsmlon,lsmlat,prcxy       ,ki,kf,temp ) ; forc_precc(ki:kf)   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,prlxy       ,ki,kf,temp ) ; forc_precl(ki:kf)   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,flwdsxy     ,ki,kf,temp ) ; clm(ki:kf)%forc_lwrad = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_pbotxy ,ki,kf,temp ) ; clm(ki:kf)%forc_pbot   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_psrfxy ,ki,kf,temp ) ; clm(ki:kf)%forc_psrf   = temp(ki:kf)       call xy2v (lsmlon,lsmlat,forc_solsxy ,ki,kf,forc_sols )        call xy2v (lsmlon,lsmlat,forc_sollxy ,ki,kf,forc_soll )        call xy2v (lsmlon,lsmlat,forc_solsdxy,ki,kf,forc_solsd)        call xy2v (lsmlon,lsmlat,forc_solldxy,ki,kf,forc_solld) !$OMP PARALLEL DO PRIVATE (k)       do k = begpatch, endpatch          clm(k)%forc_hgt_u = clm(k)%forc_hgt   !observational height of wind [m]           clm(k)%forc_hgt_t = clm(k)%forc_hgt   !observational height of temperature [m]           clm(k)%forc_hgt_q = clm(k)%forc_hgt   !observational height of humidity [m] ! Set upper limit of air temperature for snowfall at 275.65K.! This cut-off was selected based on Fig. 1, Plate 3-1, of Snow! Hydrology (1956).          if (forc_precc(k)+forc_precl(k) > 0.) then             if (clm(k)%forc_t > (tfrz + tcrit)) then                clm(k)%itypprc   = 1                clm(k)%forc_rain = forc_precc(k)+forc_precl(k)                clm(k)%forc_snow = 0.             else                clm(k)%itypprc   = 2