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 Radiation_step: IF (itimestep .eq. 1 .or. mod(itimestep,STEPRA) .eq. 1 + ra_call_offset) THEN
441
442 ! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
443 STEPABS = nint(cam_abs_freq_s/(dt*STEPRA))*STEPRA
444 IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset) THEN
445 doabsems = .true.
446 ELSE
447 doabsems = .false.
448 ENDIF
449
450 gfdl_lw = .false.
451 gfdl_sw = .false.
452
453 !---------------
454 !$OMP PARALLEL DO &
455 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
456
457 DO ij = 1 , num_tiles
458 its = i_start(ij)
459 ite = i_end(ij)
460 jts = j_start(ij)
461 jte = j_end(ij)
462
463 ! initialize data
464
465 DO j=jts,jte
466 DO i=its,ite
467 GSW(I,J)=0.
468 GLW(I,J)=0.
469 SWDOWN(I,J)=0.
470 GLAT(I,J)=XLAT(I,J)*DEGRAD
471 GLON(I,J)=XLONG(I,J)*DEGRAD
472 ENDDO
473 ENDDO
474
475 DO j=jts,jte
476 DO k=kts,kte+1
477 DO i=its,ite
478 RTHRATEN(I,K,J)=0.
479 ! SWUP(I,K,J) = 0.0
480 ! SWDN(I,K,J) = 0.0
481 ! SWUPCLEAR(I,K,J) = 0.0
482 ! SWDNCLEAR(I,K,J) = 0.0
483 ! LWUP(I,K,J) = 0.0
484 ! LWDN(I,K,J) = 0.0
485 ! LWUPCLEAR(I,K,J) = 0.0
486 ! LWDNCLEAR(I,K,J) = 0.0
487 CEMISS(I,K,J)=0.0
488 ENDDO
489 ENDDO
490 ENDDO
491
492 ! temporarily modify hydrometeors (currently only done for GD scheme and WRF-Chem)
493 !
494 IF ( PRESENT( cu_rad_feedback ) ) THEN
495 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
496 DO j=jts,jte
497 DO k=kts,kte
498 DO i=its,ite
499 qc_save(i,k,j) = qc(i,k,j)
500 qc(i,k,j) = qc(i,k,j) + qc_adjust(i,k,j)
501 ENDDO
502 ENDDO
503 ENDDO
504 ENDIF
505 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
506 DO j=jts,jte
507 DO k=kts,kte
508 DO i=its,ite
509 qi_save(i,k,j) = qi(i,k,j)
510 qi(i,k,j) = qi(i,k,j) + qi_adjust(i,k,j)
511 ENDDO
512 ENDDO
513 ENDDO
514 ENDIF
515 ENDIF
516
517
518 ! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
519 if(PRESENT(qc) .and. PRESENT(F_QC)) then
520 DO j=jts,jte
521 DO k=kts,kte
522 DO i=its,ite
523 qc_temp(I,K,J)=qc(I,K,J)
524 ENDDO
525 ENDDO
526 ENDDO
527 else
528 DO j=jts,jte
529 DO k=kts,kte
530 DO i=its,ite
531 qc_temp(I,K,J)=0.
532 ENDDO
533 ENDDO
534 ENDDO
535 endif
536 if(PRESENT(qr) .and. PRESENT(F_QR)) then
537 DO j=jts,jte
538 DO k=kts,kte
539 DO i=its,ite
540 qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J)
541 ENDDO
542 ENDDO
543 ENDDO
544 endif
545
546 !---------------
547 ! Calculate constant for short wave radiation
548
549 CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
550 DEGRAD,DPD )
551
552 if(present(DECLIN_URB))DECLIN_URB=DECLIN ! urban
553
554 lwrad_cldfra_select: SELECT CASE(lw_physics)
555
556 CASE (GFDLLWSCHEME)
557
558 !-- Do nothing, since cloud fractions (with partial cloudiness effects)
559 !-- are defined in GFDL LW/SW schemes and do not need to be initialized.
560
561 CASE (CAMLWSCHEME)
562
563 IF ( PRESENT ( CLDFRA ) .AND. &
564 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
565 ! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
566
567 CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
568 F_QV,F_QC,F_QI,F_QS,t,p, &
569 F_ICE_PHY,F_RAIN_PHY, &
570 ids,ide, jds,jde, kds,kde, &
571 ims,ime, jms,jme, kms,kme, &
572 its,ite, jts,jte, kts,kte )
573 ENDIF
574
575 CASE DEFAULT
576
577 IF ( PRESENT ( CLDFRA ) .AND. &
578 PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
579 CALL cal_cldfra(CLDFRA,qc,qi,F_QC,F_QI, &
580 ids,ide, jds,jde, kds,kde, &
581 ims,ime, jms,jme, kms,kme, &
582 its,ite, jts,jte, kts,kte )
583 ENDIF
584
585 END SELECT lwrad_cldfra_select
586
587 !pjj/cray Cray X1 cannot print from threaded region
588 #ifndef crayx1
589 WRITE(wrf_err_message,*)'SOLCON=',SOLCON,DECLIN,XTIME
590 CALL wrf_debug(50,wrf_err_message)
591 #endif
592
593 lwrad_select: SELECT CASE(lw_physics)
594
595 CASE (RRTMSCHEME)
596 CALL wrf_debug (100, 'CALL rrtm')
597
598 CALL RRTMLWRAD( &
599 RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS &
600 ,QV3D=QV &
601 ,QC3D=QC &
602 ,QR3D=QR &
603 ,QI3D=QI &
604 ,QS3D=QS &
605 ,QG3D=QG &
606 ,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,T3D=t &
607 ,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G &
608 ,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
609 ,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
610 ,ICLOUD=icloud,WARM_RAIN=warm_rain &
611 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
612 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
613 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
614 )
615
616 CASE (GFDLLWSCHEME)
617
618 CALL wrf_debug (100, 'CALL gfdllw')
619
620 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
621 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
622 PRESENT(qv) .AND. PRESENT(qc) ) THEN
623 IF ( F_QV .AND. F_QC .AND. F_QS) THEN
624 gfdl_lw = .true.
625 CALL ETARA( &
626 DT=dt,XLAND=xland &
627 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
628 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
629 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
630 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
631 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
632 ,HBOTR=hbotr, HTOPR=htopr &
633 ,ALBEDO=albedo,CUPPT=cuppt &
634 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
635 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
636 ,XTIME=xtime,JULIAN=julian &
637 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
638 ,JULYR=julyr,JULDAY=julday &
639 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
640 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
641 ,ACFRST=acfrst,NCFRST=ncfrst &
642 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
643 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
644 ,THRATEN=rthraten,THRATENLW=rthratenlw &
645 ,THRATENSW=rthratensw &
646 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
647 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
648 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
649 )
650 ELSE
651 CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
652 ENDIF
653 ELSE
654 CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
655 ENDIF
656 CASE (CAMLWSCHEME)
657 CALL wrf_debug(100, 'CALL camrad lw')
658 IF(cam_abs_dim1 .ne. 4 .or. cam_abs_dim2 .ne. kde .or. &
659 paerlev .ne. 29 .or. levsiz .ne. 59 )THEN
660 WRITE( wrf_err_message , * ) &
661 'set paerlev=29, levsiz=59, cam_abs_dim1=4, and cam_abs_dim2=number of levels (e_vert) in physics namelist for CAM radiation'
662 CALL wrf_error_fatal ( wrf_err_message )
663 ENDIF
664 IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
665 PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
666 PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
667 .AND. PRESENT(AEROSOLC_2) ) THEN
668 CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
669 SWUPT=SWUPT,SWUPTC=SWUPTC, &
670 SWDNT=SWDNT,SWDNTC=SWDNTC, &
671 LWUPT=LWUPT,LWUPTC=LWUPTC, &
672 LWDNT=LWDNT,LWDNTC=LWDNTC, &
673 SWUPB=SWUPB,SWUPBC=SWUPBC, &
674 SWDNB=SWDNB,SWDNBC=SWDNBC, &
675 LWUPB=LWUPB,LWUPBC=LWUPBC, &
676 LWDNB=LWDNB,LWDNBC=LWDNBC, &
677 SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
678 TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
679 GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
680 ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
681 ,QV3D=qv &
682 ,QC3D=qc &
683 ,QR3D=qr &
684 ,QI3D=qi &
685 ,QS3D=qs &
686 ,QG3D=qg &
687 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
688 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
689 ,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
690 ,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
691 dz8w=dz8w, &
692 CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
693 ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
694 num_months=n_ozmixm, &
695 m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
696 aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
697 paerlev=paerlev, naer_c=n_aerosolc, &
698 cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
699 GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
700 SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
701 ,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
702 ,doabsems=doabsems &
703 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
704 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
705 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
706 )
707 ELSE
708 CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
709 ENDIF
710 CASE DEFAULT
711
712 WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
713 CALL wrf_error_fatal ( wrf_err_message )
714
715 END SELECT lwrad_select
716
717 IF (lw_physics .gt. 0 .and. .not.gfdl_lw) THEN
718 DO j=jts,jte
719 DO k=kts,kte
720 DO i=its,ite
721 RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
722 ! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY)
723 IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
724 ENDDO
725 ENDDO
726 ENDDO
727 ENDIF
728 !
729
730 swrad_select: SELECT CASE(sw_physics)
731
732 CASE (SWRADSCHEME)
733 CALL wrf_debug(100, 'CALL swrad')
734 CALL SWRAD( &
735 DT=dt,RTHRATEN=rthraten,GSW=gsw &
736 ,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
737 #ifdef WRF_CHEM
738 ,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
739 ,PM2_5_DRY_EC=pm2_5_dry_ec &
740 #endif
741 ,RHO_PHY=rho,T3D=t &
742 ,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
743 ,R=r_d,CP=cp,G=g,JULDAY=julday &
744 ,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
745 ! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
746 ,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
747 ,warm_rain=warm_rain &
748 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
749 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
750 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
751 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
752 ,QV3D=qv &
753 ,QC3D=qc &
754 ,QR3D=qr &
755 ,QI3D=qi &
756 ,QS3D=qs &
757 ,QG3D=qg &
758 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
759 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
760 )
761
762 CASE (GSFCSWSCHEME)
763 CALL wrf_debug(100, 'CALL gsfcswrad')
764 CALL GSFCSWRAD( &
765 RTHRATEN=rthraten,GSW=gsw,XLAT=xlat,XLONG=xlong &
766 ,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
767 ,DZ8W=dz8w,RHO_PHY=rho &
768 ,CLDFRA3D=cldfra,RSWTOA=rswtoa &
769 ,GMT=gmt,CP=cp,G=g &
770 ! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
771 ,JULDAY=julday,XTIME=xtime &
772 ,DECLIN=declin,SOLCON=solcon &
773 ,RADFRQ=radt,DEGRAD=degrad &
774 ,TAUCLDI=taucldi,TAUCLDC=taucldc &
775 ,WARM_RAIN=warm_rain &
776 #ifdef WRF_CHEM
777 ,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
778 ,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
779 ,GAER300=gaer300,GAER400=gaer400 & ! jcb
780 ,GAER600=gaer600,GAER999=gaer999 & ! jcb
781 ,WAER300=waer300,WAER400=waer400 & ! jcb
782 ,WAER600=waer600,WAER999=waer999 & ! jcb
783 ,aer_ra_feedback=aer_ra_feedback &
784 #endif
785 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
786 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
787 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
788 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
789 ,QV3D=qv &
790 ,QC3D=qc &
791 ,QR3D=qr &
792 ,QI3D=qi &
793 ,QS3D=qs &
794 ,QG3D=qg &
795 ,QNDROP3D=qndrop &
796 ,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
797 ,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
798 ,F_QNDROP=f_qndrop &
799 )
800 CASE (CAMSWSCHEME)
801 ! Temporarily lw switch already calculates sw CAM tendency, so inactive here
802
803 DO j=jts,jte
804 DO k=kts,kte
805 DO i=its,ite
806 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
807 ENDDO
808 ENDDO
809 ENDDO
810
811 CASE (GFDLSWSCHEME)
812
813 CALL wrf_debug (100, 'CALL gfdlsw')
814
815 IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
816 PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
817 PRESENT(qv) .AND. PRESENT(qc) ) THEN
818 IF ( F_QV .AND. F_QC .AND. F_QS ) THEN
819 gfdl_sw = .true.
820 CALL ETARA( &
821 DT=dt,XLAND=xland &
822 ,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
823 ,QV=qv,QW=qc_temp,QI=qi,QS=qs &
824 ,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
825 ,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
826 ,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
827 ,HBOTR=hbotr, HTOPR=htopr &
828 ,ALBEDO=albedo,CUPPT=cuppt &
829 ,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
830 ,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
831 ,XTIME=xtime,JULIAN=julian &
832 ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
833 ,JULYR=julyr,JULDAY=julday &
834 ,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
835 ,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
836 ,ACFRST=acfrst,NCFRST=ncfrst &
837 ,ACFRCV=acfrcv,NCFRCV=ncfrcv &
838 ,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
839 ,THRATEN=rthraten,THRATENLW=rthratenlw &
840 ,THRATENSW=rthratensw &
841 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
842 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
843 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
844 )
845 ELSE
846 CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
847 ENDIF
848 ELSE
849 CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
850 ENDIF
851
852 CASE DEFAULT
853
854 WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
855 CALL wrf_error_fatal ( wrf_err_message )
856
857 END SELECT swrad_select
858
859 IF (sw_physics .gt. 0 .and. .not.gfdl_sw) THEN
860 DO j=jts,jte
861 DO k=kts,kte
862 DO i=its,ite
863 RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
864 ENDDO
865 ENDDO
866 ENDDO
867
868 DO j=jts,jte
869 DO i=its,ite
870 SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
871 ENDDO
872 ENDDO
873
874 ENDIF
875
876 IF ( PRESENT( cu_rad_feedback ) ) THEN
877 IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
878 DO j=jts,jte
879 DO k=kts,kte
880 DO i=its,ite
881 qc(i,k,j) = qc_save(i,k,j)
882 ENDDO
883 ENDDO
884 ENDDO
885 ENDIF
886 IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
887 DO j=jts,jte
888 DO k=kts,kte
889 DO i=its,ite
890 qi(i,k,j) = qi_save(i,k,j)
891 ENDDO
892 ENDDO
893 ENDDO
894 ENDIF
895 ENDIF
896
897 ENDDO
898 !$OMP END PARALLEL DO
899
900 ENDIF Radiation_step
901
902 accumulate_lw_select: SELECT CASE(lw_physics)
903
904 CASE (CAMLWSCHEME)
905 IF(PRESENT(LWUPTC))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 ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DT
918 ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DT
919 ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DT
920 ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DT
921 ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DT
922 ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DT
923 ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DT
924 ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DT
925 ENDDO
926 ENDDO
927 ENDDO
928 !$OMP END PARALLEL DO
929 ENDIF
930 CASE DEFAULT
931 END SELECT accumulate_lw_select
932
933 accumulate_sw_select: SELECT CASE(sw_physics)
934
935 CASE (CAMSWSCHEME)
936 IF(PRESENT(SWUPTC))THEN
937 !$OMP PARALLEL DO &
938 !$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
939
940 DO ij = 1 , num_tiles
941 its = i_start(ij)
942 ite = i_end(ij)
943 jts = j_start(ij)
944 jte = j_end(ij)
945
946 DO j=jts,jte
947 DO i=its,ite
948 ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DT
949 ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DT
950 ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DT
951 ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DT
952 ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DT
953 ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DT
954 ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DT
955 ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DT
956 ENDDO
957 ENDDO
958 ENDDO
959 !$OMP END PARALLEL DO
960 ENDIF
961
962 CASE DEFAULT
963 END SELECT accumulate_sw_select
964
965 END SUBROUTINE radiation_driver
966
967 !---------------------------------------------------------------------
968 !BOP
969 ! !IROUTINE: radconst - compute radiation terms
970 ! !INTERFAC:
971 SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
972 DEGRAD,DPD )
973 !---------------------------------------------------------------------
974 USE module_wrf_error
975 IMPLICIT NONE
976 !---------------------------------------------------------------------
977
978 ! !ARGUMENTS:
979 REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
980 REAL, INTENT(OUT ) :: DECLIN,SOLCON
981 REAL :: OBECL,SINOB,SXLONG,ARG, &
982 DECDEG,DJUL,RJUL,ECCFAC
983 !
984 ! !DESCRIPTION:
985 ! Compute terms used in radiation physics
986 !EOP
987
988 ! for short wave radiation
989
990 DECLIN=0.
991 SOLCON=0.
992
993 !-----OBECL : OBLIQUITY = 23.5 DEGREE.
994
995 OBECL=23.5*DEGRAD
996 SINOB=SIN(OBECL)
997
998 !-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
999
1000 IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
1001 IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
1002 SXLONG=SXLONG*DEGRAD
1003 ARG=SINOB*SIN(SXLONG)
1004 DECLIN=ASIN(ARG)
1005 DECDEG=DECLIN/DEGRAD
1006 !----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
1007 DJUL=JULIAN*360./365.
1008 RJUL=DJUL*DEGRAD
1009 ECCFAC=1.000110+0.034221*COS(RJUL)+0.001280*SIN(RJUL)+0.000719* &
1010 COS(2*RJUL)+0.000077*SIN(2*RJUL)
1011 SOLCON=1370.*ECCFAC
1012
1013 !pjj/cray Cray X1 cannot print from threaded region
1014 #ifndef crayx1
1015 write(wrf_err_message,10)DECDEG,SOLCON
1016 10 FORMAT(1X,'*** SOLAR DECLINATION ANGLE = ',F6.2,' DEGREES.', &
1017 ' SOLAR CONSTANT = ',F8.2,' W/M**2 ***')
1018 CALL wrf_debug (50, wrf_err_message)
1019 #endif
1020
1021 END SUBROUTINE radconst
1022
1023 !---------------------------------------------------------------------
1024 !BOP
1025 ! !IROUTINE: cal_cldfra - Compute cloud fraction
1026 ! !INTERFACE:
1027 SUBROUTINE cal_cldfra(CLDFRA,QC,QI,F_QC,F_QI, &
1028 ids,ide, jds,jde, kds,kde, &
1029 ims,ime, jms,jme, kms,kme, &
1030 its,ite, jts,jte, kts,kte )
1031 !---------------------------------------------------------------------
1032 IMPLICIT NONE
1033 !---------------------------------------------------------------------
1034 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1035 ims,ime, jms,jme, kms,kme, &
1036 its,ite, jts,jte, kts,kte
1037
1038 !
1039 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1040 CLDFRA
1041
1042 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1043 QI, &
1044 QC
1045
1046 LOGICAL,INTENT(IN) :: F_QC,F_QI
1047
1048 REAL thresh
1049 INTEGER:: i,j,k
1050 ! !DESCRIPTION:
1051 ! Compute cloud fraction from input ice and cloud water fields
1052 ! if provided.
1053 !
1054 ! Whether QI or QC is active or not is determined from the indices of
1055 ! the fields into the 4D scalar arrays in WRF. These indices are
1056 ! P_QI and P_QC, respectively, and they are passed in to the routine
1057 ! to enable testing to see if QI and QC represent active fields in
1058 ! the moisture 4D scalar array carried by WRF.
1059 !
1060 ! If a field is active its index will have a value greater than or
1061 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1062 ! this routine.
1063 !EOP
1064 !---------------------------------------------------------------------
1065 thresh=1.0e-6
1066
1067 IF ( f_qi .AND. f_qc ) THEN
1068 DO j = jts,jte
1069 DO k = kts,kte
1070 DO i = its,ite
1071 IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
1072 CLDFRA(i,k,j)=1.
1073 ELSE
1074 CLDFRA(i,k,j)=0.
1075 ENDIF
1076 ENDDO
1077 ENDDO
1078 ENDDO
1079 ELSE IF ( f_qc ) THEN
1080 DO j = jts,jte
1081 DO k = kts,kte
1082 DO i = its,ite
1083 IF ( QC(i,k,j) .gt. thresh) THEN
1084 CLDFRA(i,k,j)=1.
1085 ELSE
1086 CLDFRA(i,k,j)=0.
1087 ENDIF
1088 ENDDO
1089 ENDDO
1090 ENDDO
1091 ELSE
1092 DO j = jts,jte
1093 DO k = kts,kte
1094 DO i = its,ite
1095 CLDFRA(i,k,j)=0.
1096 ENDDO
1097 ENDDO
1098 ENDDO
1099 ENDIF
1100
1101 END SUBROUTINE cal_cldfra
1102
1103 !BOP
1104 ! !IROUTINE: cal_cldfra2 - Compute cloud fraction
1105 ! !INTERFACE:
1106 ! cal_cldfra_xr - Compute cloud fraction.
1107 ! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
1108 !!
1109 !!--- Cloud fraction parameterization follows Randall, 1994
1110 !! (see Hong et al., 1998)
1111 !! (modified by Ferrier, Feb '02)
1112 !
1113 SUBROUTINE cal_cldfra2(CLDFRA, QV, QC, QI, QS, &
1114 F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
1115 F_ICE_PHY,F_RAIN_PHY, &
1116 ids,ide, jds,jde, kds,kde, &
1117 ims,ime, jms,jme, kms,kme, &
1118 its,ite, jts,jte, kts,kte )
1119 !---------------------------------------------------------------------
1120 IMPLICIT NONE
1121 !---------------------------------------------------------------------
1122 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1123 ims,ime, jms,jme, kms,kme, &
1124 its,ite, jts,jte, kts,kte
1125
1126 !
1127 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1128 CLDFRA
1129
1130 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1131 QV, &
1132 QI, &
1133 QC, &
1134 QS, &
1135 t_phy, &
1136 p_phy, &
1137 F_ICE_PHY, &
1138 F_RAIN_PHY
1139
1140 LOGICAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
1141
1142 ! REAL thresh
1143 INTEGER:: i,j,k
1144 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
1145
1146 REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
1147 PEXP=0.25, RHGRID=1.0
1148 REAL , PARAMETER :: SVP1=0.61078
1149 REAL , PARAMETER :: SVP2=17.2693882
1150 REAL , PARAMETER :: SVPI2=21.8745584
1151 REAL , PARAMETER :: SVP3=35.86
1152 REAL , PARAMETER :: SVPI3=7.66
1153 REAL , PARAMETER :: SVPT0=273.15
1154 REAL , PARAMETER :: r_d = 287.
1155 REAL , PARAMETER :: r_v = 461.6
1156 REAL , PARAMETER :: ep_2=r_d/r_v
1157 ! !DESCRIPTION:
1158 ! Compute cloud fraction from input ice and cloud water fields
1159 ! if provided.
1160 !
1161 ! Whether QI or QC is active or not is determined from the indices of
1162 ! the fields into the 4D scalar arrays in WRF. These indices are
1163 ! P_QI and P_QC, respectively, and they are passed in to the routine
1164 ! to enable testing to see if QI and QC represent active fields in
1165 ! the moisture 4D scalar array carried by WRF.
1166 !
1167 ! If a field is active its index will have a value greater than or
1168 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1169 ! this routine.
1170 !EOP
1171
1172
1173 !-----------------------------------------------------------------------
1174 !--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
1175 ! (modified by Ferrier, Feb '02)
1176 !
1177 !--- Cloud fraction parameterization follows Randall, 1994
1178 ! (see Hong et al., 1998)
1179 !-----------------------------------------------------------------------
1180 ! Note: ep_2=287./461.6 Rd/Rv
1181 ! Note: R_D=287.
1182
1183 ! Alternative calculation for critical RH for grid saturation
1184 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
1185
1186 ! Calculate saturation mixing ratio weighted according to the fractions of
1187 ! water and ice.
1188 ! Following:
1189 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
1190 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
1191 !
1192 ! over ice over water
1193 ! a = 21.8745584 17.2693882
1194 ! b = 7.66 35.86
1195
1196 !---------------------------------------------------------------------
1197
1198 DO j = jts,jte
1199 DO k = kts,kte
1200 DO i = its,ite
1201 tc = t_phy(i,k,j) - SVPT0
1202 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
1203 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
1204 QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
1205 QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
1206
1207 IF ( F_QI .and. F_QC .and. F_QS) THEN
1208 QCLD=QI(i,k,j)+QC(i,k,j)+QS(I,k,j)
1209 IF (QCLD .LT. QCLDMIN) THEN
1210 weight = 0.
1211 ELSE
1212 weight = (QI(i,k,j)+QS(I,k,j)) / QCLD
1213 ENDIF
1214 ELSE IF ( F_QC ) THEN
1215
1216 ! Mixing ratios of cloud water & total ice (cloud ice + snow).
1217 ! Mixing ratios of rain are not considered in this scheme.
1218 ! F_ICE is fraction of ice
1219 ! F_RAIN is fraction of rain
1220
1221 QIMID=QC(i,k,j)*F_ICE_PHY(i,k,j)
1222 QWMID=(QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
1223
1224
1225 !
1226 !--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
1227 ! only cloud water + cloud ice + snow
1228 !
1229 QCLD=QWMID+QIMID
1230 IF (QCLD .LT. QCLDMIN) THEN
1231 weight = 0.
1232 ELSE
1233 weight = F_ICE_PHY(i,k,j)
1234 ENDIF
1235
1236 ELSE
1237 CLDFRA(i,k,j)=0.
1238 ENDIF ! IF ( F_QI .and. F_QC )
1239
1240
1241 QVS_WEIGHT = (1-weight)*QVSW + weight*QVSI
1242 RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
1243 !
1244 !--- Determine cloud fraction (modified from original algorithm)
1245 !
1246 IF (QCLD .LT. QCLDMIN) THEN
1247 !
1248 !--- Assume zero cloud fraction if there is no cloud mixing ratio
1249 !
1250 CLDFRA(i,k,j)=0.
1251 ELSEIF(RHUM.GE.RHGRID)THEN
1252 !
1253 !--- Assume cloud fraction of unity if near saturation and the cloud
1254 ! mixing ratio is at or above the minimum threshold
1255 !
1256 CLDFRA(i,k,j)=1.
1257 ELSE
1258 !
1259 !--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
1260 ! modified based on assumed grid-scale saturation at RH=RHgrid.
1261 !
1262 SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
1263 DENOM=(SUBSAT)**GAMMA
1264 ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
1265 ! prevent negative values (new)
1266 RHUM=MAX(1.E-10, RHUM)
1267 CLDFRA(i,k,j)=(RHUM/RHGRID)**PEXP*(1.-EXP(ARG))
1268 !! ARG=-1000*QCLD/(RHUM-RHGRID)
1269 !! ARG=MAX(ARG, ARGMIN)
1270 !! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
1271 IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
1272 ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
1273 ENDDO !--- End DO i
1274 ENDDO !--- End DO k
1275 ENDDO !--- End DO j
1276
1277 END SUBROUTINE cal_cldfra2
1278
1279 END MODULE module_radiation_driver