| blob | 924c820086aab14ed3a1a22ea57e3a61f53bc5ac |
2 !WRF:MEDIATION_LAYER:PHYSICS
3 !
4 MODULE module_radiation_driver
5 real, private, parameter :: ALBEDO_MIN = 0.0, ALBEDO_MAX = 1.0, ANGEXP_MIN = 0.0, AOD_MIN = 0.0, &
6 AERASY_MIN = -1.0, AERASY_MAX = 1.0, QVAPOR_MIN = 0.0, QCLOUD_MIN = 0.0, QSNOW_MIN = 0.0
7 CONTAINS
8 !BOP
9 ! !IROUTINE: radiation_driver - interface to radiation physics options
11 ! !INTERFACE:
12 SUBROUTINE radiation_driver ( &
13 ACFRCV ,ACFRST ,ALBEDO &
14 ,CFRACH ,CFRACL ,CFRACM &
15 ,CUPPT ,CZMEAN ,DT &
16 ,DZ8W ,EMISS ,GLW &
17 ,GMT ,GSW ,HBOT &
18 ,HTOP ,HBOTR ,HTOPR &
19 ,ICLOUD &
20 ,cldovrlp ,idcor &
21 ,ITIMESTEP,JULDAY, JULIAN &
22 ,JULYR ,LW_PHYSICS &
23 ,NCFRCV ,NCFRST ,NPHS &
24 ,O3INPUT, O3RAD &
25 ,AER_OPT, aerod &
26 ,swint_opt &
27 ,solar_opt &
28 ,P8W ,P ,PI &
29 ,p_top &
30 ,RADT ,RA_CALL_OFFSET &
31 ,RHO ,RLWTOA &
32 ,RSWTOA ,RTHRATEN &
33 ,RTHRATENLW ,RTHRATENSW &
34 ,RTHRATENLWC ,RTHRATENSWC &
35 ,SNOW ,STEPRA ,SWDOWN &
36 ,SWDOWNC ,SW_PHYSICS &
37 ,SW_ECLIPSE &
38 ,T8W ,T ,TAUCLDC &
39 ,TAUCLDI ,TSK ,VEGFRA &
40 ,WARM_RAIN ,XICE ,XLAND &
41 ,XLAT ,XLONG ,YR &
42 !Optional solar variables
43 ,DECLINX ,SOLCONX ,COSZEN ,HRANG &
44 , CEN_LAT &
45 ,Z &
46 ,ALEVSIZ, no_src_types &
47 ,LEVSIZ, N_OZMIXM &
48 ,N_AEROSOLC &
49 ,PAERLEV ,ID &
50 ,CAM_ABS_DIM1, CAM_ABS_DIM2 &
51 ,CAM_ABS_FREQ_S &
52 ,XTIME &
53 ,CURR_SECS, ADAPT_STEP_FLAG &
54 ,SWDOWN2, SWDDNI2, SWDDIF2, SWDDIR2 &
55 ,SWDOWNC2, SWDDNIC2 &
56 ,couple_farms &
57 ! WRF-Solar EPS
58 ,multi_perturb &
59 ,pert_farms &
60 ,perts_albedo, perts_aod, perts_angstrom, perts_assymfac &
61 ,perts_qvapor, perts_qcloud, perts_qsnow &
62 ,pert_farms_albedo, pert_farms_aod &
63 ,pert_farms_angexp, pert_farms_aerasy &
64 ,pert_farms_qv, pert_farms_qc &
65 ,pert_farms_qs &
66 ,pert_cld3 &
67 ,perts_th &
68 ,pert_cld3_qv, pert_cld3_t &
69 !BSINGH - For WRFCuP scheme (optional args)
70 ,cu_physics,shallowcu_forced_ra &
71 ,cubot,cutop,cldfra_cup &
72 ,shall &
73 !BSINGH - ENDS
74 ! indexes
75 ,IDS,IDE, JDS,JDE, KDS,KDE &
76 ,IMS,IME, JMS,JME, KMS,KME &
77 ,i_start,i_end &
78 ,j_start,j_end &
79 ,kts, kte &
80 ,num_tiles &
81 ! Optional
82 , TLWDN, TLWUP &
83 , SLWDN, SLWUP &
84 , TSWDN, TSWUP &
85 , SSWDN, SSWUP &
86 , RE_CLOUD_GSFC &
87 , RE_RAIN_GSFC &
88 , RE_ICE_GSFC &
89 , RE_SNOW_GSFC &
90 , RE_GRAUPEL_GSFC &
91 , RE_HAIL_GSFC &
92 , COD2D_OUT &
93 , CTOP2D_OUT &
94 , CLDFRA,CLDFRA_MP_ALL,CLDT,ZNU &
95 , CCLDFRA, QCCONV, QICONV &
96 , bmj_rad_feedback &
97 #if (EM_CORE == 1)
98 , lradius,iradius &
99 #endif
100 , cldfra_dp, cldfra_sh &
101 , re_cloud, re_ice, re_snow &
102 , has_reqc, has_reqi, has_reqs &
103 , PB &
104 , F_ICE_PHY,F_RAIN_PHY &
105 , F_QNC &
106 , QNC_CURR &
107 , QV, F_QV &
108 , QC, F_QC &
109 , QR, F_QR &
110 , QI, F_QI &
111 , QI2, F_QI2 &
112 , QI3, F_QI3 &
113 , QS, F_QS &
114 , QG, F_QG &
115 , QH, F_QH &
116 , QNDROP, F_QNDROP &
117 ,QNIFA,F_QNIFA &
118 ,QNWFA,F_QNWFA &
119 ,QNBCA &
120 ,wif_input_opt &
121 ,qc_tot, qi_tot &
122 ,ACSWUPT ,ACSWUPTC &
123 ,ACSWDNT ,ACSWDNTC &
124 ,ACSWUPB ,ACSWUPBC &
125 ,ACSWDNB ,ACSWDNBC &
126 ,ACLWUPT ,ACLWUPTC &
127 ,ACLWDNT ,ACLWDNTC &
128 ,ACLWUPB ,ACLWUPBC &
129 ,ACLWDNB ,ACLWDNBC &
130 ,SWUPT ,SWUPTC, SWUPTCLN &
131 ,SWDNT ,SWDNTC, SWDNTCLN &
132 ,SWUPB ,SWUPBC, SWUPBCLN &
133 ,SWDNB ,SWDNBC, SWDNBCLN &
134 ,LWUPT ,LWUPTC, LWUPTCLN &
135 ,LWDNT ,LWDNTC, LWDNTCLN &
136 ,LWUPB ,LWUPBC, LWUPBCLN &
137 ,LWDNB ,LWDNBC, LWDNBCLN &
138 ,LWCF &
139 ,SWCF &
140 ,OLR &
141 ,aerodm, PINA, aodtot &
142 ,OZMIXM, PIN &
143 ,M_PS_1, M_PS_2, AEROSOLC_1 &
144 ,AEROSOLC_2, M_HYBI0 &
145 ,ABSTOT, ABSNXT, EMSTOT &
146 ,ICLOUD_CU &
147 ,CALC_CLEAN_ATM_DIAG &
148 ,AER_RA_FEEDBACK &
149 ,QC_CU , QI_CU &
150 ,icloud_bl,qc_bl,qi_bl,cldfra_bl &
151 ,PM2_5_DRY, PM2_5_WATER &
152 ,PM2_5_DRY_EC &
153 ,TAUAER300, TAUAER400 &
154 ,TAUAER600, TAUAER999 &
155 ,GAER300, GAER400, GAER600, GAER999 &
156 ,WAER300, WAER400, WAER600, WAER999 &
157 ,TAUAERlw1, TAUAERlw2 &
158 ,TAUAERlw3, TAUAERlw4 &
159 ,TAUAERlw5, TAUAERlw6 &
160 ,TAUAERlw7, TAUAERlw8 &
161 ,TAUAERlw9, TAUAERlw10 &
162 ,TAUAERlw11, TAUAERlw12 &
163 ,TAUAERlw13, TAUAERlw14 &
164 ,TAUAERlw15, TAUAERlw16 &
165 ,progn &
166 ,slope_rad,topo_shading &
167 ,shadowmask,ht,dx,dy,dx2d,area2d &
168 ,dxkm &
169 ,diffuse_frac &
170 ,SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC &
171 ,LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC &
172 ,radtacttime &
173 ,ALSWVISDIR, ALSWVISDIF, ALSWNIRDIR, ALSWNIRDIF &
174 ,SWVISDIR, SWVISDIF, SWNIRDIR, SWNIRDIF &
175 ,SF_SURFACE_PHYSICS, IS_CAMMGMP_USED &
176 ,EXPLICIT_CONVECTION &
177 ,swddir,swddni,swddif &
178 ,swddirc,swddnic &
179 ,swdown_ref,swddir_ref,coszen_ref,Gx,gg,Bx,bb &
180 ,aer_type &
181 ,aer_aod550_opt, aer_aod550_val &
182 ,aer_angexp_opt, aer_angexp_val &
183 ,aer_ssa_opt, aer_ssa_val &
184 ,aer_asy_opt, aer_asy_val &
185 ,aod5502d, angexp2d, aerssa2d, aerasy2d &
186 ,aod5503d &
187 ,cw_rad, shcu_physics &
188 ,obscur,mask,elat_track,elon_track &
189 ,taod5502d, taod5503d &
190 ,mp_physics &
191 ,EFCG,EFCS,EFIG,EFIS,EFSG,aercu_opt &
192 ,EFSS,QS_CU &
193 ,GHG_INPUT &
194 #if (WRF_CHEM == 1)
195 ,chem &
196 ,aod_out &
197 , AOD2D_OUT &
198 , ATOP2D_OUT &
199 ,chem_opt &
200 ,gsfcrad_gocart_coupling &
201 #endif
202 ,feedback_is_ready & ! WRF-CMAQ twoway coupled model
203 ,mass_ws_i & ! WRF-CMAQ twoway coupled model
204 ,mass_ws_j & ! WRF-CMAQ twoway coupled model
205 ,mass_ws_k & ! WRF-CMAQ twoway coupled model
206 ,mass_in_i & ! WRF-CMAQ twoway coupled model
207 ,mass_in_j & ! WRF-CMAQ twoway coupled model
208 ,mass_in_k & ! WRF-CMAQ twoway coupled model
209 ,mass_ec_i & ! WRF-CMAQ twoway coupled model
210 ,mass_ec_j & ! WRF-CMAQ twoway coupled model
211 ,mass_ec_k & ! WRF-CMAQ twoway coupled model
212 ,mass_ss_i & ! WRF-CMAQ twoway coupled model
213 ,mass_ss_j & ! WRF-CMAQ twoway coupled model
214 ,mass_ss_k & ! WRF-CMAQ twoway coupled model
215 ,mass_h2o_i & ! WRF-CMAQ twoway coupled model
216 ,mass_h2o_j & ! WRF-CMAQ twoway coupled model
217 ,mass_h2o_k & ! WRF-CMAQ twoway coupled model
218 ,dgn_i & ! WRF-CMAQ twoway coupled model
219 ,dgn_j & ! WRF-CMAQ twoway coupled model
220 ,dgn_k & ! WRF-CMAQ twoway coupled model
221 ,sig_i & ! WRF-CMAQ twoway coupled model
222 ,sig_j & ! WRF-CMAQ twoway coupled model
223 ,sig_k & ! WRF-CMAQ twoway coupled model
224 ,sw_gtauxar_01 & ! WRF-CMAQ twoway coupled model
225 ,sw_gtauxar_02 & ! WRF-CMAQ twoway coupled model
226 ,sw_gtauxar_03 & ! WRF-CMAQ twoway coupled model
227 ,sw_gtauxar_04 & ! WRF-CMAQ twoway coupled model
228 ,sw_gtauxar_05 & ! WRF-CMAQ twoway coupled model
229 ,sw_ttauxar_01 & ! WRF-CMAQ twoway coupled model
230 ,sw_ttauxar_02 & ! WRF-CMAQ twoway coupled model
231 ,sw_ttauxar_03 & ! WRF-CMAQ twoway coupled model
232 ,sw_ttauxar_04 & ! WRF-CMAQ twoway coupled model
233 ,sw_ttauxar_05 & ! WRF-CMAQ twoway coupled model
234 ,sw_asy_fac_01 & ! WRF-CMAQ twoway coupled model
235 ,sw_asy_fac_02 & ! WRF-CMAQ twoway coupled model
236 ,sw_asy_fac_03 & ! WRF-CMAQ twoway coupled model
237 ,sw_asy_fac_04 & ! WRF-CMAQ twoway coupled model
238 ,sw_asy_fac_05 & ! WRF-CMAQ twoway coupled model
239 ,sw_ssa_01 & ! WRF-CMAQ twoway coupled model
240 ,sw_ssa_02 & ! WRF-CMAQ twoway coupled model
241 ,sw_ssa_03 & ! WRF-CMAQ twoway coupled model
242 ,sw_ssa_04 & ! WRF-CMAQ twoway coupled model
243 ,sw_ssa_05 & ! WRF-CMAQ twoway coupled model
244 ,ozone & ! WRF-CMAQ twoway coupled model
245 ,sw_zbbcddir & ! WRF-CMAQ twoway coupled model
246 ,sw_dirdflux & ! WRF-CMAQ twoway coupled model
247 ,sw_difdflux & ! WRF-CMAQ twoway coupled model
248 )
251 !-------------------------------------------------------------------------
253 ! !USES:
254 USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME &
255 ,RRTMG_LWSCHEME, RRTMG_SWSCHEME &
256 #if( BUILD_RRTMG_FAST == 1)
257 ,RRTMG_LWSCHEME_FAST, RRTMG_SWSCHEME_FAST &
258 #endif
259 #if( BUILD_RRTMK == 1)
260 ,RRTMK_LWSCHEME, RRTMK_SWSCHEME &
261 #endif
262 ,SWRADSCHEME, GSFCSWSCHEME &
263 ,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME &
264 ,HELDSUAREZ &
265 ,goddardswscheme & !NUWRF
266 ,goddardlwscheme & !NUWRF
267 # if (EM_CORE == 1)
268 ,CAMMGMPSCHEME &
269 #if (WRF_CHEM == 1)
271 ,num_chem & !NUWRF
272 ,p_bc1, p_bc2, p_oc1, p_oc2 & !NUWRF
273 ,p_dust_1, p_dust_2, p_dust_3 & !NUWRF
274 ,p_dust_4, p_dust_5 & !NUWRF
275 ,p_sulf, p_seas_1, p_seas_2 & !NUWRF
276 ,p_seas_3, p_seas_4 & !NUWRF
277 #endif
278 #endif
279 ,KFCUPSCHEME, BMJSCHEME & !BSINGH - Added KFCUPSCHEME for WRFCuP scheme
280 ,FLGLWSCHEME, FLGSWSCHEME &
281 ,ECLIPSESCHEME
283 USE module_model_constants
284 USE module_wrf_error , ONLY : wrf_err_message
285 USE module_state_description, ONLY : nuwrf4icescheme
287 ! *** add new modules of schemes here
289 #if ( WRF_CMAQ == 1)
290 USE module_ra_rrtmg_sw_cmaq
291 #endif
293 USE module_ra_sw , ONLY : swrad
294 USE module_ra_gsfcsw , ONLY : gsfcswrad
295 USE module_ra_rrtm , ONLY : rrtmlwrad
296 USE module_ra_rrtmg_lw , ONLY : rrtmg_lwrad, rrtmg_lwinit
297 USE module_ra_rrtmg_sw , ONLY : rrtmg_swrad
298 #if( BUILD_RRTMG_FAST == 1)
299 USE module_ra_rrtmg_lwf , ONLY : rrtmg_lwrad_fast
300 USE module_ra_rrtmg_swf , ONLY : rrtmg_swrad_fast
301 #endif
302 #if( BUILD_RRTMK == 1)
303 USE module_ra_rrtmg_swk , ONLY : rad_rrtmg_driver
304 #endif
305 USE module_ra_cam , ONLY : camrad
306 USE module_ra_gfdleta , ONLY : etara
307 USE module_ra_hs , ONLY : hsrad
309 USE module_ra_goddard , ONLY : goddardrad
310 USE module_ra_flg , ONLY : RAD_FLG
312 USE module_ra_aerosol , ONLY : calc_aerosol_goddard_sw, &
313 calc_aerosol_rrtmg_sw
314 USE module_ra_farms , ONLY : farms_driver
316 ! amontornes-bcodina 2015/09 solar eclipses
317 USE module_ra_eclipse , ONLY : solar_eclipse
319 ! This driver calls subroutines for the radiation parameterizations.
320 !
321 ! short wave radiation choices:
322 ! 1. swrad (19??)
323 ! 4. rrtmg_sw - Added November 2008, MJIacono, AER, Inc.
324 !
325 ! long wave radiation choices:
326 ! 1. rrtmlwrad
327 ! 4. rrtmg_lw - Added November 2008, MJIacono, AER, Inc.
328 !
329 !----------------------------------------------------------------------
330 IMPLICIT NONE
331 !<DESCRIPTION>
332 !
333 ! Radiation_driver is the WRF mediation layer routine that provides the interface to
334 ! to radiation physics packages in the WRF model layer. The radiation
335 ! physics packages to call are chosen by setting the namelist variable
336 ! (Rconfig entry in Registry) to the integer value assigned to the
337 ! particular package (package entry in Registry). For example, if the
338 ! namelist variable ra_lw_physics is set to 1, this corresponds to the
339 ! Registry Package entry for swradscheme. Note that the Package
340 ! names in the Registry are defined constants (frame/module_state_description.F)
341 ! in the CASE statements in this routine.
342 !
343 ! Among the arguments is moist, a four-dimensional scalar array storing
344 ! a variable number of moisture tracers, depending on the physics
345 ! configuration for the WRF run, as determined in the namelist. The
346 ! highest numbered index of active moisture tracers the integer argument
347 ! n_moist (note: the number of tracers at run time is the quantity
348 ! <tt>n_moist - PARAM_FIRST_SCALAR + 1</tt> , not n_moist. Individual tracers
349 ! may be indexed from moist by the Registry name of the tracer prepended
350 ! with P_; for example P_QC is the index of cloud water. An index
351 ! represents a valid, active field only if the index is greater than
352 ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual
353 ! indices for each tracer is defined in module_state_description and
354 ! set in <a href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a> defined in frame/module_configure.F.
355 !
356 ! Physics drivers in WRF 2.0 and higher, originally model-layer
357 ! routines, have been promoted to mediation layer routines and they
358 ! contain OpenMP threaded loops over tiles. Thus, physics drivers
359 ! are called from single-threaded regions in the solver. The physics
360 ! routines that are called from the physics drivers are model-layer
361 ! routines and fully tile-callable and thread-safe.
362 !</DESCRIPTION>
363 !
364 !======================================================================
365 ! Grid structure in physics part of WRF
366 !----------------------------------------------------------------------
367 ! The horizontal velocities used in the physics are unstaggered
368 ! relative to temperature/moisture variables. All predicted
369 ! variables are carried at half levels except w, which is at full
370 ! levels. Some arrays with names (*8w) are at w (full) levels.
371 !
372 !----------------------------------------------------------------------
373 ! In WRF, kms (smallest number) is the bottom level and kme (largest
374 ! number) is the top level. In your scheme, if 1 is at the top level,
375 ! then you have to reverse the order in the k direction.
376 !
377 ! kme - half level (no data at this level)
378 ! kme ----- full level
379 ! kme-1 - half level
380 ! kme-1 ----- full level
381 ! .
382 ! .
383 ! .
384 ! kms+2 - half level
385 ! kms+2 ----- full level
386 ! kms+1 - half level
387 ! kms+1 ----- full level
388 ! kms - half level
389 ! kms ----- full level
390 !
391 !======================================================================
392 ! Grid structure in physics part of WRF
393 !
394 !-------------------------------------
395 ! The horizontal velocities used in the physics are unstaggered
396 ! relative to temperature/moisture variables. All predicted
397 ! variables are carried at half levels except w, which is at full
398 ! levels. Some arrays with names (*8w) are at w (full) levels.
399 !
400 !==================================================================
401 ! Definitions
402 !-----------
403 ! Theta potential temperature (K)
404 ! Qv water vapor mixing ratio (kg/kg)
405 ! Qc cloud water mixing ratio (kg/kg)
406 ! Qr rain water mixing ratio (kg/kg)
407 ! Qi cloud ice mixing ratio (kg/kg)
408 ! Qs snow mixing ratio (kg/kg)
409 ! QCCONV convective cloud mixing ratio (kg/kg)
410 ! QICONV convective ice mixing ratio (kg/kg)
411 !-----------------------------------------------------------------
412 !-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3)
413 !-- PM2_5_WATER PM2.5 water mass (ug m^-3)
414 !-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3)
415 !-- RTHRATEN Theta tendency
416 ! due to radiation (K/s)
417 !-- RTHRATENLW Theta tendency
418 ! due to long wave radiation (K/s)
419 !-- RTHRATENLWC Theta tendency
420 ! due to clear-sky long wave radiation (K/s)
421 !-- RTHRATENSW Theta temperature tendency
422 ! due to short wave radiation (K/s)
423 !-- RTHRATENSWC Theta tendency
424 ! due to clear-sky short wave radiation (K/s)
425 !-- dt time step (s)
426 !-- itimestep number of time steps
427 !-- GLW downward long wave flux at ground surface (W/m^2)
428 !-- GSW net short wave flux at ground surface (W/m^2)
429 !-- SWDOWN downward short wave flux at ground surface (W/m^2)
430 !-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ)
431 !-- RLWTOA upward long wave at top of atmosphere (w/m2)
432 !-- RSWTOA upward short wave at top of atmosphere (w/m2)
433 !-- XLAT latitude, south is negative (degree)
434 !-- XLONG longitude, west is negative (degree)
435 !-- ALBEDO albedo (between 0 and 1)
436 !-- CLDFRA cloud fraction (between 0 and 1)
437 !-- CLDFRA_DP cloud fraction from deep cloud in a cumulus scheme
438 !-- CLDFRA_SH cloud fraction from shallow cloud in a cumulus scheme
439 !-- CLDFRA_MP_ALL cloud fraction from CAMMGMP microphysics scheme
440 !-- CCLDFRA convective cloud fraction (between 0 and 1)
441 !-- EMISS surface emissivity (between 0 and 1)
442 !-- rho_phy density (kg/m^3)
443 !-- rr dry air density (kg/m^3)
444 !-- moist moisture array (4D - last index is species) (kg/kg)
445 !-- n_moist number of moisture species
446 !-- qndrop Cloud droplet number (#/kg)
447 !-- p8w pressure at full levels (Pa)
448 !-- p_phy pressure (Pa)
449 !-- Pb base-state pressure (Pa)
450 !-- pi_phy exner function (dimensionless)
451 !-- dz8w dz between full levels (m)
452 !-- t_phy temperature (K)
453 !-- t8w temperature at full levels (K)
454 !-- GMT Greenwich Mean Time Hour of model start (hour)
455 !-- JULDAY the initial day (Julian day)
456 !-- RADT time for calling radiation (min)
457 !-- ra_call_offset -1 (old) means usually just before output, 0 after
458 !-- DEGRAD conversion factor for
459 ! degrees to radians (pi/180.) (rad/deg)
460 !-- DPD degrees per day for earth's
461 ! orbital position (deg/day)
462 !-- R_d gas constant for dry air (J/kg/K)
463 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
464 !-- G acceleration due to gravity (m/s^2)
465 !-- rvovrd R_v divided by R_d (dimensionless)
466 !-- XTIME time since simulation start (min)
467 !-- DECLIN solar declination angle (rad)
468 !-- SOLCON solar constant (W/m^2)
469 !-- ids start index for i in domain
470 !-- ide end index for i in domain
471 !-- jds start index for j in domain
472 !-- jde end index for j in domain
473 !-- kds start index for k in domain
474 !-- kde end index for k in domain
475 !-- ims start index for i in memory
476 !-- ime end index for i in memory
477 !-- jms start index for j in memory
478 !-- jme end index for j in memory
479 !-- kms start index for k in memory
480 !-- kme end index for k in memory
481 !-- i_start start indices for i in tile
482 !-- i_end end indices for i in tile
483 !-- j_start start indices for j in tile
484 !-- j_end end indices for j in tile
485 !-- kts start index for k in tile
486 !-- kte end index for k in tile
487 !-- num_tiles number of tiles
488 !
489 !==================================================================
490 !
491 LOGICAL, OPTIONAL, INTENT(IN) :: explicit_convection
492 LOGICAL,INTENT(IN) :: bmj_rad_feedback
494 LOGICAL :: expl_conv
495 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
496 ims,ime, jms,jme, kms,kme, &
497 kts,kte, &
498 num_tiles
500 INTEGER, INTENT(IN) :: lw_physics, sw_physics, mp_physics, sw_eclipse
501 INTEGER, INTENT(IN) :: ghg_input
502 INTEGER, INTENT(IN) :: o3input, aer_opt
503 INTEGER, INTENT(IN) :: id
504 integer, intent(in) :: swint_opt
505 integer, intent(in), OPTIONAL :: solar_opt
506 integer :: solar_opt_local
507 ! WRF-Solar EPS
508 integer, intent (in) :: multi_perturb
509 logical, intent (in) :: pert_farms, pert_cld3, couple_farms
510 real, intent (in) :: pert_farms_albedo, pert_farms_aod, pert_farms_angexp, pert_farms_aerasy, &
511 pert_farms_qv, pert_farms_qc, pert_farms_qs, pert_cld3_qv, pert_cld3_t
512 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_albedo, perts_aod, &
513 perts_angstrom, perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow, perts_th
515 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
516 i_start,i_end,j_start,j_end
518 INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset
519 INTEGER, INTENT(IN ) :: cldovrlp, idcor ! J. Henderson AER
520 INTEGER, INTENT(IN ) :: alevsiz, no_src_types
521 INTEGER, INTENT(IN ) :: levsiz, n_ozmixm
522 INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2
523 REAL, INTENT(IN ) :: cam_abs_freq_s
525 LOGICAL, INTENT(IN ) :: warm_rain
526 INTEGER, INTENT(IN ) :: cu_physics !CuP, wig 5-Oct-2006 !BSINGH - For WRFCuP scheme
527 !BSINGH - For WRFCuP scheme
528 LOGICAL, OPTIONAL, INTENT(IN) :: shallowcu_forced_ra !CuP, wig
529 !BSINGH -ENDS
530 LOGICAL, INTENT(IN ) :: is_CAMMGMP_used !BSINGH:01/31/2013: Added for CAM5 RRTMG
532 REAL, INTENT(IN ) :: RADT
534 REAL, DIMENSION( ims:ime, jms:jme ), &
535 INTENT(IN ) :: XLAND, &
536 XICE, &
537 TSK, &
538 VEGFRA, &
539 SNOW
540 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, &
541 INTENT(IN ) :: OZMIXM
542 REAL, DIMENSION( ims:ime, alevsiz, jms:jme, no_src_types, n_ozmixm-1 ), OPTIONAL, &
543 INTENT(IN ) :: AERODM
544 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, no_src_types ), OPTIONAL, &
545 INTENT(INOUT) :: AEROD
546 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
547 INTENT(INOUT) :: AODTOT
549 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ) :: OZFLG
551 REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN
552 REAL, DIMENSION(alevsiz), OPTIONAL, INTENT(IN ) :: PINA
554 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2
555 REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, &
556 INTENT(IN ) :: aerosolc_1, aerosolc_2
557 REAL, DIMENSION(paerlev), OPTIONAL, &
558 INTENT(IN ) :: m_hybi0
560 REAL, DIMENSION( ims:ime, jms:jme ), &
561 INTENT(INOUT) :: HTOP, &
562 HBOT, &
563 HTOPR, &
564 HBOTR, &
565 CUPPT
567 !BSINGH - For WRFCuP scheme
568 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
569 INTENT(INOUT) :: &
570 cutop, & !CuP, wig 1-Oct-2006
571 cubot, & !CuP, wig 1-Oct-2006
572 shall !CuP, wig 4-Feb-2008
573 !BSINGH -ENDS
576 INTEGER, INTENT(IN ) :: julyr
577 !
578 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
579 INTENT(IN ) :: dz8w, &
580 z, &
581 p8w, &
582 p, &
583 pi, &
584 t, &
585 t8w, &
586 rho
588 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
589 OPTIONAL, &
590 INTENT(IN ) :: cw_rad
591 INTEGER, OPTIONAL, INTENT(IN) :: shcu_physics
593 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
594 INTENT(IN) :: EFCG, &
595 EFCS, &
596 EFIG, &
597 EFIS, &
598 EFSG, &
599 EFSS
601 !BSINGH - For WRFCuP scheme
602 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
603 INTENT(INOUT ) :: cldfra_cup !CuP, wig 1-Oct-2006
606 !BSINGH -ENDS
608 !
609 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
610 INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
611 gaer300,gaer400,gaer600,gaer999, & ! jcb
612 waer300,waer400,waer600,waer999
614 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
615 INTENT(IN ) :: qc_cu, qi_cu
617 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
618 INTENT(IN ) :: qc_bl, qi_bl, qs_cu
620 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
621 INTENT(IN ) :: tauaerlw1,tauaerlw2,tauaerlw3,tauaerlw4, & ! czhao
622 tauaerlw5,tauaerlw6,tauaerlw7,tauaerlw8, & ! czhao
623 tauaerlw9,tauaerlw10,tauaerlw11,tauaerlw12, & ! czhao
624 tauaerlw13,tauaerlw14,tauaerlw15,tauaerlw16
626 INTEGER, INTENT(IN) :: icloud_cu
628 INTEGER, INTENT(IN ), OPTIONAL :: icloud_bl
629 INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
630 INTEGER, INTENT(IN ) :: calc_clean_atm_diag
632 !jdfcz INTEGER, OPTIONAL, INTENT(IN ) :: progn,prescribe
633 INTEGER, OPTIONAL, INTENT(IN ) :: progn
634 !
635 ! variables for aerosols (only if running with chemistry)
636 !
637 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
638 INTENT(IN ) :: pm2_5_dry, &
639 pm2_5_water, &
640 pm2_5_dry_ec
641 !
642 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
643 INTENT(INOUT) :: RTHRATEN, &
644 RTHRATENLW, &
645 RTHRATENLWC,&
646 RTHRATENSW, &
647 RTHRATENSWC
649 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
650 ! INTENT(INOUT) :: SWUP, &
651 ! SWDN, &
652 ! SWUPCLEAR, &
653 ! SWDNCLEAR, &
654 ! LWUP, &
655 ! LWDN, &
656 ! LWUPCLEAR, &
657 ! LWDNCLEAR
659 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
660 ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
661 ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
662 ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
663 ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
665 ! TOA and surface, upward and downward, total, clear (no cloud), and clean (no aerosol) fluxes
666 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
667 SWUPT, SWUPTC, SWUPTCLN, SWDNT, SWDNTC, SWDNTCLN,&
668 SWUPB, SWUPBC, SWUPBCLN, SWDNB, SWDNBC, SWDNBCLN,&
669 LWUPT, LWUPTC, LWUPTCLN, LWDNT, LWDNTC, LWDNTCLN,&
670 LWUPB, LWUPBC, LWUPBCLN, LWDNB, LWDNBC, LWDNBCLN
673 ! Upward and downward, total and clear sky layer fluxes (W m-2)
674 REAL, DIMENSION( ims:ime, kms:kme+2, jms:jme ), &
675 OPTIONAL, INTENT(INOUT) :: &
676 SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC, &
677 LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC
679 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
680 INTENT(INOUT) :: SWCF, &
681 LWCF, &
682 OLR
683 ! ---- fds (06/2010) ssib alb components ------------
684 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
685 INTENT(IN ) :: ALSWVISDIR, &
686 ALSWVISDIF, &
687 ALSWNIRDIR, &
688 ALSWNIRDIF
689 ! ---- fds (06/2010) ssib swr components ------------
690 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
691 INTENT(OUT ) :: SWVISDIR, &
692 SWVISDIF, &
693 SWNIRDIR, &
694 SWNIRDIF
695 INTEGER, OPTIONAL, INTENT(IN ) :: sf_surface_physics
696 !
697 REAL, DIMENSION( ims:ime, jms:jme ), &
698 INTENT(IN ) :: XLAT, &
699 XLONG, &
700 EMISS
701 real, dimension (ims:ime, jms:jme), intent(inout) :: albedo
702 !
703 REAL, DIMENSION( ims:ime, jms:jme ), &
704 INTENT(INOUT) :: GSW, &
705 GLW
707 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN
708 ! PAJ: FARMS coupling
709 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT), OPTIONAL :: SWDOWN2, SWDDNI2, SWDDIF2, SWDDIR2, SWDOWNC2, SWDDNIC2
711 ! ------------------------------------------------------------------------------ jararias 2013/08/10 -----------
712 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: swddir, & ! All-sky SW broadband surface direct irradiance
713 swddni, & ! All-sky SW broadband surface direct normal irradiance
714 swddif ! All-sky SW broadband surface diffuse irradiance
715 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: Gx,Bx,gg,bb, & ! For SW sza-interpolation
716 swdown_ref, &
717 swddir_ref, &
718 coszen_ref
719 ! ------------------------------------------------------------------------------ jararias 2013/11 -----------
720 INTEGER, INTENT(IN) :: aer_type, & ! rural, urban, maritime, ...
721 aer_aod550_opt, & ! input option for AOD at 550 nm
722 aer_angexp_opt, & ! input option for aerosol Angstrom exponent
723 aer_ssa_opt, & ! input option for aerosol ssa
724 aer_asy_opt, & ! input option for aerosol asy
725 aercu_opt !
726 REAL, INTENT(IN) :: aer_aod550_val, & ! AOD at 550 nm if aer_aod550_opt=1
727 aer_angexp_val, & ! aerosol Angstrom exponent if aer_angexp_opt=1
728 aer_ssa_val, & ! aerosol ssa if aer_ssa_opt=1
729 aer_asy_val ! aerosol asy if aer_asy_opt=1
730 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
731 INTENT(INOUT) :: aod5502d, & ! gridded AOD at 550 nm from auxinput
732 angexp2d, & ! gridded aerosol Angstrom exponent from auxinput
733 aerssa2d, & ! gridded aerosol ssa from auxinput
734 aerasy2d ! gridded aerosol asy from auxinput
735 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
736 INTENT(OUT) :: aod5503d ! 3D AOD at 550 nm
738 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL:: taod5503d ! Trude
739 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL:: taod5502d ! Trude
740 !
741 REAL, INTENT(IN ) :: GMT,dt, &
742 julian, xtime
743 INTEGER, INTENT(IN ),OPTIONAL :: YR
744 !
745 INTEGER, INTENT(IN ) :: JULDAY, itimestep
746 REAL, INTENT(IN ),OPTIONAL :: CURR_SECS
747 LOGICAL, INTENT(IN ),OPTIONAL :: ADAPT_STEP_FLAG
749 INTEGER,INTENT(IN) :: NPHS
750 REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: &
751 CFRACH, & !Added
752 CFRACL, & !Added
753 CFRACM, & !Added
754 CZMEAN !Added
755 REAL, DIMENSION( ims:ime, jms:jme ), &
756 INTENT(INOUT) :: &
757 RLWTOA, & !Added
758 RSWTOA, & !Added
759 ACFRST, & !Added
760 ACFRCV !Added
762 INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: &
763 NCFRST, & !Added
764 NCFRCV !Added
766 ! NUWRF JJS 20101021 vvvvv
767 ! for inline Gocart coupling
768 #if( WRF_CHEM == 1)
770 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem), &
771 INTENT(IN) :: chem
772 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(out) :: aod_out !Aeorosol Optical Depth
774 real, dimension( ims:ime, jms:jme ), OPTIONAL, INTENT(OUT) :: & !goddardrad
775 aod2d_out ,& ! column aerosol optical depth
776 atop2d_out ! aerosol top pressure [mb]
779 integer :: i24h
780 INTEGER, PARAMETER :: num_go = 14 ! number of the gocart aerosol species
781 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_go) :: aero
782 REAL, PARAMETER :: frac(4)=(/ 0.01053,0.08421,0.25263,0.65263 /) !fraction for fine dust
783 integer,intent(in) :: chem_opt ! EMK
784 integer,intent(in) :: gsfcrad_gocart_coupling ! EMK
786 #endif
787 ! NUWRF JJS 20101021 ^^^^^
789 ! JJS 20090623 vvvvv
790 ! Optional, only for Goddard LW and SW
791 REAL, DIMENSION(IMS:IME, JMS:JME, 1:8) :: ERBE_out !extra output for SDSU
792 REAL, DIMENSION(IMS:IME, JMS:JME), OPTIONAL, INTENT(INOUT) :: & !BSINGH(PNNL)- Lahey compiler forced this specification to be intent-inout
793 TLWDN, TLWUP, &
794 SLWDN, SLWUP, &
795 TSWDN, TSWUP, &
796 SSWDN, SSWUP ! for Goddard schemes
797 ! NUWRF JJS 20090623 ^^^^^
799 ! NUWRF JJS 20140225 vvvvv
800 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
801 INTENT(INOUT) :: re_cloud_gsfc, re_rain_gsfc, re_ice_gsfc, &
802 re_snow_gsfc, re_graupel_gsfc, re_hail_gsfc
803 ! NUWRF JJS 20140225 ^^^^^
805 real, dimension( ims:ime, jms:jme, 1:4 ) :: sflxd !NUWRF SW only for LIS
807 ! REAL, DIMENSION(IMS:IME, JMS:JME, 1:4) :: flxd !NUWRF extra radiation output for LIS (CLM)
808 ! 1-surface downward UV+VIS beam radiation [W/m2]
809 ! 2-surface downward UV+VIS diffuse radiation [W/m2]
810 ! 3-surface downward NIR beam radiation [W/m2]
811 ! 4-surface downward NIR diffuse radiation [W/m2]
814 real, dimension( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) :: & !goddardrad
815 cod2d_out ,& ! column optical depth
816 ctop2d_out ! cloud top pressure [mb]
818 ! Added by ZCX for low and total cloud fraction
819 REAL, DIMENSION( kms:kme ), OPTIONAL, INTENT(IN) :: znu ! eta values on half (mass)levels
820 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) :: &
821 cldt
823 ! Optional (only used by CAM lw scheme)
825 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,&
826 INTENT(INOUT) :: abstot
827 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,&
828 INTENT(INOUT) :: absnxt
829 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
830 INTENT(INOUT) :: emstot
832 !
833 ! Optional
834 !
835 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
836 OPTIONAL, &
837 INTENT(INOUT) :: CLDFRA, &
838 CCLDFRA,&
839 QCCONV, &
840 QICONV
842 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & ! ckay for sub-grid cloud fraction
843 OPTIONAL, &
844 INTENT(INOUT) :: cldfra_dp, &
845 cldfra_sh, &
846 cldfra_bl
848 !..G. Thompson
849 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT):: re_cloud, re_ice, re_snow
850 INTEGER, INTENT(INOUT):: has_reqc, has_reqi, has_reqs
852 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
853 OPTIONAL, &
854 INTENT(IN ) :: &
855 F_ICE_PHY, &
856 F_RAIN_PHY, &
857 CLDFRA_MP_ALL
859 #if (EM_CORE == 1)
860 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
861 OPTIONAL, &
862 INTENT(IN ) :: &
863 LRADIUS, &
864 IRADIUS
865 #endif
867 REAL, DIMENSION( ims:ime, jms:jme ), &
868 OPTIONAL, &
869 INTENT(OUT) :: SWDOWNC, SWDDIRC, SWDDNIC
870 !
871 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
872 OPTIONAL, &
873 INTENT(INOUT ) :: &
874 pb &
875 ,qv,qc,qr,qi,qs,qg,qh,qndrop, &
876 qnifa,qnwfa,qnbca, & ! Trude
877 qi2,qi3 ! for P3
879 LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qh,f_qndrop,&
880 f_qnifa,f_qnwfa, & ! trude
881 f_qi2,f_qi3 ! for P3
882 INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt
883 ! Solar diag
884 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT), OPTIONAL :: qc_tot, qi_tot
885 !shbaek
886 real, dimension ( ims:ime, kms:kme, jms:jme ), optional, intent(in) :: qnc_curr
887 LOGICAL, OPTIONAL :: f_qnc
888 !
889 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
890 OPTIONAL, &
891 INTENT(INOUT) :: taucldi,taucldc
893 REAL, OPTIONAL, INTENT(IN) :: dxkm
895 ! Variables for slope-dependent radiation
897 REAL, OPTIONAL, INTENT(IN) :: dx,dy
898 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN), OPTIONAL :: dx2d, area2d
899 INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,topo_shading
900 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: ht
901 INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: shadowmask
902 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(OUT) :: diffuse_frac
904 REAL , OPTIONAL, INTENT(INOUT) :: radtacttime ! Storing the time in s when radiation is called next
905 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
906 INTENT(INOUT) :: o3rad
908 ! begin WRF-CMAQ coupled model block
909 REAL, DIMENSION (ims:ime, kms:kme, jms:jme ), optional, &
910 INTENT(INOUT) :: mass_ws_i, mass_ws_j, mass_ws_k, &
911 mass_in_i, mass_in_j, mass_in_k, &
912 mass_ec_i, mass_ec_j, mass_ec_k, &
913 mass_ss_i, mass_ss_j, mass_ss_k, &
914 mass_h2o_i, mass_h2o_j, mass_h2o_k, &
915 dgn_i, dgn_j, dgn_k, &
916 sig_i, sig_j, sig_k
918 REAL, DIMENSION (ims:ime, kms:kme, jms:jme ), optional, INTENT(OUT) :: sw_gtauxar_01, &
919 sw_gtauxar_02, &
920 sw_gtauxar_03, &
921 sw_gtauxar_04, &
922 sw_gtauxar_05, &
923 sw_asy_fac_01, &
924 sw_asy_fac_02, &
925 sw_asy_fac_03, &
926 sw_asy_fac_04, &
927 sw_asy_fac_05, &
928 sw_ssa_01, &
929 sw_ssa_02, &
930 sw_ssa_03, &
931 sw_ssa_04, &
932 sw_ssa_05
934 REAL, DIMENSION( ims:ime, jms:jme ), optional, INTENT(OUT) :: sw_ttauxar_01, &
935 sw_ttauxar_02, &
936 sw_ttauxar_03, &
937 sw_ttauxar_04, &
938 sw_ttauxar_05
940 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), optional, INTENT(IN) :: ozone
941 REAL, DIMENSION( ims:ime, jms:jme ), optional, INTENT(OUT) :: sw_zbbcddir, &
942 sw_dirdflux, &
943 sw_difdflux
945 LOGICAL, INTENT(IN) :: feedback_is_ready
946 ! end WRF-CMAQ coupled model block
948 ! vert nesting
949 REAL, OPTIONAL , INTENT(IN ) :: p_top
950 REAL :: p_top_dummy
952 ! LOCAL VAR
953 INTEGER, DIMENSION( ims:ime, kms:kme, jms:jme ) :: cldfra1_flag
954 REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON
955 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS
956 REAL, DIMENSION( ims:ime, jms:jme ) :: coszr
957 REAL, DIMENSION( ims:ime, levsiz, jms:jme ) :: ozmixt
958 REAL, DIMENSION( ims:ime, alevsiz, jms:jme, 1:no_src_types ) :: aerodt
960 REAL :: DECLIN,SOLCON,XXLAT,TLOCTM,XT24, CEN_LAT, cldfra_cup_mod
961 INTEGER :: i,j,k,its,ite,jts,jte,ij
962 INTEGER :: STEPABS
963 LOGICAL :: gfdl_lw,gfdl_sw, compute_cldfra_cup
964 LOGICAL :: doabsems
965 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
966 INTEGER :: s
967 REAL :: ITRMX, &
968 ITRMN
969 REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, &
970 DJUL,RJUL,ECCFAC
971 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qc_temp
972 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save
973 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qs_save
975 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qc_cu_weight
976 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_cu_weight
977 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qs_cu_weight
979 REAL :: gridkm, Wice,Wh2o
980 REAL, DIMENSION(kms:kme):: t_1d, p_1d, Dz_1d, qv_1d, qc_1d, qi_1d, qs_1d, cf_1d
982 REAL :: next_rad_time, DTaccum
983 LOGICAL :: run_param , doing_adapt_dt , decided
984 LOGICAL :: flg_lw, flg_sw
985 !ZCX
986 REAL :: cldji,cldlji
987 !ckay
988 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: cldfra_cu
989 !------------------------------------------------------------------
990 ! solar related variables are added to declaration
991 !-------------------------------------------------
992 REAL, OPTIONAL, INTENT(OUT) :: DECLINX,SOLCONX
993 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZEN
994 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: HRANG
995 !------------------------------------------------------------------
997 ! jararias, 2013/08/10
998 real :: ioh,kt,airmass,kd
999 real, dimension(ims:ime,jms:jme) :: coszen_loc,hrang_loc
1000 ! jararias 2013/11 but modified on 2016/2/10
1001 ! the following three arrays may be dimensioned by (ims,ime,kms,kme,jms,jme,aerosol_vars)
1002 real, dimension(:,:,:,:), pointer :: tauaer_sw=>null(), ssaaer_sw=>null(), asyaer_sw=>null()
1003 !
1004 ! Montornes-Cordian eclipse variables
1005 real, dimension(ims:ime,jms:jme), optional, INTENT(OUT) :: obscur
1006 integer, dimension(ims:ime,jms:jme), optional, INTENT(OUT) :: mask
1007 real, optional, INTENT(OUT) :: elat_track, elon_track
1008 ! Local variables
1009 INTEGER, DIMENSION( ims:ime, jms:jme ) :: mask_loc
1010 REAL, DIMENSION( ims:ime, jms:jme ) :: obscur_loc
1011 REAL :: elat_loc, elon_loc
1013 ! Trude AOD variables
1014 INTEGER, PARAMETER:: taer_type = 1 ! rural, urban, maritime, ...
1015 INTEGER, PARAMETER:: taer_aod550_opt = 2 ! input option for AOD at 550 nm
1016 INTEGER, PARAMETER:: taer_angexp_opt = 3 ! input option for aerosol Angstrom exponent
1017 INTEGER, PARAMETER:: taer_ssa_opt = 3 ! input option for aerosol ssa
1018 INTEGER, PARAMETER:: taer_asy_opt = 3 ! input option for aerosol asy
1020 ! WRF-Solar EPS
1021 real, dimension(:, :, :), allocatable :: qv_tmp, qc_tmp, qs_tmp
1024 !---------- local test vars
1025 ! real, dimension(ims:ime, kms:kme, jms:jme) :: hlw, hsw
1027 LOGICAL :: proceed_cmaq_sw
1029 logical, save :: firstime = .true.
1030 logical, save :: feedback_restart, direct_sw_feedback
1032 #if ( WRF_CMAQ == 1 )
1033 if (firstime) then
1034 firstime = .false.
1035 CALL nl_get_direct_sw_feedback ( .false. , direct_sw_feedback )
1036 CALL nl_get_feedback_restart ( .false. , feedback_restart )
1037 end if
1038 #else
1039 direct_sw_feedback = .false.
1040 feedback_restart = .false.
1041 #endif
1043 ! This allows radiation schemes to correctly
1044 ! interface with the convection scheme when convection is only
1045 ! enabled in some domains:
1046 if(present(explicit_convection)) then
1047 expl_conv=explicit_convection
1048 else
1049 expl_conv=.true. ! backward compatibility for ARW
1050 endif
1052 IF ( ICLOUD == 3 ) THEN
1053 IF (PRESENT(dxkm)) then
1054 gridkm = 1.414*SQRT(dxkm*dxkm + dy*0.001*dy*0.001)
1055 ELSE IF (PRESENT(dx)) then
1056 gridkm = SQRT(dx*0.001*dx*0.001 + dy*0.001*dy*0.001)
1057 endif
1059 if (itimestep .LE. 100) then
1060 WRITE ( wrf_err_message , * ) 'Grid spacing in km ', dx, dy, gridkm
1061 CALL wrf_debug (100, wrf_err_message)
1062 endif
1063 END IF
1065 CALL wrf_debug (1, 'Top of Radiation Driver')
1066 ! WRITE ( wrf_err_message , * ) 'itimestep = ',itimestep,', dt = ',dt,', lw and sw options = ',lw_physics,sw_physics
1067 ! CALL wrf_debug (1, wrf_err_message )
1068 if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return
1070 ! amontornes-bcodina (2014-05-02) :: improving the namelist settings consistency for the FLG scheme
1071 ! if (lw_physics .ne. FLGLWSCHEME .and. sw_physics .eq. FLGSWSCHEME) then
1072 ! call wrf_error_fatal('SW and LW schemes are in conflict. SW is FLG and LW is a different scheme!')
1073 ! end if
1074 ! if (lw_physics .eq. FLGLWSCHEME .and. sw_physics .ne. FLGSWSCHEME) then
1075 ! call wrf_error_fatal('SW and LW schemes are in conflict. LW is FLG and SW is a different scheme!')
1076 ! end if
1078 ! ra_call_offset = -1 gives old method where radiation may be called just before output
1079 ! ra_call_offset = 0 gives new method where radiation may be called just after output
1080 ! and is also consistent with removal of offset in new XTIME
1081 ! also need to account for stepra=1 which always has zero modulo output
1083 doing_adapt_dt = .FALSE.
1084 IF ( PRESENT(adapt_step_flag) ) THEN
1085 IF ( adapt_step_flag ) THEN
1086 doing_adapt_dt = .TRUE.
1087 IF ( radtacttime .eq. 0. ) THEN
1088 radtacttime = CURR_SECS + radt*60.
1089 END IF
1090 END IF
1091 END IF
1093 ! Do we run through this scheme or not?
1095 ! Test 1: If this is the initial model time, then yes.
1096 ! ITIMESTEP=1
1097 ! Test 2: If the user asked for the radiation to be run every time step, then yes.
1098 ! RADT=0 or STEPRA=1
1099 ! Test 3: If not adaptive dt, and this is on the requested radiation frequency, then yes.
1100 ! MOD(ITIMESTEP,STEPRA)=0 (or 1, depending on the offset setting)
1101 ! Test 4: If using adaptive dt and the current time is past the last requested activate radiation time, then yes.
1102 ! CURR_SECS >= RADTACTTIME
1104 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
1105 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
1106 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
1108 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
1109 ! radiation run.
1111 run_param = .FALSE.
1112 decided = .FALSE.
1113 IF ( ( .NOT. decided ) .AND. &
1114 ( itimestep .EQ. 1 ) ) THEN
1115 run_param = .TRUE.
1116 decided = .TRUE.
1117 END IF
1119 IF ( ( .NOT. decided ) .AND. &
1120 ( ( radt .EQ. 0. ) .OR. ( stepra .EQ. 1 ) ) ) THEN
1121 run_param = .TRUE.
1122 decided = .TRUE.
1123 END IF
1125 IF ( ( .NOT. decided ) .AND. &
1126 ( .NOT. doing_adapt_dt ) .AND. &
1127 ( MOD(itimestep,stepra) .EQ. 1+ra_call_offset ) ) THEN
1128 run_param = .TRUE.
1129 decided = .TRUE.
1130 END IF
1132 IF ( ( .NOT. decided ) .AND. &
1133 ( doing_adapt_dt ) .AND. &
1134 ( curr_secs .GE. radtacttime ) ) THEN
1135 run_param = .TRUE.
1136 decided = .TRUE.
1137 radtacttime = curr_secs + radt*60
1138 END IF
1140 IF ( mp_physics .EQ. nuwrf4icescheme ) THEN
1141 DO ij = 1 , num_tiles
1142 its = i_start(ij)
1143 ite = i_end(ij)
1144 jts = j_start(ij)
1145 jte = j_end(ij)
1146 DO j=jts,jte
1147 DO k=kts,kte
1148 DO i=its,ite
1149 re_cloud(i,k,j) = re_cloud_gsfc(i,k,j) * 1.E-6
1150 re_ice(i,k,j) = re_ice_gsfc(i,k,j) * 1.E-6
1151 re_snow(i,k,j) = re_snow_gsfc(i,k,j) * 1.E-6
1152 ENDDO
1153 ENDDO
1154 ENDDO
1155 ENDDO
1156 END IF
1158 if(swint_opt.eq.1 .or. swint_opt == 2) then
1159 DO ij = 1 , num_tiles
1160 its = i_start(ij)
1161 ite = i_end(ij)
1162 jts = j_start(ij)
1163 jte = j_end(ij)
1164 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
1165 DEGRAD,DPD )
1166 call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, &
1167 julian,xtime,gmt,declin,degrad, &
1168 xlong,xlat,coszen_loc,hrang_loc)
1169 end do
1170 end if
1172 solar_opt_local = 0
1173 IF ( PRESENT(solar_opt) ) THEN
1174 solar_opt_local = solar_opt
1175 END IF
1176 Solar_step: IF (run_param .or. solar_opt_local == 1 .or. swint_opt == 2) THEN
1178 !---------------
1179 ! Calculate constant for short wave radiation
1180 ! moved up and out of OMP loop because it only needs to be computed once
1181 ! and because it is not entirely thread-safe (XT24, TOLOCTM and XXLAT need
1182 ! their thread-privacy) JM 20100217
1183 DO ij = 1 , num_tiles
1184 its = i_start(ij)
1185 ite = i_end(ij)
1186 jts = j_start(ij)
1187 jte = j_end(ij)
1188 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
1189 DEGRAD,DPD )
1191 ! amontornes-bcodina 2015/09 solar eclipses
1192 ! Solar eclipse prediction based on the Bessel's method
1193 ! outputs: obscur, mask, elat_track, elon_track
1194 CALL solar_eclipse(ims,ime,jms,jme,its,ite,jts,jte, &
1195 julian,gmt,YR,xtime,radt, &
1196 degrad,xlong,xlat,obscur_loc,mask_loc, &
1197 elat_loc,elon_loc,sw_eclipse )
1199 IF(PRESENT(declinx).AND.PRESENT(solconx))THEN
1200 ! saved to state arrays used in surface driver
1201 declinx=declin
1202 solconx=solcon
1203 ENDIF
1204 ! added coszen subroutine : jararias, 2013/08/10
1205 ! outputs: coszen, hrang
1206 call calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, &
1207 julian,xtime+radt*0.5,gmt, &
1208 declin,degrad,xlong,xlat,coszen,hrang)
1210 ! backup the incoming hydrometeors
1212 IF ( F_QC ) THEN
1213 DO j=jts,jte
1214 DO k=kts,kte
1215 DO i=its,ite
1216 qc_save(i,k,j) = qc(i,k,j)
1217 ENDDO
1218 ENDDO
1219 ENDDO
1220 ENDIF
1221 IF ( F_QI ) THEN
1222 DO j=jts,jte
1223 DO k=kts,kte
1224 DO i=its,ite
1225 qi_save(i,k,j) = qi(i,k,j)
1226 ENDDO
1227 ENDDO
1228 ENDDO
1229 ENDIF
1231 IF(aercu_opt.gt.0.0) THEN
1232 IF ( F_QI ) THEN
1233 DO j=jts,jte
1234 DO k=kts,kte
1235 DO i=its,ite
1236 qs_save(i,k,j) = qs(i,k,j)
1237 ENDDO
1238 ENDDO
1239 ENDDO
1240 ENDIF
1241 END IF
1243 ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
1244 if( F_QC ) then
1245 DO j=jts,jte
1246 DO k=kts,kte
1247 DO i=its,ite
1248 qc_temp(I,K,J)=qc(I,K,J)
1249 ENDDO
1250 ENDDO
1251 ENDDO
1252 else
1253 DO j=jts,jte
1254 DO k=kts,kte
1255 DO i=its,ite
1256 qc_temp(I,K,J)=0.
1257 ENDDO
1258 ENDDO
1259 ENDDO
1260 endif
1261 !
1262 ! temporarily modify hydrometeors (this is for GD scheme and WRF-Chem)
1263 !
1264 IF ( F_QC .AND. icloud_cu .EQ. 1 ) THEN
1265 DO j=jts,jte
1266 DO k=kts,kte
1267 DO i=its,ite
1268 qc(i,k,j) = qc(i,k,j) + qc_cu(i,k,j)
1269 ENDDO
1270 ENDDO
1271 ENDDO
1272 ENDIF
1273 IF ( F_QI .AND. icloud_cu .EQ. 1 ) THEN
1274 DO j=jts,jte
1275 DO k=kts,kte
1276 DO i=its,ite
1277 qi(i,k,j) = qi(i,k,j) + qi_cu(i,k,j)
1278 ENDDO
1279 ENDDO
1280 ENDDO
1281 ENDIF
1283 #if (EM_CORE == 1)
1284 ! temporarily modify hydrometeors (for P3, if 2 cat then add ice from both categories)
1285 !
1286 IF ( F_QI2) THEN
1287 DO j=jts,jte
1288 DO k=kts,kte
1289 DO i=its,ite
1290 qi(i,k,j) = qi(i,k,j) + qi2(i,k,j)
1291 ENDDO
1292 ENDDO
1293 ENDDO
1294 ENDIF
1296 ! for Jensen ISHMAEL, add the third ice species
1297 IF ( F_QI3) THEN
1298 DO j=jts,jte
1299 DO k=kts,kte
1300 DO i=its,ite
1301 qi(i,k,j) = qi(i,k,j) + qi3(i,k,j)
1302 ENDDO
1303 ENDDO
1304 ENDDO
1305 ENDIF
1306 #endif
1308 ! Choose how to compute cloud fraction (since 3.6)
1309 ! Initialize to zero
1310 DO j=jts,jte
1311 DO k=kts,kte
1312 DO i=its,ite
1313 CLDFRA(i,k,j) = 0.
1314 END DO
1315 END DO
1316 END DO
1317 !---------------
1318 ! Calculate constant for short wave radiation
1320 IF ( ICLOUD == 1 ) THEN
1322 IF ( F_QC .OR. F_QI ) THEN
1323 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
1325 CALL wrf_debug (1, 'CALL cldfra1')
1326 CALL cal_cldfra1(CLDFRA,qv,qc,qi,qs, &
1327 F_QV,F_QC,F_QI,F_QS,t,p, &
1328 F_ICE_PHY,F_RAIN_PHY, &
1329 mp_physics,cldfra1_flag, &
1330 ids,ide, jds,jde, kds,kde, &
1331 ims,ime, jms,jme, kms,kme, &
1332 its,ite, jts,jte, kts,kte )
1333 IF ( PRESENT ( CLDFRA_DP ) ) THEN
1334 ! this is for Kain-Fritsch schemes
1335 IF ( icloud_cu .EQ. 2 .OR. aercu_opt .GT. 0 ) THEN
1336 CALL wrf_debug (1, 'use kf cldfra')
1337 DO j = jts,jte
1338 DO k = kts,kte
1339 DO i = its,ite
1340 cldfra_cu(i,k,j)=cldfra_dp(i,k,j)+cldfra_sh(i,k,j) ! Cu cloud fraction
1341 CLDFRA(i,k,j)=(1.-cldfra_cu(i,k,j))*CLDFRA(i,k,j) ! Update resolved cloud fraction for Cu punch-through
1342 CLDFRA(i,k,j)=CLDFRA(i,k,j)+cldfra_cu(i,k,j) ! New total cloud fraction
1343 CLDFRA(i,k,j)=AMIN1(1.0,CLDFRA(i,k,j))
1344 qc(i,k,j) = qc(i,k,j)+qc_cu(i,k,j)*cldfra_cu(i,k,j)
1345 qi(i,k,j) = qi(i,k,j)+qi_cu(i,k,j)*cldfra_cu(i,k,j)
1346 ENDDO
1347 ENDDO
1348 ENDDO
1349 IF (aercu_opt.gt.0.0) THEN
1350 DO j = jts,jte
1351 DO k = kts,kte
1352 DO i = its,ite
1353 IF (qc(i,k,j).eq.0.0.and.qc_cu(i,k,j).gt.0.0) THEN
1354 qc_cu_weight(i,k,j) = 1.0
1355 ELSE IF (qc(i,k,j).gt.0.0.and.qc_cu(i,k,j).eq.0.0) THEN
1356 qc_cu_weight(i,k,j) = 0.0
1357 ELSE IF (qc(i,k,j).eq.0.0.and.qc_cu(i,k,j).eq.0.0) THEN
1358 qc_cu_weight(i,k,j) = 0.0
1359 ELSE
1360 qc_cu_weight(i,k,j) = (qc_cu(i,k,j)*cldfra_cu(i,k,j))/(qc(i,k,j) + qc_cu(i,k,j)*cldfra_cu(i,k,j))
1361 END IF
1362 IF (qi(i,k,j).eq.0.0.and.qi_cu(i,k,j).gt.0.0) THEN
1363 qi_cu_weight(i,k,j) = 1.0
1364 ELSE IF (qi(i,k,j).gt.0.0.and.qi_cu(i,k,j).eq.0.0) THEN
1365 qi_cu_weight(i,k,j) = 0.0
1366 ELSE IF (qi(i,k,j).eq.0.0.and.qi_cu(i,k,j).eq.0.0) THEN
1367 qi_cu_weight(i,k,j) = 0.0
1368 ELSE
1369 qi_cu_weight(i,k,j) =(qi_cu(i,k,j)*cldfra_cu(i,k,j))/(qi(i,k,j) + qi_cu(i,k,j)*cldfra_cu(i,k,j))
1370 END IF
1371 IF (qs(i,k,j).eq.0.0.and.qs_cu(i,k,j).gt.0.0) THEN
1372 qs_cu_weight(i,k,j) = 1.0
1373 ELSE IF (qs(i,k,j).gt.0.0.and.qs_cu(i,k,j).eq.0.0) THEN
1374 qs_cu_weight(i,k,j) = 0.0
1375 ELSE IF (qs(i,k,j).eq.0.0.and.qs_cu(i,k,j).eq.0.0) THEN
1376 qs_cu_weight(i,k,j) = 0.0
1377 ELSE
1378 qs_cu_weight(i,k,j)=(qs_cu(i,k,j)*cldfra_cu(i,k,j))/(qs(i,k,j) + qs_cu(i,k,j)*cldfra_cu(i,k,j))
1379 END IF
1381 ! use re_cloud, re_ice and re_snow to store combined effective radii from MSKF and Morrison microphysics
1382 re_cloud(i,k,j) = EFCS(I,K,J)*qc_cu_weight(I,K,J) &
1383 + EFCG(I,K,J)*(1-qc_cu_weight(I,K,J))
1384 re_cloud(i,k,j) = re_cloud(i,k,j) * 1.E-6
1385 re_ice(i,k,j) = EFIS(I,K,J)*qi_cu_weight(I,K,J) &
1386 + EFIG(I,K,J)*(1-qi_cu_weight(I,K,J))
1387 re_ice(i,k,j) = re_ice(i,k,j) * 1.E-6
1388 re_snow(i,k,j) = EFSS(I,K,J)*qs_cu_weight(I,K,J) &
1389 + EFSG(I,K,J)*(1-qs_cu_weight(I,K,J))
1390 re_snow(i,k,j) = re_snow(i,k,j) * 1.E-6
1391 has_reqc = 1
1392 has_reqi = 1
1393 has_reqs = 1
1394 qs(i,k,j) = qs(i,k,j)+qs_cu(i,k,j)*cldfra_cu(i,k,j)
1396 ENDDO
1397 ENDDO
1398 ENDDO
1399 END IF
1400 ENDIF
1401 ENDIF
1403 IF ( PRESENT ( CLDFRA_BL ) .AND. PRESENT ( QC_BL ) ) THEN
1404 IF ( icloud_bl > 0 ) THEN
1405 CALL wrf_debug (1, 'in rad driver; use BL clouds')
1406 IF (itimestep .NE. 1) THEN
1407 DO j = jts,jte
1408 DO i = its,ite
1409 DO k = kts,kte
1410 CLDFRA(i,k,j)=CLDFRA_BL(i,k,j)
1411 ENDDO
1412 ENDDO
1413 ENDDO
1414 ENDIF
1416 DO j = jts,jte
1417 DO i = its,ite
1418 DO k = kts,kte
1419 IF (qc(i,k,j) < 1.E-6 .AND. CLDFRA_BL(i,k,j) > 0.001) THEN
1420 qc(i,k,j)=qc(i,k,j) + QC_BL(i,k,j)*CLDFRA_BL(i,k,j)
1421 ENDIF
1422 IF (qi(i,k,j) < 1.E-8 .AND. CLDFRA_BL(i,k,j) > 0.001) THEN
1423 qi(i,k,j)=qi(i,k,j) + QI_BL(i,k,j)*CLDFRA_BL(i,k,j)
1424 ENDIF
1425 ENDDO
1426 ENDDO
1427 ENDDO
1428 ENDIF
1429 ENDIF
1431 IF ( PRESENT (cldfra_mp_all) ) THEN
1432 IF (is_CAMMGMP_used) THEN
1433 !BSINGH: cloud fraction from CAMMGMP is being used (Mods by Po-Lun)
1434 CALL wrf_debug (1, 'use cammgmp')
1435 IF (itimestep .NE. 1) THEN
1436 DO j=jts,jte
1437 DO k=kts,kte
1438 DO i=its,ite
1439 CLDFRA(i,k,j) = CLDFRA_MP_ALL(I,K,J) !PMA
1440 if (CLDFRA(i,k,j) .lt. 0.01) CLDFRA(i,k,j) = 0.
1441 ENDDO
1442 ENDDO
1443 ENDDO
1444 ENDIF
1445 ENDIF
1446 ENDIF
1447 ENDIF
1449 ELSE IF ( ICLOUD == 2 ) THEN
1450 IF ( F_QC .OR. F_QI ) THEN
1451 CALL wrf_debug (1, 'CALL cldfra2')
1452 CALL cal_cldfra2(CLDFRA,qc,qi,F_QC,F_QI, &
1453 ids,ide, jds,jde, kds,kde, &
1454 ims,ime, jms,jme, kms,kme, &
1455 its,ite, jts,jte, kts,kte )
1456 ENDIF
1458 !+---+-----------------------------------------------------------------+
1459 !..New cloud fraction scheme added by G. Thompson (2014Oct31)
1460 !+---+-----------------------------------------------------------------+
1462 ELSEIF (ICLOUD == 3) THEN
1463 IF ( F_QC .AND. F_QI ) THEN
1465 CALL wrf_debug (150, 'DEBUG: using gthompsn cloud fraction scheme')
1467 DO j = jts,jte
1468 DO i = its,ite
1470 DO k = kts,kte
1471 p_1d(k) = p(i,k,j)
1472 if (pert_cld3 .and. multi_perturb == 1) then
1473 t_1d(k) = (1.0 + perts_th(i, k, j) * pert_cld3_t) * t(i, k, j)
1474 qv_1d(k) = max (0.0, (1.0 + perts_qvapor(i, k, j) * pert_cld3_qv) * qv(i, k, j))
1475 else
1476 t_1d(k) = t(i,k,j)
1477 qv_1d(k) = qv(i,k,j)
1478 end if
1479 qc_1d(k) = qc(i,k,j)
1480 qi_1d(k) = qi(i,k,j)
1481 qs_1d(k) = qs(i,k,j)
1482 Dz_1d(k) = dz8w(i,k,j)
1483 cf_1d(k) = cldfra(i,k,j)
1484 ENDDO
1486 WRITE (wrf_err_message,*) 'DEBUG: calling cal_cldfra3 at (i,j): ', i,j, kms,kme,kts,kte
1487 CALL wrf_debug (150, wrf_err_message)
1489 CALL cal_cldfra3(cf_1d, qv_1d, qc_1d, qi_1d, qs_1d, Dz_1d, &
1490 & p_1d, t_1d, XLAND(i,j), gridkm, &
1491 & .false., 1.5, kts, kte, .false.)
1493 DO k = kts,kte
1494 qc(i,k,j) = qc_1d(k)
1495 qi(i,k,j) = qi_1d(k)
1496 qs(i,k,j) = qs_1d(k)
1497 cldfra(i,k,j) = cf_1d(k)
1498 ENDDO
1500 ENDDO
1501 ENDDO
1503 ELSE
1504 CALL wrf_error_fatal('Can not use icloud = 3 option, missing QC or QI field.')
1505 ENDIF
1507 END IF
1509 !Modify CLDFRA and QC for kfcupscheme cumulus scheme
1510 if(present(cldfra_cup)) then
1511 !BSINGH - overwrite cldfra with the cloud fraction computed in module_cu_kfcup.F
1512 !Force cloud fraction if cumulus triggered.
1513 if( cu_physics == KFCUPSCHEME ) then
1514 do j = jts,jte
1515 do k = kts,kte
1516 do i = its,ite
1518 !Find whether to overwrite cldfra or not (ONLY if ICLOUD == 1)
1519 compute_cldfra_cup = .true.
1520 if (icloud == 1 ) then
1521 compute_cldfra_cup = .false. !-- LK Berg, 4/26/17
1522 if(cldfra1_flag(i,k,j) == 1 .and. shall(i,j) .gt. 1) then
1523 CLDFRA(i,k,j)=0.
1524 elseif(cldfra1_flag(i,k,j) == 1 .and. shall(i,j) .le. 1) then
1525 CLDFRA(i,k,j)=0.
1526 compute_cldfra_cup = .true. ! No resolved clouds, but check of shallow clouds. -- LK Berg, 4/26/17
1527 elseif(cldfra1_flag(i,k,j) == 2 .and. shall(i,j) .gt. 1) then
1528 CLDFRA(i,k,j)=1.
1529 elseif(cldfra1_flag(i,k,j) == 3) then
1530 compute_cldfra_cup = .true.
1531 endif
1532 endif
1535 if(compute_cldfra_cup) then
1536 if( (int(shall(i,j)) .le.1) .and. k >= int(cubot(i,j)) .and. k <= int(cutop(i,j)) ) then ! LD Mar232012
1537 CLDFRA(i,k,j) = cldfra_cup(i,k,j)
1538 else if( shall(i,j) .gt. 1) then !!LD
1539 cldfra_cup(i,k,j) = 0.0
1540 end if
1541 endif
1542 if( shall(i,j) <= 1 .and. k >= cubot(i,j) .and. k <= cutop(i,j) ) then ! 1=Shallow Cu -- Changed to use for both deep and shallow -- LK Berg 4/26/17
1543 ! Begin: wig, 4-Feb-2008
1544 !
1545 ! Override the cloud condensate values if shallow convection triggered.
1546 ! For shallow convection, use a representative condensate value based on
1547 ! observations from CHAPS (Oklahoma area) and Florida (Blyth et al. 2005)
1548 ! or the predicted value if it is greater.
1550 cldfra_cup_mod = cldfra_cup(i,k,j) * 1.0e-3 !modified cloud fraction, assume QCLOUD is 1 g/kg -- LK Berg 4/26/17
1551 qc(i,k,j) = max(cldfra_cup_mod, qc(i,k,j) )!DE+LB 2012-Feb
1553 ! Override the cloud fraction values calculated above if shallow
1554 ! convection triggered. For shallow convection, use a representative
1555 ! cloud fraction. In this case, the median value for shallow Cu cases
1556 ! at the ARM SGP site, 36%, Berg et al. 2008, J. Clim.
1557 if( shallowcu_forced_ra )cldfra(i,k,j) = max(0.36, cldfra(i,k,j)) ! Median shallow Cu frac at ARM SGP
1558 endif
1559 ENDDO
1560 ENDDO
1561 ENDDO
1562 end if
1563 endif
1565 #if (EM_CORE==1)
1566 IF( shcu_physics == 5 ) THEN
1567 DO j=jts,jte
1568 DO k=kts,kte
1569 DO i=its,ite
1570 cldfra(I,K,J) = max(cldfra_sh(I,K,J), cldfra(I,K,J))
1571 qc(I,K,J)=cw_rad(I,K,J)+qc(I,K,J)
1572 ENDDO
1573 ENDDO
1574 ENDDO
1575 ENDIF
1576 #endif
1578 IF ( cu_physics == BMJSCHEME .AND. bmj_rad_feedback .EQV. .TRUE. ) THEN
1579 ! hydrometeors from microphysics scheme have saved in q*_save
1580 ! Modify cloud fraction and temporarily hydrometeors (PCC scheme)
1581 DO j=jts,jte
1582 DO k=kts,kte
1583 DO i=its,ite
1584 qc(i,k,j) = qc(i,k,j) + QCCONV(i,k,j)
1585 qi(i,k,j) = qi(i,k,j) + QICONV(i,k,j)
1586 ITRMX=MIN(cldfra(i,k,j),ccldfra(i,k,j))
1587 ITRMN=MAX(0.,cldfra(i,k,j)+ccldfra(i,k,j)-1.)
1588 cldfra(i,k,j)=cldfra(i,k,j)+ccldfra(i,k,j)-0.5*(ITRMX+ITRMN)
1589 ENDDO
1590 ENDDO
1591 ENDDO
1592 ENDIF
1595 END DO
1597 ENDIF Solar_step
1599 Radiation_step: IF ( run_param ) then
1600 ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
1601 STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA
1602 IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset &
1603 .or. STEPABS .eq. 1 ) THEN
1604 doabsems = .true.
1605 ELSE
1606 doabsems = .false.
1607 ENDIF
1608 IF (PRESENT(adapt_step_flag)) THEN
1609 IF ((adapt_step_flag)) THEN
1610 IF ( (itimestep .EQ. 1) .OR. (cam_abs_freq_s .EQ. 0) .OR. &
1611 ( CURR_SECS + dt >= ( INT( CURR_SECS / ( cam_abs_freq_s ) + 1 ) * cam_abs_freq_s) ) ) THEN
1612 doabsems = .true.
1613 ELSE
1614 doabsems = .false.
1615 ENDIF
1616 ENDIF
1617 ENDIF
1619 gfdl_lw = .false.
1620 gfdl_sw = .false.
1621 flg_lw = .false.
1622 flg_sw = .false.
1624 ! Allocate aerosol arrays used by aer_opt = 2 option
1625 IF ( PRESENT( AOD5502D ) ) THEN
1626 ! jararias, 2013/11
1627 IF ( aer_opt .EQ. 2 ) THEN
1628 swrad_aerosol_select: select case(sw_physics)
1630 case(GODDARDSWSCHEME)
1631 allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:11))
1632 allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:11))
1633 allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:11))
1635 case(RRTMG_SWSCHEME &
1636 #if( BUILD_RRTMG_FAST == 1)
1637 ,RRTMG_SWSCHEME_FAST &
1638 #endif
1639 #if( BUILD_RRTMK == 1)
1640 ,RRTMK_SWSCHEME &
1641 #endif
1642 )
1643 allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1644 allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1645 allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1647 end select swrad_aerosol_select
1648 ENDIF
1649 ENDIF
1651 ! Allocate aerosol arrays used by aer_opt = 3 option, and explicit AOD from QNWFA+QNIFA(+QNBCA) (Trude)
1652 IF (PRESENT(f_qnwfa) .AND. PRESENT(f_qnifa) .AND. PRESENT(taod5503d) .AND. PRESENT(taod5502d)) THEN
1653 IF (F_QNWFA .AND. aer_opt.eq.3 .AND. &
1654 (sw_physics.eq.RRTMG_SWSCHEME &
1655 #if( BUILD_RRTMG_FAST == 1)
1656 .OR. sw_physics.eq.RRTMG_SWSCHEME_FAST &
1657 #endif
1658 #if( BUILD_RRTMK == 1)
1659 .OR. sw_physics.eq.RRTMK_SWSCHEME &
1660 #endif
1661 )) THEN
1662 CALL wrf_debug (150, 'DEBUG-GT: computing 3D AOD from QNWFA+QNIFA(+QNBCA)')
1664 allocate(tauaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1665 allocate(ssaaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1666 allocate(asyaer_sw(ims:ime, kms:kme, jms:jme, 1:14))
1668 !$OMP PARALLEL DO &
1669 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1670 DO ij = 1 , num_tiles
1671 its = i_start(ij)
1672 ite = i_end(ij)
1673 jts = j_start(ij)
1674 jte = j_end(ij)
1676 do j=jts,jte
1677 do i=its,ite
1678 taod5502d(i,j) = 0.0
1679 end do
1680 end do
1682 call gt_aod (p, DZ8W, t, qv, qnwfa, qnifa, qnbca, taod5503d, &
1683 wif_input_opt, ims,ime, jms,jme, kms,kme,its,ite, jts,jte, kts,kte)
1685 do j=jts,jte
1686 do i=its,ite
1687 do k=kts,kte
1688 taod5502d(i,j) = taod5502d(i,j) + taod5503d(i,k,j)
1689 end do
1690 end do
1691 end do
1692 ENDDO
1693 !$OMP END PARALLEL DO
1694 ENDIF
1695 ENDIF
1697 !$OMP PARALLEL DO &
1698 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1700 DO ij = 1 , num_tiles
1701 its = i_start(ij)
1702 ite = i_end(ij)
1703 jts = j_start(ij)
1704 jte = j_end(ij)
1706 ! initialize data
1708 if ((itimestep.eq.1).and.(swint_opt.eq.1)) then
1709 do j=jts,jte
1710 do i=its,ite
1711 Bx(i,j)=0.
1712 bb(i,j)=0.
1713 Gx(i,j)=0.
1714 gg(i,j)=0.
1715 end do
1716 end do
1717 end if
1719 DO j=jts,jte
1720 DO i=its,ite
1721 GSW(I,J)=0.
1722 GLW(I,J)=0.
1723 SWDOWN(I,J)=0.
1724 swddir(i,j)=0. ! jararias, 2013/08/10
1725 swddni(i,j)=0. ! jararias, 2013/08/10
1726 swddif(i,j)=0. ! jararias, 2013/08/10
1727 GLAT(I,J)=XLAT(I,J)*DEGRAD
1728 GLON(I,J)=XLONG(I,J)*DEGRAD
1729 ENDDO
1730 ENDDO
1732 DO j=jts,jte
1733 DO k=kts,kte+1
1734 DO i=its,ite
1735 RTHRATEN(I,K,J)=0.
1736 RTHRATENLW(I,K,J)=0.
1737 RTHRATENLWC(I,K,J)=0.
1738 RTHRATENSW(I,K,J)=0.
1739 RTHRATENSWC(I,K,J)=0.
1740 CEMISS(I,K,J)=0.0
1741 ENDDO
1742 ENDDO
1743 ENDDO
1745 IF ( PRESENT( SWUPFLX ) ) THEN
1746 DO j=jts,jte
1747 DO k=kts,kte+2
1748 DO i=its,ite
1749 SWUPFLX(I,K,J) = 0.0
1750 SWDNFLX(I,K,J) = 0.0
1751 SWUPFLXC(I,K,J) = 0.0
1752 SWDNFLXC(I,K,J) = 0.0
1753 LWUPFLX(I,K,J) = 0.0
1754 LWDNFLX(I,K,J) = 0.0
1755 LWUPFLXC(I,K,J) = 0.0
1756 LWDNFLXC(I,K,J) = 0.0
1757 ENDDO
1758 ENDDO
1759 ENDDO
1760 ENDIF
1762 ! these are half level parameters.....so, it should start from kts to kte - 1
1763 !NUWRF JJS 20101021 vvvvv
1764 #if ( WRF_CHEM == 1)
1765 ! Pack gocart aerosol species
1766 ! All aerosol species in chem are in "ug/kg-dryair"
1767 ! and conerted to (g/m**3)
1769 aero(:,:,:,:) = 0.
1770 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
1771 chem_opt == 302 .or. chem_opt == 303) .and. &
1772 (gsfcrad_gocart_coupling == 1) ) then
1773 do j = jts, jte
1774 do k = kts, kte !corrected memory order
1775 do i = its, ite
1776 ! aero(i,k,j, 1) = max(0.0, chem(i,k,j,p_sulf)*1.0e-6*rho(i,k,j)) ! 1 = SO4
1777 aero(i,k,j, 1) = max(0.0, chem(i,k,j,p_sulf)*1.0e-6*p(i,k,j)*96.0/(8.314*t(i,k,j))) ! 1 = SO4
1778 aero(i,k,j, 2) = max(0.0, (chem(i,k,j,p_bc1)+chem(i,k,j,p_bc2))*1.0e-6*rho(i,k,j)) ! 2 = BC1+BC2
1779 aero(i,k,j, 3) = max(0.0, chem(i,k,j,p_oc1)*1.0e-6*rho(i,k,j)*1.4e0) ! 3 = OC1
1780 aero(i,k,j, 4) = max(0.0, chem(i,k,j,p_oc2)*1.0e-6*rho(i,k,j)*1.4e0) ! 4 = OC2
1781 aero(i,k,j, 5) = max(0.0, chem(i,k,j,p_seas_1)*1.0e-6*rho(i,k,j)) ! 5 = SS1
1782 aero(i,k,j, 6) = max(0.0, (chem(i,k,j,p_seas_2)+chem(i,k,j,p_seas_3)+ &
1783 chem(i,k,j,p_seas_4))*1.0e-6*rho(i,k,j)) ! 6 = SS2+SS3+SS4
1784 aero(i,k,j, 7) = max(0.0, chem(i,k,j,p_dust_1)*1.0e-6*rho(i,k,j)*frac(1)) ! 7 = DU1 dust mode 1
1785 aero(i,k,j, 8) = max(0.0, chem(i,k,j,p_dust_1)*1.0e-6*rho(i,k,j)*frac(2)) ! 8 = DU1 dust mode 2
1786 aero(i,k,j, 9) = max(0.0, chem(i,k,j,p_dust_1)*1.0e-6*rho(i,k,j)*frac(3)) ! 9 = DU1 dust mode 3
1787 aero(i,k,j,10) = max(0.0, chem(i,k,j,p_dust_1)*1.0e-6*rho(i,k,j)*frac(4)) ! 10 = DU1 dust mode 4
1788 aero(i,k,j,11) = max(0.0, chem(i,k,j,p_dust_2)*1.0e-6*rho(i,k,j)) ! 11 = DU2 dust mode 5
1789 aero(i,k,j,12) = max(0.0, chem(i,k,j,p_dust_3)*1.0e-6*rho(i,k,j)) ! 11 = DU3 dust mode 6
1790 aero(i,k,j,13) = max(0.0, chem(i,k,j,p_dust_4)*1.0e-6*rho(i,k,j)) ! 11 = DU4 dust mode 7
1791 aero(i,k,j,14) = max(0.0, chem(i,k,j,p_dust_5)*1.0e-6*rho(i,k,j)) ! 11 = DU5 dust mode 8
1792 enddo ! i
1793 enddo ! k
1794 enddo ! j
1795 end if ! if (gsfcrad_gocart_coupling == 1)
1796 #endif
1797 !NUWRF JJS 20101021 ^^^^^
1799 ! Interpolating climatological ozone and aerosol to model time and levels
1800 ! Adapted from camrad code
1801 #if (EM_CORE==1)
1802 IF ( o3input .EQ. 2 .AND. id .EQ. 1 ) THEN
1803 #else
1804 IF ( o3input .EQ. 2 ) THEN
1805 #endif
1806 ! ! Find the current month (adapted from Cavallo)
1807 ! CALL cam_time_interp( ozmixm, pin, levsiz, date_str, &
1808 ! ids , ide , jds , jde , kds , kde , &
1809 ! ims , ime , jms , jme , kms , kme , &
1810 ! its , ite , jts , jte , kts , kte )
1811 ! this is the CAM version
1812 ! interpolate to model time-step
1813 call ozn_time_int(julday,julian,ozmixm,ozmixt,levsiz,n_ozmixm, &
1814 ids , ide , jds , jde , kds , kde , &
1815 ims , ime , jms , jme , kms , kme , &
1816 its , ite , jts , jte , kts , kte )
1818 ! interpolate to model model levels
1819 call ozn_p_int(p ,pin, levsiz, ozmixt, o3rad, &
1820 ids , ide , jds , jde , kds , kde , &
1821 ims , ime , jms , jme , kms , kme , &
1822 its , ite , jts , jte , kts , kte )
1823 ENDIF
1825 IF ( PRESENT( AEROD ) ) THEN
1826 IF ( aer_opt .EQ. 1 .AND. id .EQ. 1 ) THEN
1827 call aer_time_int(julday,julian,aerodm,aerodt,alevsiz,n_ozmixm-1,no_src_types, &
1828 ids , ide , jds , jde , kds , kde , &
1829 ims , ime , jms , jme , kms , kme , &
1830 its , ite , jts , jte , kts , kte )
1832 call aer_p_int(p ,pina, alevsiz, aerodt, aerod, no_src_types, p8w, AODTOT, &
1833 ids , ide , jds , jde , kds , kde , &
1834 ims , ime , jms , jme , kms , kme , &
1835 its , ite , jts , jte , kts , kte )
1836 ENDIF
1837 ENDIF
1839 lwrad_select: SELECT CASE(lw_physics)
1841 CASE (RRTMSCHEME)
1842 CALL wrf_debug (100, 'CALL rrtm')
1844 IF ( PRESENT(p_top) ) THEN
1845 p_top_dummy = p_top
1846 END IF
1848 CALL RRTMLWRAD( &
1849 P_TOP=p_top_dummy &
1850 ,RTHRATEN=RTHRATEN,RTHRATENC=RTHRATENLWC,GLW=GLW &
1851 ,OLR=RLWTOA,EMISS=EMISS &
1852 ,QV3D=QV &
1853 ,QC3D=QC &
1854 ,QR3D=QR &
1855 ,QI3D=QI &
1856 ,QS3D=QS &
1857 ,QG3D=QG &
1858 ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t &
1859 ,T8W=t8w,RHO3D=rho,CLDFRA3D=CLDFRA,R=R_d,G=G &
1860 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
1861 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
1862 ,ICLOUD=icloud,WARM_RAIN=warm_rain &
1863 !ccc Added for time-varying trace gases.
1864 ,YR=YR,JULIAN=JULIAN,GHG_INPUT=GHG_INPUT &
1865 !ccc
1866 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1867 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1868 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1869 )
1872 ! NUWRF Version by Toshihisa Matsui and Jainn Shi 20090623
1873 case (goddardlwscheme)
1875 CALL wrf_debug(100, 'CALL NUWRF goddard longwave radiation scheme 2017')
1877 IF ( mp_physics .NE. nuwrf4icescheme ) THEN
1878 IF ( has_reqc .EQ. 1 ) THEN
1879 DO j=jts,jte
1880 DO k=kts,kte
1881 DO i=its,ite
1882 re_cloud_gsfc(i,k,j) = re_cloud(i,k,j) * 1.E+6
1883 re_ice_gsfc(i,k,j) = re_ice(i,k,j) * 1.E+6
1884 re_snow_gsfc(i,k,j) = re_snow(i,k,j) * 1.E+6
1885 re_rain_gsfc(i,k,j) = 0.
1886 re_graupel_gsfc(i,k,j) = 0.
1887 re_hail_gsfc(i,k,j) = 0.
1888 ENDDO
1889 ENDDO
1890 ENDDO
1891 ELSE
1892 WRITE ( wrf_err_message , * ) 'Must choose a microphysics that provides effective radii.'
1893 CALL wrf_debug (0, wrf_err_message)
1894 END IF
1895 END IF
1897 CALL goddardrad(sw_or_lw='lw',dx=dx &
1898 ,rthraten=rthraten,gsf=glw,xlat=xlat,xlong=xlong &
1899 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
1900 ,dz8w=dz8w,rho_phy=rho,emiss=emiss,tsk=tsk &
1901 ,cldfra3d=cldfra &
1902 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
1903 ,julday=julday,xtime=xtime &
1904 ,declin=declin,solcon=solcon &
1905 ,radfrq=radt,degrad=degrad &
1906 ,warm_rain=warm_rain &
1907 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
1908 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
1909 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
1910 ,qv=qv,qc=qc,qi=qi,qr=qr,qs=qs,qg=qg,qh=qh &
1911 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr &
1912 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg ,f_qh=f_qh &
1913 ,rec3d=re_cloud_gsfc,rei3d=re_ice_gsfc &
1914 ,rer3d=re_rain_gsfc,res3d=re_snow_gsfc & !optional
1915 ,reg3d=re_graupel_gsfc,reh3d=re_hail_gsfc & !optional
1916 ,erbe_out=erbe_out &
1917 ,itimestep=itimestep, dt_in = dt &
1918 ,obscur=obscur_loc &
1919 #if (WRF_CHEM == 1)
1920 ,AERO=aero &
1921 ,CHEM_OPT=chem_opt &
1922 ,GSFCRAD_GOCART_COUPLING=gsfcrad_gocart_coupling &
1923 #endif
1924 )
1927 CASE (GFDLLWSCHEME)
1929 CALL wrf_debug (100, 'CALL gfdllw')
1931 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
1932 PRESENT(F_QI) .AND. (PRESENT(qi) .OR. PRESENT(qs)) .AND. &
1933 PRESENT(qv) .AND. PRESENT(qc) ) THEN
1934 IF ( F_QV .AND. F_QC .AND. (F_QI .OR. F_QS)) THEN
1935 gfdl_lw = .true.
1936 CALL ETARA( &
1937 DT=dt,XLAND=xland &
1938 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
1939 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
1940 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
1941 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
1942 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
1943 ,HBOTR=hbotr, HTOPR=htopr &
1944 ,ALBEDO=albedo,CUPPT=cuppt &
1945 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
1946 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
1947 ,XTIME=xtime,JULIAN=julian &
1948 ,JULYR=julyr,JULDAY=julday &
1949 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
1950 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
1951 ,ACFRST=acfrst,NCFRST=ncfrst &
1952 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
1953 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
1954 ,THRATEN=rthraten,THRATENLW=rthratenlw &
1955 ,THRATENSW=rthratensw &
1956 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1957 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1958 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1959 )
1960 ELSE
1961 CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
1962 ENDIF
1963 ELSE
1964 CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
1965 ENDIF
1967 CASE (CAMLWSCHEME)
1969 CALL wrf_debug(100, 'CALL camrad lw')
1971 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
1972 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
1973 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
1974 .AND. PRESENT(AEROSOLC_2).AND. PRESENT(ALSWVISDIR) ) THEN
1975 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
1976 RTHRATENLWC=RTHRATENLWC,RTHRATENSWC=RTHRATENSWC, &
1977 dolw=.true.,dosw=.false., &
1978 SWUPT=SWUPT,SWUPTC=SWUPTC, &
1979 SWDNT=SWDNT,SWDNTC=SWDNTC, &
1980 LWUPT=LWUPT,LWUPTC=LWUPTC, &
1981 LWDNT=LWDNT,LWDNTC=LWDNTC, &
1982 SWUPB=SWUPB,SWUPBC=SWUPBC, &
1983 SWDNB=SWDNB,SWDNBC=SWDNBC, &
1984 LWUPB=LWUPB,LWUPBC=LWUPBC, &
1985 LWDNB=LWDNB,LWDNBC=LWDNBC, &
1986 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
1987 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
1988 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
1989 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
1990 ,QV3D=qv &
1991 ,QC3D=qc &
1992 ,QR3D=qr &
1993 ,QI3D=qi &
1994 ,QS3D=qs &
1995 ,QG3D=qg &
1996 ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010)
1997 ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010)
1998 ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010)
1999 ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010)
2000 ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds
2001 ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni & !amontornes-bcodina (2014-04-20)
2002 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
2003 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
2004 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
2005 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
2006 dz8w=dz8w, &
2007 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2008 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
2009 num_months=n_ozmixm, &
2010 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
2011 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
2012 paerlev=paerlev, naer_c=n_aerosolc, &
2013 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
2014 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
2015 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
2016 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
2017 ,doabsems=doabsems, ghg_input=ghg_input &
2018 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
2019 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
2020 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
2021 ,coszen=coszen )
2022 ELSE
2023 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
2024 ENDIF
2026 CASE (RRTMG_LWSCHEME)
2028 CALL wrf_debug (100, 'CALL rrtmg_lw')
2030 !Need to reset NLAYERS if vertical nesting is used.
2031 !This code follows that for case RRTMSCHEME within
2032 !subroutine RRTMLWRAD.
2033 IF ( PRESENT(p_top) ) THEN
2034 p_top_dummy = p_top
2035 END IF
2036 IF ( p_top_dummy .GT. 0 ) THEN !
2037 !NLAYERS is recalculated
2038 !every time the radiation scheme is called. This is
2039 !necessary if e_vert parent .NE. e_vert nest since
2040 !NLAYERS could then be different for each domain.
2041 CALL RRTMG_LWINIT( &
2042 p_top, .FALSE. , &
2043 ids, ide, jds, jde, kds, kde, &
2044 ims, ime, jms, jme, kms, kme, &
2045 its, ite, jts, jte, kts, kte )
2046 ENDIF
2048 CALL RRTMG_LWRAD( &
2049 RTHRATENLW=RTHRATEN, &
2050 RTHRATENLWC=RTHRATENLWC, &
2051 LWUPT=LWUPT,LWUPTC=LWUPTC,LWUPTCLN=LWUPTCLN, &
2052 LWDNT=LWDNT,LWDNTC=LWDNTC,LWDNTCLN=LWDNTCLN, &
2053 LWUPB=LWUPB,LWUPBC=LWUPBC,LWUPBCLN=LWUPBCLN, &
2054 LWDNB=LWDNB,LWDNBC=LWDNBC,LWDNBCLN=LWDNBCLN, &
2055 GLW=GLW,OLR=RLWTOA,LWCF=LWCF, &
2056 EMISS=EMISS, &
2057 P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, &
2058 T8W=t8w,RHO3D=rho,R=R_d,G=G, &
2059 ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,&
2060 cldovrlp=cldovrlp,idcor=idcor,XLAT=XLAT, &
2061 #if (EM_CORE == 1)
2062 LRADIUS=lradius, IRADIUS=iradius, &
2063 #endif
2064 IS_CAMMGMP_USED=is_cammgmp_used, &
2066 !ckay
2067 ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,&
2068 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
2069 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2070 QV3D=QV,QC3D=QC,QR3D=QR, &
2071 QI3D=QI,QS3D=QS,QG3D=QG, &
2072 O3INPUT=O3INPUT,O33D=O3RAD, &
2073 F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, &
2074 F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, &
2075 RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson
2076 has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson
2077 #if ( WRF_CHEM == 1 )
2078 TAUAERLW1=tauaerlw1,TAUAERLW2=tauaerlw2, & ! jcb
2079 TAUAERLW3=tauaerlw3,TAUAERLW4=tauaerlw4, & ! jcb
2080 TAUAERLW5=tauaerlw5,TAUAERLW6=tauaerlw6, & ! jcb
2081 TAUAERLW7=tauaerlw7,TAUAERLW8=tauaerlw8, & ! jcb
2082 TAUAERLW9=tauaerlw9,TAUAERLW10=tauaerlw10, & ! jcb
2083 TAUAERLW11=tauaerlw11,TAUAERLW12=tauaerlw12, & ! jcb
2084 TAUAERLW13=tauaerlw13,TAUAERLW14=tauaerlw14, & ! jcb
2085 TAUAERLW15=tauaerlw15,TAUAERLW16=tauaerlw16, & ! jcb
2086 aer_ra_feedback=aer_ra_feedback, &
2087 !jdfcz progn=progn,prescribe=prescribe, &
2088 progn=progn, &
2089 #endif
2090 calc_clean_atm_diag=calc_clean_atm_diag, &
2091 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
2092 !ccc Added for time-varying trace gases.
2093 YR=YR,JULIAN=JULIAN,GHG_INPUT=GHG_INPUT, &
2094 !ccc
2095 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
2096 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
2097 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
2098 LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, &
2099 LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC, &
2100 mp_physics=mp_physics )
2103 #if( BUILD_RRTMK == 1)
2104 CASE (RRTMK_LWSCHEME)
2106 IF ( PRESENT(F_QNC) .AND. PRESENT(QNC_CURR) ) THEN
2107 call rad_rrtmg_driver( &
2108 RTHRATEN, RTHRATENSW, &
2109 RTHRATENLWC, RTHRATENSWC, &
2110 lwupflx, lwupflxc, lwdnflx, lwdnflxc, &
2111 swupflx, swupflxc, swdnflx, swdnflxc, &
2112 lwupt, lwuptc, lwdnt, lwdntc, &
2113 lwupb, lwupbc, lwdnb, lwdnbc, &
2114 glw, olr, lwcf, &
2115 swupt, swuptc, swdnt, swdntc, &
2116 swupb, swupbc, swdnb, swdnbc, &
2117 gsw, swcf, &
2118 coszen, solcon, albedo, emiss, &
2119 t,t8w, tsk, rho, p, p8w, cldfra, &
2120 r_d, g, qnc_curr, xland, &
2121 f_qnc, f_qv, f_qc, f_qr, f_qi, &
2122 f_qs, f_qg, &
2123 qv, qc, qr, qi, qs, qg, &
2124 o3input, o3rad, &
2125 aer_opt, aerod, no_src_types, &
2126 ids,ide, jds,jde, kds,kde, &
2127 ims,ime, jms,jme, kms,kme, &
2128 its,ite, jts,jte, kts,kte)
2129 ELSE
2130 call wrf_error_fatal('Can not call RRTMG-K. Missing QNC field.')
2131 ENDIF
2132 #endif
2135 #if( BUILD_RRTMG_FAST == 1)
2136 CASE (RRTMG_LWSCHEME_FAST)
2137 CALL wrf_debug (100, 'CALL rrtmg_lw')
2139 CALL RRTMG_LWRAD_FAST( &
2140 RTHRATENLW=RTHRATEN, &
2141 RTHRATENLWC=RTHRATENLWC, &
2142 LWUPT=LWUPT,LWUPTC=LWUPTC, &
2143 LWDNT=LWDNT,LWDNTC=LWDNTC, &
2144 LWUPB=LWUPB,LWUPBC=LWUPBC, &
2145 LWDNB=LWDNB,LWDNBC=LWDNBC, &
2146 GLW=GLW,OLR=RLWTOA,LWCF=LWCF, &
2147 EMISS=EMISS, &
2148 P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, &
2149 T8W=t8w,RHO3D=rho,R=R_d,G=G, &
2150 ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,&
2151 #if (EM_CORE == 1)
2152 LRADIUS=lradius, IRADIUS=iradius, &
2153 #endif
2154 IS_CAMMGMP_USED=is_cammgmp_used, &
2156 !ckay
2157 ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,&
2158 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
2159 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2160 QV3D=QV,QC3D=QC,QR3D=QR, &
2161 QI3D=QI,QS3D=QS,QG3D=QG, &
2162 O3INPUT=O3INPUT,O33D=O3RAD, &
2163 F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, &
2164 F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, &
2165 RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson
2166 has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson
2167 #if ( WRF_CHEM == 1 )
2168 TAUAERLW1=tauaerlw1,TAUAERLW2=tauaerlw2, & ! jcb
2169 TAUAERLW3=tauaerlw3,TAUAERLW4=tauaerlw4, & ! jcb
2170 TAUAERLW5=tauaerlw5,TAUAERLW6=tauaerlw6, & ! jcb
2171 TAUAERLW7=tauaerlw7,TAUAERLW8=tauaerlw8, & ! jcb
2172 TAUAERLW9=tauaerlw9,TAUAERLW10=tauaerlw10, & ! jcb
2173 TAUAERLW11=tauaerlw11,TAUAERLW12=tauaerlw12, & ! jcb
2174 TAUAERLW13=tauaerlw13,TAUAERLW14=tauaerlw14, & ! jcb
2175 TAUAERLW15=tauaerlw15,TAUAERLW16=tauaerlw16, & ! jcb
2176 aer_ra_feedback=aer_ra_feedback, &
2177 !jdfcz progn=progn,prescribe=prescribe, &
2178 progn=progn, &
2179 #endif
2180 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
2181 !ccc Added for time-varying trace gases.
2182 YR=YR,JULIAN=JULIAN,GHG_INPUT=GHG_INPUT, &
2183 !ccc
2184 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
2185 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
2186 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
2187 LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, &
2188 LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC &
2189 )
2191 #endif
2193 CASE (HELDSUAREZ)
2194 CALL wrf_debug (100, 'CALL heldsuarez')
2196 CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, &
2197 t8w, rho, R_d,G,CP, dt, xlat, degrad, &
2198 ids,ide, jds,jde, kds,kde, &
2199 ims,ime, jms,jme, kms,kme, &
2200 its,ite, jts,jte, kts,kte )
2202 ! -- add by Jin Kong 10/2011
2203 CASE (FLGLWSCHEME)
2204 CALL wrf_debug (100, 'CALL Fu-Liou-Gu')
2205 flg_lw = .true.
2206 !test Jin Kong 11/01/2011 for ozone
2207 ozflg = 0.
2208 !test over
2209 CALL RAD_FLG( &
2210 peven=p8w,podd=p,t8w=t8w,degrees=t &
2211 ,pi3d=pi,o3=ozflg &
2212 ,G=G,CP=CP &
2213 ,albedo=ALBEDO,tskin=tsk &
2214 ,h2o=QV,cld_iccld=QI,cld_wlcld=QC &
2215 ,cld_prwc=QR,cld_pgwc=QG,cld_snow=QS &
2216 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
2217 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
2218 ,warm_rain=warm_rain &
2219 ,cloudstrf=CLDFRA &
2220 ,emiss=EMISS &
2221 ,air_den=rho &
2222 ,dz3d=dz8w &
2223 ,SOLCON=SOLCON &
2224 ,declin=DECLIN &
2225 ,xtime=xtime, xlong=xlong, xlat=xlat &
2226 ,JULDAY=JULDAY &
2227 ,gmt=gmt, radt=radt, degrad=degrad &
2228 ,dtcolumn=dt &
2229 ,ids=ids,ide=ide,jds=jds,jde=jde &
2230 ,kds=kds,kde=kde &
2231 ,ims=ims,idim=ime,jms=jms,jdim=jme &
2232 ,kms=kms,kmax=kme &
2233 ,its=its,ite=ite,jts=jts,jte=jte &
2234 ,kts=kts,kte=kte &
2235 ,uswtop=RSWTOA,ulwtop=RLWTOA &
2236 ,NETSWBOT=GSW,DLWBOT=GLW &
2237 ,DSWBOT=SWDOWN,deltat=RTHRATEN &
2238 ,dtshort=RTHRATENSW,dtlongwv=RTHRATENLW &
2239 ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni &
2240 )
2242 CALL wrf_debug(100, 'a4 Fu_Liou-Gu')
2243 ! -- end
2245 CASE DEFAULT
2247 WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
2248 CALL wrf_error_fatal ( wrf_err_message )
2250 END SELECT lwrad_select
2252 IF (lw_physics .gt. 0 .and. .not.gfdl_lw .and. .not.flg_lw) THEN
2253 DO j=jts,jte
2254 DO k=kts,kte
2255 DO i=its,ite
2256 RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
2257 ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM
2258 IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
2259 ENDDO
2260 ENDDO
2261 ENDDO
2262 ENDIF
2264 !NUWRF JJS 20090623 vvvvv
2265 IF (lw_physics .eq. goddardlwscheme) THEN
2266 IF ( PRESENT (tlwdn) ) THEN
2267 DO j=jts,jte
2268 DO i=its,ite
2269 tlwdn(i,j) = erbe_out(i,j,1) ! TOA LW downwelling flux [W/m2]
2270 tlwup(i,j) = erbe_out(i,j,2) ! TOA LW upwelling flux [W/m2]
2271 slwdn(i,j) = erbe_out(i,j,3) ! surface LW downwelling flux [W/m2]
2272 slwup(i,j) = erbe_out(i,j,4) ! surface LW upwelling flux [W/m2]
2273 olr(i,j) = -erbe_out(i,j,2)
2274 ENDDO
2275 ENDDO
2276 ENDIF
2277 ENDIF
2278 !NUWRF JJS 20090623 ^^^^^
2280 IF ( PRESENT( AOD5502D ) ) THEN
2281 ! jararias, 2013/11
2282 IF ( aer_opt .EQ. 2 ) THEN
2283 swrad_aerosol_select2: select case(sw_physics)
2285 case(RRTMG_SWSCHEME &
2286 #if( BUILD_RRTMG_FAST == 1)
2287 ,RRTMG_SWSCHEME_FAST &
2288 #endif
2289 #if( BUILD_RRTMK == 1)
2290 ,RRTMK_SWSCHEME &
2291 #endif
2292 )
2293 call wrf_debug(100, 'call calc_aerosol_rrtmg_sw')
2294 call calc_aerosol_rrtmg_sw(ht,dz8w,p,t,qv,aer_type,aer_aod550_opt,aer_angexp_opt, &
2295 aer_ssa_opt,aer_asy_opt,aer_aod550_val,aer_angexp_val, &
2296 aer_ssa_val,aer_asy_val,aod5502d,angexp2d,aerssa2d, &
2297 aerasy2d,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte, &
2298 tauaer_sw,ssaaer_sw,asyaer_sw )
2299 do j=jts,jte
2300 do i=its,ite
2301 do k=kts,kte
2302 aod5503d(i,k,j)=tauaer_sw(i,k,j,10) ! band at 550 nm
2303 end do
2304 end do
2305 end do
2307 case default
2308 end select swrad_aerosol_select2
2309 ENDIF
2310 ENDIF
2312 !..Different treatment for aer_opt=3 using QNWFA+QNIFA aerosol species (Trude)
2314 IF (PRESENT(f_qnwfa) .AND. PRESENT(f_qnifa)) THEN
2315 IF (F_QNWFA .AND. aer_opt.eq.3 .AND. &
2316 (sw_physics.eq.RRTMG_SWSCHEME &
2317 #if( BUILD_RRTMG_FAST == 1)
2318 .OR. sw_physics.eq.RRTMG_SWSCHEME_FAST &
2319 #endif
2320 #if( BUILD_RRTMK == 1)
2321 .OR. sw_physics.eq.RRTMK_SWSCHEME &
2322 #endif
2323 )) THEN
2324 call wrf_debug(100, 'call calc_aerosol_rrtmg_sw with 3D AOD values')
2325 call calc_aerosol_rrtmg_sw(ht,dz8w,p,t,qv,taer_type,taer_aod550_opt,taer_angexp_opt, &
2326 taer_ssa_opt,taer_asy_opt,aer_aod550_val,aer_angexp_val, &
2327 aer_ssa_val,aer_asy_val,taod5502d,angexp2d,aerssa2d, &
2328 aerasy2d,ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte, &
2329 tauaer_sw,ssaaer_sw,asyaer_sw, taod5503d)
2331 do j=jts,jte
2332 do i=its,ite
2333 aod5502d(i,j)= taod5502d(i, j)
2334 end do
2335 end do
2336 ENDIF
2337 ENDIF
2339 swrad_select: SELECT CASE(sw_physics)
2341 CASE (SWRADSCHEME)
2342 CALL wrf_debug(100, 'CALL swrad')
2343 CALL SWRAD( &
2344 DT=dt,RTHRATEN=rthraten,GSW=gsw &
2345 ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
2346 #if ( WRF_CHEM == 1 )
2347 ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
2348 ,PM2_5_DRY_EC=pm2_5_dry_ec &
2349 #endif
2350 ,RHO_PHY=rho,T3D=t &
2351 ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
2352 ,R=r_d,CP=cp,G=g,JULDAY=julday &
2353 ,GHG_INPUT=ghg_input &
2354 ,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
2355 ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
2356 ,warm_rain=warm_rain &
2357 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
2358 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
2359 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
2360 ,QV3D=qv &
2361 ,QC3D=qc &
2362 ,QR3D=qr &
2363 ,QI3D=qi &
2364 ,QS3D=qs &
2365 ,QG3D=qg &
2366 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
2367 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
2368 ,coszen=coszen,julian=julian &
2369 ,obscur=obscur_loc )
2371 CASE (GSFCSWSCHEME)
2372 CALL wrf_debug(100, 'CALL gsfcswrad')
2373 CALL GSFCSWRAD( &
2374 RTHRATEN=rthraten,GSW=gsw & ! PAJ: xlat and xlong removed.
2375 ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
2376 ,DZ8W=dz8w,RHO_PHY=rho &
2377 ,CLDFRA3D=cldfra,RSWTOA=rswtoa &
2378 ,CP=cp,G=g & ! PAJ: GMT removed.
2379 ,JULDAY=julday & ! PAJ: XTIME removed.
2380 ,SOLCON=solcon & ! PAJ: declin removed
2381 ! ,RADFRQ=radt,DEGRAD=degrad & ! PAJ: degrad and radfrq removed
2382 ,TAUCLDI=taucldi,TAUCLDC=taucldc &
2383 ,WARM_RAIN=warm_rain &
2385 #if ( WRF_CHEM == 1 )
2386 ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
2387 ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
2388 ,GAER300=gaer300,GAER400=gaer400 & ! jcb
2389 ,GAER600=gaer600,GAER999=gaer999 & ! jcb
2390 ,WAER300=waer300,WAER400=waer400 & ! jcb
2391 ,WAER600=waer600,WAER999=waer999 & ! jcb
2392 ,aer_ra_feedback=aer_ra_feedback &
2393 #endif
2394 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
2395 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
2396 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
2397 ,QV3D=qv &
2398 ,QC3D=qc &
2399 ,QR3D=qr &
2400 ,QI3D=qi &
2401 ,QS3D=qs &
2402 ,QG3D=qg &
2403 ,QNDROP3D=qndrop &
2404 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
2405 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
2406 ,F_QNDROP=f_qndrop &
2407 ,COSZEN=coszen &
2408 ,OBSCUR=obscur_loc &
2409 )
2411 case (goddardswscheme)
2413 CALL wrf_debug(100, 'CALL NUWRF goddard shortwave radiation scheme 2017')
2415 IF ( mp_physics .NE. nuwrf4icescheme ) THEN
2416 IF ( has_reqc .EQ. 1 ) THEN
2417 DO j=jts,jte
2418 DO k=kts,kte
2419 DO i=its,ite
2420 re_cloud_gsfc(i,k,j) = re_cloud(i,k,j) * 1.E+6
2421 re_ice_gsfc(i,k,j) = re_ice(i,k,j) * 1.E+6
2422 re_snow_gsfc(i,k,j) = re_snow(i,k,j) * 1.E+6
2423 re_rain_gsfc(i,k,j) = 0.
2424 re_graupel_gsfc(i,k,j) = 0.
2425 re_hail_gsfc(i,k,j) = 0.
2426 ENDDO
2427 ENDDO
2428 ENDDO
2429 ELSE
2430 WRITE ( wrf_err_message , * ) 'Must choose a microphysics that provides effective radii.'
2431 CALL wrf_debug (0, wrf_err_message)
2432 END IF
2433 END IF
2435 CALL goddardrad(sw_or_lw='sw',dx=dx &
2436 ,rthraten=rthraten,gsf=gsw,xlat=xlat,xlong=xlong &
2437 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
2438 ,dz8w=dz8w,rho_phy=rho,emiss=emiss,tsk=tsk &
2439 ,cldfra3d=cldfra &
2440 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
2441 ,julday=julday,xtime=xtime &
2442 ,declin=declin,solcon=solcon &
2443 ,radfrq=radt,degrad=degrad &
2444 ,warm_rain=warm_rain &
2445 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
2446 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
2447 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
2448 ,qv=qv,qc=qc,qr=qr,qi=qi,qs=qs,qg=qg,qh=qh & !optional
2449 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr & !optional
2450 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg ,f_qh=f_qh & !optional
2451 ,rec3d=re_cloud_gsfc,rei3d=re_ice_gsfc & !optional
2452 ,rer3d=re_rain_gsfc,res3d=re_snow_gsfc & !optional
2453 ,reg3d=re_graupel_gsfc,reh3d=re_hail_gsfc & !optional
2454 ,erbe_out=erbe_out &
2455 ,cod2d_out=cod2d_out,ctop2d_out=ctop2d_out & !optional
2456 ,sflxd=sflxd & !optional
2457 ,swddir=swddir,swddni=swddni,swddif=swddif & !optional
2458 ,coszen=coszen & !optional
2459 ,obscur=obscur_loc &
2460 ,itimestep=itimestep, dt_in = dt &
2461 #if (WRF_CHEM == 1)
2462 ,aod2d_out=aod2d_out, atop2d_out=atop2d_out & ! optional
2463 ,AERO=aero &
2464 ,CHEM_OPT=chem_opt &
2465 ,GSFCRAD_GOCART_COUPLING=gsfcrad_gocart_coupling &
2466 #endif
2467 )
2469 CASE (CAMSWSCHEME)
2470 CALL wrf_debug(100, 'CALL camrad sw')
2471 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
2472 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
2473 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
2474 .AND. PRESENT(AEROSOLC_2) .AND. PRESENT(ALSWVISDIR)) THEN
2475 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
2476 RTHRATENLWC=RTHRATENLWC,RTHRATENSWC=RTHRATENSWC, &
2477 dolw=.false.,dosw=.true., &
2478 SWUPT=SWUPT,SWUPTC=SWUPTC, &
2479 SWDNT=SWDNT,SWDNTC=SWDNTC, &
2480 LWUPT=LWUPT,LWUPTC=LWUPTC, &
2481 LWDNT=LWDNT,LWDNTC=LWDNTC, &
2482 SWUPB=SWUPB,SWUPBC=SWUPBC, &
2483 SWDNB=SWDNB,SWDNBC=SWDNBC, &
2484 LWUPB=LWUPB,LWUPBC=LWUPBC, &
2485 LWDNB=LWDNB,LWDNBC=LWDNBC, &
2486 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
2487 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
2488 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
2489 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
2490 ,QV3D=qv &
2491 ,QC3D=qc &
2492 ,QR3D=qr &
2493 ,QI3D=qi &
2494 ,QS3D=qs &
2495 ,QG3D=qg &
2496 ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010)
2497 ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010)
2498 ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010)
2499 ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010)
2500 ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds
2501 ,SWDDIR=swddir,SWDDIF=swddif,SWDDNI=swddni & !amontornes-bcodina (2014-04-20)
2502 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
2503 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
2504 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
2505 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
2506 dz8w=dz8w, &
2507 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2508 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
2509 num_months=n_ozmixm, &
2510 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
2511 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
2512 paerlev=paerlev, naer_c=n_aerosolc, &
2513 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
2514 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
2515 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
2516 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
2517 ,doabsems=doabsems,ghg_input=ghg_input &
2518 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
2519 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
2520 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
2521 ,coszen=coszen )
2522 ELSE
2523 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
2524 ENDIF
2525 DO j=jts,jte
2526 DO k=kts,kte
2527 DO i=its,ite
2528 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
2529 ENDDO
2530 ENDDO
2531 ENDDO
2533 CASE (RRTMG_SWSCHEME)
2535 CALL wrf_debug(100, 'CALL rrtmg_sw')
2537 if ( direct_sw_feedback .and. feedback_is_ready ) then
2538 proceed_cmaq_sw = .true.
2540 if (present(mass_ws_i)) then
2542 mass_ws_i(:,kte+1,:) = mass_ws_i(:,kte,:)
2543 mass_ws_j(:,kte+1,:) = mass_ws_j(:,kte,:)
2544 mass_ws_k(:,kte+1,:) = mass_ws_k(:,kte,:)
2546 mass_in_i(:,kte+1,:) = mass_in_i(:,kte,:)
2547 mass_in_j(:,kte+1,:) = mass_in_j(:,kte,:)
2548 mass_in_k(:,kte+1,:) = mass_in_k(:,kte,:)
2550 mass_ec_i(:,kte+1,:) = mass_ec_i(:,kte,:)
2551 mass_ec_j(:,kte+1,:) = mass_ec_j(:,kte,:)
2552 mass_ec_k(:,kte+1,:) = mass_ec_k(:,kte,:)
2554 mass_ss_i(:,kte+1,:) = mass_ss_i(:,kte,:)
2555 mass_ss_j(:,kte+1,:) = mass_ss_j(:,kte,:)
2556 mass_ss_k(:,kte+1,:) = mass_ss_k(:,kte,:)
2558 mass_h2o_i(:,kte+1,:) = mass_h2o_i(:,kte,:)
2559 mass_h2o_j(:,kte+1,:) = mass_h2o_j(:,kte,:)
2560 mass_h2o_k(:,kte+1,:) = mass_h2o_k(:,kte,:)
2562 dgn_i(:,kte+1,:) = dgn_i(:,kte,:)
2563 dgn_j(:,kte+1,:) = dgn_j(:,kte,:)
2564 dgn_k(:,kte+1,:) = dgn_k(:,kte,:)
2566 sig_i(:,kte+1,:) = sig_i(:,kte,:)
2567 sig_j(:,kte+1,:) = sig_j(:,kte,:)
2568 sig_k(:,kte+1,:) = sig_k(:,kte,:)
2569 end if
2570 else
2571 proceed_cmaq_sw = .false.
2572 end if
2574 CALL RRTMG_SWRAD( &
2575 RTHRATENSW=RTHRATENSW, &
2576 RTHRATENSWC=RTHRATENSWC, &
2577 SWUPT=SWUPT,SWUPTC=SWUPTC,SWUPTCLN=SWUPTCLN, &
2578 SWDNT=SWDNT,SWDNTC=SWDNTC,SWDNTCLN=SWDNTCLN, &
2579 SWUPB=SWUPB,SWUPBC=SWUPBC,SWUPBCLN=SWUPBCLN, &
2580 SWDNB=SWDNB,SWDNBC=SWDNBC,SWDNBCLN=SWDNBCLN, &
2581 SWCF=SWCF,GSW=GSW, &
2582 XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, &
2583 RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, &
2584 COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, &
2585 ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, &
2586 p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, &
2587 dz8w=dz8w,CLDFRA3D=CLDFRA,GHG_INPUT=GHG_INPUT, &
2588 #if (EM_CORE == 1)
2589 LRADIUS=lradius, IRADIUS=iradius, &
2590 #endif
2591 IS_CAMMGMP_USED=is_cammgmp_used, &
2592 R=R_D,G=G, &
2593 !ckay
2594 ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,&
2595 ICLOUD=icloud,WARM_RAIN=warm_rain, &
2596 cldovrlp=cldovrlp,idcor=idcor, &
2597 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
2598 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2599 QV3D=qv,QC3D=qc,QR3D=qr, &
2600 QI3D=qi,QS3D=qs,QG3D=qg, &
2601 O3INPUT=O3INPUT,O33D=O3RAD, &
2602 AER_OPT=AER_OPT,aerod=aerod,no_src=no_src_types, &
2603 ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif, & !Zhenxin ssib alb comp (06/2010)
2604 ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif, & !Zhenxin ssib alb comp (06/2010)
2605 SWVISDIR=swvisdir ,SWVISDIF=swvisdif, & !Zhenxin ssib swr comp (06/2010)
2606 SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif, & !Zhenxin ssib swr comp (06/2010)
2607 SF_SURFACE_PHYSICS=sf_surface_physics, & !Zhenxin ssib sw_phy (06/2010)
2608 F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, &
2609 F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, &
2610 RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson
2611 has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson
2612 #if ( WRF_CHEM == 1 )
2613 TAUAER300=tauaer300,TAUAER400=tauaer400, & ! jcb
2614 TAUAER600=tauaer600,TAUAER999=tauaer999, & ! jcb
2615 GAER300=gaer300,GAER400=gaer400, & ! jcb
2616 GAER600=gaer600,GAER999=gaer999, & ! jcb
2617 WAER300=waer300,WAER400=waer400, & ! jcb
2618 WAER600=waer600,WAER999=waer999, & ! jcb
2619 aer_ra_feedback=aer_ra_feedback, &
2620 !jdfcz progn=progn,prescribe=prescribe, &
2621 progn=progn, &
2622 #endif
2623 calc_clean_atm_diag=calc_clean_atm_diag, &
2624 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
2625 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
2626 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
2627 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
2628 SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, &
2629 SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC, &
2630 tauaer3d_sw=tauaer_sw, & ! jararias 2013/11
2631 ssaaer3d_sw=ssaaer_sw, & ! jararias 2013/11
2632 asyaer3d_sw=asyaer_sw, & ! jararias 2013/11
2633 swddir=swddir,swddni=swddni,swddif=swddif, & ! jararias 2013/08/10
2634 swdownc=swdownc, swddnic=swddnic, swddirc=swddirc, & ! PAJ
2635 obscur=obscur_loc, &
2636 xcoszen=coszen,yr=yr,julian=julian,mp_physics=mp_physics, & ! jararias 2013/08/14
2637 proceed_cmaq_sw=proceed_cmaq_sw, & ! WRF-CMAQ coupled model
2638 nmode=3, & ! WRF-CMAQ coupled model
2639 mass_ws_i=mass_ws_i, & ! WRF-CMAQ coupled model
2640 mass_ws_j=mass_ws_j, & ! WRF-CMAQ coupled model
2641 mass_ws_k=mass_ws_k, & ! WRF-CMAQ coupled model
2642 mass_in_i=mass_in_i, & ! WRF-CMAQ coupled model
2643 mass_in_j=mass_in_j, & ! WRF-CMAQ coupled model
2644 mass_in_k=mass_in_k, & ! WRF-CMAQ coupled model
2645 mass_ec_i=mass_ec_i, & ! WRF-CMAQ coupled model
2646 mass_ec_j=mass_ec_j, & ! WRF-CMAQ coupled model
2647 mass_ec_k=mass_ec_k, & ! WRF-CMAQ coupled model
2648 mass_ss_i=mass_ss_i, & ! WRF-CMAQ coupled model
2649 mass_ss_j=mass_ss_j, & ! WRF-CMAQ coupled model
2650 mass_ss_k=mass_ss_k, & ! WRF-CMAQ coupled model
2651 mass_h2o_i=mass_h2o_i, & ! WRF-CMAQ coupled model
2652 mass_h2o_j=mass_h2o_j, & ! WRF-CMAQ coupled model
2653 mass_h2o_k=mass_h2o_k, & ! WRF-CMAQ coupled model
2654 dgn_i=dgn_i, & ! WRF-CMAQ coupled model
2655 dgn_j=dgn_j, & ! WRF-CMAQ coupled model
2656 dgn_k=dgn_k, & ! WRF-CMAQ coupled model
2657 sig_i=sig_i, & ! WRF-CMAQ coupled model
2658 sig_j=sig_j, & ! WRF-CMAQ coupled model
2659 sig_k=sig_k, & ! WRF-CMAQ coupled model
2660 gtauxar_01=sw_gtauxar_01, & ! WRF-CMAQ coupled model
2661 gtauxar_02=sw_gtauxar_02, & ! WRF-CMAQ coupled model
2662 gtauxar_03=sw_gtauxar_03, & ! WRF-CMAQ coupled model
2663 gtauxar_04=sw_gtauxar_04, & ! WRF-CMAQ coupled model
2664 gtauxar_05=sw_gtauxar_05, & ! WRF-CMAQ coupled model
2665 asy_fac_01=sw_asy_fac_01, & ! WRF-CMAQ coupled model
2666 asy_fac_02=sw_asy_fac_02, & ! WRF-CMAQ coupled model
2667 asy_fac_03=sw_asy_fac_03, & ! WRF-CMAQ coupled model
2668 asy_fac_04=sw_asy_fac_04, & ! WRF-CMAQ coupled model
2669 asy_fac_05=sw_asy_fac_05, & ! WRF-CMAQ coupled model
2670 ssa_01=sw_ssa_01, & ! WRF-CMAQ coupled model
2671 ssa_02=sw_ssa_02, & ! WRF-CMAQ coupled model
2672 ssa_03=sw_ssa_03, & ! WRF-CMAQ coupled model
2673 ssa_04=sw_ssa_04, & ! WRF-CMAQ coupled model
2674 ssa_05=sw_ssa_05, & ! WRF-CMAQ coupled model
2675 sw_zbbcddir=sw_zbbcddir, & ! WRF-CMAQ coupled model
2676 sw_dirdflux=sw_dirdflux, & ! WRF-CMAQ coupled model
2677 sw_difdflux=sw_difdflux & ! WRF-CMAQ coupled model
2678 )
2680 ! = WRF-CMAQ twoway coupled model
2681 if (proceed_cmaq_sw) then
2682 DO j=jts,jte
2683 DO i=its,ite
2684 sw_ttauxar_01(i,j) = sum(sw_gtauxar_01(i,:,j))
2685 sw_ttauxar_02(i,j) = sum(sw_gtauxar_02(i,:,j))
2686 sw_ttauxar_03(i,j) = sum(sw_gtauxar_03(i,:,j))
2687 sw_ttauxar_04(i,j) = sum(sw_gtauxar_04(i,:,j))
2688 sw_ttauxar_05(i,j) = sum(sw_gtauxar_05(i,:,j))
2689 END DO
2690 END DO
2691 end if
2693 DO j=jts,jte
2694 DO k=kts,kte
2695 DO i=its,ite
2696 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
2697 ENDDO
2698 ENDDO
2699 ENDDO
2701 #if( BUILD_RRTMK == 1)
2702 CASE (RRTMK_SWSCHEME)
2704 DO j=jts,jte
2705 DO k=kts,kte
2706 DO i=its,ite
2707 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
2708 ENDDO
2709 ENDDO
2710 ENDDO
2711 #endif
2713 #if( BUILD_RRTMG_FAST == 1)
2714 CASE (RRTMG_SWSCHEME_FAST)
2715 CALL wrf_debug(100, 'CALL rrtmg_sw_fast')
2716 CALL RRTMG_SWRAD_FAST( &
2717 RTHRATENSW=RTHRATENSW, &
2718 RTHRATENSWC=RTHRATENSWC, &
2719 SWUPT=SWUPT,SWUPTC=SWUPTC, &
2720 SWDNT=SWDNT,SWDNTC=SWDNTC, &
2721 SWUPB=SWUPB,SWUPBC=SWUPBC, &
2722 SWDNB=SWDNB,SWDNBC=SWDNBC, &
2723 SWCF=SWCF,GSW=GSW, &
2724 XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, &
2725 RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, &
2726 COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, &
2727 ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, &
2728 p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, &
2729 dz8w=dz8w,CLDFRA3D=CLDFRA, &
2730 GHG_INPUT=GHG_INPUT, &
2731 #if (EM_CORE == 1)
2732 LRADIUS=lradius, IRADIUS=iradius, &
2733 #endif
2734 IS_CAMMGMP_USED=is_cammgmp_used, &
2735 R=R_D,G=G, &
2736 !ckay
2737 ! CLDFRA_KF3D=cldfra_KF,QC_KF3D=qc_KF,QI_KF3D=qi_KF,&
2738 ICLOUD=icloud,WARM_RAIN=warm_rain, &
2739 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
2740 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
2741 QV3D=qv,QC3D=qc,QR3D=qr, &
2742 QI3D=qi,QS3D=qs,QG3D=qg, &
2743 O3INPUT=O3INPUT,O33D=O3RAD, &
2744 AER_OPT=AER_OPT,aerod=aerod,no_src=no_src_types, &
2745 ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif, & !Zhenxin ssib alb comp (06/2010)
2746 ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif, & !Zhenxin ssib alb comp (06/2010)
2747 SWVISDIR=swvisdir ,SWVISDIF=swvisdif, & !Zhenxin ssib swr comp (06/2010)
2748 SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif, & !Zhenxin ssib swr comp (06/2010)
2749 SF_SURFACE_PHYSICS=sf_surface_physics, & !Zhenxin ssib sw_phy (06/2010)
2750 F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, &
2751 F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, &
2752 RE_CLOUD=re_cloud,RE_ICE=re_ice,RE_SNOW=re_snow, & ! G. Thompson
2753 has_reqc=has_reqc,has_reqi=has_reqi,has_reqs=has_reqs, & ! G. Thompson
2754 #if ( WRF_CHEM == 1 )
2755 TAUAER300=tauaer300,TAUAER400=tauaer400, & ! jcb
2756 TAUAER600=tauaer600,TAUAER999=tauaer999, & ! jcb
2757 GAER300=gaer300,GAER400=gaer400, & ! jcb
2758 GAER600=gaer600,GAER999=gaer999, & ! jcb
2759 WAER300=waer300,WAER400=waer400, & ! jcb
2760 WAER600=waer600,WAER999=waer999, & ! jcb
2761 aer_ra_feedback=aer_ra_feedback, &
2762 !jdfcz progn=progn,prescribe=prescribe, &
2763 progn=progn, &
2764 #endif
2765 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
2766 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
2767 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
2768 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
2769 SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, &
2770 SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC, &
2771 tauaer3d_sw=tauaer_sw, & ! jararias 2013/11
2772 ssaaer3d_sw=ssaaer_sw, & ! jararias 2013/11
2773 asyaer3d_sw=asyaer_sw, & ! jararias 2013/11
2774 swddir=swddir,swddni=swddni,swddif=swddif, & ! jararias 2013/08/10
2775 swdownc=swdownc, swddnic=swddnic, swddirc=swddirc, & ! PAJ
2776 xcoszen=coszen,yr=yr,julian=julian ) ! jararias 2013/08/14
2778 DO j=jts,jte
2779 DO k=kts,kte
2780 DO i=its,ite
2781 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
2782 ENDDO
2783 ENDDO
2784 ENDDO
2785 #endif
2787 CASE (GFDLSWSCHEME)
2789 CALL wrf_debug (100, 'CALL gfdlsw')
2791 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
2792 PRESENT(F_QI) .AND. (PRESENT(qi) .OR. PRESENT(qs)) .AND. &
2793 PRESENT(qv) .AND. PRESENT(qc) ) THEN
2794 IF ( F_QV .AND. F_QC .AND. (F_QI .OR. F_QS)) THEN
2795 gfdl_sw = .true.
2796 CALL ETARA( &
2797 DT=dt,XLAND=xland &
2798 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
2799 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
2800 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
2801 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
2802 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
2803 ,HBOTR=hbotr, HTOPR=htopr &
2804 ,ALBEDO=albedo,CUPPT=cuppt &
2805 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
2806 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
2807 ,XTIME=xtime,JULIAN=julian &
2808 ,JULYR=julyr,JULDAY=julday &
2809 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
2810 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
2811 ,ACFRST=acfrst,NCFRST=ncfrst &
2812 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
2813 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
2814 ,THRATEN=rthraten,THRATENLW=rthratenlw &
2815 ,THRATENSW=rthratensw &
2816 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
2817 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
2818 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
2819 )
2820 ELSE
2821 CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
2822 ENDIF
2823 ELSE
2824 CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
2825 ENDIF
2827 CASE (0)
2829 ! Here in case we don't want to call a sw radiation scheme
2830 ! For example, the Held-Suarez idealized test case
2831 IF (lw_physics /= HELDSUAREZ) THEN
2832 WRITE( wrf_err_message , * ) &
2833 'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')'
2834 CALL wrf_error_fatal ( wrf_err_message )
2835 END IF
2837 ! -- add by Jin Kong 10/2011
2838 !--- the following FLGSWSCHEME is for testing only
2839 CASE (FLGSWSCHEME)
2840 flg_sw = .true.
2841 !--- No need to do anything since the short and long wave is calculted in one program
2842 ! -- end
2844 CASE DEFAULT
2846 WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
2847 CALL wrf_error_fatal ( wrf_err_message )
2849 END SELECT swrad_select
2851 !NUWRF JJS 20090623 vvvvv
2852 IF (sw_physics .eq. goddardswscheme) THEN
2853 IF ( PRESENT (tswdn) ) THEN
2854 DO j=jts,jte
2855 DO i=its,ite
2856 tswdn(i,j) = erbe_out(i,j,5) ! TOA SW downwelling flux [W/m2]
2857 tswup(i,j) = erbe_out(i,j,6) ! TOA SW upwelling flux [W/m2]
2858 sswdn(i,j) = erbe_out(i,j,7) ! surface SW downwelling flux [W/m2]
2859 sswup(i,j) = erbe_out(i,j,8) ! surface SW upwelling flux [W/m2]
2860 ENDDO
2861 ENDDO
2862 ENDIF
2863 ENDIF
2864 !NUWRF JJS 20090623 ^^^^^
2866 IF (sw_physics .gt. 0 .and. .not.gfdl_sw .and. .not.flg_sw) THEN
2867 DO j=jts,jte
2868 DO k=kts,kte
2869 DO i=its,ite
2870 RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
2871 ENDDO
2872 ENDDO
2873 ENDDO
2875 DO j=jts,jte
2876 DO i=its,ite
2877 SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
2878 ENDDO
2879 ENDDO
2880 ENDIF
2882 ! jararias, 14/08/2013
2883 ! surface direct and diffuse SW fluxes computation. Only for schemes other than RRTMG and Goddard
2884 ! Backup method in case sw scheme in use does not provide surface SW direct and diffuse irradiances
2885 IF ((sw_physics .NE. RRTMG_SWSCHEME) &
2886 #if( BUILD_RRTMG_FAST == 1)
2887 .AND. (sw_physics .NE. RRTMG_SWSCHEME_FAST) &
2888 #endif
2889 .AND. (sw_physics .NE. FLGSWSCHEME) .AND. (sw_physics .NE. CAMSWSCHEME) & ! amontornes-bcodina (2014-04-20)
2890 #if( BUILD_RRTMK == 1)
2891 .AND. (sw_physics .NE. RRTMK_SWSCHEME) &
2892 #endif
2893 .AND. (sw_physics .ne. GODDARDSWSCHEME)) THEN
2894 DO j=jts,jte
2895 DO i=its,ite
2896 IF (coszen(i,j).GT.1e-3) THEN
2897 ioh=solcon*coszen(i,j) ! TOA irradiance
2898 kt=swdown(i,j)/max(ioh,1e-3) ! clearness index
2899 ! Optical air mass: Rigollier et al. (2000) doi:
2900 ! 10.1016/S0038-092X(99)00055-9
2901 airmass=exp(-ht(i,j)/8434.5)/(coszen(i,j)+ &
2902 0.50572*(asin(coszen(i,j))*57.295779513082323+6.07995)**(-1.6364))
2903 ! kt correction for air-mass at large sza: Perez et al. (1990)
2904 ! doi: 10.1016/0038-092X(90)90036-C
2905 kt=kt/(0.1+1.031*exp(-1.4/(0.9+(9.4/max(airmass,1e-3)))))
2906 ! Diffuse fraction: Ruiz-Arias et al. (2010) (Eq 33) doi:
2907 ! 10.1016/j.enconman.2009.11.024
2908 kd=0.952-1.041*exp(-exp(2.300-4.702*kt))
2909 swddif(i,j)=kd*swdown(i,j)
2910 swddir(i,j)=(1.-kd)*swdown(i,j)
2911 swddni(i,j)=swddir(i,j)/max(coszen(i,j),1e-4)
2912 ENDIF
2913 ENDDO
2914 ENDDO
2915 ENDIF
2917 IF ( PRESENT( diffuse_frac ) ) THEN
2918 DO j=jts,jte
2919 DO i=its,ite
2920 if (swdown(i,j).gt.0.001) then
2921 diffuse_frac(i,j) = swddif(i,j)/swdown(i,j)
2922 diffuse_frac(i,j) = min(diffuse_frac(i,j),1.0)
2923 else
2924 diffuse_frac(i,j) = 0.
2925 endif
2926 ENDDO
2927 ENDDO
2928 ENDIF
2930 ! jararias, aug 2013, updated 2013/11
2931 ! parameters update for SW surface fluxes interpolation
2932 IF (swint_opt.EQ.1) THEN
2933 ! interpolation applies on all-sky fluxes (swddir, swdown)
2934 CALL update_swinterp_parameters(ims,ime,jms,jme,its,ite,jts,jte, &
2935 coszen,coszen_loc,swddir,swdown, &
2936 swddir_ref,bb,Bx,swdown_ref,gg,Gx, &
2937 coszen_ref )
2938 ENDIF
2940 IF ( PRESENT( obscur ) ) THEN
2941 IF ( sw_eclipse == eclipsescheme ) THEN
2942 DO j=jts,jte
2943 DO i=its,ite
2944 obscur(i,j) = obscur_loc(i,j)
2945 mask(i,j) = mask_loc(i,j)
2946 ENDDO
2947 ENDDO
2948 elat_track = elat_loc
2949 elon_track = elon_loc
2950 ENDIF
2951 ENDIF
2953 ENDDO
2954 !$OMP END PARALLEL DO
2956 IF ( associated(tauaer_sw) ) deallocate(tauaer_sw)
2957 IF ( associated(ssaaer_sw) ) deallocate(ssaaer_sw)
2958 IF ( associated(asyaer_sw) ) deallocate(asyaer_sw)
2960 ENDIF Radiation_step
2962 ! jararias, aug 2013
2963 ! SW surface fluxes interpolation (meaningful when not in a Radiation_step)
2964 if (swint_opt .eq. 1) then
2965 call wrf_debug(100,'SW surface irradiance interpolation')
2967 !---------------
2968 !$OMP PARALLEL DO &
2969 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
2970 do ij = 1,num_tiles
2971 its = i_start(ij)
2972 ite = i_end(ij)
2973 jts = j_start(ij)
2974 jte = j_end(ij)
2975 call interp_sw_radiation(ims,ime,jms,jme,its,ite,jts,jte, &
2976 coszen_ref,coszen_loc,swddir_ref, &
2977 bb,Bx,swdown_ref,gg,Gx,albedo, &
2978 swdown,swddir,swddni,swddif,gsw )
2979 enddo
2980 !$OMP END PARALLEL DO
2981 end if
2983 ! Coupling with FARMS
2984 if( present(swdown2 ) .and. &
2985 present(swddir2 ) .and. &
2986 present(swddni2 ) .and. &
2987 present(swddif2 ) .and. &
2988 present(swdownc2 ) .and. &
2989 present(swddnic2 ) ) then
2990 if (swint_opt == 2) then
2991 call wrf_debug(100,'SW surface irradiance calculated with FARMS')
2993 select case (aer_opt)
2994 case (0)
2995 !$OMP PARALLEL DO &
2996 !$OMP PRIVATE (ij ,i, j, its, ite, jts, jte)
2997 do ij = 1,num_tiles
2998 its = i_start(ij)
2999 ite = i_end(ij)
3000 jts = j_start(ij)
3001 jte = j_end(ij)
3002 DO j=jts,jte
3003 DO i=its,ite
3004 aod5502d(i, j) = 0.0
3005 angexp2d(i, j) = 0.0
3006 aerssa2d(i, j) = 0.0
3007 aerasy2d(i, j) = 0.0
3008 ENDDO
3009 ENDDO
3010 end do
3011 !$OMP END PARALLEL DO
3013 case (1)
3014 !$OMP PARALLEL DO &
3015 !$OMP PRIVATE (ij ,i, j, its, ite, jts, jte)
3016 do ij = 1,num_tiles
3017 its = i_start(ij)
3018 ite = i_end(ij)
3019 jts = j_start(ij)
3020 jte = j_end(ij)
3021 DO j=jts,jte
3022 DO i=its,ite
3023 aod5502d(i, j) = aodtot(i, j)
3024 angexp2d(i, j) = aer_angexp_val
3025 aerssa2d(i, j) = aer_ssa_val
3026 aerasy2d(i, j) = aer_asy_val
3027 ENDDO
3028 ENDDO
3029 end do
3030 !$OMP END PARALLEL DO
3032 case (2)
3033 continue
3035 case (3)
3036 continue
3038 case default
3039 CALL wrf_error_fatal('Unkonown aer_opt to FARMS. Only options 0, 1, 2, and 3 are supported.')
3041 end select
3043 !---------------
3044 !$OMP PARALLEL DO &
3045 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
3046 do ij = 1,num_tiles
3047 its = i_start(ij)
3048 ite = i_end(ij)
3049 jts = j_start(ij)
3050 jte = j_end(ij)
3052 if (pert_farms .and. multi_perturb == 1) then
3053 allocate (qv_tmp(its:ite, kts:kte, jts:jte))
3054 allocate (qc_tmp(its:ite, kts:kte, jts:jte))
3055 allocate (qs_tmp(its:ite, kts:kte, jts:jte))
3057 call Add_multi_perturb_swrad_perturbations (perts_albedo, perts_aod, perts_angstrom, &
3058 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow, pert_farms_albedo, pert_farms_aod, &
3059 pert_farms_angexp, pert_farms_aerasy, pert_farms_qv, pert_farms_qc, pert_farms_qs, &
3060 albedo, aod5502d, angexp2d, aerasy2d, qv, qc, qs, qv_tmp, qc_tmp, qs_tmp, its, ite, jts, jte, &
3061 ims, ime, jms, jme, kms, kme, kts, kte)
3062 end if
3064 call farms_driver (ims, ime, jms, jme, its, ite, jts, jte, kms, kme, kts, kte, &
3065 p8w, rho, dz8w, albedo, aer_opt, aerssa2d, aerasy2d, aod5502d, angexp2d, &
3066 coszen_loc, qv, qi, qs, qc, re_cloud, re_ice, re_snow, &
3067 julian, swdown2, swddir2, swddni2, swddif2, swdownc2, swddnic2, &
3068 has_reqc, has_reqi, has_reqs, CLDFRA)
3070 if (pert_farms .and. multi_perturb == 1) then
3071 call Remove_multi_perturb_swrad_perturbations (perts_albedo, perts_aod, perts_angstrom, &
3072 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow, pert_farms_albedo, pert_farms_aod, &
3073 pert_farms_angexp, pert_farms_aerasy, pert_farms_qv, pert_farms_qc, pert_farms_qs, &
3074 albedo, aod5502d, angexp2d, aerasy2d, qv, qc, qs, qv_tmp, qc_tmp, qs_tmp, its, ite, jts, jte, &
3075 ims, ime, jms, jme, kms, kme, kts, kte)
3077 deallocate (qv_tmp)
3078 deallocate (qc_tmp)
3079 deallocate (qs_tmp)
3080 end if
3081 enddo
3082 !$OMP END PARALLEL DO
3084 if (couple_farms) then
3085 !$OMP PARALLEL DO &
3086 !$OMP PRIVATE (ij, i, j, its, ite, jts, jte)
3087 do ij = 1,num_tiles
3088 its = i_start(ij)
3089 ite = i_end(ij)
3090 jts = j_start(ij)
3091 jte = j_end(ij)
3092 do j = jts, jte
3093 do i = its, ite
3094 swdown(i, j) = swdown2(i, j)
3095 end do
3096 end do
3097 end do
3098 !$OMP END PARALLEL DO
3099 end if
3100 end if
3101 end if
3103 solar_opt_local = 0
3104 IF ( PRESENT(solar_opt) ) THEN
3105 solar_opt_local = solar_opt
3106 END IF
3107 IF (run_param .or. solar_opt_local == 1 .or. swint_opt == 2) THEN
3108 do ij = 1,num_tiles
3109 its = i_start(ij)
3110 ite = i_end(ij)
3111 jts = j_start(ij)
3112 jte = j_end(ij)
3113 ! Save resolved + unresolved hydrometeor mass mixing ratios for Solar diag
3114 IF ( solar_opt_local == 1 ) THEN
3115 IF ( PRESENT(qc_tot) .AND. PRESENT(qi_tot) ) THEN
3116 IF ( F_QC ) THEN
3117 DO j=jts,jte
3118 DO k=kts,kte
3119 DO i=its,ite
3120 qc_tot(i,k,j) = qc(i,k,j)
3121 ENDDO
3122 ENDDO
3123 ENDDO
3124 ENDIF
3125 IF ( F_QI ) THEN
3126 DO j=jts,jte
3127 DO k=kts,kte
3128 DO i=its,ite
3129 qi_tot(i,k,j) = qi(i,k,j)
3130 ENDDO
3131 ENDDO
3132 ENDDO
3133 ENDIF
3134 END IF
3135 ENDIF
3136 ! Restore qc & qi for any model physics configuration
3137 IF ( F_QC ) THEN
3138 DO j=jts,jte
3139 DO k=kts,kte
3140 DO i=its,ite
3141 qc(i,k,j) = qc_save(i,k,j)
3142 ENDDO
3143 ENDDO
3144 ENDDO
3145 ENDIF
3146 IF ( F_QI ) THEN
3147 DO j=jts,jte
3148 DO k=kts,kte
3149 DO i=its,ite
3150 qi(i,k,j) = qi_save(i,k,j)
3151 ENDDO
3152 ENDDO
3153 ENDDO
3154 ENDIF
3156 IF (aercu_opt.gt.0.0) THEN
3157 IF ( F_QS ) THEN
3158 DO j=jts,jte
3159 DO k=kts,kte
3160 DO i=its,ite
3161 qs(i,k,j) = qs_save(i,k,j)
3162 ENDDO
3163 ENDDO
3164 ENDDO
3165 ENDIF
3166 END IF
3167 end do
3168 END IF
3170 accumulate_lw_select: SELECT CASE(lw_physics)
3172 CASE (CAMLWSCHEME,&
3173 RRTMG_LWSCHEME &
3174 #if( BUILD_RRTMG_FAST == 1)
3175 ,RRTMG_LWSCHEME_FAST &
3176 #endif
3177 #if( BUILD_RRTMK == 1)
3178 ,RRTMK_LWSCHEME &
3179 #endif
3180 )
3181 IF(PRESENT(LWUPTC))THEN
3182 ! EM calls every DT
3183 DTaccum = DT
3184 !$OMP PARALLEL DO &
3185 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
3187 DO ij = 1 , num_tiles
3188 its = i_start(ij)
3189 ite = i_end(ij)
3190 jts = j_start(ij)
3191 jte = j_end(ij)
3193 DO j=jts,jte
3194 DO i=its,ite
3195 ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DTaccum
3196 ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DTaccum
3197 ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DTaccum
3198 ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DTaccum
3199 ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DTaccum
3200 ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DTaccum
3201 ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DTaccum
3202 ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DTaccum
3203 ENDDO
3204 ENDDO
3205 ENDDO
3206 !$OMP END PARALLEL DO
3207 ENDIF
3208 CASE DEFAULT
3209 END SELECT accumulate_lw_select
3211 accumulate_sw_select: SELECT CASE(sw_physics)
3213 CASE (CAMSWSCHEME,&
3214 RRTMG_SWSCHEME &
3215 #if( BUILD_RRTMG_FAST == 1)
3216 ,RRTMG_SWSCHEME_FAST &
3217 #endif
3218 #if( BUILD_RRTMK == 1)
3219 ,RRTMK_SWSCHEME &
3220 #endif
3221 )
3222 IF(PRESENT(SWUPTC))THEN
3223 ! EM calls the driver every DT
3224 DTaccum = DT
3225 !$OMP PARALLEL DO &
3226 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
3228 DO ij = 1 , num_tiles
3229 its = i_start(ij)
3230 ite = i_end(ij)
3231 jts = j_start(ij)
3232 jte = j_end(ij)
3234 DO j=jts,jte
3235 DO i=its,ite
3236 ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DTaccum
3237 ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DTaccum
3238 ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DTaccum
3239 ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DTaccum
3240 ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DTaccum
3241 ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DTaccum
3242 ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DTaccum
3243 ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DTaccum
3244 ENDDO
3245 ENDDO
3246 ENDDO
3247 !$OMP END PARALLEL DO
3248 ENDIF
3250 CASE DEFAULT
3251 END SELECT accumulate_sw_select
3253 ! compute cloud diagnosis (random overlapping)
3254 ! by ZCX
3255 IF ( PRESENT ( CLDFRA ) .AND. PRESENT ( CLDT ) .AND. &
3256 PRESENT ( F_QC ) .AND. PRESENT ( F_QI ) ) THEN
3258 DO ij = 1 , num_tiles
3259 its = i_start(ij)
3260 ite = i_end(ij)
3261 jts = j_start(ij)
3262 jte = j_end(ij)
3264 DO j=jts,jte
3265 DO i=its,ite
3266 cldji=1.0
3267 do k=kte-1,kts,-1
3268 cldji=cldji*(1.0-cldfra(i,k,j))
3269 enddo
3270 cldt(i,j)=1.0-cldji
3271 ! cldlji=1.0
3272 ! do k=kte-1,kts,-1
3273 ! if(znu(k).ge.0.69) then
3274 ! cldlji=cldlji*(1.0-cldfra(i,k,j))
3275 ! endif
3276 ! enddo
3277 ! cldl(i,j)=1.0-cldlji
3278 END DO
3279 END DO
3280 END DO
3281 END IF
3283 END SUBROUTINE radiation_driver
3285 SUBROUTINE pre_radiation_driver ( grid, config_flags &
3286 ,itimestep, ra_call_offset &
3287 ,XLAT, XLONG, GMT, julian, xtime, RADT, STEPRA &
3288 ,ht,dx,dy,dx2d,area2d &
3289 ,sina,cosa,shadowmask,slope_rad ,topo_shading &
3290 ,shadlen,ht_shad,ht_loc &
3291 ,ht_shad_bxs, ht_shad_bxe &
3292 ,ht_shad_bys, ht_shad_bye &
3293 ,nested, min_ptchsz &
3294 ,spec_bdy_width &
3295 ,ids, ide, jds, jde, kds, kde &
3296 ,ims, ime, jms, jme, kms, kme &
3297 ,ips, ipe, jps, jpe, kps, kpe &
3298 ,i_start, i_end &
3299 ,j_start, j_end &
3300 ,kts, kte &
3301 ,num_tiles )
3303 USE module_domain , ONLY : domain
3304 #ifdef DM_PARALLEL
3305 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks,wrf_dm_minval_integer
3306 # if (EM_CORE == 1)
3307 USE module_comm_dm , ONLY : halo_toposhad_sub
3308 # endif
3309 #endif
3310 USE module_bc
3311 USE module_model_constants
3313 IMPLICIT NONE
3315 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3316 ims,ime, jms,jme, kms,kme, &
3317 ips,ipe, jps,jpe, kps,kpe, &
3318 kts,kte, &
3319 num_tiles
3321 TYPE(domain) , INTENT(INOUT) :: grid
3322 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
3324 INTEGER, INTENT(IN ) :: itimestep, ra_call_offset, stepra, &
3325 slope_rad, topo_shading, &
3326 spec_bdy_width
3328 INTEGER, INTENT(INOUT) :: min_ptchsz
3330 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
3331 i_start,i_end,j_start,j_end
3333 REAL, INTENT(IN ) :: GMT, radt, julian, xtime, dx, dy, shadlen
3335 REAL, DIMENSION( ims:ime, jms:jme ), &
3336 INTENT(IN ) :: XLAT, &
3337 XLONG, &
3338 HT, &
3339 SINA, &
3340 COSA
3341 REAL, DIMENSION( ims:ime, jms:jme ), &
3342 INTENT(IN ), OPTIONAL :: DX2D, &
3343 AREA2D
3345 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
3347 REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width ), &
3348 INTENT(IN ) :: ht_shad_bxs, ht_shad_bxe
3349 REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width ), &
3350 INTENT(IN ) :: ht_shad_bys, ht_shad_bye
3352 INTEGER, DIMENSION( ims:ime, jms:jme ), &
3353 INTENT(INOUT) :: shadowmask
3355 LOGICAL, INTENT(IN ) :: nested
3357 !Local
3358 ! For orographic shading
3359 INTEGER :: niter,ni,psx,psy,idum,jdum,i,j,ij
3360 REAL :: DECLIN,SOLCON
3362 ! Determine minimum patch size for slope-dependent radiation
3363 if (itimestep .eq. 1) then
3364 psx = ipe-ips+1
3365 psy = jpe-jps+1
3366 min_ptchsz = min(psx,psy)
3367 idum = 0
3368 jdum = 0
3369 endif
3371 # ifdef DM_PARALLEL
3372 if (itimestep .eq. 1) then
3373 call wrf_dm_minval_integer (psx,idum,jdum)
3374 call wrf_dm_minval_integer (psy,idum,jdum)
3375 min_ptchsz = min(psx,psy)
3376 endif
3377 # endif
3379 ! Topographic shading
3381 if ((topo_shading.eq.1).and.(itimestep .eq. 1 .or. &
3382 mod(itimestep,STEPRA) .eq. 1 + ra_call_offset)) then
3384 !---------------
3385 ! Calculate constants for short wave radiation
3387 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,DEGRAD,DPD)
3389 ! Make a local copy of terrain height field
3390 do j=jms,jme
3391 do i=ims,ime
3392 ht_loc(i,j) = ht(i,j)
3393 enddo
3394 enddo
3395 ! Determine if iterations are necessary for shadows to propagate from one patch to another
3396 if ((ids.eq.ips).and.(ide.eq.ipe).and.(jds.eq.jps).and.(jde.eq.jpe)) then
3397 niter = 1
3398 else
3399 niter = int(shadlen/(dx*min_ptchsz)+3)
3400 endif
3403 IF( nested ) THEN
3405 !$OMP PARALLEL DO &
3406 !$OMP PRIVATE ( ij )
3408 DO ij = 1 , num_tiles
3410 CALL spec_bdyfield(ht_shad, &
3411 ht_shad_bxs, ht_shad_bxe, &
3412 ht_shad_bys, ht_shad_bye, &
3413 'm', config_flags, spec_bdy_width, 2,&
3414 ids,ide, jds,jde, 1 ,1 , & ! domain dims
3415 ims,ime, jms,jme, 1 ,1 , & ! memory dims
3416 ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
3417 i_start(ij), i_end(ij), &
3418 j_start(ij), j_end(ij), &
3419 1 , 1 )
3420 ENDDO
3421 ENDIF
3423 do ni = 1, niter
3425 !$OMP PARALLEL DO &
3426 !$OMP PRIVATE ( ij,i,j )
3427 do ij = 1 , num_tiles
3429 call toposhad_init (ht_shad,ht_loc, &
3430 shadowmask,nested,ni, &
3431 ids,ide, jds,jde, kds,kde, &
3432 ims,ime, jms,jme, kms,kme, &
3433 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
3434 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
3435 min(j_end(ij), jde-1), kts, kte )
3437 enddo
3438 !$OMP END PARALLEL DO
3441 !$OMP PARALLEL DO &
3442 !$OMP PRIVATE ( ij,i,j )
3443 do ij = 1 , num_tiles
3445 call toposhad (xlat,xlong,sina,cosa,xtime,gmt,radt,declin, &
3446 dx,dy,ht_shad,ht_loc,ni, &
3447 shadowmask,shadlen, &
3448 ids,ide, jds,jde, kds,kde, &
3449 ims,ime, jms,jme, kms,kme, &
3450 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
3451 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
3452 min(j_end(ij), jde-1), kts, kte )
3454 enddo
3455 !$OMP END PARALLEL DO
3457 #if defined( DM_PARALLEL ) && (EM_CORE == 1)
3458 # include "HALO_TOPOSHAD.inc"
3459 #endif
3460 enddo
3461 endif
3463 END SUBROUTINE pre_radiation_driver
3465 !---------------------------------------------------------------------
3466 !BOP
3467 ! !IROUTINE: radconst - compute radiation terms
3468 ! !INTERFAC:
3469 SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
3470 DEGRAD,DPD )
3471 !---------------------------------------------------------------------
3472 USE module_wrf_error
3473 IMPLICIT NONE
3474 !---------------------------------------------------------------------
3476 ! !ARGUMENTS:
3477 REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
3478 REAL, INTENT(OUT ) :: DECLIN,SOLCON
3479 REAL :: OBECL,SINOB,SXLONG,ARG, &
3480 DECDEG,DJUL,RJUL,ECCFAC
3481 !
3482 ! !DESCRIPTION:
3483 ! Compute terms used in radiation physics
3484 !EOP
3486 ! for short wave radiation
3488 DECLIN=0.
3489 SOLCON=0.
3491 !-----OBECL : OBLIQUITY = 23.5 DEGREE.
3493 OBECL=23.5*DEGRAD
3494 SINOB=SIN(OBECL)
3496 !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
3498 IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
3499 IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
3500 SXLONG=SXLONG*DEGRAD
3501 ARG=SINOB*SIN(SXLONG)
3502 DECLIN=ASIN(ARG)
3503 DECDEG=DECLIN/DEGRAD
3504 !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
3505 DJUL=JULIAN*360./365.
3506 RJUL=DJUL*DEGRAD
3507 ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* &
3508 COS(2*RJUL)+0.000077*SIN(2*RJUL)
3509 SOLCON=1370.*ECCFAC
3511 END SUBROUTINE radconst
3514 SUBROUTINE calc_coszen(ims,ime,jms,jme,its,ite,jts,jte, &
3515 julian,xtime,gmt, &
3516 declin,degrad,xlon,xlat,coszen,hrang)
3517 ! Added Equation of Time correction : jararias, 2013/08/10
3518 implicit none
3519 integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte
3520 real, intent(in) :: julian,declin,xtime,gmt,degrad
3521 real, dimension(ims:ime,jms:jme), intent(in) :: xlat,xlon
3522 real, dimension(ims:ime,jms:jme), intent(inout) :: coszen,hrang
3524 integer :: i,j
3525 real :: da,eot,xt24,tloctm,xxlat
3527 da=6.2831853071795862*(julian-1)/365.
3528 eot=(0.000075+0.001868*cos(da)-0.032077*sin(da) &
3529 -0.014615*cos(2*da)-0.04089*sin(2*da))*(229.18)
3530 xt24=mod(xtime,1440.)+eot
3531 do j=jts,jte
3532 do i=its,ite
3533 tloctm=gmt+xt24/60.+xlon(i,j)/15.
3534 hrang(i,j)=15.*(tloctm-12.)*degrad
3535 xxlat=xlat(i,j)*degrad
3536 coszen(i,j)=sin(xxlat)*sin(declin) &
3537 +cos(xxlat)*cos(declin) *cos(hrang(i,j))
3538 coszen(i, j) = min (max (coszen(i, j), -1.0), 1.0)
3539 enddo
3540 enddo
3541 END SUBROUTINE calc_coszen
3543 subroutine update_swinterp_parameters(ims,ime,jms,jme,its,ite,jts,jte, &
3544 coszen,coszen_loc,swddir,swdown, &
3545 swddir_ref,bb,Bx, &
3546 swdown_ref,gg,Gx, &
3547 coszen_ref )
3548 ! Author: jararias 2013/11
3549 implicit None
3550 integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte
3551 real, dimension(ims:ime,jms:jme), intent(in) :: coszen,coszen_loc,swddir,swdown
3552 real, dimension(ims:ime,jms:jme), intent(inout) :: swddir_ref,bb,Bx, &
3553 swdown_ref,gg,Gx, &
3554 coszen_ref
3556 integer :: i,j
3557 real :: swddir_0,swdown_0,coszen_0
3558 real, parameter :: coszen_min=1e-4
3560 do j=jts,jte
3561 do i=its,ite
3562 if ((coszen(i,j).gt.coszen_min) .and. (coszen_loc(i,j).gt.coszen_min)) then
3563 ! parameters update for DIR
3564 if (Bx(i,j).le.0) then
3565 swddir_0 =(coszen_loc(i,j)/coszen(i,j))*swddir(i,j) ! linear first guess estimation
3566 coszen_0 =coszen_loc(i,j)
3567 else
3568 swddir_0 =swddir_ref(i,j)
3569 coszen_0 =coszen_ref(i,j)
3570 end if
3571 if ((coszen(i,j)/coszen_0).lt.1.) then
3572 bb(i,j) =log(max(1.,swddir(i,j))/max(1.,swddir_0)) / log(min(1.-1e-4,coszen(i,j)/coszen_0))
3573 elseif ((coszen(i,j)/coszen_0).gt.1) then
3574 bb(i,j) =log(max(1.,swddir(i,j))/max(1.,swddir_0)) / log(max(1.+1e-4,coszen(i,j)/coszen_0))
3575 else
3576 bb(i,j) =0.
3577 end if
3578 bb(i,j) =max(-.5,min(2.5,bb(i,j)))
3579 Bx(i,j) =swddir(i,j)/(coszen(i,j)**bb(i,j))
3581 !write(wrf_err_message,*) 'XXX I=',i,' J=',j,' Bx=',Bx(i,j),' bb=',bb(i,j),' swddir=',swddir(i,j), &
3582 ! ' swddir_0=',swddir_0,' coszen=',coszen(i,j),' coszen_0=',coszen_0
3583 !call wrf_debug(1,wrf_err_message)
3585 ! parameters update for GHI
3586 if (Gx(i,j).le.0) then
3587 swdown_0 =(coszen_loc(i,j)/coszen(i,j))*swdown(i,j) ! linear first guess estimation
3588 coszen_0 =coszen_loc(i,j)
3589 else
3590 swdown_0 =swdown_ref(i,j)
3591 coszen_0 =coszen_ref(i,j)
3592 end if
3593 if ((coszen(i,j)/coszen_0).lt.1.) then
3594 gg(i,j) =log(max(1.,swdown(i,j))/max(1.,swdown_0)) / log(min(1.-1e-4,coszen(i,j)/coszen_0))
3595 elseif ((coszen(i,j)/coszen_0).gt.1) then
3596 gg(i,j) =log(max(1.,swdown(i,j))/max(1.,swdown_0)) / log(max(1.+1e-4,coszen(i,j)/coszen_0))
3597 else
3598 gg(i,j) =0.
3599 end if
3600 gg(i,j) =max(-.5,min(2.5,gg(i,j)))
3601 Gx(i,j) =swdown(i,j)/(coszen(i,j)**gg(i,j))
3602 else
3603 Bx(i,j) =0.
3604 bb(i,j) =0.
3605 Gx(i,j) =0.
3606 gg(i,j) =0.
3607 end if
3609 ! saving last SW run in state variables
3610 coszen_ref(i,j) =coszen(i,j)
3611 swdown_ref(i,j) =swdown(i,j)
3612 swddir_ref(i,j) =swddir(i,j)
3614 !if ((i.eq.20).and.(j.eq.20)) then
3615 ! write(wrf_err_message,'(" RADSTEP : tn=",I4," csz_0=",F9.6," csz=",F9.6," csz_1=",F9.6," Gx=",F14.2," gg=",F9.5, &
3616 ! " Bx=",F14.2," bb=",F9.5)') itimestep,coszen_0,coszen_loc(i,j),coszen(i,j),Gx(i,j),gg(i,j), &
3617 ! Bx(i,j),bb(i,j)
3618 ! call wrf_debug(1,wrf_err_message)
3619 !end if
3621 end do
3622 end do
3624 end subroutine update_swinterp_parameters
3626 subroutine interp_sw_radiation(ims,ime,jms,jme,its,ite,jts,jte, &
3627 coszen_ref,coszen_loc,swddir_ref, &
3628 bb,Bx,swdown_ref,gg,Gx,albedo, &
3629 swdown,swddir,swddni,swddif,gsw )
3630 ! Author: jararias 2013/11
3631 implicit None
3632 integer, intent(in) :: ims,ime,jms,jme,its,ite,jts,jte
3633 real, dimension(ims:ime,jms:jme), intent(in) :: coszen_ref,coszen_loc, &
3634 swddir_ref,Bx,bb, &
3635 swdown_ref,Gx,gg, &
3636 albedo
3638 real, dimension(ims:ime,jms:jme), intent(inout) :: swddir,swdown, &
3639 swddif,swddni, gsw
3641 integer :: i,j
3642 real, parameter :: coszen_min=1e-4
3644 do j=jts,jte
3645 do i=its,ite
3646 ! sza interpolation of surface fluxes
3647 if ((coszen_ref(i,j).gt.coszen_min) .and. (coszen_loc(i,j).gt.coszen_min)) then
3648 if ((bb(i,j).eq.-0.5).or.(bb(i,j).eq.2.5).or.(bb(i,j).eq.0.0)) then
3649 swddir(i,j) =(coszen_loc(i,j)/coszen_ref(i,j))*swddir_ref(i,j)
3650 else
3651 swddir(i,j) =Bx(i,j)*(coszen_loc(i,j)**bb(i,j))
3652 end if
3653 if ((gg(i,j).eq.-0.5).or.(gg(i,j).eq.2.5).or.(gg(i,j).eq.0.0)) then
3654 swdown(i,j) =(coszen_loc(i,j)/coszen_ref(i,j))*swdown_ref(i,j)
3655 else
3656 swdown(i,j) =Gx(i,j)*(coszen_loc(i,j)**gg(i,j))
3657 end if
3658 swddif(i,j) =swdown(i,j)-swddir(i,j)
3659 swddni(i,j) =swddir(i,j)/coszen_loc(i,j)
3660 gsw(i,j) =swdown(i,j)*(1.-albedo(i,j))
3661 else
3662 swddir(i,j) =0.
3663 swdown(i,j) =0.
3664 swddif(i,j) =0.
3665 swddni(i,j) =0.
3666 gsw(i,j) =0.
3667 end if
3668 end do
3669 end do
3670 end subroutine interp_sw_radiation
3672 !---------------------------------------------------------------------
3673 !BOP
3674 ! !IROUTINE: cal_cldfra2 - Compute cloud fraction
3675 ! !INTERFACE:
3676 SUBROUTINE cal_cldfra2(CLDFRA,QC,QI,F_QC,F_QI, &
3677 ids,ide, jds,jde, kds,kde, &
3678 ims,ime, jms,jme, kms,kme, &
3679 its,ite, jts,jte, kts,kte )
3680 USE module_state_description, ONLY : KFCUPSCHEME, KFETASCHEME !CuP, wig 5-Oct-2006 !BSINGH - For WRFCuP scheme
3681 !---------------------------------------------------------------------
3682 IMPLICIT NONE
3683 !---------------------------------------------------------------------
3684 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3685 ims,ime, jms,jme, kms,kme, &
3686 its,ite, jts,jte, kts,kte
3688 !
3689 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
3690 CLDFRA
3692 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
3693 QI, &
3694 QC
3695 LOGICAL,INTENT(IN) :: F_QC,F_QI
3697 REAL thresh
3698 INTEGER:: i,j,k
3699 ! !DESCRIPTION:
3700 ! Compute cloud fraction from input ice and cloud water fields
3701 ! if provided.
3702 !
3703 ! Whether QI or QC is active or not is determined from the indices of
3704 ! the fields into the 4D scalar arrays in WRF. These indices are
3705 ! P_QI and P_QC, respectively, and they are passed in to the routine
3706 ! to enable testing to see if QI and QC represent active fields in
3707 ! the moisture 4D scalar array carried by WRF.
3708 !
3709 ! If a field is active its index will have a value greater than or
3710 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
3711 ! this routine.
3712 !EOP
3713 !---------------------------------------------------------------------
3714 thresh=1.0e-6
3716 IF ( f_qi .AND. f_qc ) THEN
3717 DO j = jts,jte
3718 DO k = kts,kte
3719 DO i = its,ite
3720 IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
3721 CLDFRA(i,k,j)=1.
3722 ELSE
3723 CLDFRA(i,k,j)=0.
3724 ENDIF
3725 ENDDO
3726 ENDDO
3727 ENDDO
3728 ELSE IF ( f_qc ) THEN
3729 DO j = jts,jte
3730 DO k = kts,kte
3731 DO i = its,ite
3732 IF ( QC(i,k,j) .gt. thresh) THEN
3733 CLDFRA(i,k,j)=1.
3734 ELSE
3735 CLDFRA(i,k,j)=0.
3736 ENDIF
3737 ENDDO
3738 ENDDO
3739 ENDDO
3740 ELSE
3741 DO j = jts,jte
3742 DO k = kts,kte
3743 DO i = its,ite
3744 CLDFRA(i,k,j)=0.
3745 ENDDO
3746 ENDDO
3747 ENDDO
3748 ENDIF
3749 END SUBROUTINE cal_cldfra2
3751 !BOP
3752 ! !IROUTINE: cal_cldfra1 - Compute cloud fraction
3753 ! !INTERFACE:
3754 ! cal_cldfra_xr - Compute cloud fraction.
3755 ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
3756 !!
3757 !!--- Cloud fraction parameterization follows Xu and Randall (JAS), 1996
3758 !! (see Hong et al., 1998)
3759 !! (modified by Ferrier, Feb '02)
3760 !
3761 SUBROUTINE cal_cldfra1(CLDFRA, QV, QC, QI, QS, &
3762 F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
3763 F_ICE_PHY,F_RAIN_PHY, &
3764 mp_physics, cldfra1_flag, &
3765 ids,ide, jds,jde, kds,kde, &
3766 ims,ime, jms,jme, kms,kme, &
3767 its,ite, jts,jte, kts,kte )
3768 USE module_state_description, ONLY : KFCUPSCHEME, KFETASCHEME !wig, CuP 4-Fb-2008 !BSINGH - For WRFCuP scheme
3770 USE module_state_description, ONLY : FER_MP_HIRES, FER_MP_HIRES_ADVECT
3771 !---------------------------------------------------------------------
3772 IMPLICIT NONE
3773 !---------------------------------------------------------------------
3774 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3775 ims,ime, jms,jme, kms,kme, &
3776 its,ite, jts,jte, kts,kte
3778 !
3779 INTEGER, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: cldfra1_flag
3780 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
3781 CLDFRA
3783 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
3784 QV, &
3785 QI, &
3786 QC, &
3787 QS, &
3788 t_phy, &
3789 p_phy
3790 ! p_phy, &
3791 ! F_ICE_PHY, &
3792 ! F_RAIN_PHY
3794 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
3795 OPTIONAL, &
3796 INTENT(IN ) :: &
3797 F_ICE_PHY, &
3798 F_RAIN_PHY
3799 LOGICAL,OPTIONAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
3800 INTEGER :: mp_physics
3802 ! REAL thresh
3803 INTEGER:: i,j,k
3804 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
3806 REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
3807 PEXP=0.25, RHGRID=1.0
3808 REAL , PARAMETER :: SVP1=0.61078
3809 REAL , PARAMETER :: SVP2=17.2693882
3810 REAL , PARAMETER :: SVPI2=21.8745584
3811 REAL , PARAMETER :: SVP3=35.86
3812 REAL , PARAMETER :: SVPI3=7.66
3813 REAL , PARAMETER :: SVPT0=273.15
3814 REAL , PARAMETER :: r_d = 287.
3815 REAL , PARAMETER :: r_v = 461.6
3816 REAL , PARAMETER :: ep_2=r_d/r_v
3817 ! !DESCRIPTION:
3818 ! Compute cloud fraction from input ice and cloud water fields
3819 ! if provided.
3820 !
3821 ! Whether QI or QC is active or not is determined from the indices of
3822 ! the fields into the 4D scalar arrays in WRF. These indices are
3823 ! P_QI and P_QC, respectively, and they are passed in to the routine
3824 ! to enable testing to see if QI and QC represent active fields in
3825 ! the moisture 4D scalar array carried by WRF.
3826 !
3827 ! If a field is active its index will have a value greater than or
3828 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
3829 ! this routine.
3830 !EOP
3833 !-----------------------------------------------------------------------
3834 !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
3835 ! (modified by Ferrier, Feb '02)
3836 !
3837 !--- Cloud fraction parameterization follows Randall, 1994
3838 ! (see Hong et al., 1998)
3839 !-----------------------------------------------------------------------
3840 ! Note: ep_2=287./461.6 Rd/Rv
3841 ! Note: R_D=287.
3843 ! Alternative calculation for critical RH for grid saturation
3844 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
3846 ! Calculate saturation mixing ratio weighted according to the fractions of
3847 ! water and ice.
3848 ! Following:
3849 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
3850 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
3851 !
3852 ! over ice over water
3853 ! a = 21.8745584 17.2693882
3854 ! b = 7.66 35.86
3856 !---------------------------------------------------------------------
3858 DO j = jts,jte
3859 DO k = kts,kte
3860 DO i = its,ite
3861 tc = t_phy(i,k,j) - SVPT0
3862 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
3863 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
3864 QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
3865 QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
3867 ifouter: IF ( PRESENT(F_QI) .and. PRESENT(F_QC) .and. PRESENT(F_QS) ) THEN
3869 ! mji - For MP options 2, 4, 6, 7, 8, etc. (qc = liquid, qi = ice, qs = snow)
3870 IF ( F_QI .and. F_QC .and. F_QS) THEN
3871 QCLD = QI(i,k,j)+QC(i,k,j)+QS(i,k,j)
3872 IF (QCLD .LT. QCLDMIN) THEN
3873 weight = 0.
3874 ELSE
3875 weight = (QI(i,k,j)+QS(i,k,j)) / QCLD
3876 ENDIF
3877 ENDIF
3879 ! for P3, mp option 50 or 51
3880 IF ( F_QI .and. F_QC .and. .not. F_QS) THEN
3881 QCLD = QI(i,k,j)+QC(i,k,j)
3882 IF (QCLD .LT. QCLDMIN) THEN
3883 weight = 0.
3884 ELSE
3885 weight = (QI(i,k,j)) / QCLD
3886 ENDIF
3887 ENDIF
3889 ! mji - For MP options 1 and 3, (qc only)
3890 ! For MP=1, qc = liquid, for MP=3, qc = liquid or ice depending on temperature
3891 IF ( F_QC .and. .not. F_QI .and. .not. F_QS ) THEN
3892 QCLD = QC(i,k,j)
3893 IF (QCLD .LT. QCLDMIN) THEN
3894 weight = 0.
3895 ELSE
3896 if (t_phy(i,k,j) .gt. 273.15) weight = 0.
3897 if (t_phy(i,k,j) .le. 273.15) weight = 1.
3898 ENDIF
3899 ENDIF
3901 ! mji - For MP option 5; (qc = liquid, qs = ice)
3902 IF ( F_QC .and. .not. F_QI .and. F_QS .and. PRESENT(F_ICE_PHY) ) THEN
3904 ! Mixing ratios of cloud water & total ice (cloud ice + snow).
3905 ! Mixing ratios of rain are not considered in this scheme.
3906 ! F_ICE is fraction of ice
3907 ! F_RAIN is fraction of rain
3909 QIMID = QS(i,k,j)
3910 QWMID = QC(i,k,j)
3911 ! old method
3912 ! QIMID = QC(i,k,j)*F_ICE_PHY(i,k,j)
3913 ! QWMID = (QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
3914 !
3915 !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
3916 ! only cloud water + cloud ice + snow
3917 !
3918 QCLD=QWMID+QIMID
3919 IF (QCLD .LT. QCLDMIN) THEN
3920 weight = 0.
3921 ELSE
3922 weight = F_ICE_PHY(i,k,j)
3923 ENDIF
3924 ENDIF
3925 ! IF ( F_QC .and. F_QI .and. .not. F_QS ) THEN
3926 IF ( mp_physics .eq. FER_MP_HIRES .or. &
3927 mp_physics==fer_mp_hires_advect) THEN
3928 QIMID = QI(i,k,j) !- total ice (cloud ice + snow)
3929 QWMID = QC(i,k,j) !- cloud water
3930 QCLD=QWMID+QIMID !- cloud water + total ice
3931 IF (QCLD .LT. QCLDMIN) THEN
3932 weight = 0.
3933 ELSE
3934 weight = QIMID/QCLD
3935 if (tc<-40.) weight=1.
3936 ENDIF
3937 ENDIF
3939 ELSE
3940 CLDFRA(i,k,j)=0.
3942 ENDIF ifouter ! IF ( F_QI .and. F_QC .and. F_QS)
3945 QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI
3946 RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
3947 !
3948 !--- Determine cloud fraction (modified from original algorithm)
3949 !
3950 cldfra1_flag(i,k,j) = 0
3951 IF (QCLD .LT. QCLDMIN) THEN
3952 !
3953 !--- Assume zero cloud fraction if there is no cloud mixing ratio
3954 !
3955 CLDFRA(i,k,j)=0.
3956 cldfra1_flag(i,k,j) = 1
3957 ELSEIF(RHUM.GE.RHGRID)THEN
3958 !
3959 !--- Assume cloud fraction of unity if near saturation and the cloud
3960 ! mixing ratio is at or above the minimum threshold
3961 !
3962 CLDFRA(i,k,j)=1.
3963 cldfra1_flag(i,k,j) = 2
3964 ELSE
3965 cldfra1_flag(i,k,j) = 3
3966 !
3967 !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
3968 ! modified based on assumed grid-scale saturation at RH=RHgrid.
3969 !
3970 SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
3971 DENOM=(SUBSAT)**GAMMA
3972 ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
3973 ! prevent negative values (new)
3974 RHUM=MAX(1.E-10, RHUM)
3975 CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG))
3976 !! ARG=-1000*QCLD/(RHUM-RHGRID)
3977 !! ARG=MAX(ARG, ARGMIN)
3978 !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
3979 IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
3981 ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
3982 ENDDO !--- End DO i
3983 ENDDO !--- End DO k
3984 ENDDO !--- End DO j
3986 END SUBROUTINE cal_cldfra1
3988 !+---+-----------------------------------------------------------------+
3989 !..Cloud fraction scheme by G. Thompson (NCAR-RAL), not intended for
3990 !.. combining with any cumulus or shallow cumulus parameterization
3991 !.. scheme cloud fractions. This is intended as a stand-alone for
3992 !.. cloud fraction and is relatively good at getting widespread stratus
3993 !.. and stratoCu without caring whether any deep/shallow Cu param schemes
3994 !.. is making sub-grid-spacing clouds/precip. Under the hood, this
3995 !.. scheme follows Mocko and Cotton (1995) in applicaiton of the
3996 !.. Sundqvist et al (1989) scheme but using a grid-scale dependent
3997 !.. RH threshold, one each for land v. ocean points based on
3998 !.. experiences with HWRF testing.
3999 !+---+-----------------------------------------------------------------+
4000 !
4001 !+---+-----------------------------------------------------------------+
4003 SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, &
4004 & p, t, XLAND, gridkm, &
4005 & modify_qvapor, max_relh, &
4006 & kts,kte, debug_flag)
4007 !
4008 USE module_mp_thompson , ONLY : rsif, rslf
4009 IMPLICIT NONE
4010 !
4011 INTEGER, INTENT(IN):: kts, kte
4012 LOGICAL, INTENT(IN):: modify_qvapor
4013 REAL, DIMENSION(kts:kte), INTENT(INOUT):: qv, qc, qi, cldfra
4014 REAL, DIMENSION(kts:kte), INTENT(IN):: p, t, dz, qs
4015 REAL, INTENT(IN):: gridkm, XLAND, max_relh
4016 LOGICAL, INTENT(IN):: debug_flag
4018 !..Local vars.
4019 REAL:: RH_00L, RH_00O, RH_00
4020 REAL:: entrmnt=0.5
4021 INTEGER:: k
4022 REAL:: TC, qvsi, qvsw, RHUM, delz
4023 REAL, DIMENSION(kts:kte):: qvs, rh, rhoa
4025 character*512 dbg_msg
4027 !+---+
4029 !..Initialize cloud fraction, compute RH, and rho-air.
4031 DO k = kts,kte
4032 CLDFRA(K) = 0.0
4034 qvsw = rslf(P(k), t(k))
4035 qvsi = rsif(P(k), t(k))
4037 tc = t(k) - 273.15
4038 if (tc .ge. -12.0) then
4039 qvs(k) = qvsw
4040 elseif (tc .lt. -35.0) then
4041 qvs(k) = qvsi
4042 else
4043 qvs(k) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+35.)
4044 endif
4046 if (modify_qvapor) then
4047 if (qc(k).gt.1.E-8) then
4048 qv(k) = MAX(qv(k), qvsw)
4049 qvs(k) = qvsw
4050 endif
4051 if (qc(k).le.1.E-8 .and. qi(k).ge.1.E-9) then
4052 qv(k) = MAX(qv(k), qvsi*1.005) !..To ensure a tiny bit ice supersaturation
4053 qvs(k) = qvsi
4054 endif
4055 endif
4057 rh(k) = MAX(0.01, qv(k)/qvs(k))
4058 rhoa(k) = p(k)/(287.0*t(k))
4059 ENDDO
4062 !..First cut scale-aware. Higher resolution should require closer to
4063 !.. saturated grid box for higher cloud fraction. Simple functions
4064 !.. chosen based on Mocko and Cotton (1995) starting point and desire
4065 !.. to get near 100% RH as grid spacing moves toward 1.0km, but higher
4066 !.. RH over ocean required as compared to over land.
4068 DO k = kts,kte
4070 delz = MAX(100., dz(k))
4071 RH_00L = 0.77 + MIN(0.22,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
4072 RH_00O = 0.85 + MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
4073 RHUM = rh(k)
4075 if (qc(k).gt.1.E-6 .or. qi(k).ge.1.E-7 &
4076 & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then
4077 CLDFRA(K) = 1.0
4078 elseif (((qc(k)+qi(k)).gt.1.E-10) .and. &
4079 & ((qc(k)+qi(k)).lt.1.E-6)) then
4080 CLDFRA(K) = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k))))
4081 else
4083 IF ((XLAND-1.5).GT.0.) THEN !--- Ocean
4084 RH_00 = RH_00O
4085 ELSE !--- Land
4086 RH_00 = RH_00L
4087 ENDIF
4089 tc = MAX(-80.0, t(k) - 273.15)
4090 if (tc .lt. -12.0) RH_00 = RH_00L
4092 if (tc .ge. 25.0) then
4093 CLDFRA(K) = 0.0
4094 elseif (tc .ge. -12.0) then
4095 RHUM = MIN(rh(k), 1.0)
4096 CLDFRA(K) = MAX(0., 1.0-SQRT((1.001-RHUM)/(1.001-RH_00)))
4097 else
4098 if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then
4099 !..For HRRR model, the following look OK.
4100 RHUM = MIN(rh(k), 1.45)
4101 RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.)
4102 CLDFRA(K) = MAX(0., 1.0-SQRT((1.46-RHUM)/(1.46-RH_00)))
4103 else
4104 !..but for the GFS model, RH is way lower.
4105 RHUM = MIN(rh(k), 1.05)
4106 RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.)
4107 CLDFRA(K) = MAX(0., 1.0-SQRT((1.06-RHUM)/(1.06-RH_00)))
4108 endif
4109 endif
4110 if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.99))
4112 if (debug_flag) then
4113 WRITE (dbg_msg,*) 'DEBUG-GT: cloud fraction (k,RH_00, RHUM, CF): ',k,RH_00,RHUM,CLDFRA(K)
4114 CALL wrf_debug (150, dbg_msg)
4115 endif
4117 endif
4118 if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0
4119 ENDDO
4122 call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt, &
4123 & debug_flag, qc, qi, qs, kts,kte)
4125 !..Do a final total column adjustment since we may have added more than 1mm
4126 !.. LWP/IWP for multiple cloud decks.
4128 call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte)
4130 if (debug_flag) then
4131 WRITE (dbg_msg,*) 'DEBUG-GT: Made-up fake profile of clouds'
4132 CALL wrf_debug (150, dbg_msg)
4133 do k = kte, kts, -1
4134 write(dbg_msg,'(f7.2,2x,f7.2,2x,f6.4,2x,f7.3,x,f15.7,x,f15.7,x,f15.7)') &
4135 & T(k)-273.15, P(k)*0.01, rh(k), cldfra(k)*100., qc(k)*1000.,qi(k)*1000., qs(k)*1000.
4136 CALL wrf_debug (150, dbg_msg)
4137 enddo
4138 endif
4141 if (modify_qvapor) then
4142 DO k = kts,kte
4143 if (cldfra(k).gt.0.20 .and. cldfra(k).lt.1.0) then
4144 qv(k) = MAX(qv(k),qvs(k))
4145 endif
4146 ENDDO
4147 endif
4150 END SUBROUTINE cal_cldfra3
4152 !+---+-----------------------------------------------------------------+
4153 !..From cloud fraction array, find clouds of multi-level depth and compute
4154 !.. a reasonable value of LWP or IWP that might be contained in that depth,
4155 !.. unless existing LWC/IWC is already there.
4157 SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,&
4158 & debugfl, qc1d, qi1d, qs1d, kts,kte)
4159 !
4160 IMPLICIT NONE
4161 !
4162 INTEGER, INTENT(IN):: kts, kte
4163 LOGICAL, INTENT(IN):: debugfl
4164 REAL, INTENT(IN):: entrmnt
4165 REAL, DIMENSION(kts:kte), INTENT(IN):: qs1d,qvs1d,T1d,P1d,Dz1d
4166 REAL, DIMENSION(kts:kte), INTENT(INOUT):: cfr1d, qc1d, qi1d
4168 !..Local vars.
4169 REAL, DIMENSION(kts:kte):: theta
4170 REAL:: theta1, theta2, delz
4171 INTEGER:: k, k2, k_tropo, k_m12C, k_cldb, k_cldt, kbot, k_p150
4172 LOGICAL:: in_cloud
4173 character*512 dbg_msg
4175 !+---+
4177 k_m12C = 0
4178 k_p150 = 0
4179 DO k = kte, kts, -1
4180 theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.))
4181 if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10100.0) k_m12C = MAX(k_m12C, k)
4182 if (P1d(k).gt.14999.0 .and. k_p150.eq.0) k_p150 = k
4183 ENDDO
4184 if (k_m12C .le. kts) k_m12C = kts
4186 if (k_m12C.gt.kte-3) then
4187 WRITE (dbg_msg,*) 'DEBUG-GT: WARNING, no possible way neg12C can occur this high up: ', k_m12C
4188 CALL wrf_debug (0, dbg_msg)
4189 do k = kte, kts, -1
4190 WRITE (dbg_msg,*) 'DEBUG-GT, k, P, T : ', k,P1d(k)*0.01,T1d(k)-273.15
4191 CALL wrf_debug (0, dbg_msg)
4192 enddo
4193 call wrf_error_fatal ('FATAL ERROR, problem in temperature profile.')
4194 endif
4196 !..Find tropopause height, best surrogate, because we would not really
4197 !.. wish to put fake clouds into the stratosphere. The 10/1500 ratio
4198 !.. d(Theta)/d(Z) approximates a vertical line on typical SkewT chart
4199 !.. near typical (mid-latitude) tropopause height. Since messy data
4200 !.. could give us a false signal of such a transition, do the check over
4201 !.. three K-level change, not just a level-to-level check. This method
4202 !.. has potential failure in arctic-like conditions with extremely low
4203 !.. tropopause height, as would any other diagnostic, so ensure resulting
4204 !.. k_tropo level is above 700hPa.
4206 if ( (kte-k_p150) .lt. 3) k_p150 = kte-3
4207 DO k = k_p150-2, kts, -1
4208 theta1 = theta(k)
4209 theta2 = theta(k+2)
4210 delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2)
4211 if ( (((theta2-theta1)/delz).lt.10./1500.) .OR. P1d(k).gt.70000.) EXIT
4212 ENDDO
4213 k_tropo = MAX(kts+2, MIN(k+2, kte-1))
4215 if (k_tropo .gt. k_p150) then
4216 DO k = kte-3, k_p150-2, -1
4217 theta1 = theta(k)
4218 theta2 = theta(k+2)
4219 delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2)
4220 if ( (((theta2-theta1)/delz).lt.10./1500.) .AND. P1d(k).gt.9500.) EXIT
4221 ENDDO
4222 k_tropo = MAX(k_p150-1, MIN(k+2, kte-1))
4223 endif
4225 if (k_tropo.gt.kte-2) then
4226 WRITE (dbg_msg,*) 'DEBUG-GT: CAUTION, tropopause appears to be very high up: ', k_tropo
4227 CALL wrf_debug (150, dbg_msg)
4228 do k = kte, kts, -1
4229 WRITE (dbg_msg,*) 'DEBUG-GT, P, T : ', k,P1d(k)*0.01,T1d(k)-273.16
4230 CALL wrf_debug (150, dbg_msg)
4231 enddo
4232 elseif (debugfl) then
4233 WRITE (dbg_msg,*) 'DEBUG-GT: FOUND TROPOPAUSE k=', k_tropo
4234 CALL wrf_debug (150, dbg_msg)
4235 endif
4237 !..Eliminate possible fractional clouds above supposed tropopause.
4238 DO k = k_tropo+1, kte
4239 if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) then
4240 cfr1d(k) = 0.
4241 endif
4242 ENDDO
4244 !..Be a bit more conservative with lower cloud fraction in scenario with
4245 !.. well-mixed convective boundary layer below LCL.
4247 kbot = kts+1
4248 DO k = kbot, k_m12C
4249 if ( (theta(k)-theta(k-1)) .gt. 0.010E-3*Dz1d(k)) EXIT
4250 ENDDO
4251 kbot = MAX(kts+1, k-2)
4252 DO k = kts, kbot
4253 if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) cfr1d(k) = MAX(0.01,0.5*cfr1d(k))
4254 ENDDO
4255 DO k = kts,k_tropo
4256 if (cfr1d(k).gt.0.0) kbot = MIN(k,kbot)
4257 ENDDO
4259 !..Starting below tropo height, if cloud fraction greater than 1 percent,
4260 !.. compute an approximate total layer depth of cloud, determine a total
4261 !.. liquid water/ice path (LWP/IWP), then reduce that amount with tuning
4262 !.. parameter to represent entrainment factor, then divide up LWP/IWP
4263 !.. into delta-Z weighted amounts for individual levels per cloud layer.
4265 k_cldb = k_tropo
4266 in_cloud = .false.
4267 k = k_tropo
4268 DO WHILE (.not. in_cloud .AND. k.gt.k_m12C+1)
4269 k_cldt = 0
4270 if (cfr1d(k).ge.0.01) then
4271 in_cloud = .true.
4272 k_cldt = MAX(k_cldt, k)
4273 endif
4274 if (in_cloud) then
4275 DO k2 = k_cldt-1, k_m12C, -1
4276 if (cfr1d(k2).lt.0.01 .or. k2.eq.k_m12C) then
4277 k_cldb = k2+1
4278 goto 87
4279 endif
4280 ENDDO
4281 87 continue
4282 in_cloud = .false.
4283 endif
4284 if ((k_cldt - k_cldb + 1) .ge. 2) then
4285 if (debugfl) then
4286 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
4287 CALL wrf_debug (150, dbg_msg)
4288 endif
4289 call adjust_cloudIce(cfr1d, qi1d, qs1d, qvs1d, T1d, Dz1d, &
4290 & entrmnt, k_cldb,k_cldt,kts,kte)
4291 k = k_cldb
4292 elseif ((k_cldt - k_cldb + 1) .eq. 1) then
4293 if (debugfl) then
4294 WRITE (dbg_msg,*) 'DEBUG-GT: A single-layer ice cloud layer is found on ', k_cldb, P1d(k_cldb)*0.01
4295 CALL wrf_debug (150, dbg_msg)
4296 endif
4297 if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) &
4298 & qi1d(k_cldb)=qi1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb)
4299 k = k_cldb
4300 endif
4301 k = k - 1
4302 ENDDO
4305 k_cldb = k_m12C + 3
4306 in_cloud = .false.
4307 k = k_m12C + 2
4308 DO WHILE (.not. in_cloud .AND. k.gt.kbot)
4309 k_cldt = 0
4310 if (cfr1d(k).ge.0.01) then
4311 in_cloud = .true.
4312 k_cldt = MAX(k_cldt, k)
4313 endif
4314 if (in_cloud) then
4315 DO k2 = k_cldt-1, kbot, -1
4316 if (cfr1d(k2).lt.0.01 .or. k2.eq.kbot) then
4317 k_cldb = k2+1
4318 goto 88
4319 endif
4320 ENDDO
4321 88 continue
4322 in_cloud = .false.
4323 endif
4324 if ((k_cldt - k_cldb + 1) .ge. 2) then
4325 if (debugfl) then
4326 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
4327 CALL wrf_debug (150, dbg_msg)
4328 endif
4329 call adjust_cloudH2O(cfr1d, qc1d, qvs1d, T1d, Dz1d, &
4330 & entrmnt, k_cldb,k_cldt,kts,kte)
4331 k = k_cldb
4332 elseif ((k_cldt - k_cldb + 1) .eq. 1) then
4333 if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) &
4334 & qc1d(k_cldb)=qc1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb)
4335 k = k_cldb
4336 endif
4337 k = k - 1
4338 ENDDO
4341 END SUBROUTINE find_cloudLayers
4343 !+---+-----------------------------------------------------------------+
4345 SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte)
4346 !
4347 IMPLICIT NONE
4348 !
4349 INTEGER, INTENT(IN):: k1,k2, kts,kte
4350 REAL, INTENT(IN):: entr
4351 REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qs, qvs, T, dz
4352 REAL, DIMENSION(kts:kte), INTENT(INOUT):: qi
4353 REAL:: iwc, max_iwc, tdz, this_iwc, this_dz
4354 INTEGER:: k
4356 tdz = 0.
4357 do k = k1, k2
4358 tdz = tdz + dz(k)
4359 enddo
4360 ! max_iwc = ABS(qvs(k2)-qvs(k1))
4361 max_iwc = MAX(0.0, qvs(k1)-qvs(k2))
4362 ! print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz
4364 do k = k1, k2
4365 max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k)))
4366 enddo
4367 max_iwc = MIN(1.E-4, max_iwc)
4369 this_dz = 0.0
4370 do k = k1, k2
4371 if (k.eq.k1) then
4372 this_dz = this_dz + 0.5*dz(k)
4373 else
4374 this_dz = this_dz + dz(k)
4375 endif
4376 this_iwc = max_iwc*this_dz/tdz
4377 iwc = MAX(1.E-6, this_iwc*(1.-entr))
4378 if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then
4379 qi(k) = qi(k) + cfr(k)*cfr(k)*iwc
4380 endif
4381 enddo
4383 END SUBROUTINE adjust_cloudIce
4385 !+---+-----------------------------------------------------------------+
4387 SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte)
4388 !
4389 IMPLICIT NONE
4390 !
4391 INTEGER, INTENT(IN):: k1,k2, kts,kte
4392 REAL, INTENT(IN):: entr
4393 REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qvs, T, dz
4394 REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc
4395 REAL:: lwc, max_lwc, tdz, this_lwc, this_dz
4396 INTEGER:: k
4398 tdz = 0.
4399 do k = k1, k2
4400 tdz = tdz + dz(k)
4401 enddo
4402 ! max_lwc = ABS(qvs(k2)-qvs(k1))
4403 max_lwc = MAX(0.0, qvs(k1)-qvs(k2))
4404 ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz
4406 do k = k1, k2
4407 max_lwc = MAX(1.E-6, max_lwc - qc(k))
4408 enddo
4409 max_lwc = MIN(1.E-4, max_lwc)
4411 this_dz = 0.0
4412 do k = k1, k2
4413 if (k.eq.k1) then
4414 this_dz = this_dz + 0.5*dz(k)
4415 else
4416 this_dz = this_dz + dz(k)
4417 endif
4418 this_lwc = max_lwc*this_dz/tdz
4419 lwc = MAX(1.E-6, this_lwc*(1.-entr))
4420 if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.258.16) then
4421 qc(k) = qc(k) + cfr(k)*cfr(k)*lwc
4422 endif
4423 enddo
4425 END SUBROUTINE adjust_cloudH2O
4427 !+---+-----------------------------------------------------------------+
4429 !..Do not alter any grid-explicitly resolved hydrometeors, rather only
4430 !.. the supposed amounts due to the cloud fraction scheme.
4432 SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte)
4433 !
4434 IMPLICIT NONE
4435 !
4436 INTEGER, INTENT(IN):: kts,kte
4437 REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, Rho, dz
4438 REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc, qi
4439 REAL:: lwp, iwp, xfac
4440 INTEGER:: k
4442 lwp = 0.
4443 iwp = 0.
4444 do k = kts, kte
4445 if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
4446 lwp = lwp + qc(k)*Rho(k)*dz(k)
4447 iwp = iwp + qi(k)*Rho(k)*dz(k)
4448 endif
4449 enddo
4451 if (lwp .gt. 1.0) then
4452 xfac = 1.0/lwp
4453 do k = kts, kte
4454 if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
4455 qc(k) = qc(k)*xfac
4456 endif
4457 enddo
4458 endif
4460 if (iwp .gt. 1.0) then
4461 xfac = 1.0/iwp
4462 do k = kts, kte
4463 if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
4464 qi(k) = qi(k)*xfac
4465 endif
4466 enddo
4467 endif
4469 END SUBROUTINE adjust_cloudFinal
4472 !+---+-----------------------------------------------------------------+
4474 SUBROUTINE toposhad_init(ht_shad,ht_loc,shadowmask,nested,iter, &
4475 ids,ide, jds,jde, kds,kde, &
4476 ims,ime, jms,jme, kms,kme, &
4477 ips,ipe, jps,jpe, kps,kpe, &
4478 its,ite, jts,jte, kts,kte )
4480 USE module_model_constants
4482 implicit none
4484 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4485 ims,ime, jms,jme, kms,kme, &
4486 ips,ipe, jps,jpe, kps,kpe, &
4487 its,ite, jts,jte, kts,kte
4489 LOGICAL, INTENT(IN) :: nested
4491 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad, ht_loc
4493 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
4494 INTEGER, INTENT(IN) :: iter
4496 ! Local variables
4498 INTEGER :: i, j
4500 if (iter.eq.1) then
4502 ! Initialize shadow mask
4503 do j=jts,jte
4504 do i=its,ite
4505 shadowmask(i,j) = 0
4506 ENDDO
4507 ENDDO
4509 ! Initialize shading height
4511 IF ( nested ) THEN ! Do not overwrite input from parent domain
4512 do j=max(jts,jds+2),min(jte,jde-3)
4513 do i=max(its,ids+2),min(ite,ide-3)
4514 ht_shad(i,j) = ht_loc(i,j)-0.001
4515 ENDDO
4516 ENDDO
4517 ELSE
4518 do j=jts,jte
4519 do i=its,ite
4520 ht_shad(i,j) = ht_loc(i,j)-0.001
4521 ENDDO
4522 ENDDO
4523 ENDIF
4525 IF ( nested ) THEN ! Check if a shadow exceeding the topography height is available at the lateral domain edge from nesting
4526 if (its.eq.ids) then
4527 do j=jts,jte
4528 if (ht_shad(its,j) .gt. ht_loc(its,j)) then
4529 shadowmask(its,j) = 1
4530 ht_loc(its,j) = ht_shad(its,j)
4531 endif
4532 if (ht_shad(its+1,j) .gt. ht_loc(its+1,j)) then
4533 shadowmask(its+1,j) = 1
4534 ht_loc(its+1,j) = ht_shad(its+1,j)
4535 endif
4536 enddo
4537 endif
4538 if (ite.eq.ide-1) then
4539 do j=jts,jte
4540 if (ht_shad(ite,j) .gt. ht_loc(ite,j)) then
4541 shadowmask(ite,j) = 1
4542 ht_loc(ite,j) = ht_shad(ite,j)
4543 endif
4544 if (ht_shad(ite-1,j) .gt. ht_loc(ite-1,j)) then
4545 shadowmask(ite-1,j) = 1
4546 ht_loc(ite-1,j) = ht_shad(ite-1,j)
4547 endif
4548 enddo
4549 endif
4550 if (jts.eq.jds) then
4551 do i=its,ite
4552 if (ht_shad(i,jts) .gt. ht_loc(i,jts)) then
4553 shadowmask(i,jts) = 1
4554 ht_loc(i,jts) = ht_shad(i,jts)
4555 endif
4556 if (ht_shad(i,jts+1) .gt. ht_loc(i,jts+1)) then
4557 shadowmask(i,jts+1) = 1
4558 ht_loc(i,jts+1) = ht_shad(i,jts+1)
4559 endif
4560 enddo
4561 endif
4562 if (jte.eq.jde-1) then
4563 do i=its,ite
4564 if (ht_shad(i,jte) .gt. ht_loc(i,jte)) then
4565 shadowmask(i,jte) = 1
4566 ht_loc(i,jte) = ht_shad(i,jte)
4567 endif
4568 if (ht_shad(i,jte-1) .gt. ht_loc(i,jte-1)) then
4569 shadowmask(i,jte-1) = 1
4570 ht_loc(i,jte-1) = ht_shad(i,jte-1)
4571 endif
4572 enddo
4573 endif
4574 ENDIF
4576 else
4578 ! Fill the local topography field at the points next to internal tile boundaries with ht_shad values
4579 ! A 2-pt halo has been applied to the ht_shad before the repeated call of this subroutine
4581 if ((its.ne.ids).and.(its.eq.ips)) then
4582 do j=jts-2,jte+2
4583 ht_loc(its-1,j) = max(ht_loc(its-1,j),ht_shad(its-1,j))
4584 ht_loc(its-2,j) = max(ht_loc(its-2,j),ht_shad(its-2,j))
4585 enddo
4586 endif
4587 if ((ite.ne.ide-1).and.(ite.eq.ipe)) then
4588 do j=jts-2,jte+2
4589 ht_loc(ite+1,j) = max(ht_loc(ite+1,j),ht_shad(ite+1,j))
4590 ht_loc(ite+2,j) = max(ht_loc(ite+2,j),ht_shad(ite+2,j))
4591 enddo
4592 endif
4593 if ((jts.ne.jds).and.(jts.eq.jps)) then
4594 do i=its-2,ite+2
4595 ht_loc(i,jts-1) = max(ht_loc(i,jts-1),ht_shad(i,jts-1))
4596 ht_loc(i,jts-2) = max(ht_loc(i,jts-2),ht_shad(i,jts-2))
4597 enddo
4598 endif
4599 if ((jte.ne.jde-1).and.(jte.eq.jpe)) then
4600 do i=its-2,ite+2
4601 ht_loc(i,jte+1) = max(ht_loc(i,jte+1),ht_shad(i,jte+1))
4602 ht_loc(i,jte+2) = max(ht_loc(i,jte+2),ht_shad(i,jte+2))
4603 enddo
4604 endif
4606 endif
4608 END SUBROUTINE toposhad_init
4610 SUBROUTINE toposhad(xlat,xlong,sina,cosa,xtime,gmt,radfrq,declin, &
4611 dx,dy,ht_shad,ht_loc,iter, &
4612 shadowmask,shadlen, &
4613 ids,ide, jds,jde, kds,kde, &
4614 ims,ime, jms,jme, kms,kme, &
4615 ips,ipe, jps,jpe, kps,kpe, &
4616 its,ite, jts,jte, kts,kte )
4619 USE module_model_constants
4621 implicit none
4623 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4624 ims,ime, jms,jme, kms,kme, &
4625 ips,ipe, jps,jpe, kps,kpe, &
4626 its,ite, jts,jte, kts,kte
4628 INTEGER, INTENT(IN) :: iter
4630 REAL, INTENT(IN) :: RADFRQ,XTIME,DECLIN,dx,dy,gmt,shadlen
4632 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT, XLONG, sina, cosa
4634 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
4636 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
4638 ! Local variables
4640 REAL :: pi, xt24, wgt, ri, rj, argu, sol_azi, topoelev, dxabs, tloctm, hrang, xxlat, csza
4641 INTEGER :: gpshad, ii, jj, i1, i2, j1, j2, i, j
4645 XT24=MOD(XTIME+RADFRQ*0.5,1440.)
4646 pi = 4.*atan(1.)
4647 gpshad = int(shadlen/dx+1.)
4649 if (iter.eq.1) then
4652 j_loop1: DO J=jts,jte
4653 i_loop1: DO I=its,ite
4655 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
4656 HRANG=15.*(TLOCTM-12.)*DEGRAD
4657 XXLAT=XLAT(i,j)*DEGRAD
4658 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
4660 if (csza.lt.1.e-2) then ! shadow mask does not need to be computed
4661 shadowmask(i,j) = 0
4662 ht_shad(i,j) = ht_loc(i,j)-0.001
4663 goto 120
4664 endif
4666 ! Solar azimuth angle
4668 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
4669 if (argu.gt.1) argu = 1
4670 if (argu.lt.-1) argu = -1
4671 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
4672 if (cosa(i,j).ge.0) then
4673 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
4674 else
4675 sol_azi = sol_azi + pi - asin(sina(i,j))
4676 endif
4678 ! Scan for higher surrounding topography
4680 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
4682 do jj = j+1,j+gpshad
4683 ri = i + (jj-j)*tan(sol_azi)
4684 i1 = int(ri)
4685 i2 = i1+1
4686 wgt = ri-i1
4687 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
4688 if ((jj.ge.jpe+3).or.(i1.le.ips-3).or.(i2.ge.ipe+3)) then
4689 ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
4690 goto 120
4691 endif
4692 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
4693 if (sin(topoelev).ge.csza) then
4694 shadowmask(i,j) = 1
4695 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4696 endif
4697 enddo
4699 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
4700 do ii = i+1,i+gpshad
4701 rj = j - (ii-i)*tan(pi/2.+sol_azi)
4702 j1 = int(rj)
4703 j2 = j1+1
4704 wgt = rj-j1
4705 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
4706 if ((ii.ge.ipe+3).or.(j1.le.jps-3).or.(j2.ge.jpe+3)) then
4707 ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
4708 goto 120
4709 endif
4710 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
4711 if (sin(topoelev).ge.csza) then
4712 shadowmask(i,j) = 1
4713 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4714 endif
4715 enddo
4717 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
4718 do jj = j-1,j-gpshad,-1
4719 ri = i + (jj-j)*tan(sol_azi)
4720 i1 = int(ri)
4721 i2 = i1+1
4722 wgt = ri-i1
4723 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
4724 if ((jj.le.jps-3).or.(i1.le.ips-3).or.(i2.ge.ipe+3)) then
4725 ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
4726 goto 120
4727 endif
4728 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
4729 if (sin(topoelev).ge.csza) then
4730 shadowmask(i,j) = 1
4731 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4732 endif
4733 enddo
4735 else ! sun is in the western quarter
4736 do ii = i-1,i-gpshad,-1
4737 rj = j - (ii-i)*tan(pi/2.+sol_azi)
4738 j1 = int(rj)
4739 j2 = j1+1
4740 wgt = rj-j1
4741 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
4742 if ((ii.le.ips-3).or.(j1.le.jps-3).or.(j2.ge.jpe+3)) then
4743 ! if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
4744 goto 120
4745 endif
4746 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
4747 if (sin(topoelev).ge.csza) then
4748 shadowmask(i,j) = 1
4749 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4750 endif
4751 enddo
4752 endif
4754 120 continue
4756 ENDDO i_loop1
4757 ENDDO j_loop1
4759 else ! iteration > 1
4762 j_loop2: DO J=jts,jte
4763 i_loop2: DO I=its,ite
4765 ! if (shadowmask(i,j).eq.-1) then ! this indicates that the search ended at a lateral boundary during iteration 1
4767 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
4768 HRANG=15.*(TLOCTM-12.)*DEGRAD
4769 XXLAT=XLAT(i,j)*DEGRAD
4770 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
4772 if (csza.lt.1.e-2) then ! shadow mask does not need to be computed
4773 shadowmask(i,j) = 0
4774 ht_shad(i,j) = ht_loc(i,j)-0.001
4775 goto 220
4776 endif
4778 ! Solar azimuth angle
4780 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
4781 if (argu.gt.1) argu = 1
4782 if (argu.lt.-1) argu = -1
4783 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
4784 if (cosa(i,j).ge.0) then
4785 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
4786 else
4787 sol_azi = sol_azi + pi - asin(sina(i,j))
4788 endif
4790 ! Scan for higher surrounding topography
4792 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
4794 do jj = j+1,j+gpshad
4795 ri = i + (jj-j)*tan(sol_azi)
4796 i1 = int(ri)
4797 i2 = i1+1
4798 wgt = ri-i1
4799 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
4800 if ((jj.ge.min(jde,jpe+3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220
4801 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
4802 if (sin(topoelev).ge.csza) then
4803 shadowmask(i,j) = 1
4804 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4805 endif
4806 enddo
4808 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
4809 do ii = i+1,i+gpshad
4810 rj = j - (ii-i)*tan(pi/2.+sol_azi)
4811 j1 = int(rj)
4812 j2 = j1+1
4813 wgt = rj-j1
4814 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
4815 if ((ii.ge.min(ide,ipe+3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220
4816 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
4817 if (sin(topoelev).ge.csza) then
4818 shadowmask(i,j) = 1
4819 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4820 endif
4821 enddo
4823 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
4824 do jj = j-1,j-gpshad,-1
4825 ri = i + (jj-j)*tan(sol_azi)
4826 i1 = int(ri)
4827 i2 = i1+1
4828 wgt = ri-i1
4829 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
4830 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
4831 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
4832 if (sin(topoelev).ge.csza) then
4833 shadowmask(i,j) = 1
4834 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4835 endif
4836 enddo
4838 else ! sun is in the western quarter
4839 do ii = i-1,i-gpshad,-1
4840 rj = j - (ii-i)*tan(pi/2.+sol_azi)
4841 j1 = int(rj)
4842 j2 = j1+1
4843 wgt = rj-j1
4844 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
4845 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
4846 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
4847 if (sin(topoelev).ge.csza) then
4848 shadowmask(i,j) = 1
4849 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
4850 endif
4851 enddo
4852 endif
4854 220 continue
4855 ! endif
4857 ENDDO i_loop2
4858 ENDDO j_loop2
4860 endif ! iteration
4862 END SUBROUTINE toposhad
4864 SUBROUTINE ozn_time_int(julday,julian,ozmixm,ozmixt,levsiz,num_months, &
4865 ids , ide , jds , jde , kds , kde , &
4866 ims , ime , jms , jme , kms , kme , &
4867 its , ite , jts , jte , kts , kte )
4869 ! adapted from oznint from CAM module
4870 ! input: ozmixm - read from physics_init
4871 ! output: ozmixt - time interpolated
4873 ! USE module_ra_cam_support, ONLY : getfactors
4875 IMPLICIT NONE
4877 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4878 ims,ime, jms,jme, kms,kme, &
4879 its,ite, jts,jte, kts,kte
4881 INTEGER, INTENT(IN ) :: levsiz, num_months
4883 REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months ), &
4884 INTENT(IN ) :: ozmixm
4886 INTEGER, INTENT(IN ) :: JULDAY
4887 REAL, INTENT(IN ) :: JULIAN
4889 REAL, DIMENSION( ims:ime, levsiz, jms:jme ), &
4890 INTENT(OUT ) :: ozmixt
4892 !Local
4893 REAL :: intJULIAN
4894 integer :: np1,np,nm,m,k,i,j
4895 integer :: IJUL
4896 integer, dimension(12) :: date_oz
4897 data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/
4898 real, parameter :: daysperyear = 365. ! number of days in a year
4899 real :: cdayozp, cdayozm
4900 real :: fact1, fact2, deltat
4901 logical :: finddate
4902 logical :: ozncyc
4903 CHARACTER(LEN=256) :: msgstr
4905 ozncyc = .true.
4906 ! JULIAN starts from 0.0 at 0Z on 1 Jan.
4907 intJULIAN = JULIAN + 1.0 ! offset by one day
4908 ! jan 1st 00z is julian=1.0 here
4909 IJUL=INT(intJULIAN)
4910 ! Note that following will drift.
4911 ! Need to use actual month/day info to compute julian.
4912 intJULIAN=intJULIAN-FLOAT(IJUL)
4913 IJUL=MOD(IJUL,365)
4914 IF(IJUL.EQ.0)IJUL=365
4915 intJULIAN=intJULIAN+IJUL
4916 np1=1
4917 finddate=.false.
4919 ! do m=1,num_months
4920 do m=1,12
4921 if(date_oz(m).gt.intjulian.and..not.finddate) then
4922 np1=m
4923 finddate=.true.
4924 endif
4925 enddo
4926 cdayozp=date_oz(np1)
4928 if(np1.gt.1) then
4929 cdayozm=date_oz(np1-1)
4930 np=np1
4931 nm=np-1
4932 else
4933 cdayozm=date_oz(12)
4934 np=np1
4935 nm=12
4936 endif
4938 ! call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, &
4939 ! fact1, fact2)
4940 !
4941 ! Determine time interpolation factors. Account for December-January
4942 ! interpolation if dataset is being cycled yearly.
4943 !
4944 if (ozncyc .and. np1 == 1) then ! Dec-Jan interpolation
4945 deltat = cdayozp + daysperyear - cdayozm
4946 if (intjulian > cdayozp) then ! We are in December
4947 fact1 = (cdayozp + daysperyear - intjulian)/deltat
4948 fact2 = (intjulian - cdayozm)/deltat
4949 else ! We are in January
4950 fact1 = (cdayozp - intjulian)/deltat
4951 fact2 = (intjulian + daysperyear - cdayozm)/deltat
4952 end if
4953 else
4954 deltat = cdayozp - cdayozm
4955 fact1 = (cdayozp - intjulian)/deltat
4956 fact2 = (intjulian - cdayozm)/deltat
4957 end if
4958 !
4959 ! Time interpolation.
4960 !
4961 do j=jts,jte
4962 do k=1,levsiz
4963 do i=its,ite
4964 ozmixt(i,k,j) = ozmixm(i,k,j,nm+1)*fact1 + ozmixm(i,k,j,np+1)*fact2
4965 end do
4966 end do
4967 end do
4969 END SUBROUTINE ozn_time_int
4971 SUBROUTINE ozn_p_int(p ,pin, levsiz, ozmixt, o3vmr, &
4972 ids , ide , jds , jde , kds , kde , &
4973 ims , ime , jms , jme , kms , kme , &
4974 its , ite , jts , jte , kts , kte )
4976 !-----------------------------------------------------------------------
4977 !
4978 ! Purpose: Interpolate ozone from current time-interpolated values to model levels
4979 !
4980 ! Method: Use pressure values to determine interpolation levels
4981 !
4982 ! Author: Bruce Briegleb
4983 ! WW: Adapted for general use
4984 !
4985 !--------------------------------------------------------------------------
4986 implicit none
4987 !--------------------------------------------------------------------------
4988 !
4989 ! Arguments
4990 !
4991 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4992 ims,ime, jms,jme, kms,kme, &
4993 its,ite, jts,jte, kts,kte
4995 integer, intent(in) :: levsiz ! number of ozone layers
4997 real, intent(in) :: p(ims:ime,kms:kme,jms:jme) ! level pressures (mks, bottom-up)
4998 real, intent(in) :: pin(levsiz) ! ozone data level pressures (mks, top-down)
4999 real, intent(in) :: ozmixt(ims:ime,levsiz,jms:jme) ! ozone mixing ratio
5001 real, intent(out) :: o3vmr(ims:ime,kms:kme,jms:jme) ! ozone volume mixing ratio
5002 !
5003 ! local storage
5004 !
5005 real pmid(its:ite,kts:kte)
5006 integer i,j ! longitude index
5007 integer k, kk, kkstart, kout! level indices
5008 integer kupper(its:ite) ! Level indices for interpolation
5009 integer kount ! Counter
5010 integer ncol, pver
5012 real dpu ! upper level pressure difference
5013 real dpl ! lower level pressure difference
5015 ncol = ite - its + 1
5016 pver = kte - kts + 1
5018 do j=jts,jte
5019 !
5020 ! Initialize index array
5021 !
5022 ! do i=1, ncol
5023 do i=its, ite
5024 kupper(i) = 1
5025 end do
5026 !
5027 ! Reverse the pressure array, and pin is in Pa, the same as model pmid
5028 !
5029 do k = kts,kte
5030 kk = kte - k + kts
5031 do i = its,ite
5032 pmid(i,kk) = p(i,k,j)
5033 enddo
5034 enddo
5036 do k=1,pver
5038 kout = pver - k + 1
5039 ! kout = k
5040 !
5041 ! Top level we need to start looking is the top level for the previous k
5042 ! for all longitude points
5043 !
5044 kkstart = levsiz
5045 ! do i=1,ncol
5046 do i=its,ite
5047 kkstart = min0(kkstart,kupper(i))
5048 end do
5049 kount = 0
5050 !
5051 ! Store level indices for interpolation
5052 !
5053 do kk=kkstart,levsiz-1
5054 ! do i=1,ncol
5055 do i=its,ite
5056 if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then
5057 kupper(i) = kk
5058 kount = kount + 1
5059 end if
5060 end do
5061 !
5062 ! If all indices for this level have been found, do the interpolation and
5063 ! go to the next level
5064 !
5065 if (kount.eq.ncol) then
5066 ! do i=1,ncol
5067 do i=its,ite
5068 dpu = pmid(i,k) - pin(kupper(i))
5069 dpl = pin(kupper(i)+1) - pmid(i,k)
5070 o3vmr(i,kout,j) = (ozmixt(i,kupper(i),j)*dpl + &
5071 ozmixt(i,kupper(i)+1,j)*dpu)/(dpl + dpu)
5072 end do
5073 goto 35
5074 end if
5075 end do
5076 !
5077 ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
5078 ! must extrapolate from the bottom or top ozone data level for at least some
5079 ! of the longitude points.
5080 !
5081 ! do i=1,ncol
5082 do i=its,ite
5083 if (pmid(i,k) .lt. pin(1)) then
5084 o3vmr(i,kout,j) = ozmixt(i,1,j)*pmid(i,k)/pin(1)
5085 else if (pmid(i,k) .gt. pin(levsiz)) then
5086 o3vmr(i,kout,j) = ozmixt(i,levsiz,j)
5087 else
5088 dpu = pmid(i,k) - pin(kupper(i))
5089 dpl = pin(kupper(i)+1) - pmid(i,k)
5090 o3vmr(i,kout,j) = (ozmixt(i,kupper(i),j)*dpl + &
5091 ozmixt(i,kupper(i)+1,j)*dpu)/(dpl + dpu)
5092 end if
5093 end do
5095 if (kount.gt.ncol) then
5096 ! call endrun ('OZN_P_INT: Bad ozone data: non-monotonicity suspected')
5097 call wrf_error_fatal ('OZN_P_INT: Bad ozone data: non-monotonicity suspected')
5098 end if
5099 35 continue
5101 end do
5102 end do
5104 return
5105 END SUBROUTINE ozn_p_int
5107 SUBROUTINE aer_time_int(julday,julian,aerodm,aerodt,levsiz,num_months,no_src, &
5108 ids , ide , jds , jde , kds , kde , &
5109 ims , ime , jms , jme , kms , kme , &
5110 its , ite , jts , jte , kts , kte )
5112 ! adapted from oznint from CAM module
5113 ! input: aerodm - read from physics_init
5114 ! output: aerodt - time interpolated
5116 ! USE module_ra_cam_support, ONLY : getfactors
5118 IMPLICIT NONE
5120 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
5121 ims,ime, jms,jme, kms,kme, &
5122 its,ite, jts,jte, kts,kte
5124 INTEGER, INTENT(IN ) :: levsiz, num_months, no_src
5126 REAL, DIMENSION( ims:ime, levsiz, jms:jme, num_months, no_src ), &
5127 INTENT(IN ) :: aerodm
5129 INTEGER, INTENT(IN ) :: JULDAY
5130 REAL, INTENT(IN ) :: JULIAN
5132 REAL, DIMENSION( ims:ime, levsiz, jms:jme, no_src ), &
5133 INTENT(OUT ) :: aerodt
5135 !Local
5136 REAL :: intJULIAN
5137 integer :: np1,np,nm,m,k,i,j,s
5138 integer :: IJUL
5139 integer, dimension(12) :: date_oz
5140 data date_oz/16, 45, 75, 105, 136, 166, 197, 228, 258, 289, 319, 350/
5141 real, parameter :: daysperyear = 365. ! number of days in a year
5142 real :: cdayozp, cdayozm
5143 real :: fact1, fact2, deltat
5144 logical :: finddate
5145 logical :: ozncyc
5146 CHARACTER(LEN=256) :: msgstr
5148 ozncyc = .true.
5149 ! JULIAN starts from 0.0 at 0Z on 1 Jan.
5150 intJULIAN = JULIAN + 1.0 ! offset by one day
5151 ! jan 1st 00z is julian=1.0 here
5152 IJUL=INT(intJULIAN)
5153 ! Note that following will drift.
5154 ! Need to use actual month/day info to compute julian.
5155 intJULIAN=intJULIAN-FLOAT(IJUL)
5156 IJUL=MOD(IJUL,365)
5157 IF(IJUL.EQ.0)IJUL=365
5158 intJULIAN=intJULIAN+IJUL
5159 np1=1
5160 finddate=.false.
5162 ! do m=1,num_months
5163 do m=1,12
5164 if(date_oz(m).gt.intjulian.and..not.finddate) then
5165 np1=m
5166 finddate=.true.
5167 endif
5168 enddo
5169 cdayozp=date_oz(np1)
5171 if(np1.gt.1) then
5172 cdayozm=date_oz(np1-1)
5173 np=np1
5174 nm=np-1
5175 else
5176 cdayozm=date_oz(12)
5177 np=np1
5178 nm=12
5179 endif
5181 ! call getfactors(ozncyc,np1, cdayozm, cdayozp,intjulian, &
5182 ! fact1, fact2)
5183 !
5184 ! Determine time interpolation factors. Account for December-January
5185 ! interpolation if dataset is being cycled yearly.
5186 !
5187 if (ozncyc .and. np1 == 1) then ! Dec-Jan interpolation
5188 deltat = cdayozp + daysperyear - cdayozm
5189 if (intjulian > cdayozp) then ! We are in December
5190 fact1 = (cdayozp + daysperyear - intjulian)/deltat
5191 fact2 = (intjulian - cdayozm)/deltat
5192 else ! We are in January
5193 fact1 = (cdayozp - intjulian)/deltat
5194 fact2 = (intjulian + daysperyear - cdayozm)/deltat
5195 end if
5196 else
5197 deltat = cdayozp - cdayozm
5198 fact1 = (cdayozp - intjulian)/deltat
5199 fact2 = (intjulian - cdayozm)/deltat
5200 end if
5201 !
5202 ! Time interpolation.
5203 !
5204 do s=1, no_src
5205 do j=jts,jte
5206 do k=1,levsiz
5207 do i=its,ite
5208 aerodt(i,k,j,s) = aerodm(i,k,j,nm,s)*fact1 + aerodm(i,k,j,np,s)*fact2
5209 end do
5210 end do
5211 end do
5212 end do
5214 END SUBROUTINE aer_time_int
5216 SUBROUTINE aer_p_int(p ,pin, levsiz, aerodt, aerod, no_src, pf, totaod, &
5217 ids , ide , jds , jde , kds , kde , &
5218 ims , ime , jms , jme , kms , kme , &
5219 its , ite , jts , jte , kts , kte )
5221 !-----------------------------------------------------------------------
5222 !
5223 ! Purpose: Interpolate aerosol from current time-interpolated values to model levels
5224 !
5225 ! Method: Use pressure values to determine interpolation levels
5226 !
5227 ! Author: Bruce Briegleb
5228 ! WW: Adapted for general use
5229 !
5230 ! p: model level pressure at half levels (Pa, bottom-up)
5231 ! pf: model level pressure at full levles (Pa, bottom-up)
5232 !
5233 !--------------------------------------------------------------------------
5234 implicit none
5235 !--------------------------------------------------------------------------
5236 !
5237 ! Arguments
5238 !
5239 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
5240 ims,ime, jms,jme, kms,kme, &
5241 its,ite, jts,jte, kts,kte
5243 integer, intent(in) :: levsiz ! number of aerosol layers
5244 integer, intent(in) :: no_src ! types of aerosol
5246 real, intent(in) :: p(ims:ime,kms:kme,jms:jme)
5247 real, intent(in) :: pf(ims:ime,kms:kme,jms:jme)
5248 real, intent(in) :: pin(levsiz) ! aerosol data level pressures (mks, top-down)
5249 real, intent(in) :: aerodt(ims:ime,levsiz,jms:jme,1:no_src) ! aerosol optical depth
5251 real, intent(out) :: aerod(ims:ime,kms:kme,jms:jme,1:no_src) ! aerosol optical depth
5252 real, intent(out) :: totaod(ims:ime,jms:jme) ! total aerosol optical depth
5253 !
5254 ! local storage
5255 !
5256 real pmid(its:ite,kts:kte)
5257 integer i,j ! longitude index
5258 integer k, kk, kkstart, kout! level indices
5259 integer kupper(its:ite) ! Level indices for interpolation
5260 integer kount ! Counter
5261 integer ncol, pver, s
5263 real dpu ! upper level pressure difference
5264 real dpl ! lower level pressure difference
5265 real dpm ! pressure difference in a model layer surrounding half p
5267 ncol = ite - its + 1
5268 pver = kte - kts + 1
5270 do s=1,no_src
5271 do j=jts,jte
5272 !
5273 ! Initialize index array
5274 !
5275 do i=its, ite
5276 kupper(i) = 1
5277 end do
5278 !
5279 ! The pressure from incoming data is in hPa and top-down,
5280 ! while model pressure is in Pa and bottom-up
5281 !
5282 do k = kts,kte
5283 kk = kte - k + kts
5284 do i = its,ite
5285 pmid(i,kk) = p(i,k,j)*0.01
5286 enddo
5287 enddo
5289 do k=1,pver
5291 kout = pver - k + 1
5292 !
5293 ! Top level we need to start looking is the top level for the previous k
5294 ! for all longitude points
5295 !
5296 kkstart = levsiz
5297 do i=its,ite
5298 kkstart = min0(kkstart,kupper(i))
5299 end do
5300 kount = 0
5301 !
5302 ! Store level indices for interpolation
5303 !
5304 do kk=kkstart,levsiz-1
5305 do i=its,ite
5306 if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then
5307 kupper(i) = kk
5308 kount = kount + 1
5309 end if
5310 end do
5311 !
5312 ! If all indices for this level have been found, do the interpolation and
5313 ! go to the next level
5314 !
5315 if (kount.eq.ncol) then
5316 do i=its,ite
5317 dpu = pmid(i,k) - pin(kupper(i))
5318 dpl = pin(kupper(i)+1) - pmid(i,k)
5319 dpm = pf(i,kout,j) - pf(i,kout+1,j)
5320 aerod(i,kout,j,s) = (aerodt(i,kupper(i),j,s)*dpl + &
5321 aerodt(i,kupper(i)+1,j,s)*dpu)/(dpl + dpu)
5322 aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm
5323 end do
5324 goto 35
5325 end if
5326 end do
5327 !
5328 ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
5329 ! must extrapolate from the bottom or top aerosol data level for at least some
5330 ! of the longitude points.
5331 !
5332 do i=its,ite
5333 if (pmid(i,k) .lt. pin(1)) then
5334 dpm = pf(i,kout,j) - pf(i,kout+1,j)
5335 aerod(i,kout,j,s) = aerodt(i,1,j,s)*pmid(i,k)/pin(1)
5336 aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm
5337 else if (pmid(i,k) .gt. pin(levsiz)) then
5338 dpm = pf(i,kout,j) - pf(i,kout+1,j)
5339 aerod(i,kout,j,s) = aerodt(i,levsiz,j,s)
5340 aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm
5341 else
5342 dpu = pmid(i,k) - pin(kupper(i))
5343 dpl = pin(kupper(i)+1) - pmid(i,k)
5344 dpm = pf(i,kout,j) - pf(i,kout+1,j)
5345 aerod(i,kout,j,s) = (aerodt(i,kupper(i),j,s)*dpl + &
5346 aerodt(i,kupper(i)+1,j,s)*dpu)/(dpl + dpu)
5347 aerod(i,kout,j,s) = aerod(i,kout,j,s)*dpm
5348 end if
5349 end do
5351 if (kount.gt.ncol) then
5352 call wrf_error_fatal ('AER_P_INT: Bad aerosol data: non-monotonicity suspected')
5353 end if
5354 35 continue
5356 end do
5357 end do
5358 end do
5360 do j=jts,jte
5361 do i=its,ite
5362 totaod(i,j) = 0.
5363 end do
5364 end do
5366 do s=1,no_src
5367 do j=jts,jte
5368 do k=1,pver
5369 do i=its,ite
5370 totaod(i,j) = totaod(i,j) + aerod(i,k,j,s)
5371 end do
5372 end do
5373 end do
5374 end do
5376 return
5377 END SUBROUTINE aer_p_int
5380 !+---+-----------------------------------------------------------------+
5382 SUBROUTINE gt_aod(p_phy,DZ8W,t_phy,qvapor, nwfa,nifa,nbca, taod5503d, &
5383 & wif_input_opt, ims,ime, jms,jme, kms,kme, its,ite, jts,jte, kts,kte)
5385 USE module_mp_thompson, only: RSLF
5387 IMPLICIT NONE
5389 INTEGER, INTENT(IN):: ims,ime, jms,jme, kms,kme, &
5390 & its,ite, jts,jte, kts,kte
5392 REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: &
5393 & t_phy,p_phy, DZ8W, &
5394 & qvapor, nifa, nwfa, nbca
5395 REAL,DIMENSION(ims:ime,kms:kme,jms:jme),INTENT(INOUT):: taod5503d
5396 INTEGER, INTENT(IN):: wif_input_opt
5398 !..Local variables.
5400 REAL, DIMENSION(its:ite,kts:kte,jts:jte):: AOD_wfa, AOD_ifa, AOD_bca
5401 REAL:: RH, a_RH,b_RH, rh_d,rh_f, rhoa,qvsat, unit_bext1,unit_bext3,unit_bext4
5402 REAL:: ntemp
5403 INTEGER :: i, k, j, RH_idx, RH_idx1, RH_idx2, t_idx
5404 INTEGER, PARAMETER:: rind=8
5405 REAL, DIMENSION(rind), PARAMETER:: rh_arr = &
5406 & (/10., 60., 70., 80., 85., 90., 95., 99.8/)
5407 REAL, DIMENSION(rind,4,3) :: lookup_tabl ! RH, temp, water-friendly, ice-friendly
5409 lookup_tabl(1,1,1) = 5.73936E-15
5410 lookup_tabl(1,1,2) = 2.63577E-12
5411 lookup_tabl(1,1,3) = 7.83852E-16
5412 lookup_tabl(1,2,1) = 5.73936E-15
5413 lookup_tabl(1,2,2) = 2.63577E-12
5414 lookup_tabl(1,2,3) = 7.83803E-16
5415 lookup_tabl(1,3,1) = 5.73936E-15
5416 lookup_tabl(1,3,2) = 2.63577E-12
5417 lookup_tabl(1,3,3) = 7.83770E-16
5418 lookup_tabl(1,4,1) = 5.73936E-15
5419 lookup_tabl(1,4,2) = 2.63577E-12
5420 lookup_tabl(1,4,3) = 7.83692E-16
5422 lookup_tabl(2,1,1) = 6.93515E-15
5423 lookup_tabl(2,1,2) = 2.72095E-12
5424 lookup_tabl(2,1,3) = 9.96291E-16
5425 lookup_tabl(2,2,1) = 6.93168E-15
5426 lookup_tabl(2,2,2) = 2.72092E-12
5427 lookup_tabl(2,2,3) = 9.94950E-16
5428 lookup_tabl(2,3,1) = 6.92570E-15
5429 lookup_tabl(2,3,2) = 2.72091E-12
5430 lookup_tabl(2,3,3) = 9.93238E-16
5431 lookup_tabl(2,4,1) = 6.91833E-15
5432 lookup_tabl(2,4,2) = 2.72087E-12
5433 lookup_tabl(2,4,3) = 9.91346E-16
5435 lookup_tabl(3,1,1) = 7.24707E-15
5436 lookup_tabl(3,1,2) = 2.77219E-12
5437 lookup_tabl(3,1,3) = 1.11936E-15
5438 lookup_tabl(3,2,1) = 7.23809E-15
5439 lookup_tabl(3,2,2) = 2.77222E-12
5440 lookup_tabl(3,2,3) = 1.11687E-15
5441 lookup_tabl(3,3,1) = 7.23108E-15
5442 lookup_tabl(3,3,2) = 2.77201E-12
5443 lookup_tabl(3,3,3) = 1.11447E-15
5444 lookup_tabl(3,4,1) = 7.21800E-15
5445 lookup_tabl(3,4,2) = 2.77111E-12
5446 lookup_tabl(3,4,3) = 1.11061E-15
5448 lookup_tabl(4,1,1) = 8.95130E-15
5449 lookup_tabl(4,1,2) = 2.87263E-12
5450 lookup_tabl(4,1,3) = 1.36116E-15
5451 lookup_tabl(4,2,1) = 9.01582E-15
5452 lookup_tabl(4,2,2) = 2.87252E-12
5453 lookup_tabl(4,2,3) = 1.35479E-15
5454 lookup_tabl(4,3,1) = 9.13216E-15
5455 lookup_tabl(4,3,2) = 2.87241E-12
5456 lookup_tabl(4,3,3) = 1.34787E-15
5457 lookup_tabl(4,4,1) = 9.16219E-15
5458 lookup_tabl(4,4,2) = 2.87211E-12
5459 lookup_tabl(4,4,3) = 1.33910E-15
5461 lookup_tabl(5,1,1) = 1.06695E-14
5462 lookup_tabl(5,1,2) = 2.96752E-12
5463 lookup_tabl(5,1,3) = 1.58848E-15
5464 lookup_tabl(5,2,1) = 1.06370E-14
5465 lookup_tabl(5,2,2) = 2.96726E-12
5466 lookup_tabl(5,2,3) = 1.57854E-15
5467 lookup_tabl(5,3,1) = 1.05999E-14
5468 lookup_tabl(5,3,2) = 2.96702E-12
5469 lookup_tabl(5,3,3) = 1.56648E-15
5470 lookup_tabl(5,4,1) = 1.05443E-14
5471 lookup_tabl(5,4,2) = 2.96603E-12
5472 lookup_tabl(5,4,3) = 1.55057E-15
5474 lookup_tabl(6,1,1) = 1.37908E-14
5475 lookup_tabl(6,1,2) = 3.15081E-12
5476 lookup_tabl(6,1,3) = 2.02033E-15
5477 lookup_tabl(6,2,1) = 1.37172E-14
5478 lookup_tabl(6,2,2) = 3.15020E-12
5479 lookup_tabl(6,2,3) = 1.99948E-15
5480 lookup_tabl(6,3,1) = 1.36362E-14
5481 lookup_tabl(6,3,2) = 3.14927E-12
5482 lookup_tabl(6,3,3) = 1.97488E-15
5483 lookup_tabl(6,4,1) = 1.35287E-14
5484 lookup_tabl(6,4,2) = 3.14817E-12
5485 lookup_tabl(6,4,3) = 1.94523E-15
5487 lookup_tabl(7,1,1) = 2.26019E-14
5488 lookup_tabl(7,1,2) = 3.66798E-12
5489 lookup_tabl(7,1,3) = 3.14156E-15
5490 lookup_tabl(7,2,1) = 2.24435E-14
5491 lookup_tabl(7,2,2) = 3.66540E-12
5492 lookup_tabl(7,2,3) = 3.08114E-15
5493 lookup_tabl(7,3,1) = 2.23254E-14
5494 lookup_tabl(7,3,2) = 3.66173E-12
5495 lookup_tabl(7,3,3) = 3.01021E-15
5496 lookup_tabl(7,4,1) = 2.20496E-14
5497 lookup_tabl(7,4,2) = 3.65796E-12
5498 lookup_tabl(7,4,3) = 2.92456E-15
5500 lookup_tabl(8,1,1) = 4.41983E-13
5501 lookup_tabl(8,1,2) = 7.50091E-11
5502 lookup_tabl(8,1,3) = 1.95503E-14
5503 lookup_tabl(8,2,1) = 3.93335E-13
5504 lookup_tabl(8,2,2) = 6.79097E-11
5505 lookup_tabl(8,2,3) = 1.74308E-14
5506 lookup_tabl(8,3,1) = 3.45569E-13
5507 lookup_tabl(8,3,2) = 6.07845E-11
5508 lookup_tabl(8,3,3) = 1.53194E-14
5509 lookup_tabl(8,4,1) = 2.96971E-13
5510 lookup_tabl(8,4,2) = 5.36085E-11
5511 lookup_tabl(8,4,3) = 1.32479E-14
5513 DO j=jts,jte
5514 DO k=kts,kte
5515 DO i=its,ite
5516 AOD_wfa(i,k,j) = 0.
5517 AOD_ifa(i,k,j) = 0.
5518 if ( wif_input_opt .eq. 2 ) then
5519 AOD_bca(i,k,j) = 0.
5520 end if
5521 END DO
5522 END DO
5523 END DO
5525 DO j=jts,jte
5526 DO k=kts,kte
5527 DO i=its,ite
5528 rhoa = p_phy(i,k,j)/(287.*t_phy(i,k,j))
5529 t_idx = MAX(1, MIN(nint(10.999-0.0333*t_phy(i,k,j)),4))
5530 qvsat = rslf(p_phy(i,k,j),t_phy(i,k,j))
5531 RH = MIN(99.1, MAX(10.1, qvapor(i,k,j)/qvsat*100.))
5533 !..Get the index for the RH array element
5535 if (RH .lt. 60) then
5536 RH_idx1 = 1
5537 RH_idx2 = 2
5538 elseif (RH .ge. 60 .AND. RH.lt.80) then
5539 a_RH = 0.1
5540 b_RH = -4
5541 RH_idx = nint(a_RH*RH+b_RH)
5542 rh_d = rh-rh_arr(rh_idx)
5543 if (rh_d .lt. 0) then
5544 RH_idx1 = RH_idx-1
5545 RH_idx2 = RH_idx
5546 else
5547 RH_idx1 = RH_idx
5548 RH_idx2 = RH_idx+1
5549 if (RH_idx2.gt.rind) then
5550 RH_idx2 = rind
5551 RH_idx1 = rind-1
5552 endif
5553 endif
5554 else
5555 a_RH = 0.2
5556 b_RH = -12.
5557 RH_idx = MIN(rind, nint(a_RH*RH+b_RH))
5558 rh_d = rh-rh_arr(rh_idx)
5559 if (rh_d .lt. 0) then
5560 RH_idx1 = RH_idx-1
5561 RH_idx2 = RH_idx
5562 else
5563 RH_idx1 = RH_idx
5564 RH_idx2 = RH_idx+1
5565 if (RH_idx2.gt.rind) then
5566 RH_idx2 = rind
5567 RH_idx1 = rind-1
5568 endif
5569 endif
5570 endif
5572 !..RH fraction to be used
5574 rh_f = MAX(0., MIN(1.0, (rh/(100-rh)-rh_arr(rh_idx1) &
5575 & /(100-rh_arr(rh_idx1))) &
5576 & /(rh_arr(rh_idx2)/(100-rh_arr(rh_idx2)) &
5577 & -rh_arr(rh_idx1)/(100-rh_arr(rh_idx1))) ))
5580 unit_bext1 = lookup_tabl(RH_idx1,t_idx,1) &
5581 & + (lookup_tabl(RH_idx2,t_idx,1) &
5582 & - lookup_tabl(RH_idx1,t_idx,1))*rh_f
5583 unit_bext3 = lookup_tabl(RH_idx1,t_idx,2) &
5584 & + (lookup_tabl(RH_idx2,t_idx,2) &
5585 & - lookup_tabl(RH_idx1,t_idx,2))*rh_f
5586 if ( wif_input_opt .eq. 2 ) then
5587 unit_bext4 = lookup_tabl(RH_idx1,t_idx,3) &
5588 & + (lookup_tabl(RH_idx2,t_idx,3) &
5589 & - lookup_tabl(RH_idx1,t_idx,3))*rh_f
5590 end if
5592 ntemp = MAX(1., MIN(99999.E6, nwfa(i,k,j)))
5593 AOD_wfa(i,k,j) = unit_bext1*ntemp*dz8w(i,k,j)*rhoa
5595 ntemp = MAX(0.01, MIN(9999.E6, nifa(i,k,j)))
5596 AOD_ifa(i,k,j) = unit_bext3*ntemp*dz8w(i,k,j)*rhoa
5598 if ( wif_input_opt .eq. 2 ) then
5599 ntemp = MAX(1., MIN(9999.E9, nbca(i,k,j)))
5600 AOD_bca(i,k,j) = unit_bext4*ntemp*dz8w(i,k,j)*rhoa
5601 end if
5603 END DO
5604 END DO
5605 END DO
5607 DO j=jts,jte
5608 DO k=kts,kte
5609 DO i=its,ite
5610 if ( wif_input_opt .eq. 2 ) then
5611 taod5503d(i,k,j) = MAX(1.E-3, aod_wfa(i,k,j) + aod_ifa(i,k,j) + aod_bca(i,k,j))
5612 else
5613 taod5503d(i,k,j) = MAX(1.E-3, aod_wfa(i,k,j) + aod_ifa(i,k,j))
5614 end if
5615 END DO
5616 END DO
5617 END DO
5619 END SUBROUTINE gt_aod
5621 !+---+-----------------------------------------------------------------+
5623 subroutine Add_multi_perturb_swrad_perturbations (perts_albedo, perts_aod, perts_angstrom, &
5624 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow, pert_farms_albedo, pert_farms_aod, &
5625 pert_farms_angexp, pert_farms_aerasy, pert_farms_qv, pert_farms_qc, pert_farms_qs, &
5626 albedo, aod5502d, angexp2d, aerasy2d, qv, qc, qs, qv_tmp, qc_tmp, qs_tmp, its, ite, jts, jte, &
5627 ims, ime, jms, jme, kms, kme, kts, kte)
5629 implicit none
5631 integer, intent(in) :: its, ite, jts, jte, ims, ime, jms, jme, kms, kme, kts, kte
5632 real, intent(in) :: pert_farms_albedo, pert_farms_aod, pert_farms_angexp, pert_farms_aerasy, &
5633 pert_farms_qv, pert_farms_qc, pert_farms_qs
5634 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_albedo, perts_aod, perts_angstrom, &
5635 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow
5636 real, dimension(ims:ime, kms:kme, jms:jme), intent (inout) :: qv, qc, qs
5637 real, dimension(ims:ime, jms:jme), intent (inout) :: albedo, aod5502d, angexp2d, aerasy2d
5638 real, dimension (its:ite, kts:kte, jts:jte), intent(out) :: qv_tmp, qc_tmp, qs_tmp
5640 integer :: i, j, k
5643 k = 1
5644 do j = jts, jte
5645 do i = its, ite
5646 albedo(i, j) = min (max (ALBEDO_MIN, (1.0 + perts_albedo(i, k, j) * pert_farms_albedo) * albedo(i, j)), AERASY_MAX)
5647 angexp2d(i, j) = max (ANGEXP_MIN, (1.0 + perts_angstrom(i, k, j) * pert_farms_angexp) * angexp2d(i, j))
5648 aod5502d(i, j) = max (AOD_MIN, (1.0 + perts_aod(i, k, j) * pert_farms_aod) * aod5502d(i, j))
5649 aerasy2d(i, j) = min (max (AERASY_MIN, (1.0 + perts_assymfac(i, k, j) * pert_farms_aerasy) * aerasy2d(i, j)), AERASY_MAX)
5650 end do
5651 end do
5653 do j = jts, jte
5654 do k = kts, kte
5655 do i = its, ite
5656 qv_tmp(i, k, j) = qv(i, k, j)
5657 qv(i, k, j) = max (QVAPOR_MIN, (1.0 + perts_qvapor(i, k, j) * pert_farms_qv) * qv(i, k, j))
5658 qc_tmp(i, k, j) = qc(i, k, j)
5659 qc(i, k, j) = max (QCLOUD_MIN, (1.0 + perts_qcloud(i, k, j) * pert_farms_qc) * qc(i, k, j))
5660 qs_tmp(i, k, j) = qs(i, k, j)
5661 qs(i, k, j) = max (QSNOW_MIN, (1.0 + perts_qsnow(i, k, j) * pert_farms_qs) * qs(i, k, j))
5662 end do
5663 end do
5664 end do
5666 end subroutine Add_multi_perturb_swrad_perturbations
5668 subroutine Remove_multi_perturb_swrad_perturbations (perts_albedo, perts_aod, perts_angstrom, &
5669 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow, pert_farms_albedo, pert_farms_aod, &
5670 pert_farms_angexp, pert_farms_aerasy, pert_farms_qv, pert_farms_qc, pert_farms_qs, &
5671 albedo, aod5502d, angexp2d, aerasy2d, qv, qc, qs, qv_tmp, qc_tmp, qs_tmp, its, ite, jts, jte, &
5672 ims, ime, jms, jme, kms, kme, kts, kte)
5674 implicit none
5676 integer, intent(in) :: its, ite, jts, jte, ims, ime, jms, jme, kms, kme, kts, kte
5677 real, intent(in) :: pert_farms_albedo, pert_farms_aod, pert_farms_angexp, pert_farms_aerasy, &
5678 pert_farms_qv, pert_farms_qc, pert_farms_qs
5679 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_albedo, perts_aod, perts_angstrom, &
5680 perts_assymfac, perts_qvapor, perts_qcloud, perts_qsnow
5681 real, dimension(ims:ime, kms:kme, jms:jme), intent (inout) :: qv, qc, qs
5682 real, dimension(ims:ime, jms:jme), intent (inout) :: albedo, aod5502d, angexp2d, aerasy2d
5683 real, dimension (its:ite, kts:kte, jts:jte), intent(in) :: qv_tmp, qc_tmp, qs_tmp
5685 integer :: i, j, k
5688 k = 1
5689 do j = jts, jte
5690 do i = its, ite
5691 albedo(i, j) = min (max (ALBEDO_MIN, albedo(i, j) / (1.0 + perts_albedo(i, k, j) * pert_farms_albedo)), ALBEDO_MAX)
5692 angexp2d(i, j) = max (ANGEXP_MIN, angexp2d(i, j) / (1.0 + perts_angstrom(i, k, j) * pert_farms_angexp))
5693 aod5502d(i, j) = max (AOD_MIN, aod5502d(i, j) / (1.0 + perts_aod(i, k, j) * pert_farms_aod))
5694 aerasy2d(i, j) = min (max (AERASY_MIN, aerasy2d(i, j) / (1.0 + perts_assymfac(i, k, j) * pert_farms_aerasy)), AERASY_MAX)
5695 end do
5696 end do
5698 do j = jts, jte
5699 do k = kts, kte
5700 do i = its, ite
5701 qv(i, k, j) = max (QVAPOR_MIN, qv(i, k, j) - perts_qvapor(i, k, j) * pert_farms_qv * qv_tmp(i, k, j))
5702 qc(i, k, j) = max (QCLOUD_MIN, qc(i, k, j) - perts_qcloud(i, k, j) * pert_farms_qc * qc_tmp(i, k, j))
5703 qs(i, k, j) = max (QSNOW_MIN, qs(i, k, j) - perts_qsnow(i, k, j) * pert_farms_qs * qs_tmp(i, k, j))
5704 end do
5705 end do
5706 end do
5708 end subroutine Remove_multi_perturb_swrad_perturbations
5710 END MODULE module_radiation_driver