| blob | 743272b9a0b07a11ca6d1e101dfb9cc6536503fa |
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 itimestep,dt ,lw_physics,sw_physics ,NPHS &
11 ,RTHRATENLW ,RTHRATENSW ,RTHRATEN &
12 ,ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC & ! Optional
13 ,ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC & ! Optional
14 ,ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC & ! Optional
15 ,ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC & ! Optional
16 , SWUPT, SWUPTC, SWDNT, SWDNTC & ! Optional
17 , SWUPB, SWUPBC, SWDNB, SWDNBC & ! Optional
18 , LWUPT, LWUPTC, LWDNT, LWDNTC & ! Optional
19 , LWUPB, LWUPBC, LWDNB, LWDNBC & ! Optional
20 ,LWCF,SWCF,OLR & ! Optional
21 ,SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC & ! Optional
22 ,LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC & ! Optional
23 ,GLW, GSW, SWDOWN, XLAT, XLONG, ALBEDO &
24 ,EMISS, rho, p8w, p , pi , dz8w ,t, t8w, GMT &
25 ,XLAND, XICE, TSK, HTOP,HBOT,HTOPR,HBOTR, CUPPT, VEGFRA, SNOW &
26 ,julyr, JULDAY, julian, YR, xtime, RADT, STEPRA, ICLOUD, warm_rain &
27 ,declinx,solconx,COSZEN, HRANG & !Optional
28 ,cen_lat &
29 ,ra_call_offset,RSWTOA,RLWTOA, CZMEAN &
30 ,CFRACL, CFRACM, CFRACH &
31 ,ACFRST,NCFRST,ACFRCV,NCFRCV,SWDOWNC &
32 ,z &
33 ,levsiz, n_ozmixm, n_aerosolc, paerlev &
34 ,cam_abs_dim1, cam_abs_dim2, cam_abs_freq_s &
35 ,ozmixm,pin & ! Optional
36 ,m_ps_1,m_ps_2,aerosolc_1,aerosolc_2,m_hybi0 & ! Optional
37 ,abstot, absnxt, emstot & ! Optional
38 ,taucldi, taucldc & ! Optional
39 ,ids, ide, jds, jde, kds, kde &
40 ,ims, ime, jms, jme, kms, kme &
41 ,i_start, i_end &
42 ,j_start, j_end &
43 ,kts, kte &
44 ,num_tiles, CURR_SECS, adapt_step_flag &
45 ,qv,qc,qr,qi,qs,qg,qndrop &
46 ,f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop &
47 ,CLDFRA ,Pb &
48 ,f_ice_phy,f_rain_phy &
49 ,pm2_5_dry, pm2_5_water, pm2_5_dry_ec &
50 ,tauaer300, tauaer400, tauaer600, tauaer999 & ! jcb
51 ,gaer300, gaer400, gaer600, gaer999 & ! jcb
52 ,waer300, waer400, waer600, waer999 & ! jcb
53 ,qc_adjust ,qi_adjust & ! jm
54 ,cu_rad_feedback, aer_ra_feedback & ! jm
55 ,ht,dx,dy,shadowmask,slope_rad ,topo_shading & ! slope-dependent radiation
56 !JJS 20101021 vvvvv
57 ,tlwdn, tlwup & ! optoinal for Goddard schemes
58 ,slwdn, slwup & ! optoinal for Goddard schemes
59 ,tswdn, tswup & ! optoinal for Goddard schemes
60 ,sswdn, sswup & ! optoinal for Goddard schemes
61 !JJS 20101021 ^^^^^
62 )
65 !-------------------------------------------------------------------------
67 ! !USES:
68 USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME &
69 ,RRTMG_LWSCHEME, RRTMG_SWSCHEME &
70 ,SWRADSCHEME, GSFCSWSCHEME &
71 ,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME &
72 ,HELDSUAREZ &
73 #ifdef HWRF
74 ,HWRFSWSCHEME, HWRFLWSCHEME &
75 #endif
76 ,goddardlwscheme &
77 ,goddardswscheme
80 USE module_model_constants
81 USE module_wrf_error , ONLY : wrf_err_message
83 ! *** add new modules of schemes here
85 USE module_ra_sw , ONLY : swrad
86 USE module_ra_gsfcsw , ONLY : gsfcswrad
87 USE module_ra_rrtm , ONLY : rrtmlwrad
88 USE module_ra_rrtmg_lw , ONLY : rrtmg_lwrad
89 USE module_ra_rrtmg_sw , ONLY : rrtmg_swrad
90 USE module_ra_cam , ONLY : camrad
91 USE module_ra_gfdleta , ONLY : etara
92 #ifdef HWRF
93 USE module_ra_hwrf
94 #endif
95 USE module_ra_hs , ONLY : hsrad
97 USE module_ra_goddard , ONLY : goddardrad !Goddard Radiation Scheme by Toshihisa Matsui and Jainn Shi
101 ! This driver calls subroutines for the radiation parameterizations.
102 !
103 ! short wave radiation choices:
104 ! 1. swrad (19??)
105 ! 4. rrtmg_sw - Added November 2008, MJIacono, AER, Inc.
106 !
107 ! long wave radiation choices:
108 ! 1. rrtmlwrad
109 ! 4. rrtmg_lw - Added November 2008, MJIacono, AER, Inc.
110 !
111 !----------------------------------------------------------------------
112 IMPLICIT NONE
113 !<DESCRIPTION>
114 !
115 ! Radiation_driver is the WRF mediation layer routine that provides the interface to
116 ! to radiation physics packages in the WRF model layer. The radiation
117 ! physics packages to call are chosen by setting the namelist variable
118 ! (Rconfig entry in Registry) to the integer value assigned to the
119 ! particular package (package entry in Registry). For example, if the
120 ! namelist variable ra_lw_physics is set to 1, this corresponds to the
121 ! Registry Package entry for swradscheme. Note that the Package
122 ! names in the Registry are defined constants (frame/module_state_description.F)
123 ! in the CASE statements in this routine.
124 !
125 ! Among the arguments is moist, a four-dimensional scalar array storing
126 ! a variable number of moisture tracers, depending on the physics
127 ! configuration for the WRF run, as determined in the namelist. The
128 ! highest numbered index of active moisture tracers the integer argument
129 ! n_moist (note: the number of tracers at run time is the quantity
130 ! <tt>n_moist - PARAM_FIRST_SCALAR + 1</tt> , not n_moist. Individual tracers
131 ! may be indexed from moist by the Registry name of the tracer prepended
132 ! with P_; for example P_QC is the index of cloud water. An index
133 ! represents a valid, active field only if the index is greater than
134 ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual
135 ! indices for each tracer is defined in module_state_description and
136 ! set in <a href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a> defined in frame/module_configure.F.
137 !
138 ! Physics drivers in WRF 2.0 and higher, originally model-layer
139 ! routines, have been promoted to mediation layer routines and they
140 ! contain OpenMP threaded loops over tiles. Thus, physics drivers
141 ! are called from single-threaded regions in the solver. The physics
142 ! routines that are called from the physics drivers are model-layer
143 ! routines and fully tile-callable and thread-safe.
144 !</DESCRIPTION>
145 !
146 !======================================================================
147 ! Grid structure in physics part of WRF
148 !----------------------------------------------------------------------
149 ! The horizontal velocities used in the physics are unstaggered
150 ! relative to temperature/moisture variables. All predicted
151 ! variables are carried at half levels except w, which is at full
152 ! levels. Some arrays with names (*8w) are at w (full) levels.
153 !
154 !----------------------------------------------------------------------
155 ! In WRF, kms (smallest number) is the bottom level and kme (largest
156 ! number) is the top level. In your scheme, if 1 is at the top level,
157 ! then you have to reverse the order in the k direction.
158 !
159 ! kme - half level (no data at this level)
160 ! kme ----- full level
161 ! kme-1 - half level
162 ! kme-1 ----- full level
163 ! .
164 ! .
165 ! .
166 ! kms+2 - half level
167 ! kms+2 ----- full level
168 ! kms+1 - half level
169 ! kms+1 ----- full level
170 ! kms - half level
171 ! kms ----- full level
172 !
173 !======================================================================
174 ! Grid structure in physics part of WRF
175 !
176 !-------------------------------------
177 ! The horizontal velocities used in the physics are unstaggered
178 ! relative to temperature/moisture variables. All predicted
179 ! variables are carried at half levels except w, which is at full
180 ! levels. Some arrays with names (*8w) are at w (full) levels.
181 !
182 !==================================================================
183 ! Definitions
184 !-----------
185 ! Theta potential temperature (K)
186 ! Qv water vapor mixing ratio (kg/kg)
187 ! Qc cloud water mixing ratio (kg/kg)
188 ! Qr rain water mixing ratio (kg/kg)
189 ! Qi cloud ice mixing ratio (kg/kg)
190 ! Qs snow mixing ratio (kg/kg)
191 !-----------------------------------------------------------------
192 !-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3)
193 !-- PM2_5_WATER PM2.5 water mass (ug m^-3)
194 !-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3)
195 !-- RTHRATEN Theta tendency
196 ! due to radiation (K/s)
197 !-- RTHRATENLW Theta tendency
198 ! due to long wave radiation (K/s)
199 !-- RTHRATENSW Theta temperature tendency
200 ! due to short wave radiation (K/s)
201 !-- dt time step (s)
202 !-- itimestep number of time steps
203 !-- GLW downward long wave flux at ground surface (W/m^2)
204 !-- GSW net short wave flux at ground surface (W/m^2)
205 !-- SWDOWN downward short wave flux at ground surface (W/m^2)
206 !-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ)
207 !-- RLWTOA upward long wave at top of atmosphere (w/m2)
208 !-- RSWTOA upward short wave at top of atmosphere (w/m2)
209 !-- XLAT latitude, south is negative (degree)
210 !-- XLONG longitude, west is negative (degree)
211 !-- ALBEDO albedo (between 0 and 1)
212 !-- CLDFRA cloud fraction (between 0 and 1)
213 !-- EMISS surface emissivity (between 0 and 1)
214 !-- rho_phy density (kg/m^3)
215 !-- rr dry air density (kg/m^3)
216 !-- moist moisture array (4D - last index is species) (kg/kg)
217 !-- n_moist number of moisture species
218 !-- qndrop Cloud droplet number (#/kg)
219 !-- p8w pressure at full levels (Pa)
220 !-- p_phy pressure (Pa)
221 !-- Pb base-state pressure (Pa)
222 !-- pi_phy exner function (dimensionless)
223 !-- dz8w dz between full levels (m)
224 !-- t_phy temperature (K)
225 !-- t8w temperature at full levels (K)
226 !-- GMT Greenwich Mean Time Hour of model start (hour)
227 !-- JULDAY the initial day (Julian day)
228 !-- RADT time for calling radiation (min)
229 !-- ra_call_offset -1 (old) means usually just before output, 0 after
230 !-- DEGRAD conversion factor for
231 ! degrees to radians (pi/180.) (rad/deg)
232 !-- DPD degrees per day for earth's
233 ! orbital position (deg/day)
234 !-- R_d gas constant for dry air (J/kg/K)
235 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
236 !-- G acceleration due to gravity (m/s^2)
237 !-- rvovrd R_v divided by R_d (dimensionless)
238 !-- XTIME time since simulation start (min)
239 !-- DECLIN solar declination angle (rad)
240 !-- SOLCON solar constant (W/m^2)
241 !-- ids start index for i in domain
242 !-- ide end index for i in domain
243 !-- jds start index for j in domain
244 !-- jde end index for j in domain
245 !-- kds start index for k in domain
246 !-- kde end index for k in domain
247 !-- ims start index for i in memory
248 !-- ime end index for i in memory
249 !-- jms start index for j in memory
250 !-- jme end index for j in memory
251 !-- kms start index for k in memory
252 !-- kme end index for k in memory
253 !-- i_start start indices for i in tile
254 !-- i_end end indices for i in tile
255 !-- j_start start indices for j in tile
256 !-- j_end end indices for j in tile
257 !-- kts start index for k in tile
258 !-- kte end index for k in tile
259 !-- num_tiles number of tiles
260 !
261 !==================================================================
262 !
263 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
264 ims,ime, jms,jme, kms,kme, &
265 kts,kte, &
266 num_tiles
268 INTEGER, INTENT(IN) :: lw_physics, sw_physics
270 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
271 i_start,i_end,j_start,j_end
273 INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset
274 INTEGER, INTENT(IN ) :: levsiz, n_ozmixm
275 INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2
276 REAL, INTENT(IN ) :: cam_abs_freq_s
278 LOGICAL, INTENT(IN ) :: warm_rain
280 REAL, INTENT(IN ) :: RADT
282 REAL, DIMENSION( ims:ime, jms:jme ), &
283 INTENT(IN ) :: XLAND, &
284 XICE, &
285 TSK, &
286 VEGFRA, &
287 SNOW
288 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, &
289 INTENT(IN ) :: OZMIXM
291 REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN
293 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2
294 REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, &
295 INTENT(IN ) :: aerosolc_1, aerosolc_2
296 REAL, DIMENSION(paerlev), OPTIONAL, &
297 INTENT(IN ) :: m_hybi0
299 REAL, DIMENSION( ims:ime, jms:jme ), &
300 INTENT(INOUT) :: HTOP, &
301 HBOT, &
302 HTOPR, &
303 HBOTR, &
304 CUPPT
306 INTEGER, INTENT(IN ) :: julyr
307 !
308 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
309 INTENT(IN ) :: dz8w, &
310 z, &
311 p8w, &
312 p, &
313 pi, &
314 t, &
315 t8w, &
316 rho
317 !
318 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
319 INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
320 gaer300,gaer400,gaer600,gaer999, & ! jcb
321 waer300,waer400,waer600,waer999, & ! jcb
322 qc_adjust, qi_adjust
324 LOGICAL, OPTIONAL :: cu_rad_feedback
326 INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
328 !
329 ! variables for aerosols (only if running with chemistry)
330 !
331 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
332 INTENT(IN ) :: pm2_5_dry, &
333 pm2_5_water, &
334 pm2_5_dry_ec
335 !
336 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
337 INTENT(INOUT) :: RTHRATEN, &
338 RTHRATENLW, &
339 RTHRATENSW
341 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
342 ! INTENT(INOUT) :: SWUP, &
343 ! SWDN, &
344 ! SWUPCLEAR, &
345 ! SWDNCLEAR, &
346 ! LWUP, &
347 ! LWDN, &
348 ! LWUPCLEAR, &
349 ! LWDNCLEAR
351 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
352 ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
353 ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
354 ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
355 ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
357 ! TOA and surface, upward and downward, total and clear fluxes
358 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
359 SWUPT, SWUPTC, SWDNT, SWDNTC, &
360 SWUPB, SWUPBC, SWDNB, SWDNBC, &
361 LWUPT, LWUPTC, LWDNT, LWDNTC, &
362 LWUPB, LWUPBC, LWDNB, LWDNBC
364 ! Upward and downward, total and clear sky layer fluxes (W m-2)
365 REAL, DIMENSION( ims:ime, kms:kme+2, jms:jme ), &
366 OPTIONAL, INTENT(INOUT) :: &
367 SWUPFLX,SWUPFLXC,SWDNFLX,SWDNFLXC, &
368 LWUPFLX,LWUPFLXC,LWDNFLX,LWDNFLXC
370 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
371 INTENT(INOUT) :: SWCF, &
372 LWCF, &
373 OLR
376 !
377 REAL, DIMENSION( ims:ime, jms:jme ), &
378 INTENT(IN ) :: XLAT, &
379 XLONG, &
380 ALBEDO, &
381 EMISS
382 !
383 REAL, DIMENSION( ims:ime, jms:jme ), &
384 INTENT(INOUT) :: GSW, &
385 GLW
387 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN
388 !
389 REAL, INTENT(IN ) :: GMT,dt, &
390 julian, xtime
391 INTEGER, INTENT(IN ),OPTIONAL :: YR
392 !
393 INTEGER, INTENT(IN ) :: JULDAY, itimestep
394 REAL, INTENT(IN ),OPTIONAL :: CURR_SECS
395 LOGICAL, INTENT(IN ),OPTIONAL :: ADAPT_STEP_FLAG
397 INTEGER,INTENT(IN) :: NPHS
398 REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: &
399 CFRACH, & !Added
400 CFRACL, & !Added
401 CFRACM, & !Added
402 CZMEAN !Added
403 REAL, DIMENSION( ims:ime, jms:jme ), &
404 INTENT(INOUT) :: &
405 RLWTOA, & !Added
406 RSWTOA, & !Added
407 ACFRST, & !Added
408 ACFRCV !Added
410 INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: &
411 NCFRST, & !Added
412 NCFRCV !Added
414 ! Optional, only for Goddard LW and SW
415 REAL, DIMENSION(IMS:IME, JMS:JME, 1:8) :: ERBE_out !extra output for SDSU
416 REAL, DIMENSION(IMS:IME, JMS:JME), OPTIONAL, INTENT(OUT) :: &
417 TLWDN, TLWUP, &
418 SLWDN, SLWUP, &
419 TSWDN, TSWUP, &
420 SSWDN, SSWUP ! for Goddard schemes
422 ! Optional (only used by CAM lw scheme)
424 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,&
425 INTENT(INOUT) :: abstot
426 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,&
427 INTENT(INOUT) :: absnxt
428 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
429 INTENT(INOUT) :: emstot
431 !
432 ! Optional
433 !
434 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
435 OPTIONAL, &
436 INTENT(INOUT) :: CLDFRA
438 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
439 OPTIONAL, &
440 INTENT(IN ) :: &
441 F_ICE_PHY, &
442 F_RAIN_PHY
444 REAL, DIMENSION( ims:ime, jms:jme ), &
445 OPTIONAL, &
446 INTENT(OUT) :: SWDOWNC
447 !
448 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
449 OPTIONAL, &
450 INTENT(INOUT ) :: &
451 pb &
452 ,qv,qc,qr,qi,qs,qg,qndrop
454 LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop
455 !
456 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
457 OPTIONAL, &
458 INTENT(INOUT) :: taucldi,taucldc
460 ! Variables for slope-dependent radiation
462 REAL, OPTIONAL, INTENT(IN) :: dx,dy
463 INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,topo_shading
464 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: ht
465 INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: shadowmask
468 ! LOCAL VAR
470 REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON
471 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS
472 REAL, DIMENSION( ims:ime, jms:jme ) :: coszr
474 REAL :: DECLIN,SOLCON,XXLAT,TLOCTM,XT24, CEN_LAT
475 INTEGER :: i,j,k,its,ite,jts,jte,ij
476 INTEGER :: STEPABS
477 LOGICAL :: gfdl_lw,gfdl_sw
478 LOGICAL :: doabsems
479 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
481 REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, &
482 DJUL,RJUL,ECCFAC
483 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_temp,qc_temp
484 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save
486 REAL :: next_rad_time
487 LOGICAL :: run_param
488 !------------------------------------------------------------------
489 ! solar related variables are added to declaration
490 !-------------------------------------------------
491 REAL, OPTIONAL, INTENT(OUT) :: DECLINX,SOLCONX
492 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZEN
493 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: HRANG
494 !------------------------------------------------------------------
496 if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return
498 ! ra_call_offset = -1 gives old method where radiation may be called just before output
499 ! ra_call_offset = 0 gives new method where radiation may be called just after output
500 ! and is also consistent with removal of offset in new XTIME
501 ! also need to account for stepra=1 which always has zero modulo output
503 !
504 ! Calculate whether or not to run the radiation step.
505 ! If CURR_SECS is passed in, we will calculate based upon that. If it
506 ! is not passed in, we'll do the old method of using the time STEP number
507 !
509 IF ( (itimestep .EQ. 1) .OR. (MOD(itimestep,STEPRA) .EQ. 1 + ra_call_offset) .OR. &
510 (STEPRA .EQ. 1) ) THEN
511 run_param = .TRUE.
512 ELSE
513 run_param = .FALSE.
514 ENDIF
515 IF (PRESENT(adapt_step_flag)) THEN
516 IF ((adapt_step_flag)) THEN
517 IF ( (itimestep .EQ. 1) .OR. (radt .EQ. 0) .OR. &
518 ( CURR_SECS + dt >= ( INT( CURR_SECS / ( radt * 60 ) + 1 ) * radt * 60) ) ) THEN
519 run_param = .TRUE.
520 ELSE
521 run_param = .FALSE.
522 ENDIF
523 ENDIF
524 ENDIF
526 Radiation_step: IF ( run_param ) then
528 ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
529 STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA
530 IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset &
531 .or. STEPABS .eq. 1 ) THEN
532 doabsems = .true.
533 ELSE
534 doabsems = .false.
535 ENDIF
536 IF (PRESENT(adapt_step_flag)) THEN
537 IF ((adapt_step_flag)) THEN
538 IF ( (itimestep .EQ. 1) .OR. (cam_abs_freq_s .EQ. 0) .OR. &
539 ( CURR_SECS + dt >= ( INT( CURR_SECS / ( cam_abs_freq_s ) + 1 ) * cam_abs_freq_s) ) ) THEN
540 doabsems = .true.
541 ELSE
542 doabsems = .false.
543 ENDIF
544 ENDIF
545 ENDIF
547 gfdl_lw = .false.
548 gfdl_sw = .false.
551 ! moved up and out of OMP loop because it only needs to be computed once
552 ! and because it is not entirely thread-safe (XT24, TOLOCTM and XXLAT need
553 ! their thread-privacy) JM 20100217
554 DO ij = 1 , num_tiles
555 its = i_start(ij)
556 ite = i_end(ij)
557 jts = j_start(ij)
558 jte = j_end(ij)
559 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
560 DEGRAD,DPD )
562 IF(PRESENT(declinx).AND.PRESENT(solconx))THEN
563 ! saved to state arrays used in surface driver
564 declinx=declin
565 solconx=solcon
566 ENDIF
568 IF(PRESENT(coszen).AND.PRESENT(hrang))THEN
569 ! state arrays of hrang and coszen used in surface driver
570 XT24=MOD(XTIME+RADT*0.5,1440.)
571 DO j=jts,jte
572 DO i=its,ite
573 TLOCTM=GMT+XT24/60.+XLONG(I,J)/15.
574 HRANG(I,J)=15.*(TLOCTM-12.)*DEGRAD
575 XXLAT=XLAT(I,J)*DEGRAD
576 COSZEN(I,J)=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG(I,J))
577 ENDDO
578 ENDDO
579 ENDIF
580 ENDDO
582 !---------------
583 !$OMP PARALLEL DO &
584 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
586 DO ij = 1 , num_tiles
587 its = i_start(ij)
588 ite = i_end(ij)
589 jts = j_start(ij)
590 jte = j_end(ij)
592 ! initialize data
594 DO j=jts,jte
595 DO i=its,ite
596 GSW(I,J)=0.
597 GLW(I,J)=0.
598 SWDOWN(I,J)=0.
599 GLAT(I,J)=XLAT(I,J)*DEGRAD
600 GLON(I,J)=XLONG(I,J)*DEGRAD
601 ENDDO
602 ENDDO
604 DO j=jts,jte
605 DO k=kts,kte+1
606 DO i=its,ite
607 RTHRATEN(I,K,J)=0.
608 RTHRATENLW(I,K,J)=0.
609 RTHRATENSW(I,K,J)=0.
610 ! SWUP(I,K,J) = 0.0
611 ! SWDN(I,K,J) = 0.0
612 ! SWUPCLEAR(I,K,J) = 0.0
613 ! SWDNCLEAR(I,K,J) = 0.0
614 ! LWUP(I,K,J) = 0.0
615 ! LWDN(I,K,J) = 0.0
616 ! LWUPCLEAR(I,K,J) = 0.0
617 ! LWDNCLEAR(I,K,J) = 0.0
618 CEMISS(I,K,J)=0.0
619 ENDDO
620 ENDDO
621 ENDDO
623 IF ( PRESENT( SWUPFLX ) ) THEN
624 DO j=jts,jte
625 DO k=kts,kte+2
626 DO i=its,ite
627 SWUPFLX(I,K,J) = 0.0
628 SWDNFLX(I,K,J) = 0.0
629 SWUPFLXC(I,K,J) = 0.0
630 SWDNFLXC(I,K,J) = 0.0
631 LWUPFLX(I,K,J) = 0.0
632 LWDNFLX(I,K,J) = 0.0
633 LWUPFLXC(I,K,J) = 0.0
634 LWDNFLXC(I,K,J) = 0.0
635 ENDDO
636 ENDDO
637 ENDDO
638 ENDIF
640 ! temporarily modify hydrometeors (currently only done for GD scheme and WRF-Chem)
641 !
642 IF ( PRESENT( cu_rad_feedback ) ) THEN
643 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
644 DO j=jts,jte
645 DO k=kts,kte
646 DO i=its,ite
647 qc_save(i,k,j) = qc(i,k,j)
648 qc(i,k,j) = qc(i,k,j) + qc_adjust(i,k,j)
649 ENDDO
650 ENDDO
651 ENDDO
652 ENDIF
653 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
654 DO j=jts,jte
655 DO k=kts,kte
656 DO i=its,ite
657 qi_save(i,k,j) = qi(i,k,j)
658 qi(i,k,j) = qi(i,k,j) + qi_adjust(i,k,j)
659 ENDDO
660 ENDDO
661 ENDDO
662 ENDIF
663 ENDIF
666 ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
667 if(PRESENT(qc) .and. PRESENT(F_QC)) then
668 DO j=jts,jte
669 DO k=kts,kte
670 DO i=its,ite
671 qc_temp(I,K,J)=qc(I,K,J)
672 ENDDO
673 ENDDO
674 ENDDO
675 else
676 DO j=jts,jte
677 DO k=kts,kte
678 DO i=its,ite
679 qc_temp(I,K,J)=0.
680 ENDDO
681 ENDDO
682 ENDDO
683 endif
684 if(PRESENT(qr) .and. PRESENT(F_QR)) then
685 DO j=jts,jte
686 DO k=kts,kte
687 DO i=its,ite
688 qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J)
689 ENDDO
690 ENDDO
691 ENDDO
692 endif
694 !---------------
695 ! Calculate constant for short wave radiation
697 lwrad_cldfra_select: SELECT CASE(lw_physics)
699 CASE (GFDLLWSCHEME)
701 !-- Do nothing, since cloud fractions (with partial cloudiness effects)
702 !-- are defined in GFDL LW/SW schemes and do not need to be initialized.
704 CASE (CAMLWSCHEME)
706 IF ( PRESENT ( CLDFRA ) .AND. &
707 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
708 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
710 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
711 F_QV,F_QC,F_QI,F_QS,t,p, &
712 F_ICE_PHY,F_RAIN_PHY, &
713 ids,ide, jds,jde, kds,kde, &
714 ims,ime, jms,jme, kms,kme, &
715 its,ite, jts,jte, kts,kte )
717 ENDIF
719 CASE (RRTMG_LWSCHEME)
721 IF ( PRESENT ( CLDFRA ) .AND. &
722 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
723 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
725 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
726 F_QV,F_QC,F_QI,F_QS,t,p, &
727 F_ICE_PHY,F_RAIN_PHY, &
728 ids,ide, jds,jde, kds,kde, &
729 ims,ime, jms,jme, kms,kme, &
730 its,ite, jts,jte, kts,kte )
732 ENDIF
734 CASE DEFAULT
736 IF ( PRESENT ( CLDFRA ) .AND. &
737 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
738 CALL cal_cldfra(CLDFRA,qc,qi,F_QC,F_QI, &
739 ids,ide, jds,jde, kds,kde, &
740 ims,ime, jms,jme, kms,kme, &
741 its,ite, jts,jte, kts,kte )
742 ENDIF
744 END SELECT lwrad_cldfra_select
746 lwrad_select: SELECT CASE(lw_physics)
749 CASE (RRTMSCHEME)
750 CALL wrf_debug (100, 'CALL rrtm')
752 CALL RRTMLWRAD( &
753 RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS &
754 ,QV3D=QV &
755 ,QC3D=QC &
756 ,QR3D=QR &
757 ,QI3D=QI &
758 ,QS3D=QS &
759 ,QG3D=QG &
760 ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t &
761 ,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G &
762 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
763 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
764 ,ICLOUD=icloud,WARM_RAIN=warm_rain &
765 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
766 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
767 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
768 )
770 ! by Toshihisa Matsui and Jainn Shi 20090623
771 CASE (goddardlwscheme)
773 CALL wrf_debug (100, 'CALL goddard longwave radiation scheme ')
774 IF (itimestep.eq.1) then
775 call wrf_message('running goddard lw radiation')
776 ENDIF
777 CALL goddardrad(sw_or_lw='lw' &
778 ,rthraten=rthraten,gsf=glw,xlat=xlat,xlong=xlong &
779 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
780 ,dz8w=dz8w,rho_phy=rho,emiss=emiss &
781 ,cldfra3d=cldfra &
782 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
783 ,julday=julday,xtime=xtime &
784 ,declin=declin,solcon=solcon &
785 , center_lat = cen_lat &
786 ,radfrq=radt,degrad=degrad &
787 ,taucldi=taucldi,taucldc=taucldc &
788 ,warm_rain=warm_rain &
789 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
790 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
791 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
792 ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban
793 ,qv3d=qv &
794 ,qc3d=qc &
795 ,qr3d=qr &
796 ,qi3d=qi &
797 ,qs3d=qs &
798 ,qg3d=qg &
799 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr &
800 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg &
801 ,erbe_out=erbe_out & !optional
802 )
804 CASE (GFDLLWSCHEME)
806 CALL wrf_debug (100, 'CALL gfdllw')
808 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
809 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
810 PRESENT(qv) .AND. PRESENT(qc) ) THEN
811 IF ( F_QV .AND. F_QC .AND. F_QS) THEN
812 gfdl_lw = .true.
813 CALL ETARA( &
814 DT=dt,XLAND=xland &
815 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
816 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
817 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
818 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
819 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
820 ,HBOTR=hbotr, HTOPR=htopr &
821 ,ALBEDO=albedo,CUPPT=cuppt &
822 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
823 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
824 ,XTIME=xtime,JULIAN=julian &
825 ,JULYR=julyr,JULDAY=julday &
826 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
827 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
828 ,ACFRST=acfrst,NCFRST=ncfrst &
829 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
830 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
831 ,THRATEN=rthraten,THRATENLW=rthratenlw &
832 ,THRATENSW=rthratensw &
833 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
834 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
835 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
836 )
837 ELSE
838 CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
839 ENDIF
840 ELSE
841 CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
842 ENDIF
844 #ifdef HWRF
845 CASE (HWRFLWSCHEME)
847 CALL wrf_debug (100, 'CALL hwrflw')
849 gfdl_lw = .true.
851 CALL HWRFRA(DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, &
852 XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, &
853 QV=qv,QW=qc_temp,QI=Qi, &
854 TSK2D=tsk,GLW=GLW,GSW=GSW, &
855 TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added
856 GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt,&
857 VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified
858 NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified
859 julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, &
860 CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added
861 ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added
862 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
863 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
864 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
867 #endif
869 CASE (CAMLWSCHEME)
871 CALL wrf_debug(100, 'CALL camrad lw')
872 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
873 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
874 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
875 .AND. PRESENT(AEROSOLC_2) ) THEN
876 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
877 dolw=.true.,dosw=.false., &
878 SWUPT=SWUPT,SWUPTC=SWUPTC, &
879 SWDNT=SWDNT,SWDNTC=SWDNTC, &
880 LWUPT=LWUPT,LWUPTC=LWUPTC, &
881 LWDNT=LWDNT,LWDNTC=LWDNTC, &
882 SWUPB=SWUPB,SWUPBC=SWUPBC, &
883 SWDNB=SWDNB,SWDNBC=SWDNBC, &
884 LWUPB=LWUPB,LWUPBC=LWUPBC, &
885 LWDNB=LWDNB,LWDNBC=LWDNBC, &
886 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
887 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
888 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
889 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
890 ,QV3D=qv &
891 ,QC3D=qc &
892 ,QR3D=qr &
893 ,QI3D=qi &
894 ,QS3D=qs &
895 ,QG3D=qg &
896 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
897 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
898 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
899 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
900 dz8w=dz8w, &
901 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
902 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
903 num_months=n_ozmixm, &
904 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
905 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
906 paerlev=paerlev, naer_c=n_aerosolc, &
907 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
908 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
909 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
910 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
911 ,doabsems=doabsems &
912 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
913 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
914 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
915 )
916 ELSE
917 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
918 ENDIF
920 CASE (RRTMG_LWSCHEME)
921 CALL wrf_debug (100, 'CALL rrtmg_lw')
923 CALL RRTMG_LWRAD( &
924 RTHRATENLW=RTHRATEN, &
925 LWUPT=LWUPT,LWUPTC=LWUPTC, &
926 LWDNT=LWDNT,LWDNTC=LWDNTC, &
927 LWUPB=LWUPB,LWUPBC=LWUPBC, &
928 LWDNB=LWDNB,LWDNBC=LWDNBC, &
929 GLW=GLW,OLR=RLWTOA,LWCF=LWCF, &
930 EMISS=EMISS, &
931 P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t, &
932 T8W=t8w,RHO3D=rho,R=R_d,G=G, &
933 ICLOUD=icloud,WARM_RAIN=warm_rain,CLDFRA3D=CLDFRA,&
934 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
935 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
936 QV3D=QV,QC3D=QC,QR3D=QR, &
937 QI3D=QI,QS3D=QS,QG3D=QG, &
938 F_QV=F_QV,F_QC=F_QC,F_QR=F_QR, &
939 F_QI=F_QI,F_QS=F_QS,F_QG=F_QG, &
940 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
941 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
942 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
943 LWUPFLX=LWUPFLX,LWUPFLXC=LWUPFLXC, &
944 LWDNFLX=LWDNFLX,LWDNFLXC=LWDNFLXC &
945 )
947 CASE (HELDSUAREZ)
948 CALL wrf_debug (100, 'CALL heldsuarez')
950 CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, &
951 t8w, rho, R_d,G,CP, dt, xlat, degrad, &
952 ids,ide, jds,jde, kds,kde, &
953 ims,ime, jms,jme, kms,kme, &
954 its,ite, jts,jte, kts,kte )
956 CASE DEFAULT
958 WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
959 CALL wrf_error_fatal ( wrf_err_message )
961 END SELECT lwrad_select
963 IF (lw_physics .gt. 0 .and. .not.gfdl_lw) THEN
964 DO j=jts,jte
965 DO k=kts,kte
966 DO i=its,ite
967 RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
968 ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY)
969 IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
970 ENDDO
971 ENDDO
972 ENDDO
973 ENDIF
975 !JJS 20090623 vvvvv
976 IF (lw_physics .eq. goddardlwscheme) THEN
977 DO j=jts,jte
978 DO i=its,ite
979 tlwdn(i,j) = erbe_out(i,j,1) ! TOA LW downwelling flux [W/m2]
980 tlwup(i,j) = erbe_out(i,j,2) ! TOA LW upwelling flux [W/m2]
981 slwdn(i,j) = erbe_out(i,j,3) ! surface LW downwelling flux [W/m2]
982 slwup(i,j) = erbe_out(i,j,4) ! surface LW upwelling flux [W/m2]
983 ENDDO
984 ENDDO
985 ENDIF
986 !JJS 20090623 ^^^^^
988 !
989 swrad_cldfra_select: SELECT CASE(sw_physics)
991 CASE (CAMSWSCHEME)
993 IF ( PRESENT ( CLDFRA ) .AND. &
994 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
995 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
997 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
998 F_QV,F_QC,F_QI,F_QS,t,p, &
999 F_ICE_PHY,F_RAIN_PHY, &
1000 ids,ide, jds,jde, kds,kde, &
1001 ims,ime, jms,jme, kms,kme, &
1002 its,ite, jts,jte, kts,kte )
1003 ENDIF
1005 CASE (RRTMG_SWSCHEME)
1007 IF ( PRESENT ( CLDFRA ) .AND. &
1008 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
1009 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
1011 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
1012 F_QV,F_QC,F_QI,F_QS,t,p, &
1013 F_ICE_PHY,F_RAIN_PHY, &
1014 ids,ide, jds,jde, kds,kde, &
1015 ims,ime, jms,jme, kms,kme, &
1016 its,ite, jts,jte, kts,kte )
1017 ENDIF
1019 CASE DEFAULT
1021 END SELECT swrad_cldfra_select
1023 swrad_select: SELECT CASE(sw_physics)
1025 CASE (SWRADSCHEME)
1026 CALL wrf_debug(100, 'CALL swrad')
1027 CALL SWRAD( &
1028 DT=dt,RTHRATEN=rthraten,GSW=gsw &
1029 ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
1030 #ifdef WRF_CHEM
1031 ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
1032 ,PM2_5_DRY_EC=pm2_5_dry_ec &
1033 #endif
1034 ,RHO_PHY=rho,T3D=t &
1035 ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
1036 ,R=r_d,CP=cp,G=g,JULDAY=julday &
1037 ,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
1038 ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
1039 ,warm_rain=warm_rain &
1040 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1041 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1042 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1043 ,QV3D=qv &
1044 ,QC3D=qc &
1045 ,QR3D=qr &
1046 ,QI3D=qi &
1047 ,QS3D=qs &
1048 ,QG3D=qg &
1049 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1050 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg )
1052 CASE (GSFCSWSCHEME)
1053 CALL wrf_debug(100, 'CALL gsfcswrad')
1054 CALL GSFCSWRAD( &
1055 RTHRATEN=rthraten,GSW=gsw,XLAT=xlat,XLONG=xlong &
1056 ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
1057 ,DZ8W=dz8w,RHO_PHY=rho &
1058 ,CLDFRA3D=cldfra,RSWTOA=rswtoa &
1059 ,GMT=gmt,CP=cp,G=g &
1060 ,JULDAY=julday,XTIME=xtime &
1061 ,DECLIN=declin,SOLCON=solcon &
1062 ,RADFRQ=radt,DEGRAD=degrad &
1063 ,TAUCLDI=taucldi,TAUCLDC=taucldc &
1064 ,WARM_RAIN=warm_rain &
1065 #ifdef WRF_CHEM
1066 ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
1067 ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
1068 ,GAER300=gaer300,GAER400=gaer400 & ! jcb
1069 ,GAER600=gaer600,GAER999=gaer999 & ! jcb
1070 ,WAER300=waer300,WAER400=waer400 & ! jcb
1071 ,WAER600=waer600,WAER999=waer999 & ! jcb
1072 ,aer_ra_feedback=aer_ra_feedback &
1073 #endif
1074 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1075 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1076 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1077 ,QV3D=qv &
1078 ,QC3D=qc &
1079 ,QR3D=qr &
1080 ,QI3D=qi &
1081 ,QS3D=qs &
1082 ,QG3D=qg &
1083 ,QNDROP3D=qndrop &
1084 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1085 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
1086 ,F_QNDROP=f_qndrop &
1087 )
1089 ! by Toshihisa Matsui and Jainn Shi 20090624
1091 CASE (goddardswscheme)
1093 CALL wrf_debug(100, 'CALL goddard shortwave radiation scheme ')
1094 IF (itimestep.eq.1) then
1095 call wrf_message('running goddard sw radiation')
1096 ENDIF
1097 CALL goddardrad(sw_or_lw='sw' &
1098 ,rthraten=rthraten,gsf=gsw,xlat=xlat,xlong=xlong &
1099 ,alb=albedo,t3d=t,p3d=p,p8w3d=p8w,pi3d=pi &
1100 ,dz8w=dz8w,rho_phy=rho,emiss=emiss &
1101 ,cldfra3d=cldfra &
1102 ,gmt=gmt,cp=cp,g=g,t8w=t8w &
1103 ,julday=julday,xtime=xtime &
1104 ,declin=declin,solcon=solcon &
1105 , center_lat = cen_lat &
1106 ,radfrq=radt,degrad=degrad &
1107 ,taucldi=taucldi,taucldc=taucldc &
1108 ,warm_rain=warm_rain &
1109 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
1110 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
1111 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte &
1112 ! ,cosz_urb2d=cosz_urb2d ,omg_urb2d=omg_urb2d & !urban
1113 ,qv3d=qv &
1114 ,qc3d=qc &
1115 ,qr3d=qr &
1116 ,qi3d=qi &
1117 ,qs3d=qs &
1118 ,qg3d=qg &
1119 ,f_qv=f_qv,f_qc=f_qc,f_qr=f_qr &
1120 ,f_qi=f_qi,f_qs=f_qs,f_qg=f_qg &
1121 ,erbe_out=erbe_out & !optional
1122 )
1124 CASE (CAMSWSCHEME)
1125 CALL wrf_debug(100, 'CALL camrad sw')
1126 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
1127 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
1128 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
1129 .AND. PRESENT(AEROSOLC_2) ) THEN
1130 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
1131 dolw=.false.,dosw=.true., &
1132 SWUPT=SWUPT,SWUPTC=SWUPTC, &
1133 SWDNT=SWDNT,SWDNTC=SWDNTC, &
1134 LWUPT=LWUPT,LWUPTC=LWUPTC, &
1135 LWDNT=LWDNT,LWDNTC=LWDNTC, &
1136 SWUPB=SWUPB,SWUPBC=SWUPBC, &
1137 SWDNB=SWDNB,SWDNBC=SWDNBC, &
1138 LWUPB=LWUPB,LWUPBC=LWUPBC, &
1139 LWDNB=LWDNB,LWDNBC=LWDNBC, &
1140 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
1141 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
1142 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
1143 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
1144 ,QV3D=qv &
1145 ,QC3D=qc &
1146 ,QR3D=qr &
1147 ,QI3D=qi &
1148 ,QS3D=qs &
1149 ,QG3D=qg &
1150 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
1151 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
1152 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
1153 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
1154 dz8w=dz8w, &
1155 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1156 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
1157 num_months=n_ozmixm, &
1158 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
1159 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
1160 paerlev=paerlev, naer_c=n_aerosolc, &
1161 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
1162 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,YR=YR,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
1163 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
1164 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
1165 ,doabsems=doabsems &
1166 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1167 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1168 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1169 )
1170 ELSE
1171 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
1172 ENDIF
1173 DO j=jts,jte
1174 DO k=kts,kte
1175 DO i=its,ite
1176 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
1177 ENDDO
1178 ENDDO
1179 ENDDO
1181 CASE (RRTMG_SWSCHEME)
1182 CALL wrf_debug(100, 'CALL rrtmg_sw')
1183 CALL RRTMG_SWRAD( &
1184 RTHRATENSW=RTHRATENSW, &
1185 SWUPT=SWUPT,SWUPTC=SWUPTC, &
1186 SWDNT=SWDNT,SWDNTC=SWDNTC, &
1187 SWUPB=SWUPB,SWUPBC=SWUPBC, &
1188 SWDNB=SWDNB,SWDNBC=SWDNBC, &
1189 SWCF=SWCF,GSW=GSW, &
1190 XTIME=XTIME,GMT=GMT,XLAT=XLAT,XLONG=XLONG, &
1191 RADT=RADT,DEGRAD=DEGRAD,DECLIN=DECLIN, &
1192 COSZR=COSZR,JULDAY=JULDAY,SOLCON=SOLCON, &
1193 ALBEDO=ALBEDO,t3d=t,t8w=t8w,TSK=TSK, &
1194 p3d=p,p8w=p8w,pi3d=pi,rho3d=rho, &
1195 dz8w=dz8w,CLDFRA3D=CLDFRA,R=R_D,G=G, &
1196 ICLOUD=icloud,WARM_RAIN=warm_rain, &
1197 F_ICE_PHY=F_ICE_PHY,F_RAIN_PHY=F_RAIN_PHY, &
1198 XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
1199 QV3D=qv,QC3D=qc,QR3D=qr, &
1200 QI3D=qi,QS3D=qs,QG3D=qg, &
1201 F_QV=f_qv,F_QC=f_qc,F_QR=f_qr, &
1202 F_QI=f_qi,F_QS=f_qs,F_QG=f_qg, &
1203 IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde,&
1204 IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme,&
1205 ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte,&
1206 SWUPFLX=SWUPFLX,SWUPFLXC=SWUPFLXC, &
1207 SWDNFLX=SWDNFLX,SWDNFLXC=SWDNFLXC &
1208 )
1209 DO j=jts,jte
1210 DO k=kts,kte
1211 DO i=its,ite
1212 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
1213 ENDDO
1214 ENDDO
1215 ENDDO
1217 CASE (GFDLSWSCHEME)
1219 CALL wrf_debug (100, 'CALL gfdlsw')
1221 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
1222 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
1223 PRESENT(qv) .AND. PRESENT(qc) ) THEN
1224 IF ( F_QV .AND. F_QC .AND. F_QS ) THEN
1225 gfdl_sw = .true.
1226 CALL ETARA( &
1227 DT=dt,XLAND=xland &
1228 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
1229 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
1230 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
1231 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
1232 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
1233 ,HBOTR=hbotr, HTOPR=htopr &
1234 ,ALBEDO=albedo,CUPPT=cuppt &
1235 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
1236 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
1237 ,XTIME=xtime,JULIAN=julian &
1238 ,JULYR=julyr,JULDAY=julday &
1239 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
1240 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
1241 ,ACFRST=acfrst,NCFRST=ncfrst &
1242 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
1243 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
1244 ,THRATEN=rthraten,THRATENLW=rthratenlw &
1245 ,THRATENSW=rthratensw &
1246 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1247 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1248 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1249 )
1250 ELSE
1251 CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
1252 ENDIF
1253 ELSE
1254 CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
1255 ENDIF
1257 #ifdef HWRF
1258 CASE (HWRFSWSCHEME)
1260 CALL wrf_debug (100, 'CALL hwrfsw')
1262 gfdl_sw = .true.
1263 CALL HWRFRA(DT=dt,thraten=RTHRATEN,thratenlw=RTHRATENLW,thratensw=RTHRATENSW,pi3d=pi, &
1264 XLAND=xland,P8w=p8w,DZ8w=dz8w,RHO_PHY=rho,P_PHY=p,T=t, &
1265 QV=qv,QW=qc_temp,QI=Qi, &
1266 TSK2D=tsk,GLW=GLW,GSW=GSW, &
1267 TOTSWDN=swdown,TOTLWDN=glw,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean, & !Added
1268 GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot,htopr=htopr,hbotr=hbotr,ALBEDO=albedo,CUPPT=cuppt, &
1269 VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt, & !Modified
1270 NSTEPRA=stepra,NPHS=nphs,itimestep=itimestep, & !Modified
1271 julyr=julyr,julday=julday,gfdl_lw=gfdl_lw,gfdl_sw=gfdl_sw, &
1272 CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach, & !Added
1273 ACFRST=acfrst,NCFRST=ncfrst,ACFRCV=acfrcv,NCFRCV=ncfrcv, & !Added
1274 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
1275 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
1276 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
1277 #endif
1279 CASE (0)
1281 ! Here in case we don't want to call a sw radiation scheme
1282 ! For example, the Held-Suarez idealized test case
1283 IF (lw_physics /= HELDSUAREZ) THEN
1284 WRITE( wrf_err_message , * ) &
1285 'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')'
1286 CALL wrf_error_fatal ( wrf_err_message )
1287 END IF
1289 CASE DEFAULT
1291 WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
1292 CALL wrf_error_fatal ( wrf_err_message )
1294 END SELECT swrad_select
1296 !JJS 20090623 vvvvv
1297 IF (sw_physics .eq. goddardswscheme) THEN
1298 DO j=jts,jte
1299 DO i=its,ite
1300 tswdn(i,j) = erbe_out(i,j,5) ! TOA SW downwelling flux [W/m2]
1301 tswup(i,j) = erbe_out(i,j,6) ! TOA SW upwelling flux [W/m2]
1302 sswdn(i,j) = erbe_out(i,j,7) ! surface SW downwelling flux [W/m2]
1303 sswup(i,j) = erbe_out(i,j,8) ! surface SW upwelling flux [W/m2]
1304 ENDDO
1305 ENDDO
1306 ENDIF
1307 !JJS 20090623 ^^^^^
1309 IF (sw_physics .gt. 0 .and. .not.gfdl_sw) THEN
1310 DO j=jts,jte
1311 DO k=kts,kte
1312 DO i=its,ite
1313 RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
1314 ENDDO
1315 ENDDO
1316 ENDDO
1318 DO j=jts,jte
1319 DO i=its,ite
1320 SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
1321 ENDDO
1322 ENDDO
1324 ENDIF
1326 IF ( PRESENT( cu_rad_feedback ) ) THEN
1327 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
1328 DO j=jts,jte
1329 DO k=kts,kte
1330 DO i=its,ite
1331 qc(i,k,j) = qc_save(i,k,j)
1332 ENDDO
1333 ENDDO
1334 ENDDO
1335 ENDIF
1336 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
1337 DO j=jts,jte
1338 DO k=kts,kte
1339 DO i=its,ite
1340 qi(i,k,j) = qi_save(i,k,j)
1341 ENDDO
1342 ENDDO
1343 ENDDO
1344 ENDIF
1345 ENDIF
1347 ENDDO
1348 !$OMP END PARALLEL DO
1350 ENDIF Radiation_step
1352 accumulate_lw_select: SELECT CASE(lw_physics)
1354 CASE (CAMLWSCHEME,RRTMG_LWSCHEME)
1355 IF(PRESENT(LWUPTC))THEN
1356 !$OMP PARALLEL DO &
1357 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1359 DO ij = 1 , num_tiles
1360 its = i_start(ij)
1361 ite = i_end(ij)
1362 jts = j_start(ij)
1363 jte = j_end(ij)
1365 DO j=jts,jte
1366 DO i=its,ite
1367 ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DT
1368 ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DT
1369 ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DT
1370 ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DT
1371 ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DT
1372 ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DT
1373 ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DT
1374 ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DT
1375 ENDDO
1376 ENDDO
1377 ENDDO
1378 !$OMP END PARALLEL DO
1379 ENDIF
1380 CASE DEFAULT
1381 END SELECT accumulate_lw_select
1383 accumulate_sw_select: SELECT CASE(sw_physics)
1385 CASE (CAMSWSCHEME,RRTMG_SWSCHEME)
1386 IF(PRESENT(SWUPTC))THEN
1387 !$OMP PARALLEL DO &
1388 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1390 DO ij = 1 , num_tiles
1391 its = i_start(ij)
1392 ite = i_end(ij)
1393 jts = j_start(ij)
1394 jte = j_end(ij)
1396 DO j=jts,jte
1397 DO i=its,ite
1398 ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DT
1399 ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DT
1400 ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DT
1401 ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DT
1402 ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DT
1403 ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DT
1404 ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DT
1405 ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DT
1406 ENDDO
1407 ENDDO
1408 ENDDO
1409 !$OMP END PARALLEL DO
1410 ENDIF
1412 CASE DEFAULT
1413 END SELECT accumulate_sw_select
1415 END SUBROUTINE radiation_driver
1417 SUBROUTINE pre_radiation_driver ( grid, config_flags &
1418 ,itimestep, ra_call_offset &
1419 ,XLAT, XLONG, GMT, julian, xtime, RADT, STEPRA &
1420 ,ht,dx,dy,sina,cosa,shadowmask,slope_rad ,topo_shading &
1421 ,shadlen,ht_shad,ht_loc &
1422 ,ht_shad_bxs, ht_shad_bxe &
1423 ,ht_shad_bys, ht_shad_bye &
1424 ,nested, min_ptchsz &
1425 ,spec_bdy_width &
1426 ,ids, ide, jds, jde, kds, kde &
1427 ,ims, ime, jms, jme, kms, kme &
1428 ,ips, ipe, jps, jpe, kps, kpe &
1429 ,i_start, i_end &
1430 ,j_start, j_end &
1431 ,kts, kte &
1432 ,num_tiles )
1434 USE module_domain , ONLY : domain
1435 #ifdef DM_PARALLEL
1436 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks,wrf_dm_minval_integer
1437 # if (EM_CORE == 1)
1438 USE module_comm_dm , ONLY : halo_toposhad_sub
1439 # endif
1440 #endif
1441 USE module_bc
1442 USE module_model_constants
1444 IMPLICIT NONE
1446 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1447 ims,ime, jms,jme, kms,kme, &
1448 ips,ipe, jps,jpe, kps,kpe, &
1449 kts,kte, &
1450 num_tiles
1452 TYPE(domain) , INTENT(INOUT) :: grid
1453 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1455 INTEGER, INTENT(IN ) :: itimestep, ra_call_offset, stepra, &
1456 slope_rad, topo_shading, &
1457 spec_bdy_width
1459 INTEGER, INTENT(INOUT) :: min_ptchsz
1461 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
1462 i_start,i_end,j_start,j_end
1464 REAL, INTENT(IN ) :: GMT, radt, julian, xtime, dx, dy, shadlen
1466 REAL, DIMENSION( ims:ime, jms:jme ), &
1467 INTENT(IN ) :: XLAT, &
1468 XLONG, &
1469 HT, &
1470 SINA, &
1471 COSA
1473 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
1475 REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width ), &
1476 INTENT(IN ) :: ht_shad_bxs, ht_shad_bxe
1477 REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width ), &
1478 INTENT(IN ) :: ht_shad_bys, ht_shad_bye
1480 INTEGER, DIMENSION( ims:ime, jms:jme ), &
1481 INTENT(INOUT) :: shadowmask
1483 LOGICAL, INTENT(IN ) :: nested
1485 !Local
1486 ! For orographic shading
1487 INTEGER :: niter,ni,psx,psy,idum,jdum,i,j,ij
1488 REAL :: DECLIN,SOLCON
1490 ! Determine minimum patch size for slope-dependent radiation
1491 if (itimestep .eq. 1) then
1492 psx = ipe-ips+1
1493 psy = jpe-jps+1
1494 min_ptchsz = min(psx,psy)
1495 idum = 0
1496 jdum = 0
1497 endif
1499 # ifdef DM_PARALLEL
1500 if (itimestep .eq. 1) then
1501 call wrf_dm_minval_integer (psx,idum,jdum)
1502 call wrf_dm_minval_integer (psy,idum,jdum)
1503 min_ptchsz = min(psx,psy)
1504 endif
1505 # endif
1507 ! Topographic shading
1509 if ((topo_shading.eq.1).and.(itimestep .eq. 1 .or. &
1510 mod(itimestep,STEPRA) .eq. 1 + ra_call_offset)) then
1512 !---------------
1513 ! Calculate constants for short wave radiation
1515 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,DEGRAD,DPD)
1517 ! Make a local copy of terrain height field
1518 do j=jms,jme
1519 do i=ims,ime
1520 ht_loc(i,j) = ht(i,j)
1521 enddo
1522 enddo
1523 ! Determine if iterations are necessary for shadows to propagate from one patch to another
1524 if ((ids.eq.ips).and.(ide.eq.ipe).and.(jds.eq.jps).and.(jde.eq.jpe)) then
1525 niter = 1
1526 else
1527 niter = int(shadlen/(dx*min_ptchsz)+3)
1528 endif
1532 IF( nested ) THEN
1534 !$OMP PARALLEL DO &
1535 !$OMP PRIVATE ( ij )
1537 DO ij = 1 , num_tiles
1539 CALL spec_bdyfield(ht_shad, &
1540 ht_shad_bxs, ht_shad_bxe, &
1541 ht_shad_bys, ht_shad_bye, &
1542 'm', config_flags, spec_bdy_width, 2,&
1543 ids,ide, jds,jde, 1 ,1 , & ! domain dims
1544 ims,ime, jms,jme, 1 ,1 , & ! memory dims
1545 ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
1546 i_start(ij), i_end(ij), &
1547 j_start(ij), j_end(ij), &
1548 1 , 1 )
1549 ENDDO
1550 ENDIF
1552 do ni = 1, niter
1554 !$OMP PARALLEL DO &
1555 !$OMP PRIVATE ( ij,i,j )
1556 do ij = 1 , num_tiles
1558 call toposhad_init (ht_shad,ht_loc, &
1559 shadowmask,nested,ni, &
1560 ids,ide, jds,jde, kds,kde, &
1561 ims,ime, jms,jme, kms,kme, &
1562 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1563 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1564 min(j_end(ij), jde-1), kts, kte )
1566 enddo
1567 !$OMP END PARALLEL DO
1570 !$OMP PARALLEL DO &
1571 !$OMP PRIVATE ( ij,i,j )
1572 do ij = 1 , num_tiles
1574 call toposhad (xlat,xlong,sina,cosa,xtime,gmt,radt,declin, &
1575 dx,dy,ht_shad,ht_loc,ni, &
1576 shadowmask,shadlen, &
1577 ids,ide, jds,jde, kds,kde, &
1578 ims,ime, jms,jme, kms,kme, &
1579 ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1580 i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1581 min(j_end(ij), jde-1), kts, kte )
1583 enddo
1584 !$OMP END PARALLEL DO
1586 #if defined( DM_PARALLEL ) && (EM_CORE == 1)
1587 # include "HALO_TOPOSHAD.inc"
1588 #endif
1589 enddo
1590 endif
1592 END SUBROUTINE pre_radiation_driver
1594 !---------------------------------------------------------------------
1595 !BOP
1596 ! !IROUTINE: radconst - compute radiation terms
1597 ! !INTERFAC:
1598 SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
1599 DEGRAD,DPD )
1600 !---------------------------------------------------------------------
1601 USE module_wrf_error
1602 IMPLICIT NONE
1603 !---------------------------------------------------------------------
1605 ! !ARGUMENTS:
1606 REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
1607 REAL, INTENT(OUT ) :: DECLIN,SOLCON
1608 REAL :: OBECL,SINOB,SXLONG,ARG, &
1609 DECDEG,DJUL,RJUL,ECCFAC
1610 !
1611 ! !DESCRIPTION:
1612 ! Compute terms used in radiation physics
1613 !EOP
1615 ! for short wave radiation
1617 DECLIN=0.
1618 SOLCON=0.
1620 !-----OBECL : OBLIQUITY = 23.5 DEGREE.
1622 OBECL=23.5*DEGRAD
1623 SINOB=SIN(OBECL)
1625 !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
1627 IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
1628 IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
1629 SXLONG=SXLONG*DEGRAD
1630 ARG=SINOB*SIN(SXLONG)
1631 DECLIN=ASIN(ARG)
1632 DECDEG=DECLIN/DEGRAD
1633 !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
1634 DJUL=JULIAN*360./365.
1635 RJUL=DJUL*DEGRAD
1636 ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* &
1637 COS(2*RJUL)+0.000077*SIN(2*RJUL)
1638 SOLCON=1370.*ECCFAC
1640 END SUBROUTINE radconst
1642 !---------------------------------------------------------------------
1643 !BOP
1644 ! !IROUTINE: cal_cldfra - Compute cloud fraction
1645 ! !INTERFACE:
1646 SUBROUTINE cal_cldfra(CLDFRA,QC,QI,F_QC,F_QI, &
1647 ids,ide, jds,jde, kds,kde, &
1648 ims,ime, jms,jme, kms,kme, &
1649 its,ite, jts,jte, kts,kte )
1650 !---------------------------------------------------------------------
1651 IMPLICIT NONE
1652 !---------------------------------------------------------------------
1653 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1654 ims,ime, jms,jme, kms,kme, &
1655 its,ite, jts,jte, kts,kte
1657 !
1658 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1659 CLDFRA
1661 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1662 QI, &
1663 QC
1665 LOGICAL,INTENT(IN) :: F_QC,F_QI
1667 REAL thresh
1668 INTEGER:: i,j,k
1669 ! !DESCRIPTION:
1670 ! Compute cloud fraction from input ice and cloud water fields
1671 ! if provided.
1672 !
1673 ! Whether QI or QC is active or not is determined from the indices of
1674 ! the fields into the 4D scalar arrays in WRF. These indices are
1675 ! P_QI and P_QC, respectively, and they are passed in to the routine
1676 ! to enable testing to see if QI and QC represent active fields in
1677 ! the moisture 4D scalar array carried by WRF.
1678 !
1679 ! If a field is active its index will have a value greater than or
1680 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1681 ! this routine.
1682 !EOP
1683 !---------------------------------------------------------------------
1684 thresh=1.0e-6
1686 IF ( f_qi .AND. f_qc ) THEN
1687 DO j = jts,jte
1688 DO k = kts,kte
1689 DO i = its,ite
1690 IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
1691 CLDFRA(i,k,j)=1.
1692 ELSE
1693 CLDFRA(i,k,j)=0.
1694 ENDIF
1695 ENDDO
1696 ENDDO
1697 ENDDO
1698 ELSE IF ( f_qc ) THEN
1699 DO j = jts,jte
1700 DO k = kts,kte
1701 DO i = its,ite
1702 IF ( QC(i,k,j) .gt. thresh) THEN
1703 CLDFRA(i,k,j)=1.
1704 ELSE
1705 CLDFRA(i,k,j)=0.
1706 ENDIF
1707 ENDDO
1708 ENDDO
1709 ENDDO
1710 ELSE
1711 DO j = jts,jte
1712 DO k = kts,kte
1713 DO i = its,ite
1714 CLDFRA(i,k,j)=0.
1715 ENDDO
1716 ENDDO
1717 ENDDO
1718 ENDIF
1720 END SUBROUTINE cal_cldfra
1722 !BOP
1723 ! !IROUTINE: cal_cldfra2 - Compute cloud fraction
1724 ! !INTERFACE:
1725 ! cal_cldfra_xr - Compute cloud fraction.
1726 ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
1727 !!
1728 !!--- Cloud fraction parameterization follows Randall, 1994
1729 !! (see Hong et al., 1998)
1730 !! (modified by Ferrier, Feb '02)
1731 !
1732 SUBROUTINE cal_cldfra2(CLDFRA, QV, QC, QI, QS, &
1733 F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
1734 F_ICE_PHY,F_RAIN_PHY, &
1735 ids,ide, jds,jde, kds,kde, &
1736 ims,ime, jms,jme, kms,kme, &
1737 its,ite, jts,jte, kts,kte )
1738 !---------------------------------------------------------------------
1739 IMPLICIT NONE
1740 !---------------------------------------------------------------------
1741 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1742 ims,ime, jms,jme, kms,kme, &
1743 its,ite, jts,jte, kts,kte
1745 !
1746 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1747 CLDFRA
1749 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1750 QV, &
1751 QI, &
1752 QC, &
1753 QS, &
1754 t_phy, &
1755 p_phy
1756 ! p_phy, &
1757 ! F_ICE_PHY, &
1758 ! F_RAIN_PHY
1760 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1761 OPTIONAL, &
1762 INTENT(IN ) :: &
1763 F_ICE_PHY, &
1764 F_RAIN_PHY
1766 LOGICAL,OPTIONAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
1768 ! REAL thresh
1769 INTEGER:: i,j,k
1770 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
1772 REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
1773 PEXP=0.25, RHGRID=1.0
1774 REAL , PARAMETER :: SVP1=0.61078
1775 REAL , PARAMETER :: SVP2=17.2693882
1776 REAL , PARAMETER :: SVPI2=21.8745584
1777 REAL , PARAMETER :: SVP3=35.86
1778 REAL , PARAMETER :: SVPI3=7.66
1779 REAL , PARAMETER :: SVPT0=273.15
1780 REAL , PARAMETER :: r_d = 287.
1781 REAL , PARAMETER :: r_v = 461.6
1782 REAL , PARAMETER :: ep_2=r_d/r_v
1783 ! !DESCRIPTION:
1784 ! Compute cloud fraction from input ice and cloud water fields
1785 ! if provided.
1786 !
1787 ! Whether QI or QC is active or not is determined from the indices of
1788 ! the fields into the 4D scalar arrays in WRF. These indices are
1789 ! P_QI and P_QC, respectively, and they are passed in to the routine
1790 ! to enable testing to see if QI and QC represent active fields in
1791 ! the moisture 4D scalar array carried by WRF.
1792 !
1793 ! If a field is active its index will have a value greater than or
1794 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1795 ! this routine.
1796 !EOP
1799 !-----------------------------------------------------------------------
1800 !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
1801 ! (modified by Ferrier, Feb '02)
1802 !
1803 !--- Cloud fraction parameterization follows Randall, 1994
1804 ! (see Hong et al., 1998)
1805 !-----------------------------------------------------------------------
1806 ! Note: ep_2=287./461.6 Rd/Rv
1807 ! Note: R_D=287.
1809 ! Alternative calculation for critical RH for grid saturation
1810 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
1812 ! Calculate saturation mixing ratio weighted according to the fractions of
1813 ! water and ice.
1814 ! Following:
1815 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
1816 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
1817 !
1818 ! over ice over water
1819 ! a = 21.8745584 17.2693882
1820 ! b = 7.66 35.86
1822 !---------------------------------------------------------------------
1824 DO j = jts,jte
1825 DO k = kts,kte
1826 DO i = its,ite
1827 tc = t_phy(i,k,j) - SVPT0
1828 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
1829 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
1830 QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
1831 QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
1833 IF ( PRESENT(F_QI) .and. PRESENT(F_QC) .and. PRESENT(F_QS) ) THEN
1835 ! mji - For MP options 2, 4, 6, 7, 8, etc. (qc = liquid, qi = ice, qs = snow)
1836 IF ( F_QI .and. F_QC .and. F_QS) THEN
1837 QCLD = QI(i,k,j)+QC(i,k,j)+QS(i,k,j)
1838 IF (QCLD .LT. QCLDMIN) THEN
1839 weight = 0.
1840 ELSE
1841 weight = (QI(i,k,j)+QS(i,k,j)) / QCLD
1842 ENDIF
1843 ENDIF
1845 ! mji - For MP options 1 and 3, (qc only)
1846 ! For MP=1, qc = liquid, for MP=3, qc = liquid or ice depending on temperature
1847 IF ( F_QC .and. .not. F_QI .and. .not. F_QS ) THEN
1848 QCLD = QC(i,k,j)
1849 IF (QCLD .LT. QCLDMIN) THEN
1850 weight = 0.
1851 ELSE
1852 if (t_phy(i,k,j) .gt. 273.15) weight = 0.
1853 if (t_phy(i,k,j) .le. 273.15) weight = 1.
1854 ENDIF
1855 ENDIF
1857 ! mji - For MP option 5; (qc = liquid, qs = ice)
1858 IF ( F_QC .and. .not. F_QI .and. F_QS .and. PRESENT(F_ICE_PHY) ) THEN
1860 ! Mixing ratios of cloud water & total ice (cloud ice + snow).
1861 ! Mixing ratios of rain are not considered in this scheme.
1862 ! F_ICE is fraction of ice
1863 ! F_RAIN is fraction of rain
1865 QIMID = QS(i,k,j)
1866 QWMID = QC(i,k,j)
1867 ! old method
1868 ! QIMID = QC(i,k,j)*F_ICE_PHY(i,k,j)
1869 ! QWMID = (QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
1870 !
1871 !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
1872 ! only cloud water + cloud ice + snow
1873 !
1874 QCLD=QWMID+QIMID
1875 IF (QCLD .LT. QCLDMIN) THEN
1876 weight = 0.
1877 ELSE
1878 weight = F_ICE_PHY(i,k,j)
1879 ENDIF
1880 ENDIF
1882 ELSE
1883 CLDFRA(i,k,j)=0.
1885 ENDIF ! IF ( F_QI .and. F_QC .and. F_QS)
1888 QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI
1889 RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
1890 !
1891 !--- Determine cloud fraction (modified from original algorithm)
1892 !
1893 IF (QCLD .LT. QCLDMIN) THEN
1894 !
1895 !--- Assume zero cloud fraction if there is no cloud mixing ratio
1896 !
1897 CLDFRA(i,k,j)=0.
1898 ELSEIF(RHUM.GE.RHGRID)THEN
1899 !
1900 !--- Assume cloud fraction of unity if near saturation and the cloud
1901 ! mixing ratio is at or above the minimum threshold
1902 !
1903 CLDFRA(i,k,j)=1.
1904 ELSE
1905 !
1906 !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
1907 ! modified based on assumed grid-scale saturation at RH=RHgrid.
1908 !
1909 SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
1910 DENOM=(SUBSAT)**GAMMA
1911 ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
1912 ! prevent negative values (new)
1913 RHUM=MAX(1.E-10, RHUM)
1914 CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG))
1915 !! ARG=-1000*QCLD/(RHUM-RHGRID)
1916 !! ARG=MAX(ARG, ARGMIN)
1917 !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
1918 IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
1919 ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
1920 ENDDO !--- End DO i
1921 ENDDO !--- End DO k
1922 ENDDO !--- End DO j
1924 END SUBROUTINE cal_cldfra2
1927 SUBROUTINE toposhad_init(ht_shad,ht_loc,shadowmask,nested,iter, &
1928 ids,ide, jds,jde, kds,kde, &
1929 ims,ime, jms,jme, kms,kme, &
1930 ips,ipe, jps,jpe, kps,kpe, &
1931 its,ite, jts,jte, kts,kte )
1933 USE module_model_constants
1935 implicit none
1937 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
1938 ims,ime, jms,jme, kms,kme, &
1939 ips,ipe, jps,jpe, kps,kpe, &
1940 its,ite, jts,jte, kts,kte
1942 LOGICAL, INTENT(IN) :: nested
1944 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad, ht_loc
1946 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
1947 INTEGER, INTENT(IN) :: iter
1949 ! Local variables
1951 INTEGER :: i, j
1953 if (iter.eq.1) then
1955 ! Initialize shadow mask
1956 do j=jts,jte
1957 do i=its,ite
1958 shadowmask(i,j) = 0
1959 ENDDO
1960 ENDDO
1962 ! Initialize shading height
1964 IF ( nested ) THEN ! Do not overwrite input from parent domain
1965 do j=max(jts,jds+2),min(jte,jde-3)
1966 do i=max(its,ids+2),min(ite,ide-3)
1967 ht_shad(i,j) = ht_loc(i,j)-0.001
1968 ENDDO
1969 ENDDO
1970 ELSE
1971 do j=jts,jte
1972 do i=its,ite
1973 ht_shad(i,j) = ht_loc(i,j)-0.001
1974 ENDDO
1975 ENDDO
1976 ENDIF
1978 IF ( nested ) THEN ! Check if a shadow exceeding the topography height is available at the lateral domain edge from nesting
1979 if (its.eq.ids) then
1980 do j=jts,jte
1981 if (ht_shad(its,j) .gt. ht_loc(its,j)) then
1982 shadowmask(its,j) = 1
1983 ht_loc(its,j) = ht_shad(its,j)
1984 endif
1985 if (ht_shad(its+1,j) .gt. ht_loc(its+1,j)) then
1986 shadowmask(its+1,j) = 1
1987 ht_loc(its+1,j) = ht_shad(its+1,j)
1988 endif
1989 enddo
1990 endif
1991 if (ite.eq.ide-1) then
1992 do j=jts,jte
1993 if (ht_shad(ite,j) .gt. ht_loc(ite,j)) then
1994 shadowmask(ite,j) = 1
1995 ht_loc(ite,j) = ht_shad(ite,j)
1996 endif
1997 if (ht_shad(ite-1,j) .gt. ht_loc(ite-1,j)) then
1998 shadowmask(ite-1,j) = 1
1999 ht_loc(ite-1,j) = ht_shad(ite-1,j)
2000 endif
2001 enddo
2002 endif
2003 if (jts.eq.jds) then
2004 do i=its,ite
2005 if (ht_shad(i,jts) .gt. ht_loc(i,jts)) then
2006 shadowmask(i,jts) = 1
2007 ht_loc(i,jts) = ht_shad(i,jts)
2008 endif
2009 if (ht_shad(i,jts+1) .gt. ht_loc(i,jts+1)) then
2010 shadowmask(i,jts+1) = 1
2011 ht_loc(i,jts+1) = ht_shad(i,jts+1)
2012 endif
2013 enddo
2014 endif
2015 if (jte.eq.jde-1) then
2016 do i=its,ite
2017 if (ht_shad(i,jte) .gt. ht_loc(i,jte)) then
2018 shadowmask(i,jte) = 1
2019 ht_loc(i,jte) = ht_shad(i,jte)
2020 endif
2021 if (ht_shad(i,jte-1) .gt. ht_loc(i,jte-1)) then
2022 shadowmask(i,jte-1) = 1
2023 ht_loc(i,jte-1) = ht_shad(i,jte-1)
2024 endif
2025 enddo
2026 endif
2027 ENDIF
2029 else
2031 ! Fill the local topography field at the points next to internal tile boundaries with ht_shad values
2032 ! A 2-pt halo has been applied to the ht_shad before the repeated call of this subroutine
2034 if ((its.ne.ids).and.(its.eq.ips)) then
2035 do j=jts-2,jte+2
2036 ht_loc(its-1,j) = max(ht_loc(its-1,j),ht_shad(its-1,j))
2037 ht_loc(its-2,j) = max(ht_loc(its-2,j),ht_shad(its-2,j))
2038 enddo
2039 endif
2040 if ((ite.ne.ide-1).and.(ite.eq.ipe)) then
2041 do j=jts-2,jte+2
2042 ht_loc(ite+1,j) = max(ht_loc(ite+1,j),ht_shad(ite+1,j))
2043 ht_loc(ite+2,j) = max(ht_loc(ite+2,j),ht_shad(ite+2,j))
2044 enddo
2045 endif
2046 if ((jts.ne.jds).and.(jts.eq.jps)) then
2047 do i=its-2,ite+2
2048 ht_loc(i,jts-1) = max(ht_loc(i,jts-1),ht_shad(i,jts-1))
2049 ht_loc(i,jts-2) = max(ht_loc(i,jts-2),ht_shad(i,jts-2))
2050 enddo
2051 endif
2052 if ((jte.ne.jde-1).and.(jte.eq.jpe)) then
2053 do i=its-2,ite+2
2054 ht_loc(i,jte+1) = max(ht_loc(i,jte+1),ht_shad(i,jte+1))
2055 ht_loc(i,jte+2) = max(ht_loc(i,jte+2),ht_shad(i,jte+2))
2056 enddo
2057 endif
2059 endif
2061 END SUBROUTINE toposhad_init
2066 SUBROUTINE toposhad(xlat,xlong,sina,cosa,xtime,gmt,radfrq,declin, &
2067 dx,dy,ht_shad,ht_loc,iter, &
2068 shadowmask,shadlen, &
2069 ids,ide, jds,jde, kds,kde, &
2070 ims,ime, jms,jme, kms,kme, &
2071 ips,ipe, jps,jpe, kps,kpe, &
2072 its,ite, jts,jte, kts,kte )
2075 USE module_model_constants
2077 implicit none
2079 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
2080 ims,ime, jms,jme, kms,kme, &
2081 ips,ipe, jps,jpe, kps,kpe, &
2082 its,ite, jts,jte, kts,kte
2084 INTEGER, INTENT(IN) :: iter
2086 REAL, INTENT(IN) :: RADFRQ,XTIME,DECLIN,dx,dy,gmt,shadlen
2088 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT, XLONG, sina, cosa
2090 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
2092 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
2094 ! Local variables
2096 REAL :: pi, xt24, wgt, ri, rj, argu, sol_azi, topoelev, dxabs, tloctm, hrang, xxlat, csza
2097 INTEGER :: gpshad, ii, jj, i1, i2, j1, j2, i, j
2101 XT24=MOD(XTIME+RADFRQ*0.5,1440.)
2102 pi = 4.*atan(1.)
2103 gpshad = int(shadlen/dx+1.)
2105 if (iter.eq.1) then
2108 j_loop1: DO J=jts,jte
2109 i_loop1: DO I=its,ite
2111 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
2112 HRANG=15.*(TLOCTM-12.)*DEGRAD
2113 XXLAT=XLAT(i,j)*DEGRAD
2114 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
2116 if (csza.lt.1.e-2) then ! shadow mask does not need to be computed
2117 shadowmask(i,j) = 0
2118 ht_shad(i,j) = ht_loc(i,j)-0.001
2119 goto 120
2120 endif
2122 ! Solar azimuth angle
2124 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
2125 if (argu.gt.1) argu = 1
2126 if (argu.lt.-1) argu = -1
2127 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
2128 if (cosa(i,j).ge.0) then
2129 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
2130 else
2131 sol_azi = sol_azi + pi - asin(sina(i,j))
2132 endif
2134 ! Scan for higher surrounding topography
2136 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
2138 do jj = j+1,j+gpshad
2139 ri = i + (jj-j)*tan(sol_azi)
2140 i1 = int(ri)
2141 i2 = i1+1
2142 wgt = ri-i1
2143 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2144 if ((jj.ge.jpe+1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
2145 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2146 goto 120
2147 endif
2148 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2149 if (sin(topoelev).ge.csza) then
2150 shadowmask(i,j) = 1
2151 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2152 endif
2153 enddo
2155 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
2156 do ii = i+1,i+gpshad
2157 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2158 j1 = int(rj)
2159 j2 = j1+1
2160 wgt = rj-j1
2161 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2162 if ((ii.ge.ipe+1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
2163 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2164 goto 120
2165 endif
2166 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2167 if (sin(topoelev).ge.csza) then
2168 shadowmask(i,j) = 1
2169 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2170 endif
2171 enddo
2173 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
2174 do jj = j-1,j-gpshad,-1
2175 ri = i + (jj-j)*tan(sol_azi)
2176 i1 = int(ri)
2177 i2 = i1+1
2178 wgt = ri-i1
2179 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2180 if ((jj.le.jps-1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
2181 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2182 goto 120
2183 endif
2184 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2185 if (sin(topoelev).ge.csza) then
2186 shadowmask(i,j) = 1
2187 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2188 endif
2189 enddo
2191 else ! sun is in the western quarter
2192 do ii = i-1,i-gpshad,-1
2193 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2194 j1 = int(rj)
2195 j2 = j1+1
2196 wgt = rj-j1
2197 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2198 if ((ii.le.ips-1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
2199 if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
2200 goto 120
2201 endif
2202 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2203 if (sin(topoelev).ge.csza) then
2204 shadowmask(i,j) = 1
2205 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2206 endif
2207 enddo
2208 endif
2210 120 continue
2212 ENDDO i_loop1
2213 ENDDO j_loop1
2215 else ! iteration > 1
2218 j_loop2: DO J=jts,jte
2219 i_loop2: DO I=its,ite
2221 ! if (shadowmask(i,j).eq.-1) then ! this indicates that the search ended at a lateral boundary during iteration 1
2223 TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
2224 HRANG=15.*(TLOCTM-12.)*DEGRAD
2225 XXLAT=XLAT(i,j)*DEGRAD
2226 CSZA=SIN(XXLAT)*SIN(DECLIN)+COS(XXLAT)*COS(DECLIN)*COS(HRANG)
2228 ! Solar azimuth angle
2230 argu=(csza*sin(XXLAT)-sin(DECLIN))/(sin(acos(csza))*cos(XXLAT))
2231 if (argu.gt.1) argu = 1
2232 if (argu.lt.-1) argu = -1
2233 sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
2234 if (cosa(i,j).ge.0) then
2235 sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
2236 else
2237 sol_azi = sol_azi + pi - asin(sina(i,j))
2238 endif
2240 ! Scan for higher surrounding topography
2242 if ((sol_azi.gt.1.75*pi).or.(sol_azi.lt.0.25*pi)) then ! sun is in the northern quarter
2244 do jj = j+1,j+gpshad
2245 ri = i + (jj-j)*tan(sol_azi)
2246 i1 = int(ri)
2247 i2 = i1+1
2248 wgt = ri-i1
2249 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2250 if ((jj.ge.min(jde,jpe+3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220
2251 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2252 if (sin(topoelev).ge.csza) then
2253 shadowmask(i,j) = 1
2254 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2255 endif
2256 enddo
2258 else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
2259 do ii = i+1,i+gpshad
2260 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2261 j1 = int(rj)
2262 j2 = j1+1
2263 wgt = rj-j1
2264 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2265 if ((ii.ge.min(ide,ipe+3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220
2266 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2267 if (sin(topoelev).ge.csza) then
2268 shadowmask(i,j) = 1
2269 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2270 endif
2271 enddo
2273 else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
2274 do jj = j-1,j-gpshad,-1
2275 ri = i + (jj-j)*tan(sol_azi)
2276 i1 = int(ri)
2277 i2 = i1+1
2278 wgt = ri-i1
2279 dxabs = sqrt((dy*(jj-j))**2+(dx*(ri-i))**2)
2280 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
2281 topoelev=atan((wgt*ht_loc(i2,jj)+(1.-wgt)*ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
2282 if (sin(topoelev).ge.csza) then
2283 shadowmask(i,j) = 1
2284 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2285 endif
2286 enddo
2288 else ! sun is in the western quarter
2289 do ii = i-1,i-gpshad,-1
2290 rj = j - (ii-i)*tan(pi/2.+sol_azi)
2291 j1 = int(rj)
2292 j2 = j1+1
2293 wgt = rj-j1
2294 dxabs = sqrt((dx*(ii-i))**2+(dy*(rj-j))**2)
2295 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
2296 topoelev=atan((wgt*ht_loc(ii,j2)+(1.-wgt)*ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
2297 if (sin(topoelev).ge.csza) then
2298 shadowmask(i,j) = 1
2299 ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
2300 endif
2301 enddo
2302 endif
2304 220 continue
2305 ! endif
2307 ENDDO i_loop2
2308 ENDDO j_loop2
2310 endif ! iteration
2312 END SUBROUTINE toposhad
2314 END MODULE module_radiation_driver