| blob | d4f77dfa8f18f7f79a206e306475afb7b5352512 |
1 !WRF:MEDIATION_LAYER:PHYSICS
2 !
3 MODULE module_radiation_driver
4 CONTAINS
5 !BOP
6 ! !IROUTINE: radiation_driver - interface to radiation physics options
8 ! !INTERFACE:
9 SUBROUTINE radiation_driver ( &
10 ACFRCV ,ACFRST ,ALBEDO &
11 ,CFRACH ,CFRACL ,CFRACM &
12 ,CUPPT ,CZMEAN ,DT &
13 ,DZ8W ,EMISS ,GLW &
14 ,GMT ,GSW ,HBOT &
15 ,HTOP ,HBOTR ,HTOPR &
16 ,ICLOUD &
17 ,ITIMESTEP,JULDAY, JULIAN &
18 ,JULYR ,LW_PHYSICS &
19 ,NCFRCV ,NCFRST ,NPHS &
20 ,P8W ,P ,PI &
21 ,RADT ,RA_CALL_OFFSET &
22 ,RHO ,RLWTOA &
23 ,RSWTOA ,RTHRATEN &
24 ,RTHRATENLW ,RTHRATENSW &
25 ,SNOW ,STEPRA ,SWDOWN &
26 ,SWDOWNC ,SW_PHYSICS &
27 ,T8W ,T ,TAUCLDC &
28 ,TAUCLDI ,TSK ,VEGFRA &
29 ,WARM_RAIN ,XICE ,XLAND &
30 ,XLAT ,XLONG ,YR &
31 !Optional solar variables
32 ,DECLINX ,SOLCONX ,COSZEN ,HRANG &
33 , CEN_LAT &
34 ,Z &
35 ,LEVSIZ, N_OZMIXM &
36 ,N_AEROSOLC &
37 ,PAERLEV &
38 ,CAM_ABS_DIM1, CAM_ABS_DIM2 &
39 ,CAM_ABS_FREQ_S &
40 ,XTIME &
41 ,CURR_SECS, ADAPT_STEP_FLAG &
42 ! indexes
43 ,IDS,IDE, JDS,JDE, KDS,KDE &
44 ,IMS,IME, JMS,JME, KMS,KME &
45 ,i_start,i_end &
46 ,j_start,j_end &
47 ,kts, kte &
48 ,num_tiles &
49 ! Optional
50 , TLWDN, TLWUP & ! goddard schemes
51 , SLWDN, SLWUP & ! goddard schemes
52 , TSWDN, TSWUP & ! goddard schemes
53 , SSWDN, SSWUP & ! goddard schemes
54 , CLDFRA &
55 , PB &
56 , F_ICE_PHY,F_RAIN_PHY &
57 , QV, F_QV &
58 , QC, F_QC &
59 , QR, F_QR &
60 , QI, F_QI &
61 , QS, F_QS &
62 , QG, F_QG &
63 , QNDROP, F_QNDROP &
64 ,ACSWUPT ,ACSWUPTC &
65 ,ACSWDNT ,ACSWDNTC &
66 ,ACSWUPB ,ACSWUPBC &
67 ,ACSWDNB ,ACSWDNBC &
68 ,ACLWUPT ,ACLWUPTC &
69 ,ACLWDNT ,ACLWDNTC &
70 ,ACLWUPB ,ACLWUPBC &
71 ,ACLWDNB ,ACLWDNBC &
72 ,SWUPT ,SWUPTC &
73 ,SWDNT ,SWDNTC &
74 ,SWUPB ,SWUPBC &
75 ,SWDNB ,SWDNBC &
76 ,LWUPT ,LWUPTC &
77 ,LWDNT ,LWDNTC &
78 ,LWUPB ,LWUPBC &
79 ,LWDNB ,LWDNBC &
80 ,LWCF &
81 ,SWCF &
82 ,OLR &
83 ,OZMIXM, PIN &
84 ,M_PS_1, M_PS_2, AEROSOLC_1 &
85 ,AEROSOLC_2, M_HYBI0 &
86 ,ABSTOT, ABSNXT, EMSTOT &
87 ,CU_RAD_FEEDBACK &
88 ,AER_RA_FEEDBACK &
89 ,QC_ADJUST , QI_ADJUST &
90 ,PM2_5_DRY, PM2_5_WATER &
91 ,PM2_5_DRY_EC &
92 ,TAUAER300, TAUAER400 & ! jcb
93 ,TAUAER600, TAUAER999 & ! jcb
94 ,GAER300, GAER400, GAER600, GAER999 & ! jcb
95 ,WAER300, WAER400, WAER600, WAER999 & ! jcb
96 ,TAUAERlw1, TAUAERlw2 &
97 ,TAUAERlw3, TAUAERlw4 &
98 ,TAUAERlw5, TAUAERlw6 &
99 ,TAUAERlw7, TAUAERlw8 &
100 ,TAUAERlw9, TAUAERlw10 &
101 ,TAUAERlw11, TAUAERlw12 &
102 ,TAUAERlw13, TAUAERlw14 &
103 ,TAUAERlw15, TAUAERlw16 &
104 ,progn &
105 ,slope_rad,topo_shading &
106 ,shadowmask,ht,dx,dy &
107 ,SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC & ! Optional
108 ,LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC & ! Optional
109 ,radtacttime &
110 ,ALSWVISDIR, ALSWVISDIF, ALSWNIRDIR, ALSWNIRDIF & !fds ssib alb comp (06/2010)
111 ,SWVISDIR, SWVISDIF, SWNIRDIR, SWNIRDIF & !fds ssib swr comp (06/2010)
112 ,SF_SURFACE_PHYSICS & !fds
113 )
116 !-------------------------------------------------------------------------
118 ! !USES:
119 USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME &
120 ,RRTMG_LWSCHEME, RRTMG_SWSCHEME &
121 ,SWRADSCHEME, GSFCSWSCHEME &
122 ,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME &
123 ,HELDSUAREZ &
124 #ifdef HWRF
125 ,HWRFSWSCHEME, HWRFLWSCHEME &
126 #endif
127 ,goddardlwscheme &
128 ,goddardswscheme &
129 ,FLGLWSCHEME, FLGSWSCHEME
131 USE module_model_constants
132 #ifndef HWRF
133 USE module_wrf_error , ONLY : wrf_err_message
134 #endif
136 ! *** add new modules of schemes here
138 USE module_ra_sw , ONLY : swrad
139 USE module_ra_gsfcsw , ONLY : gsfcswrad
140 USE module_ra_rrtm , ONLY : rrtmlwrad
141 USE module_ra_rrtmg_lw , ONLY : rrtmg_lwrad
142 USE module_ra_rrtmg_sw , ONLY : rrtmg_swrad
143 USE module_ra_cam , ONLY : camrad
144 USE module_ra_gfdleta , ONLY : etara
145 #ifdef HWRF
146 USE module_ra_hwrf
147 #endif
148 USE module_ra_hs , ONLY : hsrad
150 USE module_ra_goddard , ONLY : goddardrad
151 USE module_ra_flg , ONLY : RAD_FLG
155 ! This driver calls subroutines for the radiation parameterizations.
156 !
157 ! short wave radiation choices:
158 ! 1. swrad (19??)
159 ! 4. rrtmg_sw - Added November 2008, MJIacono, AER, Inc.
160 !
161 ! long wave radiation choices:
162 ! 1. rrtmlwrad
163 ! 4. rrtmg_lw - Added November 2008, MJIacono, AER, Inc.
164 !
165 !----------------------------------------------------------------------
166 IMPLICIT NONE
167 !<DESCRIPTION>
168 !
169 ! Radiation_driver is the WRF mediation layer routine that provides the interface to
170 ! to radiation physics packages in the WRF model layer. The radiation
171 ! physics packages to call are chosen by setting the namelist variable
172 ! (Rconfig entry in Registry) to the integer value assigned to the
173 ! particular package (package entry in Registry). For example, if the
174 ! namelist variable ra_lw_physics is set to 1, this corresponds to the
175 ! Registry Package entry for swradscheme. Note that the Package
176 ! names in the Registry are defined constants (frame/module_state_description.F)
177 ! in the CASE statements in this routine.
178 !
179 ! Among the arguments is moist, a four-dimensional scalar array storing
180 ! a variable number of moisture tracers, depending on the physics
181 ! configuration for the WRF run, as determined in the namelist. The
182 ! highest numbered index of active moisture tracers the integer argument
183 ! n_moist (note: the number of tracers at run time is the quantity
184 ! <tt>n_moist - PARAM_FIRST_SCALAR + 1</tt> , not n_moist. Individual tracers
185 ! may be indexed from moist by the Registry name of the tracer prepended
186 ! with P_; for example P_QC is the index of cloud water. An index
187 ! represents a valid, active field only if the index is greater than
188 ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual
189 ! indices for each tracer is defined in module_state_description and
190 ! set in <a href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a> defined in frame/module_configure.F.
191 !
192 ! Physics drivers in WRF 2.0 and higher, originally model-layer
193 ! routines, have been promoted to mediation layer routines and they
194 ! contain OpenMP threaded loops over tiles. Thus, physics drivers
195 ! are called from single-threaded regions in the solver. The physics
196 ! routines that are called from the physics drivers are model-layer
197 ! routines and fully tile-callable and thread-safe.
198 !</DESCRIPTION>
199 !
200 !======================================================================
201 ! Grid structure in physics part of WRF
202 !----------------------------------------------------------------------
203 ! The horizontal velocities used in the physics are unstaggered
204 ! relative to temperature/moisture variables. All predicted
205 ! variables are carried at half levels except w, which is at full
206 ! levels. Some arrays with names (*8w) are at w (full) levels.
207 !
208 !----------------------------------------------------------------------
209 ! In WRF, kms (smallest number) is the bottom level and kme (largest
210 ! number) is the top level. In your scheme, if 1 is at the top level,
211 ! then you have to reverse the order in the k direction.
212 !
213 ! kme - half level (no data at this level)
214 ! kme ----- full level
215 ! kme-1 - half level
216 ! kme-1 ----- full level
217 ! .
218 ! .
219 ! .
220 ! kms+2 - half level
221 ! kms+2 ----- full level
222 ! kms+1 - half level
223 ! kms+1 ----- full level
224 ! kms - half level
225 ! kms ----- full level
226 !
227 !======================================================================
228 ! Grid structure in physics part of WRF
229 !
230 !-------------------------------------
231 ! The horizontal velocities used in the physics are unstaggered
232 ! relative to temperature/moisture variables. All predicted
233 ! variables are carried at half levels except w, which is at full
234 ! levels. Some arrays with names (*8w) are at w (full) levels.
235 !
236 !==================================================================
237 ! Definitions
238 !-----------
239 ! Theta potential temperature (K)
240 ! Qv water vapor mixing ratio (kg/kg)
241 ! Qc cloud water mixing ratio (kg/kg)
242 ! Qr rain water mixing ratio (kg/kg)
243 ! Qi cloud ice mixing ratio (kg/kg)
244 ! Qs snow mixing ratio (kg/kg)
245 !-----------------------------------------------------------------
246 !-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3)
247 !-- PM2_5_WATER PM2.5 water mass (ug m^-3)
248 !-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3)
249 !-- RTHRATEN Theta tendency
250 ! due to radiation (K/s)
251 !-- RTHRATENLW Theta tendency
252 ! due to long wave radiation (K/s)
253 !-- RTHRATENSW Theta temperature tendency
254 ! due to short wave radiation (K/s)
255 !-- dt time step (s)
256 !-- itimestep number of time steps
257 !-- GLW downward long wave flux at ground surface (W/m^2)
258 !-- GSW net short wave flux at ground surface (W/m^2)
259 !-- SWDOWN downward short wave flux at ground surface (W/m^2)
260 !-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ)
261 !-- RLWTOA upward long wave at top of atmosphere (w/m2)
262 !-- RSWTOA upward short wave at top of atmosphere (w/m2)
263 !-- XLAT latitude, south is negative (degree)
264 !-- XLONG longitude, west is negative (degree)
265 !-- ALBEDO albedo (between 0 and 1)
266 !-- CLDFRA cloud fraction (between 0 and 1)
267 !-- EMISS surface emissivity (between 0 and 1)
268 !-- rho_phy density (kg/m^3)
269 !-- rr dry air density (kg/m^3)
270 !-- moist moisture array (4D - last index is species) (kg/kg)
271 !-- n_moist number of moisture species
272 !-- qndrop Cloud droplet number (#/kg)
273 !-- p8w pressure at full levels (Pa)
274 !-- p_phy pressure (Pa)
275 !-- Pb base-state pressure (Pa)
276 !-- pi_phy exner function (dimensionless)
277 !-- dz8w dz between full levels (m)
278 !-- t_phy temperature (K)
279 !-- t8w temperature at full levels (K)
280 !-- GMT Greenwich Mean Time Hour of model start (hour)
281 !-- JULDAY the initial day (Julian day)
282 !-- RADT time for calling radiation (min)
283 !-- ra_call_offset -1 (old) means usually just before output, 0 after
284 !-- DEGRAD conversion factor for
285 ! degrees to radians (pi/180.) (rad/deg)
286 !-- DPD degrees per day for earth's
287 ! orbital position (deg/day)
288 !-- R_d gas constant for dry air (J/kg/K)
289 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
290 !-- G acceleration due to gravity (m/s^2)
291 !-- rvovrd R_v divided by R_d (dimensionless)
292 !-- XTIME time since simulation start (min)
293 !-- DECLIN solar declination angle (rad)
294 !-- SOLCON solar constant (W/m^2)
295 !-- ids start index for i in domain
296 !-- ide end index for i in domain
297 !-- jds start index for j in domain
298 !-- jde end index for j in domain
299 !-- kds start index for k in domain
300 !-- kde end index for k in domain
301 !-- ims start index for i in memory
302 !-- ime end index for i in memory
303 !-- jms start index for j in memory
304 !-- jme end index for j in memory
305 !-- kms start index for k in memory
306 !-- kme end index for k in memory
307 !-- i_start start indices for i in tile
308 !-- i_end end indices for i in tile
309 !-- j_start start indices for j in tile
310 !-- j_end end indices for j in tile
311 !-- kts start index for k in tile
312 !-- kte end index for k in tile
313 !-- num_tiles number of tiles
314 !
315 !==================================================================
316 !
317 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
318 ims,ime, jms,jme, kms,kme, &
319 kts,kte, &
320 num_tiles
322 INTEGER, INTENT(IN) :: lw_physics, sw_physics
324 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
325 i_start,i_end,j_start,j_end
327 INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset
328 INTEGER, INTENT(IN ) :: levsiz, n_ozmixm
329 INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2
330 REAL, INTENT(IN ) :: cam_abs_freq_s
332 LOGICAL, INTENT(IN ) :: warm_rain
334 REAL, INTENT(IN ) :: RADT
336 REAL, DIMENSION( ims:ime, jms:jme ), &
337 INTENT(IN ) :: XLAND, &
338 XICE, &
339 TSK, &
340 VEGFRA, &
341 SNOW
342 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, &
343 INTENT(IN ) :: OZMIXM
345 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ) :: OZFLG
347 REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN
349 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2
350 REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, &
351 INTENT(IN ) :: aerosolc_1, aerosolc_2
352 REAL, DIMENSION(paerlev), OPTIONAL, &
353 INTENT(IN ) :: m_hybi0
355 REAL, DIMENSION( ims:ime, jms:jme ), &
356 INTENT(INOUT) :: HTOP, &
357 HBOT, &
358 HTOPR, &
359 HBOTR, &
360 CUPPT
362 INTEGER, INTENT(IN ) :: julyr
363 !
364 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
365 INTENT(IN ) :: dz8w, &
366 z, &
367 p8w, &
368 p, &
369 pi, &
370 t, &
371 t8w, &
372 rho
373 !
374 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
375 INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
376 gaer300,gaer400,gaer600,gaer999, & ! jcb
377 waer300,waer400,waer600,waer999, & ! jcb
378 qc_adjust, qi_adjust
380 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
381 INTENT(IN ) :: tauaerlw1,tauaerlw2,tauaerlw3,tauaerlw4, & ! czhao
382 tauaerlw5,tauaerlw6,tauaerlw7,tauaerlw8, & ! czhao
383 tauaerlw9,tauaerlw10,tauaerlw11,tauaerlw12, & ! czhao
384 tauaerlw13,tauaerlw14,tauaerlw15,tauaerlw16
386 LOGICAL, INTENT(IN) :: cu_rad_feedback
388 INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
390 !jdfcz INTEGER, OPTIONAL, INTENT(IN ) :: progn,prescribe
391 INTEGER, OPTIONAL, INTENT(IN ) :: progn
392 !
393 ! variables for aerosols (only if running with chemistry)
394 !
395 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
396 INTENT(IN ) :: pm2_5_dry, &
397 pm2_5_water, &
398 pm2_5_dry_ec
399 !
400 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
401 INTENT(INOUT) :: RTHRATEN, &
402 RTHRATENLW, &
403 RTHRATENSW
405 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
406 ! INTENT(INOUT) :: SWUP, &
407 ! SWDN, &
408 ! SWUPCLEAR, &
409 ! SWDNCLEAR, &
410 ! LWUP, &
411 ! LWDN, &
412 ! LWUPCLEAR, &
413 ! LWDNCLEAR
415 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
416 ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
417 ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
418 ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
419 ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
421 ! TOA and surface, upward and downward, total and clear fluxes
422 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
423 SWUPT, SWUPTC, SWDNT, SWDNTC, &
424 SWUPB, SWUPBC, SWDNB, SWDNBC, &
425 LWUPT, LWUPTC, LWDNT, LWDNTC, &
426 LWUPB, LWUPBC, LWDNB, LWDNBC
428 ! Upward and downward, total and clear sky layer fluxes (W m-2)
429 REAL, DIMENSION( ims:ime, kms:kme+2, jms:jme ), &
430 OPTIONAL, INTENT(INOUT) :: &
431 SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC, &
432 LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC
434 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
435 INTENT(INOUT) :: SWCF, &
436 LWCF, &
437 OLR
438 ! ---- fds (06/2010) ssib alb components ------------
439 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
440 INTENT(IN ) :: ALSWVISDIR, &
441 ALSWVISDIF, &
442 ALSWNIRDIR, &
443 ALSWNIRDIF
444 ! ---- fds (06/2010) ssib swr components ------------
445 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
446 INTENT(OUT ) :: SWVISDIR, &
447 SWVISDIF, &
448 SWNIRDIR, &
449 SWNIRDIF
450 INTEGER, OPTIONAL, INTENT(IN ) :: sf_surface_physics
451 !
452 REAL, DIMENSION( ims:ime, jms:jme ), &
453 INTENT(IN ) :: XLAT, &
454 XLONG, &
455 ALBEDO, &
456 EMISS
457 !
458 REAL, DIMENSION( ims:ime, jms:jme ), &
459 INTENT(INOUT) :: GSW, &
460 GLW
462 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN
463 !
464 REAL, INTENT(IN ) :: GMT,dt, &
465 julian, xtime
466 INTEGER, INTENT(IN ),OPTIONAL :: YR
467 !
468 INTEGER, INTENT(IN ) :: JULDAY, itimestep
469 REAL, INTENT(IN ),OPTIONAL :: CURR_SECS
470 LOGICAL, INTENT(IN ),OPTIONAL :: ADAPT_STEP_FLAG
472 INTEGER,INTENT(IN) :: NPHS
473 REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: &
474 CFRACH, & !Added
475 CFRACL, & !Added
476 CFRACM, & !Added
477 CZMEAN !Added
478 REAL, DIMENSION( ims:ime, jms:jme ), &
479 INTENT(INOUT) :: &
480 RLWTOA, & !Added
481 RSWTOA, & !Added
482 ACFRST, & !Added
483 ACFRCV !Added
485 INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: &
486 NCFRST, & !Added
487 NCFRCV !Added
489 ! Optional, only for Goddard LW and SW
490 REAL, DIMENSION(IMS:IME, JMS:JME, 1:8) :: ERBE_out !extra output for SDSU
491 REAL, DIMENSION(IMS:IME, JMS:JME), OPTIONAL, INTENT(OUT) :: &
492 TLWDN, TLWUP, &
493 SLWDN, SLWUP, &
494 TSWDN, TSWUP, &
495 SSWDN, SSWUP ! for Goddard schemes
497 ! Optional (only used by CAM lw scheme)
499 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,&
500 INTENT(INOUT) :: abstot
501 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,&
502 INTENT(INOUT) :: absnxt
503 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
504 INTENT(INOUT) :: emstot
506 !
507 ! Optional
508 !
509 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
510 OPTIONAL, &
511 INTENT(INOUT) :: CLDFRA
513 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
514 OPTIONAL, &
515 INTENT(IN ) :: &
516 F_ICE_PHY, &
517 F_RAIN_PHY
519 REAL, DIMENSION( ims:ime, jms:jme ), &
520 OPTIONAL, &
521 INTENT(OUT) :: SWDOWNC
522 !
523 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
524 OPTIONAL, &
525 INTENT(INOUT ) :: &
526 pb &
527 ,qv,qc,qr,qi,qs,qg,qndrop
529 LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop
530 !
531 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
532 OPTIONAL, &
533 INTENT(INOUT) :: taucldi,taucldc
535 ! Variables for slope-dependent radiation
537 REAL, OPTIONAL, INTENT(IN) :: dx,dy
538 INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,topo_shading
539 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: ht
540 INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: shadowmask
542 REAL , OPTIONAL, INTENT(INOUT) :: radtacttime ! Storing the time in s when radiation is called next
543 ! LOCAL VAR
545 REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON
546 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS
547 REAL, DIMENSION( ims:ime, jms:jme ) :: coszr
549 REAL :: DECLIN,SOLCON,XXLAT,TLOCTM,XT24, CEN_LAT
550 INTEGER :: i,j,k,its,ite,jts,jte,ij
551 INTEGER :: STEPABS
552 LOGICAL :: gfdl_lw,gfdl_sw
553 LOGICAL :: doabsems
554 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
556 REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, &
557 DJUL,RJUL,ECCFAC
558 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_temp,qc_temp
559 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save
561 REAL :: next_rad_time
562 LOGICAL :: run_param , doing_adapt_dt , decided
563 LOGICAL :: flg_lw, flg_sw
564 !------------------------------------------------------------------
565 ! solar related variables are added to declaration
566 !-------------------------------------------------
567 REAL, OPTIONAL, INTENT(OUT) :: DECLINX,SOLCONX
568 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZEN
569 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: HRANG
570 !------------------------------------------------------------------
571 #ifdef HWRF
572 CHARACTER(len=265) :: wrf_err_message
573 #endif
575 CALL wrf_debug (1, 'Top of Radiation Driver')
576 ! WRITE ( wrf_err_message , * ) 'itimestep = ',itimestep,', dt = ',dt,', lw and sw options = ',lw_physics,sw_physics
577 ! CALL wrf_debug (1, wrf_err_message )
578 if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return
580 ! ra_call_offset = -1 gives old method where radiation may be called just before output
581 ! ra_call_offset = 0 gives new method where radiation may be called just after output
582 ! and is also consistent with removal of offset in new XTIME
583 ! also need to account for stepra=1 which always has zero modulo output
585 doing_adapt_dt = .FALSE.
586 IF ( PRESENT(adapt_step_flag) ) THEN
587 IF ( adapt_step_flag ) THEN
588 doing_adapt_dt = .TRUE.
589 IF ( radtacttime .eq. 0. ) THEN
590 radtacttime = CURR_SECS + radt*60.
591 END IF
592 END IF
593 END IF
595 ! Do we run through this scheme or not?
597 ! Test 1: If this is the initial model time, then yes.
598 ! ITIMESTEP=1
599 ! Test 2: If the user asked for the radiation to be run every time step, then yes.
600 ! RADT=0 or STEPRA=1
601 ! Test 3: If not adaptive dt, and this is on the requested radiation frequency, then yes.
602 ! MOD(ITIMESTEP,STEPRA)=0 (or 1, depending on the offset setting)
603 ! Test 4: If using adaptive dt and the current time is past the last requested activate radiation time, then yes.
604 ! CURR_SECS >= RADTACTTIME
606 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
607 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
608 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
610 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
611 ! radiation run.
613 run_param = .FALSE.
614 decided = .FALSE.
615 IF ( ( .NOT. decided ) .AND. &
616 ( itimestep .EQ. 1 ) ) THEN
617 run_param = .TRUE.
618 decided = .TRUE.
619 END IF
621 IF ( ( .NOT. decided ) .AND. &
622 ( ( radt .EQ. 0. ) .OR. ( stepra .EQ. 1 ) ) ) THEN
623 run_param = .TRUE.
624 decided = .TRUE.
625 END IF
627 IF ( ( .NOT. decided ) .AND. &
628 ( .NOT. doing_adapt_dt ) .AND. &
629 ( MOD(itimestep,stepra) .EQ. 1+ra_call_offset ) ) THEN
630 run_param = .TRUE.
631 decided = .TRUE.
632 END IF
634 IF ( ( .NOT. decided ) .AND. &
635 ( doing_adapt_dt ) .AND. &
636 ( curr_secs .GE. radtacttime ) ) THEN
637 run_param = .TRUE.
638 decided = .TRUE.
639 radtacttime = curr_secs + radt*60
640 END IF
642 Radiation_step: IF ( run_param ) then
644 ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
645 STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA
646 IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset &
647 .or. STEPABS .eq. 1 ) THEN
648 doabsems = .true.
649 ELSE
650 doabsems = .false.
651 ENDIF
652 IF (PRESENT(adapt_step_flag)) THEN
653 IF ((adapt_step_flag)) THEN
654 IF ( (itimestep .EQ. 1) .OR. (cam_abs_freq_s .EQ. 0) .OR. &
655 ( CURR_SECS + dt >= ( INT( CURR_SECS / ( cam_abs_freq_s ) + 1 ) * cam_abs_freq_s) ) ) THEN
656 doabsems = .true.
657 ELSE
658 doabsems = .false.
659 ENDIF
660 ENDIF
661 ENDIF
663 gfdl_lw = .false.
664 gfdl_sw = .false.
667 ! moved up and out of OMP loop because it only needs to be computed once
668 ! and because it is not entirely thread-safe (XT24, TOLOCTM and XXLAT need
669 ! their thread-privacy) JM 20100217
670 DO ij = 1 , num_tiles
671 its = i_start(ij)
672 ite = i_end(ij)
673 jts = j_start(ij)
674 jte = j_end(ij)
675 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
676 DEGRAD,DPD )
678 IF(PRESENT(declinx).AND.PRESENT(solconx))THEN
679 ! saved to state arrays used in surface driver
680 declinx=declin
681 solconx=solcon
682 ENDIF
684 IF(PRESENT(coszen).AND.PRESENT(hrang))THEN
685 ! state arrays of hrang and coszen used in surface driver
686 XT24=MOD(XTIME+RADT*0.5,1440.)
687 DO j=jts,jte
688 DO i=its,ite
689 TLOCTM=GMT+XT24/60.+XLONG(I,J)/15.
690 HRANG(I,J)=15.*(TLOCTM-12.)*DEGRAD
691 XXLAT=XLAT(I,J)*DEGRAD
692 COSZEN(I,J)=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG(I,J))
693 ENDDO
694 ENDDO
695 ENDIF
696 ENDDO
698 !---------------
699 !$OMP PARALLEL DO &
700 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
702 DO ij = 1 , num_tiles
703 its = i_start(ij)
704 ite = i_end(ij)
705 jts = j_start(ij)
706 jte = j_end(ij)
708 ! initialize data
710 DO j=jts,jte
711 DO i=its,ite
712 GSW(I,J)=0.
713 GLW(I,J)=0.
714 SWDOWN(I,J)=0.
715 GLAT(I,J)=XLAT(I,J)*DEGRAD
716 GLON(I,J)=XLONG(I,J)*DEGRAD
717 ENDDO
718 ENDDO
720 DO j=jts,jte
721 DO k=kts,kte+1
722 DO i=its,ite
723 RTHRATEN(I,K,J)=0.
724 RTHRATENLW(I,K,J)=0.
725 RTHRATENSW(I,K,J)=0.
726 ! SWUP(I,K,J) = 0.0
727 ! SWDN(I,K,J) = 0.0
728 ! SWUPCLEAR(I,K,J) = 0.0
729 ! SWDNCLEAR(I,K,J) = 0.0
730 ! LWUP(I,K,J) = 0.0
731 ! LWDN(I,K,J) = 0.0
732 ! LWUPCLEAR(I,K,J) = 0.0
733 ! LWDNCLEAR(I,K,J) = 0.0
734 CEMISS(I,K,J)=0.0
735 ENDDO
736 ENDDO
737 ENDDO
739 IF ( PRESENT( SWUPFLX ) ) THEN
740 DO j=jts,jte
741 DO k=kts,kte+2
742 DO i=its,ite
743 SWUPFLX(I,K,J) = 0.0
744 SWDNFLX(I,K,J) = 0.0
745 SWUPFLXC(I,K,J) = 0.0
746 SWDNFLXC(I,K,J) = 0.0
747 LWUPFLX(I,K,J) = 0.0
748 LWDNFLX(I,K,J) = 0.0
749 LWUPFLXC(I,K,J) = 0.0
750 LWDNFLXC(I,K,J) = 0.0
751 ENDDO
752 ENDDO
753 ENDDO
754 ENDIF
756 ! temporarily modify hydrometeors (currently only done for GD scheme and WRF-Chem)
757 !
758 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
759 DO j=jts,jte
760 DO k=kts,kte
761 DO i=its,ite
762 qc_save(i,k,j) = qc(i,k,j)
763 qc(i,k,j) = qc(i,k,j) + qc_adjust(i,k,j)
764 ENDDO
765 ENDDO
766 ENDDO
767 ENDIF
768 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
769 DO j=jts,jte
770 DO k=kts,kte
771 DO i=its,ite
772 qi_save(i,k,j) = qi(i,k,j)
773 qi(i,k,j) = qi(i,k,j) + qi_adjust(i,k,j)
774 ENDDO
775 ENDDO
776 ENDDO
777 ENDIF
779 ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
780 if(PRESENT(qc) .and. PRESENT(F_QC)) then
781 DO j=jts,jte
782 DO k=kts,kte
783 DO i=its,ite
784 qc_temp(I,K,J)=qc(I,K,J)
785 ENDDO
786 ENDDO
787 ENDDO
788 else
789 DO j=jts,jte
790 DO k=kts,kte
791 DO i=its,ite
792 qc_temp(I,K,J)=0.
793 ENDDO
794 ENDDO
795 ENDDO
796 endif
797 if(PRESENT(qr) .and. PRESENT(F_QR)) then
798 DO j=jts,jte
799 DO k=kts,kte
800 DO i=its,ite
801 qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J)
802 ENDDO
803 ENDDO
804 ENDDO
805 endif
807 !---------------
808 ! Calculate constant for short wave radiation
810 lwrad_cldfra_select: SELECT CASE(lw_physics)
812 CASE (GFDLLWSCHEME)
814 !-- Do nothing, since cloud fractions (with partial cloudiness effects)
815 !-- are defined in GFDL LW/SW schemes and do not need to be initialized.
817 CASE (CAMLWSCHEME)
819 IF ( PRESENT ( CLDFRA ) .AND. &
820 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
821 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
823 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
824 F_QV,F_QC,F_QI,F_QS,t,p, &
825 F_ICE_PHY,F_RAIN_PHY, &
826 ids,ide, jds,jde, kds,kde, &
827 ims,ime, jms,jme, kms,kme, &
828 its,ite, jts,jte, kts,kte )
830 ENDIF
832 CASE (RRTMG_LWSCHEME)
834 IF ( PRESENT ( CLDFRA ) .AND. &
835 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
836 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
838 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
839 F_QV,F_QC,F_QI,F_QS,t,p, &
840 F_ICE_PHY,F_RAIN_PHY, &
841 ids,ide, jds,jde, kds,kde, &
842 ims,ime, jms,jme, kms,kme, &
843 its,ite, jts,jte, kts,kte )
845 ENDIF
847 CASE DEFAULT
849 IF ( PRESENT ( CLDFRA ) .AND. &
850 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
851 CALL cal_cldfra(CLDFRA,qc,qi,F_QC,F_QI, &
852 ids,ide, jds,jde, kds,kde, &
853 ims,ime, jms,jme, kms,kme, &
854 its,ite, jts,jte, kts,kte )
855 ENDIF
857 END SELECT lwrad_cldfra_select
859 lwrad_select: SELECT CASE(lw_physics)
862 CASE (RRTMSCHEME)
863 CALL wrf_debug (100, 'CALL rrtm')
865 CALL RRTMLWRAD( &
866 RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS &
867 ,QV3D=QV &
868 ,QC3D=QC &
869 ,QR3D=QR &
870 ,QI3D=QI &
871 ,QS3D=QS &
872 ,QG3D=QG &
873 ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t &
874 ,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G &
875 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
876 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
877 ,ICLOUD=icloud,WARM_RAIN=warm_rain &
878 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
879 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
880 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
881 )
883 CASE (goddardlwscheme)
885 CALL wrf_debug (100, 'CALL goddard longwave radiation scheme ')
886 IF (itimestep.eq.1) then
887 call wrf_message('running goddard lw radiation')
888 ENDIF
889 CALL goddardrad(sw_or_lw='lw' &
890 ,rthraten=rthraten,gsf=glw,xlat=xlat,xlong=xlong &
891 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
892 ,dz8w=dz8w,rho_phy=rho,emiss=emiss &
893 ,cldfra3d=cldfra &
894 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
895 ,julday=julday,xtime=xtime &
896 ,declin=declin,solcon=solcon &
897 , center_lat = cen_lat &
898 ,radfrq=radt,degrad=degrad &
899 ,taucldi=taucldi,taucldc=taucldc &
900 ,warm_rain=warm_rain &
901 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
902 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
903 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
904 ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban
905 ,qv3d=qv &
906 ,qc3d=qc &
907 ,qr3d=qr &
908 ,qi3d=qi &
909 ,qs3d=qs &
910 ,qg3d=qg &
911 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr &
912 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg &
913 ,erbe_out=erbe_out & !optional
914 )
916 CASE (GFDLLWSCHEME)
918 CALL wrf_debug (100, 'CALL gfdllw')
920 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
921 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
922 PRESENT(qv) .AND. PRESENT(qc) ) THEN
923 IF ( F_QV .AND. F_QC .AND. F_QS) THEN
924 gfdl_lw = .true.
925 CALL ETARA( &
926 DT=dt,XLAND=xland &
927 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
928 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
929 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
930 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
931 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
932 ,HBOTR=hbotr, HTOPR=htopr &
933 ,ALBEDO=albedo,CUPPT=cuppt &
934 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
935 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
936 ,XTIME=xtime,JULIAN=julian &
937 ,JULYR=julyr,JULDAY=julday &
938 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
939 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
940 ,ACFRST=acfrst,NCFRST=ncfrst &
941 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
942 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
943 ,THRATEN=rthraten,THRATENLW=rthratenlw &
944 ,THRATENSW=rthratensw &
945 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
946 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
947 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
948 )
949 ELSE
950 CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
951 ENDIF
952 ELSE
953 CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
954 ENDIF
956 #ifdef HWRF
957 CASE (HWRFLWSCHEME)
959 CALL wrf_debug (100, 'CALL hwrflw')
961 gfdl_lw = .true.
963 CALL HWRFRA(DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, &
964 XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, &
965 QV=qv,QW=qc_temp,QI=Qi, &
966 TSK2D=tsk,GLW=GLW,GSW=GSW, &
967 TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added
968 GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt,&
969 VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified
970 NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified
971 julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, &
972 CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added
973 ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added
974 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
975 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
976 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
979 #endif
981 CASE (CAMLWSCHEME)
983 CALL wrf_debug(100, 'CALL camrad lw')
985 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
986 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
987 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
988 .AND. PRESENT(AEROSOLC_2).AND. PRESENT(ALSWVISDIR) ) THEN
989 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
990 dolw=.true.,dosw=.false., &
991 SWUPT=SWUPT,SWUPTC=SWUPTC, &
992 SWDNT=SWDNT,SWDNTC=SWDNTC, &
993 LWUPT=LWUPT,LWUPTC=LWUPTC, &
994 LWDNT=LWDNT,LWDNTC=LWDNTC, &
995 SWUPB=SWUPB,SWUPBC=SWUPBC, &
996 SWDNB=SWDNB,SWDNBC=SWDNBC, &
997 LWUPB=LWUPB,LWUPBC=LWUPBC, &
998 LWDNB=LWDNB,LWDNBC=LWDNBC, &
999 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
1000 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
1001 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
1002 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
1003 ,QV3D=qv &
1004 ,QC3D=qc &
1005 ,QR3D=qr &
1006 ,QI3D=qi &
1007 ,QS3D=qs &
1008 ,QG3D=qg &
1009 ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010)
1010 ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010)
1011 ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010)
1012 ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010)
1013 ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds
1014 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1015 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
1016 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
1017 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
1018 dz8w=dz8w, &
1019 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1020 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
1021 num_months=n_ozmixm, &
1022 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
1023 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
1024 paerlev=paerlev, naer_c=n_aerosolc, &
1025 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
1026 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
1027 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
1028 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
1029 ,doabsems=doabsems &
1030 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1031 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1032 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1033 )
1034 ELSE
1035 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
1036 ENDIF
1038 CASE (RRTMG_LWSCHEME)
1039 CALL wrf_debug (100, 'CALL rrtmg_lw')
1041 CALL RRTMG_LWRAD( &
1042 RTHRATENLW=RTHRATEN, &
1043 LWUPT=LWUPT,LWUPTC=LWUPTC, &
1044 LWDNT=LWDNT,LWDNTC=LWDNTC, &
1045 LWUPB=LWUPB,LWUPBC=LWUPBC, &
1046 LWDNB=LWDNB,LWDNBC=LWDNBC, &
1047 GLW=GLW,OLR=RLWTOA,LWCF=LWCF, &
1048 EMISS=EMISS, &
1049 P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, &
1050 T8W=t8w,RHO3D=rho,R=R_d,G=G, &
1051 ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,&
1052 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
1053 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1054 QV3D=QV,QC3D=QC,QR3D=QR, &
1055 QI3D=QI,QS3D=QS,QG3D=QG, &
1056 F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, &
1057 F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, &
1058 #ifdef WRF_CHEM
1059 TAUAERLW1=tauaerlw1,TAUAERLW2=tauaerlw2, & ! jcb
1060 TAUAERLW3=tauaerlw3,TAUAERLW4=tauaerlw4, & ! jcb
1061 TAUAERLW5=tauaerlw5,TAUAERLW6=tauaerlw6, & ! jcb
1062 TAUAERLW7=tauaerlw7,TAUAERLW8=tauaerlw8, & ! jcb
1063 TAUAERLW9=tauaerlw9,TAUAERLW10=tauaerlw10, & ! jcb
1064 TAUAERLW11=tauaerlw11,TAUAERLW12=tauaerlw12, & ! jcb
1065 TAUAERLW13=tauaerlw13,TAUAERLW14=tauaerlw14, & ! jcb
1066 TAUAERLW15=tauaerlw15,TAUAERLW16=tauaerlw16, & ! jcb
1067 aer_ra_feedback=aer_ra_feedback, &
1068 !jdfcz progn=progn,prescribe=prescribe, &
1069 progn=progn, &
1070 #endif
1071 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
1072 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
1073 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
1074 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
1075 LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, &
1076 LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC &
1077 )
1079 CASE (HELDSUAREZ)
1080 CALL wrf_debug (100, 'CALL heldsuarez')
1082 CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, &
1083 t8w, rho, R_d,G,CP, dt, xlat, degrad, &
1084 ids,ide, jds,jde, kds,kde, &
1085 ims,ime, jms,jme, kms,kme, &
1086 its,ite, jts,jte, kts,kte )
1088 ! -- add by Jin Kong 10/2011
1089 CASE (FLGLWSCHEME)
1090 CALL wrf_debug (100, 'CALL Fu-Liou-Gu')
1091 flg_lw = .true.
1092 !test Jin Kong 11/01/2011 for ozone
1093 ozflg = 0.
1094 !test over
1095 CALL RAD_FLG( &
1096 peven=p8w,podd=p,t8w=t8w,degrees=t &
1097 ,pi3d=pi,o3=ozflg &
1098 ,G=G,CP=CP &
1099 ,albedo=ALBEDO,tskin=tsk &
1100 ,h2o=QV,cld_iccld=QI,cld_wlcld=QC &
1101 ,cld_prwc=QR,cld_pgwc=QG,cld_snow=QS &
1102 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
1103 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
1104 ,warm_rain=warm_rain &
1105 ,cloudstrf=CLDFRA &
1106 ,emiss=EMISS &
1107 ,air_den=rho &
1108 ,dz3d=dz8w &
1109 ,SOLCON=SOLCON &
1110 ,declin=DECLIN &
1111 ,xtime=xtime, xlong=xlong, xlat=xlat &
1112 ,JULDAY=JULDAY &
1113 ,gmt=gmt, radt=radt, degrad=degrad &
1114 ,dtcolumn=dt &
1115 ,ids=ids,ide=ide,jds=jds,jde=jde &
1116 ,kds=kds,kde=kde &
1117 ,ims=ims,idim=ime,jms=jms,jdim=jme &
1118 ,kms=kms,kmax=kme &
1119 ,its=its,ite=ite,jts=jts,jte=jte &
1120 ,kts=kts,kte=kte &
1121 ,uswtop=RSWTOA,ulwtop=RLWTOA &
1122 ,NETSWBOT=GSW,DLWBOT=GLW &
1123 ,DSWBOT=SWDOWN,deltat=RTHRATEN &
1124 ,dtshort=RTHRATENSW,dtlongwv=RTHRATENLW &
1125 )
1127 CALL wrf_debug(100, 'a4 Fu_Liou-Gu')
1128 ! -- end
1131 CASE DEFAULT
1133 WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
1134 CALL wrf_error_fatal ( wrf_err_message )
1136 END SELECT lwrad_select
1138 IF (lw_physics .gt. 0 .and. .not.gfdl_lw) THEN
1139 DO j=jts,jte
1140 DO k=kts,kte
1141 DO i=its,ite
1142 RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
1143 ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY)
1144 IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
1145 ENDDO
1146 ENDDO
1147 ENDDO
1148 ENDIF
1150 IF (lw_physics .eq. goddardlwscheme) THEN
1151 IF ( PRESENT (tlwdn) ) THEN
1152 DO j=jts,jte
1153 DO i=its,ite
1154 tlwdn(i,j) = erbe_out(i,j,1) ! TOA LW downwelling flux [W/m2]
1155 tlwup(i,j) = erbe_out(i,j,2) ! TOA LW upwelling flux [W/m2]
1156 slwdn(i,j) = erbe_out(i,j,3) ! surface LW downwelling flux [W/m2]
1157 slwup(i,j) = erbe_out(i,j,4) ! surface LW upwelling flux [W/m2]
1158 ENDDO
1159 ENDDO
1160 ENDIF
1161 ENDIF
1163 !
1164 swrad_cldfra_select: SELECT CASE(sw_physics)
1166 CASE (CAMSWSCHEME)
1168 IF ( PRESENT ( CLDFRA ) .AND. &
1169 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
1170 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
1172 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
1173 F_QV,F_QC,F_QI,F_QS,t,p, &
1174 F_ICE_PHY,F_RAIN_PHY, &
1175 ids,ide, jds,jde, kds,kde, &
1176 ims,ime, jms,jme, kms,kme, &
1177 its,ite, jts,jte, kts,kte )
1178 ENDIF
1180 CASE (RRTMG_SWSCHEME)
1182 IF ( PRESENT ( CLDFRA ) .AND. &
1183 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
1184 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
1186 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
1187 F_QV,F_QC,F_QI,F_QS,t,p, &
1188 F_ICE_PHY,F_RAIN_PHY, &
1189 ids,ide, jds,jde, kds,kde, &
1190 ims,ime, jms,jme, kms,kme, &
1191 its,ite, jts,jte, kts,kte )
1192 ENDIF
1194 CASE DEFAULT
1196 END SELECT swrad_cldfra_select
1198 swrad_select: SELECT CASE(sw_physics)
1200 CASE (SWRADSCHEME)
1201 CALL wrf_debug(100, 'CALL swrad')
1202 CALL SWRAD( &
1203 DT=dt,RTHRATEN=rthraten,GSW=gsw &
1204 ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
1205 #ifdef WRF_CHEM
1206 ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
1207 ,PM2_5_DRY_EC=pm2_5_dry_ec &
1208 #endif
1209 ,RHO_PHY=rho,T3D=t &
1210 ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
1211 ,R=r_d,CP=cp,G=g,JULDAY=julday &
1212 ,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
1213 ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
1214 ,warm_rain=warm_rain &
1215 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1216 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1217 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1218 ,QV3D=qv &
1219 ,QC3D=qc &
1220 ,QR3D=qr &
1221 ,QI3D=qi &
1222 ,QS3D=qs &
1223 ,QG3D=qg &
1224 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1225 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg )
1227 CASE (GSFCSWSCHEME)
1228 CALL wrf_debug(100, 'CALL gsfcswrad')
1229 CALL GSFCSWRAD( &
1230 RTHRATEN=rthraten,GSW=gsw,XLAT=xlat,XLONG=xlong &
1231 ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
1232 ,DZ8W=dz8w,RHO_PHY=rho &
1233 ,CLDFRA3D=cldfra,RSWTOA=rswtoa &
1234 ,GMT=gmt,CP=cp,G=g &
1235 ,JULDAY=julday,XTIME=xtime &
1236 ,DECLIN=declin,SOLCON=solcon &
1237 ,RADFRQ=radt,DEGRAD=degrad &
1238 ,TAUCLDI=taucldi,TAUCLDC=taucldc &
1239 ,WARM_RAIN=warm_rain &
1240 #ifdef WRF_CHEM
1241 ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
1242 ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
1243 ,GAER300=gaer300,GAER400=gaer400 & ! jcb
1244 ,GAER600=gaer600,GAER999=gaer999 & ! jcb
1245 ,WAER300=waer300,WAER400=waer400 & ! jcb
1246 ,WAER600=waer600,WAER999=waer999 & ! jcb
1247 ,aer_ra_feedback=aer_ra_feedback &
1248 #endif
1249 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1250 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1251 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1252 ,QV3D=qv &
1253 ,QC3D=qc &
1254 ,QR3D=qr &
1255 ,QI3D=qi &
1256 ,QS3D=qs &
1257 ,QG3D=qg &
1258 ,QNDROP3D=qndrop &
1259 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1260 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
1261 ,F_QNDROP=f_qndrop &
1262 )
1264 CASE (goddardswscheme)
1266 CALL wrf_debug(100, 'CALL goddard shortwave radiation scheme ')
1267 IF (itimestep.eq.1) then
1268 call wrf_message('running goddard sw radiation')
1269 ENDIF
1270 CALL goddardrad(sw_or_lw='sw' &
1271 ,rthraten=rthraten,gsf=gsw,xlat=xlat,xlong=xlong &
1272 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
1273 ,dz8w=dz8w,rho_phy=rho,emiss=emiss &
1274 ,cldfra3d=cldfra &
1275 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
1276 ,julday=julday,xtime=xtime &
1277 ,declin=declin,solcon=solcon &
1278 , center_lat = cen_lat &
1279 ,radfrq=radt,degrad=degrad &
1280 ,taucldi=taucldi,taucldc=taucldc &
1281 ,warm_rain=warm_rain &
1282 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
1283 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
1284 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
1285 ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban
1286 ,qv3d=qv &
1287 ,qc3d=qc &
1288 ,qr3d=qr &
1289 ,qi3d=qi &
1290 ,qs3d=qs &
1291 ,qg3d=qg &
1292 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr &
1293 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg &
1294 ,erbe_out=erbe_out & !optional
1295 )
1297 CASE (CAMSWSCHEME)
1298 CALL wrf_debug(100, 'CALL camrad sw')
1299 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
1300 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
1301 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
1302 .AND. PRESENT(AEROSOLC_2) .AND. PRESENT(ALSWVISDIR)) THEN
1303 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
1304 dolw=.false.,dosw=.true., &
1305 SWUPT=SWUPT,SWUPTC=SWUPTC, &
1306 SWDNT=SWDNT,SWDNTC=SWDNTC, &
1307 LWUPT=LWUPT,LWUPTC=LWUPTC, &
1308 LWDNT=LWDNT,LWDNTC=LWDNTC, &
1309 SWUPB=SWUPB,SWUPBC=SWUPBC, &
1310 SWDNB=SWDNB,SWDNBC=SWDNBC, &
1311 LWUPB=LWUPB,LWUPBC=LWUPBC, &
1312 LWDNB=LWDNB,LWDNBC=LWDNBC, &
1313 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
1314 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
1315 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
1316 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
1317 ,QV3D=qv &
1318 ,QC3D=qc &
1319 ,QR3D=qr &
1320 ,QI3D=qi &
1321 ,QS3D=qs &
1322 ,QG3D=qg &
1323 ,ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif & !fds ssib alb comp (06/2010)
1324 ,ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif & !fds ssib alb comp (06/2010)
1325 ,SWVISDIR=swvisdir ,SWVISDIF=swvisdif & !fds ssib swr comp (06/2010)
1326 ,SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif & !fds ssib swr comp (06/2010)
1327 ,SF_SURFACE_PHYSICS=sf_surface_physics & !fds
1328 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1329 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
1330 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
1331 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
1332 dz8w=dz8w, &
1333 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1334 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
1335 num_months=n_ozmixm, &
1336 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
1337 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
1338 paerlev=paerlev, naer_c=n_aerosolc, &
1339 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
1340 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
1341 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
1342 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
1343 ,doabsems=doabsems &
1344 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1345 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1346 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1347 )
1348 ELSE
1349 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
1350 ENDIF
1351 DO j=jts,jte
1352 DO k=kts,kte
1353 DO i=its,ite
1354 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
1355 ENDDO
1356 ENDDO
1357 ENDDO
1359 CASE (RRTMG_SWSCHEME)
1360 CALL wrf_debug(100, 'CALL rrtmg_sw')
1361 CALL RRTMG_SWRAD( &
1362 RTHRATENSW=RTHRATENSW, &
1363 SWUPT=SWUPT,SWUPTC=SWUPTC, &
1364 SWDNT=SWDNT,SWDNTC=SWDNTC, &
1365 SWUPB=SWUPB,SWUPBC=SWUPBC, &
1366 SWDNB=SWDNB,SWDNBC=SWDNBC, &
1367 SWCF=SWCF,GSW=GSW, &
1368 XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, &
1369 RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, &
1370 COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, &
1371 ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, &
1372 p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, &
1373 dz8w=dz8w,CLDFRA3D=CLDFRA,R=R_D,G=G, &
1374 ICLOUD=icloud,WARM_RAIN=warm_rain, &
1375 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
1376 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1377 QV3D=qv,QC3D=qc,QR3D=qr, &
1378 QI3D=qi,QS3D=qs,QG3D=qg, &
1379 ALSWVISDIR=alswvisdir ,ALSWVISDIF=alswvisdif, & !Zhenxin ssib alb comp (06/2010)
1380 ALSWNIRDIR=alswnirdir ,ALSWNIRDIF=alswnirdif, & !Zhenxin ssib alb comp (06/2010)
1381 SWVISDIR=swvisdir ,SWVISDIF=swvisdif, & !Zhenxin ssib swr comp (06/2010)
1382 SWNIRDIR=swnirdir ,SWNIRDIF=swnirdif, & !Zhenxin ssib swr comp (06/2010)
1383 SF_SURFACE_PHYSICS=sf_surface_physics, & !Zhenxin ssib sw_phy (06/2010)
1384 F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, &
1385 F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, &
1386 #ifdef WRF_CHEM
1387 TAUAER300=tauaer300,TAUAER400=tauaer400, & ! jcb
1388 TAUAER600=tauaer600,TAUAER999=tauaer999, & ! jcb
1389 GAER300=gaer300,GAER400=gaer400, & ! jcb
1390 GAER600=gaer600,GAER999=gaer999, & ! jcb
1391 WAER300=waer300,WAER400=waer400, & ! jcb
1392 WAER600=waer600,WAER999=waer999, & ! jcb
1393 aer_ra_feedback=aer_ra_feedback, &
1394 !jdfcz progn=progn,prescribe=prescribe, &
1395 progn=progn, &
1396 #endif
1397 QNDROP3D=qndrop,F_QNDROP=f_qndrop, &
1398 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
1399 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
1400 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
1401 SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, &
1402 SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC &
1403 )
1404 DO j=jts,jte
1405 DO k=kts,kte
1406 DO i=its,ite
1407 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
1408 ENDDO
1409 ENDDO
1410 ENDDO
1412 CASE (GFDLSWSCHEME)
1414 CALL wrf_debug (100, 'CALL gfdlsw')
1416 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
1417 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
1418 PRESENT(qv) .AND. PRESENT(qc) ) THEN
1419 IF ( F_QV .AND. F_QC .AND. F_QS ) THEN
1420 gfdl_sw = .true.
1421 CALL ETARA( &
1422 DT=dt,XLAND=xland &
1423 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
1424 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
1425 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
1426 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
1427 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
1428 ,HBOTR=hbotr, HTOPR=htopr &
1429 ,ALBEDO=albedo,CUPPT=cuppt &
1430 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
1431 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
1432 ,XTIME=xtime,JULIAN=julian &
1433 ,JULYR=julyr,JULDAY=julday &
1434 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
1435 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
1436 ,ACFRST=acfrst,NCFRST=ncfrst &
1437 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
1438 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
1439 ,THRATEN=rthraten,THRATENLW=rthratenlw &
1440 ,THRATENSW=rthratensw &
1441 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1442 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1443 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1444 )
1445 ELSE
1446 CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
1447 ENDIF
1448 ELSE
1449 CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
1450 ENDIF
1452 #ifdef HWRF
1453 CASE (HWRFSWSCHEME)
1455 CALL wrf_debug (100, 'CALL hwrfsw')
1457 gfdl_sw = .true.
1458 CALL HWRFRA(DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, &
1459 XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, &
1460 QV=qv,QW=qc_temp,QI=Qi, &
1461 TSK2D=tsk,GLW=GLW,GSW=GSW, &
1462 TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added
1463 GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt, &
1464 VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified
1465 NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified
1466 julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, &
1467 CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added
1468 ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added
1469 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
1470 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
1471 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
1472 #endif
1474 CASE (0)
1476 ! Here in case we don't want to call a sw radiation scheme
1477 ! For example, the Held-Suarez idealized test case
1478 IF (lw_physics /= HELDSUAREZ) THEN
1479 WRITE( wrf_err_message , * ) &
1480 'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')'
1481 CALL wrf_error_fatal ( wrf_err_message )
1482 END IF
1484 ! -- add by Jin Kong 10/2011
1485 !--- the following FLGSWSCHEME is for testing only
1486 CASE (FLGSWSCHEME)
1487 flg_sw = .true.
1488 !--- No need to do anything since the short and long wave is calculted in one program
1489 ! -- end
1491 CASE DEFAULT
1493 WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
1494 CALL wrf_error_fatal ( wrf_err_message )
1496 END SELECT swrad_select
1498 IF (sw_physics .eq. goddardswscheme) THEN
1499 IF ( PRESENT (tswdn) ) THEN
1500 DO j=jts,jte
1501 DO i=its,ite
1502 tswdn(i,j) = erbe_out(i,j,5) ! TOA SW downwelling flux [W/m2]
1503 tswup(i,j) = erbe_out(i,j,6) ! TOA SW upwelling flux [W/m2]
1504 sswdn(i,j) = erbe_out(i,j,7) ! surface SW downwelling flux [W/m2]
1505 sswup(i,j) = erbe_out(i,j,8) ! surface SW upwelling flux [W/m2]
1506 ENDDO
1507 ENDDO
1508 ENDIF
1509 ENDIF
1511 IF (sw_physics .gt. 0 .and. .not.gfdl_sw) THEN
1512 DO j=jts,jte
1513 DO k=kts,kte
1514 DO i=its,ite
1515 RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
1516 ENDDO
1517 ENDDO
1518 ENDDO
1520 DO j=jts,jte
1521 DO i=its,ite
1522 SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
1523 ENDDO
1524 ENDDO
1526 ENDIF
1528 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
1529 DO j=jts,jte
1530 DO k=kts,kte
1531 DO i=its,ite
1532 qc(i,k,j) = qc_save(i,k,j)
1533 ENDDO
1534 ENDDO
1535 ENDDO
1536 ENDIF
1537 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
1538 DO j=jts,jte
1539 DO k=kts,kte
1540 DO i=its,ite
1541 qi(i,k,j) = qi_save(i,k,j)
1542 ENDDO
1543 ENDDO
1544 ENDDO
1545 ENDIF
1547 ENDDO
1548 !$OMP END PARALLEL DO
1550 ENDIF Radiation_step
1552 accumulate_lw_select: SELECT CASE(lw_physics)
1554 CASE (CAMLWSCHEME,RRTMG_LWSCHEME)
1555 IF(PRESENT(LWUPTC))THEN
1556 !$OMP PARALLEL DO &
1557 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1559 DO ij = 1 , num_tiles
1560 its = i_start(ij)
1561 ite = i_end(ij)
1562 jts = j_start(ij)
1563 jte = j_end(ij)
1565 DO j=jts,jte
1566 DO i=its,ite
1567 ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DT
1568 ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DT
1569 ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DT
1570 ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DT
1571 ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DT
1572 ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DT
1573 ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DT
1574 ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DT
1575 ENDDO
1576 ENDDO
1577 ENDDO
1578 !$OMP END PARALLEL DO
1579 ENDIF
1580 CASE DEFAULT
1581 END SELECT accumulate_lw_select
1583 accumulate_sw_select: SELECT CASE(sw_physics)
1585 CASE (CAMSWSCHEME,RRTMG_SWSCHEME)
1586 IF(PRESENT(SWUPTC))THEN
1587 !$OMP PARALLEL DO &
1588 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1590 DO ij = 1 , num_tiles
1591 its = i_start(ij)
1592 ite = i_end(ij)
1593 jts = j_start(ij)
1594 jte = j_end(ij)
1596 DO j=jts,jte
1597 DO i=its,ite
1598 ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DT
1599 ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DT
1600 ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DT
1601 ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DT
1602 ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DT
1603 ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DT
1604 ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DT
1605 ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DT
1606 ENDDO
1607 ENDDO
1608 ENDDO
1609 !$OMP END PARALLEL DO
1610 ENDIF
1612 CASE DEFAULT
1613 END SELECT accumulate_sw_select
1615 END SUBROUTINE radiation_driver
1617 SUBROUTINE pre_radiation_driver ( grid, config_flags &
1618 ,itimestep, ra_call_offset &
1619 ,XLAT, XLONG, GMT, julian, xtime, RADT, STEPRA &
1620 ,ht,dx,dy,sina,cosa,shadowmask,slope_rad ,topo_shading &
1621 ,shadlen,ht_shad,ht_loc &
1622 ,ht_shad_bxs, ht_shad_bxe &
1623 ,ht_shad_bys, ht_shad_bye &
1624 ,nested, min_ptchsz &
1625 ,spec_bdy_width &
1626 ,ids, ide, jds, jde, kds, kde &
1627 ,ims, ime, jms, jme, kms, kme &
1628 ,ips, ipe, jps, jpe, kps, kpe &
1629 ,i_start, i_end &
1630 ,j_start, j_end &
1631 ,kts, kte &
1632 ,num_tiles )
1634 USE module_domain , ONLY : domain
1635 #ifdef DM_PARALLEL
1636 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks,wrf_dm_minval_integer
1637 # if (EM_CORE == 1)
1638 USE module_comm_dm , ONLY : halo_toposhad_sub
1639 # endif
1640 #endif
1641 USE module_bc
1642 USE module_model_constants
1644 IMPLICIT NONE
1646 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1647 ims,ime, jms,jme, kms,kme, &
1648 ips,ipe, jps,jpe, kps,kpe, &
1649 kts,kte, &
1650 num_tiles
1652 TYPE(domain) , INTENT(INOUT) :: grid
1653 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1655 INTEGER, INTENT(IN ) :: itimestep, ra_call_offset, stepra, &
1656 slope_rad, topo_shading, &
1657 spec_bdy_width
1659 INTEGER, INTENT(INOUT) :: min_ptchsz
1661 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
1662 i_start,i_end,j_start,j_end
1664 REAL, INTENT(IN ) :: GMT, radt, julian, xtime, dx, dy, shadlen
1666 REAL, DIMENSION( ims:ime, jms:jme ), &
1667 INTENT(IN ) :: XLAT, &
1668 XLONG, &
1669 HT, &
1670 SINA, &
1671 COSA
1673 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
1675 REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width ), &
1676 INTENT(IN ) :: ht_shad_bxs, ht_shad_bxe
1677 REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width ), &
1678 INTENT(IN ) :: ht_shad_bys, ht_shad_bye
1680 INTEGER, DIMENSION( ims:ime, jms:jme ), &
1681 INTENT(INOUT) :: shadowmask
1683 LOGICAL, INTENT(IN ) :: nested
1685 !Local
1686 ! For orographic shading
1687 INTEGER :: niter,ni,psx,psy,idum,jdum,i,j,ij
1688 REAL :: DECLIN,SOLCON
1690 ! Determine minimum patch size for slope-dependent radiation
1691 if (itimestep .eq. 1) then
1692 psx = ipe-ips+1
1693 psy = jpe-jps+1
1694 min_ptchsz = min(psx,psy)
1695 idum = 0
1696 jdum = 0
1697 endif
1699 # ifdef DM_PARALLEL
1700 if (itimestep .eq. 1) then
1701 call wrf_dm_minval_integer (psx,idum,jdum)
1702 call wrf_dm_minval_integer (psy,idum,jdum)
1703 min_ptchsz = min(psx,psy)
1704 endif
1705 # endif
1707 ! Topographic shading
1709 if ((topo_shading.eq.1).and.(itimestep .eq. 1 .or. &
1710 mod(itimestep,STEPRA) .eq. 1 + ra_call_offset)) then
1712 !---------------
1713 ! Calculate constants for short wave radiation
1715 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,DEGRAD,DPD)
1717 ! Make a local copy of terrain height field
1718 do j=jms,jme
1719 do i=ims,ime
1720 ht_loc(i,j) = ht(i,j)
1721 enddo
1722 enddo
1723 ! Determine if iterations are necessary for shadows to propagate from one patch to another
1724 if ((ids.eq.ips).and.(ide.eq.ipe).and.(jds.eq.jps).and.(jde.eq.jpe)) then
1725 niter = 1
1726 else
1727 niter = int(shadlen/(dx*min_ptchsz)+3)
1728 endif
1732 IF( nested ) THEN
1734 !$OMP PARALLEL DO &
1735 !$OMP PRIVATE ( ij )
1737 DO ij = 1 , num_tiles
1739 CALL spec_bdyfield(ht_shad, &
1740 ht_shad_bxs, ht_shad_bxe, &
1741 ht_shad_bys, ht_shad_bye, &
1742 'm', config_flags, spec_bdy_width, 2,&
1743 ids,ide, jds,jde, 1 ,1 , & ! domain dims
1744 ims,ime, jms,jme, 1 ,1 , & ! memory dims
1745 ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
1746 i_start(ij), i_end(ij), &
1747 j_start(ij), j_end(ij), &
1748 1 , 1 )
1749 ENDDO
1750 ENDIF
1752 do ni = 1, niter
1754 !$OMP PARALLEL DO &
1755 !$OMP PRIVATE ( ij,i,j )
1756 do ij = 1 , num_tiles
1758 call toposhad_init (ht_shad,ht_loc, &
1759 shadowmask,nested,ni, &
1760 ids,ide, jds,jde, kds,kde, &
1761 ims,ime, jms,jme, kms,kme, &
1762 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1763 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1764 min(j_end(ij), jde-1), kts, kte )
1766 enddo
1767 !$OMP END PARALLEL DO
1770 !$OMP PARALLEL DO &
1771 !$OMP PRIVATE ( ij,i,j )
1772 do ij = 1 , num_tiles
1774 call toposhad (xlat,xlong,sina,cosa,xtime,gmt,radt,declin, &
1775 dx,dy,ht_shad,ht_loc,ni, &
1776 shadowmask,shadlen, &
1777 ids,ide, jds,jde, kds,kde, &
1778 ims,ime, jms,jme, kms,kme, &
1779 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1780 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1781 min(j_end(ij), jde-1), kts, kte )
1783 enddo
1784 !$OMP END PARALLEL DO
1786 #if defined( DM_PARALLEL ) && (EM_CORE == 1)
1787 # include "HALO_TOPOSHAD.inc"
1788 #endif
1789 enddo
1790 endif
1792 END SUBROUTINE pre_radiation_driver
1794 !---------------------------------------------------------------------
1795 !BOP
1796 ! !IROUTINE: radconst - compute radiation terms
1797 ! !INTERFAC:
1798 SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
1799 DEGRAD,DPD )
1800 !---------------------------------------------------------------------
1801 USE module_wrf_error
1802 IMPLICIT NONE
1803 !---------------------------------------------------------------------
1805 ! !ARGUMENTS:
1806 REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
1807 REAL, INTENT(OUT ) :: DECLIN,SOLCON
1808 REAL :: OBECL,SINOB,SXLONG,ARG, &
1809 DECDEG,DJUL,RJUL,ECCFAC
1810 !
1811 ! !DESCRIPTION:
1812 ! Compute terms used in radiation physics
1813 !EOP
1815 ! for short wave radiation
1817 DECLIN=0.
1818 SOLCON=0.
1820 !-----OBECL : OBLIQUITY = 23.5 DEGREE.
1822 OBECL=23.5*DEGRAD
1823 SINOB=SIN(OBECL)
1825 !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
1827 IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
1828 IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
1829 SXLONG=SXLONG*DEGRAD
1830 ARG=SINOB*SIN(SXLONG)
1831 DECLIN=ASIN(ARG)
1832 DECDEG=DECLIN/DEGRAD
1833 !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
1834 DJUL=JULIAN*360./365.
1835 RJUL=DJUL*DEGRAD
1836 ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* &
1837 COS(2*RJUL)+0.000077*SIN(2*RJUL)
1838 SOLCON=1370.*ECCFAC
1840 END SUBROUTINE radconst
1842 !---------------------------------------------------------------------
1843 !BOP
1844 ! !IROUTINE: cal_cldfra - Compute cloud fraction
1845 ! !INTERFACE:
1846 SUBROUTINE cal_cldfra(CLDFRA,QC,QI,F_QC,F_QI, &
1847 ids,ide, jds,jde, kds,kde, &
1848 ims,ime, jms,jme, kms,kme, &
1849 its,ite, jts,jte, kts,kte )
1850 !---------------------------------------------------------------------
1851 IMPLICIT NONE
1852 !---------------------------------------------------------------------
1853 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1854 ims,ime, jms,jme, kms,kme, &
1855 its,ite, jts,jte, kts,kte
1857 !
1858 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1859 CLDFRA
1861 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1862 QI, &
1863 QC
1865 LOGICAL,INTENT(IN) :: F_QC,F_QI
1867 REAL thresh
1868 INTEGER:: i,j,k
1869 ! !DESCRIPTION:
1870 ! Compute cloud fraction from input ice and cloud water fields
1871 ! if provided.
1872 !
1873 ! Whether QI or QC is active or not is determined from the indices of
1874 ! the fields into the 4D scalar arrays in WRF. These indices are
1875 ! P_QI and P_QC, respectively, and they are passed in to the routine
1876 ! to enable testing to see if QI and QC represent active fields in
1877 ! the moisture 4D scalar array carried by WRF.
1878 !
1879 ! If a field is active its index will have a value greater than or
1880 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1881 ! this routine.
1882 !EOP
1883 !---------------------------------------------------------------------
1884 thresh=1.0e-6
1886 IF ( f_qi .AND. f_qc ) THEN
1887 DO j = jts,jte
1888 DO k = kts,kte
1889 DO i = its,ite
1890 IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
1891 CLDFRA(i,k,j)=1.
1892 ELSE
1893 CLDFRA(i,k,j)=0.
1894 ENDIF
1895 ENDDO
1896 ENDDO
1897 ENDDO
1898 ELSE IF ( f_qc ) THEN
1899 DO j = jts,jte
1900 DO k = kts,kte
1901 DO i = its,ite
1902 IF ( QC(i,k,j) .gt. thresh) THEN
1903 CLDFRA(i,k,j)=1.
1904 ELSE
1905 CLDFRA(i,k,j)=0.
1906 ENDIF
1907 ENDDO
1908 ENDDO
1909 ENDDO
1910 ELSE
1911 DO j = jts,jte
1912 DO k = kts,kte
1913 DO i = its,ite
1914 CLDFRA(i,k,j)=0.
1915 ENDDO
1916 ENDDO
1917 ENDDO
1918 ENDIF
1920 END SUBROUTINE cal_cldfra
1922 !BOP
1923 ! !IROUTINE: cal_cldfra2 - Compute cloud fraction
1924 ! !INTERFACE:
1925 ! cal_cldfra_xr - Compute cloud fraction.
1926 ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
1927 !!
1928 !!--- Cloud fraction parameterization follows Randall, 1994
1929 !! (see Hong et al., 1998)
1930 !! (modified by Ferrier, Feb '02)
1931 !
1932 SUBROUTINE cal_cldfra2(CLDFRA, QV, QC, QI, QS, &
1933 F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
1934 F_ICE_PHY,F_RAIN_PHY, &
1935 ids,ide, jds,jde, kds,kde, &
1936 ims,ime, jms,jme, kms,kme, &
1937 its,ite, jts,jte, kts,kte )
1938 !---------------------------------------------------------------------
1939 IMPLICIT NONE
1940 !---------------------------------------------------------------------
1941 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1942 ims,ime, jms,jme, kms,kme, &
1943 its,ite, jts,jte, kts,kte
1945 !
1946 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1947 CLDFRA
1949 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1950 QV, &
1951 QI, &
1952 QC, &
1953 QS, &
1954 t_phy, &
1955 p_phy
1956 ! p_phy, &
1957 ! F_ICE_PHY, &
1958 ! F_RAIN_PHY
1960 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1961 OPTIONAL, &
1962 INTENT(IN ) :: &
1963 F_ICE_PHY, &
1964 F_RAIN_PHY
1966 LOGICAL,OPTIONAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
1968 ! REAL thresh
1969 INTEGER:: i,j,k
1970 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
1972 REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
1973 PEXP=0.25, RHGRID=1.0
1974 REAL , PARAMETER :: SVP1=0.61078
1975 REAL , PARAMETER :: SVP2=17.2693882
1976 REAL , PARAMETER :: SVPI2=21.8745584
1977 REAL , PARAMETER :: SVP3=35.86
1978 REAL , PARAMETER :: SVPI3=7.66
1979 REAL , PARAMETER :: SVPT0=273.15
1980 REAL , PARAMETER :: r_d = 287.
1981 REAL , PARAMETER :: r_v = 461.6
1982 REAL , PARAMETER :: ep_2=r_d/r_v
1983 ! !DESCRIPTION:
1984 ! Compute cloud fraction from input ice and cloud water fields
1985 ! if provided.
1986 !
1987 ! Whether QI or QC is active or not is determined from the indices of
1988 ! the fields into the 4D scalar arrays in WRF. These indices are
1989 ! P_QI and P_QC, respectively, and they are passed in to the routine
1990 ! to enable testing to see if QI and QC represent active fields in
1991 ! the moisture 4D scalar array carried by WRF.
1992 !
1993 ! If a field is active its index will have a value greater than or
1994 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1995 ! this routine.
1996 !EOP
1999 !-----------------------------------------------------------------------
2000 !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
2001 ! (modified by Ferrier, Feb '02)
2002 !
2003 !--- Cloud fraction parameterization follows Randall, 1994
2004 ! (see Hong et al., 1998)
2005 !-----------------------------------------------------------------------
2006 ! Note: ep_2=287./461.6 Rd/Rv
2007 ! Note: R_D=287.
2009 ! Alternative calculation for critical RH for grid saturation
2010 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
2012 ! Calculate saturation mixing ratio weighted according to the fractions of
2013 ! water and ice.
2014 ! Following:
2015 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
2016 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
2017 !
2018 ! over ice over water
2019 ! a = 21.8745584 17.2693882
2020 ! b = 7.66 35.86
2022 !---------------------------------------------------------------------
2024 DO j = jts,jte
2025 DO k = kts,kte
2026 DO i = its,ite
2027 tc = t_phy(i,k,j) - SVPT0
2028 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
2029 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
2030 QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
2031 QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
2033 IF ( PRESENT(F_QI) .and. PRESENT(F_QC) .and. PRESENT(F_QS) ) THEN
2035 ! mji - For MP options 2, 4, 6, 7, 8, etc. (qc = liquid, qi = ice, qs = snow)
2036 IF ( F_QI .and. F_QC .and. F_QS) THEN
2037 QCLD = QI(i,k,j)+QC(i,k,j)+QS(i,k,j)
2038 IF (QCLD .LT. QCLDMIN) THEN
2039 weight = 0.
2040 ELSE
2041 weight = (QI(i,k,j)+QS(i,k,j)) / QCLD
2042 ENDIF
2043 ENDIF
2045 ! mji - For MP options 1 and 3, (qc only)
2046 ! For MP=1, qc = liquid, for MP=3, qc = liquid or ice depending on temperature
2047 IF ( F_QC .and. .not. F_QI .and. .not. F_QS ) THEN
2048 QCLD = QC(i,k,j)
2049 IF (QCLD .LT. QCLDMIN) THEN
2050 weight = 0.
2051 ELSE
2052 if (t_phy(i,k,j) .gt. 273.15) weight = 0.
2053 if (t_phy(i,k,j) .le. 273.15) weight = 1.
2054 ENDIF
2055 ENDIF
2057 ! mji - For MP option 5; (qc = liquid, qs = ice)
2058 IF ( F_QC .and. .not. F_QI .and. F_QS .and. PRESENT(F_ICE_PHY) ) THEN
2060 ! Mixing ratios of cloud water & total ice (cloud ice + snow).
2061 ! Mixing ratios of rain are not considered in this scheme.
2062 ! F_ICE is fraction of ice
2063 ! F_RAIN is fraction of rain
2065 QIMID = QS(i,k,j)
2066 QWMID = QC(i,k,j)
2067 ! old method
2068 ! QIMID = QC(i,k,j)*F_ICE_PHY(i,k,j)
2069 ! QWMID = (QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
2070 !
2071 !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
2072 ! only cloud water + cloud ice + snow
2073 !
2074 QCLD=QWMID+QIMID
2075 IF (QCLD .LT. QCLDMIN) THEN
2076 weight = 0.
2077 ELSE
2078 weight = F_ICE_PHY(i,k,j)
2079 ENDIF
2080 ENDIF
2082 ELSE
2083 CLDFRA(i,k,j)=0.
2085 ENDIF ! IF ( F_QI .and. F_QC .and. F_QS)
2088 QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI
2089 RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
2090 !
2091 !--- Determine cloud fraction (modified from original algorithm)
2092 !
2093 IF (QCLD .LT. QCLDMIN) THEN
2094 !
2095 !--- Assume zero cloud fraction if there is no cloud mixing ratio
2096 !
2097 CLDFRA(i,k,j)=0.
2098 ELSEIF(RHUM.GE.RHGRID)THEN
2099 !
2100 !--- Assume cloud fraction of unity if near saturation and the cloud
2101 ! mixing ratio is at or above the minimum threshold
2102 !
2103 CLDFRA(i,k,j)=1.
2104 ELSE
2105 !
2106 !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
2107 ! modified based on assumed grid-scale saturation at RH=RHgrid.
2108 !
2109 SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
2110 DENOM=(SUBSAT)**GAMMA
2111 ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
2112 ! prevent negative values (new)
2113 RHUM=MAX(1.E-10, RHUM)
2114 CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG))
2115 !! ARG=-1000*QCLD/(RHUM-RHGRID)
2116 !! ARG=MAX(ARG, ARGMIN)
2117 !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
2118 IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
2119 ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
2120 ENDDO !--- End DO i
2121 ENDDO !--- End DO k
2122 ENDDO !--- End DO j
2124 END SUBROUTINE cal_cldfra2
2127 SUBROUTINE toposhad_init(ht_shad,ht_loc,shadowmask,nested,iter, &
2128 ids,ide, jds,jde, kds,kde, &
2129 ims,ime, jms,jme, kms,kme, &
2130 ips,ipe, jps,jpe, kps,kpe, &
2131 its,ite, jts,jte, kts,kte )
2133 USE module_model_constants
2135 implicit none
2137 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
2138 ims,ime, jms,jme, kms,kme, &
2139 ips,ipe, jps,jpe, kps,kpe, &
2140 its,ite, jts,jte, kts,kte
2142 LOGICAL, INTENT(IN) :: nested
2144 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad, ht_loc
2146 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
2147 INTEGER, INTENT(IN) :: iter
2149 ! Local variables
2151 INTEGER :: i, j
2153 if (iter.eq.1) then
2155 ! Initialize shadow mask
2156 do j=jts,jte
2157 do i=its,ite
2158 shadowmask(i,j) = 0
2159 ENDDO
2160 ENDDO
2162 ! Initialize shading height
2164 IF ( nested ) THEN ! Do not overwrite input from parent domain
2165 do j=max(jts,jds+2),min(jte,jde-3)
2166 do i=max(its,ids+2),min(ite,ide-3)
2167 ht_shad(i,j) = ht_loc(i,j)-0.001
2168 ENDDO
2169 ENDDO
2170 ELSE
2171 do j=jts,jte
2172 do i=its,ite
2173 ht_shad(i,j) = ht_loc(i,j)-0.001
2174 ENDDO
2175 ENDDO
2176 ENDIF
2178 IF ( nested ) THEN ! Check if a shadow exceeding the topography height is available at the lateral domain edge from nesting
2179 if (its.eq.ids) then
2180 do j=jts,jte
2181 if (ht_shad(its,j) .gt. ht_loc(its,j)) then
2182 shadowmask(its,j) = 1
2183 ht_loc(its,j) = ht_shad(its,j)
2184 endif
2185 if (ht_shad(its+1,j) .gt. ht_loc(its+1,j)) then
2186 shadowmask(its+1,j) = 1
2187 ht_loc(its+1,j) = ht_shad(its+1,j)
2188 endif
2189 enddo
2190 endif
2191 if (ite.eq.ide-1) then
2192 do j=jts,jte
2193 if (ht_shad(ite,j) .gt. ht_loc(ite,j)) then
2194 shadowmask(ite,j) = 1
2195 ht_loc(ite,j) = ht_shad(ite,j)
2196 endif
2197 if (ht_shad(ite-1,j) .gt. ht_loc(ite-1,j)) then
2198 shadowmask(ite-1,j) = 1
2199 ht_loc(ite-1,j) = ht_shad(ite-1,j)
2200 endif
2201 enddo
2202 endif
2203 if (jts.eq.jds) then
2204 do i=its,ite
2205 if (ht_shad(i,jts) .gt. ht_loc(i,jts)) then
2206 shadowmask(i,jts) = 1
2207 ht_loc(i,jts) = ht_shad(i,jts)
2208 endif
2209 if (ht_shad(i,jts+1) .gt. ht_loc(i,jts+1)) then
2210 shadowmask(i,jts+1) = 1
2211 ht_loc(i,jts+1) = ht_shad(i,jts+1)
2212 endif
2213 enddo
2214 endif
2215 if (jte.eq.jde-1) then
2216 do i=its,ite
2217 if (ht_shad(i,jte) .gt. ht_loc(i,jte)) then
2218 shadowmask(i,jte) = 1
2219 ht_loc(i,jte) = ht_shad(i,jte)
2220 endif
2221 if (ht_shad(i,jte-1) .gt. ht_loc(i,jte-1)) then
2222 shadowmask(i,jte-1) = 1
2223 ht_loc(i,jte-1) = ht_shad(i,jte-1)
2224 endif
2225 enddo
2226 endif
2227 ENDIF
2229 else
2231 ! Fill the local topography field at the points next to internal tile boundaries with ht_shad values
2232 ! A 2-pt halo has been applied to the ht_shad before the repeated call of this subroutine
2234 if ((its.ne.ids).and.(its.eq.ips)) then
2235 do j=jts-2,jte+2
2236 ht_loc(its-1,j) = max(ht_loc(its-1,j),ht_shad(its-1,j))
2237 ht_loc(its-2,j) = max(ht_loc(its-2,j),ht_shad(its-2,j))
2238 enddo
2239 endif
2240 if ((ite.ne.ide-1).and.(ite.eq.ipe)) then
2241 do j=jts-2,jte+2
2242 ht_loc(ite+1,j) = max(ht_loc(ite+1,j),ht_shad(ite+1,j))
2243 ht_loc(ite+2,j) = max(ht_loc(ite+2,j),ht_shad(ite+2,j))
2244 enddo
2245 endif
2246 if ((jts.ne.jds).and.(jts.eq.jps)) then
2247 do i=its-2,ite+2
2248 ht_loc(i,jts-1) = max(ht_loc(i,jts-1),ht_shad(i,jts-1))
2249 ht_loc(i,jts-2) = max(ht_loc(i,jts-2),ht_shad(i,jts-2))
2250 enddo
2251 endif
2252 if ((jte.ne.jde-1).and.(jte.eq.jpe)) then
2253 do i=its-2,ite+2
2254 ht_loc(i,jte+1) = max(ht_loc(i,jte+1),ht_shad(i,jte+1))
2255 ht_loc(i,jte+2) = max(ht_loc(i,jte+2),ht_shad(i,jte+2))
2256 enddo
2257 endif
2259 endif
2261 END SUBROUTINE toposhad_init
2266 SUBROUTINE toposhad(xlat,xlong,sina,cosa,xtime,gmt,radfrq,declin, &
2267 dx,dy,ht_shad,ht_loc,iter, &
2268 shadowmask,shadlen, &
2269 ids,ide, jds,jde, kds,kde, &
2270 ims,ime, jms,jme, kms,kme, &
2271 ips,ipe, jps,jpe, kps,kpe, &
2272 its,ite, jts,jte, kts,kte )
2275 USE module_model_constants
2277 implicit none
2279 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
2280 ims,ime, jms,jme, kms,kme, &
2281 ips,ipe, jps,jpe, kps,kpe, &
2282 its,ite, jts,jte, kts,kte
2284 INTEGER, INTENT(IN) :: iter
2286 REAL, INTENT(IN) :: RADFRQ,XTIME,DECLIN,dx,dy,gmt,shadlen
2288 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT, XLONG, sina, cosa
2290 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
2292 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
2294 ! Local variables
2296 REAL :: pi, xt24, wgt, ri, rj, argu, sol_azi, topoelev, dxabs, tloctm, hrang, xxlat, csza
2297 INTEGER :: gpshad, ii, jj, i1, i2, j1, j2, i, j
2301 XT24=MOD(XTIME+RADFRQ*0.5,1440.)
2302 pi = 4.*atan(1.)
2303 gpshad = int(shadlen/dx+1.)
2305 if (iter.eq.1) then
2308 j_loop1: DO J=jts,jte
2309 i_loop1: DO I=its,ite
2311 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
2312 HRANG=15.*(TLOCTM-12.)*DEGRAD
2313 XXLAT=XLAT(i,j)*DEGRAD
2314 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
2316 if (csza.lt.1.e-2) then ! shadow mask does not need to be computed
2317 shadowmask(i,j) = 0
2318 ht_shad(i,j) = ht_loc(i,j)-0.001
2319 goto 120
2320 endif
2322 ! Solar azimuth angle
2324 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
2325 if (argu.gt.1) argu = 1
2326 if (argu.lt.-1) argu = -1
2327 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
2328 if (cosa(i,j).ge.0) then
2329 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
2330 else
2331 sol_azi = sol_azi + pi - asin(sina(i,j))
2332 endif
2334 ! Scan for higher surrounding topography
2336 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
2338 do jj = j+1,j+gpshad
2339 ri = i + (jj-j)*tan(sol_azi)
2340 i1 = int(ri)
2341 i2 = i1+1
2342 wgt = ri-i1
2343 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2344 if ((jj.ge.jpe+1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
2345 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2346 goto 120
2347 endif
2348 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2349 if (sin(topoelev).ge.csza) then
2350 shadowmask(i,j) = 1
2351 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2352 endif
2353 enddo
2355 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
2356 do ii = i+1,i+gpshad
2357 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2358 j1 = int(rj)
2359 j2 = j1+1
2360 wgt = rj-j1
2361 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2362 if ((ii.ge.ipe+1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
2363 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2364 goto 120
2365 endif
2366 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2367 if (sin(topoelev).ge.csza) then
2368 shadowmask(i,j) = 1
2369 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2370 endif
2371 enddo
2373 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
2374 do jj = j-1,j-gpshad,-1
2375 ri = i + (jj-j)*tan(sol_azi)
2376 i1 = int(ri)
2377 i2 = i1+1
2378 wgt = ri-i1
2379 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2380 if ((jj.le.jps-1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
2381 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2382 goto 120
2383 endif
2384 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2385 if (sin(topoelev).ge.csza) then
2386 shadowmask(i,j) = 1
2387 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2388 endif
2389 enddo
2391 else ! sun is in the western quarter
2392 do ii = i-1,i-gpshad,-1
2393 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2394 j1 = int(rj)
2395 j2 = j1+1
2396 wgt = rj-j1
2397 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2398 if ((ii.le.ips-1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
2399 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2400 goto 120
2401 endif
2402 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2403 if (sin(topoelev).ge.csza) then
2404 shadowmask(i,j) = 1
2405 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2406 endif
2407 enddo
2408 endif
2410 120 continue
2412 ENDDO i_loop1
2413 ENDDO j_loop1
2415 else ! iteration > 1
2418 j_loop2: DO J=jts,jte
2419 i_loop2: DO I=its,ite
2421 ! if (shadowmask(i,j).eq.-1) then ! this indicates that the search ended at a lateral boundary during iteration 1
2423 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
2424 HRANG=15.*(TLOCTM-12.)*DEGRAD
2425 XXLAT=XLAT(i,j)*DEGRAD
2426 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
2428 ! Solar azimuth angle
2430 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
2431 if (argu.gt.1) argu = 1
2432 if (argu.lt.-1) argu = -1
2433 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
2434 if (cosa(i,j).ge.0) then
2435 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
2436 else
2437 sol_azi = sol_azi + pi - asin(sina(i,j))
2438 endif
2440 ! Scan for higher surrounding topography
2442 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
2444 do jj = j+1,j+gpshad
2445 ri = i + (jj-j)*tan(sol_azi)
2446 i1 = int(ri)
2447 i2 = i1+1
2448 wgt = ri-i1
2449 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2450 if ((jj.ge.min(jde,jpe+3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220
2451 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2452 if (sin(topoelev).ge.csza) then
2453 shadowmask(i,j) = 1
2454 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2455 endif
2456 enddo
2458 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
2459 do ii = i+1,i+gpshad
2460 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2461 j1 = int(rj)
2462 j2 = j1+1
2463 wgt = rj-j1
2464 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2465 if ((ii.ge.min(ide,ipe+3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220
2466 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2467 if (sin(topoelev).ge.csza) then
2468 shadowmask(i,j) = 1
2469 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2470 endif
2471 enddo
2473 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
2474 do jj = j-1,j-gpshad,-1
2475 ri = i + (jj-j)*tan(sol_azi)
2476 i1 = int(ri)
2477 i2 = i1+1
2478 wgt = ri-i1
2479 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2480 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
2481 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2482 if (sin(topoelev).ge.csza) then
2483 shadowmask(i,j) = 1
2484 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2485 endif
2486 enddo
2488 else ! sun is in the western quarter
2489 do ii = i-1,i-gpshad,-1
2490 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2491 j1 = int(rj)
2492 j2 = j1+1
2493 wgt = rj-j1
2494 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2495 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
2496 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2497 if (sin(topoelev).ge.csza) then
2498 shadowmask(i,j) = 1
2499 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2500 endif
2501 enddo
2502 endif
2504 220 continue
2505 ! endif
2507 ENDDO i_loop2
2508 ENDDO j_loop2
2510 endif ! iteration
2512 END SUBROUTINE toposhad
2514 END MODULE module_radiation_driver