atmdrvmod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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