!!!!! ========================================================== !!!!!
!!!!! 'module_radiation_driver' descriptions !!!!!
!!!!! ========================================================== !!!!!
! !
! this is the radiation driver module. it prepares atmospheric !
! profiles and invokes main radiation calculations. !
! !
! in module 'module_radiation_driver' there are twe externally !
! callable subroutine: !
! !
! 'radinit' -- initialization routine !
! input: !
! ( si, NLAY, iflip, idate, jdate, ICTM, ISOL, ICO2, !
! IAER, IALB, IEMS, ICWP, NP3D, isubcsw, isubclw, !
! iovrsw, iovrlw, me ) !
! output: !
! ( none ) !
! !
! 'grrad' -- setup and invoke main radiation calls !
! input: !
! ( prsi,prsl,prslk,tgrs,qgrs,oz,vvl,slmsk, !
! xlon,xlat,tsfc,snowd,sncovr,snoalb,zorl,hprim, !
! alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc, !
! solcon,coszen,coszdg,k1oz,k2oz,facoz, !
! cv,cvt,cvb,iovrsw,iovrlw,fcice,frain,rrime,flgmin, !
! icsdsw,icsdlw, np3d,ntcw,ncld,ntoz, NTRAC,NFXR, !
! dtlw,dtsw, lsswr,lslwr,lssav,sashal,norad_precip, !
! crick_proof, ccnorm, !
! IX, IM, LM, iflip, me, lprnt, ipt, kdt, !
! output: !
! htrsw,topfsw,sfcfsw,sfalb, !
! htrlw,topflw,sfcflw,tsflw,semis,cldcov,cldsa !
! input/output: !
! fluxr !
! optional output: !
! HTRSWB,HTRLWB) !
! !
! !
! external modules referenced: !
! 'module machine' in 'machine.f' !
! 'module funcphys' in 'funcphys.f' !
! 'module physcons' in 'physcons.f !
! !
! 'module module_radiation_gases' in 'radiation_gases.f' !
! 'module module_radiation_aerosols' in 'radiation_aerosols.f' !
! 'module module_radiation_surface' in 'radiation_surface.f' !
! 'module module_radiation_clouds' in 'radiation_clouds.f' !
! !
! 'module module_radsw_cntr_para' in 'radsw_xxxx_param.f' !
! 'module module_radsw_parameters' in 'radsw_xxxx_param.f' !
! 'module module_radsw_main' in 'radsw_xxxx_main.f' !
! !
! 'module module_radlw_cntr_para' in 'radlw_xxxx_param.f' !
! 'module module_radlw_parameters' in 'radlw_xxxx_param.f' !
! 'module module_radlw_main' in 'radlw_xxxx_main.f' !
! !
! where xxxx may vary according to different scheme selection !
! !
! !
! program history log: !
! mm-dd-yy ncep - created program grrad !
! 08-12-03 yu-tai hou - re-written for modulized radiations !
! 11-06-03 yu-tai hou - modified !
! 01-18-05 s. moorthi - NOAH/ICE model changes added !
! 05-10-05 yu-tai hou - modified module structure !
! 12-xx-05 s. moorthi - sfc lw flux adj by mean temperature !
! 02-20-06 yu-tai hou - add time variation for co2 data, and !
! solar const. add sfc emiss change !
! 03-21-06 s. Moorthi - added surface temp over ice !
! 07-28-06 yu-tai hou - add stratospheric vocanic aerosols !
! 03-14-07 yu-tai hou - add generalized spectral band interp !
! for aerosol optical prop. (sw and lw) !
! 04-10-07 yu-tai hou - spectral band sw/lw heating rates !
! 05-04-07 yu-tai hou - make options for clim based and modis !
! based (h. wei and c. marshall) albedo !
! 09-05-08 yu-tai hou - add the initial date and time 'idate' !
! and control param 'ICTM' to the passing param list!
! to handel different time/date requirements for !
! external data (co2, aeros, solcon, ...) !
! 10-10-08 yu-tai hou - add the ICTM=-2 option for combining !
! initial condition data with seasonal cycle from !
! climatology. !
! 03-12-09 yu-tai hou - use two time stamps to keep tracking !
! dates for init cond and fcst time. remove volcanic!
! aerosols data in climate hindcast (ICTM=-2). !
! 03-16-09 yu-tai hou - included sub-column clouds approx. !
! control flags isubcsw/isubclw in initialization !
! subroutine. passed auxiliary cloud control arrays !
! icsdsw/icsdlw (if isubcsw/isubclw =2, it will be !
! the user provided permutation seeds) to the sw/lw !
! radiation calculation programs. also moved cloud !
! overlapping control flags iovrsw/iovrlw from main !
! radiation routines to the initialization routines.!
! 04-02-09 yu-tai hou - modified surface control flag iems to !
! have additional function of if the surface-air !
! interface have the same or different temperature !
! for radiation calculations. !
! 04-03-09 yu-tai hou - modified to add lw surface emissivity !
! as output variable. changed the sign of sfcnsw to !
! be positive value denote solar flux goes into the !
! ground (this is needed to reduce sign confusion !
! in other part of model) !
! 04-20-09 carlos perez - prepare driver for nmmb. added option!
! of run the gfs's radiation on nmmb !
! 09-09-09 fanglin yang (thru s.moorthi) added QME5 QME6 to E-20!
! 01-09-10 sarah lu - added gocart option, revised grrad for!
! gocart coupling. calling argument modifed: ldiag3 !
! removed; cldcov/fluxr sequence changed; cldcov is !
! changed from accumulative to instant field and !
! from input/output to output field !
! 01-24-10 sarah lu - added aod to fluxr, added prslk and !
! oz to setaer input argument (for gocart coupling),!
! added tau_gocart to setaer output argument (for, !
! aerosol diag) !
! 07-08-10 s.moorthi - updated the NEMS version for new physics !
! 07-28-10 yu-tai hou - changed grrad interface to allow all !
! components of sw/lw toa/sfc instantaneous values !
! being passed to the calling program. moved the !
! computaion of sfc net sw flux (sfcnsw) to the !
! calling program. merged carlos' nmmb modification.!
! 07-30-10 s. moorthi - corrected some errors associated with !
! unit changes !
! 12-02-10 s. moorthi/y. hou - removed the use of aerosol flags !
! 'iaersw' 'iaerlw' from radiations and replaced !
! them by using the runtime variable iaerflg and !
! laswflg defined in module radiation_aerosols. !
! also replaced param nspc in grrad with the use of !
! max_num_gridcomp in module radiation_aerosols. !
! 01-03-11 y. hou - added sea/land madk 'slmsk' to the !
! argument list of subrotine setaer call for the !
! newly modified horizontal bi-linear interpolation !
! in climatological aerosols schem. !
! 09-30-11 H.M. Lin - rename "grrad" to "grrad_nmmb" for using!
! in regional to avoid the confliction with global !
! when compile the nems. Also change names of the !
! related subroutines to "astronomy_nmmb" & !
! "solinit_nmmb" !
! 03-15-12 H.M. Lin - changed grrad interface to allow 'cldsa'!
! (fraction of clouds for low, mid, hi, tot, bl) !
! passed to the calling program. !
! 09-25-12 y. hou - added optional extra top layer 'LTP' on !
! top of low ceiling models !
! (** to cover the gases atop, expecially O3) !
! !
!!!!! ========================================================== !!!!!
!!!!! end descriptions !!!!!
!!!!! ========================================================== !!!!!
!=============================================!
module module_radiation_driver_nmmb !
!.............................................!
!
use physparam
use physcons, only : con_eps, con_epsm1, &
& con_fvirt, PI=>con_pi
use funcphys, only : fpvs, fpkapx
use MODULE_CONSTANTS, only : EPSQ
use module_radiation_astronomy_nmmb,only : sol_init_nmmb, &
& sol_update_nmmb, coszmn_nmmb
use module_radiation_gases_nmmb,only : NF_VGAS, getgases, getozn, &
& gas_init, gas_update
use module_radiation_aerosols_nmmb,only : NF_AESW, NF_AELW,setaer,&
& aer_init, aer_update
! &, NSPC1 ! optn for aod output
use module_radiation_surface_nmmb, only : NF_ALBD,sfc_init,setalb,&
& setemis
use module_radiation_clouds_nmmb, only : NF_CLDS, cld_init, &
& progcld1, progcld2, diagcld1 &
& , progcld8
use module_radsw_parameters, only : topfsw_type, sfcfsw_type, &
& profsw_type,cmpfsw_type,NBDSW
use module_radsw_main_nmmb, only : rswinit, swrad
use module_radlw_parameters, only : topflw_type, sfcflw_type, &
& proflw_type, NBDLW
use module_radlw_main_nmmb, only : rlwinit, lwrad
!
implicit none
!
private
! --- version tag and last revision date
character(40), parameter :: &
& VTAGRAD='NCEP-Radiation_driver v5.2 Jan 2013 '
! & VTAGRAD='NCEP-Radiation_driver v5.1 Nov 2012 '
! & VTAGRAD='NCEP-Radiation_driver v5.0 Aug 2012 '
! --- constant values
real (kind=kind_phys) :: QMIN, QME5, QME6
! parameter (QMIN=1.0e-10, QME5=1.0e-5, QME6=1.0e-6, EPSQ=1.0e-12)
parameter (QMIN=1.0e-10, QME5=1.0e-7, QME6=1.0e-7)
! parameter (QMIN=1.0e-10, QME5=1.0e-20, QME6=1.0e-20, EPSQ=1.0e-12)
real, parameter :: prsmin = 1.0e-6 ! toa pressure minimum value in mb (hpa)
! --- control flags set in subr radinit:
integer :: itsfc =0 ! flag for lw sfc air/ground interface temp setting
! --- data input control variables set in subr radupdate:
integer :: month0=0, iyear0=0, monthd=0
logical :: loz1st =.true. ! first-time clim ozone data read flag
! --- optional extra top layer on top of low ceiling models (Y. Hou, 2012-09-25)
! integer, parameter :: LTP = 0 ! do no add an extra top layer
integer, parameter :: LTP = 1 ! add an extra top layer
logical :: lextop = (LTP > 0)
!
integer, parameter :: CHK= CHNK_RRTM ! CHNK_RRTM is macro defined in configure.nems
! --- publicly accessible module programs:
public radinit_nmmb, radupdate_nmmb, grrad_nmmb, dayparts
! =================
contains
! =================
!-----------------------------------
subroutine radinit_nmmb &
!...................................
! --- inputs:
& ( si, NLAY, me )
! --- outputs:
! ( none )
! ================= subprogram documentation block ================ !
! !
! subprogram: radinit initialization of radiation calculations !
! !
! !
! program history log: !
! 08-14-2003 yu-tai hou created !
! !
! usage: call radinit !
! !
! attributes: !
! language: fortran 90 !
! machine: ibm sp !
! !
! ==================== defination of variables ==================== !
! !
! input parameters: !
! si : model vertical sigma interface !
! NLAY : number of model vertical layers !
! iflip : control flag for direction of vertical index !
! =0: index from toa to surface !
! =1: index from surface to toa !
! idate(8) : ncep absolute date and time of initial condition !
! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
! jdate(8) : ncep absolute date and time at fcst time !
! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
! ICTM :=yyyy#, external data time/date control flag !
! = -2: same as 0, but superimpose seasonal cycle !
! from climatology data set. !
! = -1: use user provided external data for the !
! forecast time, no extrapolation. !
! = 0: use data at initial cond time, if not !
! available, use latest, no extrapolation. !
! = 1: use data at the forecast time, if not !
! available, use latest and extrapolation. !
! =yyyy0: use yyyy data for the forecast time, !
! no further data extrapolation. !
! =yyyy1: use yyyy data for the fcst. if needed, do !
! extrapolation to match the fcst time. !
! ISOL :=0: use a fixed solar constant value !
! =1: use 11-year cycle solar constant table !
! ICO2 :=0: use prescribed global mean co2 (old oper) !
! =1: use observed co2 annual mean value only !
! =2: use obs co2 monthly data with 2-d variation !
! IAER : 3-digit aerosol flag (for volc, lw, sw) !
! = 0: turn all aeros effects off (sw,lw,volc) !
! = 1: use clim tropspheric aerosol for sw only !
! = 10: use clim tropspheric aerosol for lw only !
! = 11: use clim tropspheric aerosol for both sw/lw !
! =100: volc aerosol only for both sw and lw !
! =101: volc and clim trops aerosol for sw only !
! =110: volc and clim trops aerosol for lw only !
! =111: volc and clim trops aerosol for both sw/lw !
! = 2: gocart prognostic, without volc forcing !
! =12: gocart prognostic, with volcanic forcing !
! IALB : control flag for surface albedo schemes !
! =0: climatology, based on surface veg types !
! =1: modis retrieval based surface albedo scheme !
! =2: use externally provided albedoes directly. !
! IEMS : ab 2-digit control flag !
! a =0 set sfc air/ground t same for lw radiation !
! =1 set sfc air/ground t diff for lw radiation !
! b =0 use fixed sfc emissivity=1.0 (black-body) !
! =1 use varying climtology sfc emiss (veg based)!
! =2 future development (not yet) !
! ICWP : control flag for cloud generation schemes !
! =-1: use diagnostic cloud scheme (GFDL type !
! =0 : use diagnostic cloud scheme !
! =1 : use prognostic cloud scheme (default) !
! NP3D :=3: ferrier's microphysics cloud scheme !
! =4: zhao/carr/sundqvist microphysics cloud !
! =5: nmmb ferrier+bmj microphysics scheme !
! =8: Thompson microphysics scheme !
! isubcsw/isubclw : sub-column cloud approx control flag (sw/lw rad) !
! =0: with out sub-column cloud approximation !
! =1: mcica sub-col approx. prescribed random seed !
! =2: mcica sub-col approx. provided random seed !
! iovrsw/iovrlw : control flag for cloud overlap (sw/lw rad) !
! =0: random overlapping clouds !
! =1: max/ran overlapping clouds !
! me : print control flag !
! !
! outputs: (none) !
! !
! usage: call radinit !
! !
! subroutines called: cldinit, aerinit, rlwinit, rswinit, gasinit !
! !
! =================================================================== !
!
implicit none
! --- inputs:
integer, intent(in) :: NLAY, me
real (kind=kind_phys), intent(in) :: si(:)
! --- outputs: (none, to module variables)
! --- locals:
integer :: RICWP, icldflg_org
!
!===> ... begin here
!
! --- ... set up control variables
itsfc = iemsflg / 10 ! sfc air/ground temp control
loz1st = (ioznflg == 0) ! first-time clim ozone data read flag
if (me == 0) then
! print *,' NEW RADIATION PROGRAM STRUCTURES -- SEP 01 2004'
print *,' NEW RADIATION PROGRAM STRUCTURES BECAME OPER. ', &
& ' May 01 2007'
print *, VTAGRAD !print out version tag
print *,' - Selected Control Flag settings: ICTMflg=',ictmflg, &
& ' ISOLar =',isolar, ' ICO2flg=',ico2flg,' IAERflg=',iaerflg, &
& ' IALBflg=',ialbflg,' IEMSflg=',iemsflg,' ICLDflg=',icldflg, &
& ' ICMPHYS=',icmphys,' IOZNflg=',ioznflg
print *,' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw, &
& ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw
print *,' LSASHAL=',lsashal,' LCRICK=',lcrick,' LCNORM=',lcnorm,&
& ' LNOPREC=',lnoprec
print *,' LTP =',LTP,', add extra top layer =',lextop
if ( ictmflg==0 .or. ictmflg==-2 ) then
print *,' Data usage is limited by initial condition!'
print *,' No volcanic aerosols'
endif
if ( isubclw == 0 ) then
print *,' - ISUBCLW=',isubclw,' No McICA, use grid ', &
& 'averaged cloud in LW radiation'
elseif ( isubclw == 1 ) then
print *,' - ISUBCLW=',isubclw,' Use McICA with fixed', &
& 'permutation seeds for LW random number generator'
elseif ( isubclw == 2 ) then
print *,' - ISUBCLW=',isubclw,' Use McICA with random ', &
& 'permutation seeds for LW random number generator'
else
print *,' - ERROR!!! ISUBCLW=',isubclw,' is not a ', &
& 'valid option '
stop
endif
if ( isubcsw == 0 ) then
print *,' - ISUBCSW=',isubcsw,' No McICA, use grid ', &
& 'averaged cloud in SW radiation'
elseif ( isubcsw == 1 ) then
print *,' - ISUBCSW=',isubcsw,' Use McICA with fixed', &
& 'permutation seeds for SW random number generator'
elseif ( isubcsw == 2 ) then
print *,' - ISUBCSW=',isubcsw,' Use McICA with random ', &
& 'permutation seeds for SW random number generator'
else
print *,' - ERROR!!! ISUBCSW=',isubcsw,' is not a ', &
& 'valid option '
stop
endif
if ( isubcsw /= isubclw ) then
print *,' - *** Notice *** ISUBCSW /= ISUBCLW !!!', &
& isubcsw, isubclw
endif
endif
! --- ... call astronomy initialization routine
call sol_init_nmmb ( me )
! --- ... call aerosols initialization routine
call aer_init ( NLAY, me )
! --- ... call co2 and other gases initialization routine
call gas_init ( me )
! --- ... call surface initialization routine
call sfc_init ( me )
! --- ... call cloud initialization routine
call cld_init ( si, NLAY, me)
!==========================================================
! The following use for NP3D=5 GFDL cloud in radiation
!==========================================================
if ( icmphys==5 ) then
ricwp = 0
icldflg_org = icldflg
icldflg = ricwp
endif
! --- ... call lw radiation initialization routine
call rlwinit ( me )
! --- ... call sw radiation initialization routine
call rswinit ( me )
!
!==========================================================
! The following use for NP3D=5 GFDL cloud in radiation
! return icldflg to initial setting
!==========================================================
if ( icmphys==5 ) then
icldflg = icldflg_org
endif
return
!...................................
end subroutine radinit_nmmb
!-----------------------------------
!-----------------------------------
subroutine radupdate_nmmb &
!...................................
! --- inputs:
& ( idate, jdate, deltsw, deltim, lsswr, me, &
! --- outputs:
& slag,sdec,cdec,solcon &
& )
! ================= subprogram documentation block ================ !
! !
! subprogram: radupdate calls many update subroutines to check and !
! update radiation required but time varying data sets and module !
! variables. !
! !
! usage: call radupdate !
! !
! attributes: !
! language: fortran 90 !
! machine: ibm sp !
! !
! ==================== defination of variables ==================== !
! !
! input parameters: !
! idate(8) : ncep absolute date and time of initial condition !
! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
! jdate(8) : ncep absolute date and time at fcst time !
! (yr, mon, day, t-zone, hr, min, sec, mil-sec) !
! deltsw : sw radiation calling frequency in seconds !
! deltim : model timestep in seconds !
! lsswr : logical flags for sw radiation calculations !
! me : print control flag !
! !
! outputs: !
! slag : equation of time in radians !
! sdec, cdec : sin and cos of the solar declination angle !
! solcon : sun-earth distance adjusted solar constant (w/m2) !
! !
! module variables: !
! isolar : solar constant cntrl (see ISOL description) !
! ictmfl : flag for initial condition gh-gas data source !
! loz1st : first-time clim ozone data read flag !
! !
! subroutines called: sol_update, aer_update, gas_update !
! !
! =================================================================== !
!
implicit none
! --- inputs:
integer, intent(in) :: idate(:), jdate(:), me
logical, intent(in) :: lsswr
real (kind=kind_phys), intent(in) :: deltsw, deltim
! --- outputs:
real (kind=kind_phys), intent(out) :: slag, sdec, cdec, solcon
! --- locals:
integer :: iyear, imon, iday, ihour
integer :: kyear, kmon, kday, khour
logical :: lmon_chg ! month change flag
logical :: lco2_chg ! cntrl flag for updating co2 data
logical :: lsol_chg ! cntrl flag for updating solar constant
!
!===> ... begin here
!
! --- ... time stamp at fcst time
iyear = jdate(1)
imon = jdate(2)
iday = jdate(3)
ihour = jdate(5)
! --- ... set up time stamp used for green house gases (** currently co2 only)
if ( ictmflg==0 .or. ictmflg==-2 ) then ! get external data at initial condition time
kyear = idate(1)
kmon = idate(2)
kday = idate(3)
khour = idate(5)
else ! get external data at fcst or specified time
kyear = iyear
kmon = imon
kday = iday
khour = ihour
endif ! end if_ictmflg_block
if ( month0 /= imon ) then
lmon_chg = .true.
month0 = imon
else
lmon_chg = .false.
endif
! --- ... call astronomy update routine, yearly update, no time interpolation
if (lsswr) then
if ( isolar == 0 .or. isolar == 10 ) then
lsol_chg = .false.
elseif ( iyear0/=iyear ) then
lsol_chg = .true.
else
lsol_chg = ( isolar==4 .and. lmon_chg )
endif
iyear0 = iyear
call sol_update_nmmb &
! --- input:
& ( jdate,kyear,deltsw,deltim,lsol_chg, me, &
! --- outputs:
& slag,sdec,cdec,solcon &
& )
endif ! end_if_lsswr_block
! --- ... call aerosols update routine, monthly update, no time interpolation
if ( lmon_chg ) then
call aer_update ( iyear, imon, me )
endif
! --- ... call co2 and other gases update routine
if ( monthd /= kmon ) then
monthd = kmon
lco2_chg = .true.
else
lco2_chg = .false.
endif
call gas_update ( kyear,kmon,kday,khour,loz1st,lco2_chg, me )
if ( loz1st ) loz1st = .false.
! --- ... call surface update routine (currently not needed)
! call sfc_update ( iyear, imon, me )
! --- ... call clouds update routine (currently not needed)
! call cld_update ( iyear, imon, me )
!
return
!...................................
end subroutine radupdate_nmmb
!-----------------------------------
!-----------------------------------
subroutine dayparts &
!...................................
! --- inputs:
& ( XLON,SINLAT,COSLAT,SOLHR,MYPE,dtswav,nrads, &
& veclen , &
! --- outputs:
& dp_start,dp_len,dp_day,ndayparts &
& )
implicit none
integer, intent(in) :: veclen,mype,nrads
real (kind=kind_phys), intent(in) :: dtswav, solhr
real (kind=kind_phys), dimension(veclen), intent(in) :: &
& sinlat, coslat, xlon
integer, intent(out), dimension(veclen) :: dp_start,dp_len
logical, intent(out), dimension(veclen) :: dp_day
integer, intent(out) :: ndayparts
! local
logical isday
integer i
real (kind=kind_phys), dimension(veclen) :: coszen_loc, coszdg
call coszmn_nmmb &
! --- inputs:
& ( xlon,sinlat,coslat,solhr,veclen,mype, &
& dtswav,nrads , & ! Extra input
! --- outputs:
& coszen_loc, coszdg & !jm
& )
if ( veclen < 1 ) then
ndayparts = 0
else
ndayparts = 1
dp_day(ndayparts) = (coszen_loc(1) >= 0.0001)
dp_start(ndayparts) = 1
dp_len(ndayparts) = 1
do i = 2, veclen
isday = (coszen_loc(i) >= 0.0001)
if ( dp_day(ndayparts) .NEQV. isday ) THEN
ndayparts=ndayparts+1
dp_start(ndayparts) = i
dp_day(ndayparts) = isday
dp_len(ndayparts) = 1
else
dp_len(ndayparts) = dp_len(ndayparts) + 1
endif
enddo
endif
!...................................
end subroutine dayparts
!-----------------------------------
!-----------------------------------
subroutine grrad_nmmb &
!...................................
! --- inputs:
& ( prsi,prsl,prslk,tgrs,qgrs,tracer,vvl,slmsk, &
& xlon,xlat,tsfc,snowd,sncovr,snoalb,zorl,hprim, &
! & alvsf,alnsf,alvwf,alnwf,facsf,facwf, & ! processed inside
& SALBEDO,SM, & ! input for albedo cal.
& fice,tisfc, &
& sinlat,coslat,solhr,jdate,solcon, &
& dtswav,nrads, & ! extra input
& cv,cvt,cvb,fcice,frain,rrime,flgmin, &
& icsdsw,icsdlw, ntcw,ncld,ntoz, NTRAC,NFXR, &
& CPATHFAC4LW, & ! enhance factor of cloud depth for LW
& dtlw,dtsw, lsswr,lslwr,lssav, &
& its,jts, IX, IM, LM, isday, me, lprnt, ipt, kdt, &
! --- additional inputs: ! GFDL type
& tauclouds,cldf, & ! GFDL type
& cld_fraction, & ! Thompson cloud fraction
! --- outputs:
& htrsw,topfsw,sfcfsw,sfalb,coszen,coszdg, &
& htrlw,topflw,sfcflw,tsflw,semis,cldcov,cldsa, &
! --- input/output:
& fluxr &
!! --- optional outputs:
&, HTRSWB,HTRLWB &
& )
! ================= subprogram documentation block ================ !
! !
! this program is the driver of radiation calculation subroutines. * !
! It sets up profile variables for radiation input, including * !
! clouds, surface albedos, atmospheric aerosols, ozone, etc. * !
! * !
! * !
! usage: call grrad * !
! * !
! subprograms called: * !
! setalb, setemis, setaer, getozn, getgases, * !
! progcld1, progcld2, progcld8, diagcds, * !
! swrad, lwrad, fpvs * !
! * !
! attributes: * !
! language: fortran 90 * !
! machine: ibm-sp, sgi * !
! * !
! * !
! ==================== defination of variables ==================== !
! !
! input variables: !
! prsi (IX,LM+1) : model level pressure in Pa !
! prsl (IX,LM) : model layer mean pressure in Pa !
! prslk (IX,LM) : Exner function !
! tgrs (IX,LM) : model layer mean temperature in k !
! qgrs (IX,LM) : layer specific humidity in gm/gm !
! oz (IX,LM,NTRAC):layer ozone mass mixing ratio !
! vvl (IX,LM) : layer mean vertical velocity in Pa/sec !
! slmsk (IM) : sea/land mask array (sea:0,land:1,sea-ice:2) !
! xlon,xlat (IM) : grid longitude/latitude in radians !
! tsfc (IM) : surface temperature in k !
! snowd (IM) : snow depth water equivalent in mm !
! sncovr(IM) : snow cover in fraction !
! snoalb(IM) : maximum snow albedo in fraction !
! zorl (IM) : surface roughness in cm !
! hprim (IM) : topographic standard deviation in m !
! alvsf (IM) : mean vis albedo with strong cosz dependency !
! alnsf (IM) : mean nir albedo with strong cosz dependency !
! alvwf (IM) : mean vis albedo with weak cosz dependency !
! alnwf (IM) : mean nir albedo with weak cosz dependency !
! facsf (IM) : fractional coverage with strong cosz dependen !
! facwf (IM) : fractional coverage with weak cosz dependency !
! fice (IM) : ice fraction over open water grid !
! tisfc (IM) : surface temperature over ice fraction !
! solcon : solar constant (sun-earth distant adjusted) !
! coszen(IM) : mean cos of zenith angle over rad call period !
! coszdg(IM) : daytime mean cosz over rad call period !
! k1oz,k2oz,facoz : parameters for climatological ozone !
! cv (IM) : fraction of convective cloud !
! cvt, cvb (IM) : convective cloud top/bottom pressure in cb !
! iovrsw/iovrlw : control flag for cloud overlap (sw/lw rad) !
! =0 random overlapping clouds !
! =1 max/ran overlapping clouds !
! fcice : fraction of cloud ice (in ferrier scheme) !
! frain : fraction of rain water (in ferrier scheme) !
! rrime : mass ratio of total to unrimed ice ( >= 1 ) !
! flgmin : minimim large ice fraction !
! icsdsw/icsdlw : auxiliary cloud control arrays passed to main !
! (IM) radiations. if isubcsw/isubclw (input to init) !
! are set to 2, the arrays contains provided !
! random seeds for sub-column clouds generators !
! np3d : =3 brad ferrier microphysics scheme !
! =4 zhao/carr/sundqvist microphysics scheme !
! =5 external microphysics scheme provided bulk/ !
! grey quantities of cloud fields. (nmmb vars)!
! =8 Thompson microphysics scheme !
! ntcw : =0 no cloud condensate calculated !
! >0 array index location for cloud condensate !
! ncld : only used when ntcw .gt. 0 !
! ntoz : =0 climatological ozone profile !
! >0 interactive ozone profile !
! NTRAC : dimension veriable for array oz !
! NFXR : second dimension of input/output array fluxr !
! dtlw, dtsw : time duration for lw/sw radiation call in sec !
! lsswr, lslwr : logical flags for sw/lw radiation calls !
! lssav : logical flag for store 3-d cloud field !
! sashal : logical flag for Jongil's shallow convection !
! norad_precip : logical flag for not using precip in radiation !
! crick_proof : logical flag for eliminating CRICK !
! ccnorm : logical flag for incloud condensate mixing ratio!
! IX,IM : horizontal dimention and num of used points !
! LM : vertical layer dimension !
! iflip : control flag for in/out vertical indexing !
! =0 index from toa to surface !
! =1 index from surface to toa !
! me : control flag for parallel process !
! lprnt : control flag for diagnostic print out !
! ipt : index for diagnostic printout point !
! kdt : time-step number !
! !
! output variables: !
! htrsw (IX,LM) : total sky sw heating rate in k/sec !
! topfsw(IM) : sw radiation fluxes at toa, components: !
! (check module_radsw_parameters for definition) !
! %upfxc - total sky upward sw flux at toa (w/m**2) !
! %dnflx - total sky downward sw flux at toa (w/m**2) !
! %upfx0 - clear sky upward sw flux at toa (w/m**2) !
! sfcfsw(IM) : sw radiation fluxes at sfc, components: !
! (check module_radsw_parameters for definition) !
! %upfxc - total sky upward sw flux at sfc (w/m**2) !
! %dnfxc - total sky downward sw flux at sfc (w/m**2) !
! %upfx0 - clear sky upward sw flux at sfc (w/m**2) !
! %dnfx0 - clear sky downward sw flux at sfc (w/m**2) !
! sfalb (IM) : mean surface diffused sw albedo !
! cldcov(IX,LM) : 3-d cloud fraction !
! cldsa(IX,5) : fraction of clouds for low, mid, hi, tot, bl !
! htrlw (IX,LM) : total sky lw heating rate in k/sec !
! topflw(IM) : lw radiation fluxes at top, component: !
! (check module_radlw_paramters for definition) !
! %upfxc - total sky upward lw flux at toa (w/m**2) !
! %upfx0 - clear sky upward lw flux at toa (w/m**2) !
! sfcflw(IM) : lw radiation fluxes at sfc, component: !
! (check module_radlw_paramters for definition) !
! %upfxc - total sky upward lw flux at sfc (w/m**2) !
! %upfx0 - clear sky upward lw flux at sfc (w/m**2) !
! %dnfxc - total sky downward lw flux at sfc (w/m**2) !
! %dnfx0 - clear sky downward lw flux at sfc (w/m**2) !
! semis (IM) : surface lw emissivity in fraction !
! tsflw (IM) : surface air temp during lw calculation in k !
! !
! input and output variables: !
! fluxr (IX,NFXR) : to save time accumulated 2-d fields defined as:!
! 1 - toa total sky upwd lw radiation flux !
! 2 - toa total sky upwd sw radiation flux !
! 3 - sfc total sky upwd sw radiation flux !
! 4 - sfc total sky dnwd sw radiation flux !
! 5 - high domain cloud fraction !
! 6 - mid domain cloud fraction !
! 7 - low domain cloud fraction !
! 8 - high domain mean cloud top pressure !
! 9 - mid domain mean cloud top pressure !
! 10 - low domain mean cloud top pressure !
! 11 - high domain mean cloud base pressure !
! 12 - mid domain mean cloud base pressure !
! 13 - low domain mean cloud base pressure !
! 14 - high domain mean cloud top temperature !
! 15 - mid domain mean cloud top temperature !
! 16 - low domain mean cloud top temperature !
! 17 - total cloud fraction !
! 18 - boundary layer domain cloud fraction !
! 19 - sfc total sky dnwd lw radiation flux !
! 20 - sfc total sky upwd lw radiation flux !
! 21 - sfc total sky dnwd sw uv-b radiation flux !
! 22 - sfc clear sky dnwd sw uv-b radiation flux !
! 23 - toa incoming solar radiation flux !
! 24 - sfc vis beam dnwd sw radiation flux !
! 25 - sfc vis diff dnwd sw radiation flux !
! 26 - sfc nir beam dnwd sw radiation flux !
! 27 - sfc nir diff dnwd sw radiation flux !
! 28 - toa clear sky upwd lw radiation flux !
! 29 - toa clear sky upwd sw radiation flux !
! 30 - sfc clear sky dnwd lw radiation flux !
! 31 - sfc clear sky upwd sw radiation flux !
! 32 - sfc clear sky dnwd sw radiation flux !
! 33 - sfc clear sky upwd lw radiation flux !
!optional 34 - aeros opt depth at 550nm (all components) !
! .... - optional for test and future use !
! !
! optional output variables: !
! htrswb(IX,LM,NBDSW) : spectral band total sky sw heating rate !
! htrlwb(IX,LM,NBDLW) : spectral band total sky lw heating rate !
! !
! !
! definitions of internal variable arrays: !
! !
! 1. fixed gases: (defined in 'module_radiation_gases') !
! gasvmr(:,:,1) - co2 volume mixing ratio !
! gasvmr(:,:,2) - n2o volume mixing ratio !
! gasvmr(:,:,3) - ch4 volume mixing ratio !
! gasvmr(:,:,4) - o2 volume mixing ratio !
! gasvmr(:,:,5) - co volume mixing ratio !
! gasvmr(:,:,6) - cf11 volume mixing ratio !
! gasvmr(:,:,7) - cf12 volume mixing ratio !
! gasvmr(:,:,8) - cf22 volume mixing ratio !
! gasvmr(:,:,9) - ccl4 volume mixing ratio !
! !
! 2. cloud profiles: (defined in 'module_radiation_clouds') !
! --- for prognostic cloud --- !
! clouds(:,:,1) - layer total cloud fraction !
! clouds(:,:,2) - layer cloud liq water path !
! clouds(:,:,3) - mean effective radius for liquid cloud !
! clouds(:,:,4) - layer cloud ice water path !
! clouds(:,:,5) - mean effective radius for ice cloud !
! clouds(:,:,6) - layer rain drop water path !
! clouds(:,:,7) - mean effective radius for rain drop !
! clouds(:,:,8) - layer snow flake water path !
! clouds(:,:,9) - mean effective radius for snow flake !
! --- for diagnostic cloud --- !
! clouds(:,:,1) - layer total cloud fraction !
! clouds(:,:,2) - layer cloud optical depth !
! clouds(:,:,3) - layer cloud single scattering albedo !
! clouds(:,:,4) - layer cloud asymmetry factor !
! !
! 3. surface albedo: (defined in 'module_radiation_surface') !
! sfcalb( :,1 ) - near ir direct beam albedo !
! sfcalb( :,2 ) - near ir diffused albedo !
! sfcalb( :,3 ) - uv+vis direct beam albedo !
! sfcalb( :,4 ) - uv+vis diffused albedo !
! !
! 4. sw aerosol profiles: (defined in 'module_radiation_aerosols') !
! faersw(:,:,:,1)- sw aerosols optical depth !
! faersw(:,:,:,2)- sw aerosols single scattering albedo !
! faersw(:,:,:,3)- sw aerosols asymmetry parameter !
! !
! 5. lw aerosol profiles: (defined in 'module_radiation_aerosols') !
! faerlw(:,:,:,1)- lw aerosols optical depth !
! faerlw(:,:,:,2)- lw aerosols single scattering albedo !
! faerlw(:,:,:,3)- lw aerosols asymmetry parameter !
! !
! 6. sw fluxes at toa: (defined in 'module_radsw_main') !
! (topfsw_type -- derived data type for toa rad fluxes) !
! topfsw(:)%upfxc - total sky upward flux at toa !
! topfsw(:)%dnfxc - total sky downward flux at toa !
! topfsw(:)%upfx0 - clear sky upward flux at toa !
! !
! 7. lw fluxes at toa: (defined in 'module_radlw_main') !
! (topflw_type -- derived data type for toa rad fluxes) !
! topflw(:)%upfxc - total sky upward flux at toa !
! topflw(:)%upfx0 - clear sky upward flux at toa !
! !
! 8. sw fluxes at sfc: (defined in 'module_radsw_main') !
! (sfcfsw_type -- derived data type for sfc rad fluxes) !
! sfcfsw(:)%upfxc - total sky upward flux at sfc !
! sfcfsw(:)%dnfxc - total sky downward flux at sfc !
! sfcfsw(:)%upfx0 - clear sky upward flux at sfc !
! sfcfsw(:)%dnfx0 - clear sky downward flux at sfc !
! !
! 9. lw fluxes at sfc: (defined in 'module_radlw_main') !
! (sfcflw_type -- derived data type for sfc rad fluxes) !
! sfcflw(:)%upfxc - total sky upward flux at sfc !
! sfcflw(:)%dnfxc - total sky downward flux at sfc !
! sfcflw(:)%dnfx0 - clear sky downward flux at sfc !
! !
!! optional radiation outputs: !
!! 10. sw flux profiles: (defined in 'module_radsw_main') !
!! (profsw_type -- derived data type for rad vertical profiles) !
!! fswprf(:,:)%upfxc - total sky upward flux !
!! fswprf(:,:)%dnfxc - total sky downward flux !
!! fswprf(:,:)%upfx0 - clear sky upward flux !
!! fswprf(:,:)%dnfx0 - clear sky downward flux !
!! !
!! 11. lw flux profiles: (defined in 'module_radlw_main') !
!! (proflw_type -- derived data type for rad vertical profiles) !
!! flwprf(:,:)%upfxc - total sky upward flux !
!! flwprf(:,:)%dnfxc - total sky downward flux !
!! flwprf(:,:)%upfx0 - clear sky upward flux !
!! flwprf(:,:)%dnfx0 - clear sky downward flux !
!! !
!! 12. sw sfc components: (defined in 'module_radsw_main') !
!! (cmpfsw_type -- derived data type for component sfc fluxes) !
!! scmpsw(:)%uvbfc - total sky downward uv-b flux at sfc !
!! scmpsw(:)%uvbf0 - clear sky downward uv-b flux at sfc !
!! scmpsw(:)%nirbm - total sky sfc downward nir direct flux !
!! scmpsw(:)%nirdf - total sky sfc downward nir diffused flux !
!! scmpsw(:)%visbm - total sky sfc downward uv+vis direct flx !
!! scmpsw(:)%visdf - total sky sfc downward uv+vis diff flux !
! !
! ====================== end of definations ======================= !
!
implicit none
! --- constant parameter
! The following need to be commented out after the newer version
! (Hsin-mu Lin, 2011-04-25)
! integer, parameter :: NSPC = 6
! --- inputs: (for rank>1 arrays, horizontal dimensioned by IX)
integer, intent(in) :: IX,IM, LM, NTRAC, NFXR, me, &
& ntoz, ntcw, ncld, ipt, kdt, &
& nrads, cld_fraction
integer, intent(in) :: its,jts
logical, intent(in) :: isday ! true = day, false = night
integer, intent(in) :: icsdsw(IX), icsdlw(IX), jdate(8)
logical, intent(in) :: lsswr, lslwr, lssav, lprnt
real (kind=kind_phys), dimension(IX,LM+1), intent(in) :: prsi
real (kind=kind_phys), dimension(IX,LM), intent(in) :: prsl, &
& prslk, tgrs, qgrs, vvl, fcice, frain, rrime
real (kind=kind_phys), dimension(IX), intent(in) :: flgmin
real (kind=kind_phys), dimension(IX), intent(in) :: slmsk, &
& xlon, xlat, tsfc, snowd, zorl, hprim, &
! & alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
& cv, cvt, cvb, fice, tisfc, &
& sncovr, snoalb, sinlat, coslat
real (kind=kind_phys), intent(in) :: solcon, dtlw, dtsw, solhr, &
& tracer(IX,LM,NTRAC) &
& , dtswav, SALBEDO(IX), SM(IX) ! extra input
real (kind=kind_phys), dimension(IX,LM), intent(in) :: tauclouds, & ! GFDL type
& cldf
! --- outputs: (horizontal dimensioned by IX)
real (kind=kind_phys), dimension(IX,LM),intent(out):: htrsw,htrlw,&
& cldcov
real (kind=kind_phys), dimension(IX,5),intent(out):: cldsa
real (kind=kind_phys), dimension(IX), intent(out):: tsflw, &
& sfalb, semis, coszen, coszdg
real, intent(in) :: CPATHFAC4LW
type (topfsw_type), dimension(IX), intent(out) :: topfsw
type (sfcfsw_type), dimension(IX), intent(out) :: sfcfsw
type (topflw_type), dimension(IX), intent(out) :: topflw
type (sfcflw_type), dimension(IX), intent(out) :: sfcflw
! --- variables are for both input and output:
real (kind=kind_phys), intent(inout) :: fluxr(IX,NFXR)
!! --- optional outputs:
real (kind=kind_phys), dimension(IX,LM,NBDSW), optional, &
& intent(out) :: htrswb
real (kind=kind_phys), dimension(IX,LM,NBDLW), optional, &
& intent(out) :: htrlwb
! --- local variables: (horizontal dimensioned by CHK)
type (topfsw_type), dimension(CHK) :: topfsw_loc
type (sfcfsw_type), dimension(CHK) :: sfcfsw_loc
real (kind=kind_phys), dimension(CHK) :: coszen_loc !jm
integer, dimension(CHK) :: icsdsw_loc !jm
real (kind=kind_phys), dimension(CHK,LM+1+LTP)::plvl, tlvl
real (kind=kind_phys), dimension(CHK,LM+LTP) :: &
& plyr, tlyr, qlyr, &
& olyr, rhly, qstl, vvel, clw, prslk1, tem2da, tem2db, tvly, &
& cldfk1
real (kind=kind_phys), dimension(CHK,LM+LTP) :: &
& qc2d, qi2d, qs2d, ni2d
real (kind=kind_phys), dimension(CHK) :: tsfa, cvt1, cvb1, tem1d, &
& sfcemis, tsfg, tskn
real (kind=kind_phys), dimension(CHK) :: &
& alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
& CZMEAN
real (kind=kind_phys), dimension(CHK,LM+LTP,NF_CLDS) :: clouds
real (kind=kind_phys), dimension(CHK,LM+LTP,NF_CLDS) :: dummys
real (kind=kind_phys), dimension(CHK,LM+LTP,NF_VGAS) :: gasvmr
real (kind=kind_phys), dimension(CHK, NF_ALBD) :: sfcalb
! real (kind=kind_phys), dimension(CHK, NSPC1) :: aerodp ! optn for aod output
real (kind=kind_phys), dimension(CHK,LM+LTP,NTRAC) :: tracer1
real (kind=kind_phys), dimension(CHK,LM+LTP,NBDSW,NF_AESW):: faersw
real (kind=kind_phys), dimension(CHK,LM+LTP,NBDLW,NF_AELW):: faerlw
! real (kind=kind_phys), dimension(CHK,LM+LTP,NSPC-1) :: tau_gocart
real (kind=kind_phys), dimension(CHK,LM+LTP) :: htswc
real (kind=kind_phys), dimension(CHK,LM+LTP) :: htlwc
real (kind=kind_phys), dimension(CHK,LM+LTP) :: gcice, grain, grime
!! --- may be used for optional sw/lw outputs:
!! take out "!!" as needed
!! real (kind=kind_phys), dimension(CHK,LM+LTP) :: htsw0
!! type (profsw_type), dimension(CHK,LM+1+LTP) :: fswprf
type (cmpfsw_type), dimension(CHK) :: scmpsw
real (kind=kind_phys), dimension(CHK,LM+LTP,NBDSW) :: htswb
!! real (kind=kind_phys), dimension(CHK,LM+LTP) :: htlw0
!! type (proflw_type), dimension(CHK,LM+1+LTP) :: flwprf
real (kind=kind_phys), dimension(CHK,LM+LTP,NBDLW) :: htlwb
!jm real (kind=kind_phys)::raddt,es(CHK),qs, delt, tem0d, cldsa(CHK,5)&
!jm & ,SFCALBEDO(CHK), SMX(CHK)
real (kind=kind_phys)::raddt,es(CHK),qs, delt, tem0d &
& ,SFCALBEDO(CHK), SMX(CHK)
integer j2
integer :: i, j, k, k1, lv, icec, itop, ibtc, nday, idxday(CHK), &
& mbota(CHK,3), mtopa(CHK,3), LP1, nb, LMK, LMP, kd, lla, llb, &
& lya, lyb, kt, kb, np3d
! --- for debug test use
! real (kind=kind_phys) :: temlon, temlat, alon, alat
! integer :: ipt
! logical :: lprnt1
integer ipass,iorder !jm
!!==========================================================================
! --- Albedo (ETA era) calculation from 2012 RRTM (Lin, 20130225)
!
INTEGER :: IQ,JX(CHK)
! INTEGER,EXTERNAL :: omp_get_thread_num
! REAL (kind=kind_phys) :: ALBD0, ALVD1, ALND1
REAL (kind=kind_phys) :: ZEN, DZEN, ALB1, ALB2
REAL (kind=kind_phys), PARAMETER :: TWENTY=20.0, &
HP537=0.537, &
ONE=1., &
DEGRAD1=180.0/PI, &
H74E1=74.0, &
HAF=0.5, &
HNINETY=90., &
FIFTY=50., &
QUARTR=0.25, &
HNINE=9.0, &
HP1=0.1, &
H15E1=15.0
REAL (kind=kind_phys), DIMENSION(IM) :: ALVB,ALNB,ALVD,ALND
REAL (kind=kind_phys), DIMENSION(20) :: ZA
REAL (kind=kind_phys), DIMENSION(19) :: DZA
REAL (kind=kind_phys), DIMENSION(21) :: TRN
REAL (kind=kind_phys), DIMENSION(21,20) :: ALBD
DATA TRN/.00,.05,.10,.15,.20,.25,.30,.35,.40,.45,.50,.55,.60,.65, &
.70,.75,.80,.85,.90,.95,1.00/
DATA ALBD/.061,.062,.072,.087,.115,.163,.235,.318,.395,.472,.542, &
.604,.655,.693,.719,.732,.730,.681,.581,.453,.425,.061,.062,.070, &
.083,.108,.145,.198,.263,.336,.415,.487,.547,.595,.631,.656,.670, &
.652,.602,.494,.398,.370,.061,.061,.068,.079,.098,.130,.174,.228, &
.290,.357,.424,.498,.556,.588,.603,.592,.556,.488,.393,.342,.325, &
.061,.061,.065,.073,.086,.110,.150,.192,.248,.306,.360,.407,.444, &
.469,.480,.474,.444,.386,.333,.301,.290,.061,.061,.065,.070,.082, &
.101,.131,.168,.208,.252,.295,.331,.358,.375,.385,.377,.356,.320, &
.288,.266,.255,.061,.061,.063,.068,.077,.092,.114,.143,.176,.210, &
.242,.272,.288,.296,.300,.291,.273,.252,.237,.266,.220,.061,.061, &
.062,.066,.072,.084,.103,.127,.151,.176,.198,.219,.236,.245,.250, &
.246,.235,.222,.211,.205,.200, &
.061,.061,.061,.065,.071,.079,.094,.113,.134,.154,.173, &
.185,.190,.193,.193,.190,.188,.185,.182,.180,.178,.061,.061,.061, &
.064,.067,.072,.083,.099,.117,.135,.150,.160,.164,.165,.164,.162, &
.160,.159,.158,.157,.157,.061,.061,.061,.062,.065,.068,.074,.084, &
.097,.111,.121,.127,.130,.131,.131,.130,.129,.127,.126,.125,.122, &
.061,.061,.061,.061,.062,.064,.070,.076,.085,.094,.101,.105,.107, &
.106,.103,.100,.097,.096,.095,.095,.095,.061,.061,.061,.060,.061, &
.062,.065,.070,.075,.081,.086,.089,.090,.088,.084,.080,.077,.075, &
.074,.074,.074,.061,.061,.060,.060,.060,.061,.063,.065,.068,.072, &
.076,.077,.076,.074,.071,.067,.064,.062,.061,.061,.061,.061,.061, &
.060,.060,.060,.060,.061,.062,.065,.068,.069,.069,.068,.065,.061, &
.058,.055,.054,.053,.052,.052, &
.061,.061,.060,.060,.060,.060,.060,.060,.062,.065,.065, &
.063,.060,.057,.054,.050,.047,.046,.045,.044,.044,.061,.061,.060, &
.060,.060,.059,.059,.059,.059,.059,.058,.055,.051,.047,.043,.039, &
.035,.033,.032,.031,.031,.061,.061,.060,.060,.060,.059,.059,.058, &
.057,.056,.054,.051,.047,.043,.039,.036,.033,.030,.028,.027,.026, &
.061,.061,.060,.060,.060,.059,.059,.058,.057,.055,.052,.049,.045, &
.040,.036,.032,.029,.027,.026,.025,.025,.061,.061,.060,.060,.060, &
.059,.059,.058,.056,.053,.050,.046,.042,.038,.034,.031,.028,.026, &
.025,.025,.025,.061,.061,.060,.060,.059,.058,.058,.057,.055,.053, &
.050,.046,.042,.038,.034,.030,.028,.029,.025,.025,.025/
!
DATA ZA/90.,88.,86.,84.,82.,80.,78.,76.,74.,70.,66.,62.,58.,54., &
50.,40.,30.,20.,10.,0.0/
!
DATA DZA/8*2.0,6*4.0,5*10.0/
!
! --- end of Albedo (ETA era) calculation from 2012 RRTM
!!==========================================================================
!
!===> ... begin here
!
np3d = icmphys
LP1 = LM + 1 ! num of in/out levels
!===========================================================================
! --- ... set local /level/layer indexes corresponding to in/out variables
! ( for optional extra top layer)
LMK = LM + LTP ! num of local layer
LMP = LMK + 1 ! num of local level
if ( lextop ) then
if ( ivflip == 1 ) then ! vertical from sfc upward
kd = 0 ! index diff between in/out and local
kt = 1 ! index diff between lyr and upper bound
kb = 0 ! index diff between lyr and lower bound
lla = LMK ! local index at the 2nd level from top
llb = LMP ! local index at toa level
lya = LM ! local index for the 2nd layer from top
lyb = LP1 ! local index for the top layer
else ! vertical from toa downward
kd = 1 ! index diff between in/out and local
kt = 0 ! index diff between lyr and upper bound
kb = 1 ! index diff between lyr and lower bound
lla = 2 ! local index at the 2nd level from top
llb = 1 ! local index at toa level
lya = 2 ! local index for the 2nd layer from top
lyb = 1 ! local index for the top layer
endif ! end if_ivflip_block
else
kd = 0
if ( ivflip == 1 ) then ! vertical from sfc upward
kt = 1 ! index diff between lyr and upper bound
kb = 0 ! index diff between lyr and lower bound
else ! vertical from toa downward
kt = 0 ! index diff between lyr and upper bound
kb = 1 ! index diff between lyr and lower bound
endif ! end if_ivflip_block
endif ! end if_lextop_block
!
! --- End of index setting for optional extra top layer
!===========================================================================
raddt = min(dtsw, dtlw)
! --- ... for debug test
! alon = 120.0
! alat = 29.5
! ipt = 0
! do i = 1, IM
! temlon = xlon(i) * 57.29578
! if (temlon < 0.0) temlon = temlon + 360.0
! temlat = xlat(i) * 57.29578
! lprnt1 = abs(temlon-alon) < 1.1 .and. abs(temlat-alat) < 1.1
! if ( lprnt1 ) then
! ipt = i
! exit
! endif
! enddo
! print *,' in grrad : raddt=',raddt
! --- ... setup surface ground temp and ground/air skin temp if required
if ( itsfc == 0 ) then ! use same sfc skin-air/ground temp
do i = 1, IM
tskn(i) = tsfc(i)
tsfg(i) = tsfc(i)
enddo
else ! use diff sfc skin-air/ground temp
do i = 1, IM
!! tskn(i) = ta (i) ! not yet
!! tsfg(i) = tg (i) ! not yet
tskn(i) = tsfc(i)
tsfg(i) = tsfc(i)
enddo
endif
! --- ... prepare atmospheric profiles for radiation input
! convert pressure unit from Pa to mb
! if (im > ipt) then
! write(0,*)' prsi=',prsi(ipt,1:10)
! write(0,*)' prsi=',prsl(ipt,1:10)
! write(0,*)' tgrs=',tgrs(ipt,1:10)
! endif
!=============================================================
! ** NOTE: For NMMB, the conversion is "from cb to mb"
!=============================================================
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
!DEC$ SIMD
do i = 1, IM
plvl(i,k1) = 10.0 * prsi(i,k) ! cb to mb
plyr(i,k1) = 10.0 * prsl(i,k) ! cb to mb
! plvl(i,k1) = 0.01 * prsi(i,k) ! Pa to mb
! plyr(i,k1) = 0.01 * prsl(i,k) ! Pa to mb
tlyr(i,k1) = tgrs(i,k)
prslk1(i,k1) = prslk(i,k)
ENDDO
! --- ... compute relative humidity
! Note for qs (**Saturation specific humidity):
! ** prsl unit in NMMB is kPa
! it can be done in es (**Saturation vapor pressure) for 0.001*fpvs
! or multiple every prsl by 1000
! (Hsin-mu Lin, 20110927)
! note this loop won't vectorize because of the call to fpvs
DO i = 1, IM
es(i) = min( prsl(i,k), 0.001 * fpvs( tgrs(i,k) ) ) ! fpvs in pa
ENDDO
!DEC$ SIMD
DO i = 1, IM
qs = max( QMIN,con_eps*es(i)/(prsl(i,k) + con_epsm1*es(i)) )
rhly(i,k1) = max( 0.0, min( 1.0, max(QMIN, qgrs(i,k))/qs ) )
qstl(i,k1) = qs
enddo
do j = 1, NTRAC
do i = 1, IM
tracer1(i,k1,j) = tracer(i,k,j)
enddo
enddo
enddo
do i = 1, IM
plvl(i,LP1+kd) = 10.0 * prsi(i,LP1) ! cb to mb
! plvl(i,LP1+kd) = 0.01 * prsi(i,LP1) ! Pa to mb
enddo
!======================================================================
! --- ... values for extra top layer
if ( lextop ) then ! values for extra top layer
!DEC$ SIMD
do i = 1, IM
plvl(i,llb) = prsmin
if ( plvl(i,lla)<=prsmin ) plvl(i,lla) = 2.0*prsmin
plyr(i,lyb) = 0.5 * plvl(i,lla)
tlyr(i,lyb) = tlyr(i,lya)
rhly(i,lyb) = rhly(i,lya)
qstl(i,lyb) = qstl(i,lya)
enddo
! loop won't vectorize because of call to fpkapx
do i = 1, IM
prslk1(i,lyb) = fpkapx(plyr(i,lyb)*100.0) ! fpkapx in pa
enddo
! --- note: may need to take care the top layer mount
do i = 1, IM
tracer1(i,lyb,1) = tracer1(i,lya,1) ! 1st tracer element is ozone
enddo
do j = 2, NTRAC ! This loop variable should now be correct
do i = 1, IM
tracer1(i,lyb,j) = 0. ! enforce no cloud above domain top
enddo
enddo
endif
!
!======================================================================
!======================================================================
! --- ... extra variables needed for using Thompson cloud fraction in
! progcld2 & progcld8 for rrtm. For non Thompson cloud fraction
! the value will be "0" as initialization in module_RA_RRTM
!
do k = 1, LM
k1 = k + kd
do i = 1, IM
cldfk1(i,k1)= cldf(i,k)
enddo
enddo
!========================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
cldfk1(i,lyb) = cldf(i,lya)
enddo
endif
!
!======================================================================
!======================================================================
! --- ... extra variables needed for ferrier's microphysics
if (np3d == 3) then
do k = 1, LM
k1 = k + kd
do i = 1, IM
gcice(i,k1)= fcice(i,k)
grain(i,k1)= frain(i,k)
grime(i,k1)= rrime(i,k)
enddo
enddo
!========================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
gcice(i,lyb) = 1.0 ! fcice(i,lya)
grain(i,lyb) = 0.0 ! frain(i,lya)
grime(i,lyb) = 1.0 ! rrime(i,lya)
enddo
endif
endif ! if_np3d
!
!======================================================================
! --- ... get layer ozone mass mixing ratio
if (ntoz > 0) then ! interactive ozone generation
do k = 1, LMK
do i = 1, IM
olyr(i,k) = max( QMIN, tracer1(i,k,ntoz) )
enddo
enddo
else ! climatological ozone
! print *,' in grrad : calling getozn'
call getozn &
! --- inputs:
& ( prslk1,xlat, &
& IM, LMK, &
! --- outputs:
& olyr &
& )
endif ! end_if_ntoz
! --- ... compute cosin of zenith angle
call coszmn_nmmb &
! --- inputs:
& ( xlon,sinlat,coslat,solhr,IM, me, &
& dtswav,nrads , & ! Extra input
! --- outputs:
& coszen_loc, coszdg & !jm
& )
! --- ... set up non-prognostic gas volume mixing ratioes
call getgases &
! --- inputs:
& ( plvl, xlon, xlat, &
& IM, LMK, &
! --- outputs:
& gasvmr &
& )
! --- ... get temperature at layer interface, and layer moisture
do k = 2, LMK
do i = 1, IM
tem2da(i,k) = log( plyr(i,k) )
tem2db(i,k) = log( plvl(i,k) )
enddo
enddo
if (ivflip == 0) then ! input data from toa to sfc
do i = 1, IM
tem1d (i) = QME6
tem2da(i,1) = log( plyr(i,1) )
!---------------------------------------------
! take care for the model top (20130410, Lin)
!---------------------------------------------
if ( plvl(i,1) == 0.0 ) then
tem2db(i,1) = -14.0
else
tem2db(i,1) = log( plvl(i,1) )
endif
! tem2db(i,1) = 1.0
tsfa (i) = tlyr(i,LMK) ! sfc layer air temp
tlvl(i,1) = tlyr(i,1)
tlvl(i,LMP) = tskn(i)
enddo
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM
qlyr(i,k1) = max( tem1d(i), qgrs(i,k) )
tem1d(i) = min( QME5, qlyr(i,k1) )
tvly(i,k1) = tgrs(i,k) * (1.0 + con_fvirt*qlyr(i,k1)) ! virtual temp in K
enddo
enddo
!============================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
qlyr(i,lyb) = qlyr(i,lya)
tvly(i,lyb) = tvly(i,lya)
enddo
endif
!============================================
do k = 2, LMK
do i = 1, IM
tlvl(i,k) = tlyr(i,k) + (tlyr(i,k-1) - tlyr(i,k)) &
& * (tem2db(i,k) - tem2da(i,k)) &
& / (tem2da(i,k-1) - tem2da(i,k))
enddo
enddo
else ! input data from sfc to toa
do i = 1, IM
tem1d (i) = QME6
tem2da(i,1) = log( plyr(i,1) )
tem2db(i,1) = log( plvl(i,1) )
tsfa (i) = tlyr(i,1) ! sfc layer air temp
tlvl(i,1) = tskn(i)
tlvl(i,LMP) = tlyr(i,LMK)
enddo
do k = LM, 1, -1
do i = 1, IM
qlyr(i,k) = max( tem1d(i), qgrs(i,k) )
tem1d(i) = min( QME5, qlyr(i,k) )
tvly(i,k) = tgrs(i,k) * (1.0 + con_fvirt*qlyr(i,k)) ! virtual temp in K
enddo
enddo
!============================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
qlyr(i,lyb) = qlyr(i,lya)
tvly(i,lyb) = tvly(i,lya)
enddo
endif
!============================================
do k = 1, LMK-1
do i = 1, IM
tlvl(i,k+1) = tlyr(i,k) + (tlyr(i,k+1) - tlyr(i,k)) &
& * (tem2db(i,k+1) - tem2da(i,k)) &
& / (tem2da(i,k+1) - tem2da(i,k))
enddo
enddo
endif ! end_if_ivflip
! --- check for daytime points
!jm with modification in RA_RRTM, the grrad_nmm routine will
!jm only be called for all day or all night points, the number
!jm of which is passed in as IM.
nday = im
! --- ... setup aerosols property profile for radiation
!check print *,' in grrad : calling setaer '
call setaer &
! --- inputs:
& ( plvl,plyr,prslk1,tvly,rhly,slmsk,tracer1,xlon,xlat, &
& IM,LMK,LMP, lsswr,lslwr, &
! --- extra 2 vars for old code that has been replaced by 'tvly'
! can be eliminated after aggrements being reached
& tlyr,qlyr, &
! --- outputs:
& faersw,faerlw &
! &, tau_gocart &
& )
! --- ... obtain cloud information for radiation calculations
if (ntcw > 0) then ! prognostic cloud scheme
do k = 1, LMK
do i = 1, IM
clw(i,k) = 0.0
enddo
do j = 1, ncld
lv = ntcw + j - 1
do i = 1, IM
clw(i,k) = clw(i,k) + tracer1(i,k,lv) ! cloud condensate amount
enddo
enddo
enddo
do k = 1, LMK
do i = 1, IM
if ( clw(i,k) < EPSQ ) clw(i,k) = 0.0
enddo
enddo
if (np3d == 4) then ! zhao/moorthi's prognostic cloud scheme
call progcld1 &
! --- inputs:
& ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, &
& xlat,xlon,slmsk, &
& IM, LMK, LMP, &
! --- outputs:
& clouds,cldsa,mtopa,mbota &
& )
elseif (np3d == 3) then ! ferrier's microphysics
! print *,' in grrad : calling progcld2'
call progcld2 &
! --- inputs:
& ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, &
& xlat,xlon,slmsk, gcice,grain,grime,flgmin, &
& IM, LMK, LMP, &
& cldfk1, cld_fraction, &
! --- outputs:
& clouds,cldsa,mtopa,mbota &
& )
elseif (np3d == 5) then ! nmmb ferrier+bmj (GFDL type diagnostic)
!
! ======================================================================
! The original GFS RRTM will use "call progcld3".
! ** After merging the nmmb & GFS RRTM, "np3d=5" is temporarily used
! for the GFDL diagnostic cloud type which is in the original meso
! nmmb (Hsin-mu Lin 12/28/2011)
! ======================================================================
!
! call progcld3 &
! --- inputs:
! & ( plyr,plvl,tlyr,qlyr,qstl,rhly,clw, &
! & xlat,xlon,slmsk, fcice,frain,rrime,flgmin, &
! & IM, LMK, LMP, iflip, iovrsw, &
! --- outputs:
! & clouds,cldsa,mtopa,mbota &
! & )
! print *,' in grrad : !!! need to develop progcld3 for nmmb', &
! & ' temporarily using progcld2 !!!'
do k = 1, LM ! GFDL type
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM ! GFDL type
clouds(i,k1,1) = cldf(i,k) ! GFDL type
clouds(i,k1,2) = tauclouds(i,k) ! GFDL type
clouds(i,k1,3) = 0.99 ! GFDL type
clouds(i,k1,4) = 0.84 ! GFDL type
enddo ! GFDL type
enddo ! GFDL type
!============================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
clouds(i,lyb,1) = clouds(i,lya,1)
clouds(i,lyb,2) = clouds(i,lya,2)
clouds(i,lyb,3) = clouds(i,lya,3)
clouds(i,lyb,4) = clouds(i,lya,4)
enddo
endif
!============================================
!..This is far from optimal! The specific array element numbers are
!.. potentially flexible for cloud water, cloud ice, snow mixing ratios
!.. and cloud ice number concentration. Therefore, numerical values of
!.. 4, 5, 6, 9 are hard-wired when they should be made flexible.
elseif (np3d == 8) then ! Thompson microphysics
! write(*,*) ' DEBUG-GT, double-check radiation flags ',iswcliq , iswcice, ilwcliq, ilwcice
! write(*,*) ' DEBUG-GT, double-check NTRAC number ', ntrac
! write(*,*) ' DEBUG-GT, double-check LMK (LM+LTP) ', LMK
do k = 1, LMK
do i = 1, IM
qc2d(i,k) = tracer1(i,k,4)
if (qc2d(i,k) .LE. EPSQ) qc2d(i,k)=0.0
qi2d(i,k) = tracer1(i,k,5)
if (qi2d(i,k) .LE. EPSQ) qi2d(i,k)=0.0
qs2d(i,k) = tracer1(i,k,6)
if (qs2d(i,k) .LE. EPSQ) qs2d(i,k)=0.0
ni2d(i,k) = MAX(1.E-8, tracer1(i,k,9))
enddo
enddo
call progcld8 (plyr,plvl, tlyr, qlyr, qc2d, qi2d, qs2d, ni2d, &
& xlat, IM, LMK, LMP, &
& cldfk1, cld_fraction, &
& clouds, cldsa, mtopa, mbota)
endif ! end if_np3d
!
else ! diagnostic cloud scheme
do i = 1, IM
cvt1(i) = 10.0 * cvt(i)
cvb1(i) = 10.0 * cvb(i)
! cvt1(i) = 0.01 * cvt(i)
! cvb1(i) = 0.01 * cvb(i)
enddo
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM
vvel(i,k1) = 10.0 * vvl (i,k)
! vvel(i,k1) = 0.01 * vvl (i,k)
enddo
enddo
!============================================
! --- ... values for extra top layer
if ( lextop ) then
do i = 1, IM
vvel(i,lyb) = vvel(i,lya)
enddo
endif
!============================================
! --- compute diagnostic cloud related quantities
call diagcld1 &
! --- inputs:
& ( plyr,plvl,tlyr,rhly,vvel,cv,cvt1,cvb1, &
& xlat,xlon,slmsk, &
& IM, LMK, LMP, &
! --- outputs:
& clouds,cldsa,mtopa,mbota &
& )
endif ! end_if_ntcw
!==========================================================================
! Albedo (ETA era) calculation from 2012 RRTM (Lin, 20130225)
!==========================================================================
IF (ialbflg == 2) THEN
do i = 1, IM
SMX(i) = SM(i)
SFCALBEDO(i) = SALBEDO(i)
!..... THE FOLLOWING CODE GETS ALBEDO FROM PAYNE,1972 TABLES IF
! 1) OPEN SEA POINT (SLMSK=1); 2) KALB=0
IQ=INT(TWENTY*HP537+ONE)
CZMEAN(i)=coszen_loc(i) !jm
IF(CZMEAN(i).GT.0.0 .AND. SMX(i).GT.0.5) THEN
ZEN=DEGRAD1*ACOS(MAX(CZMEAN(i),0.0))
IF(ZEN.GE.H74E1) JX(I)=INT(HAF*(HNINETY-ZEN)+ONE)
IF(ZEN.LT.H74E1 .AND. ZEN.GE.FIFTY) &
JX(I)=INT(QUARTR*(H74E1-ZEN)+HNINE)
IF(ZEN .LT. FIFTY) JX(I)=INT(HP1*(FIFTY-ZEN)+H15E1)
DZEN=-(ZEN-ZA(JX(I)))/DZA(JX(I))
ALB1=ALBD(IQ,JX(I))+DZEN*(ALBD(IQ,JX(I)+1)&
-ALBD(IQ,JX(I)))
ALB2=ALBD(IQ+1,JX(I))+DZEN*(ALBD(IQ+1,JX(I)+1)&
-ALBD(IQ+1,JX(I)))
SFCALBEDO(I)=ALB1+TWENTY*(ALB2-ALB1)*(HP537-TRN(IQ)) ! BSF
ENDIF
!..... VISIBLE AND NEAR IR DIFFUSE ALBEDO
ALVD(I) = SFCALBEDO(I) ! BSF
ALND(I) = SFCALBEDO(I) ! BSF
!..... VISIBLE AND NEAR IR DIRECT BEAM ALBEDO
ALVB(I) = SFCALBEDO(I) ! BSF
ALNB(I) = SFCALBEDO(I) ! BSF
!
!--- Remove diurnal variation of land surface albedos (Ferrier, 6/28/05)
!--- Turn back on to mimic NAM 8/17/05
!
!..... VISIBLE AND NEAR IR DIRECT BEAM ALBEDO, IF NOT OCEAN NOR SNOW
! ..FUNCTION OF COSINE SOLAR ZENITH ANGLE..
!
!=== The following line commented out by the "fixed" are used instead of the
! original version (09/2012)
!
!fixed IF (SMX .LT. 0.5) THEN
!fixed IF (SFCALBEDO .LE. 0.5) THEN
!fixed ALBD0=-18.0 * (0.5 - ACOS(CZMEAN(I))/PI)
!fixed ALBD0=EXP (ALBD0)
!fixed ALVD1=(ALVD(I) - 0.054313) / 0.945687
!fixed ALND1=(ALND(I) - 0.054313) / 0.945687
!fixed ALVB(I)=ALVD1 + (1.0 - ALVD1) * ALBD0
!fixed ALNB(I)=ALND1 + (1.0 - ALND1) * ALBD0
!fixed! !-- Put in an upper limit on beam albedos
!fixed ALVB(I)=MIN(0.5,ALVB(I))
!fixed ALNB(I)=MIN(0.5,ALNB(I))
!fixed ENDIF
!fixed ENDIF
!!! WE INTRODUCE HERE DIRECT AND DIFFUSE ALBEDO... FOR THIS OPTION, THERE IS A CHANGE IN GRRAD.f
alvsf(i) = ALVB(I) !For this option ALVSF1 is direct visible albedo
alnsf(i) = ALNB(I) !For this option ALNSF1 is direct nir albedo
alvwf(i) = ALVD(I) !For this option ALVWF1 is diffuse visible albedo
alnwf(i) = ALND(I) !For this option ALNWF1 is diffuse nir albedo
facsf(i) = 0. !not used with this option
facwf(i) = 0. !not used for this option
ENDDO
ENDIF !---- end of IALB=2 GFDL TYPE RADIATION
!
! --- End of Albedo (ETA era) calculation from 2012 RRTM ------------------
!==========================================================================
! --- ... start radiation calculations
! remember to set heating rate unit to k/sec!
!#ifdef ORIG_LWSW_ORDERING
iorder = 0
!#else
! iorder = mod(omp_get_thread_num(),2) ! every other thread
! iorder = mod(omp_get_thread_num()/2,2) ! alternate by groups of 4
!#endif
DO ipass=0,1
IF( ipass .eq. iorder ) then
if (lsswr) then
! --- setup surface albedo for sw radiation, incl xw (nov04) sea-ice
call setalb &
! --- inputs:
& ( slmsk,snowd,sncovr,snoalb,zorl,coszen_loc,tsfg,tsfa,hprim, & !jm
& alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc, &
& IM, &
! --- outputs:
& sfcalb &
& )
! --- lu [+4L]: derive SFALB from vis- and nir- diffuse surface albedo
do i = 1, IM
sfalb(i) = max(0.01, 0.5 * (sfcalb(i,2) + sfcalb(i,4)))
enddo
! if (nday > 0) then
if (isday) then
! print *,' in grrad : calling swrad'
!jm pad out vector
#if 1
if ( nday < CHK ) then
coszen_loc(nday+1:CHK) = coszen_loc(nday)
do k = 1,lmp
plvl (nday+1:CHK,k) = plvl (nday,k)
tlvl (nday+1:CHK,k) = tlvl (nday,k)
enddo
do k = 1,lmk
plyr (nday+1:CHK,k) = plyr (nday,k)
tlyr (nday+1:CHK,k) = tlyr (nday,k)
qlyr (nday+1:CHK,k) = qlyr (nday,k)
olyr (nday+1:CHK,k) = olyr (nday,k)
enddo
do j = 1,4
sfcalb (nday+1:CHK,j) = sfcalb (nday,j)
enddo
do j = 1,9
do k = 1,lmk
gasvmr (nday+1:CHK,k,j) = gasvmr (nday,k,j)
clouds (nday+1:CHK,k,j) = clouds (nday,k,j)
enddo
enddo
do j = 1,3
do j2 = 1,nbdsw
do k = 1,lmk
faersw (nday+1:CHK,k,j2,j) = faersw (nday,k,j2,j)
enddo
enddo
enddo
endif
#endif
if ( present(htrswb) ) then
call swrad &
! --- inputs:
& ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
& clouds,icsdsw_loc,faersw,sfcalb, &
& coszen_loc,solcon, nday,idxday, &
& IM, LMK, LMP, lprnt, &
! --- outputs:
& htswc,topfsw_loc,sfcfsw_loc &
!! --- optional:
!! &, HSW0=htsw0,FLXPRF=fswprf &
&, HSWB=htswb,FDNCMP=scmpsw &
&, ITS=its,JTS=jts &
& )
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do j = 1, NBDSW
do i = 1, IM
htrswb(i,k,j) = htswb(i,k1,j)
enddo
enddo
enddo
else
call swrad &
! --- inputs:
& ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
& clouds,icsdsw_loc,faersw,sfcalb, &
& coszen_loc,solcon, nday,idxday, &
& IM, LMK, LMP, lprnt, &
! --- outputs:
& htswc,topfsw_loc,sfcfsw_loc &
!! --- optional:
!! &, HSW0=htsw0,FLXPRF=fswprf,HSWB=htswb &
&, FDNCMP=scmpsw &
&, ITS=its,JTS=jts &
& )
endif
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM
htrsw(i,k) = htswc(i,k1)
enddo
enddo
else ! if_isday_block
do k = 1, LM
do i = 1, IM
htrsw(i,k) = 0.0
enddo
enddo
sfcfsw_loc = sfcfsw_type( 0.0, 0.0, 0.0, 0.0 )
topfsw_loc = topfsw_type( 0.0, 0.0, 0.0 )
scmpsw = cmpfsw_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 )
!! --- optional:
!! fswprf= profsw_type( 0.0, 0.0, 0.0, 0.0 )
if ( present(htrswb) ) then
do j = 1, NBDSW
do k = 1, LM
do i = 1, IM
htrswb(i,k,j) = 0.0
enddo
enddo
enddo
endif
endif ! end_if_isday
topfsw(1:IM) = topfsw_loc(1:IM)
sfcfsw(1:IM) = sfcfsw_loc(1:IM)
endif ! end_if_lsswr
ELSE ! ipass == iorder
if (lslwr) then
! --- setup surface emissivity for lw radiation
call setemis &
! --- inputs:
& ( xlon,xlat,slmsk,snowd,sncovr,zorl,tsfg,tsfa,hprim, &
& IM, &
! --- outputs:
& sfcemis &
& )
! print *,' in grrad : calling lwrad'
!==========================================================================
! Special for lwrad to enhence the emissivity
! it is similar to *CPATHFAC4LW to odcld in radlw (Hsin-Mu Lin, 20140520)
!==========================================================================
if (np3d == 3) then ! ferrier's microphysics
do k = 1, LMK
do i = 1, IM
dummys(i,k,2) = clouds(i,k,2)
dummys(i,k,4) = clouds(i,k,4)
clouds(i,k,2) = dummys(i,k,2)*CPATHFAC4LW ! cloud liquid path
clouds(i,k,4) = dummys(i,k,4)*CPATHFAC4LW ! cloud ice path
enddo
enddo
endif
!==========================================================================
if ( present(htrlwb) ) then
call lwrad &
! --- inputs:
& ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
& clouds,icsdlw,faerlw,sfcemis,tsfg, &
& IM, LMK, LMP, lprnt, &
! --- outputs:
& htlwc,topflw,sfcflw &
!! --- optional:
!! &, HLW0=htlw0,FLXPRF=flwprf &
&, HLWB=htlwb &
&, ITS=its,JTS=jts &
& )
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do j = 1, NBDLW
do i = 1, IM
htrlwb(i,k,j) = htlwb(i,k1,j)
enddo
enddo
enddo
else
! do i = 1, IM
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k sfcemis",3i7,e25.15)')its+i-1,jts,k,sfcemis(i)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tsfg",3i7,e25.15)')its+i-1,jts,k,tsfg(i)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k icsdlw",4i7)')its+i-1,jts,k,icsdlw(i)
! do k = 1, LM
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k plyr",3i7,e25.15)')its+i-1,jts,k,plyr(i,k)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k plvl",3i7,e25.15)')its+i-1,jts,k,plvl(i,k)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tlyr",3i7,e25.15)')its+i-1,jts,k,tlyr(i,k)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k tlvl",3i7,e25.15)')its+i-1,jts,k,tlvl(i,k)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k qlyr",3i7,e25.15)')its+i-1,jts,k,qlyr(i,k)
! if (its+i-1.eq.140.and.jts.eq.47)write(0,'("two i,j,k olyr",3i7,e25.15)')its+i-1,jts,k,olyr(i,k)
! enddo
! enddo
call lwrad &
! --- inputs:
& ( plyr,plvl,tlyr,tlvl,qlyr,olyr,gasvmr, &
& clouds,icsdlw,faerlw,sfcemis,tsfg, &
& IM, LMK, LMP, lprnt, &
! --- outputs
& htlwc,topflw,sfcflw &
!! --- optional:
!! &, HLW0=htlw0,FLXPRF=flwprf,HLWB=htlwb &
&, ITS=its,JTS=jts &
& )
endif
do i = 1, IM
semis (i) = sfcemis(i)
! --- save surface air temp for diurnal adjustment at model t-steps
tsflw (i) = tsfa(i)
enddo
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM
htrlw(i,k) = htlwc(i,k1)
! if (k.eq.1.and.its+i-1.eq.140.and.jts.eq.47)write(0,'("two IM,its,jts",2i7,e25.15)')its+i-1,jts,htrlw(i,1)
enddo
enddo
!==========================================================================
! return *1.5 clouds to original value (Hsin-Mu Lin, 20140520)
!==========================================================================
if (np3d == 3) then ! ferrier's microphysics
do k = 1, LMK
do i = 1, IM
clouds(i,k,2) = dummys(i,k,2) ! cloud liquid path
clouds(i,k,4) = dummys(i,k,4) ! cloud ice path
enddo
enddo
endif
!==========================================================================
endif ! end_if_lslwr
ENDIF ! iorder
ENDDO ! jm ipass
! --- ... collect the fluxr data for wrtsfc
if (lssav) then
! --- save lw toa and sfc fluxes
if (lslwr) then
do i = 1, IM
! --- lw total-sky fluxes
fluxr(i,1 ) = fluxr(i,1 ) + dtlw * topflw(i)%upfxc ! total sky top lw up
fluxr(i,19) = fluxr(i,19) + dtlw * sfcflw(i)%dnfxc ! total sky sfc lw dn
fluxr(i,20) = fluxr(i,20) + dtlw * sfcflw(i)%upfxc ! total sky sfc lw up
! --- lw clear-sky fluxes
fluxr(i,28) = fluxr(i,28) + dtlw * topflw(i)%upfx0 ! clear sky top lw up
fluxr(i,30) = fluxr(i,30) + dtlw * sfcflw(i)%dnfx0 ! clear sky sfc lw dn
fluxr(i,33) = fluxr(i,33) + dtlw * sfcflw(i)%upfx0 ! clear sky sfc lw up
enddo
endif
! --- save sw toa and sfc fluxes with proper diurnal sw wgt. coszen=mean cosz over daylight
! part of sw calling interval, while coszdg= mean cosz over entire interval
if (lsswr) then
do i = 1, IM
if (coszen_loc(i) > 0.) then
! --- sw total-sky fluxes
tem0d = dtsw * coszdg(i) / coszen_loc(i)
fluxr(i,2 ) = fluxr(i,2) + topfsw(i)%upfxc * tem0d ! total sky top sw up
fluxr(i,3 ) = fluxr(i,3) + sfcfsw(i)%upfxc * tem0d ! total sky sfc sw up
fluxr(i,4 ) = fluxr(i,4) + sfcfsw(i)%dnfxc * tem0d ! total sky sfc sw dn
! --- sw uv-b fluxes
fluxr(i,21) = fluxr(i,21) + scmpsw(i)%uvbfc * tem0d ! total sky uv-b sw dn
fluxr(i,22) = fluxr(i,22) + scmpsw(i)%uvbf0 * tem0d ! clear sky uv-b sw dn
! --- sw toa incoming fluxes
fluxr(i,23) = fluxr(i,23) + topfsw(i)%dnfxc * tem0d ! top sw dn
! --- sw sfc flux components
fluxr(i,24) = fluxr(i,24) + scmpsw(i)%visbm * tem0d ! uv/vis beam sw dn
fluxr(i,25) = fluxr(i,25) + scmpsw(i)%visdf * tem0d ! uv/vis diff sw dn
fluxr(i,26) = fluxr(i,26) + scmpsw(i)%nirbm * tem0d ! nir beam sw dn !?
fluxr(i,27) = fluxr(i,27) + scmpsw(i)%nirdf * tem0d ! nir diff sw dn !?
! --- sw clear-sky fluxes
fluxr(i,29) = fluxr(i,29) + topfsw(i)%upfx0 * tem0d ! clear sky top sw up
fluxr(i,31) = fluxr(i,31) + sfcfsw(i)%upfx0 * tem0d ! clear sky sfc sw up
fluxr(i,32) = fluxr(i,32) + sfcfsw(i)%dnfx0 * tem0d ! clear sky sfc sw dn
endif
enddo
endif
! --- save total cloud and bl cloud
if (lsswr .or. lslwr) then
do i = 1, IM
fluxr(i,17) = fluxr(i,17) + raddt * cldsa(i,4) !?
fluxr(i,18) = fluxr(i,18) + raddt * cldsa(i,5) !?
enddo
! --- save cld frac,toplyr,botlyr and top temp, note that the order
! of h,m,l cloud is reversed for the fluxr output.
! --- save interface pressure (Pa) of top/bot
do j = 1, 3
do i = 1, IM
tem0d = raddt * cldsa(i,j)
itop = mtopa(i,j) - kd
ibtc = mbota(i,j) - kd
!==========================================
!=== Zavisa mentioned on 01/15/2013
! for our current run, there are no problem for np3d=3
! without this modification, while np3d=5 need this
itop=min(max(itop,1),lm) !zj quick fix
ibtc=min(max(ibtc,1),lm) !zj quick fix
!==========================================
fluxr(i, 8-j) = fluxr(i, 8-j) + tem0d
fluxr(i,11-j) = fluxr(i,11-j) + prsi(i,itop+kt) * tem0d
fluxr(i,14-j) = fluxr(i,14-j) + prsi(i,ibtc+kb) * tem0d
fluxr(i,17-j) = fluxr(i,17-j) + tgrs(i,itop) * tem0d
enddo
enddo
endif
do k = 1, LM
k1 = k + kd ! kd: index diff between in/out and local
do i = 1, IM
cldcov(i,k) = clouds(i,k1,1)
enddo
enddo
! --- save optional vertically integrated aerosol optical depth at
! wavelenth of 550nm aerodp(:,1), and other optional aod for
! individual species aerodp(:,2:NSPC1)
! if ( laswflg ) then
! if ( NFXR > 33 ) then
! do i = 1, IM
! fluxr(i,34) = fluxr(i,34) + dtsw*aerodp(i,1) ! total aod at 550nm (all species)
! enddo
! if ( lspcodp ) then
! do j = 2, NSPC1
! k = 33 + j
! do i = 1, IM
! fluxr(i,k) = fluxr(i,k) + dtsw*aerodp(i,j) ! aod at 550nm for indiv species
! enddo
! enddo
! endif ! end_if_lspcodp
! else
! print *,' !Error! Need to increase array fluxr size NFXR ',&
! & ' to be able to output aerosol optical depth'
! stop
! endif ! end_if_nfxr
! endif ! end_if_laswflg
endif ! end_if_lssav
do i = 1, IM
coszen(i) = coszen_loc(i) ! jm
enddo
!
return
!...................................
end subroutine grrad_nmmb
!-----------------------------------
!
!.............................................!
end module module_radiation_driver_nmmb !
!=============================================!