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