module_radiation_driver.F
References to this file elsewhere.
1 !WRF:MEDIATION_LAYER:PHYSICS
2 !
3 MODULE module_radiation_driver
4 CONTAINS
5 !BOP
6 ! !IROUTINE: radiation_driver - interface to radiation physics options
7
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 ,GLW, GSW, SWDOWN, XLAT, XLONG, ALBEDO &
22 ,EMISS, rho, p8w, p , pi , dz8w ,t, t8w, GMT &
23 ,XLAND, XICE, TSK, HTOP,HBOT,HTOPR,HBOTR, CUPPT, VEGFRA, SNOW &
24 ,julyr, JULDAY, julian, xtime, RADT, STEPRA, ICLOUD, warm_rain &
25 ,declin_urb,COSZ_URB2D, omg_urb2d & !Optional urban
26 ,ra_call_offset,RSWTOA,RLWTOA, CZMEAN &
27 ,CFRACL, CFRACM, CFRACH &
28 ,ACFRST,NCFRST,ACFRCV,NCFRCV,SWDOWNC &
29 ,z &
30 ,levsiz, n_ozmixm, n_aerosolc, paerlev &
31 ,cam_abs_dim1, cam_abs_dim2, cam_abs_freq_s &
32 ,ozmixm,pin & ! Optional
33 ,m_ps_1,m_ps_2,aerosolc_1,aerosolc_2,m_hybi0 & ! Optional
34 ,abstot, absnxt, emstot & ! Optional
35 ,taucldi, taucldc & ! Optional
36 ,ids, ide, jds, jde, kds, kde &
37 ,ims, ime, jms, jme, kms, kme &
38 ,i_start, i_end &
39 ,j_start, j_end &
40 ,kts, kte &
41 ,num_tiles &
42 ,qv,qc,qr,qi,qs,qg,qndrop &
43 ,f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop &
44 ,CLDFRA ,Pb &
45 ,f_ice_phy,f_rain_phy &
46 ,pm2_5_dry, pm2_5_water, pm2_5_dry_ec &
47 ,tauaer300, tauaer400, tauaer600, tauaer999 & ! jcb
48 ,gaer300, gaer400, gaer600, gaer999 & ! jcb
49 ,waer300, waer400, waer600, waer999 & ! jcb
50 ,qc_adjust ,qi_adjust & ! jm
51 ,cu_rad_feedback, aer_ra_feedback & ! jm
52
53 )
54
55 !-------------------------------------------------------------------------
56
57 ! !USES:
58 USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME &
59 ,SWRADSCHEME, GSFCSWSCHEME &
60 ,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME &
61 ,HELDSUAREZ
62 USE module_model_constants
63 USE module_wrf_error
64
65 ! *** add new modules of schemes here
66
67 USE module_ra_sw
68 USE module_ra_gsfcsw
69 USE module_ra_rrtm
70 USE module_ra_cam
71 USE module_ra_gfdleta
72 USE module_ra_hs
73
74 ! This driver calls subroutines for the radiation parameterizations.
75 !
76 ! short wave radiation choices:
77 ! 1. swrad (19??)
78 !
79 ! long wave radiation choices:
80 ! 1. rrtmlwrad
81 !
82 !----------------------------------------------------------------------
83 IMPLICIT NONE
84 !<DESCRIPTION>
85 !
86 ! Radiation_driver is the WRF mediation layer routine that provides the interface to
87 ! to radiation physics packages in the WRF model layer. The radiation
88 ! physics packages to call are chosen by setting the namelist variable
89 ! (Rconfig entry in Registry) to the integer value assigned to the
90 ! particular package (package entry in Registry). For example, if the
91 ! namelist variable ra_lw_physics is set to 1, this corresponds to the
92 ! Registry Package entry for swradscheme. Note that the Package
93 ! names in the Registry are defined constants (frame/module_state_description.F)
94 ! in the CASE statements in this routine.
95 !
96 ! Among the arguments is moist, a four-dimensional scalar array storing
97 ! a variable number of moisture tracers, depending on the physics
98 ! configuration for the WRF run, as determined in the namelist. The
99 ! highest numbered index of active moisture tracers the integer argument
100 ! n_moist (note: the number of tracers at run time is the quantity
101 ! <tt>n_moist - PARAM_FIRST_SCALAR + 1</tt> , not n_moist. Individual tracers
102 ! may be indexed from moist by the Registry name of the tracer prepended
103 ! with P_; for example P_QC is the index of cloud water. An index
104 ! represents a valid, active field only if the index is greater than
105 ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual
106 ! indices for each tracer is defined in module_state_description and
107 ! set in <a href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a> defined in frame/module_configure.F.
108 !
109 ! Physics drivers in WRF 2.0 and higher, originally model-layer
110 ! routines, have been promoted to mediation layer routines and they
111 ! contain OpenMP threaded loops over tiles. Thus, physics drivers
112 ! are called from single-threaded regions in the solver. The physics
113 ! routines that are called from the physics drivers are model-layer
114 ! routines and fully tile-callable and thread-safe.
115 !</DESCRIPTION>
116 !
117 !======================================================================
118 ! Grid structure in physics part of WRF
119 !----------------------------------------------------------------------
120 ! The horizontal velocities used in the physics are unstaggered
121 ! relative to temperature/moisture variables. All predicted
122 ! variables are carried at half levels except w, which is at full
123 ! levels. Some arrays with names (*8w) are at w (full) levels.
124 !
125 !----------------------------------------------------------------------
126 ! In WRF, kms (smallest number) is the bottom level and kme (largest
127 ! number) is the top level. In your scheme, if 1 is at the top level,
128 ! then you have to reverse the order in the k direction.
129 !
130 ! kme - half level (no data at this level)
131 ! kme ----- full level
132 ! kme-1 - half level
133 ! kme-1 ----- full level
134 ! .
135 ! .
136 ! .
137 ! kms+2 - half level
138 ! kms+2 ----- full level
139 ! kms+1 - half level
140 ! kms+1 ----- full level
141 ! kms - half level
142 ! kms ----- full level
143 !
144 !======================================================================
145 ! Grid structure in physics part of WRF
146 !
147 !-------------------------------------
148 ! The horizontal velocities used in the physics are unstaggered
149 ! relative to temperature/moisture variables. All predicted
150 ! variables are carried at half levels except w, which is at full
151 ! levels. Some arrays with names (*8w) are at w (full) levels.
152 !
153 !==================================================================
154 ! Definitions
155 !-----------
156 ! Theta potential temperature (K)
157 ! Qv water vapor mixing ratio (kg/kg)
158 ! Qc cloud water mixing ratio (kg/kg)
159 ! Qr rain water mixing ratio (kg/kg)
160 ! Qi cloud ice mixing ratio (kg/kg)
161 ! Qs snow mixing ratio (kg/kg)
162 !-----------------------------------------------------------------
163 !-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3)
164 !-- PM2_5_WATER PM2.5 water mass (ug m^-3)
165 !-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3)
166 !-- RTHRATEN Theta tendency
167 ! due to radiation (K/s)
168 !-- RTHRATENLW Theta tendency
169 ! due to long wave radiation (K/s)
170 !-- RTHRATENSW Theta temperature tendency
171 ! due to short wave radiation (K/s)
172 !-- dt time step (s)
173 !-- itimestep number of time steps
174 !-- GLW downward long wave flux at ground surface (W/m^2)
175 !-- GSW net short wave flux at ground surface (W/m^2)
176 !-- SWDOWN downward short wave flux at ground surface (W/m^2)
177 !-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ)
178 !-- RLWTOA upward long wave at top of atmosphere (w/m2)
179 !-- RSWTOA upward short wave at top of atmosphere (w/m2)
180 !-- XLAT latitude, south is negative (degree)
181 !-- XLONG longitude, west is negative (degree)
182 !-- ALBEDO albedo (between 0 and 1)
183 !-- CLDFRA cloud fraction (between 0 and 1)
184 !-- EMISS surface emissivity (between 0 and 1)
185 !-- rho_phy density (kg/m^3)
186 !-- rr dry air density (kg/m^3)
187 !-- moist moisture array (4D - last index is species) (kg/kg)
188 !-- n_moist number of moisture species
189 !-- qndrop Cloud droplet number (#/kg)
190 !-- p8w pressure at full levels (Pa)
191 !-- p_phy pressure (Pa)
192 !-- Pb base-state pressure (Pa)
193 !-- pi_phy exner function (dimensionless)
194 !-- dz8w dz between full levels (m)
195 !-- t_phy temperature (K)
196 !-- t8w temperature at full levels (K)
197 !-- GMT Greenwich Mean Time Hour of model start (hour)
198 !-- JULDAY the initial day (Julian day)
199 !-- RADT time for calling radiation (min)
200 !-- ra_call_offset -1 (old) means usually just before output, 0 after
201 !-- DEGRAD conversion factor for
202 ! degrees to radians (pi/180.) (rad/deg)
203 !-- DPD degrees per day for earth's
204 ! orbital position (deg/day)
205 !-- R_d gas constant for dry air (J/kg/K)
206 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
207 !-- G acceleration due to gravity (m/s^2)
208 !-- rvovrd R_v divided by R_d (dimensionless)
209 !-- XTIME time since simulation start (min)
210 !-- DECLIN solar declination angle (rad)
211 !-- SOLCON solar constant (W/m^2)
212 !-- ids start index for i in domain
213 !-- ide end index for i in domain
214 !-- jds start index for j in domain
215 !-- jde end index for j in domain
216 !-- kds start index for k in domain
217 !-- kde end index for k in domain
218 !-- ims start index for i in memory
219 !-- ime end index for i in memory
220 !-- jms start index for j in memory
221 !-- jme end index for j in memory
222 !-- kms start index for k in memory
223 !-- kme end index for k in memory
224 !-- i_start start indices for i in tile
225 !-- i_end end indices for i in tile
226 !-- j_start start indices for j in tile
227 !-- j_end end indices for j in tile
228 !-- kts start index for k in tile
229 !-- kte end index for k in tile
230 !-- num_tiles number of tiles
231 !
232 !==================================================================
233 !
234 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
235 ims,ime, jms,jme, kms,kme, &
236 kts,kte, &
237 num_tiles
238
239 INTEGER, INTENT(IN) :: lw_physics, sw_physics
240
241 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
242 & i_start,i_end,j_start,j_end
243
244 INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset
245 INTEGER, INTENT(IN ) :: levsiz, n_ozmixm
246 INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2
247 REAL, INTENT(IN ) :: cam_abs_freq_s
248
249 LOGICAL, INTENT(IN ) :: warm_rain
250
251 REAL, INTENT(IN ) :: RADT
252
253 REAL, DIMENSION( ims:ime, jms:jme ), &
254 INTENT(IN ) :: XLAND, &
255 XICE, &
256 TSK, &
257 VEGFRA, &
258 SNOW
259 REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, &
260 INTENT(IN ) :: OZMIXM
261
262 REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN
263
264 REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2
265 REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, &
266 INTENT(IN ) :: aerosolc_1, aerosolc_2
267 REAL, DIMENSION(paerlev), OPTIONAL, &
268 INTENT(IN ) :: m_hybi0
269
270 REAL, DIMENSION( ims:ime, jms:jme ), &
271 INTENT(INOUT) :: HTOP, &
272 HBOT, &
273 HTOPR, &
274 HBOTR, &
275 CUPPT
276
277 INTEGER, INTENT(IN ) :: julyr
278 !
279 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
280 INTENT(IN ) :: dz8w, &
281 z, &
282 p8w, &
283 p, &
284 pi, &
285 t, &
286 t8w, &
287 rho
288 !
289 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
290 INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
291 gaer300,gaer400,gaer600,gaer999, & ! jcb
292 waer300,waer400,waer600,waer999, & ! jcb
293 qc_adjust, qi_adjust
294
295 LOGICAL, OPTIONAL :: cu_rad_feedback
296
297 INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
298
299 !
300 ! variables for aerosols (only if running with chemistry)
301 !
302 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
303 INTENT(IN ) :: pm2_5_dry, &
304 pm2_5_water, &
305 pm2_5_dry_ec
306 !
307 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
308 INTENT(INOUT) :: RTHRATEN, &
309 RTHRATENLW, &
310 RTHRATENSW
311
312 ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
313 ! INTENT(INOUT) :: SWUP, &
314 ! SWDN, &
315 ! SWUPCLEAR, &
316 ! SWDNCLEAR, &
317 ! LWUP, &
318 ! LWDN, &
319 ! LWUPCLEAR, &
320 ! LWDNCLEAR
321
322 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
323 ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
324 ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
325 ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
326 ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
327 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
328 SWUPT, SWUPTC, SWDNT, SWDNTC, &
329 SWUPB, SWUPBC, SWDNB, SWDNBC, &
330 LWUPT, LWUPTC, LWDNT, LWDNTC, &
331 LWUPB, LWUPBC, LWDNB, LWDNBC
332
333 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
334 INTENT(INOUT) :: SWCF, &
335 LWCF, &
336 OLR
337
338
339 !
340 REAL, DIMENSION( ims:ime, jms:jme ), &
341 INTENT(IN ) :: XLAT, &
342 XLONG, &
343 ALBEDO, &
344 EMISS
345 !
346 REAL, DIMENSION( ims:ime, jms:jme ), &
347 INTENT(INOUT) :: GSW, &
348 GLW
349
350 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN
351 !
352 REAL, INTENT(IN ) :: GMT,dt, &
353 julian, xtime
354 !
355 INTEGER, INTENT(IN ) :: JULDAY, itimestep
356
357 INTEGER,INTENT(IN) :: NPHS
358 REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: &
359 CFRACH, & !Added
360 CFRACL, & !Added
361 CFRACM, & !Added
362 CZMEAN !Added
363 REAL, DIMENSION( ims:ime, jms:jme ), &
364 INTENT(INOUT) :: &
365 RLWTOA, & !Added
366 RSWTOA, & !Added
367 ACFRST, & !Added
368 ACFRCV !Added
369
370 INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: &
371 NCFRST, & !Added
372 NCFRCV !Added
373 ! Optional (only used by CAM lw scheme)
374
375 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,&
376 INTENT(INOUT) :: abstot
377 REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,&
378 INTENT(INOUT) :: absnxt
379 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
380 INTENT(INOUT) :: emstot
381
382 !
383 ! Optional
384 !
385 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
386 OPTIONAL, &
387 INTENT(INOUT) :: CLDFRA
388
389 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
390 OPTIONAL, &
391 INTENT(IN ) :: &
392 F_ICE_PHY, &
393 F_RAIN_PHY
394
395 REAL, DIMENSION( ims:ime, jms:jme ), &
396 OPTIONAL, &
397 INTENT(OUT) :: SWDOWNC
398 !
399 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
400 OPTIONAL, &
401 INTENT(INOUT ) :: &
402 pb &
403 ,qv,qc,qr,qi,qs,qg,qndrop
404
405 LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop
406 !
407 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
408 OPTIONAL, &
409 INTENT(INOUT) :: taucldi,taucldc
410
411 ! LOCAL VAR
412
413 REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON
414 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS
415 REAL, DIMENSION( ims:ime, jms:jme ) :: coszr
416
417 REAL :: DECLIN,SOLCON
418 INTEGER :: i,j,k,its,ite,jts,jte,ij
419 INTEGER :: STEPABS
420 LOGICAL :: gfdl_lw,gfdl_sw
421 LOGICAL :: doabsems
422 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
423
424 REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, &
425 DJUL,RJUL,ECCFAC
426 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_temp,qc_temp
427 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save
428
429 !------------------------------------------------------------------
430 ! urban related variables are added to declaration
431 !-------------------------------------------------
432 REAL, OPTIONAL, INTENT(OUT) :: DECLIN_URB !urban
433 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZ_URB2D !urban
434 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: omg_urb2d !urban
435 !------------------------------------------------------------------
436
437 if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return
438
439 ! ra_call_offset = -1 gives old method where radiation may be called just before output
440 ! ra_call_offset = 0 gives new method where radiation may be called just after output
441 ! and is also consistent with removal of offset in new XTIME
442 ! also need to account for stepra=1 which always has zero modulo output
443 Radiation_step: IF (itimestep .eq. 1 .or. mod(itimestep,STEPRA) .eq. 1 + ra_call_offset &
444 .or. STEPRA .eq. 1 ) THEN
445
446 ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
447 STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA
448 IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset &
449 .or. STEPABS .eq. 1 ) THEN
450 doabsems = .true.
451 ELSE
452 doabsems = .false.
453 ENDIF
454
455 gfdl_lw = .false.
456 gfdl_sw = .false.
457
458 !---------------
459 !$OMP PARALLEL DO &
460 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
461
462 DO ij = 1 , num_tiles
463 its = i_start(ij)
464 ite = i_end(ij)
465 jts = j_start(ij)
466 jte = j_end(ij)
467
468 ! initialize data
469
470 DO j=jts,jte
471 DO i=its,ite
472 GSW(I,J)=0.
473 GLW(I,J)=0.
474 SWDOWN(I,J)=0.
475 GLAT(I,J)=XLAT(I,J)*DEGRAD
476 GLON(I,J)=XLONG(I,J)*DEGRAD
477 ENDDO
478 ENDDO
479
480 DO j=jts,jte
481 DO k=kts,kte+1
482 DO i=its,ite
483 RTHRATEN(I,K,J)=0.
484 ! SWUP(I,K,J) = 0.0
485 ! SWDN(I,K,J) = 0.0
486 ! SWUPCLEAR(I,K,J) = 0.0
487 ! SWDNCLEAR(I,K,J) = 0.0
488 ! LWUP(I,K,J) = 0.0
489 ! LWDN(I,K,J) = 0.0
490 ! LWUPCLEAR(I,K,J) = 0.0
491 ! LWDNCLEAR(I,K,J) = 0.0
492 CEMISS(I,K,J)=0.0
493 ENDDO
494 ENDDO
495 ENDDO
496
497 ! temporarily modify hydrometeors (currently only done for GD scheme and WRF-Chem)
498 !
499 IF ( PRESENT( cu_rad_feedback ) ) THEN
500 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
501 DO j=jts,jte
502 DO k=kts,kte
503 DO i=its,ite
504 qc_save(i,k,j) = qc(i,k,j)
505 qc(i,k,j) = qc(i,k,j) + qc_adjust(i,k,j)
506 ENDDO
507 ENDDO
508 ENDDO
509 ENDIF
510 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
511 DO j=jts,jte
512 DO k=kts,kte
513 DO i=its,ite
514 qi_save(i,k,j) = qi(i,k,j)
515 qi(i,k,j) = qi(i,k,j) + qi_adjust(i,k,j)
516 ENDDO
517 ENDDO
518 ENDDO
519 ENDIF
520 ENDIF
521
522
523 ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
524 if(PRESENT(qc) .and. PRESENT(F_QC)) then
525 DO j=jts,jte
526 DO k=kts,kte
527 DO i=its,ite
528 qc_temp(I,K,J)=qc(I,K,J)
529 ENDDO
530 ENDDO
531 ENDDO
532 else
533 DO j=jts,jte
534 DO k=kts,kte
535 DO i=its,ite
536 qc_temp(I,K,J)=0.
537 ENDDO
538 ENDDO
539 ENDDO
540 endif
541 if(PRESENT(qr) .and. PRESENT(F_QR)) then
542 DO j=jts,jte
543 DO k=kts,kte
544 DO i=its,ite
545 qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J)
546 ENDDO
547 ENDDO
548 ENDDO
549 endif
550
551 !---------------
552 ! Calculate constant for short wave radiation
553
554 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
555 DEGRAD,DPD )
556
557 if(present(DECLIN_URB))DECLIN_URB=DECLIN ! urban
558
559 lwrad_cldfra_select: SELECT CASE(lw_physics)
560
561 CASE (GFDLLWSCHEME)
562
563 !-- Do nothing, since cloud fractions (with partial cloudiness effects)
564 !-- are defined in GFDL LW/SW schemes and do not need to be initialized.
565
566 CASE (CAMLWSCHEME)
567
568 IF ( PRESENT ( CLDFRA ) .AND. &
569 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
570 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
571
572 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
573 F_QV,F_QC,F_QI,F_QS,t,p, &
574 F_ICE_PHY,F_RAIN_PHY, &
575 ids,ide, jds,jde, kds,kde, &
576 ims,ime, jms,jme, kms,kme, &
577 its,ite, jts,jte, kts,kte )
578 ENDIF
579
580 CASE DEFAULT
581
582 IF ( PRESENT ( CLDFRA ) .AND. &
583 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
584 CALL cal_cldfra(CLDFRA,qc,qi,F_QC,F_QI, &
585 ids,ide, jds,jde, kds,kde, &
586 ims,ime, jms,jme, kms,kme, &
587 its,ite, jts,jte, kts,kte )
588 ENDIF
589
590 END SELECT lwrad_cldfra_select
591
592 !pjj/cray Cray X1 cannot print from threaded region
593 #ifndef crayx1
594 WRITE(wrf_err_message,*)'SOLCON=',SOLCON,DECLIN,XTIME
595 CALL wrf_debug(50,wrf_err_message)
596 #endif
597
598 lwrad_select: SELECT CASE(lw_physics)
599
600 CASE (RRTMSCHEME)
601 CALL wrf_debug (100, 'CALL rrtm')
602
603 CALL RRTMLWRAD( &
604 RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS &
605 ,QV3D=QV &
606 ,QC3D=QC &
607 ,QR3D=QR &
608 ,QI3D=QI &
609 ,QS3D=QS &
610 ,QG3D=QG &
611 ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t &
612 ,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G &
613 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
614 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
615 ,ICLOUD=icloud,WARM_RAIN=warm_rain &
616 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
617 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
618 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
619 )
620
621 CASE (GFDLLWSCHEME)
622
623 CALL wrf_debug (100, 'CALL gfdllw')
624
625 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
626 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
627 PRESENT(qv) .AND. PRESENT(qc) ) THEN
628 IF ( F_QV .AND. F_QC .AND. F_QS) THEN
629 gfdl_lw = .true.
630 CALL ETARA( &
631 DT=dt,XLAND=xland &
632 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
633 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
634 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
635 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
636 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
637 ,HBOTR=hbotr, HTOPR=htopr &
638 ,ALBEDO=albedo,CUPPT=cuppt &
639 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
640 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
641 ,XTIME=xtime,JULIAN=julian &
642 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
643 ,JULYR=julyr,JULDAY=julday &
644 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
645 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
646 ,ACFRST=acfrst,NCFRST=ncfrst &
647 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
648 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
649 ,THRATEN=rthraten,THRATENLW=rthratenlw &
650 ,THRATENSW=rthratensw &
651 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
652 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
653 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
654 )
655 ELSE
656 CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
657 ENDIF
658 ELSE
659 CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
660 ENDIF
661 CASE (CAMLWSCHEME)
662 CALL wrf_debug(100, 'CALL camrad lw')
663 IF(cam_abs_dim1 .ne. 4 .or. cam_abs_dim2 .ne. kde .or. &
664 paerlev .ne. 29 .or. levsiz .ne. 59 )THEN
665 WRITE( wrf_err_message , * ) &
666 'set paerlev=29, levsiz=59, cam_abs_dim1=4, and cam_abs_dim2=number of levels (e_vert) in physics namelist for CAM radiation'
667 CALL wrf_error_fatal ( wrf_err_message )
668 ENDIF
669 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
670 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
671 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
672 .AND. PRESENT(AEROSOLC_2) ) THEN
673 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
674 SWUPT=SWUPT,SWUPTC=SWUPTC, &
675 SWDNT=SWDNT,SWDNTC=SWDNTC, &
676 LWUPT=LWUPT,LWUPTC=LWUPTC, &
677 LWDNT=LWDNT,LWDNTC=LWDNTC, &
678 SWUPB=SWUPB,SWUPBC=SWUPBC, &
679 SWDNB=SWDNB,SWDNBC=SWDNBC, &
680 LWUPB=LWUPB,LWUPBC=LWUPBC, &
681 LWDNB=LWDNB,LWDNBC=LWDNBC, &
682 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
683 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
684 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
685 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
686 ,QV3D=qv &
687 ,QC3D=qc &
688 ,QR3D=qr &
689 ,QI3D=qi &
690 ,QS3D=qs &
691 ,QG3D=qg &
692 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
693 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
694 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
695 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
696 dz8w=dz8w, &
697 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
698 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
699 num_months=n_ozmixm, &
700 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
701 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
702 paerlev=paerlev, naer_c=n_aerosolc, &
703 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
704 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
705 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
706 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
707 ,doabsems=doabsems &
708 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
709 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
710 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
711 )
712 ELSE
713 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
714 ENDIF
715
716
717 CASE (HELDSUAREZ)
718 CALL wrf_debug (100, 'CALL heldsuarez')
719
720 CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, &
721 t8w, rho, R_d,G,CP, dt, xlat, degrad, &
722 ids,ide, jds,jde, kds,kde, &
723 ims,ime, jms,jme, kms,kme, &
724 its,ite, jts,jte, kts,kte )
725
726 CASE DEFAULT
727
728 WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
729 CALL wrf_error_fatal ( wrf_err_message )
730
731 END SELECT lwrad_select
732
733 IF (lw_physics .gt. 0 .and. .not.gfdl_lw) THEN
734 DO j=jts,jte
735 DO k=kts,kte
736 DO i=its,ite
737 RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
738 ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY)
739 IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
740 ENDDO
741 ENDDO
742 ENDDO
743 ENDIF
744 !
745
746 swrad_select: SELECT CASE(sw_physics)
747
748 CASE (SWRADSCHEME)
749 CALL wrf_debug(100, 'CALL swrad')
750 CALL SWRAD( &
751 DT=dt,RTHRATEN=rthraten,GSW=gsw &
752 ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
753 #ifdef WRF_CHEM
754 ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
755 ,PM2_5_DRY_EC=pm2_5_dry_ec &
756 #endif
757 ,RHO_PHY=rho,T3D=t &
758 ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
759 ,R=r_d,CP=cp,G=g,JULDAY=julday &
760 ,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
761 ! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
762 ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
763 ,warm_rain=warm_rain &
764 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
765 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
766 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
767 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
768 ,QV3D=qv &
769 ,QC3D=qc &
770 ,QR3D=qr &
771 ,QI3D=qi &
772 ,QS3D=qs &
773 ,QG3D=qg &
774 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
775 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
776 )
777
778 CASE (GSFCSWSCHEME)
779 CALL wrf_debug(100, 'CALL gsfcswrad')
780 CALL GSFCSWRAD( &
781 RTHRATEN=rthraten,GSW=gsw,XLAT=xlat,XLONG=xlong &
782 ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
783 ,DZ8W=dz8w,RHO_PHY=rho &
784 ,CLDFRA3D=cldfra,RSWTOA=rswtoa &
785 ,GMT=gmt,CP=cp,G=g &
786 ! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
787 ,JULDAY=julday,XTIME=xtime &
788 ,DECLIN=declin,SOLCON=solcon &
789 ,RADFRQ=radt,DEGRAD=degrad &
790 ,TAUCLDI=taucldi,TAUCLDC=taucldc &
791 ,WARM_RAIN=warm_rain &
792 #ifdef WRF_CHEM
793 ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
794 ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
795 ,GAER300=gaer300,GAER400=gaer400 & ! jcb
796 ,GAER600=gaer600,GAER999=gaer999 & ! jcb
797 ,WAER300=waer300,WAER400=waer400 & ! jcb
798 ,WAER600=waer600,WAER999=waer999 & ! jcb
799 ,aer_ra_feedback=aer_ra_feedback &
800 #endif
801 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
802 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
803 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
804 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
805 ,QV3D=qv &
806 ,QC3D=qc &
807 ,QR3D=qr &
808 ,QI3D=qi &
809 ,QS3D=qs &
810 ,QG3D=qg &
811 ,QNDROP3D=qndrop &
812 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
813 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
814 ,F_QNDROP=f_qndrop &
815 )
816 CASE (CAMSWSCHEME)
817 ! Temporarily lw switch already calculates sw CAM tendency, so inactive here
818
819 DO j=jts,jte
820 DO k=kts,kte
821 DO i=its,ite
822 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
823 ENDDO
824 ENDDO
825 ENDDO
826
827 CASE (GFDLSWSCHEME)
828
829 CALL wrf_debug (100, 'CALL gfdlsw')
830
831 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
832 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
833 PRESENT(qv) .AND. PRESENT(qc) ) THEN
834 IF ( F_QV .AND. F_QC .AND. F_QS ) THEN
835 gfdl_sw = .true.
836 CALL ETARA( &
837 DT=dt,XLAND=xland &
838 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
839 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
840 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
841 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
842 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
843 ,HBOTR=hbotr, HTOPR=htopr &
844 ,ALBEDO=albedo,CUPPT=cuppt &
845 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
846 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
847 ,XTIME=xtime,JULIAN=julian &
848 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
849 ,JULYR=julyr,JULDAY=julday &
850 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
851 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
852 ,ACFRST=acfrst,NCFRST=ncfrst &
853 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
854 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
855 ,THRATEN=rthraten,THRATENLW=rthratenlw &
856 ,THRATENSW=rthratensw &
857 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
858 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
859 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
860 )
861 ELSE
862 CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
863 ENDIF
864 ELSE
865 CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
866 ENDIF
867
868 CASE (0)
869
870 ! Here in case we don't want to call a sw radiation scheme
871 ! For example, the Held-Suarez idealized test case
872 IF (lw_physics /= HELDSUAREZ) THEN
873 WRITE( wrf_err_message , * ) 'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')'
874 CALL wrf_error_fatal ( wrf_err_message )
875 END IF
876
877 CASE DEFAULT
878
879 WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
880 CALL wrf_error_fatal ( wrf_err_message )
881
882 END SELECT swrad_select
883
884 IF (sw_physics .gt. 0 .and. .not.gfdl_sw) THEN
885 DO j=jts,jte
886 DO k=kts,kte
887 DO i=its,ite
888 RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
889 ENDDO
890 ENDDO
891 ENDDO
892
893 DO j=jts,jte
894 DO i=its,ite
895 SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
896 ENDDO
897 ENDDO
898
899 ENDIF
900
901 IF ( PRESENT( cu_rad_feedback ) ) THEN
902 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
903 DO j=jts,jte
904 DO k=kts,kte
905 DO i=its,ite
906 qc(i,k,j) = qc_save(i,k,j)
907 ENDDO
908 ENDDO
909 ENDDO
910 ENDIF
911 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
912 DO j=jts,jte
913 DO k=kts,kte
914 DO i=its,ite
915 qi(i,k,j) = qi_save(i,k,j)
916 ENDDO
917 ENDDO
918 ENDDO
919 ENDIF
920 ENDIF
921
922 ENDDO
923 !$OMP END PARALLEL DO
924
925 ENDIF Radiation_step
926
927 accumulate_lw_select: SELECT CASE(lw_physics)
928
929 CASE (CAMLWSCHEME)
930 IF(PRESENT(LWUPTC))THEN
931 !$OMP PARALLEL DO &
932 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
933
934 DO ij = 1 , num_tiles
935 its = i_start(ij)
936 ite = i_end(ij)
937 jts = j_start(ij)
938 jte = j_end(ij)
939
940 DO j=jts,jte
941 DO i=its,ite
942 ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DT
943 ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DT
944 ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DT
945 ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DT
946 ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DT
947 ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DT
948 ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DT
949 ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DT
950 ENDDO
951 ENDDO
952 ENDDO
953 !$OMP END PARALLEL DO
954 ENDIF
955 CASE DEFAULT
956 END SELECT accumulate_lw_select
957
958 accumulate_sw_select: SELECT CASE(sw_physics)
959
960 CASE (CAMSWSCHEME)
961 IF(PRESENT(SWUPTC))THEN
962 !$OMP PARALLEL DO &
963 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
964
965 DO ij = 1 , num_tiles
966 its = i_start(ij)
967 ite = i_end(ij)
968 jts = j_start(ij)
969 jte = j_end(ij)
970
971 DO j=jts,jte
972 DO i=its,ite
973 ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DT
974 ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DT
975 ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DT
976 ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DT
977 ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DT
978 ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DT
979 ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DT
980 ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DT
981 ENDDO
982 ENDDO
983 ENDDO
984 !$OMP END PARALLEL DO
985 ENDIF
986
987 CASE DEFAULT
988 END SELECT accumulate_sw_select
989
990 END SUBROUTINE radiation_driver
991
992 !---------------------------------------------------------------------
993 !BOP
994 ! !IROUTINE: radconst - compute radiation terms
995 ! !INTERFAC:
996 SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
997 DEGRAD,DPD )
998 !---------------------------------------------------------------------
999 USE module_wrf_error
1000 IMPLICIT NONE
1001 !---------------------------------------------------------------------
1002
1003 ! !ARGUMENTS:
1004 REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
1005 REAL, INTENT(OUT ) :: DECLIN,SOLCON
1006 REAL :: OBECL,SINOB,SXLONG,ARG, &
1007 DECDEG,DJUL,RJUL,ECCFAC
1008 !
1009 ! !DESCRIPTION:
1010 ! Compute terms used in radiation physics
1011 !EOP
1012
1013 ! for short wave radiation
1014
1015 DECLIN=0.
1016 SOLCON=0.
1017
1018 !-----OBECL : OBLIQUITY = 23.5 DEGREE.
1019
1020 OBECL=23.5*DEGRAD
1021 SINOB=SIN(OBECL)
1022
1023 !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
1024
1025 IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
1026 IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
1027 SXLONG=SXLONG*DEGRAD
1028 ARG=SINOB*SIN(SXLONG)
1029 DECLIN=ASIN(ARG)
1030 DECDEG=DECLIN/DEGRAD
1031 !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
1032 DJUL=JULIAN*360./365.
1033 RJUL=DJUL*DEGRAD
1034 ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* &
1035 COS(2*RJUL)+0.000077*SIN(2*RJUL)
1036 SOLCON=1370.*ECCFAC
1037
1038 !pjj/cray Cray X1 cannot print from threaded region
1039 #ifndef crayx1
1040 write(wrf_err_message,10)DECDEG,SOLCON
1041 10 FORMAT(1X,'*** SOLAR DECLINATION ANGLE = ',F6.2,' DEGREES.', &
1042 ' SOLAR CONSTANT = ',F8.2,' W/M**2 ***')
1043 CALL wrf_debug (50, wrf_err_message)
1044 #endif
1045
1046 END SUBROUTINE radconst
1047
1048 !---------------------------------------------------------------------
1049 !BOP
1050 ! !IROUTINE: cal_cldfra - Compute cloud fraction
1051 ! !INTERFACE:
1052 SUBROUTINE cal_cldfra(CLDFRA,QC,QI,F_QC,F_QI, &
1053 ids,ide, jds,jde, kds,kde, &
1054 ims,ime, jms,jme, kms,kme, &
1055 its,ite, jts,jte, kts,kte )
1056 !---------------------------------------------------------------------
1057 IMPLICIT NONE
1058 !---------------------------------------------------------------------
1059 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1060 ims,ime, jms,jme, kms,kme, &
1061 its,ite, jts,jte, kts,kte
1062
1063 !
1064 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1065 CLDFRA
1066
1067 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1068 QI, &
1069 QC
1070
1071 LOGICAL,INTENT(IN) :: F_QC,F_QI
1072
1073 REAL thresh
1074 INTEGER:: i,j,k
1075 ! !DESCRIPTION:
1076 ! Compute cloud fraction from input ice and cloud water fields
1077 ! if provided.
1078 !
1079 ! Whether QI or QC is active or not is determined from the indices of
1080 ! the fields into the 4D scalar arrays in WRF. These indices are
1081 ! P_QI and P_QC, respectively, and they are passed in to the routine
1082 ! to enable testing to see if QI and QC represent active fields in
1083 ! the moisture 4D scalar array carried by WRF.
1084 !
1085 ! If a field is active its index will have a value greater than or
1086 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1087 ! this routine.
1088 !EOP
1089 !---------------------------------------------------------------------
1090 thresh=1.0e-6
1091
1092 IF ( f_qi .AND. f_qc ) THEN
1093 DO j = jts,jte
1094 DO k = kts,kte
1095 DO i = its,ite
1096 IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
1097 CLDFRA(i,k,j)=1.
1098 ELSE
1099 CLDFRA(i,k,j)=0.
1100 ENDIF
1101 ENDDO
1102 ENDDO
1103 ENDDO
1104 ELSE IF ( f_qc ) THEN
1105 DO j = jts,jte
1106 DO k = kts,kte
1107 DO i = its,ite
1108 IF ( QC(i,k,j) .gt. thresh) THEN
1109 CLDFRA(i,k,j)=1.
1110 ELSE
1111 CLDFRA(i,k,j)=0.
1112 ENDIF
1113 ENDDO
1114 ENDDO
1115 ENDDO
1116 ELSE
1117 DO j = jts,jte
1118 DO k = kts,kte
1119 DO i = its,ite
1120 CLDFRA(i,k,j)=0.
1121 ENDDO
1122 ENDDO
1123 ENDDO
1124 ENDIF
1125
1126 END SUBROUTINE cal_cldfra
1127
1128 !BOP
1129 ! !IROUTINE: cal_cldfra2 - Compute cloud fraction
1130 ! !INTERFACE:
1131 ! cal_cldfra_xr - Compute cloud fraction.
1132 ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
1133 !!
1134 !!--- Cloud fraction parameterization follows Randall, 1994
1135 !! (see Hong et al., 1998)
1136 !! (modified by Ferrier, Feb '02)
1137 !
1138 SUBROUTINE cal_cldfra2(CLDFRA, QV, QC, QI, QS, &
1139 F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
1140 F_ICE_PHY,F_RAIN_PHY, &
1141 ids,ide, jds,jde, kds,kde, &
1142 ims,ime, jms,jme, kms,kme, &
1143 its,ite, jts,jte, kts,kte )
1144 !---------------------------------------------------------------------
1145 IMPLICIT NONE
1146 !---------------------------------------------------------------------
1147 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1148 ims,ime, jms,jme, kms,kme, &
1149 its,ite, jts,jte, kts,kte
1150
1151 !
1152 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1153 CLDFRA
1154
1155 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1156 QV, &
1157 QI, &
1158 QC, &
1159 QS, &
1160 t_phy, &
1161 p_phy, &
1162 F_ICE_PHY, &
1163 F_RAIN_PHY
1164
1165 LOGICAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
1166
1167 ! REAL thresh
1168 INTEGER:: i,j,k
1169 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
1170
1171 REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
1172 PEXP=0.25, RHGRID=1.0
1173 REAL , PARAMETER :: SVP1=0.61078
1174 REAL , PARAMETER :: SVP2=17.2693882
1175 REAL , PARAMETER :: SVPI2=21.8745584
1176 REAL , PARAMETER :: SVP3=35.86
1177 REAL , PARAMETER :: SVPI3=7.66
1178 REAL , PARAMETER :: SVPT0=273.15
1179 REAL , PARAMETER :: r_d = 287.
1180 REAL , PARAMETER :: r_v = 461.6
1181 REAL , PARAMETER :: ep_2=r_d/r_v
1182 ! !DESCRIPTION:
1183 ! Compute cloud fraction from input ice and cloud water fields
1184 ! if provided.
1185 !
1186 ! Whether QI or QC is active or not is determined from the indices of
1187 ! the fields into the 4D scalar arrays in WRF. These indices are
1188 ! P_QI and P_QC, respectively, and they are passed in to the routine
1189 ! to enable testing to see if QI and QC represent active fields in
1190 ! the moisture 4D scalar array carried by WRF.
1191 !
1192 ! If a field is active its index will have a value greater than or
1193 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1194 ! this routine.
1195 !EOP
1196
1197
1198 !-----------------------------------------------------------------------
1199 !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
1200 ! (modified by Ferrier, Feb '02)
1201 !
1202 !--- Cloud fraction parameterization follows Randall, 1994
1203 ! (see Hong et al., 1998)
1204 !-----------------------------------------------------------------------
1205 ! Note: ep_2=287./461.6 Rd/Rv
1206 ! Note: R_D=287.
1207
1208 ! Alternative calculation for critical RH for grid saturation
1209 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
1210
1211 ! Calculate saturation mixing ratio weighted according to the fractions of
1212 ! water and ice.
1213 ! Following:
1214 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
1215 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
1216 !
1217 ! over ice over water
1218 ! a = 21.8745584 17.2693882
1219 ! b = 7.66 35.86
1220
1221 !---------------------------------------------------------------------
1222
1223 DO j = jts,jte
1224 DO k = kts,kte
1225 DO i = its,ite
1226 tc = t_phy(i,k,j) - SVPT0
1227 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
1228 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
1229 QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
1230 QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
1231
1232 IF ( F_QI .and. F_QC .and. F_QS) THEN
1233 QCLD=QI(i,k,j)+QC(i,k,j)+QS(I,k,j)
1234 IF (QCLD .LT. QCLDMIN) THEN
1235 weight = 0.
1236 ELSE
1237 weight = (QI(i,k,j)+QS(I,k,j)) / QCLD
1238 ENDIF
1239 ELSE IF ( F_QC ) THEN
1240
1241 ! Mixing ratios of cloud water & total ice (cloud ice + snow).
1242 ! Mixing ratios of rain are not considered in this scheme.
1243 ! F_ICE is fraction of ice
1244 ! F_RAIN is fraction of rain
1245
1246 QIMID=QC(i,k,j)*F_ICE_PHY(i,k,j)
1247 QWMID=(QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
1248
1249
1250 !
1251 !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
1252 ! only cloud water + cloud ice + snow
1253 !
1254 QCLD=QWMID+QIMID
1255 IF (QCLD .LT. QCLDMIN) THEN
1256 weight = 0.
1257 ELSE
1258 weight = F_ICE_PHY(i,k,j)
1259 ENDIF
1260
1261 ELSE
1262 CLDFRA(i,k,j)=0.
1263 ENDIF ! IF ( F_QI .and. F_QC )
1264
1265
1266 QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI
1267 RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
1268 !
1269 !--- Determine cloud fraction (modified from original algorithm)
1270 !
1271 IF (QCLD .LT. QCLDMIN) THEN
1272 !
1273 !--- Assume zero cloud fraction if there is no cloud mixing ratio
1274 !
1275 CLDFRA(i,k,j)=0.
1276 ELSEIF(RHUM.GE.RHGRID)THEN
1277 !
1278 !--- Assume cloud fraction of unity if near saturation and the cloud
1279 ! mixing ratio is at or above the minimum threshold
1280 !
1281 CLDFRA(i,k,j)=1.
1282 ELSE
1283 !
1284 !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
1285 ! modified based on assumed grid-scale saturation at RH=RHgrid.
1286 !
1287 SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
1288 DENOM=(SUBSAT)**GAMMA
1289 ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
1290 ! prevent negative values (new)
1291 RHUM=MAX(1.E-10, RHUM)
1292 CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG))
1293 !! ARG=-1000*QCLD/(RHUM-RHGRID)
1294 !! ARG=MAX(ARG, ARGMIN)
1295 !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
1296 IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
1297 ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
1298 ENDDO !--- End DO i
1299 ENDDO !--- End DO k
1300 ENDDO !--- End DO j
1301
1302 END SUBROUTINE cal_cldfra2
1303
1304 END MODULE module_radiation_driver