!!!!! ========================================================== !!!!! !!!!! '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 ! !=============================================!