!WRF:MEDIATION_LAYER:PHYSICS ! MODULE module_radiation_driver CONTAINS !BOP ! !IROUTINE: radiation_driver - interface to radiation physics options ! !INTERFACE: SUBROUTINE radiation_driver ( & ACFRCV ,ACFRST ,ALBEDO & ,CFRACH ,CFRACL ,CFRACM & ,CUPPT ,CZMEAN ,DT & ,DZ8W ,EMISS ,GLW & ,GMT ,GSW ,HBOT & ,HTOP ,HBOTR ,HTOPR & ,ICLOUD & ,ITIMESTEP,JULDAY, JULIAN & ,JULYR ,LW_PHYSICS & ,NCFRCV ,NCFRST ,NPHS & ,O3INPUT, O3RAD & ,AER_OPT, aerod & ,swint_opt & ,P8W ,P ,PI & ,p_top & ,RADT ,RA_CALL_OFFSET & ,RHO ,RLWTOA & ,RSWTOA ,RTHRATEN & ,RTHRATENLW ,RTHRATENSW & ,SNOW ,STEPRA ,SWDOWN & ,SWDOWNC ,SW_PHYSICS & ,SW_ECLIPSE & ! amontornes-bcodina 2015/09 solar eclipses ,T8W ,T ,TAUCLDC & ,TAUCLDI ,TSK ,VEGFRA & ,WARM_RAIN ,XICE ,XLAND & ,XLAT ,XLONG ,YR & ! tgs - needed for albedo solar angle correction & ,IVGTYP,alb_sol & & ,ALBBCK,ALBSOL,ALBBCKSOL & !Optional solar variables ,DECLINX ,SOLCONX ,COSZEN ,HRANG & , CEN_LAT & ,Z & ,ALEVSIZ, no_src_types & ,LEVSIZ, N_OZMIXM & ,N_AEROSOLC & ,PAERLEV ,ID & ,CAM_ABS_DIM1, CAM_ABS_DIM2 & ,CAM_ABS_FREQ_S & ,XTIME & ,CURR_SECS, ADAPT_STEP_FLAG & !BSINGH - For WRFCuP scheme (optional args) ,cu_physics,shallowcu_forced_ra & !CuP, wig 1-Oct-2006 ,cubot,cutop,cldfra_cup & !CuP, wig 1-Oct-2006 ,shall & !CuP, wig 4-Feb-2008 !BSINGH - ENDS ! indexes ,IDS,IDE, JDS,JDE, KDS,KDE & ,IMS,IME, JMS,JME, KMS,KME & ,i_start,i_end & ,j_start,j_end & ,kts, kte & ,num_tiles & ! Optional , TLWDN, TLWUP & ! goddard schemes , SLWDN, SLWUP & ! goddard schemes , TSWDN, TSWUP & ! goddard schemes , SSWDN, SSWUP & ! goddard schemes , CLDFRA,CLDFRA_MP_ALL,CLDT,ZNU & #if (EM_CORE == 1) , lradius,iradius & #endif , cldfra_dp, cldfra_sh & ! ckay for sub-grid cloud fraction , re_cloud, re_ice, re_snow & ! G. Thompson , has_reqc, has_reqi, has_reqs & ! G. Thompson , PB & , F_ICE_PHY,F_RAIN_PHY & , QV, F_QV & , QC, F_QC & , QR, F_QR & , QI, F_QI & , QS, F_QS & , QG, F_QG & , QNDROP, F_QNDROP & ,QNIFA,F_QNIFA & ! trude ,QNWFA,F_QNWFA & ! trude ,ACSWUPT ,ACSWUPTC & ,ACSWDNT ,ACSWDNTC & ,ACSWUPB ,ACSWUPBC & ,ACSWDNB ,ACSWDNBC & ,ACLWUPT ,ACLWUPTC & ,ACLWDNT ,ACLWDNTC & ,ACLWUPB ,ACLWUPBC & ,ACLWDNB ,ACLWDNBC & ,SWUPT ,SWUPTC & ,SWDNT ,SWDNTC & ,SWUPB ,SWUPBC & ,SWDNB ,SWDNBC & ,LWUPT ,LWUPTC & ,LWDNT ,LWDNTC & ,LWUPB ,LWUPBC & ,LWDNB ,LWDNBC & ,LWCF & ,SWCF & ,OLR & ,aerodm, PINA, aodtot & ,OZMIXM, PIN & ,M_PS_1, M_PS_2, AEROSOLC_1 & ,AEROSOLC_2, M_HYBI0 & ,ABSTOT, ABSNXT, EMSTOT & ,ICLOUD_CU & ,AER_RA_FEEDBACK & ,QC_CU , QI_CU & ,icloud_bl,qc_bl,cldfra_bl & !JOE-subgrid bl clouds ,PM2_5_DRY, PM2_5_WATER & ,PM2_5_DRY_EC & ,TAUAER300, TAUAER400 & ! jcb ,TAUAER600, TAUAER999 & ! jcb ,GAER300, GAER400, GAER600, GAER999 & ! jcb ,WAER300, WAER400, WAER600, WAER999 & ! jcb ,TAUAERlw1, TAUAERlw2 & ,TAUAERlw3, TAUAERlw4 & ,TAUAERlw5, TAUAERlw6 & ,TAUAERlw7, TAUAERlw8 & ,TAUAERlw9, TAUAERlw10 & ,TAUAERlw11, TAUAERlw12 & ,TAUAERlw13, TAUAERlw14 & ,TAUAERlw15, TAUAERlw16 & ,progn & ,smkfdb_opt, chem_opt, AOD3D_SMOKE & ! RAR: optional arguments for smoke feedback ,slope_rad,topo_shading & ,shadowmask,ht,dx,dy & ,dxkm & ,diffuse_frac & ,SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC & ! Optional ,LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC & ! Optional ,radtacttime & ,ALSWVISDIR, ALSWVISDIF, ALSWNIRDIR, ALSWNIRDIF & !fds ssib alb comp (06/2010) ,SWVISDIR, SWVISDIF, SWNIRDIR, SWNIRDIF & !fds ssib swr comp (06/2010) ,SF_SURFACE_PHYSICS, IS_CAMMGMP_USED & !fds ,EXPLICIT_CONVECTION & ! .true.=no conv. scheme ,swddir,swddni,swddif & ! jararias 2013/08 ,swddnic,swddifc & ! J. Kenyon, 30 Sep 2016 ,swdown_ref,swddir_ref,coszen_ref,Gx,gg,Bx,bb & ,aer_type & ! jararias 2013/11 ,aer_aod550_opt, aer_aod550_val & ,aer_angexp_opt, aer_angexp_val & ,aer_ssa_opt, aer_ssa_val & ,aer_asy_opt, aer_asy_val & ,aod5502d, angexp2d, aerssa2d, aerasy2d & ,aod5503d & ,obscur,mask,elat_track,elon_track & ! amontornes-bcodina 2015/09 solar eclipses ,taod5502d, taod5503d & ! Trude ,mp_physics & ,spp_rad,pattern_spp_rad & ) !------------------------------------------------------------------------- ! !USES: USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME & ,RRTMG_LWSCHEME, RRTMG_SWSCHEME & ,RRTMG_LWSCHEME_FAST, RRTMG_SWSCHEME_FAST & ,SWRADSCHEME, GSFCSWSCHEME & ,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME & ,HELDSUAREZ & #if ( HWRF == 1 ) ,HWRFSWSCHEME, HWRFLWSCHEME & #endif ,goddardlwscheme & ,goddardswscheme & ,KFCUPSCHEME & !BSINGH - Added KFCUPSCHEME for WRFCuP scheme ,FLGLWSCHEME, FLGSWSCHEME USE module_model_constants #ifndef HWRF USE module_wrf_error , ONLY : wrf_err_message #endif ! *** add new modules of schemes here USE module_ra_sw , ONLY : swrad USE module_ra_gsfcsw , ONLY : gsfcswrad USE module_ra_rrtm , ONLY : rrtmlwrad USE module_ra_rrtmg_lw , ONLY : rrtmg_lwrad USE module_ra_rrtmg_sw , ONLY : rrtmg_swrad USE module_ra_rrtmg_lwf , ONLY : rrtmg_lwrad_fast USE module_ra_rrtmg_swf , ONLY : rrtmg_swrad_fast USE module_ra_cam , ONLY : camrad USE module_ra_gfdleta , ONLY : etara #if ( HWRF == 1 ) USE module_ra_hwrf #endif USE module_ra_hs , ONLY : hsrad USE module_ra_goddard , ONLY : goddardrad USE module_ra_flg , ONLY : RAD_FLG USE module_ra_aerosol , ONLY : calc_aerosol_goddard_sw, & calc_aerosol_rrtmg_sw ! amontornes-bcodina 2015/09 solar eclipses USE module_ra_eclipse , ONLY : solar_eclipse ! This driver calls subroutines for the radiation parameterizations. ! ! short wave radiation choices: ! 1. swrad (19??) ! 4. rrtmg_sw - Added November 2008, MJIacono, AER, Inc. ! ! long wave radiation choices: ! 1. rrtmlwrad ! 4. rrtmg_lw - Added November 2008, MJIacono, AER, Inc. ! !---------------------------------------------------------------------- IMPLICIT NONE ! ! ! Radiation_driver is the WRF mediation layer routine that provides the interface to ! to radiation physics packages in the WRF model layer. The radiation ! physics packages to call are chosen by setting the namelist variable ! (Rconfig entry in Registry) to the integer value assigned to the ! particular package (package entry in Registry). For example, if the ! namelist variable ra_lw_physics is set to 1, this corresponds to the ! Registry Package entry for swradscheme. Note that the Package ! names in the Registry are defined constants (frame/module_state_description.F) ! in the CASE statements in this routine. ! ! Among the arguments is moist, a four-dimensional scalar array storing ! a variable number of moisture tracers, depending on the physics ! configuration for the WRF run, as determined in the namelist. The ! highest numbered index of active moisture tracers the integer argument ! n_moist (note: the number of tracers at run time is the quantity ! n_moist - PARAM_FIRST_SCALAR + 1 , not n_moist. Individual tracers ! may be indexed from moist by the Registry name of the tracer prepended ! with P_; for example P_QC is the index of cloud water. An index ! represents a valid, active field only if the index is greater than ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual ! indices for each tracer is defined in module_state_description and ! set in set_scalar_indices_from_config defined in frame/module_configure.F. ! ! Physics drivers in WRF 2.0 and higher, originally model-layer ! routines, have been promoted to mediation layer routines and they ! contain OpenMP threaded loops over tiles. Thus, physics drivers ! are called from single-threaded regions in the solver. The physics ! routines that are called from the physics drivers are model-layer ! routines and fully tile-callable and thread-safe. ! ! !====================================================================== ! Grid structure in physics part of WRF !---------------------------------------------------------------------- ! The horizontal velocities used in the physics are unstaggered ! relative to temperature/moisture variables. All predicted ! variables are carried at half levels except w, which is at full ! levels. Some arrays with names (*8w) are at w (full) levels. ! !---------------------------------------------------------------------- ! In WRF, kms (smallest number) is the bottom level and kme (largest ! number) is the top level. In your scheme, if 1 is at the top level, ! then you have to reverse the order in the k direction. ! ! kme - half level (no data at this level) ! kme ----- full level ! kme-1 - half level ! kme-1 ----- full level ! . ! . ! . ! kms+2 - half level ! kms+2 ----- full level ! kms+1 - half level ! kms+1 ----- full level ! kms - half level ! kms ----- full level ! !====================================================================== ! Grid structure in physics part of WRF ! !------------------------------------- ! The horizontal velocities used in the physics are unstaggered ! relative to temperature/moisture variables. All predicted ! variables are carried at half levels except w, which is at full ! levels. Some arrays with names (*8w) are at w (full) levels. ! !================================================================== ! Definitions !----------- ! Theta potential temperature (K) ! Qv water vapor mixing ratio (kg/kg) ! Qc cloud water mixing ratio (kg/kg) ! Qr rain water mixing ratio (kg/kg) ! Qi cloud ice mixing ratio (kg/kg) ! Qs snow mixing ratio (kg/kg) !----------------------------------------------------------------- !-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3) !-- PM2_5_WATER PM2.5 water mass (ug m^-3) !-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3) !-- RTHRATEN Theta tendency ! due to radiation (K/s) !-- RTHRATENLW Theta tendency ! due to long wave radiation (K/s) !-- RTHRATENSW Theta temperature tendency ! due to short wave radiation (K/s) !-- dt time step (s) !-- itimestep number of time steps !-- GLW downward long wave flux at ground surface (W/m^2) !-- GSW net short wave flux at ground surface (W/m^2) !-- SWDOWN downward short wave flux at ground surface (W/m^2) !-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ) !-- RLWTOA upward long wave at top of atmosphere (w/m2) !-- RSWTOA upward short wave at top of atmosphere (w/m2) !-- XLAT latitude, south is negative (degree) !-- XLONG longitude, west is negative (degree) !-- ALBEDO albedo (between 0 and 1) !-- CLDFRA cloud fraction (between 0 and 1) !-- CLDFRA_DP cloud fraction from deep cloud in a cumulus scheme !-- CLDFRA_SH cloud fraction from shallow cloud in a cumulus scheme !-- CLDFRA_MP_ALL cloud fraction from CAMMGMP microphysics scheme !-- EMISS surface emissivity (between 0 and 1) !-- rho_phy density (kg/m^3) !-- rr dry air density (kg/m^3) !-- moist moisture array (4D - last index is species) (kg/kg) !-- n_moist number of moisture species !-- qndrop Cloud droplet number (#/kg) !-- p8w pressure at full levels (Pa) !-- p_phy pressure (Pa) !-- Pb base-state pressure (Pa) !-- pi_phy exner function (dimensionless) !-- dz8w dz between full levels (m) !-- t_phy temperature (K) !-- t8w temperature at full levels (K) !-- GMT Greenwich Mean Time Hour of model start (hour) !-- JULDAY the initial day (Julian day) !-- RADT time for calling radiation (min) !-- ra_call_offset -1 (old) means usually just before output, 0 after !-- DEGRAD conversion factor for ! degrees to radians (pi/180.) (rad/deg) !-- DPD degrees per day for earth's ! orbital position (deg/day) !-- R_d gas constant for dry air (J/kg/K) !-- CP heat capacity at constant pressure for dry air (J/kg/K) !-- G acceleration due to gravity (m/s^2) !-- rvovrd R_v divided by R_d (dimensionless) !-- XTIME time since simulation start (min) !-- DECLIN solar declination angle (rad) !-- SOLCON solar constant (W/m^2) !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !-- i_start start indices for i in tile !-- i_end end indices for i in tile !-- j_start start indices for j in tile !-- j_end end indices for j in tile !-- kts start index for k in tile !-- kte end index for k in tile !-- num_tiles number of tiles ! !================================================================== ! LOGICAL, OPTIONAL, INTENT(IN) :: explicit_convection LOGICAL :: expl_conv INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & kts,kte, & num_tiles INTEGER, INTENT(IN) :: lw_physics, sw_physics, mp_physics, sw_eclipse INTEGER, INTENT(IN) :: o3input, aer_opt INTEGER, INTENT(IN) :: id integer, intent(in) :: swint_opt INTEGER, OPTIONAL, INTENT(IN) :: smkfdb_opt, chem_opt ! RAR: INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & i_start,i_end,j_start,j_end INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset INTEGER, INTENT(IN ) :: alevsiz, no_src_types INTEGER, INTENT(IN ) :: levsiz, n_ozmixm INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2 REAL, INTENT(IN ) :: cam_abs_freq_s LOGICAL, INTENT(IN ) :: warm_rain INTEGER, INTENT(IN ) :: cu_physics !CuP, wig 5-Oct-2006 !BSINGH - For WRFCuP scheme !BSINGH - For WRFCuP scheme LOGICAL, OPTIONAL, INTENT(IN) :: shallowcu_forced_ra !CuP, wig !BSINGH -ENDS LOGICAL, INTENT(IN ) :: is_CAMMGMP_used !BSINGH:01/31/2013: Added for CAM5 RRTMG REAL, INTENT(IN ) :: RADT REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN ) :: XLAND, & XICE, & TSK, & VEGFRA, & SNOW !tgs - for solar angle albedo correction INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN) :: IVGTYP INTEGER, INTENT(IN) :: alb_sol REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ALBBCK, & ALBSOL,ALBBCKSOL ! REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, & INTENT(IN ) :: OZMIXM REAL, DIMENSION( ims:ime, alevsiz, jms:jme, no_src_types, n_ozmixm-1 ), OPTIONAL, & INTENT(IN ) :: AERODM REAL, DIMENSION( ims:ime, kms:kme, jms:jme, no_src_types ), OPTIONAL, & INTENT(INOUT) :: AEROD REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, & INTENT(INOUT) :: AODTOT REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ) :: OZFLG REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN REAL, DIMENSION(alevsiz), OPTIONAL, INTENT(IN ) :: PINA REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2 REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, & INTENT(IN ) :: aerosolc_1, aerosolc_2 REAL, DIMENSION(paerlev), OPTIONAL, & INTENT(IN ) :: m_hybi0 REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: HTOP, & HBOT, & HTOPR, & HBOTR, & CUPPT !BSINGH - For WRFCuP scheme REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, & INTENT(INOUT) :: & cutop, & !CuP, wig 1-Oct-2006 cubot, & !CuP, wig 1-Oct-2006 shall !CuP, wig 4-Feb-2008 !BSINGH -ENDS INTEGER, INTENT(IN ) :: julyr ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(IN ) :: dz8w, & z, & p8w, & p, & pi, & t, & t8w, & rho !BSINGH - For WRFCuP scheme REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, & INTENT(INOUT ) :: cldfra_cup !CuP, wig 1-Oct-2006 !BSINGH -ENDS ! Stochastic parameter perturbations INTEGER, INTENT(IN), OPTIONAL :: spp_rad REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN),OPTIONAL ::pattern_spp_rad ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , & INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb gaer300,gaer400,gaer600,gaer999, & ! jcb waer300,waer400,waer600,waer999 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , & INTENT(IN ) :: qc_cu, qi_cu, qc_bl REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , & INTENT(IN ) :: tauaerlw1,tauaerlw2,tauaerlw3,tauaerlw4, & ! czhao tauaerlw5,tauaerlw6,tauaerlw7,tauaerlw8, & ! czhao tauaerlw9,tauaerlw10,tauaerlw11,tauaerlw12, & ! czhao tauaerlw13,tauaerlw14,tauaerlw15,tauaerlw16 INTEGER, INTENT(IN) :: icloud_cu INTEGER, INTENT(IN ), OPTIONAL :: icloud_bl INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback !jdfcz INTEGER, OPTIONAL, INTENT(IN ) :: progn,prescribe INTEGER, OPTIONAL, INTENT(IN ) :: progn ! ! variables for aerosols (only if running with chemistry) ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , & INTENT(IN ) :: pm2_5_dry, & pm2_5_water, & pm2_5_dry_ec ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT) :: RTHRATEN, & RTHRATENLW, & RTHRATENSW ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , & ! INTENT(INOUT) :: SWUP, & ! SWDN, & ! SWUPCLEAR, & ! SWDNCLEAR, & ! LWUP, & ! LWDN, & ! LWUPCLEAR, & ! LWDNCLEAR REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::& ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, & ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, & ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, & ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC ! TOA and surface, upward and downward, total and clear fluxes REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::& SWUPT, SWUPTC, SWDNT, SWDNTC, & SWUPB, SWUPBC, SWDNB, SWDNBC, & LWUPT, LWUPTC, LWDNT, LWDNTC, & LWUPB, LWUPBC, LWDNB, LWDNBC ! Upward and downward, total and clear sky layer fluxes (W m-2) REAL, DIMENSION( ims:ime, kms:kme+2, jms:jme ), & OPTIONAL, INTENT(INOUT) :: & SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC, & LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , & INTENT(INOUT) :: SWCF, & LWCF, & OLR ! ---- fds (06/2010) ssib alb components ------------ REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , & INTENT(IN ) :: ALSWVISDIR, & ALSWVISDIF, & ALSWNIRDIR, & ALSWNIRDIF ! ---- fds (06/2010) ssib swr components ------------ REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , & INTENT(OUT ) :: SWVISDIR, & SWVISDIF, & SWNIRDIR, & SWNIRDIF INTEGER, OPTIONAL, INTENT(IN ) :: sf_surface_physics ! REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN ) :: XLAT, & XLONG, & ALBEDO, & EMISS ! REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: GSW, & GLW REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN ! ------------------------------------------------------------------------------ jararias 2013/08/10 ----------- REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: swddir, & ! All-sky SW broadband surface direct irradiance swddni, & ! All-sky SW broadband surface direct normal irradiance swddif, & ! All-sky SW broadband surface diffuse irradiance swddnic,& ! Clear-sky SW broadband surface direct normal irradiance: J. Kenyon, 30 Sep 2016 swddifc ! Clear-sky SW broadband surface diffuse irradiance: J. Kenyon, 30 Sep 2016 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: Gx,Bx,gg,bb, & ! For SW sza-interpolation swdown_ref, & swddir_ref, & coszen_ref ! ------------------------------------------------------------------------------ jararias 2013/11 ----------- INTEGER, INTENT(IN) :: aer_type, & ! rural, urban, maritime, ... aer_aod550_opt, & ! input option for AOD at 550 nm aer_angexp_opt, & ! input option for aerosol Angstrom exponent aer_ssa_opt, & ! input option for aerosol ssa aer_asy_opt ! input option for aerosol asy REAL, INTENT(IN) :: aer_aod550_val, & ! AOD at 550 nm if aer_aod550_opt=1 aer_angexp_val, & ! aerosol Angstrom exponent if aer_angexp_opt=1 aer_ssa_val, & ! aerosol ssa if aer_ssa_opt=1 aer_asy_val ! aerosol asy if aer_asy_opt=1 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, & INTENT(INOUT) :: aod5502d, & ! gridded AOD at 550 nm from auxinput angexp2d, & ! gridded aerosol Angstrom exponent from auxinput aerssa2d, & ! gridded aerosol ssa from auxinput aerasy2d ! gridded aerosol asy from auxinput REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, & INTENT(OUT) :: aod5503d ! 3D AOD at 550 nm REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL:: taod5503d ! Trude REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL:: taod5502d ! Trude ! REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: AOD3D_SMOKE ! RAR REAL, INTENT(IN ) :: GMT,dt, & julian, xtime INTEGER, INTENT(IN ),OPTIONAL :: YR ! INTEGER, INTENT(IN ) :: JULDAY, itimestep REAL, INTENT(IN ),OPTIONAL :: CURR_SECS LOGICAL, INTENT(IN ),OPTIONAL :: ADAPT_STEP_FLAG INTEGER,INTENT(IN) :: NPHS REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: & CFRACH, & !Added CFRACL, & !Added CFRACM, & !Added CZMEAN !Added REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: & RLWTOA, & !Added RSWTOA, & !Added ACFRST, & !Added ACFRCV !Added INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: & NCFRST, & !Added NCFRCV !Added ! Optional, only for Goddard LW and SW REAL, DIMENSION(IMS:IME, JMS:JME, 1:8) :: ERBE_out !extra output for SDSU REAL, DIMENSION(IMS:IME, JMS:JME), OPTIONAL, INTENT(INOUT) :: & !BSINGH(PNNL)- Lahey compiler forced this specification to be intent-inout TLWDN, TLWUP, & SLWDN, SLWUP, & TSWDN, TSWUP, & SSWDN, SSWUP ! for Goddard schemes ! Added by ZCX for low and total cloud fraction REAL, DIMENSION( kms:kme ), OPTIONAL, INTENT(IN) :: znu ! eta values on half (mass)levels REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) :: & cldt ! Optional (only used by CAM lw scheme) REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,& INTENT(INOUT) :: abstot REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,& INTENT(INOUT) :: absnxt REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,& INTENT(INOUT) :: emstot ! ! Optional ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(INOUT) :: CLDFRA REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & ! ckay for sub-grid cloud fraction OPTIONAL, & INTENT(INOUT) :: cldfra_dp, & cldfra_sh, & cldfra_bl !..G. Thompson REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: re_cloud, re_ice, re_snow INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(IN ) :: & F_ICE_PHY, & F_RAIN_PHY, & CLDFRA_MP_ALL #if (EM_CORE == 1) REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(IN ) :: & LRADIUS, & IRADIUS #endif REAL, DIMENSION( ims:ime, jms:jme ), & OPTIONAL, & INTENT(OUT) :: SWDOWNC ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(INOUT ) :: & pb & ,qv,qc,qr,qi,qs,qg,qndrop, & qnifa,qnwfa ! Trude LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop, & f_qnifa,f_qnwfa ! trude ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(INOUT) :: taucldi,taucldc REAL, OPTIONAL, INTENT(IN) :: dxkm ! Variables for slope-dependent radiation REAL, OPTIONAL, INTENT(IN) :: dx,dy INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,topo_shading REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: ht INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: shadowmask REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(OUT) :: diffuse_frac REAL , OPTIONAL, INTENT(INOUT) :: radtacttime ! Storing the time in s when radiation is called next REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(INOUT) :: o3rad ! vert nesting REAL, OPTIONAL , INTENT(IN ) :: p_top REAL :: p_top_dummy ! LOCAL VAR INTEGER, DIMENSION( ims:ime, kms:kme, jms:jme ) :: cldfra1_flag REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS REAL, DIMENSION( ims:ime, jms:jme ) :: coszr REAL, DIMENSION( ims:ime, levsiz, jms:jme ) :: ozmixt REAL, DIMENSION( ims:ime, alevsiz, jms:jme, 1:no_src_types ) :: aerodt REAL :: DECLIN,SOLCON,XXLAT,TLOCTM,XT24, CEN_LAT, cldfra_cup_mod INTEGER :: i,j,k,its,ite,jts,jte,ij INTEGER :: STEPABS LOGICAL :: gfdl_lw,gfdl_sw, compute_cldfra_cup LOGICAL :: doabsems LOGICAL, EXTERNAL :: wrf_dm_on_monitor INTEGER :: s REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, & DJUL,RJUL,ECCFAC REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qc_temp REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qs_save REAL :: gridkm, Wice,Wh2o REAL :: next_rad_time, DTaccum LOGICAL :: run_param , doing_adapt_dt , decided LOGICAL :: flg_lw, flg_sw !ZCX REAL :: cldji,cldlji !ckay - do not need this - modified code below ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: cldfra_cu !------------------------------------------------------------------ ! solar related variables are added to declaration !------------------------------------------------- REAL, OPTIONAL, INTENT(OUT) :: DECLINX,SOLCONX REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZEN REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: HRANG !------------------------------------------------------------------ ! jararias, 2013/08/10 real :: ioh,kt,airmass,kd real, dimension(ims:ime,jms:jme) :: coszen_loc,hrang_loc ! jararias 2013/11 but modified on 2016/2/10 ! the following three arrays may be dimensioned by (ims,ime,kms,kme,jms,jme,aerosol_vars) real, dimension(:,:,:,:), pointer :: tauaer_sw=>null(), ssaaer_sw=>null(), asyaer_sw=>null() !eclipse variables real, dimension(ims:ime,jms:jme), INTENT(OUT) :: obscur integer, dimension(ims:ime,jms:jme), INTENT(OUT) :: mask real, INTENT(OUT) :: elat_track, elon_track ! Trude AOD variables INTEGER, PARAMETER:: taer_type = 1 ! rural, urban, maritime, ... INTEGER, PARAMETER:: taer_aod550_opt = 2 ! input option for AOD at 550 nm INTEGER, PARAMETER:: taer_angexp_opt = 3 ! input option for aerosol Angstrom exponent INTEGER, PARAMETER:: taer_ssa_opt = 3 ! input option for aerosol ssa INTEGER, PARAMETER:: taer_asy_opt = 3 ! input option for aerosol asy #if ( HWRF == 1 ) CHARACTER(len=255) :: wrf_err_message #endif !tgs--- istwe for albedo correction according to Ivan Csiszar integer :: ilst real :: dm !tgs --- 30 July 2019 --- reduce parameter d from 0.4 to 0.25 for short ! vegetation to get more realistic albedo in morning and evening ! hours. ! real, dimension (2), parameter :: d = (/0.1,0.4/) real, dimension (2), parameter :: d = (/0.1,0.25/) integer, dimension (1:21) :: istwe DATA istwe /5*1,2,2,1,1,5*2,1,2,2,1,2,2,2/ ! MODIS 21 classes ! DATA istwe/9*2,6*1,2,2,1,2,2,1,2,2,1/ ! USGS 24 classes ! This allows radiation schemes (mainly HWRF) to correctly ! interface with the convection scheme when convection is only ! enabled in some domains: if(present(explicit_convection)) then expl_conv=explicit_convection else expl_conv=.true. ! backward compatibility for ARW endif IF ( ICLOUD == 3 ) THEN IF (PRESENT(dxkm)) then gridkm = 1.414*SQRT(dxkm*dxkm + dy*0.001*dy*0.001) ELSE IF (PRESENT(dx)) then gridkm = SQRT(dx*0.001*dx*0.001 + dy*0.001*dy*0.001) endif if (itimestep .LE. 100) then WRITE ( wrf_err_message , * ) 'Grid spacing in km ', dx, dy, gridkm CALL wrf_debug (100, wrf_err_message) endif END IF CALL wrf_debug (1, 'Top of Radiation Driver') ! WRITE ( wrf_err_message , * ) 'itimestep = ',itimestep,', dt = ',dt,', lw and sw options = ',lw_physics,sw_physics ! CALL wrf_debug (1, wrf_err_message ) if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return ! amontornes-bcodina (2014-05-02) :: improving the namelist settings consistency for the FLG scheme ! if (lw_physics .ne. FLGLWSCHEME .and. sw_physics .eq. FLGSWSCHEME) then ! call wrf_error_fatal('SW and LW schemes are in conflict. SW is FLG and LW is a different scheme!') ! end if ! if (lw_physics .eq. FLGLWSCHEME .and. sw_physics .ne. FLGSWSCHEME) then ! call wrf_error_fatal('SW and LW schemes are in conflict. LW is FLG and SW is a different scheme!') ! end if ! ra_call_offset = -1 gives old method where radiation may be called just before output ! ra_call_offset = 0 gives new method where radiation may be called just after output ! and is also consistent with removal of offset in new XTIME ! also need to account for stepra=1 which always has zero modulo output doing_adapt_dt = .FALSE. IF ( PRESENT(adapt_step_flag) ) THEN IF ( adapt_step_flag ) THEN doing_adapt_dt = .TRUE. IF ( radtacttime .eq. 0. ) THEN radtacttime = CURR_SECS + radt*60. END IF END IF END IF ! Do we run through this scheme or not? ! Test 1: If this is the initial model time, then yes. ! ITIMESTEP=1 ! Test 2: If the user asked for the radiation to be run every time step, then yes. ! RADT=0 or STEPRA=1 ! Test 3: If not adaptive dt, and this is on the requested radiation frequency, then yes. ! MOD(ITIMESTEP,STEPRA)=0 (or 1, depending on the offset setting) ! Test 4: If using adaptive dt and the current time is past the last requested activate radiation time, then yes. ! CURR_SECS >= RADTACTTIME ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme. ! We only proceed to other tests if the previous tests all have left decided as FALSE. ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next ! radiation run. run_param = .FALSE. decided = .FALSE. IF ( ( .NOT. decided ) .AND. & ( itimestep .EQ. 1 ) ) THEN run_param = .TRUE. decided = .TRUE. END IF IF ( ( .NOT. decided ) .AND. & ( ( radt .EQ. 0. ) .OR. ( stepra .EQ. 1 ) ) ) THEN run_param = .TRUE. decided = .TRUE. END IF IF ( ( .NOT. decided ) .AND. & ( .NOT. doing_adapt_dt ) .AND. & ( MOD(itimestep,stepra) .EQ. 1+ra_call_offset ) ) THEN run_param = .TRUE. decided = .TRUE. END IF IF ( ( .NOT. decided ) .AND. & ( doing_adapt_dt ) .AND. & ( curr_secs .GE. radtacttime ) ) THEN run_param = .TRUE. decided = .TRUE. radtacttime = curr_secs + radt*60 END IF if(swint_opt.eq.1) then DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, & DEGRAD,DPD ) call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, & julian,xtime,gmt,declin,degrad, & xlong,xlat,coszen_loc,hrang_loc) end do end if Radiation_step: IF ( run_param ) then ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default) STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset & .or. STEPABS .eq. 1 ) THEN doabsems = .true. ELSE doabsems = .false. ENDIF IF (PRESENT(adapt_step_flag)) THEN IF ((adapt_step_flag)) THEN IF ( (itimestep .EQ. 1) .OR. (cam_abs_freq_s .EQ. 0) .OR. & ( CURR_SECS + dt >= ( INT( CURR_SECS / ( cam_abs_freq_s ) + 1 ) * cam_abs_freq_s) ) ) THEN doabsems = .true. ELSE doabsems = .false. ENDIF ENDIF ENDIF gfdl_lw = .false. gfdl_sw = .false. flg_lw = .false. flg_sw = .false. ! Allocate aerosol arrays used by aer_opt = 2 option IF ( PRESENT( AOD5502D ) ) THEN ! jararias, 2013/11 IF ( aer_opt .EQ. 2 ) THEN swrad_aerosol_select: select case(sw_physics) case(GODDARDSWSCHEME) allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:11)) allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:11)) allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:11)) case(RRTMG_SWSCHEME,RRTMG_SWSCHEME_FAST) allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) end select swrad_aerosol_select ENDIF ENDIF ! Allocate aerosol arrays used by aer_opt = 3 option, and explicit AOD from QNWFA+QNIFA (Trude) IF (PRESENT(f_qnwfa) .AND. PRESENT(f_qnifa) .AND. PRESENT(taod5503d) .AND. PRESENT(taod5502d)) THEN IF (F_QNWFA .AND. aer_opt.eq.3 .AND. & (sw_physics.eq.RRTMG_SWSCHEME .OR. & sw_physics.eq.RRTMG_SWSCHEME_FAST)) THEN CALL wrf_debug (150, 'DEBUG-GT: computing 3D AOD from QNWFA+QNIFA') allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:14)) !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte) DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) do j=jts,jte do i=its,ite taod5502d(i,j) = 0.0 end do end do call gt_aod (p, DZ8W, t, qv, qnwfa, qnifa, taod5503d, & ims,ime, jms,jme, kms,kme,its,ite, jts,jte, kts,kte) ! RAR: #if ( WRF_SMOKE == 1 ) IF (chem_opt==18 .AND. smkfdb_opt==1) THEN do j=jts,jte do i=its,ite do k=kts,kte taod5503d(i,k,j)= taod5503d(i,k,j) + MIN(3.0,AOD3D_SMOKE(i,k,j)) enddo enddo enddo ENDIF #endif do j=jts,jte do i=its,ite do k=kts,kte taod5502d(i,j) = taod5502d(i,j) + taod5503d(i,k,j) end do end do end do ENDDO !$OMP END PARALLEL DO ENDIF ENDIF !--------------- ! Calculate constant for short wave radiation ! moved up and out of OMP loop because it only needs to be computed once ! and because it is not entirely thread-safe (XT24, TOLOCTM and XXLAT need ! their thread-privacy) JM 20100217 DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, & DEGRAD,DPD ) ! amontornes-bcodina 2015/09 solar eclipses ! Solar eclipse prediction based on the Bessel's method ! outputs: obscur, mask, elat_track, elon_track CALL solar_eclipse(ims,ime,jms,jme,its,ite,jts,jte, & julian,gmt,YR,xtime,radt, & degrad,xlong,xlat,obscur,mask, & elat_track,elon_track,sw_eclipse ) IF(PRESENT(declinx).AND.PRESENT(solconx))THEN ! saved to state arrays used in surface driver declinx=declin solconx=solcon ENDIF ! added coszen subroutine : jararias, 2013/08/10 ! outputs: coszen, hrang call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, & julian,xtime+radt*0.5,gmt, & declin,degrad,xlong,xlat,coszen,hrang) !tgs - 3 april 2018 ! normalize instantaneous surface albedos (as) for overhead sun (an) ! using the parameterization from Dickinson(1983), an=as/dm do j=jts,jte do i=its,ite if(itimestep == 1 ) then albsol(i,j) = albedo(i,j) albbcksol(i,j) = albbck(i,j) endif if(alb_sol == 1) then if(coszen(i,j) > 0.) then ! no albedo normalizing for water, snow or ice if(xland(i,j) < 1.5 .and. snow(i,j)==0. .and. xice(i,j)==0.) then ilst=istwe(ivgtyp(i,j)) ! 1 or 2 dm = (1.+2.*d(ilst))/(1.+2.*d(ilst)*coszen(i,j)) albsol(i,j) = albbck(i,j)*dm albbcksol(i,j) = albbck(i,j)*dm endif endif ! coszen > 0. endif ! alb_sol enddo enddo !tgs - end DO j=jts,jte DO i=its,ite albsol(i,j) = MIN(albsol(i,j),0.9) ENDDO ENDDO ENDDO !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte) DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) ! initialize data if ((itimestep.eq.1).and.(swint_opt.eq.1)) then do j=jts,jte do i=its,ite Bx(i,j)=0. bb(i,j)=0. Gx(i,j)=0. gg(i,j)=0. end do end do end if DO j=jts,jte DO i=its,ite GSW(I,J)=0. GLW(I,J)=0. SWDOWN(I,J)=0. swddir(i,j)=0. ! jararias, 2013/08/10 swddni(i,j)=0. ! jararias, 2013/08/10 swddif(i,j)=0. ! jararias, 2013/08/10 GLAT(I,J)=XLAT(I,J)*DEGRAD GLON(I,J)=XLONG(I,J)*DEGRAD ENDDO ENDDO DO j=jts,jte DO k=kts,kte+1 DO i=its,ite RTHRATEN(I,K,J)=0. RTHRATENLW(I,K,J)=0. RTHRATENSW(I,K,J)=0. ! SWUP(I,K,J) = 0.0 ! SWDN(I,K,J) = 0.0 ! SWUPCLEAR(I,K,J) = 0.0 ! SWDNCLEAR(I,K,J) = 0.0 ! LWUP(I,K,J) = 0.0 ! LWDN(I,K,J) = 0.0 ! LWUPCLEAR(I,K,J) = 0.0 ! LWDNCLEAR(I,K,J) = 0.0 CEMISS(I,K,J)=0.0 ENDDO ENDDO ENDDO IF ( PRESENT( SWUPFLX ) ) THEN DO j=jts,jte DO k=kts,kte+2 DO i=its,ite SWUPFLX(I,K,J) = 0.0 SWDNFLX(I,K,J) = 0.0 SWUPFLXC(I,K,J) = 0.0 SWDNFLXC(I,K,J) = 0.0 LWUPFLX(I,K,J) = 0.0 LWDNFLX(I,K,J) = 0.0 LWUPFLXC(I,K,J) = 0.0 LWDNFLXC(I,K,J) = 0.0 ENDDO ENDDO ENDDO ENDIF ! backup the incoming hydrometeors IF ( F_QC ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qc_save(i,k,j) = qc(i,k,j) ENDDO ENDDO ENDDO ENDIF IF ( F_QI ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qi_save(i,k,j) = qi(i,k,j) ENDDO ENDDO ENDDO ENDIF ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006) if( F_QC ) then DO j=jts,jte DO k=kts,kte DO i=its,ite qc_temp(I,K,J)=qc(I,K,J) ENDDO ENDDO ENDDO else DO j=jts,jte DO k=kts,kte DO i=its,ite qc_temp(I,K,J)=0. ENDDO ENDDO ENDDO endif ! Remove this - to match NAM operational (affects GFDL and HWRF schemes) ! if( F_QR ) then ! DO j=jts,jte ! DO k=kts,kte ! DO i=its,ite ! qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J) ! ENDDO ! ENDDO ! ENDDO ! endif ! ! temporarily modify hydrometeors (this is for GD scheme and WRF-Chem) ! IF ( F_QC .AND. icloud_cu .EQ. 1 ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qc(i,k,j) = qc(i,k,j) + qc_cu(i,k,j) ENDDO ENDDO ENDDO ENDIF IF ( F_QI .AND. icloud_cu .EQ. 1 ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qi(i,k,j) = qi(i,k,j) + qi_cu(i,k,j) ENDDO ENDDO ENDDO ENDIF ! Choose how to compute cloud fraction (since 3.6) ! Initialize to zero DO j=jts,jte DO k=kts,kte DO i=its,ite CLDFRA(i,k,j) = 0. END DO END DO END DO IF ( ICLOUD == 1 ) THEN IF ( F_QC .AND. F_QI ) THEN ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998) CALL wrf_debug (1, 'CALL cldfra1') CALL cal_cldfra1(CLDFRA,qv,qc,qi,qs, & F_QV,F_QC,F_QI,F_QS,t,p, & F_ICE_PHY,F_RAIN_PHY, & mp_physics,cldfra1_flag, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IF ( PRESENT ( CLDFRA_DP ) ) THEN ! this is for Kain-Fritsch scheme IF ( icloud_cu .EQ. 2 ) THEN CALL wrf_debug (1, 'use kf cldfra') DO j = jts,jte DO k = kts,kte DO i = its,ite !JOE-changed to not need cldfra_cu (save memory) !cldfra_cu(i,k,j)=cldfra_dp(i,k,j)+cldfra_sh(i,k,j) ! Cu cloud fraction !CLDFRA(i,k,j)=(1.-cldfra_cu(i,k,j))*CLDFRA(i,k,j) ! Update resolved cloud fraction for Cu punch-through CLDFRA(i,k,j)=(1.-(cldfra_dp(i,k,j)+cldfra_sh(i,k,j)))*CLDFRA(i,k,j) !CLDFRA(i,k,j)=CLDFRA(i,k,j)+cldfra_cu(i,k,j) ! New total cloud fraction CLDFRA(i,k,j)=CLDFRA(i,k,j)+(cldfra_dp(i,k,j)+cldfra_sh(i,k,j)) ! New total cloud fraction CLDFRA(i,k,j)=AMIN1(1.0,CLDFRA(i,k,j)) !qc(i,k,j) = qc(i,k,j)+qc_cu(i,k,j)*cldfra_cu(i,k,j) qc(i,k,j) = qc(i,k,j)+qc_cu(i,k,j)*(cldfra_dp(i,k,j)+cldfra_sh(i,k,j)) !qi(i,k,j) = qi(i,k,j)+qi_cu(i,k,j)*cldfra_cu(i,k,j) qi(i,k,j) = qi(i,k,j)+qi_cu(i,k,j)*(cldfra_dp(i,k,j)+cldfra_sh(i,k,j)) ENDDO ENDDO ENDDO ENDIF ENDIF IF ( PRESENT ( CLDFRA_BL ) .AND. PRESENT ( QC_BL ) ) THEN IF ( icloud_bl > 0 ) THEN CALL wrf_debug (1, 'in rad driver; use BL clouds') !tgs - 17sept19 - at the first time step, PBL hasn't ! been called yet and CLDFRA_BL = 0. Therefore, replace CLDFRA ! with CLDFRA_BL starting from the 2-nd time step. !tgs - 20sept19 - added a check for itimestep=1, which ensures ! that Xu-Randall cloud fraction is used only if CLDFRA_BL=0. IF(itimestep == 1) THEN DO j = jts,jte DO i = its,ite DO k = kts,kte if( CLDFRA_BL(i,k,j)>0.001) THEN CLDFRA(i,k,j)=CLDFRA_BL(i,k,j) endif ENDDO ENDDO ENDDO ELSE DO j = jts,jte DO i = its,ite DO k = kts,kte CLDFRA(i,k,j)=CLDFRA_BL(i,k,j) ENDDO ENDDO ENDDO ENDIF DO j = jts,jte DO i = its,ite DO k = kts,kte IF (qc(i,k,j) < 1.E-6 .AND. qi(i,k,j) < 1.E-8 .AND. CLDFRA_BL(i,k,j)>0.001) THEN !Partition the BL clouds into water & ice according to a linear !approximation of Hobbs et al. (1974). This allows us to only use !one 3D array for both cloud water & ice. ! Wice = 1. - MIN(1., MAX(0., (t(i,k,j)-254.)/15.)) ! Wh2o = 1. - Wice !CLDFRA(i,k,j)=MAX(CLDFRA(i,k,j),CLDFRA_BL(i,k,j)) !CLDFRA(i,k,j)=MAX(0.0,MIN(1.0,CLDFRA(i,k,j))) qc(i,k,j)=qc(i,k,j) + QC_BL(i,k,j)*(MIN(1., MAX(0., (t(i,k,j)-254.)/15.)))*CLDFRA_BL(i,k,j) qi(i,k,j)=qi(i,k,j) + QC_BL(i,k,j)*(1. - MIN(1., MAX(0., (t(i,k,j)-254.)/15.)))*CLDFRA_BL(i,k,j) ENDIF ENDDO ENDDO ENDDO ENDIF ENDIF IF ( PRESENT (cldfra_mp_all) ) THEN IF (is_CAMMGMP_used) THEN !BSINGH: cloud fraction from CAMMGMP is being used (Mods by Po-Lun) CALL wrf_debug (1, 'use cammgmp') IF (itimestep .NE. 1) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite CLDFRA(i,k,j) = CLDFRA_MP_ALL(I,K,J) !PMA if (CLDFRA(i,k,j) .lt. 0.01) CLDFRA(i,k,j) = 0. ENDDO ENDDO ENDDO ENDIF ENDIF ENDIF ENDIF ELSE IF ( ICLOUD == 2 ) THEN IF ( F_QC .AND. F_QI ) THEN CALL wrf_debug (1, 'CALL cldfra2') CALL cal_cldfra2(CLDFRA,qc,qi,F_QC,F_QI, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ENDIF !+---+-----------------------------------------------------------------+ !..New cloud fraction scheme added by G. Thompson (2014Oct31) !+---+-----------------------------------------------------------------+ ELSEIF (ICLOUD == 3) THEN IF ( F_QC .AND. F_QI ) THEN IF ( F_QS ) THEN DO j = jts,jte DO k = kts,kte DO i = its,ite qs_save(i,k,j) = qs(i,k,j) ENDDO ENDDO ENDDO ENDIF CALL wrf_debug (150, 'DEBUG: using gthompsn cloud fraction scheme') CALL cal_cldfra3(CLDFRA, qv, qc, qi, qs, & & p,t,rho, XLAND, gridkm, & & ids,ide, jds,jde, kds,kde, & & ims,ime, jms,jme, kms,kme, & & its,ite, jts,jte, kts,kte) ELSE CALL wrf_error_fatal('Can not use icloud = 3 option, missing QC or QI field.') ENDIF END IF ! ww: Interpolating climatological ozone and aerosol to model time and levels ! Adapted from camrad code IF ( PRESENT( O3RAD ) ) THEN call wrf_debug(1,'Have o3rad') #if (EM_CORE==1) IF ( o3input .EQ. 2 .AND. id .EQ. 1 ) THEN #else IF ( o3input .EQ. 2 ) THEN #endif ! ! Find the current month (adapted from Cavallo) ! CALL cam_time_interp( ozmixm, pin, levsiz, date_str, & ! ids , ide , jds , jde , kds , kde , & ! ims , ime , jms , jme , kms , kme , & ! its , ite , jts , jte , kts , kte ) ! this is the CAM version ! interpolate to model time-step call ozn_time_int(julday,julian,ozmixm,ozmixt,levsiz,n_ozmixm, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! interpolate to model model levels call ozn_p_int(p ,pin, levsiz, ozmixt, o3rad, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ENDIF ENDIF IF ( PRESENT( AEROD ) ) THEN IF ( aer_opt .EQ. 1 .AND. id .EQ. 1 ) THEN call aer_time_int(julday,julian,aerodm,aerodt,alevsiz,n_ozmixm-1,no_src_types, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) call aer_p_int(p ,pina, alevsiz, aerodt, aerod, no_src_types, p8w, AODTOT, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ENDIF ENDIF !Modify CLDFRA and QC for kfcupscheme cumulus scheme if(present(cldfra_cup)) then !BSINGH - overwrite cldfra with the cloud fraction computed in module_cu_kfcup.F !Force cloud fraction if cumulus triggered. if( cu_physics == KFCUPSCHEME ) then do j = jts,jte do k = kts,kte do i = its,ite !Find whether to overwrite cldfra or not (ONLY if ICLOUD == 1) compute_cldfra_cup = .true. if (icloud == 1 ) then compute_cldfra_cup = .false. if(cldfra1_flag(i,k,j) == 1 .and. shall(i,j) .gt. 1) then CLDFRA(i,k,j)=0. elseif(cldfra1_flag(i,k,j) == 2 .and. shall(i,j) .gt. 1) then CLDFRA(i,k,j)=1. elseif(cldfra1_flag(i,k,j) == 3) then compute_cldfra_cup = .true. endif endif if(compute_cldfra_cup) then if( (int(shall(i,j)) .le.1) .and. k >= int(cubot(i,j)) .and. k <= int(cutop(i,j)) ) then ! LD Mar232012 CLDFRA(i,k,j) = cldfra_cup(i,k,j) else if( shall(i,j) .gt. 1) then !!LD cldfra_cup(i,k,j) = 0.0 end if endif if( shall(i,j) == 1 .and. k >= cubot(i,j) .and. k <= cutop(i,j) ) then ! 1=Shallow Cu ! Begin: wig, 4-Feb-2008 ! ! Override the cloud condensate values if shallow convection triggered. ! For shallow convection, use a representative condensate value based on ! observations from CHAPS (Oklahoma area) and Florida (Blyth et al. 2005) ! or the predicted value if it is greater. cldfra_cup_mod = cldfra_cup(i,k,j)* 1.0e-3*(1+qv(i,k,j))!modified cloud fraction qc(i,k,j) = max( 1.0e-3*cldfra_cup_mod, qc(i,k,j) )!DE+LB 2012-Feb ! Override the cloud fraction values calculated above if shallow ! convection triggered. For shallow convection, use a representative ! cloud fraction. In this case, the median value for shallow Cu cases ! at the ARM SGP site, 36%, Berg et al. 2008, J. Clim. if( shallowcu_forced_ra )cldfra(i,k,j) = max(0.36, cldfra(i,k,j)) ! Median shallow Cu frac at ARM SGP endif ENDDO ENDDO ENDDO end if endif lwrad_select: SELECT CASE(lw_physics) CASE (RRTMSCHEME) CALL wrf_debug (100, 'CALL rrtm') IF ( PRESENT(p_top) ) THEN p_top_dummy = p_top ELSE p_top_dummy = -1. ! not used by NMM END IF CALL RRTMLWRAD( & P_TOP=p_top_dummy & ,RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS & ,QV3D=QV & ,QC3D=QC & ,QR3D=QR & ,QI3D=QI & ,QS3D=QS & ,QG3D=QG & ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t & ,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G & ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR & ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG & ,ICLOUD=icloud,WARM_RAIN=warm_rain & !ccc Added for time-varying trace gases. ,YR=YR,JULIAN=JULIAN & !ccc ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ) CASE (goddardlwscheme) CALL wrf_debug (100, 'CALL goddard longwave radiation scheme ') IF (itimestep.eq.1) then call wrf_message('running goddard lw radiation') ENDIF CALL goddardrad(sw_or_lw='lw' & ,rthraten=rthraten,gsf=glw,xlat=xlat,xlong=xlong & ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi & ,dz8w=dz8w,rho_phy=rho,emiss=emiss & ,cldfra3d=cldfra & ,gmt=gmt,cp=cp,g=g,t8w=t8w & ,julday=julday,xtime=xtime & ,declin=declin,solcon=solcon & , center_lat = cen_lat & ,radfrq=radt,degrad=degrad & ,taucldi=taucldi,taucldc=taucldc & ,warm_rain=warm_rain & ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte & ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban ,qv3d=qv & ,qc3d=qc & ,qr3d=qr & ,qi3d=qi & ,qs3d=qs & ,qg3d=qg & ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr & ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg & ,erbe_out=erbe_out & !optional ,aer_opt=aer_opt & ,tauaer3d_sw=tauaer_sw & ! jararias, 2013/11 ,ssaaer3d_sw=ssaaer_sw & ! jararias, 2013/11 ,asyaer3d_sw=asyaer_sw & ! jararias, 2012/11 ) CASE (GFDLLWSCHEME) CALL wrf_debug (100, 'CALL gfdllw') IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. & PRESENT(F_QS) .AND. PRESENT(qs) .AND. & PRESENT(qv) .AND. PRESENT(qc) ) THEN IF ( F_QV .AND. F_QC .AND. F_QS) THEN gfdl_lw = .true. CALL ETARA( & DT=dt,XLAND=xland & ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t & ,QV=qv,QW=qc_temp,QI=qi,QS=qs & ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW & ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi & ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot & ,HBOTR=hbotr, HTOPR=htopr & ,ALBEDO=albedo,CUPPT=cuppt & ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt & ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep & ,XTIME=xtime,JULIAN=julian & ,JULYR=julyr,JULDAY=julday & ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw & ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach & ,ACFRST=acfrst,NCFRST=ncfrst & ,ACFRCV=acfrcv,NCFRCV=ncfrcv & ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean & ,THRATEN=rthraten,THRATENLW=rthratenlw & ,THRATENSW=rthratensw & ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ) ELSE CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.') ENDIF ELSE CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.') ENDIF #if ( HWRF == 1 ) CASE (HWRFLWSCHEME) CALL wrf_debug (100, 'CALL hwrflw') gfdl_lw = .true. CALL HWRFRA(explicit_convection=expl_conv, & DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, & XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, & QV=qv,QW=qc_temp,QI=Qi, & TSK2D=tsk,GLW=GLW,GSW=GSW, & TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt,& VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, & CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, & ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, & its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte ) #endif CASE (CAMLWSCHEME) CALL wrf_debug(100, 'CALL camrad lw') IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. & PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. & PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) & .AND. PRESENT(AEROSOLC_2).AND. PRESENT(ALSWVISDIR) ) THEN CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, & dolw=.true.,dosw=.false., & SWUPT=SWUPT,SWUPTC=SWUPTC, & SWDNT=SWDNT,SWDNTC=SWDNTC, & LWUPT=LWUPT,LWUPTC=LWUPTC, & LWDNT=LWDNT,LWDNTC=LWDNTC, & SWUPB=SWUPB,SWUPBC=SWUPBC, & SWDNB=SWDNB,SWDNBC=SWDNBC, & LWUPB=LWUPB,LWUPBC=LWUPBC, & LWDNB=LWDNB,LWDNBC=LWDNBC, & SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, & TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, & GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, & ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS & ,QV3D=qv & ,QC3D=qc & ,QR3D=qr & ,QI3D=qi & ,QS3D=qs & ,QG3D=qg & ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010) ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010) ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010) ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010) ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni & !amontornes-bcodina (2014-04-20) ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr & ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg & ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy & ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, & dz8w=dz8w, & CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, & ozmixm=ozmixm,pin0=pin,levsiz=levsiz, & num_months=n_ozmixm, & m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, & aerosolcn=aerosolc_2,m_hybi0=m_hybi0, & paerlev=paerlev, naer_c=n_aerosolc, & cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, & GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, & SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 & ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot & ,doabsems=doabsems & ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ,coszen=coszen ) ELSE CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' ) ENDIF CASE (RRTMG_LWSCHEME) CALL wrf_debug (100, 'CALL rrtmg_lw') CALL RRTMG_LWRAD( & RTHRATENLW=RTHRATEN, & LWUPT=LWUPT,LWUPTC=LWUPTC, & LWDNT=LWDNT,LWDNTC=LWDNTC, & LWUPB=LWUPB,LWUPBC=LWUPBC, & LWDNB=LWDNB,LWDNBC=LWDNBC, & GLW=GLW,OLR=RLWTOA,LWCF=LWCF, & EMISS=EMISS, & P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, & T8W=t8w,RHO3D=rho,R=R_d,G=G, & ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,& #if (EM_CORE == 1) LRADIUS=lradius, IRADIUS=iradius, & #endif IS_CAMMGMP_USED=is_cammgmp_used, & !ckay ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,& F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, & XLAND=XLAND,XICE=XICE,SNOW=SNOW, & QV3D=QV,QC3D=QC,QR3D=QR, & QI3D=QI,QS3D=QS,QG3D=QG, & O3INPUT=O3INPUT,O33D=O3RAD, & F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, & F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, & RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson #if ( WRF_CHEM == 1 ) TAUAERLW1=tauaerlw1,TAUAERLW2=tauaerlw2, & ! jcb TAUAERLW3=tauaerlw3,TAUAERLW4=tauaerlw4, & ! jcb TAUAERLW5=tauaerlw5,TAUAERLW6=tauaerlw6, & ! jcb TAUAERLW7=tauaerlw7,TAUAERLW8=tauaerlw8, & ! jcb TAUAERLW9=tauaerlw9,TAUAERLW10=tauaerlw10, & ! jcb TAUAERLW11=tauaerlw11,TAUAERLW12=tauaerlw12, & ! jcb TAUAERLW13=tauaerlw13,TAUAERLW14=tauaerlw14, & ! jcb TAUAERLW15=tauaerlw15,TAUAERLW16=tauaerlw16, & ! jcb aer_ra_feedback=aer_ra_feedback, & !jdfcz progn=progn,prescribe=prescribe, & progn=progn, & #endif QNDROP3D=qndrop,F_QNDROP=f_qndrop, & !ccc Added for time-varying trace gases. YR=YR,JULIAN=JULIAN, & !ccc IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,& IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,& ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,& LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, & LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC, & mp_physics=mp_physics ) CASE (RRTMG_LWSCHEME_FAST) CALL wrf_debug (100, 'CALL rrtmg_lw') CALL RRTMG_LWRAD_FAST( & RTHRATENLW=RTHRATEN, & LWUPT=LWUPT,LWUPTC=LWUPTC, & LWDNT=LWDNT,LWDNTC=LWDNTC, & LWUPB=LWUPB,LWUPBC=LWUPBC, & LWDNB=LWDNB,LWDNBC=LWDNBC, & GLW=GLW,OLR=RLWTOA,LWCF=LWCF, & EMISS=EMISS, & P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, & T8W=t8w,RHO3D=rho,R=R_d,G=G, & ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,& #if (EM_CORE == 1) LRADIUS=lradius, IRADIUS=iradius, & #endif IS_CAMMGMP_USED=is_cammgmp_used, & !ckay ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,& F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, & XLAND=XLAND,XICE=XICE,SNOW=SNOW, & QV3D=QV,QC3D=QC,QR3D=QR, & QI3D=QI,QS3D=QS,QG3D=QG, & O3INPUT=O3INPUT,O33D=O3RAD, & F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, & F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, & RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson #if ( WRF_CHEM == 1 ) TAUAERLW1=tauaerlw1,TAUAERLW2=tauaerlw2, & ! jcb TAUAERLW3=tauaerlw3,TAUAERLW4=tauaerlw4, & ! jcb TAUAERLW5=tauaerlw5,TAUAERLW6=tauaerlw6, & ! jcb TAUAERLW7=tauaerlw7,TAUAERLW8=tauaerlw8, & ! jcb TAUAERLW9=tauaerlw9,TAUAERLW10=tauaerlw10, & ! jcb TAUAERLW11=tauaerlw11,TAUAERLW12=tauaerlw12, & ! jcb TAUAERLW13=tauaerlw13,TAUAERLW14=tauaerlw14, & ! jcb TAUAERLW15=tauaerlw15,TAUAERLW16=tauaerlw16, & ! jcb aer_ra_feedback=aer_ra_feedback, & !jdfcz progn=progn,prescribe=prescribe, & progn=progn, & #endif QNDROP3D=qndrop,F_QNDROP=f_qndrop, & !ccc Added for time-varying trace gases. YR=YR,JULIAN=JULIAN, & !ccc IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,& IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,& ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,& LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, & LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC & ) CASE (HELDSUAREZ) CALL wrf_debug (100, 'CALL heldsuarez') CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, & t8w, rho, R_d,G,CP, dt, xlat, degrad, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! -- add by Jin Kong 10/2011 CASE (FLGLWSCHEME) CALL wrf_debug (100, 'CALL Fu-Liou-Gu') flg_lw = .true. !test Jin Kong 11/01/2011 for ozone ozflg = 0. !test over CALL RAD_FLG( & peven=p8w,podd=p,t8w=t8w,degrees=t & ,pi3d=pi,o3=ozflg & ,G=G,CP=CP & ,albedo=ALBEDO,tskin=tsk & ,h2o=QV,cld_iccld=QI,cld_wlcld=QC & ,cld_prwc=QR,cld_pgwc=QG,cld_snow=QS & ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR & ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG & ,warm_rain=warm_rain & ,cloudstrf=CLDFRA & ,emiss=EMISS & ,air_den=rho & ,dz3d=dz8w & ,SOLCON=SOLCON & ,declin=DECLIN & ,xtime=xtime, xlong=xlong, xlat=xlat & ,JULDAY=JULDAY & ,gmt=gmt, radt=radt, degrad=degrad & ,dtcolumn=dt & ,ids=ids,ide=ide,jds=jds,jde=jde & ,kds=kds,kde=kde & ,ims=ims,idim=ime,jms=jms,jdim=jme & ,kms=kms,kmax=kme & ,its=its,ite=ite,jts=jts,jte=jte & ,kts=kts,kte=kte & ,uswtop=RSWTOA,ulwtop=RLWTOA & ,NETSWBOT=GSW,DLWBOT=GLW & ,DSWBOT=SWDOWN,deltat=RTHRATEN & ,dtshort=RTHRATENSW,dtlongwv=RTHRATENLW & ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni & ! amontornes-bcodina (2014-04-20) ) CALL wrf_debug(100, 'a4 Fu_Liou-Gu') ! -- end CASE DEFAULT WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics CALL wrf_error_fatal ( wrf_err_message ) END SELECT lwrad_select IF (lw_physics .gt. 0 .and. .not.gfdl_lw .and. .not.flg_lw) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite RTHRATENLW(I,K,J)=RTHRATEN(I,K,J) ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY) IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J) ENDDO ENDDO ENDDO ENDIF IF (lw_physics .eq. goddardlwscheme) THEN IF ( PRESENT (tlwdn) ) THEN DO j=jts,jte DO i=its,ite tlwdn(i,j) = erbe_out(i,j,1) ! TOA LW downwelling flux [W/m2] tlwup(i,j) = erbe_out(i,j,2) ! TOA LW upwelling flux [W/m2] slwdn(i,j) = erbe_out(i,j,3) ! surface LW downwelling flux [W/m2] slwup(i,j) = erbe_out(i,j,4) ! surface LW upwelling flux [W/m2] ENDDO ENDDO ENDIF ENDIF IF ( PRESENT( AOD5502D ) ) THEN ! jararias, 2013/11 IF ( aer_opt .EQ. 2 ) THEN swrad_aerosol_select2: select case(sw_physics) case(GODDARDSWSCHEME) call wrf_debug(100, 'call calc_aerosol_goddard_sw') call calc_aerosol_goddard_sw(ht,dz8w,p,t,qv,aer_type,aer_aod550_opt,aer_angexp_opt, & aer_ssa_opt,aer_asy_opt,aer_aod550_val,aer_angexp_val, & aer_ssa_val,aer_asy_val,aod5502d,angexp2d,aerssa2d, & aerasy2d,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte, & tauaer_sw,ssaaer_sw,asyaer_sw ) do j=jts,jte do i=its,ite do k=kts,kte aod5503d(i,k,j)=tauaer_sw(i,k,j,8) ! band at 550 nm end do end do end do case(RRTMG_SWSCHEME,RRTMG_SWSCHEME_FAST) call wrf_debug(100, 'call calc_aerosol_rrtmg_sw') call calc_aerosol_rrtmg_sw(ht,dz8w,p,t,qv,aer_type,aer_aod550_opt,aer_angexp_opt, & aer_ssa_opt,aer_asy_opt,aer_aod550_val,aer_angexp_val, & aer_ssa_val,aer_asy_val,aod5502d,angexp2d,aerssa2d, & aerasy2d,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte, & tauaer_sw,ssaaer_sw,asyaer_sw ) do j=jts,jte do i=its,ite do k=kts,kte aod5503d(i,k,j)=tauaer_sw(i,k,j,10) ! band at 550 nm end do end do end do case default end select swrad_aerosol_select2 ENDIF ENDIF !..Different treatment for aer_opt=3 using QNWFA+QNIFA aerosol species (Trude) IF (PRESENT(f_qnwfa) .AND. PRESENT(f_qnifa)) THEN IF (F_QNWFA .AND. aer_opt.eq.3 .AND. & (sw_physics.eq.RRTMG_SWSCHEME .OR. & sw_physics.eq.RRTMG_SWSCHEME_FAST)) THEN call wrf_debug(100, 'call calc_aerosol_rrtmg_sw with 3D AOD values') call calc_aerosol_rrtmg_sw(ht,dz8w,p,t,qv,taer_type,taer_aod550_opt,taer_angexp_opt, & taer_ssa_opt,taer_asy_opt,aer_aod550_val,aer_angexp_val, & aer_ssa_val,aer_asy_val,taod5502d,angexp2d,aerssa2d, & aerasy2d,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte, & tauaer_sw,ssaaer_sw,asyaer_sw, taod5503d) ENDIF ENDIF swrad_select: SELECT CASE(sw_physics) CASE (SWRADSCHEME) CALL wrf_debug(100, 'CALL swrad') CALL SWRAD( & DT=dt,RTHRATEN=rthraten,GSW=gsw & ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo & #if ( WRF_CHEM == 1 ) ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water & ,PM2_5_DRY_EC=pm2_5_dry_ec & #endif ,RHO_PHY=rho,T3D=t & ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt & ,R=r_d,CP=cp,G=g,JULDAY=julday & ,XTIME=xtime,DECLIN=declin,SOLCON=solcon & ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad & ,warm_rain=warm_rain & ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ,QV3D=qv & ,QC3D=qc & ,QR3D=qr & ,QI3D=qi & ,QS3D=qs & ,QG3D=qg & ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr & ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg & ,coszen=coszen,julian=julian ) CASE (GSFCSWSCHEME) CALL wrf_debug(100, 'CALL gsfcswrad') CALL GSFCSWRAD( & RTHRATEN=rthraten,GSW=gsw & ! PAJ: xlat and xlong removed. ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi & ,DZ8W=dz8w,RHO_PHY=rho & ,CLDFRA3D=cldfra,RSWTOA=rswtoa & ,CP=cp,G=g & ! PAJ: GMT removed. ,JULDAY=julday & ! PAJ: XTIME removed. ,SOLCON=solcon & ! PAJ: declin removed ! ,RADFRQ=radt,DEGRAD=degrad & ! PAJ: degrad and radfrq removed ,TAUCLDI=taucldi,TAUCLDC=taucldc & ,WARM_RAIN=warm_rain & #if ( WRF_CHEM == 1 ) ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb ,GAER300=gaer300,GAER400=gaer400 & ! jcb ,GAER600=gaer600,GAER999=gaer999 & ! jcb ,WAER300=waer300,WAER400=waer400 & ! jcb ,WAER600=waer600,WAER999=waer999 & ! jcb ,aer_ra_feedback=aer_ra_feedback & #endif ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ,QV3D=qv & ,QC3D=qc & ,QR3D=qr & ,QI3D=qi & ,QS3D=qs & ,QG3D=qg & ,QNDROP3D=qndrop & ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr & ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg & ,F_QNDROP=f_qndrop & ,COSZEN=coszen & ) CASE (goddardswscheme) CALL wrf_debug(100, 'CALL goddard shortwave radiation scheme ') IF (itimestep.eq.1) then call wrf_message('running goddard sw radiation') ENDIF CALL goddardrad(sw_or_lw='sw' & ,rthraten=rthraten,gsf=gsw,xlat=xlat,xlong=xlong & ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi & ,dz8w=dz8w,rho_phy=rho,emiss=emiss & ,cldfra3d=cldfra & ,gmt=gmt,cp=cp,g=g,t8w=t8w & ,julday=julday,xtime=xtime & ,declin=declin,solcon=solcon & ,center_lat = cen_lat & ,radfrq=radt,degrad=degrad & ,taucldi=taucldi,taucldc=taucldc & ,warm_rain=warm_rain & ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde & ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme & ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte & ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban ,qv3d=qv & ,qc3d=qc & ,qr3d=qr & ,qi3d=qi & ,qs3d=qs & ,qg3d=qg & ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr & ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg & ,erbe_out=erbe_out & !optional ,swddir=swddir,swddni=swddni,swddif=swddif & ! jararias, 14/08/2013 ,coszen=coszen,julian=julian & ! jararias, 14/08/2013 ,tauaer3d_sw=tauaer_sw & ! jararias, 2013/11 ,ssaaer3d_sw=ssaaer_sw & ! jararias, 2013/11 ,asyaer3d_sw=asyaer_sw & ! jararias, 2012/11 ,aer_opt=aer_opt & ) CASE (CAMSWSCHEME) CALL wrf_debug(100, 'CALL camrad sw') IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. & PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. & PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) & .AND. PRESENT(AEROSOLC_2) .AND. PRESENT(ALSWVISDIR)) THEN CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, & dolw=.false.,dosw=.true., & SWUPT=SWUPT,SWUPTC=SWUPTC, & SWDNT=SWDNT,SWDNTC=SWDNTC, & LWUPT=LWUPT,LWUPTC=LWUPTC, & LWDNT=LWDNT,LWDNTC=LWDNTC, & SWUPB=SWUPB,SWUPBC=SWUPBC, & SWDNB=SWDNB,SWDNBC=SWDNBC, & LWUPB=LWUPB,LWUPBC=LWUPBC, & LWDNB=LWDNB,LWDNBC=LWDNBC, & SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, & TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, & GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, & ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS & ,QV3D=qv & ,QC3D=qc & ,QR3D=qr & ,QI3D=qi & ,QS3D=qs & ,QG3D=qg & ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010) ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010) ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010) ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010) ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni & !amontornes-bcodina (2014-04-20) ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr & ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg & ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy & ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, & dz8w=dz8w, & CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, & ozmixm=ozmixm,pin0=pin,levsiz=levsiz, & num_months=n_ozmixm, & m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, & aerosolcn=aerosolc_2,m_hybi0=m_hybi0, & paerlev=paerlev, naer_c=n_aerosolc, & cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, & GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, & SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 & ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot & ,doabsems=doabsems & ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ,coszen=coszen ) ELSE CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' ) ENDIF DO j=jts,jte DO k=kts,kte DO i=its,ite RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J) ENDDO ENDDO ENDDO CASE (RRTMG_SWSCHEME) CALL wrf_debug(100, 'CALL rrtmg_sw') CALL RRTMG_SWRAD( & RTHRATENSW=RTHRATENSW, & SWUPT=SWUPT,SWUPTC=SWUPTC, & SWDNT=SWDNT,SWDNTC=SWDNTC, & SWUPB=SWUPB,SWUPBC=SWUPBC, & SWDNB=SWDNB,SWDNBC=SWDNBC, & SWCF=SWCF,GSW=GSW, & XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, & RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, & COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, & ALBEDO=ALBSOL,t3d=t,t8w=t8w,TSK=TSK, & !tgs ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, & p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, & dz8w=dz8w,CLDFRA3D=CLDFRA, & #if (EM_CORE == 1) LRADIUS=lradius, IRADIUS=iradius, & #endif IS_CAMMGMP_USED=is_cammgmp_used, & R=R_D,G=G, & !ckay ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,& ICLOUD=icloud,WARM_RAIN=warm_rain, & F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, & XLAND=XLAND,XICE=XICE,SNOW=SNOW, & QV3D=qv,QC3D=qc,QR3D=qr, & QI3D=qi,QS3D=qs,QG3D=qg, & O3INPUT=O3INPUT,O33D=O3RAD, & AER_OPT=AER_OPT,aerod=aerod,no_src=no_src_types, & ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif, & !Zhenxin ssib alb comp (06/2010) ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif, & !Zhenxin ssib alb comp (06/2010) SWVISDIR=swvisdir ,SWVISDIF=swvisdif, & !Zhenxin ssib swr comp (06/2010) SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif, & !Zhenxin ssib swr comp (06/2010) SF_SURFACE_PHYSICS=sf_surface_physics, & !Zhenxin ssib sw_phy (06/2010) F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, & F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, & RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson #if ( WRF_CHEM == 1 ) TAUAER300=tauaer300,TAUAER400=tauaer400, & ! jcb TAUAER600=tauaer600,TAUAER999=tauaer999, & ! jcb GAER300=gaer300,GAER400=gaer400, & ! jcb GAER600=gaer600,GAER999=gaer999, & ! jcb WAER300=waer300,WAER400=waer400, & ! jcb WAER600=waer600,WAER999=waer999, & ! jcb aer_ra_feedback=aer_ra_feedback, & !jdfcz progn=progn,prescribe=prescribe, & progn=progn, & #endif QNDROP3D=qndrop,F_QNDROP=f_qndrop, & IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,& IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,& ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,& SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, & SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC, & tauaer3d_sw=tauaer_sw, & ! jararias 2013/11 ssaaer3d_sw=ssaaer_sw, & ! jararias 2013/11 asyaer3d_sw=asyaer_sw, & ! jararias 2013/11 swddir=swddir,swddni=swddni,swddif=swddif, & ! jararias 2013/08/10 swddnic=swddnic,swddifc=swddifc, & ! J. Kenyon, 30 Sep 2016 obscur=obscur, & ! amontornes-bcodina 2015/09 solar eclipses xcoszen=coszen,julian=julian,mp_physics=mp_physics,& ! jararias 2013/08/14 spp_rad=spp_rad,pattern_spp_rad=pattern_spp_rad ) DO j=jts,jte DO k=kts,kte DO i=its,ite RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J) ENDDO ENDDO ENDDO CASE (RRTMG_SWSCHEME_FAST) CALL wrf_debug(100, 'CALL rrtmg_sw_fast') CALL RRTMG_SWRAD_FAST( & RTHRATENSW=RTHRATENSW, & SWUPT=SWUPT,SWUPTC=SWUPTC, & SWDNT=SWDNT,SWDNTC=SWDNTC, & SWUPB=SWUPB,SWUPBC=SWUPBC, & SWDNB=SWDNB,SWDNBC=SWDNBC, & SWCF=SWCF,GSW=GSW, & XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, & RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, & COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, & ALBEDO=ALBSOL,t3d=t,t8w=t8w,TSK=TSK, & !tgs ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, & p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, & dz8w=dz8w,CLDFRA3D=CLDFRA, & #if (EM_CORE == 1) LRADIUS=lradius, IRADIUS=iradius, & #endif IS_CAMMGMP_USED=is_cammgmp_used, & R=R_D,G=G, & !ckay ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,& ICLOUD=icloud,WARM_RAIN=warm_rain, & F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, & XLAND=XLAND,XICE=XICE,SNOW=SNOW, & QV3D=qv,QC3D=qc,QR3D=qr, & QI3D=qi,QS3D=qs,QG3D=qg, & O3INPUT=O3INPUT,O33D=O3RAD, & AER_OPT=AER_OPT,aerod=aerod,no_src=no_src_types, & ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif, & !Zhenxin ssib alb comp (06/2010) ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif, & !Zhenxin ssib alb comp (06/2010) SWVISDIR=swvisdir ,SWVISDIF=swvisdif, & !Zhenxin ssib swr comp (06/2010) SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif, & !Zhenxin ssib swr comp (06/2010) SF_SURFACE_PHYSICS=sf_surface_physics, & !Zhenxin ssib sw_phy (06/2010) F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, & F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, & RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson #if ( WRF_CHEM == 1 ) TAUAER300=tauaer300,TAUAER400=tauaer400, & ! jcb TAUAER600=tauaer600,TAUAER999=tauaer999, & ! jcb GAER300=gaer300,GAER400=gaer400, & ! jcb GAER600=gaer600,GAER999=gaer999, & ! jcb WAER300=waer300,WAER400=waer400, & ! jcb WAER600=waer600,WAER999=waer999, & ! jcb aer_ra_feedback=aer_ra_feedback, & !jdfcz progn=progn,prescribe=prescribe, & progn=progn, & #endif QNDROP3D=qndrop,F_QNDROP=f_qndrop, & IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,& IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,& ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,& SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, & SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC, & tauaer3d_sw=tauaer_sw, & ! jararias 2013/11 ssaaer3d_sw=ssaaer_sw, & ! jararias 2013/11 asyaer3d_sw=asyaer_sw, & ! jararias 2013/11 swddir=swddir,swddni=swddni,swddif=swddif, & ! jararias 2013/08/10 xcoszen=coszen,julian=julian ) ! jararias 2013/08/14 DO j=jts,jte DO k=kts,kte DO i=its,ite RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J) ENDDO ENDDO ENDDO CASE (GFDLSWSCHEME) CALL wrf_debug (100, 'CALL gfdlsw') IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. & PRESENT(F_QS) .AND. PRESENT(qs) .AND. & PRESENT(qv) .AND. PRESENT(qc) ) THEN IF ( F_QV .AND. F_QC .AND. F_QS ) THEN gfdl_sw = .true. CALL ETARA( & DT=dt,XLAND=xland & ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t & ,QV=qv,QW=qc_temp,QI=qi,QS=qs & ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW & ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi & ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot & ,HBOTR=hbotr, HTOPR=htopr & ,ALBEDO=albedo,CUPPT=cuppt & ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt & ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep & ,XTIME=xtime,JULIAN=julian & ,JULYR=julyr,JULDAY=julday & ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw & ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach & ,ACFRST=acfrst,NCFRST=ncfrst & ,ACFRCV=acfrcv,NCFRCV=ncfrcv & ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean & ,THRATEN=rthraten,THRATENLW=rthratenlw & ,THRATENSW=rthratensw & ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde & ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme & ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte & ) ELSE CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.') ENDIF ELSE CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.') ENDIF #if ( HWRF == 1 ) CASE (HWRFSWSCHEME) CALL wrf_debug (100, 'CALL hwrfsw') gfdl_sw = .true. CALL HWRFRA(explicit_convection=expl_conv, & DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, & XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, & QV=qv,QW=qc_temp,QI=Qi, & TSK2D=tsk,GLW=GLW,GSW=GSW, & TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt, & VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, & CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, & ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, & its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte ) #endif CASE (0) ! Here in case we don't want to call a sw radiation scheme ! For example, the Held-Suarez idealized test case IF (lw_physics /= HELDSUAREZ) THEN WRITE( wrf_err_message , * ) & 'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')' CALL wrf_error_fatal ( wrf_err_message ) END IF ! -- add by Jin Kong 10/2011 !--- the following FLGSWSCHEME is for testing only CASE (FLGSWSCHEME) flg_sw = .true. !--- No need to do anything since the short and long wave is calculted in one program ! -- end CASE DEFAULT WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics CALL wrf_error_fatal ( wrf_err_message ) END SELECT swrad_select IF (sw_physics .eq. goddardswscheme) THEN IF ( PRESENT (tswdn) ) THEN DO j=jts,jte DO i=its,ite tswdn(i,j) = erbe_out(i,j,5) ! TOA SW downwelling flux [W/m2] tswup(i,j) = erbe_out(i,j,6) ! TOA SW upwelling flux [W/m2] sswdn(i,j) = erbe_out(i,j,7) ! surface SW downwelling flux [W/m2] sswup(i,j) = erbe_out(i,j,8) ! surface SW upwelling flux [W/m2] ENDDO ENDDO ENDIF ENDIF IF (sw_physics .gt. 0 .and. .not.gfdl_sw .and. .not.flg_sw) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J) ENDDO ENDDO ENDDO DO j=jts,jte DO i=its,ite SWDOWN(I,J)=GSW(I,J)/(1.-ALBSOL(I,J)) !tgs SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J)) ENDDO ENDDO ENDIF ! jararias, 14/08/2013 ! surface direct and diffuse SW fluxes computation. Only for schemes other than RRTMG and Goddard ! Backup method in case sw scheme in use does not provide surface SW direct and diffuse irradiances IF ((sw_physics .NE. RRTMG_SWSCHEME) .AND. (sw_physics .NE. RRTMG_SWSCHEME_FAST) & .AND. (sw_physics .NE. FLGSWSCHEME) .AND. (sw_physics .NE. CAMSWSCHEME) & ! amontornes-bcodina (2014-04-20) .AND. (sw_physics .ne. GODDARDSWSCHEME)) THEN DO j=jts,jte DO i=its,ite IF (coszen(i,j).GT.1e-3) THEN ioh=solcon*coszen(i,j) ! TOA irradiance kt=swdown(i,j)/max(ioh,1e-3) ! clearness index ! Optical air mass: Rigollier et al. (2000) doi: ! 10.1016/S0038-092X(99)00055-9 airmass=exp(-ht(i,j)/8434.5)/(coszen(i,j)+ & 0.50572*(asin(coszen(i,j))*57.295779513082323+6.07995)**(-1.6364)) ! kt correction for air-mass at large sza: Perez et al. (1990) ! doi: 10.1016/0038-092X(90)90036-C kt=kt/(0.1+1.031*exp(-1.4/(0.9+(9.4/max(airmass,1e-3))))) ! Diffuse fraction: Ruiz-Arias et al. (2010) (Eq 33) doi: ! 10.1016/j.enconman.2009.11.024 kd=0.952-1.041*exp(-exp(2.300-4.702*kt)) swddif(i,j)=kd*swdown(i,j) swddir(i,j)=(1.-kd)*swdown(i,j) swddni(i,j)=swddir(i,j)/max(coszen(i,j),1e-4) ENDIF ENDDO ENDDO ENDIF IF ( PRESENT( diffuse_frac ) ) THEN DO j=jts,jte DO i=its,ite if (swdown(i,j).gt.0.001) then diffuse_frac(i,j) = swddif(i,j)/swdown(i,j) diffuse_frac(i,j) = min(diffuse_frac(i,j),1.0) else diffuse_frac(i,j) = 0. endif ENDDO ENDDO ENDIF ! Restore qc & qi for any model physics configuration IF ( F_QC ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qc(i,k,j) = qc_save(i,k,j) ENDDO ENDDO ENDDO ENDIF IF ( F_QI ) THEN DO j=jts,jte DO k=kts,kte DO i=its,ite qi(i,k,j) = qi_save(i,k,j) ENDDO ENDDO ENDDO ENDIF IF (ICLOUD == 3 .AND. F_QS ) THEN DO j = jts,jte DO k = kts,kte DO i = its,ite qs(i,k,j) = qs_save(i,k,j) ENDDO ENDDO ENDDO ENDIF ! jararias, aug 2013, updated 2013/11 ! parameters update for SW surface fluxes interpolation IF (swint_opt.EQ.1) THEN ! interpolation applies on all-sky fluxes (swddir, swdown) CALL update_swinterp_parameters(ims,ime,jms,jme,its,ite,jts,jte, & coszen,coszen_loc,swddir,swdown, & swddir_ref,bb,Bx,swdown_ref,gg,Gx, & coszen_ref ) ENDIF ENDDO !$OMP END PARALLEL DO IF ( associated(tauaer_sw) ) deallocate(tauaer_sw) IF ( associated(ssaaer_sw) ) deallocate(ssaaer_sw) IF ( associated(asyaer_sw) ) deallocate(asyaer_sw) ENDIF Radiation_step ! jararias, aug 2013 ! SW surface fluxes interpolation (meaningful when not in a Radiation_step) if (swint_opt .eq. 1) then call wrf_debug(100,'SW surface irradiance interpolation') !--------------- !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte) do ij = 1,num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) call interp_sw_radiation(ims,ime,jms,jme,its,ite,jts,jte, & coszen_ref,coszen_loc,swddir_ref, & bb,Bx,swdown_ref,gg,Gx,albsol, & ! bb,Bx,swdown_ref,gg,Gx,albedo, & swdown,swddir,swddni,swddif,gsw ) enddo !$OMP END PARALLEL DO end if accumulate_lw_select: SELECT CASE(lw_physics) CASE (CAMLWSCHEME,RRTMG_LWSCHEME,RRTMG_LWSCHEME_FAST) IF(PRESENT(LWUPTC))THEN ! NMM calls the driver every RADT time steps, EM calls every DT #if (EM_CORE == 1) DTaccum = DT #else DTaccum = RADT*60 #endif !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte) DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) DO j=jts,jte DO i=its,ite ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DTaccum ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DTaccum ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DTaccum ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DTaccum ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DTaccum ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DTaccum ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DTaccum ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DTaccum ENDDO ENDDO ENDDO !$OMP END PARALLEL DO ENDIF CASE DEFAULT END SELECT accumulate_lw_select accumulate_sw_select: SELECT CASE(sw_physics) CASE (CAMSWSCHEME,RRTMG_SWSCHEME,RRTMG_SWSCHEME_FAST) IF(PRESENT(SWUPTC))THEN ! NMM calls the driver every RADT time steps, EM calls every DT #if (EM_CORE == 1) DTaccum = DT #else DTaccum = RADT*60 #endif !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte) DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) DO j=jts,jte DO i=its,ite ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DTaccum ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DTaccum ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DTaccum ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DTaccum ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DTaccum ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DTaccum ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DTaccum ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DTaccum ENDDO ENDDO ENDDO !$OMP END PARALLEL DO ENDIF CASE DEFAULT END SELECT accumulate_sw_select ! compute cloud diagnosis (random overlapping) ! by ZCX IF ( PRESENT ( CLDFRA ) .AND. PRESENT ( CLDT ) .AND. & PRESENT ( F_QC ) .AND. PRESENT ( F_QI ) ) THEN DO ij = 1 , num_tiles its = i_start(ij) ite = i_end(ij) jts = j_start(ij) jte = j_end(ij) DO j=jts,jte DO i=its,ite cldji=1.0 do k=kte-1,kts,-1 cldji=cldji*(1.0-cldfra(i,k,j)) enddo cldt(i,j)=1.0-cldji ! cldlji=1.0 ! do k=kte-1,kts,-1 ! if(znu(k).ge.0.69) then ! cldlji=cldlji*(1.0-cldfra(i,k,j)) ! endif ! enddo ! cldl(i,j)=1.0-cldlji END DO END DO END DO END IF END SUBROUTINE radiation_driver SUBROUTINE pre_radiation_driver ( grid, config_flags & ,itimestep, ra_call_offset & ,XLAT, XLONG, GMT, julian, xtime, RADT, STEPRA & ,ht,dx,dy,sina,cosa,shadowmask,slope_rad ,topo_shading & ,shadlen,ht_shad,ht_loc & ,ht_shad_bxs, ht_shad_bxe & ,ht_shad_bys, ht_shad_bye & ,nested, min_ptchsz & ,spec_bdy_width & ,ids, ide, jds, jde, kds, kde & ,ims, ime, jms, jme, kms, kme & ,ips, ipe, jps, jpe, kps, kpe & ,i_start, i_end & ,j_start, j_end & ,kts, kte & ,num_tiles ) USE module_domain , ONLY : domain #ifdef DM_PARALLEL USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks,wrf_dm_minval_integer # if (EM_CORE == 1) USE module_comm_dm , ONLY : halo_toposhad_sub # endif #endif USE module_bc USE module_model_constants IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & kts,kte, & num_tiles TYPE(domain) , INTENT(INOUT) :: grid TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: itimestep, ra_call_offset, stepra, & slope_rad, topo_shading, & spec_bdy_width INTEGER, INTENT(INOUT) :: min_ptchsz INTEGER, DIMENSION(num_tiles), INTENT(IN) :: & i_start,i_end,j_start,j_end REAL, INTENT(IN ) :: GMT, radt, julian, xtime, dx, dy, shadlen REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN ) :: XLAT, & XLONG, & HT, & SINA, & COSA REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width ), & INTENT(IN ) :: ht_shad_bxs, ht_shad_bxe REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width ), & INTENT(IN ) :: ht_shad_bys, ht_shad_bye INTEGER, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: shadowmask LOGICAL, INTENT(IN ) :: nested !Local ! For orographic shading INTEGER :: niter,ni,psx,psy,idum,jdum,i,j,ij REAL :: DECLIN,SOLCON ! Determine minimum patch size for slope-dependent radiation if (itimestep .eq. 1) then psx = ipe-ips+1 psy = jpe-jps+1 min_ptchsz = min(psx,psy) idum = 0 jdum = 0 endif # ifdef DM_PARALLEL if (itimestep .eq. 1) then call wrf_dm_minval_integer (psx,idum,jdum) call wrf_dm_minval_integer (psy,idum,jdum) min_ptchsz = min(psx,psy) endif # endif ! Topographic shading if ((topo_shading.eq.1).and.(itimestep .eq. 1 .or. & mod(itimestep,STEPRA) .eq. 1 + ra_call_offset)) then !--------------- ! Calculate constants for short wave radiation CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,DEGRAD,DPD) ! Make a local copy of terrain height field do j=jms,jme do i=ims,ime ht_loc(i,j) = ht(i,j) enddo enddo ! Determine if iterations are necessary for shadows to propagate from one patch to another if ((ids.eq.ips).and.(ide.eq.ipe).and.(jds.eq.jps).and.(jde.eq.jpe)) then niter = 1 else niter = int(shadlen/(dx*min_ptchsz)+3) endif IF( nested ) THEN !$OMP PARALLEL DO & !$OMP PRIVATE ( ij ) DO ij = 1 , num_tiles CALL spec_bdyfield(ht_shad, & ht_shad_bxs, ht_shad_bxe, & ht_shad_bys, ht_shad_bye, & 'm', config_flags, spec_bdy_width, 2,& ids,ide, jds,jde, 1 ,1 , & ! domain dims ims,ime, jms,jme, 1 ,1 , & ! memory dims ips,ipe, jps,jpe, 1 ,1 , & ! patch dims i_start(ij), i_end(ij), & j_start(ij), j_end(ij), & 1 , 1 ) ENDDO ENDIF do ni = 1, niter !$OMP PARALLEL DO & !$OMP PRIVATE ( ij,i,j ) do ij = 1 , num_tiles call toposhad_init (ht_shad,ht_loc, & shadowmask,nested,ni, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, & i_start(ij),min(i_end(ij), ide-1),j_start(ij),& min(j_end(ij), jde-1), kts, kte ) enddo !$OMP END PARALLEL DO !$OMP PARALLEL DO & !$OMP PRIVATE ( ij,i,j ) do ij = 1 , num_tiles call toposhad (xlat,xlong,sina,cosa,xtime,gmt,radt,declin, & dx,dy,ht_shad,ht_loc,ni, & shadowmask,shadlen, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, & i_start(ij),min(i_end(ij), ide-1),j_start(ij),& min(j_end(ij), jde-1), kts, kte ) enddo !$OMP END PARALLEL DO #if defined( DM_PARALLEL ) && (EM_CORE == 1) # include "HALO_TOPOSHAD.inc" #endif enddo endif END SUBROUTINE pre_radiation_driver !--------------------------------------------------------------------- !BOP ! !IROUTINE: radconst - compute radiation terms ! !INTERFAC: SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, & DEGRAD,DPD ) !--------------------------------------------------------------------- USE module_wrf_error IMPLICIT NONE !--------------------------------------------------------------------- ! !ARGUMENTS: REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN REAL, INTENT(OUT ) :: DECLIN,SOLCON REAL :: OBECL,SINOB,SXLONG,ARG, & DECDEG,DJUL,RJUL,ECCFAC ! ! !DESCRIPTION: ! Compute terms used in radiation physics !EOP ! for short wave radiation DECLIN=0. SOLCON=0. !-----OBECL : OBLIQUITY = 23.5 DEGREE. OBECL=23.5*DEGRAD SINOB=SIN(OBECL) !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX: IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.) IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.) SXLONG=SXLONG*DEGRAD ARG=SINOB*SIN(SXLONG) DECLIN=ASIN(ARG) DECDEG=DECLIN/DEGRAD !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976) DJUL=JULIAN*360./365. RJUL=DJUL*DEGRAD ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* & COS(2*RJUL)+0.000077*SIN(2*RJUL) SOLCON=1370.*ECCFAC END SUBROUTINE radconst SUBROUTINE calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, & julian,xtime,gmt, & declin,degrad,xlon,xlat,coszen,hrang) ! Added Equation of Time correction : jararias, 2013/08/10 implicit none integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte real, intent(in) :: julian,declin,xtime,gmt,degrad real, dimension(ims:ime,jms:jme), intent(in) :: xlat,xlon real, dimension(ims:ime,jms:jme), intent(inout) :: coszen,hrang integer :: i,j real :: da,eot,xt24,tloctm,xxlat da=6.2831853071795862*(julian-1)/365. eot=(0.000075+0.001868*cos(da)-0.032077*sin(da) & -0.014615*cos(2*da)-0.04089*sin(2*da))*(229.18) xt24=mod(xtime,1440.)+eot do j=jts,jte do i=its,ite tloctm=gmt+xt24/60.+xlon(i,j)/15. hrang(i,j)=15.*(tloctm-12.)*degrad xxlat=xlat(i,j)*degrad coszen(i,j)=sin(xxlat)*sin(declin) & +cos(xxlat)*cos(declin) *cos(hrang(i,j)) enddo enddo END SUBROUTINE calc_coszen subroutine update_swinterp_parameters(ims,ime,jms,jme,its,ite,jts,jte, & coszen,coszen_loc,swddir,swdown, & swddir_ref,bb,Bx, & swdown_ref,gg,Gx, & coszen_ref ) ! Author: jararias 2013/11 implicit None integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte real, dimension(ims:ime,jms:jme), intent(in) :: coszen,coszen_loc,swddir,swdown real, dimension(ims:ime,jms:jme), intent(inout) :: swddir_ref,bb,Bx, & swdown_ref,gg,Gx, & coszen_ref integer :: i,j real :: swddir_0,swdown_0,coszen_0 real, parameter :: coszen_min=1e-4 do j=jts,jte do i=its,ite if ((coszen(i,j).gt.coszen_min) .and. (coszen_loc(i,j).gt.coszen_min)) then ! parameters update for DIR if (Bx(i,j).le.0) then swddir_0 =(coszen_loc(i,j)/coszen(i,j))*swddir(i,j) ! linear first guess estimation coszen_0 =coszen_loc(i,j) else swddir_0 =swddir_ref(i,j) coszen_0 =coszen_ref(i,j) end if if ((coszen(i,j)/coszen_0).lt.1.) then bb(i,j) =log(max(1.,swddir(i,j))/max(1.,swddir_0)) / log(min(1.-1e-4,coszen(i,j)/coszen_0)) elseif ((coszen(i,j)/coszen_0).gt.1) then bb(i,j) =log(max(1.,swddir(i,j))/max(1.,swddir_0)) / log(max(1.+1e-4,coszen(i,j)/coszen_0)) else bb(i,j) =0. end if bb(i,j) =max(-.5,min(2.5,bb(i,j))) Bx(i,j) =swddir(i,j)/(coszen(i,j)**bb(i,j)) !write(wrf_err_message,*) 'XXX I=',i,' J=',j,' Bx=',Bx(i,j),' bb=',bb(i,j),' swddir=',swddir(i,j), & ! ' swddir_0=',swddir_0,' coszen=',coszen(i,j),' coszen_0=',coszen_0 !call wrf_debug(1,wrf_err_message) ! parameters update for GHI if (Gx(i,j).le.0) then swdown_0 =(coszen_loc(i,j)/coszen(i,j))*swdown(i,j) ! linear first guess estimation coszen_0 =coszen_loc(i,j) else swdown_0 =swdown_ref(i,j) coszen_0 =coszen_ref(i,j) end if if ((coszen(i,j)/coszen_0).lt.1.) then gg(i,j) =log(max(1.,swdown(i,j))/max(1.,swdown_0)) / log(min(1.-1e-4,coszen(i,j)/coszen_0)) elseif ((coszen(i,j)/coszen_0).gt.1) then gg(i,j) =log(max(1.,swdown(i,j))/max(1.,swdown_0)) / log(max(1.+1e-4,coszen(i,j)/coszen_0)) else gg(i,j) =0. end if gg(i,j) =max(-.5,min(2.5,gg(i,j))) Gx(i,j) =swdown(i,j)/(coszen(i,j)**gg(i,j)) else Bx(i,j) =0. bb(i,j) =0. Gx(i,j) =0. gg(i,j) =0. end if ! saving last SW run in state variables coszen_ref(i,j) =coszen(i,j) swdown_ref(i,j) =swdown(i,j) swddir_ref(i,j) =swddir(i,j) !if ((i.eq.20).and.(j.eq.20)) then ! write(wrf_err_message,'(" RADSTEP : tn=",I4," csz_0=",F9.6," csz=",F9.6," csz_1=",F9.6," Gx=",F14.2," gg=",F9.5, & ! " Bx=",F14.2," bb=",F9.5)') itimestep,coszen_0,coszen_loc(i,j),coszen(i,j),Gx(i,j),gg(i,j), & ! Bx(i,j),bb(i,j) ! call wrf_debug(1,wrf_err_message) !end if end do end do end subroutine update_swinterp_parameters subroutine interp_sw_radiation(ims,ime,jms,jme,its,ite,jts,jte, & coszen_ref,coszen_loc,swddir_ref, & bb,Bx,swdown_ref,gg,Gx,albedo, & swdown,swddir,swddni,swddif,gsw ) ! Author: jararias 2013/11 implicit None integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte real, dimension(ims:ime,jms:jme), intent(in) :: coszen_ref,coszen_loc, & swddir_ref,Bx,bb, & swdown_ref,Gx,gg, & albedo real, dimension(ims:ime,jms:jme), intent(inout) :: swddir,swdown, & swddif,swddni, gsw integer :: i,j real, parameter :: coszen_min=1e-4 do j=jts,jte do i=its,ite ! sza interpolation of surface fluxes if ((coszen_ref(i,j).gt.coszen_min) .and. (coszen_loc(i,j).gt.coszen_min)) then if ((bb(i,j).eq.-0.5).or.(bb(i,j).eq.2.5)) then swddir(i,j) =(coszen_loc(i,j)/coszen_ref(i,j))*swddir_ref(i,j) else swddir(i,j) =Bx(i,j)*(coszen_loc(i,j)**bb(i,j)) end if if ((gg(i,j).eq.-0.5).or.(gg(i,j).eq.2.5)) then swdown(i,j) =(coszen_loc(i,j)/coszen_ref(i,j))*swdown_ref(i,j) else swdown(i,j) =Gx(i,j)*(coszen_loc(i,j)**gg(i,j)) end if swddif(i,j) =swdown(i,j)-swddir(i,j) swddni(i,j) =swddir(i,j)/coszen_loc(i,j) gsw(i,j) =swdown(i,j)*(1.-albedo(i,j)) else swddir(i,j) =0. swdown(i,j) =0. swddif(i,j) =0. swddni(i,j) =0. gsw(i,j) =0. end if end do end do end subroutine interp_sw_radiation !--------------------------------------------------------------------- !BOP ! !IROUTINE: cal_cldfra2 - Compute cloud fraction ! !INTERFACE: SUBROUTINE cal_cldfra2(CLDFRA,QC,QI,F_QC,F_QI, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_state_description, ONLY : KFCUPSCHEME, KFETASCHEME !CuP, wig 5-Oct-2006 !BSINGH - For WRFCuP scheme !--------------------------------------------------------------------- IMPLICIT NONE !--------------------------------------------------------------------- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: & CLDFRA REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: & QI, & QC LOGICAL,INTENT(IN) :: F_QC,F_QI REAL thresh INTEGER:: i,j,k ! !DESCRIPTION: ! Compute cloud fraction from input ice and cloud water fields ! if provided. ! ! Whether QI or QC is active or not is determined from the indices of ! the fields into the 4D scalar arrays in WRF. These indices are ! P_QI and P_QC, respectively, and they are passed in to the routine ! to enable testing to see if QI and QC represent active fields in ! the moisture 4D scalar array carried by WRF. ! ! If a field is active its index will have a value greater than or ! equal to PARAM_FIRST_SCALAR, which is also an input argument to ! this routine. !EOP !--------------------------------------------------------------------- thresh=1.0e-6 IF ( f_qi .AND. f_qc ) THEN DO j = jts,jte DO k = kts,kte DO i = its,ite IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN CLDFRA(i,k,j)=1. ELSE CLDFRA(i,k,j)=0. ENDIF ENDDO ENDDO ENDDO ELSE IF ( f_qc ) THEN DO j = jts,jte DO k = kts,kte DO i = its,ite IF ( QC(i,k,j) .gt. thresh) THEN CLDFRA(i,k,j)=1. ELSE CLDFRA(i,k,j)=0. ENDIF ENDDO ENDDO ENDDO ELSE DO j = jts,jte DO k = kts,kte DO i = its,ite CLDFRA(i,k,j)=0. ENDDO ENDDO ENDDO ENDIF END SUBROUTINE cal_cldfra2 !BOP ! !IROUTINE: cal_cldfra1 - Compute cloud fraction ! !INTERFACE: ! cal_cldfra_xr - Compute cloud fraction. ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done !! !!--- Cloud fraction parameterization follows Xu and Randall (JAS), 1996 !! (see Hong et al., 1998) !! (modified by Ferrier, Feb '02) ! SUBROUTINE cal_cldfra1(CLDFRA, QV, QC, QI, QS, & F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, & F_ICE_PHY,F_RAIN_PHY, & mp_physics, cldfra1_flag, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_state_description, ONLY : KFCUPSCHEME, KFETASCHEME !wig, CuP 4-Fb-2008 !BSINGH - For WRFCuP scheme #if (HWRF == 1) USE module_state_description, ONLY : FER_MP_HIRES, FER_MP_HIRES_ADVECT, ETAMP_HWRF #else USE module_state_description, ONLY : FER_MP_HIRES, FER_MP_HIRES_ADVECT #endif !--------------------------------------------------------------------- IMPLICIT NONE !--------------------------------------------------------------------- INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ! INTEGER, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: cldfra1_flag REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: & CLDFRA REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: & QV, & QI, & QC, & QS, & t_phy, & p_phy ! p_phy, & ! F_ICE_PHY, & ! F_RAIN_PHY REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & OPTIONAL, & INTENT(IN ) :: & F_ICE_PHY, & F_RAIN_PHY LOGICAL,OPTIONAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS INTEGER :: mp_physics ! REAL thresh INTEGER:: i,j,k REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, & PEXP=0.25, RHGRID=1.0 REAL , PARAMETER :: SVP1=0.61078 REAL , PARAMETER :: SVP2=17.2693882 REAL , PARAMETER :: SVPI2=21.8745584 REAL , PARAMETER :: SVP3=35.86 REAL , PARAMETER :: SVPI3=7.66 REAL , PARAMETER :: SVPT0=273.15 REAL , PARAMETER :: r_d = 287. REAL , PARAMETER :: r_v = 461.6 REAL , PARAMETER :: ep_2=r_d/r_v ! !DESCRIPTION: ! Compute cloud fraction from input ice and cloud water fields ! if provided. ! ! Whether QI or QC is active or not is determined from the indices of ! the fields into the 4D scalar arrays in WRF. These indices are ! P_QI and P_QC, respectively, and they are passed in to the routine ! to enable testing to see if QI and QC represent active fields in ! the moisture 4D scalar array carried by WRF. ! ! If a field is active its index will have a value greater than or ! equal to PARAM_FIRST_SCALAR, which is also an input argument to ! this routine. !EOP !----------------------------------------------------------------------- !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION ! (modified by Ferrier, Feb '02) ! !--- Cloud fraction parameterization follows Randall, 1994 ! (see Hong et al., 1998) !----------------------------------------------------------------------- ! Note: ep_2=287./461.6 Rd/Rv ! Note: R_D=287. ! Alternative calculation for critical RH for grid saturation ! RHGRID=0.90+.08*((100.-DX)/95.)**.5 ! Calculate saturation mixing ratio weighted according to the fractions of ! water and ice. ! Following: ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ] ! ! over ice over water ! a = 21.8745584 17.2693882 ! b = 7.66 35.86 !--------------------------------------------------------------------- DO j = jts,jte DO k = kts,kte DO i = its,ite tc = t_phy(i,k,j) - SVPT0 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) ) esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) ) QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw ) QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi ) ifouter: IF ( PRESENT(F_QI) .and. PRESENT(F_QC) .and. PRESENT(F_QS) ) THEN ! mji - For MP options 2, 4, 6, 7, 8, etc. (qc = liquid, qi = ice, qs = snow) IF ( F_QI .and. F_QC .and. F_QS) THEN QCLD = QI(i,k,j)+QC(i,k,j)+QS(i,k,j) IF (QCLD .LT. QCLDMIN) THEN weight = 0. ELSE weight = (QI(i,k,j)+QS(i,k,j)) / QCLD ENDIF ENDIF ! for P3, mp option 50 or 51 IF ( F_QI .and. F_QC .and. .not. F_QS) THEN QCLD = QI(i,k,j)+QC(i,k,j) IF (QCLD .LT. QCLDMIN) THEN weight = 0. ELSE weight = (QI(i,k,j)) / QCLD ENDIF ENDIF ! mji - For MP options 1 and 3, (qc only) ! For MP=1, qc = liquid, for MP=3, qc = liquid or ice depending on temperature IF ( F_QC .and. .not. F_QI .and. .not. F_QS ) THEN QCLD = QC(i,k,j) IF (QCLD .LT. QCLDMIN) THEN weight = 0. ELSE if (t_phy(i,k,j) .gt. 273.15) weight = 0. if (t_phy(i,k,j) .le. 273.15) weight = 1. ENDIF ENDIF ! mji - For MP option 5; (qc = liquid, qs = ice) IF ( F_QC .and. .not. F_QI .and. F_QS .and. PRESENT(F_ICE_PHY) ) THEN ! Mixing ratios of cloud water & total ice (cloud ice + snow). ! Mixing ratios of rain are not considered in this scheme. ! F_ICE is fraction of ice ! F_RAIN is fraction of rain QIMID = QS(i,k,j) QWMID = QC(i,k,j) ! old method ! QIMID = QC(i,k,j)*F_ICE_PHY(i,k,j) ! QWMID = (QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j)) ! !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud, ! only cloud water + cloud ice + snow ! QCLD=QWMID+QIMID IF (QCLD .LT. QCLDMIN) THEN weight = 0. ELSE weight = F_ICE_PHY(i,k,j) ENDIF ENDIF !BSF - For HWRF MP option; (qc = liquid, qi = cloud ice+snow) ! IF ( F_QC .and. F_QI .and. .not. F_QS ) THEN #if (HWRF == 1) IF ( mp_physics .eq. FER_MP_HIRES .or. & mp_physics .eq. FER_MP_HIRES_ADVECT .or. & mp_physics .eq. ETAMP_HWRF) THEN #else IF ( mp_physics .eq. FER_MP_HIRES .or. & mp_physics==fer_mp_hires_advect) THEN #endif QIMID = QI(i,k,j) !- total ice (cloud ice + snow) QWMID = QC(i,k,j) !- cloud water QCLD=QWMID+QIMID !- cloud water + total ice IF (QCLD .LT. QCLDMIN) THEN weight = 0. ELSE weight = QIMID/QCLD if (tc<-40.) weight=1. ENDIF ENDIF ELSE CLDFRA(i,k,j)=0. ENDIF ifouter ! IF ( F_QI .and. F_QC .and. F_QS) QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity ! !--- Determine cloud fraction (modified from original algorithm) ! cldfra1_flag(i,k,j) = 0 IF (QCLD .LT. QCLDMIN) THEN ! !--- Assume zero cloud fraction if there is no cloud mixing ratio ! CLDFRA(i,k,j)=0. cldfra1_flag(i,k,j) = 1 ELSEIF(RHUM.GE.RHGRID)THEN ! !--- Assume cloud fraction of unity if near saturation and the cloud ! mixing ratio is at or above the minimum threshold ! CLDFRA(i,k,j)=1. cldfra1_flag(i,k,j) = 2 ELSE cldfra1_flag(i,k,j) = 3 ! !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995) ! modified based on assumed grid-scale saturation at RH=RHgrid. ! SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j)) DENOM=(SUBSAT)**GAMMA ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001 ! prevent negative values (new) RHUM=MAX(1.E-10, RHUM) CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG)) !! ARG=-1000*QCLD/(RHUM-RHGRID) !! ARG=MAX(ARG, ARGMIN) !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG)) IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0. ENDIF !--- End IF (QCLD .LT. QCLDMIN) ... ENDDO !--- End DO i ENDDO !--- End DO k ENDDO !--- End DO j END SUBROUTINE cal_cldfra1 !+---+-----------------------------------------------------------------+ !..Cloud fraction scheme by G. Thompson (NCAR-RAL), not intended for !.. combining with any cumulus or shallow cumulus parameterization !.. scheme cloud fractions. This is intended as a stand-alone for !.. cloud fraction and is relatively good at getting widespread stratus !.. and stratoCu without caring whether any deep/shallow Cu param schemes !.. is making sub-grid-spacing clouds/precip. Under the hood, this !.. scheme follows Mocko and Cotton (1995) in applicaiton of the !.. Sundqvist et al (1989) scheme but using a grid-scale dependent !.. RH threshold, one each for land v. ocean points based on !.. experiences with HWRF testing. !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, & & p,t,rho, XLAND, gridkm, & ! & rand_perturb_on, kme_stoch, rand_pert, & & ids,ide, jds,jde, kds,kde, & & ims,ime, jms,jme, kms,kme, & & its,ite, jts,jte, kts,kte) ! USE module_mp_thompson , ONLY : rsif, rslf IMPLICIT NONE ! INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & & ims,ime, jms,jme, kms,kme, & ! & kme_stoch, & & its,ite, jts,jte, kts,kte ! INTEGER, INTENT(IN):: rand_perturb_on REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: qv,p,t,rho REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT):: qc,qi,qs ! REAL, DIMENSION(ims:ime,kms:kme_stoch,jms:jme), INTENT(IN):: rand_pert REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN):: XLAND REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT):: cldfra REAL, INTENT(IN):: gridkm !..Local vars. REAL:: RH_00L, RH_00O, RH_00, RHI_max, entrmnt REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: qvsat INTEGER:: i,j,k REAL:: TK, TC, qvsi, qvsw, RHUM, xx, yy REAL, DIMENSION(kts:kte):: qvs1d, cfr1d, T1d, & & P1d, R1d, qc1d, qi1d, qs1d character*512 dbg_msg LOGICAL, parameter :: debug_flag=.false. !+---+ !..First cut scale-aware. Higher resolution should require closer to !.. saturated grid box for higher cloud fraction. Simple functions !.. chosen based on Mocko and Cotton (1995) starting point and desire !.. to get near 100% RH as grid spacing moves toward 1.0km, but higher !.. RH over ocean required as compared to over land. RH_00L = 0.7 + SQRT(1./(25.0+gridkm*gridkm*gridkm)) RH_00O = 0.81 + SQRT(1./(50.0+gridkm*gridkm*gridkm)) DO j = jts,jte DO k = kts,kte DO i = its,ite CLDFRA(I,K,J) = 0.0 if (qc(i,k,j).gt.1.E-6 .or. qi(i,k,j).ge.1.E-7 .or. qs(i,k,j).gt.1.E-5) then CLDFRA(I,K,J) = 1.0 qvsat(i,k,j) = qv(i,k,j) else TK = t(i,k,j) TC = TK - 273.16 qvsw = rslf(P(i,k,j), TK) qvsi = rsif(P(i,k,j), TK) if (tc .ge. -12.0) then qvsat(i,k,j) = qvsw elseif (tc .lt. -20.0) then qvsat(i,k,j) = qvsi else qvsat(i,k,j) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+20.) endif RHUM = MAX(0.01, MIN(qv(i,k,j)/qvsat(i,k,j), 0.9999)) IF ((XLAND(I,J)-1.5).GT.0.) THEN !--- Ocean RH_00 = RH_00O ELSE !--- Land RH_00 = RH_00L ENDIF if (tc .ge. -12.0) then RHUM = MIN(0.999, RHUM) CLDFRA(I,K,J) = MAX(0.0, 1.0-SQRT((1.0-RHUM)/(1.-RH_00))) elseif (tc.lt.-12..and.tc.gt.-70. .and. RHUM.gt.RH_00L) then RHUM = MAX(0.01, MIN(qv(i,k,j)/qvsat(i,k,j), 1.0 - 1.E-6)) CLDFRA(I,K,J) = MAX(0., 1.0-SQRT((1.0-RHUM)/(1.0-RH_00L))) endif CLDFRA(I,K,J) = MIN(0.90, CLDFRA(I,K,J)) endif ENDDO ENDDO ENDDO !..Prepare for a 1-D column to find various cloud layers. DO j = jts,jte DO i = its,ite ! if (i.gt.10.and.i.le.20 .and. j.gt.10.and.j.le.20) then ! debug_flag = .true. ! else ! debug_flag = .false. ! endif ! if (rand_perturb_on .eq. 1) then ! entrmnt = MAX(0.01, MIN(0.99, 0.5 + rand_pert(i,1,j)*0.5)) ! else entrmnt = 0.5 ! endif DO k = kts,kte qvs1d(k) = qvsat(i,k,j) cfr1d(k) = cldfra(i,k,j) T1d(k) = t(i,k,j) P1d(k) = p(i,k,j) R1d(k) = rho(i,k,j) qc1d(k) = qc(i,k,j) qi1d(k) = qi(i,k,j) qs1d(k) = qs(i,k,j) ENDDO ! if (debug_flag) then ! WRITE (dbg_msg,*) 'DEBUG-GT: finding cloud layers at point (', i, ', ', j, ')' ! CALL wrf_debug (150, dbg_msg) ! endif call find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & & debug_flag, qc1d, qi1d, qs1d, kts,kte) DO k = kts,kte cldfra(i,k,j) = cfr1d(k) qc(i,k,j) = qc1d(k) qi(i,k,j) = qi1d(k) ENDDO ENDDO ENDDO END SUBROUTINE cal_cldfra3 !+---+-----------------------------------------------------------------+ !..From cloud fraction array, find clouds of multi-level depth and compute !.. a reasonable value of LWP or IWP that might be contained in that depth, !.. unless existing LWC/IWC is already there. SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & & debugfl, qc1d, qi1d, qs1d, kts,kte) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: kts, kte LOGICAL, INTENT(IN):: debugfl REAL, INTENT(IN):: entrmnt REAL, DIMENSION(kts:kte), INTENT(IN):: qvs1d,T1d,P1d,R1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: cfr1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc1d, qi1d, qs1d !..Local vars. REAL, DIMENSION(kts:kte):: theta, dz REAL:: Z1, Z2, theta1, theta2, ht1, ht2 INTEGER:: k, k2, k_tropo, k_m12C, k_m40C, k_cldb, k_cldt, kbot LOGICAL:: in_cloud character*512 dbg_msg !+---+ k_m12C = 0 k_m40C = 0 DO k = kte, kts, -1 theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.)) if (T1d(k)-273.16 .gt. -40.0 .and. P1d(k).gt.7000.0) k_m40C = MAX(k_m40C, k) if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10000.0) k_m12C = MAX(k_m12C, k) ENDDO if (k_m40C .le. kts) k_m40C = kts if (k_m12C .le. kts) k_m12C = kts Z2 = 44307.692 * (1.0 - (P1d(kte)/101325.)**0.190) DO k = kte-1, kts, -1 Z1 = 44307.692 * (1.0 - (P1d(k)/101325.)**0.190) dz(k+1) = Z2 - Z1 Z2 = Z1 ENDDO dz(kts) = dz(kts+1) !..Find tropopause height, best surrogate, because we would not really !.. wish to put fake clouds into the stratosphere. The 10/1500 ratio !.. d(Theta)/d(Z) approximates a vertical line on typical SkewT chart !.. near typical (mid-latitude) tropopause height. Since messy data !.. could give us a false signal of such a transition, do the check over !.. three K-level change, not just a level-to-level check. This method !.. has potential failure in arctic-like conditions with extremely low !.. tropopause height, as would any other diagnostic, so ensure resulting !.. k_tropo level is above 4km. DO k = kte-3, kts, -1 theta1 = theta(k) theta2 = theta(k+2) ht1 = 44307.692 * (1.0 - (P1d(k)/101325.)**0.190) ht2 = 44307.692 * (1.0 - (P1d(k+2)/101325.)**0.190) if ( (((theta2-theta1)/(ht2-ht1)) .lt. 10./1500. ) .AND. & & (ht1.lt.19000.) .and. (ht1.gt.4000.) ) then goto 86 endif ENDDO 86 continue k_tropo = MAX(kts+2, k+2) ! if (debugfl) then ! print*, ' FOUND TROPOPAUSE ', k_tropo, ' near ', ht2, ' m' ! WRITE (dbg_msg,*) 'DEBUG-GT: FOUND TROPOPAUSE ', k_tropo, ' near ', ht2, ' m' ! CALL wrf_debug (150, dbg_msg) ! endif !..Eliminate possible fractional clouds above supposed tropopause. DO k = k_tropo+1, kte if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.0.999) then cfr1d(k) = 0. endif ENDDO !..We would like to prevent fractional clouds below LCL in idealized !.. situation with deep well-mixed convective PBL, that otherwise is !.. likely to get clouds in more realistic capping inversion layer. kbot = kts+2 DO k = kbot, k_m12C if ( (theta(k)-theta(k-1)) .gt. 0.05E-3*dz(k)) EXIT ENDDO kbot = MAX(kts+1, k-2) DO k = kts, kbot if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.0.999) cfr1d(k) = 0. ENDDO !..Starting below tropo height, if cloud fraction greater than 1 percent, !.. compute an approximate total layer depth of cloud, determine a total !.. liquid water/ice path (LWP/IWP), then reduce that amount with tuning !.. parameter to represent entrainment factor, then divide up LWP/IWP !.. into delta-Z weighted amounts for individual levels per cloud layer. k_cldb = k_tropo in_cloud = .false. k = k_tropo DO WHILE (.not. in_cloud .AND. k.gt.k_m12C) k_cldt = 0 if (cfr1d(k).ge.0.01) then in_cloud = .true. k_cldt = MAX(k_cldt, k) endif if (in_cloud) then DO k2 = k_cldt-1, k_m12C, -1 if (cfr1d(k2).lt.0.01 .or. k2.eq.k_m12C) then k_cldb = k2+1 goto 87 endif ENDDO 87 continue in_cloud = .false. endif if ((k_cldt - k_cldb + 1) .ge. 2) then ! if (debugfl) then ! print*, 'An ice cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 ! WRITE (dbg_msg,*) 'DEBUG-GT: An ice cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 ! CALL wrf_debug (150, dbg_msg) ! endif call adjust_cloudIce(cfr1d, qi1d, qs1d, qvs1d, T1d,R1d,dz, & & entrmnt, k_cldb,k_cldt,kts,kte) k = k_cldb else if (cfr1d(k_cldb).gt.0.and.qi1d(k_cldb).lt.1.E-6) & & qi1d(k_cldb)=1.E-5*cfr1d(k_cldb) endif k = k - 1 ENDDO k_cldb = k_tropo in_cloud = .false. k = k_m12C + 2 DO WHILE (.not. in_cloud .AND. k.gt.kbot) k_cldt = 0 if (cfr1d(k).ge.0.01) then in_cloud = .true. k_cldt = MAX(k_cldt, k) endif if (in_cloud) then DO k2 = k_cldt-1, kbot, -1 if (cfr1d(k2).lt.0.01 .or. k2.eq.kbot) then k_cldb = k2+1 goto 88 endif ENDDO 88 continue in_cloud = .false. endif if ((k_cldt - k_cldb + 1) .ge. 2) then ! if (debugfl) then ! print*, 'A water cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 ! WRITE (dbg_msg,*) 'DEBUG-GT: A water cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 ! CALL wrf_debug (150, dbg_msg) ! endif call adjust_cloudH2O(cfr1d, qc1d, qvs1d, T1d,R1d,dz, & & entrmnt, k_cldb,k_cldt,kts,kte) k = k_cldb else if (cfr1d(k_cldb).gt.0.and.qc1d(k_cldb).lt.1.E-6) & & qc1d(k_cldb)=1.E-5*cfr1d(k_cldb) endif k = k - 1 ENDDO !..Do a final total column adjustment since we may have added more than 1mm !.. LWP/IWP for multiple cloud decks. call adjust_cloudFinal(cfr1d, qc1d, qi1d, R1d,dz, kts,kte,k_tropo) ! if (debugfl) then ! print*, ' Made-up fake profile of clouds' ! do k = kte, kts, -1 ! write(*,'(i3, 2x, f8.2, 2x, f9.2, 2x, f6.2, 2x, f15.7, 2x, f15.7)') & ! & K, T1d(k)-273.15, P1d(k)*0.01, cfr1d(k)*100., qc1d(k)*1000.,qi1d(k)*1000. ! enddo ! WRITE (dbg_msg,*) 'DEBUG-GT: Made-up fake profile of clouds' ! CALL wrf_debug (150, dbg_msg) ! do k = kte, kts, -1 ! write(dbg_msg,'(f8.2, 2x, f9.2, 2x, f6.2, 2x, f15.7, 2x, f15.7)') & ! & T1d(k)-273.15, P1d(k)*0.01, cfr1d(k)*100., qc1d(k)*1000.,qi1d(k)*1000. ! CALL wrf_debug (150, dbg_msg) ! enddo ! endif END SUBROUTINE find_cloudLayers !+---+-----------------------------------------------------------------+ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs, T,Rho,dz, entr, k1,k2,kts,kte) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: k1,k2, kts,kte REAL, INTENT(IN):: entr REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qvs, T, Rho, dz REAL, DIMENSION(kts:kte), INTENT(INOUT):: qi, qs REAL:: iwc, max_iwc, tdz, this_iwc, this_dz, iwp_exists INTEGER:: k, kmid tdz = 0. do k = k1, k2 tdz = tdz + dz(k) enddo kmid = NINT(0.5*(k1+k2)) max_iwc = ABS(qvs(k2-1)-qvs(k1)) ! print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz iwp_exists = 0. do k = k1, k2 iwp_exists = iwp_exists + (qi(k)+qs(k))*Rho(k)*dz(k) enddo this_dz = 0.0 do k = k1, k2 if (k.eq.k1) then this_dz = this_dz + 0.5*dz(k) else this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz iwc = MAX(1.E-6, this_iwc*(1.-entr)) if (cfr(k).gt.0.01.and.cfr(k).lt.0.99.and.T(k).ge.203.16) then qi(k) = qi(k) + 0.1*cfr(k)*iwc elseif (qi(k).lt.1.E-5.and.cfr(k).ge.0.99.and.T(k).ge.203.16) then qi(k) = qi(k) + 0.01*iwc endif enddo END SUBROUTINE adjust_cloudIce !+---+-----------------------------------------------------------------+ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs, T,Rho,dz, entr, k1,k2,kts,kte) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: k1,k2, kts,kte REAL, INTENT(IN):: entr REAL, DIMENSION(kts:kte):: cfr, qc, qvs, T, Rho, dz REAL:: lwc, max_lwc, tdz, this_lwc, this_dz, lwp_exists INTEGER:: k, kmid tdz = 0. do k = k1, k2 tdz = tdz + dz(k) enddo kmid = NINT(0.5*(k1+k2)) max_lwc = ABS(qvs(k2-1)-qvs(k1)) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz lwp_exists = 0. do k = k1, k2 lwp_exists = lwp_exists + qc(k)*Rho(k)*dz(k) enddo this_dz = 0.0 do k = k1, k2 if (k.eq.k1) then this_dz = this_dz + 0.5*dz(k) else this_dz = this_dz + dz(k) endif this_lwc = max_lwc*this_dz/tdz lwc = MAX(1.E-6, this_lwc*(1.-entr)) if (cfr(k).gt.0.01.and.cfr(k).lt.0.99.and.T(k).lt.298.16.and.T(k).ge.253.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc elseif (cfr(k).ge.0.99.and.qc(k).lt.1.E-5.and.T(k).lt.298.16.and.T(k).ge.253.16) then qc(k) = qc(k) + 0.1*lwc endif enddo END SUBROUTINE adjust_cloudH2O !+---+-----------------------------------------------------------------+ !..Do not alter any grid-explicitly resolved hydrometeors, rather only !.. the supposed amounts due to the cloud fraction scheme. SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte,k_tropo) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: kts,kte,k_tropo REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, Rho, dz REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc, qi REAL:: lwp, iwp, xfac INTEGER:: k lwp = 0. iwp = 0. do k = kts, k_tropo if (cfr(k).gt.0.0) then lwp = lwp + qc(k)*Rho(k)*dz(k) iwp = iwp + qi(k)*Rho(k)*dz(k) endif enddo if (lwp .gt. 1.5) then xfac = 1./lwp do k = kts, k_tropo if (cfr(k).gt.0.01 .and. cfr(k).lt.0.99) then qc(k) = qc(k)*xfac endif enddo endif if (iwp .gt. 1.5) then xfac = 1./iwp do k = kts, k_tropo if (cfr(k).gt.0.01 .and. cfr(k).lt.0.99) then qi(k) = qi(k)*xfac endif enddo endif END SUBROUTINE adjust_cloudFinal !+---+-----------------------------------------------------------------+ SUBROUTINE toposhad_init(ht_shad,ht_loc,shadowmask,nested,iter, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & its,ite, jts,jte, kts,kte ) USE module_model_constants implicit none INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & its,ite, jts,jte, kts,kte LOGICAL, INTENT(IN) :: nested REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad, ht_loc INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask INTEGER, INTENT(IN) :: iter ! Local variables INTEGER :: i, j if (iter.eq.1) then ! Initialize shadow mask do j=jts,jte do i=its,ite shadowmask(i,j) = 0 ENDDO ENDDO ! Initialize shading height IF ( nested ) THEN ! Do not overwrite input from parent domain do j=max(jts,jds+2),min(jte,jde-3) do i=max(its,ids+2),min(ite,ide-3) ht_shad(i,j) = ht_loc(i,j)-0.001 ENDDO ENDDO ELSE do j=jts,jte do i=its,ite ht_shad(i,j) = ht_loc(i,j)-0.001 ENDDO ENDDO ENDIF IF ( nested ) THEN ! Check if a shadow exceeding the topography height is available at the lateral domain edge from nesting if (its.eq.ids) then do j=jts,jte if (ht_shad(its,j) .gt. ht_loc(its,j)) then shadowmask(its,j) = 1 ht_loc(its,j) = ht_shad(its,j) endif if (ht_shad(its+1,j) .gt. ht_loc(its+1,j)) then shadowmask(its+1,j) = 1 ht_loc(its+1,j) = ht_shad(its+1,j) endif enddo endif if (ite.eq.ide-1) then do j=jts,jte if (ht_shad(ite,j) .gt. ht_loc(ite,j)) then shadowmask(ite,j) = 1 ht_loc(ite,j) = ht_shad(ite,j) endif if (ht_shad(ite-1,j) .gt. ht_loc(ite-1,j)) then shadowmask(ite-1,j) = 1 ht_loc(ite-1,j) = ht_shad(ite-1,j) endif enddo endif if (jts.eq.jds) then do i=its,ite if (ht_shad(i,jts) .gt. ht_loc(i,jts)) then shadowmask(i,jts) = 1 ht_loc(i,jts) = ht_shad(i,jts) endif if (ht_shad(i,jts+1) .gt. ht_loc(i,jts+1)) then shadowmask(i,jts+1) = 1 ht_loc(i,jts+1) = ht_shad(i,jts+1) endif enddo endif if (jte.eq.jde-1) then do i=its,ite if (ht_shad(i,jte) .gt. ht_loc(i,jte)) then shadowmask(i,jte) = 1 ht_loc(i,jte) = ht_shad(i,jte) endif if (ht_shad(i,jte-1) .gt. ht_loc(i,jte-1)) then shadowmask(i,jte-1) = 1 ht_loc(i,jte-1) = ht_shad(i,jte-1) endif enddo endif ENDIF else ! Fill the local topography field at the points next to internal tile boundaries with ht_shad values ! A 2-pt halo has been applied to the ht_shad before the repeated call of this subroutine if ((its.ne.ids).and.(its.eq.ips)) then do j=jts-2,jte+2 ht_loc(its-1,j) = max(ht_loc(its-1,j),ht_shad(its-1,j)) ht_loc(its-2,j) = max(ht_loc(its-2,j),ht_shad(its-2,j)) enddo endif if ((ite.ne.ide-1).and.(ite.eq.ipe)) then do j=jts-2,jte+2 ht_loc(ite+1,j) = max(ht_loc(ite+1,j),ht_shad(ite+1,j)) ht_loc(ite+2,j) = max(ht_loc(ite+2,j),ht_shad(ite+2,j)) enddo endif if ((jts.ne.jds).and.(jts.eq.jps)) then do i=its-2,ite+2 ht_loc(i,jts-1) = max(ht_loc(i,jts-1),ht_shad(i,jts-1)) ht_loc(i,jts-2) = max(ht_loc(i,jts-2),ht_shad(i,jts-2)) enddo endif if ((jte.ne.jde-1).and.(jte.eq.jpe)) then do i=its-2,ite+2 ht_loc(i,jte+1) = max(ht_loc(i,jte+1),ht_shad(i,jte+1)) ht_loc(i,jte+2) = max(ht_loc(i,jte+2),ht_shad(i,jte+2)) enddo endif endif END SUBROUTINE toposhad_init SUBROUTINE toposhad(xlat,xlong,sina,cosa,xtime,gmt,radfrq,declin, & dx,dy,ht_shad,ht_loc,iter, & shadowmask,shadlen, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & its,ite, jts,jte, kts,kte ) USE module_model_constants implicit none INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & ips,ipe, jps,jpe, kps,kpe, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN) :: iter REAL, INTENT(IN) :: RADFRQ,XTIME,DECLIN,dx,dy,gmt,shadlen REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT, XLONG, sina, cosa REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask ! Local variables REAL :: pi, xt24, wgt, ri, rj, argu, sol_azi, topoelev, dxabs, tloctm, hrang, xxlat, csza INTEGER :: gpshad, ii, jj, i1, i2, j1, j2, i, j XT24=MOD(XTIME+RADFRQ*0.5,1440.) pi = 4.*atan(1.) gpshad = int(shadlen/dx+1.) if (iter.eq.1) then j_loop1: DO J=jts,jte i_loop1: DO I=its,ite TLOCTM=GMT+XT24/60.+XLONG(i,j)/15. HRANG=15.*(TLOCTM-12.)*DEGRAD XXLAT=XLAT(i,j)*DEGRAD CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG) if (csza.lt.1.e-2) then ! shadow mask does not need to be computed shadowmask(i,j) = 0 ht_shad(i,j) = ht_loc(i,j)-0.001 goto 120 endif ! Solar azimuth angle argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT)) if (argu.gt.1) argu = 1 if (argu.lt.-1) argu = -1 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun if (cosa(i,j).ge.0) then sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid else sol_azi = sol_azi + pi - asin(sina(i,j)) endif ! Scan for higher surrounding topography if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter do jj = j+1,j+gpshad ri = i + (jj-j)*tan(sol_azi) i1 = int(ri) i2 = i1+1 wgt = ri-i1 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2) if ((jj.ge.jpe+3).or.(i1.le.ips-3).or.(i2.ge.ipe+3)) then ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1 goto 120 endif topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter do ii = i+1,i+gpshad rj = j - (ii-i)*tan(pi/2.+sol_azi) j1 = int(rj) j2 = j1+1 wgt = rj-j1 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2) if ((ii.ge.ipe+3).or.(j1.le.jps-3).or.(j2.ge.jpe+3)) then ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1 goto 120 endif topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter do jj = j-1,j-gpshad,-1 ri = i + (jj-j)*tan(sol_azi) i1 = int(ri) i2 = i1+1 wgt = ri-i1 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2) if ((jj.le.jps-3).or.(i1.le.ips-3).or.(i2.ge.ipe+3)) then ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1 goto 120 endif topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else ! sun is in the western quarter do ii = i-1,i-gpshad,-1 rj = j - (ii-i)*tan(pi/2.+sol_azi) j1 = int(rj) j2 = j1+1 wgt = rj-j1 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2) if ((ii.le.ips-3).or.(j1.le.jps-3).or.(j2.ge.jpe+3)) then ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1 goto 120 endif topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo endif 120 continue ENDDO i_loop1 ENDDO j_loop1 else ! iteration > 1 j_loop2: DO J=jts,jte i_loop2: DO I=its,ite ! if (shadowmask(i,j).eq.-1) then ! this indicates that the search ended at a lateral boundary during iteration 1 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15. HRANG=15.*(TLOCTM-12.)*DEGRAD XXLAT=XLAT(i,j)*DEGRAD CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG) if (csza.lt.1.e-2) then ! shadow mask does not need to be computed shadowmask(i,j) = 0 ht_shad(i,j) = ht_loc(i,j)-0.001 goto 220 endif ! Solar azimuth angle argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT)) if (argu.gt.1) argu = 1 if (argu.lt.-1) argu = -1 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun if (cosa(i,j).ge.0) then sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid else sol_azi = sol_azi + pi - asin(sina(i,j)) endif ! Scan for higher surrounding topography if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter do jj = j+1,j+gpshad ri = i + (jj-j)*tan(sol_azi) i1 = int(ri) i2 = i1+1 wgt = ri-i1 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2) if ((jj.ge.min(jde,jpe+3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter do ii = i+1,i+gpshad rj = j - (ii-i)*tan(pi/2.+sol_azi) j1 = int(rj) j2 = j1+1 wgt = rj-j1 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2) if ((ii.ge.min(ide,ipe+3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter do jj = j-1,j-gpshad,-1 ri = i + (jj-j)*tan(sol_azi) i1 = int(ri) i2 = i1+1 wgt = ri-i1 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2) if ((jj.le.max(jds-1,jps-3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo else ! sun is in the western quarter do ii = i-1,i-gpshad,-1 rj = j - (ii-i)*tan(pi/2.+sol_azi) j1 = int(rj) j2 = j1+1 wgt = rj-j1 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2) if ((ii.le.max(ids-1,ips-3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs) if (sin(topoelev).ge.csza) then shadowmask(i,j) = 1 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza)))) endif enddo endif 220 continue ! endif ENDDO i_loop2 ENDDO j_loop2 endif ! iteration END SUBROUTINE toposhad SUBROUTINE ozn_time_int(julday,julian,ozmixm,ozmixt,levsiz,num_months, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! adapted from oznint from CAM module ! input: ozmixm - read from physics_init ! output: ozmixt - time interpolated ! USE module_ra_cam_support, ONLY : getfactors IMPLICIT NONE INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: levsiz, num_months REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), & INTENT(IN ) :: ozmixm INTEGER, INTENT(IN ) :: JULDAY REAL, INTENT(IN ) :: JULIAN REAL, DIMENSION( ims:ime, levsiz, jms:jme ), & INTENT(OUT ) :: ozmixt !Local REAL :: intJULIAN integer :: np1,np,nm,m,k,i,j integer :: IJUL integer, dimension(12) :: date_oz data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/ real, parameter :: daysperyear = 365. ! number of days in a year real :: cdayozp, cdayozm real :: fact1, fact2, deltat logical :: finddate logical :: ozncyc CHARACTER(LEN=256) :: msgstr ozncyc = .true. ! JULIAN starts from 0.0 at 0Z on 1 Jan. intJULIAN = JULIAN + 1.0 ! offset by one day ! jan 1st 00z is julian=1.0 here IJUL=INT(intJULIAN) ! Note that following will drift. ! Need to use actual month/day info to compute julian. intJULIAN=intJULIAN-FLOAT(IJUL) IJUL=MOD(IJUL,365) IF(IJUL.EQ.0)IJUL=365 intJULIAN=intJULIAN+IJUL np1=1 finddate=.false. ! do m=1,num_months do m=1,12 if(date_oz(m).gt.intjulian.and..not.finddate) then np1=m finddate=.true. endif enddo cdayozp=date_oz(np1) if(np1.gt.1) then cdayozm=date_oz(np1-1) np=np1 nm=np-1 else cdayozm=date_oz(12) np=np1 nm=12 endif ! call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, & ! fact1, fact2) ! ! Determine time interpolation factors. Account for December-January ! interpolation if dataset is being cycled yearly. ! if (ozncyc .and. np1 == 1) then ! Dec-Jan interpolation deltat = cdayozp + daysperyear - cdayozm if (intjulian > cdayozp) then ! We are in December fact1 = (cdayozp + daysperyear - intjulian)/deltat fact2 = (intjulian - cdayozm)/deltat else ! We are in January fact1 = (cdayozp - intjulian)/deltat fact2 = (intjulian + daysperyear - cdayozm)/deltat end if else deltat = cdayozp - cdayozm fact1 = (cdayozp - intjulian)/deltat fact2 = (intjulian - cdayozm)/deltat end if ! ! Time interpolation. ! do j=jts,jte do k=1,levsiz do i=its,ite ozmixt(i,k,j) = ozmixm(i,k,j,nm)*fact1 + ozmixm(i,k,j,np)*fact2 end do end do end do END SUBROUTINE ozn_time_int SUBROUTINE ozn_p_int(p ,pin, levsiz, ozmixt, o3vmr, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) !----------------------------------------------------------------------- ! ! Purpose: Interpolate ozone from current time-interpolated values to model levels ! ! Method: Use pressure values to determine interpolation levels ! ! Author: Bruce Briegleb ! WW: Adapted for general use ! !-------------------------------------------------------------------------- implicit none !-------------------------------------------------------------------------- ! ! Arguments ! INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, intent(in) :: levsiz ! number of ozone layers real, intent(in) :: p(ims:ime,kms:kme,jms:jme) ! level pressures (mks, bottom-up) real, intent(in) :: pin(levsiz) ! ozone data level pressures (mks, top-down) real, intent(in) :: ozmixt(ims:ime,levsiz,jms:jme) ! ozone mixing ratio real, intent(out) :: o3vmr(ims:ime,kms:kme,jms:jme) ! ozone volume mixing ratio ! ! local storage ! real pmid(its:ite,kts:kte) integer i,j ! longitude index integer k, kk, kkstart, kout! level indices integer kupper(its:ite) ! Level indices for interpolation integer kount ! Counter integer ncol, pver real dpu ! upper level pressure difference real dpl ! lower level pressure difference ncol = ite - its + 1 pver = kte - kts + 1 do j=jts,jte ! ! Initialize index array ! ! do i=1, ncol do i=its, ite kupper(i) = 1 end do ! ! Reverse the pressure array, and pin is in Pa, the same as model pmid ! do k = kts,kte kk = kte - k + kts do i = its,ite pmid(i,kk) = p(i,k,j) enddo enddo do k=1,pver kout = pver - k + 1 ! kout = k ! ! Top level we need to start looking is the top level for the previous k ! for all longitude points ! kkstart = levsiz ! do i=1,ncol do i=its,ite kkstart = min0(kkstart,kupper(i)) end do kount = 0 ! ! Store level indices for interpolation ! do kk=kkstart,levsiz-1 ! do i=1,ncol do i=its,ite if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then kupper(i) = kk kount = kount + 1 end if end do ! ! If all indices for this level have been found, do the interpolation and ! go to the next level ! if (kount.eq.ncol) then ! do i=1,ncol do i=its,ite dpu = pmid(i,k) - pin(kupper(i)) dpl = pin(kupper(i)+1) - pmid(i,k) o3vmr(i,kout,j) = (ozmixt(i,kupper(i),j)*dpl + & ozmixt(i,kupper(i)+1,j)*dpu)/(dpl + dpu) end do goto 35 end if end do ! ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and ! must extrapolate from the bottom or top ozone data level for at least some ! of the longitude points. ! ! do i=1,ncol do i=its,ite if (pmid(i,k) .lt. pin(1)) then o3vmr(i,kout,j) = ozmixt(i,1,j)*pmid(i,k)/pin(1) else if (pmid(i,k) .gt. pin(levsiz)) then o3vmr(i,kout,j) = ozmixt(i,levsiz,j) else dpu = pmid(i,k) - pin(kupper(i)) dpl = pin(kupper(i)+1) - pmid(i,k) o3vmr(i,kout,j) = (ozmixt(i,kupper(i),j)*dpl + & ozmixt(i,kupper(i)+1,j)*dpu)/(dpl + dpu) end if end do if (kount.gt.ncol) then ! call endrun ('OZN_P_INT: Bad ozone data: non-monotonicity suspected') call wrf_error_fatal ('OZN_P_INT: Bad ozone data: non-monotonicity suspected') end if 35 continue end do end do return END SUBROUTINE ozn_p_int SUBROUTINE aer_time_int(julday,julian,aerodm,aerodt,levsiz,num_months,no_src, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) ! adapted from oznint from CAM module ! input: aerodm - read from physics_init ! output: aerodt - time interpolated ! USE module_ra_cam_support, ONLY : getfactors IMPLICIT NONE INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: levsiz, num_months, no_src REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months, no_src ), & INTENT(IN ) :: aerodm INTEGER, INTENT(IN ) :: JULDAY REAL, INTENT(IN ) :: JULIAN REAL, DIMENSION( ims:ime, levsiz, jms:jme, no_src ), & INTENT(OUT ) :: aerodt !Local REAL :: intJULIAN integer :: np1,np,nm,m,k,i,j,s integer :: IJUL integer, dimension(12) :: date_oz data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/ real, parameter :: daysperyear = 365. ! number of days in a year real :: cdayozp, cdayozm real :: fact1, fact2, deltat logical :: finddate logical :: ozncyc CHARACTER(LEN=256) :: msgstr ozncyc = .true. ! JULIAN starts from 0.0 at 0Z on 1 Jan. intJULIAN = JULIAN + 1.0 ! offset by one day ! jan 1st 00z is julian=1.0 here IJUL=INT(intJULIAN) ! Note that following will drift. ! Need to use actual month/day info to compute julian. intJULIAN=intJULIAN-FLOAT(IJUL) IJUL=MOD(IJUL,365) IF(IJUL.EQ.0)IJUL=365 intJULIAN=intJULIAN+IJUL np1=1 finddate=.false. ! do m=1,num_months do m=1,12 if(date_oz(m).gt.intjulian.and..not.finddate) then np1=m finddate=.true. endif enddo cdayozp=date_oz(np1) if(np1.gt.1) then cdayozm=date_oz(np1-1) np=np1 nm=np-1 else cdayozm=date_oz(12) np=np1 nm=12 endif ! call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, & ! fact1, fact2) ! ! Determine time interpolation factors. Account for December-January ! interpolation if dataset is being cycled yearly. ! if (ozncyc .and. np1 == 1) then ! Dec-Jan interpolation deltat = cdayozp + daysperyear - cdayozm if (intjulian > cdayozp) then ! We are in December fact1 = (cdayozp + daysperyear - intjulian)/deltat fact2 = (intjulian - cdayozm)/deltat else ! We are in January fact1 = (cdayozp - intjulian)/deltat fact2 = (intjulian + daysperyear - cdayozm)/deltat end if else deltat = cdayozp - cdayozm fact1 = (cdayozp - intjulian)/deltat fact2 = (intjulian - cdayozm)/deltat end if ! ! Time interpolation. ! do s=1, no_src do j=jts,jte do k=1,levsiz do i=its,ite aerodt(i,k,j,s) = aerodm(i,k,j,nm,s)*fact1 + aerodm(i,k,j,np,s)*fact2 end do end do end do end do END SUBROUTINE aer_time_int SUBROUTINE aer_p_int(p ,pin, levsiz, aerodt, aerod, no_src, pf, totaod, & ids , ide , jds , jde , kds , kde , & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) !----------------------------------------------------------------------- ! ! Purpose: Interpolate aerosol from current time-interpolated values to model levels ! ! Method: Use pressure values to determine interpolation levels ! ! Author: Bruce Briegleb ! WW: Adapted for general use ! ! p: model level pressure at half levels (Pa, bottom-up) ! pf: model level pressure at full levles (Pa, bottom-up) ! !-------------------------------------------------------------------------- implicit none !-------------------------------------------------------------------------- ! ! Arguments ! INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte integer, intent(in) :: levsiz ! number of aerosol layers integer, intent(in) :: no_src ! types of aerosol real, intent(in) :: p(ims:ime,kms:kme,jms:jme) real, intent(in) :: pf(ims:ime,kms:kme,jms:jme) real, intent(in) :: pin(levsiz) ! aerosol data level pressures (mks, top-down) real, intent(in) :: aerodt(ims:ime,levsiz,jms:jme,1:no_src) ! aerosol optical depth real, intent(out) :: aerod(ims:ime,kms:kme,jms:jme,1:no_src) ! aerosol optical depth real, intent(out) :: totaod(ims:ime,jms:jme) ! total aerosol optical depth ! ! local storage ! real pmid(its:ite,kts:kte) integer i,j ! longitude index integer k, kk, kkstart, kout! level indices integer kupper(its:ite) ! Level indices for interpolation integer kount ! Counter integer ncol, pver, s real dpu ! upper level pressure difference real dpl ! lower level pressure difference real dpm ! pressure difference in a model layer surrounding half p ncol = ite - its + 1 pver = kte - kts + 1 do s=1,no_src do j=jts,jte ! ! Initialize index array ! do i=its, ite kupper(i) = 1 end do ! ! The pressure from incoming data is in hPa and top-down, ! while model pressure is in Pa and bottom-up ! do k = kts,kte kk = kte - k + kts do i = its,ite pmid(i,kk) = p(i,k,j)*0.01 enddo enddo do k=1,pver kout = pver - k + 1 ! ! Top level we need to start looking is the top level for the previous k ! for all longitude points ! kkstart = levsiz do i=its,ite kkstart = min0(kkstart,kupper(i)) end do kount = 0 ! ! Store level indices for interpolation ! do kk=kkstart,levsiz-1 do i=its,ite if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then kupper(i) = kk kount = kount + 1 end if end do ! ! If all indices for this level have been found, do the interpolation and ! go to the next level ! if (kount.eq.ncol) then do i=its,ite dpu = pmid(i,k) - pin(kupper(i)) dpl = pin(kupper(i)+1) - pmid(i,k) dpm = pf(i,kout,j) - pf(i,kout+1,j) aerod(i,kout,j,s) = (aerodt(i,kupper(i),j,s)*dpl + & aerodt(i,kupper(i)+1,j,s)*dpu)/(dpl + dpu) aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm end do goto 35 end if end do ! ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and ! must extrapolate from the bottom or top aerosol data level for at least some ! of the longitude points. ! do i=its,ite if (pmid(i,k) .lt. pin(1)) then dpm = pf(i,kout,j) - pf(i,kout+1,j) aerod(i,kout,j,s) = aerodt(i,1,j,s)*pmid(i,k)/pin(1) aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm else if (pmid(i,k) .gt. pin(levsiz)) then dpm = pf(i,kout,j) - pf(i,kout+1,j) aerod(i,kout,j,s) = aerodt(i,levsiz,j,s) aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm else dpu = pmid(i,k) - pin(kupper(i)) dpl = pin(kupper(i)+1) - pmid(i,k) dpm = pf(i,kout,j) - pf(i,kout+1,j) aerod(i,kout,j,s) = (aerodt(i,kupper(i),j,s)*dpl + & aerodt(i,kupper(i)+1,j,s)*dpu)/(dpl + dpu) aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm end if end do if (kount.gt.ncol) then call wrf_error_fatal ('AER_P_INT: Bad aerosol data: non-monotonicity suspected') end if 35 continue end do end do end do do j=jts,jte do i=its,ite totaod(i,j) = 0. end do end do do s=1,no_src do j=jts,jte do k=1,pver do i=its,ite totaod(i,j) = totaod(i,j) + aerod(i,k,j,s) end do end do end do end do return END SUBROUTINE aer_p_int !+---+-----------------------------------------------------------------+ SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa, taod5503d, & & ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte) USE module_mp_thompson, only: RSLF IMPLICIT NONE INTEGER, INTENT(IN):: ims,ime, jms,jme, kms,kme, & & its,ite, jts,jte, kts,kte REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: & & t_phy,p_phy, DZ8W, & & qvapor, nifa, nwfa REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT):: taod5503d !..Local variables. REAL, DIMENSION(its:ite,kts:kte,jts:jte):: AOD_wfa, AOD_ifa REAL:: RH, a_RH,b_RH, rh_d,rh_f, rhoa,qvsat, unit_bext1,unit_bext3 REAL:: ntemp INTEGER :: i, k, j, RH_idx, RH_idx1, RH_idx2, t_idx INTEGER, PARAMETER:: rind=8 REAL, DIMENSION(rind), PARAMETER:: rh_arr = & & (/10., 60., 70., 80., 85., 90., 95., 99.8/) REAL, DIMENSION(rind,4,2) :: lookup_tabl ! RH, temp, water-friendly, ice-friendly lookup_tabl(1,1,1) = 5.73936E-15 lookup_tabl(1,1,2) = 2.63577E-12 lookup_tabl(1,2,1) = 5.73936E-15 lookup_tabl(1,2,2) = 2.63577E-12 lookup_tabl(1,3,1) = 5.73936E-15 lookup_tabl(1,3,2) = 2.63577E-12 lookup_tabl(1,4,1) = 5.73936E-15 lookup_tabl(1,4,2) = 2.63577E-12 lookup_tabl(2,1,1) = 6.93515E-15 lookup_tabl(2,1,2) = 2.72095E-12 lookup_tabl(2,2,1) = 6.93168E-15 lookup_tabl(2,2,2) = 2.72092E-12 lookup_tabl(2,3,1) = 6.92570E-15 lookup_tabl(2,3,2) = 2.72091E-12 lookup_tabl(2,4,1) = 6.91833E-15 lookup_tabl(2,4,2) = 2.72087E-12 lookup_tabl(3,1,1) = 7.24707E-15 lookup_tabl(3,1,2) = 2.77219E-12 lookup_tabl(3,2,1) = 7.23809E-15 lookup_tabl(3,2,2) = 2.77222E-12 lookup_tabl(3,3,1) = 7.23108E-15 lookup_tabl(3,3,2) = 2.77201E-12 lookup_tabl(3,4,1) = 7.21800E-15 lookup_tabl(3,4,2) = 2.77111E-12 lookup_tabl(4,1,1) = 8.95130E-15 lookup_tabl(4,1,2) = 2.87263E-12 lookup_tabl(4,2,1) = 9.01582E-15 lookup_tabl(4,2,2) = 2.87252E-12 lookup_tabl(4,3,1) = 9.13216E-15 lookup_tabl(4,3,2) = 2.87241E-12 lookup_tabl(4,4,1) = 9.16219E-15 lookup_tabl(4,4,2) = 2.87211E-12 lookup_tabl(5,1,1) = 1.06695E-14 lookup_tabl(5,1,2) = 2.96752E-12 lookup_tabl(5,2,1) = 1.06370E-14 lookup_tabl(5,2,2) = 2.96726E-12 lookup_tabl(5,3,1) = 1.05999E-14 lookup_tabl(5,3,2) = 2.96702E-12 lookup_tabl(5,4,1) = 1.05443E-14 lookup_tabl(5,4,2) = 2.96603E-12 lookup_tabl(6,1,1) = 1.37908E-14 lookup_tabl(6,1,2) = 3.15081E-12 lookup_tabl(6,2,1) = 1.37172E-14 lookup_tabl(6,2,2) = 3.15020E-12 lookup_tabl(6,3,1) = 1.36362E-14 lookup_tabl(6,3,2) = 3.14927E-12 lookup_tabl(6,4,1) = 1.35287E-14 lookup_tabl(6,4,2) = 3.14817E-12 lookup_tabl(7,1,1) = 2.26019E-14 lookup_tabl(7,1,2) = 3.66798E-12 lookup_tabl(7,2,1) = 2.24435E-14 lookup_tabl(7,2,2) = 3.66540E-12 lookup_tabl(7,3,1) = 2.23254E-14 lookup_tabl(7,3,2) = 3.66173E-12 lookup_tabl(7,4,1) = 2.20496E-14 lookup_tabl(7,4,2) = 3.65796E-12 lookup_tabl(8,1,1) = 4.41983E-13 lookup_tabl(8,1,2) = 7.50091E-11 lookup_tabl(8,2,1) = 3.93335E-13 lookup_tabl(8,2,2) = 6.79097E-11 lookup_tabl(8,3,1) = 3.45569E-13 lookup_tabl(8,3,2) = 6.07845E-11 lookup_tabl(8,4,1) = 2.96971E-13 lookup_tabl(8,4,2) = 5.36085E-11 DO j=jts,jte DO k=kts,kte DO i=its,ite AOD_wfa(i,k,j) = 0. AOD_ifa(i,k,j) = 0. END DO END DO END DO DO j=jts,jte DO k=kts,kte DO i=its,ite rhoa = p_phy(i,k,j)/(287.*t_phy(i,k,j)) t_idx = MAX(1, MIN(nint(10.999-0.0333*t_phy(i,k,j)),4)) qvsat = rslf(p_phy(i,k,j),t_phy(i,k,j)) RH = MIN(98., MAX(10.1, qvapor(i,k,j)/qvsat*100.)) !..Get the index for the RH array element if (RH .lt. 60) then RH_idx1 = 1 RH_idx2 = 2 elseif (RH .ge. 60 .AND. RH.lt.80) then a_RH = 0.1 b_RH = -4 RH_idx = nint(a_RH*RH+b_RH) rh_d = rh-rh_arr(rh_idx) if (rh_d .lt. 0) then RH_idx1 = RH_idx-1 RH_idx2 = RH_idx else RH_idx1 = RH_idx RH_idx2 = RH_idx+1 if (RH_idx2.gt.rind) then RH_idx2 = rind RH_idx1 = rind-1 endif endif else a_RH = 0.2 b_RH = -12. RH_idx = MIN(rind, nint(a_RH*RH+b_RH)) rh_d = rh-rh_arr(rh_idx) if (rh_d .lt. 0) then RH_idx1 = RH_idx-1 RH_idx2 = RH_idx else RH_idx1 = RH_idx RH_idx2 = RH_idx+1 if (RH_idx2.gt.rind) then RH_idx2 = rind RH_idx1 = rind-1 endif endif endif !..RH fraction to be used rh_f = MAX(0., MIN(1.0, (rh/(100-rh)-rh_arr(rh_idx1) & & /(100-rh_arr(rh_idx1))) & & /(rh_arr(rh_idx2)/(100-rh_arr(rh_idx2)) & & -rh_arr(rh_idx1)/(100-rh_arr(rh_idx1))) )) unit_bext1 = lookup_tabl(RH_idx1,t_idx,1) & & + (lookup_tabl(RH_idx2,t_idx,1) & & - lookup_tabl(RH_idx1,t_idx,1))*rh_f unit_bext3 = lookup_tabl(RH_idx1,t_idx,2) & & + (lookup_tabl(RH_idx2,t_idx,2) & & - lookup_tabl(RH_idx1,t_idx,2))*rh_f ntemp = MAX(1., MIN(99999.E6, nwfa(i,k,j))) AOD_wfa(i,k,j) = unit_bext1*ntemp*dz8w(i,k,j)*rhoa ntemp = MAX(0.01, MIN(9999.E6, nifa(i,k,j))) AOD_ifa(i,k,j) = unit_bext3*ntemp*dz8w(i,k,j)*rhoa END DO END DO END DO DO j=jts,jte DO k=kts,kte DO i=its,ite taod5503d(i,k,j) = aod_wfa(i,k,j) + aod_ifa(i,k,j) END DO END DO END DO END SUBROUTINE gt_aod !+---+-----------------------------------------------------------------+ END MODULE module_radiation_driver