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