#ifdef AIX
@PROCESS ALIAS_SIZE(805306368)
#endif
module datm_comp_mod 1,16
! !USES:
use shr_const_mod
use shr_sys_mod
use shr_kind_mod , only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, &
CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
use shr_file_mod , only: shr_file_getunit, shr_file_getlogunit, shr_file_getloglevel, &
shr_file_setlogunit, shr_file_setloglevel, shr_file_setio, &
shr_file_freeunit
use shr_cal_mod , only: shr_cal_date2eday
use shr_mpi_mod , only: shr_mpi_bcast
use mct_mod
use esmf_mod
use perf_mod
use shr_strdata_mod
use shr_dmodel_mod
use shr_pcdf_mod
use datm_shr_mod
use seq_cdata_mod
use seq_infodata_mod
use seq_timemgr_mod
use seq_flds_indices , only: nflds_a2x, &
nflds_x2a
use seq_flds_mod , only: seq_flds_a2x_fields, &
seq_flds_x2a_fields
!
! !PUBLIC TYPES:
implicit none
private ! except
!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------
public :: datm_comp_init
public :: datm_comp_run
public :: datm_comp_final
!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------
!--- other ---
character(CS) :: myModelName = 'atm' ! user defined model name
integer(IN) :: mpicom
integer(IN) :: COMPID ! mct comp id
integer(IN) :: my_task ! my task in mpi communicator mpicom
integer(IN) :: npes ! total number of tasks
integer(IN),parameter :: master_task=0 ! task number of master task
integer(IN) :: logunit ! logging unit number
character(CL) :: atm_mode ! mode
integer(IN) :: dbug = 0 ! debug level (higher is more)
logical :: firstcall = .true. ! first call logical
logical :: scmMode = .false. ! single column mode
real(R8) :: scmLat = shr_const_SPVAL ! single column lat
real(R8) :: scmLon = shr_const_SPVAL ! single column lon
integer :: phase ! phase of method
logical :: read_restart ! start from restart
real(R8) :: orbEccen ! orb eccentricity (unit-less)
real(R8) :: orbMvelpp ! orb moving vernal eq (radians)
real(R8) :: orbLambm0 ! orb mean long of perhelion (radians)
real(R8) :: orbObliqr ! orb obliquity (radians)
real(R8) :: tbotmax ! units detector
real(R8) :: tdewmax ! units detector
real(R8) :: anidrmax ! existance detector
integer(IN) :: iradsw ! radiation logical
character(CL) :: factorFn ! file containing correction factors
character(len=*),parameter :: rpfile = 'rpointer.atm'
character(len=*),parameter :: nullstr = 'undefined'
real(R8),parameter :: aerodep_spval = 1.e29_r8 ! special aerosol deposition
real(R8),parameter :: tKFrz = SHR_CONST_TKFRZ
real(R8),parameter :: degtorad = SHR_CONST_PI/180.0_R8
real(R8),parameter :: pstd = SHR_CONST_PSTD ! standard pressure ~ Pa
real(R8),parameter :: stebol = SHR_CONST_STEBOL ! Stefan-Boltzmann constant ~ W/m^2/K^4
real(R8),parameter :: rdair = SHR_CONST_RDAIR ! dry air gas constant ~ J/K/kg
real(R8),parameter :: avg_c0 = 61.846_R8
real(R8),parameter :: avg_c1 = 1.107_R8
real(R8),parameter :: amp_c0 = -21.841_R8
real(R8),parameter :: amp_c1 = -0.447_R8
real(R8),parameter :: phs_c0 = 0.298_R8
real(R8),parameter :: dLWarc = -5.000_R8
real(R8) ,save :: dTarc(12)
data dTarc / 0.49_R8, 0.06_R8,-0.73_R8, -0.89_R8,-0.77_R8,-1.02_R8, &
& -1.99_R8,-0.91_R8, 1.72_R8, 2.30_R8, 1.81_R8, 1.06_R8/
integer(IN) :: kz,ku,kv,ktbot,kptem,kshum,kdens,kpbot,kpslv,klwdn
integer(IN) :: krc,krl,ksc,ksl,kswndr,kswndf,kswvdr,kswvdf,kswnet,kco2p,kco2d
integer(IN) :: kbid,kbod,kbiw,koid,kood,koiw,kdw1,kdw2,kdw3,kdw4,kdd1,kdd2,kdd3,kdd4
integer(IN) :: kanidr,kanidf,kavsdr,kavsdf
integer(IN) :: stbot,swind,sz,spbot,sshum,stdew,srh,slwdn,sswdn,sswdndf,sswdndr
integer(IN) :: sprecc,sprecl,sprecn,sco2p,sco2d,sswup,sprec
type(shr_strdata_type) :: SDATM
type(mct_rearr) :: rearr
type(mct_avect) :: avstrm ! av of data from stream
integer(IN), pointer :: imask(:)
real(R8), pointer :: yc(:)
real(R8), pointer :: windFactor(:)
real(R8), pointer :: winddFactor(:)
real(R8), pointer :: qsatFactor(:)
integer(IN),parameter :: ktrans = 56
character(16),parameter :: avofld(1:ktrans) = &
(/"Sa_z ","Sa_u ","Sa_v ","Sa_tbot ", &
"Sa_ptem ","Sa_shum ","Sa_dens ","Sa_pbot ", &
"Sa_pslv ","Faxa_lwdn ","Faxa_rainc ","Faxa_rainl ", &
"Faxa_snowc ","Faxa_snowl ","Faxa_swndr ","Faxa_swvdr ", &
"Faxa_swndf ","Faxa_swvdf ","Faxa_swnet ","Sa_co2prog ", &
"Sa_co2diag ","Faxa_bcphidry ","Faxa_bcphodry ","Faxa_bcphiwet ", &
"Faxa_ocphidry ","Faxa_ocphodry ","Faxa_ocphiwet ","Faxa_dstwet1 ", &
"Faxa_dstwet2 ","Faxa_dstwet3 ","Faxa_dstwet4 ","Faxa_dstdry1 ", &
"Faxa_dstdry2 ","Faxa_dstdry3 ","Faxa_dstdry4 ", &
"Sx_tref ","Sx_qref ","Sx_avsdr ","Sx_anidr ", &
"Sx_avsdf ","Sx_anidf ","Sx_t ","So_t ", &
"Sx_snowh ","Sx_lfrac ","Sx_ifrac ","Sx_ofrac ", &
"Faxx_taux ","Faxx_tauy ","Faxx_lat ","Faxx_sen ", &
"Faxx_lwup ","Faxx_evap ","Faxx_fco2_lnd ","Faxx_fco2_ocn ", &
"Faxx_fdms " /)
character(16),parameter :: avifld(1:ktrans) = &
(/"z ","u ","v ","tbot ", &
"ptem ","shum ","dens ","pbot ", &
"pslv ","lwdn ","rainc ","rainl ", &
"snowc ","snowl ","swndr ","swvdr ", &
"swndf ","swvdf ","swnet ","co2prog ", &
"co2diag ","bcphidry ","bcphodry ","bcphiwet ", &
"ocphidry ","ocphodry ","ocphiwet ","dstwet1 ", &
"dstwet2 ","dstwet3 ","dstwet4 ","dstdry1 ", &
"dstdry2 ","dstdry3 ","dstdry4 ", &
"tref ","qref ","avsdr ","anidr ", &
"avsdf ","anidf ","ts ","to ", &
"snowh ","lfrac ","ifrac ","ofrac ", &
"taux ","tauy ","lat ","sen ", &
"lwup ","evap ","co2lnd ","co2ocn ", &
"dms " /)
integer(IN),parameter :: ktranss = 18
character(16),parameter :: stofld(1:ktranss) = &
(/"strm_tbot ","strm_wind ","strm_z ","strm_pbot ", &
"strm_shum ","strm_tdew ","strm_rh ","strm_lwdn ", &
"strm_swdn ","strm_swdndf ","strm_swdndr ","strm_precc ", &
"strm_precl ","strm_precn ","strm_co2prog ","strm_co2diag ", &
"strm_swup ","strm_prec " /)
character(16),parameter :: stifld(1:ktranss) = &
(/"tbot ","wind ","z ","pbot ", &
"shum ","tdew ","rh ","lwdn ", &
"swdn ","swdndf ","swdndr ","precc ", &
"precl ","precn ","co2prog ","co2diag ", &
"swup ","prec " /)
save
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: datm_comp_init
!
! !DESCRIPTION:
! initialize data atm model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
subroutine datm_comp_init( EClock, cdata, x2a, a2x, NLFilename ) 2,49
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock) , intent(in) :: EClock
type(seq_cdata) , intent(inout) :: cdata
type(mct_aVect) , intent(inout) :: x2a, a2x
character(len=*), optional , intent(in) :: NLFilename ! Namelist filename
!EOP
!--- local variables ---
integer(IN) :: n,k ! generic counters
integer(IN) :: ierr ! error code
integer(IN) :: gsize ! global size
integer(IN) :: lsize ! local size
integer(IN) :: shrlogunit, shrloglev ! original log unit and level
integer(IN) :: nunit ! unit number
integer(IN) :: kmask ! field reference
integer(IN) :: klat ! field reference
integer(IN) :: kfld ! fld index
integer(IN) :: cnt ! counter
logical :: atm_present ! flag
logical :: atm_prognostic ! flag
type(seq_infodata_type), pointer :: infodata
type(mct_gsMap) , pointer :: gsMap
type(mct_gGrid) , pointer :: ggrid
character(CL) :: filePath ! generic file path
character(CL) :: fileName ! generic file name
character(CS) :: timeName ! domain file: time variable name
character(CS) :: lonName ! domain file: lon variable name
character(CS) :: latName ! domain file: lat variable name
character(CS) :: maskName ! domain file: mask variable name
character(CS) :: areaName ! domain file: area variable name
integer(IN) :: yearFirst ! first year to use in data stream
integer(IN) :: yearLast ! last year to use in data stream
integer(IN) :: yearAlign ! data year that aligns with yearFirst
character(CL) :: atm_in ! dshr atm namelist
character(CL) :: decomp ! decomp strategy
character(CL) :: rest_file ! restart filename
character(CL) :: rest_file_strm ! restart filename for streams
character(CL) :: restfilm ! model restart file namelist
character(CL) :: restfils ! stream restart file namelist
logical :: exists ! filename existance
integer(IN) :: nu ! unit number
integer(IN) :: idt ! integer timestep
integer(IN) :: CurrentYMD ! model date
integer(IN) :: CurrentTOD ! model sec into model date
integer(IN) :: stepno ! step number
real(R8) :: nextsw_cday ! calendar of next atm sw
character(CL) :: flds_strm
logical :: presaero ! true => send valid prescribe aero fields to coupler
!----- define namelist -----
namelist / datm_nml / &
atm_in, decomp, iradsw, factorFn, restfilm, restfils, presaero
!--- formats ---
character(*), parameter :: F00 = "('(datm_comp_init) ',8a)"
character(*), parameter :: F0L = "('(datm_comp_init) ',a, l)"
character(*), parameter :: F01 = "('(datm_comp_init) ',a,5i8)"
character(*), parameter :: F02 = "('(datm_comp_init) ',a,4es13.6)"
character(*), parameter :: F03 = "('(datm_comp_init) ',a,i8,a)"
character(*), parameter :: F04 = "('(datm_comp_init) ',2a,2i8,'s')"
character(*), parameter :: F05 = "('(datm_comp_init) ',a,2f10.4)"
character(*), parameter :: F90 = "('(datm_comp_init) ',73('='))"
character(*), parameter :: F91 = "('(datm_comp_init) ',73('-'))"
character(*), parameter :: subName = "(datm_comp_init) "
!-------------------------------------------------------------------------------
call t_startf('DATM_INIT')
! Set cdata pointers
call seq_cdata_setptrs(cdata, ID=COMPID, mpicom=mpicom, &
gsMap=gsmap, dom=ggrid, infodata=infodata)
call seq_infodata_getData(infodata,atm_phase=phase)
if (phase == 1) then
! Determine communicator groups and sizes
call mpi_comm_rank(mpicom, my_task, ierr)
call mpi_comm_size(mpicom, npes, ierr)
!--- open log file ---
if (my_task == master_task) then
logUnit = shr_file_getUnit()
call shr_file_setIO('atm_modelio.nml',logUnit)
else
logUnit = 6
endif
endif
!----------------------------------------------------------------------------
! Reset shr logging to my log file
!----------------------------------------------------------------------------
call shr_file_getLogUnit (shrlogunit)
call shr_file_getLogLevel(shrloglev)
call shr_file_setLogUnit (logUnit)
if (phase == 1) then
!----------------------------------------------------------------------------
! Set a Few Defaults
!----------------------------------------------------------------------------
call seq_infodata_getData(infodata,single_column=scmMode, &
& scmlat=scmlat, scmlon=scmLon)
call seq_infodata_GetData(infodata,orb_eccen=orbEccen,orb_mvelpp=orbMvelpp, &
orb_lambm0=orbLambm0,orb_obliqr=orbObliqr )
atm_present = .false.
atm_prognostic = .false.
call seq_infodata_GetData(infodata,read_restart=read_restart)
!----------------------------------------------------------------------------
! Read datm_in
!----------------------------------------------------------------------------
call t_startf('datm_readnml')
filename = "datm_in"
atm_in = "unset"
decomp = "1d"
iradsw = 0
factorFn = 'null'
restfilm = trim(nullstr)
restfils = trim(nullstr)
presaero = .false.
if (my_task == master_task) then
nunit = shr_file_getUnit() ! get unused unit number
open (nunit,file=trim(filename),status="old",action="read")
read (nunit,nml=datm_nml,iostat=ierr)
close(nunit)
call shr_file_freeUnit(nunit)
if (ierr > 0) then
write(logunit,F01) 'ERROR: reading input namelist, '//trim(filename)//' iostat=',ierr
call shr_sys_abort(subName//': namelist read error '//trim(filename))
end if
write(logunit,F00)' atm_in = ',trim(atm_in)
write(logunit,F00)' decomp = ',trim(decomp)
write(logunit,F01)' iradsw = ',iradsw
write(logunit,F00)' factorFn = ',trim(factorFn)
write(logunit,F00)' restfilm = ',trim(restfilm)
write(logunit,F00)' restfils = ',trim(restfils)
write(logunit,F0L)' presaero = ',presaero
endif
call shr_mpi_bcast(atm_in,mpicom,'atm_in')
call shr_mpi_bcast(decomp,mpicom,'decomp')
call shr_mpi_bcast(iradsw,mpicom,'iradsw')
call shr_mpi_bcast(factorFn,mpicom,'factorFn')
call shr_mpi_bcast(restfilm,mpicom,'restfilm')
call shr_mpi_bcast(restfils,mpicom,'restfils')
call shr_mpi_bcast(presaero,mpicom,'presaero')
rest_file = trim(restfilm)
rest_file_strm = trim(restfils)
!----------------------------------------------------------------------------
! Read dshr namelist
!----------------------------------------------------------------------------
call shr_strdata_readnml(SDATM,trim(atm_in),mpicom=mpicom)
!----------------------------------------------------------------------------
! Validate mode
!----------------------------------------------------------------------------
atm_mode = trim(SDATM%dataMode)
! check that we know how to handle the mode
if (trim(atm_mode) == 'NULL' .or. &
trim(atm_mode) == 'TN460' .or. &
trim(atm_mode) == 'CORE2' .or. &
trim(atm_mode) == 'CLMNCEP' .or. &
trim(atm_mode) == 'CPLHIST' .or. &
trim(atm_mode) == 'COPYALL') then
if (my_task == master_task) &
write(logunit,F00) ' atm mode = ',trim(atm_mode)
else
write(logunit,F00) ' ERROR illegal atm mode = ',trim(atm_mode)
call shr_sys_abort()
endif
call t_stopf('datm_readnml')
!----------------------------------------------------------------------------
! Initialize datasets
!----------------------------------------------------------------------------
call t_startf('datm_strdata_init')
if (trim(atm_mode) /= 'NULL') then
atm_present = .true.
if (scmmode) then
if (my_task == master_task) &
write(logunit,F05) ' scm lon lat = ',scmlon,scmlat
call shr_strdata_init(SDATM,mpicom,compid,name='atm', &
scmmode=scmmode,scmlon=scmlon,scmlat=scmlat)
else
call shr_strdata_init(SDATM,mpicom,compid,name='atm')
endif
!--- overwrite mask and frac ---
k = mct_aVect_indexRA(SDATM%grid%data,'mask')
SDATM%grid%data%rAttr(k,:) = 1.0_R8
k = mct_aVect_indexRA(SDATM%grid%data,'frac')
SDATM%grid%data%rAttr(k,:) = 1.0_R8
!--- set data needed for cosz t-interp method ---
call seq_timemgr_EClockGetData( EClock, dtime=idt )
call shr_strdata_setOrbs(SDATM,orbEccen,orbMvelpp,orbLambm0,orbObliqr,idt)
endif
if (my_task == master_task) then
call shr_strdata_print(SDATM,'ATM data')
endif
call t_stopf('datm_strdata_init')
!----------------------------------------------------------------------------
! Set flag to specify data components
!----------------------------------------------------------------------------
call seq_infodata_PutData(infodata, &
atm_present=atm_present, atm_prognostic=atm_prognostic, &
atm_nx=SDATM%nxg, atm_ny=SDATM%nyg )
call seq_infodata_PutData( infodata, atm_aero=presaero)
!----------------------------------------------------------------------------
! Initialize MCT global seg map, 1d decomp
!----------------------------------------------------------------------------
call t_startf('datm_initgsmaps')
if (my_task == master_task) write(logunit,F00) ' initialize gsmaps'
call shr_sys_flush(logunit)
call shr_dmodel_gsmapcreate(gsmap,SDATM%nxg*SDATM%nyg,compid,mpicom,decomp)
lsize = mct_gsmap_lsize(gsmap,mpicom)
if (atm_present) then
call mct_rearr_init(SDATM%gsmap,gsmap,mpicom,rearr)
endif
call t_stopf('datm_initgsmaps')
!----------------------------------------------------------------------------
! Initialize MCT domain
!----------------------------------------------------------------------------
call t_startf('datm_initmctdom')
if (my_task == master_task) write(logunit,F00) 'copy domains'
call shr_sys_flush(logunit)
if (atm_present) call shr_dmodel_rearrGGrid(SDATM%grid, ggrid, gsmap, rearr, mpicom)
call t_stopf('datm_initmctdom')
!----------------------------------------------------------------------------
! Initialize MCT attribute vectors
!----------------------------------------------------------------------------
call t_startf('datm_initmctavs')
if (my_task == master_task) write(logunit,F00) 'allocate AVs'
call shr_sys_flush(logunit)
call mct_aVect_init(a2x, rList=seq_flds_a2x_fields, lsize=lsize)
call mct_aVect_zero(a2x)
kz = mct_aVect_indexRA(a2x,'Sa_z')
ku = mct_aVect_indexRA(a2x,'Sa_u')
kv = mct_aVect_indexRA(a2x,'Sa_v')
ktbot = mct_aVect_indexRA(a2x,'Sa_tbot')
kptem = mct_aVect_indexRA(a2x,'Sa_ptem')
kshum = mct_aVect_indexRA(a2x,'Sa_shum')
kdens = mct_aVect_indexRA(a2x,'Sa_dens')
kpbot = mct_aVect_indexRA(a2x,'Sa_pbot')
kpslv = mct_aVect_indexRA(a2x,'Sa_pslv')
klwdn = mct_aVect_indexRA(a2x,'Faxa_lwdn')
krc = mct_aVect_indexRA(a2x,'Faxa_rainc')
krl = mct_aVect_indexRA(a2x,'Faxa_rainl')
ksc = mct_aVect_indexRA(a2x,'Faxa_snowc')
ksl = mct_aVect_indexRA(a2x,'Faxa_snowl')
kswndr= mct_aVect_indexRA(a2x,'Faxa_swndr')
kswndf= mct_aVect_indexRA(a2x,'Faxa_swndf')
kswvdr= mct_aVect_indexRA(a2x,'Faxa_swvdr')
kswvdf= mct_aVect_indexRA(a2x,'Faxa_swvdf')
kswnet= mct_aVect_indexRA(a2x,'Faxa_swnet')
kco2p = mct_aVect_indexRA(a2x,'Sa_co2prog',perrWith='quiet')
kco2d = mct_aVect_indexRA(a2x,'Sa_co2diag',perrWith='quiet')
kbid = mct_aVect_indexRA(a2x,'Faxa_bcphidry')
kbod = mct_aVect_indexRA(a2x,'Faxa_bcphodry')
kbiw = mct_aVect_indexRA(a2x,'Faxa_bcphiwet')
koid = mct_aVect_indexRA(a2x,'Faxa_ocphidry')
kood = mct_aVect_indexRA(a2x,'Faxa_ocphodry')
koiw = mct_aVect_indexRA(a2x,'Faxa_ocphiwet')
kdd1 = mct_aVect_indexRA(a2x,'Faxa_dstdry1')
kdd2 = mct_aVect_indexRA(a2x,'Faxa_dstdry2')
kdd3 = mct_aVect_indexRA(a2x,'Faxa_dstdry3')
kdd4 = mct_aVect_indexRA(a2x,'Faxa_dstdry4')
kdw1 = mct_aVect_indexRA(a2x,'Faxa_dstwet1')
kdw2 = mct_aVect_indexRA(a2x,'Faxa_dstwet2')
kdw3 = mct_aVect_indexRA(a2x,'Faxa_dstwet3')
kdw4 = mct_aVect_indexRA(a2x,'Faxa_dstwet4')
call mct_aVect_init(x2a, rList=seq_flds_x2a_fields, lsize=lsize)
call mct_aVect_zero(x2a)
kanidr = mct_aVect_indexRA(x2a,'Sx_anidr')
kanidf = mct_aVect_indexRA(x2a,'Sx_anidf')
kavsdr = mct_aVect_indexRA(x2a,'Sx_avsdr')
kavsdf = mct_aVect_indexRA(x2a,'Sx_avsdf')
!--- figure out what's on the streams ---
cnt = 0
flds_strm = ''
do n = 1,SDATM%nstreams
do k = 1,ktranss
kfld = mct_aVect_indexRA(SDATM%avs(n),trim(stifld(k)),perrWith='quiet')
if (kfld > 0) then
cnt = cnt + 1
if (cnt == 1) then
flds_strm = trim(stofld(k))
else
flds_strm = trim(flds_strm)//':'//trim(stofld(k))
endif
endif
enddo
enddo
if (my_task == master_task) write(logunit,F00) ' flds_strm = ',trim(flds_strm)
call shr_sys_flush(logunit)
call mct_aVect_init(avstrm, rList=flds_strm, lsize=lsize)
call mct_aVect_zero(avstrm)
stbot = mct_aVect_indexRA(avstrm,'strm_tbot',perrWith='quiet')
swind = mct_aVect_indexRA(avstrm,'strm_wind',perrWith='quiet')
sz = mct_aVect_indexRA(avstrm,'strm_z',perrWith='quiet')
spbot = mct_aVect_indexRA(avstrm,'strm_pbot',perrWith='quiet')
sshum = mct_aVect_indexRA(avstrm,'strm_shum',perrWith='quiet')
stdew = mct_aVect_indexRA(avstrm,'strm_tdew',perrWith='quiet')
srh = mct_aVect_indexRA(avstrm,'strm_rh',perrWith='quiet')
slwdn = mct_aVect_indexRA(avstrm,'strm_lwdn',perrWith='quiet')
sswdn = mct_aVect_indexRA(avstrm,'strm_swdn',perrWith='quiet')
sswdndf= mct_aVect_indexRA(avstrm,'strm_swdndf',perrWith='quiet')
sswdndr= mct_aVect_indexRA(avstrm,'strm_swdndr',perrWith='quiet')
sprecc = mct_aVect_indexRA(avstrm,'strm_precc',perrWith='quiet')
sprecl = mct_aVect_indexRA(avstrm,'strm_precl',perrWith='quiet')
sprecn = mct_aVect_indexRA(avstrm,'strm_precn',perrWith='quiet')
sco2p = mct_aVect_indexRA(avstrm,'strm_co2p',perrWith='quiet')
sco2d = mct_aVect_indexRA(avstrm,'strm_co2d',perrWith='quiet')
sswup = mct_aVect_indexRA(avstrm,'strm_swup',perrWith='quiet')
sprec = mct_aVect_indexRA(avstrm,'strm_prec',perrWith='quiet')
allocate(imask(lsize))
allocate(yc(lsize))
allocate(windFactor(lsize))
allocate(winddFactor(lsize))
allocate(qsatFactor(lsize))
kmask = mct_aVect_indexRA(ggrid%data,'mask')
imask(:) = nint(ggrid%data%rAttr(kmask,:))
klat = mct_aVect_indexRA(ggrid%data,'lat')
yc(:) = ggrid%data%rAttr(klat,:)
call t_stopf('datm_initmctavs')
!----------------------------------------------------------------------------
! Read restart
!----------------------------------------------------------------------------
if (read_restart) then
if (trim(rest_file) == trim(nullstr) .and. &
trim(rest_file_strm) == trim(nullstr)) then
if (my_task == master_task) then
write(logunit,F00) ' restart filenames from rpointer'
call shr_sys_flush(logunit)
inquire(file=trim(rpfile),exist=exists)
if (.not.exists) then
write(logunit,F00) ' ERROR: rpointer file does not exist'
call shr_sys_abort(trim(subname)//' ERROR: rpointer file missing')
endif
nu = shr_file_getUnit()
open(nu,file=trim(rpfile),form='formatted')
read(nu,'(a)') rest_file
read(nu,'(a)') rest_file_strm
close(nu)
call shr_file_freeUnit(nu)
inquire(file=trim(rest_file_strm),exist=exists)
endif
call shr_mpi_bcast(rest_file,mpicom,'rest_file')
call shr_mpi_bcast(rest_file_strm,mpicom,'rest_file_strm')
else
! use namelist already read
if (my_task == master_task) then
write(logunit,F00) ' restart filenames from namelist '
call shr_sys_flush(logunit)
inquire(file=trim(rest_file_strm),exist=exists)
endif
endif
call shr_mpi_bcast(exists,mpicom,'exists')
! if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file)
! call shr_pcdf_readwrite('read',trim(rest_file),mpicom,gsmap,rf1=somtp,rf1n='somtp')
if (exists) then
if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file_strm)
call shr_strdata_restRead(trim(rest_file_strm),SDATM,mpicom)
else
if (my_task == master_task) write(logunit,F00) ' file not found, skipping ',trim(rest_file_strm)
endif
call shr_sys_flush(logunit)
endif
if (read_restart) then
call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt )
nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw )
else
call seq_timemgr_EClockGetData( EClock, curr_cday=nextsw_cday, stepno=stepno )
endif
call seq_infodata_PutData(infodata, nextsw_cday=nextsw_cday )
else
call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt)
nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw )
call seq_infodata_PutData(infodata, nextsw_cday=nextsw_cday )
endif
!----------------------------------------------------------------------------
! Set initial atm state, needed for CCSM atm initialization
!----------------------------------------------------------------------------
call t_adj_detailf(+2)
call datm_comp_run( EClock, cdata, x2a, a2x)
call t_adj_detailf(-2)
!----------------------------------------------------------------------------
! Reset shr logging to original values
!----------------------------------------------------------------------------
if (my_task == master_task) write(logunit,F00) 'datm_comp_init done'
call shr_sys_flush(logunit)
call shr_file_setLogUnit (shrlogunit)
call shr_file_setLogLevel(shrloglev)
call shr_sys_flush(logunit)
call t_stopf('DATM_INIT')
end subroutine datm_comp_init
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: datm_comp_run
!
! !DESCRIPTION:
! run method for dead atm model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
subroutine datm_comp_run( EClock, cdata, x2a, a2x) 2,37
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(ESMF_Clock) ,intent(in) :: EClock
type(seq_cdata) ,intent(inout) :: cdata
type(mct_aVect) ,intent(inout) :: x2a ! driver -> dead
type(mct_aVect) ,intent(inout) :: a2x ! dead -> driver
!EOP
!--- local ---
type(mct_gsMap) , pointer :: gsMap
type(mct_gGrid) , pointer :: ggrid
integer(IN) :: CurrentYMD ! model date
integer(IN) :: CurrentTOD ! model sec into model date
integer(IN) :: yy,mm,dd ! year month day
integer(IN) :: n ! indices
integer(IN) :: lsize ! size of attr vect
integer(IN) :: shrlogunit, shrloglev ! original log unit and level
logical :: mssrmlf ! remove old data
integer(IN) :: idt ! integer timestep
real(R8) :: dt ! timestep
logical :: write_restart ! restart now
character(CL) :: case_name ! case name
character(CL) :: rest_file ! restart_file
character(CL) :: rest_file_strm ! restart_file
integer(IN) :: nu ! unit number
integer(IN) :: stepno ! step number
real(R8) :: nextsw_cday ! calendar of next atm sw
integer(IN) :: eday ! elapsed day
real(R8) :: rday ! elapsed day
real(R8) :: cosFactor ! cosine factor
real(R8) :: factor ! generic/temporary correction factor
real(R8) :: avg_alb ! average albedo
real(R8) :: tMin ! minimum temperature
!--- temporaries
real(R8) :: uprime,vprime,swndr,swndf,swvdr,swvdf,ratio_rvrf
real(R8) :: tbot,pbot,rtmp,vp,ea,e,qsat,frac,qsatT
type(seq_infodata_type), pointer :: infodata
character(*), parameter :: F00 = "('(datm_comp_run) ',8a)"
character(*), parameter :: F04 = "('(datm_comp_run) ',2a,2i8,'s')"
character(*), parameter :: subName = "(datm_comp_run) "
!-------------------------------------------------------------------------------
call t_startf('DATM_RUN')
call t_startf('datm_run1')
!----------------------------------------------------------------------------
! Reset shr logging to my log file
!----------------------------------------------------------------------------
call shr_file_getLogUnit (shrlogunit)
call shr_file_getLogLevel(shrloglev)
call shr_file_setLogUnit (logUnit)
call seq_cdata_setptrs(cdata, gsMap=gsmap, dom=ggrid)
call seq_cdata_setptrs(cdata, infodata=infodata)
call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
call seq_timemgr_EClockGetData( EClock, curr_yr=yy, curr_mon=mm, curr_day=dd)
call seq_timemgr_EClockGetData( EClock, stepno=stepno, dtime=idt)
dt = idt * 1.0_r8
write_restart = seq_timemgr_RestartAlarmIsOn(EClock)
call t_stopf('datm_run1')
!--------------------
! ADVANCE ATM
!--------------------
call t_barrierf('datm_BARRIER',mpicom)
call t_startf('datm')
nextsw_cday = datm_shr_getNextRadCDay( CurrentYMD, CurrentTOD, stepno, idt, iradsw )
call seq_infodata_PutData(infodata, nextsw_cday=nextsw_cday )
!--- copy all fields from streams to a2x as default ---
if (trim(atm_mode) /= 'NULL') then
call t_startf('datm_strdata_advance')
call shr_strdata_advance(SDATM,currentYMD,currentTOD,mpicom,'datm')
call t_stopf('datm_strdata_advance')
call t_barrierf('datm_scatter_BARRIER',mpicom)
call t_startf('datm_scatter')
if (trim(atm_mode) /= 'COPYALL') then
lsize = mct_avect_lsize(a2x)
do n = 1,lsize
a2x%rAttr(kbid,n) = aerodep_spval
a2x%rAttr(kbod,n) = aerodep_spval
a2x%rAttr(kbiw,n) = aerodep_spval
a2x%rAttr(koid,n) = aerodep_spval
a2x%rAttr(kood,n) = aerodep_spval
a2x%rAttr(koiw,n) = aerodep_spval
a2x%rAttr(kdd1,n) = aerodep_spval
a2x%rAttr(kdd2,n) = aerodep_spval
a2x%rAttr(kdd3,n) = aerodep_spval
a2x%rAttr(kdd4,n) = aerodep_spval
a2x%rAttr(kdw1,n) = aerodep_spval
a2x%rAttr(kdw2,n) = aerodep_spval
a2x%rAttr(kdw3,n) = aerodep_spval
a2x%rAttr(kdw4,n) = aerodep_spval
enddo
endif
do n = 1,SDATM%nstreams
call shr_dmodel_translateAV(SDATM%avs(n),a2x,avifld,avofld,rearr)
enddo
do n = 1,SDATM%nstreams
call shr_dmodel_translateAV(SDATM%avs(n),avstrm,stifld,stofld,rearr)
enddo
call t_stopf('datm_scatter')
else
call mct_aVect_zero(a2x)
endif
call t_startf('datm_mode')
select case (trim(atm_mode))
case('COPYALL')
! do nothing extra
case('CPLHIST')
! do nothing extra
case('TN460')
if (firstcall) then
if (sprec < 1 .or. sswdn < 1 .or. sswup < 1) then
write(logunit,F00) ' ERROR: prec swdn and swup must be in streams for TN460'
call shr_sys_abort(trim(subname)//' ERROR: prec swdn and swup must be in streams for TN460')
endif
call datm_shr_TN460getFactors(factorFn,windFactor,qsatFactor, &
mpicom,compid,gsmap,ggrid,SDATM%nxg,SDATM%nyg)
endif
call shr_cal_date2eday(currentYMD,eday)
rday = mod(eday,365) + real(currentTOD)/SHR_CONST_CDAY
cosfactor = cos((2.0_R8*SHR_CONST_PI*rday)/365 - phs_c0)
lsize = mct_avect_lsize(a2x)
do n = 1,lsize
a2x%rAttr(kz,n) = 10.0_R8
!--- correction to NCEP winds based on QSCAT ---
a2x%rAttr(ku,n) = a2x%rAttr(ku,n)*windFactor(n)
a2x%rAttr(kv,n) = a2x%rAttr(kv,n)*windFactor(n)
!--- density, tbot, & pslv taken directly from input stream, set pbot ---
! a2x%rAttr(kdens,n) = <copied from input streams bundle>
! a2x%rAttr(ktbot,n) = <copied from input streams bundle>
! a2x%rAttr(kpslv,n) = <copied from input streams bundle>
a2x%rAttr(kpbot,n) = a2x%rAttr(kpslv,n)
!--- correction to NCEP Arctic & Antarctic air T & potential T ---
if ( yc(n) < -60.0_R8 ) then
tMin = (avg_c0 + avg_c1*yc(n)) + (amp_c0 + amp_c1*yc(n))*cosFactor + tKFrz
a2x%rAttr(ktbot,n) = max(a2x%rAttr(ktbot,n), tMin)
else if ( yc(n) > 60.0_R8 ) then
factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + factor * dTarc(mm)
endif
a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)
!--- correction to NCEP relative humidity for heat budget balance ---
qsatT = 640380.0_R8/ exp(5107.4_R8/a2x%rAttr(ktbot,n))
a2x%rAttr(kshum,n) = a2x%rAttr(kshum,n) - &
(qsatT)*qsatFactor(n)/(100.0_R8*a2x%rAttr(kdens,n))
!-------------------------------------------------------------------------
! PRECIPITATION DATA
!-------------------------------------------------------------------------
avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)/86400.0_R8 ! convert mm/day to kg/m^2/s
avstrm%rAttr(sprec,n) = (avstrm%rAttr(sprec,n)*1.14168_R8)+7.44923e-07_R8 ! correction to GXGXS for water budget
a2x%rAttr(krc,n) = 0.0_R8 ! default zero
a2x%rAttr(ksc,n) = 0.0_R8
if (a2x%rAttr(ktbot,n) < tKFrz ) then ! assign precip to rain/snow components
a2x%rAttr(krl,n) = 0.0_R8
a2x%rAttr(ksl,n) = avstrm%rAttr(sprec,n)
else
a2x%rAttr(krl,n) = avstrm%rAttr(sprec,n)
a2x%rAttr(ksl,n) = 0.0_R8
endif
!-------------------------------------------------------------------------
! RADIATION DATA
!-------------------------------------------------------------------------
!--- fabricate required swdn components from net swdn ---
a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdn,n)*(0.28_R8)
a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdn,n)*(0.31_R8)
a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdn,n)*(0.24_R8)
a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdn,n)*(0.17_R8)
!--- compute net short-wave based on LY08 latitudinally-varying albedo ---
if (firstcall) then
a2x%rAttr(kswnet,n) = avstrm%rAttr(sswdn,n)-avstrm%rAttr(sswup,n)
else
avg_alb = (x2a%rAttr(kavsdr,n) + x2a%rAttr(kanidr,n) + &
x2a%rAttr(kavsdf,n) + x2a%rAttr(kanidf,n)) * 0.25_R8
a2x%rAttr(kswnet,n) = avstrm%rAttr(sswdn,n)*(1.0_R8 - avg_alb)
endif
!--- corrections to GISS sswdn for heat budget balancing ---
factor = 1.0_R8
if ( -60.0_R8 < yc(n) .and. yc(n) < -50.0_R8 ) then
factor = 1.0_R8 - (yc(n) + 60.0_R8)*(0.05_R8/10.0_R8)
else if ( -50.0_R8 < yc(n) .and. yc(n) < 30.0_R8 ) then
factor = 0.95_R8
else if ( 30.0_R8 < yc(n) .and. yc(n) < 40._R8 ) then
factor = 1.0_R8 - (40.0_R8 - yc(n))*(0.05_R8/10.0_R8)
endif
a2x%rAttr(kswnet,n) = a2x%rAttr(kswnet,n)*factor
a2x%rAttr(kswvdr,n) = a2x%rAttr(kswvdr,n)*factor
a2x%rAttr(kswndr,n) = a2x%rAttr(kswndr,n)*factor
a2x%rAttr(kswvdf,n) = a2x%rAttr(kswvdf,n)*factor
a2x%rAttr(kswndf,n) = a2x%rAttr(kswndf,n)*factor
!--- correction to GISS lwdn in Arctic ---
if ( yc(n) > 60._R8 ) then
factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
a2x%rAttr(klwdn,n) = a2x%rAttr(klwdn,n) + factor * dLWarc
endif
enddo ! lsize
case('CORE2')
if (firstcall) then
if (sprec < 1 .or. sswdn < 1) then
write(logunit,F00) ' ERROR: prec and swdn must be in streams for CORE2'
call shr_sys_abort(trim(subname)//' ERROR: prec and swdn must be in streams for CORE2')
endif
call datm_shr_CORE2getFactors(factorFn,windFactor,winddFactor,qsatFactor, &
mpicom,compid,gsmap,ggrid,SDATM%nxg,SDATM%nyg)
endif
call shr_cal_date2eday(currentYMD,eday)
rday = mod(eday,365) + real(currentTOD)/SHR_CONST_CDAY
cosfactor = cos((2.0_R8*SHR_CONST_PI*rday)/365 - phs_c0)
lsize = mct_avect_lsize(a2x)
do n = 1,lsize
a2x%rAttr(kz,n) = 10.0_R8
!--- correction to NCEP winds based on QSCAT ---
uprime = a2x%rAttr(ku,n)*windFactor(n)
vprime = a2x%rAttr(kv,n)*windFactor(n)
a2x%rAttr(ku,n) = uprime*cos(winddFactor(n)*degtorad)- &
vprime*sin(winddFactor(n)*degtorad)
a2x%rAttr(kv,n) = uprime*sin(winddFactor(n)*degtorad)+ &
vprime*cos(winddFactor(n)*degtorad)
!--- density, tbot, & pslv taken directly from input stream, set pbot ---
! a2x%rAttr(kdens,n) = <copied from input streams bundle>
! a2x%rAttr(ktbot,n) = <copied from input streams bundle>
! a2x%rAttr(kpslv,n) = <copied from input streams bundle>
a2x%rAttr(kpbot,n) = a2x%rAttr(kpslv,n)
!--- correction to NCEP Arctic & Antarctic air T & potential T ---
if ( yc(n) < -60.0_R8 ) then
tMin = (avg_c0 + avg_c1*yc(n)) + (amp_c0 + amp_c1*yc(n))*cosFactor + tKFrz
a2x%rAttr(ktbot,n) = max(a2x%rAttr(ktbot,n), tMin)
else if ( yc(n) > 60.0_R8 ) then
factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + factor * dTarc(mm)
endif
a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)
!--- correction to NCEP relative humidity for heat budget balance ---
a2x%rAttr(kshum,n) = a2x%rAttr(kshum,n) + qsatFactor(n)
!-------------------------------------------------------------------------
! PRECIPITATION DATA
!-------------------------------------------------------------------------
avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)/86400.0_R8 ! convert mm/day to kg/m^2/s
! only correct satellite products, do not correct Serreze Arctic data
if ( yc(n) < 58. ) then
avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)*1.14168_R8
endif
if ( yc(n) >= 58. .and. yc(n) < 68. ) then
factor = MAX(0.0_R8, 1.0_R8 - 0.1_R8*(yc(n)-58.0_R8) )
avstrm%rAttr(sprec,n) = avstrm%rAttr(sprec,n)*(factor*(1.14168_R8 - 1.0_R8) + 1.0_R8)
endif
a2x%rAttr(krc,n) = 0.0_R8 ! default zero
a2x%rAttr(ksc,n) = 0.0_R8
if (a2x%rAttr(ktbot,n) < tKFrz ) then ! assign precip to rain/snow components
a2x%rAttr(krl,n) = 0.0_R8
a2x%rAttr(ksl,n) = avstrm%rAttr(sprec,n)
else
a2x%rAttr(krl,n) = avstrm%rAttr(sprec,n)
a2x%rAttr(ksl,n) = 0.0_R8
endif
!-------------------------------------------------------------------------
! RADIATION DATA
!-------------------------------------------------------------------------
!--- fabricate required swdn components from net swdn ---
a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdn,n)*(0.28_R8)
a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdn,n)*(0.31_R8)
a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdn,n)*(0.24_R8)
a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdn,n)*(0.17_R8)
!--- compute net short-wave based on LY08 latitudinally-varying albedo ---
avg_alb = ( 0.069 - 0.011*cos(2.0_R8*yc(n)*degtorad ) )
a2x%rAttr(kswnet,n) = avstrm%rAttr(sswdn,n)*(1.0_R8 - avg_alb)
!--- corrections to GISS sswdn for heat budget balancing ---
factor = 1.0_R8
if ( -60.0_R8 < yc(n) .and. yc(n) < -50.0_R8 ) then
factor = 1.0_R8 - (yc(n) + 60.0_R8)*(0.05_R8/10.0_R8)
else if ( -50.0_R8 < yc(n) .and. yc(n) < 30.0_R8 ) then
factor = 0.95_R8
else if ( 30.0_R8 < yc(n) .and. yc(n) < 40._R8 ) then
factor = 1.0_R8 - (40.0_R8 - yc(n))*(0.05_R8/10.0_R8)
endif
a2x%rAttr(kswnet,n) = a2x%rAttr(kswnet,n)*factor
a2x%rAttr(kswvdr,n) = a2x%rAttr(kswvdr,n)*factor
a2x%rAttr(kswndr,n) = a2x%rAttr(kswndr,n)*factor
a2x%rAttr(kswvdf,n) = a2x%rAttr(kswvdf,n)*factor
a2x%rAttr(kswndf,n) = a2x%rAttr(kswndf,n)*factor
!--- correction to GISS lwdn in Arctic ---
if ( yc(n) > 60._R8 ) then
factor = MIN(1.0_R8, 0.1_R8*(yc(n)-60.0_R8) )
a2x%rAttr(klwdn,n) = a2x%rAttr(klwdn,n) + factor * dLWarc
endif
enddo ! lsize
case('CLMNCEP')
if (firstcall) then
if (swind < 1 .or. stbot < 1) then
write(logunit,F00) ' ERROR: wind and tbot must be in streams for CLMNCEP'
call shr_sys_abort(trim(subname)//' ERROR: wind and tbot must be in streams for CLMNCEP')
endif
rtmp = maxval(a2x%rAttr(ktbot,:))
call shr_mpi_max(rtmp,tbotmax,mpicom,'datm_tbot',all=.true.)
rtmp = maxval(x2a%rAttr(kanidr,:))
call shr_mpi_max(rtmp,anidrmax,mpicom,'datm_ani',all=.true.)
if (stdew > 0) then
rtmp = maxval(avstrm%rAttr(stdew,:))
call shr_mpi_max(rtmp,tdewmax,mpicom,'datm_tdew',all=.true.)
endif
if (my_task == master_task) &
write(logunit,*) trim(subname),' max values = ',tbotmax,tdewmax,anidrmax
endif
lsize = mct_avect_lsize(a2x)
do n = 1,SDATM%nstreams
call shr_dmodel_translateAV(SDATM%avs(n),avstrm,stifld,stofld,rearr)
enddo
do n = 1,lsize
!--- bottom layer height ---
if (sz < 1) a2x%rAttr(kz,n) = 30.0_R8
!--- temperature ---
if (tbotmax < 50.0_R8) a2x%rAttr(ktbot,n) = a2x%rAttr(ktbot,n) + tkFrz
a2x%rAttr(kptem,n) = a2x%rAttr(ktbot,n)
!--- pressure ---
if (spbot < 1) a2x%rAttr(kpbot,n) = pstd
a2x%rAttr(kpslv,n) = a2x%rAttr(kpbot,n)
!--- u, v wind velocity ---
a2x%rAttr(ku,n) = avstrm%rAttr(swind,n)/sqrt(2.0_R8)
a2x%rAttr(kv,n) = a2x%rAttr(ku,n)
!--- specific humidity ---
tbot = a2x%rAttr(ktbot,n)
pbot = a2x%rAttr(kpbot,n)
if (sshum > 0) then
e = datm_shr_esat(tbot,tbot)
qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
if (qsat < a2x%rAttr(kshum,n)) then
a2x%rAttr(kshum,n) = qsat
endif
else if (srh > 0) then
e = avstrm%rAttr(srh,n) * 0.01_R8 * datm_shr_esat(tbot,tbot)
qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
a2x%rAttr(kshum,n) = qsat
else if (stdew > 0) then
if (tdewmax < 50.0_R8) avstrm%rAttr(stdew,n) = avstrm%rAttr(stdew,n) + tkFrz
e = datm_shr_esat(avstrm%rAttr(stdew,n),tbot)
qsat = (0.622_R8 * e)/(pbot - 0.378_R8 * e)
a2x%rAttr(kshum,n) = qsat
else
call shr_sys_abort(subname//'ERROR: cannot compute shum')
endif
!--- density ---
vp = (a2x%rAttr(kshum,n)*pbot) / (0.622_R8 + 0.0378_R8 * a2x%rAttr(kshum,n))
a2x%rAttr(kdens,n) = (pbot - 0.0378_R8 * vp) / (tbot*rdair)
!--- downward longwave ---
if (slwdn < 1) then
e = a2x%rAttr(kpslv,n) * a2x%rAttr(kshum,n) / (0.622_R8 + 0.378_R8 * a2x%rAttr(kshum,n))
ea = 0.70_R8 + 5.95e-05_R8 * 0.01_R8 * e * exp(1500.0_R8/tbot)
a2x%rAttr(klwdn,n) = ea * stebol * tbot**4
endif
!--- shortwave radiation ---
if (sswdndf > 0 .and. sswdndr > 0) then
a2x%rAttr(kswndr,n) = avstrm%rAttr(sswdndr,n) * 0.50_R8
a2x%rAttr(kswvdr,n) = avstrm%rAttr(sswdndr,n) * 0.50_R8
a2x%rAttr(kswndf,n) = avstrm%rAttr(sswdndf,n) * 0.50_R8
a2x%rAttr(kswvdf,n) = avstrm%rAttr(sswdndf,n) * 0.50_R8
elseif (sswdn > 0) then
! relationship between incoming NIR or VIS radiation and ratio of
! direct to diffuse radiation calculated based on one year's worth of
! hourly CAM output from CAM version cam3_5_55
swndr = avstrm%rAttr(sswdn,n) * 0.50_R8
ratio_rvrf = min(0.99_R8,max(0.29548_R8 + 0.00504_R8*swndr &
-1.4957e-05_R8*swndr**2 + 1.4881e-08_R8*swndr**3,0.01_R8))
a2x%rAttr(kswndr,n) = ratio_rvrf*swndr
swndf = avstrm%rAttr(sswdn,n) * 0.50_R8
a2x%rAttr(kswndf,n) = (1._R8 - ratio_rvrf)*swndf
swvdr = avstrm%rAttr(sswdn,n) * 0.50_R8
ratio_rvrf = min(0.99_R8,max(0.17639_R8 + 0.00380_R8*swvdr &
-9.0039e-06_R8*swvdr**2 + 8.1351e-09_R8*swvdr**3,0.01_R8))
a2x%rAttr(kswvdr,n) = ratio_rvrf*swvdr
swvdf = avstrm%rAttr(sswdn,n) * 0.50_R8
a2x%rAttr(kswvdf,n) = (1._R8 - ratio_rvrf)*swvdf
else
call shr_sys_abort(subName//'ERROR: cannot compute short-wave down')
endif
!--- swnet: a diagnostic quantity ---
if (anidrmax < 1.0e-8 .or. anidrmax > SHR_CONST_SPVAL * 0.9_R8) then
a2x%rAttr(kswnet,n) = 0.0_R8
else
a2x%rAttr(kswnet,n) = (1.0_R8-x2a%rAttr(kanidr,n))*a2x%rAttr(kswndr,n) + &
(1.0_R8-x2a%rAttr(kavsdr,n))*a2x%rAttr(kswvdr,n) + &
(1.0_R8-x2a%rAttr(kanidf,n))*a2x%rAttr(kswndf,n) + &
(1.0_R8-x2a%rAttr(kavsdf,n))*a2x%rAttr(kswvdf,n)
endif
!--- rain and snow ---
if (sprecc > 0 .and. sprecl > 0) then
a2x%rAttr(krc,n) = avstrm%rAttr(sprecc,n)
a2x%rAttr(krl,n) = avstrm%rAttr(sprecl,n)
elseif (sprecn > 0) then
a2x%rAttr(krc,n) = avstrm%rAttr(sprecn,n)*0.1_R8
a2x%rAttr(krl,n) = avstrm%rAttr(sprecn,n)*0.9_R8
else
call shr_sys_abort(subName//'ERROR: cannot compute rain and snow')
endif
!--- split precip between rain & snow ---
!--- note: aribitrarily small negative values cause CLM to crash ---
frac = (tbot - tkFrz)*0.5_R8 ! ramp near freezing
frac = min(1.0_R8,max(0.0_R8,frac)) ! bound in [0,1]
a2x%rAttr(ksc,n) = max(0.0_R8, a2x%rAttr(krc,n)*(1.0_R8 - frac) )
a2x%rAttr(ksl,n) = max(0.0_R8, a2x%rAttr(krl,n)*(1.0_R8 - frac) )
a2x%rAttr(krc,n) = max(0.0_R8, a2x%rAttr(krc,n)*( frac) )
a2x%rAttr(krl,n) = max(0.0_R8, a2x%rAttr(krl,n)*( frac) )
enddo
end select
call t_stopf('datm_mode')
if (write_restart) then
call t_startf('datm_restart')
call seq_infodata_GetData( infodata, case_name=case_name)
write(rest_file,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
trim(case_name), '.datm.r.', yy,'-',mm,'-',dd,'-',currentTOD,'.nc'
write(rest_file_strm,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
trim(case_name), '.datm.rs1.', yy,'-',mm,'-',dd,'-',currentTOD,'.bin'
if (my_task == master_task) then
nu = shr_file_getUnit()
open(nu,file=trim(rpfile),form='formatted')
write(nu,'(a)') rest_file
write(nu,'(a)') rest_file_strm
close(nu)
call shr_file_freeUnit(nu)
endif
! if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file),currentYMD,currentTOD
! call shr_pcdf_readwrite('write',trim(rest_file),mpicom,gsmap,clobber=.true., &
! rf1=somtp,rf1n='somtp')
if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file_strm),currentYMD,currentTOD
call shr_strdata_restWrite(trim(rest_file_strm),SDATM,mpicom,trim(case_name),'SDATM strdata')
call shr_sys_flush(logunit)
call t_stopf('datm_restart')
endif
call t_stopf('datm')
!----------------------------------------------------------------------------
! Log output for model date
! Reset shr logging to original values
!----------------------------------------------------------------------------
call t_startf('datm_run2')
if (my_task == master_task) then
write(logunit,F04) trim(myModelName),': model date ', CurrentYMD,CurrentTOD
call shr_sys_flush(logunit)
end if
firstcall = .false.
call shr_file_setLogUnit (shrlogunit)
call shr_file_setLogLevel(shrloglev)
call shr_sys_flush(logunit)
call t_stopf('datm_run2')
call t_stopf('DATM_RUN')
end subroutine datm_comp_run
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: datm_comp_final
!
! !DESCRIPTION:
! finalize method for dead atm model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
!
subroutine datm_comp_final() 1
implicit none
!EOP
!--- formats ---
character(*), parameter :: F00 = "('(datm_comp_final) ',8a)"
character(*), parameter :: F91 = "('(datm_comp_final) ',73('-'))"
character(*), parameter :: subName = "(datm_comp_final) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
call t_startf('DATM_FINAL')
if (my_task == master_task) then
write(logunit,F91)
write(logunit,F00) trim(myModelName),': end of main integration loop'
write(logunit,F91)
end if
call t_stopf('DATM_FINAL')
end subroutine datm_comp_final
!===============================================================================
!===============================================================================
end module datm_comp_mod