Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_radiation_driver.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 54.2 KB

Line
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/(dtSTEPRA))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.034221COS(RJUL)+0.001280SIN(RJUL)+0.000719* &
1001	COS(2RJUL)+0.000077SIN(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/M2 *')
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 + weightQVSI
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format