Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_radiation_driver.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 80.8 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, CURR_SECS, adapt_step_flag &
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	,ht,dx,dy,sina,cosa,shadowmask,slope_rad ,topo_shading ) ! slope-dependent radiation
53
54	!-------------------------------------------------------------------------
55
56	! !USES:
57	USE module_state_description, ONLY : RRTMSCHEME, GFDLLWSCHEME &
58	,SWRADSCHEME, GSFCSWSCHEME &
59	,GFDLSWSCHEME, CAMLWSCHEME, CAMSWSCHEME &
60	,HELDSUAREZ
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	USE module_ra_hs
72
73	! This driver calls subroutines for the radiation parameterizations.
74	!
75	! short wave radiation choices:
76	! 1. swrad (19??)
77	!
78	! long wave radiation choices:
79	! 1. rrtmlwrad
80	!
81	!----------------------------------------------------------------------
82	IMPLICIT NONE
83	!<DESCRIPTION>
84	!
85	! Radiation_driver is the WRF mediation layer routine that provides the interface to
86	! to radiation physics packages in the WRF model layer. The radiation
87	! physics packages to call are chosen by setting the namelist variable
88	! (Rconfig entry in Registry) to the integer value assigned to the
89	! particular package (package entry in Registry). For example, if the
90	! namelist variable ra_lw_physics is set to 1, this corresponds to the
91	! Registry Package entry for swradscheme. Note that the Package
92	! names in the Registry are defined constants (frame/module_state_description.F)
93	! in the CASE statements in this routine.
94	!
95	! Among the arguments is moist, a four-dimensional scalar array storing
96	! a variable number of moisture tracers, depending on the physics
97	! configuration for the WRF run, as determined in the namelist. The
98	! highest numbered index of active moisture tracers the integer argument
99	! n_moist (note: the number of tracers at run time is the quantity
100	! <tt>n_moist - PARAM_FIRST_SCALAR + 1</tt> , not n_moist. Individual tracers
101	! may be indexed from moist by the Registry name of the tracer prepended
102	! with P_; for example P_QC is the index of cloud water. An index
103	! represents a valid, active field only if the index is greater than
104	! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual
105	! indices for each tracer is defined in module_state_description and
106	! set in <a href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a> defined in frame/module_configure.F.
107	!
108	! Physics drivers in WRF 2.0 and higher, originally model-layer
109	! routines, have been promoted to mediation layer routines and they
110	! contain OpenMP threaded loops over tiles. Thus, physics drivers
111	! are called from single-threaded regions in the solver. The physics
112	! routines that are called from the physics drivers are model-layer
113	! routines and fully tile-callable and thread-safe.
114	!</DESCRIPTION>
115	!
116	!======================================================================
117	! Grid structure in physics part of WRF
118	!----------------------------------------------------------------------
119	! The horizontal velocities used in the physics are unstaggered
120	! relative to temperature/moisture variables. All predicted
121	! variables are carried at half levels except w, which is at full
122	! levels. Some arrays with names (*8w) are at w (full) levels.
123	!
124	!----------------------------------------------------------------------
125	! In WRF, kms (smallest number) is the bottom level and kme (largest
126	! number) is the top level. In your scheme, if 1 is at the top level,
127	! then you have to reverse the order in the k direction.
128	!
129	! kme - half level (no data at this level)
130	! kme ----- full level
131	! kme-1 - half level
132	! kme-1 ----- full level
133	! .
134	! .
135	! .
136	! kms+2 - half level
137	! kms+2 ----- full level
138	! kms+1 - half level
139	! kms+1 ----- full level
140	! kms - half level
141	! kms ----- full level
142	!
143	!======================================================================
144	! Grid structure in physics part of WRF
145	!
146	!-------------------------------------
147	! The horizontal velocities used in the physics are unstaggered
148	! relative to temperature/moisture variables. All predicted
149	! variables are carried at half levels except w, which is at full
150	! levels. Some arrays with names (*8w) are at w (full) levels.
151	!
152	!==================================================================
153	! Definitions
154	!-----------
155	! Theta potential temperature (K)
156	! Qv water vapor mixing ratio (kg/kg)
157	! Qc cloud water mixing ratio (kg/kg)
158	! Qr rain water mixing ratio (kg/kg)
159	! Qi cloud ice mixing ratio (kg/kg)
160	! Qs snow mixing ratio (kg/kg)
161	!-----------------------------------------------------------------
162	!-- PM2_5_DRY Dry PM2.5 aerosol mass for all species (ug m^-3)
163	!-- PM2_5_WATER PM2.5 water mass (ug m^-3)
164	!-- PM2_5_DRY_EC Dry PM2.5 elemental carbon aersol mass (ug m^-3)
165	!-- RTHRATEN Theta tendency
166	! due to radiation (K/s)
167	!-- RTHRATENLW Theta tendency
168	! due to long wave radiation (K/s)
169	!-- RTHRATENSW Theta temperature tendency
170	! due to short wave radiation (K/s)
171	!-- dt time step (s)
172	!-- itimestep number of time steps
173	!-- GLW downward long wave flux at ground surface (W/m^2)
174	!-- GSW net short wave flux at ground surface (W/m^2)
175	!-- SWDOWN downward short wave flux at ground surface (W/m^2)
176	!-- SWDOWNC clear-sky downward short wave flux at ground surface (W/m^2; optional; for AQ)
177	!-- RLWTOA upward long wave at top of atmosphere (w/m2)
178	!-- RSWTOA upward short wave at top of atmosphere (w/m2)
179	!-- XLAT latitude, south is negative (degree)
180	!-- XLONG longitude, west is negative (degree)
181	!-- ALBEDO albedo (between 0 and 1)
182	!-- CLDFRA cloud fraction (between 0 and 1)
183	!-- EMISS surface emissivity (between 0 and 1)
184	!-- rho_phy density (kg/m^3)
185	!-- rr dry air density (kg/m^3)
186	!-- moist moisture array (4D - last index is species) (kg/kg)
187	!-- n_moist number of moisture species
188	!-- qndrop Cloud droplet number (#/kg)
189	!-- p8w pressure at full levels (Pa)
190	!-- p_phy pressure (Pa)
191	!-- Pb base-state pressure (Pa)
192	!-- pi_phy exner function (dimensionless)
193	!-- dz8w dz between full levels (m)
194	!-- t_phy temperature (K)
195	!-- t8w temperature at full levels (K)
196	!-- GMT Greenwich Mean Time Hour of model start (hour)
197	!-- JULDAY the initial day (Julian day)
198	!-- RADT time for calling radiation (min)
199	!-- ra_call_offset -1 (old) means usually just before output, 0 after
200	!-- DEGRAD conversion factor for
201	! degrees to radians (pi/180.) (rad/deg)
202	!-- DPD degrees per day for earth's
203	! orbital position (deg/day)
204	!-- R_d gas constant for dry air (J/kg/K)
205	!-- CP heat capacity at constant pressure for dry air (J/kg/K)
206	!-- G acceleration due to gravity (m/s^2)
207	!-- rvovrd R_v divided by R_d (dimensionless)
208	!-- XTIME time since simulation start (min)
209	!-- DECLIN solar declination angle (rad)
210	!-- SOLCON solar constant (W/m^2)
211	!-- ids start index for i in domain
212	!-- ide end index for i in domain
213	!-- jds start index for j in domain
214	!-- jde end index for j in domain
215	!-- kds start index for k in domain
216	!-- kde end index for k in domain
217	!-- ims start index for i in memory
218	!-- ime end index for i in memory
219	!-- jms start index for j in memory
220	!-- jme end index for j in memory
221	!-- kms start index for k in memory
222	!-- kme end index for k in memory
223	!-- i_start start indices for i in tile
224	!-- i_end end indices for i in tile
225	!-- j_start start indices for j in tile
226	!-- j_end end indices for j in tile
227	!-- kts start index for k in tile
228	!-- kte end index for k in tile
229	!-- num_tiles number of tiles
230	!
231	!==================================================================
232	!
233	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
234	ims,ime, jms,jme, kms,kme, &
235	kts,kte, &
236	num_tiles
237
238	INTEGER, INTENT(IN) :: lw_physics, sw_physics
239
240	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
241	i_start,i_end,j_start,j_end
242
243	INTEGER, INTENT(IN ) :: STEPRA,ICLOUD,ra_call_offset
244	INTEGER, INTENT(IN ) :: levsiz, n_ozmixm
245	INTEGER, INTENT(IN ) :: paerlev, n_aerosolc, cam_abs_dim1, cam_abs_dim2
246	REAL, INTENT(IN ) :: cam_abs_freq_s
247
248	LOGICAL, INTENT(IN ) :: warm_rain
249
250	REAL, INTENT(IN ) :: RADT
251
252	REAL, DIMENSION( ims:ime, jms:jme ), &
253	INTENT(IN ) :: XLAND, &
254	XICE, &
255	TSK, &
256	VEGFRA, &
257	SNOW
258	REAL, DIMENSION( ims:ime, levsiz, jms:jme, n_ozmixm ), OPTIONAL, &
259	INTENT(IN ) :: OZMIXM
260
261	REAL, DIMENSION(levsiz), OPTIONAL, INTENT(IN ) :: PIN
262
263	REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN ) :: m_ps_1,m_ps_2
264	REAL, DIMENSION( ims:ime, paerlev, jms:jme, n_aerosolc ), OPTIONAL, &
265	INTENT(IN ) :: aerosolc_1, aerosolc_2
266	REAL, DIMENSION(paerlev), OPTIONAL, &
267	INTENT(IN ) :: m_hybi0
268
269	REAL, DIMENSION( ims:ime, jms:jme ), &
270	INTENT(INOUT) :: HTOP, &
271	HBOT, &
272	HTOPR, &
273	HBOTR, &
274	CUPPT
275
276	INTEGER, INTENT(IN ) :: julyr
277	!
278	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
279	INTENT(IN ) :: dz8w, &
280	z, &
281	p8w, &
282	p, &
283	pi, &
284	t, &
285	t8w, &
286	rho
287	!
288	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
289	INTENT(IN ) :: tauaer300,tauaer400,tauaer600,tauaer999, & ! jcb
290	gaer300,gaer400,gaer600,gaer999, & ! jcb
291	waer300,waer400,waer600,waer999, & ! jcb
292	qc_adjust, qi_adjust
293
294	LOGICAL, OPTIONAL :: cu_rad_feedback
295
296	INTEGER, INTENT(IN ), OPTIONAL :: aer_ra_feedback
297
298	!
299	! variables for aerosols (only if running with chemistry)
300	!
301	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
302	INTENT(IN ) :: pm2_5_dry, &
303	pm2_5_water, &
304	pm2_5_dry_ec
305	!
306	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
307	INTENT(INOUT) :: RTHRATEN, &
308	RTHRATENLW, &
309	RTHRATENSW
310
311	! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL , &
312	! INTENT(INOUT) :: SWUP, &
313	! SWDN, &
314	! SWUPCLEAR, &
315	! SWDNCLEAR, &
316	! LWUP, &
317	! LWDN, &
318	! LWUPCLEAR, &
319	! LWDNCLEAR
320
321	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
322	ACSWUPT,ACSWUPTC,ACSWDNT,ACSWDNTC, &
323	ACSWUPB,ACSWUPBC,ACSWDNB,ACSWDNBC, &
324	ACLWUPT,ACLWUPTC,ACLWDNT,ACLWDNTC, &
325	ACLWUPB,ACLWUPBC,ACLWDNB,ACLWDNBC
326	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(INOUT) ::&
327	SWUPT, SWUPTC, SWDNT, SWDNTC, &
328	SWUPB, SWUPBC, SWDNB, SWDNBC, &
329	LWUPT, LWUPTC, LWDNT, LWDNTC, &
330	LWUPB, LWUPBC, LWDNB, LWDNBC
331
332	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL , &
333	INTENT(INOUT) :: SWCF, &
334	LWCF, &
335	OLR
336
337
338	!
339	REAL, DIMENSION( ims:ime, jms:jme ), &
340	INTENT(IN ) :: XLAT, &
341	XLONG, &
342	ALBEDO, &
343	EMISS
344	!
345	REAL, DIMENSION( ims:ime, jms:jme ), &
346	INTENT(INOUT) :: GSW, &
347	GLW
348
349	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: SWDOWN
350	!
351	REAL, INTENT(IN ) :: GMT,dt, &
352	julian, xtime
353	!
354	INTEGER, INTENT(IN ) :: JULDAY, itimestep
355	REAL, INTENT(IN ),OPTIONAL :: CURR_SECS
356	LOGICAL, INTENT(IN ),OPTIONAL :: ADAPT_STEP_FLAG
357
358	INTEGER,INTENT(IN) :: NPHS
359	REAL, DIMENSION( ims:ime, jms:jme ),INTENT(OUT) :: &
360	CFRACH, & !Added
361	CFRACL, & !Added
362	CFRACM, & !Added
363	CZMEAN !Added
364	REAL, DIMENSION( ims:ime, jms:jme ), &
365	INTENT(INOUT) :: &
366	RLWTOA, & !Added
367	RSWTOA, & !Added
368	ACFRST, & !Added
369	ACFRCV !Added
370
371	INTEGER,DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: &
372	NCFRST, & !Added
373	NCFRCV !Added
374	! Optional (only used by CAM lw scheme)
375
376	REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim2, jms:jme ), OPTIONAL ,&
377	INTENT(INOUT) :: abstot
378	REAL, DIMENSION( ims:ime, kms:kme, cam_abs_dim1, jms:jme ), OPTIONAL ,&
379	INTENT(INOUT) :: absnxt
380	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL ,&
381	INTENT(INOUT) :: emstot
382
383	!
384	! Optional
385	!
386	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
387	OPTIONAL, &
388	INTENT(INOUT) :: CLDFRA
389
390	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
391	OPTIONAL, &
392	INTENT(IN ) :: &
393	F_ICE_PHY, &
394	F_RAIN_PHY
395
396	REAL, DIMENSION( ims:ime, jms:jme ), &
397	OPTIONAL, &
398	INTENT(OUT) :: SWDOWNC
399	!
400	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
401	OPTIONAL, &
402	INTENT(INOUT ) :: &
403	pb &
404	,qv,qc,qr,qi,qs,qg,qndrop
405
406	LOGICAL, OPTIONAL :: f_qv,f_qc,f_qr,f_qi,f_qs,f_qg,f_qndrop
407	!
408	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
409	OPTIONAL, &
410	INTENT(INOUT) :: taucldi,taucldc
411
412	! Variables for slope-dependent radiation
413
414	REAL, OPTIONAL, INTENT(IN) :: dx,dy
415	INTEGER, OPTIONAL, INTENT(IN) :: slope_rad,topo_shading
416	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: sina,cosa,ht
417	INTEGER, DIMENSION( ims:ime, jms:jme ), OPTIONAL, INTENT(IN) :: shadowmask
418
419
420	! LOCAL VAR
421
422	REAL, DIMENSION( ims:ime, jms:jme ) :: GLAT,GLON
423	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: CEMISS
424	REAL, DIMENSION( ims:ime, jms:jme ) :: coszr
425
426	REAL :: DECLIN,SOLCON
427	INTEGER :: i,j,k,its,ite,jts,jte,ij
428	INTEGER :: STEPABS
429	LOGICAL :: gfdl_lw,gfdl_sw
430	LOGICAL :: doabsems
431	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
432
433	REAL :: OBECL,SINOB,SXLONG,ARG,DECDEG, &
434	DJUL,RJUL,ECCFAC
435	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_temp,qc_temp
436	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qi_save,qc_save
437
438	REAL :: next_rad_time
439	LOGICAL :: run_param
440	!------------------------------------------------------------------
441	! urban related variables are added to declaration
442	!-------------------------------------------------
443	REAL, OPTIONAL, INTENT(OUT) :: DECLIN_URB !urban
444	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: COSZ_URB2D !urban
445	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme), INTENT(OUT) :: omg_urb2d !urban
446	!------------------------------------------------------------------
447
448	if (lw_physics .eq. 0 .and. sw_physics .eq. 0) return
449
450	! ra_call_offset = -1 gives old method where radiation may be called just before output
451	! ra_call_offset = 0 gives new method where radiation may be called just after output
452	! and is also consistent with removal of offset in new XTIME
453	! also need to account for stepra=1 which always has zero modulo output
454
455	!
456	! Calculate whether or not to run the radiation step.
457	! If CURR_SECS is passed in, we will calculate based upon that. If it
458	! is not passed in, we'll do the old method of using the time STEP number
459	!
460
461	IF ( (itimestep .EQ. 1) .OR. (MOD(itimestep,STEPRA) .EQ. 1 + ra_call_offset) .OR. &
462	(STEPRA .EQ. 1) ) THEN
463	run_param = .TRUE.
464	ELSE
465	run_param = .FALSE.
466	ENDIF
467	IF (PRESENT(adapt_step_flag)) THEN
468	IF ((adapt_step_flag)) THEN
469	IF ( (itimestep .EQ. 1) .OR. (radt .EQ. 0) .OR. &
470	( CURR_SECS + dt >= ( INT( CURR_SECS / ( radt * 60 ) + 1 ) * radt * 60) ) ) THEN
471	run_param = .TRUE.
472	ELSE
473	run_param = .FALSE.
474	ENDIF
475	ENDIF
476	ENDIF
477
478	Radiation_step: IF ( run_param ) then
479
480	! CAM-specific additional radiation frequency - cam_abs_freq_s (=21600s by default)
481	STEPABS = nint(cam_abs_freq_s/(dtSTEPRA))STEPRA
482	IF (itimestep .eq. 1 .or. mod(itimestep,STEPABS) .eq. 1 + ra_call_offset &
483	.or. STEPABS .eq. 1 ) THEN
484	doabsems = .true.
485	ELSE
486	doabsems = .false.
487	ENDIF
488	IF (PRESENT(adapt_step_flag)) THEN
489	IF ((adapt_step_flag)) THEN
490	IF ( (itimestep .EQ. 1) .OR. (cam_abs_freq_s .EQ. 0) .OR. &
491	( CURR_SECS + dt >= ( INT( CURR_SECS / ( cam_abs_freq_s ) + 1 ) * cam_abs_freq_s) ) ) THEN
492	doabsems = .true.
493	ELSE
494	doabsems = .false.
495	ENDIF
496	ENDIF
497	ENDIF
498
499	gfdl_lw = .false.
500	gfdl_sw = .false.
501
502	!---------------
503	!$OMP PARALLEL DO &
504	!$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
505
506	DO ij = 1 , num_tiles
507	its = i_start(ij)
508	ite = i_end(ij)
509	jts = j_start(ij)
510	jte = j_end(ij)
511
512	! initialize data
513
514	DO j=jts,jte
515	DO i=its,ite
516	GSW(I,J)=0.
517	GLW(I,J)=0.
518	SWDOWN(I,J)=0.
519	GLAT(I,J)=XLAT(I,J)*DEGRAD
520	GLON(I,J)=XLONG(I,J)*DEGRAD
521	ENDDO
522	ENDDO
523
524	DO j=jts,jte
525	DO k=kts,kte+1
526	DO i=its,ite
527	RTHRATEN(I,K,J)=0.
528	! SWUP(I,K,J) = 0.0
529	! SWDN(I,K,J) = 0.0
530	! SWUPCLEAR(I,K,J) = 0.0
531	! SWDNCLEAR(I,K,J) = 0.0
532	! LWUP(I,K,J) = 0.0
533	! LWDN(I,K,J) = 0.0
534	! LWUPCLEAR(I,K,J) = 0.0
535	! LWDNCLEAR(I,K,J) = 0.0
536	CEMISS(I,K,J)=0.0
537	ENDDO
538	ENDDO
539	ENDDO
540
541	! temporarily modify hydrometeors (currently only done for GD scheme and WRF-Chem)
542	!
543	IF ( PRESENT( cu_rad_feedback ) ) THEN
544	IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
545	DO j=jts,jte
546	DO k=kts,kte
547	DO i=its,ite
548	qc_save(i,k,j) = qc(i,k,j)
549	qc(i,k,j) = qc(i,k,j) + qc_adjust(i,k,j)
550	ENDDO
551	ENDDO
552	ENDDO
553	ENDIF
554	IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
555	DO j=jts,jte
556	DO k=kts,kte
557	DO i=its,ite
558	qi_save(i,k,j) = qi(i,k,j)
559	qi(i,k,j) = qi(i,k,j) + qi_adjust(i,k,j)
560	ENDDO
561	ENDDO
562	ENDDO
563	ENDIF
564	ENDIF
565
566
567	! Fill temporary water variable depending on micro package (tgs 25 Apr 2006)
568	if(PRESENT(qc) .and. PRESENT(F_QC)) then
569	DO j=jts,jte
570	DO k=kts,kte
571	DO i=its,ite
572	qc_temp(I,K,J)=qc(I,K,J)
573	ENDDO
574	ENDDO
575	ENDDO
576	else
577	DO j=jts,jte
578	DO k=kts,kte
579	DO i=its,ite
580	qc_temp(I,K,J)=0.
581	ENDDO
582	ENDDO
583	ENDDO
584	endif
585	if(PRESENT(qr) .and. PRESENT(F_QR)) then
586	DO j=jts,jte
587	DO k=kts,kte
588	DO i=its,ite
589	qc_temp(I,K,J) = qc_temp(I,K,J) + qr(I,K,J)
590	ENDDO
591	ENDDO
592	ENDDO
593	endif
594
595	!---------------
596	! Calculate constant for short wave radiation
597
598	CALL radconst(XTIME,DECLIN,SOLCON,JULIAN, &
599	DEGRAD,DPD )
600
601
602	if(present(DECLIN_URB))DECLIN_URB=DECLIN ! urban
603
604	lwrad_cldfra_select: SELECT CASE(lw_physics)
605
606	CASE (GFDLLWSCHEME)
607
608	!-- Do nothing, since cloud fractions (with partial cloudiness effects)
609	!-- are defined in GFDL LW/SW schemes and do not need to be initialized.
610
611	CASE (CAMLWSCHEME)
612
613	IF ( PRESENT ( CLDFRA ) .AND. &
614	PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
615	! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
616
617	CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
618	F_QV,F_QC,F_QI,F_QS,t,p, &
619	F_ICE_PHY,F_RAIN_PHY, &
620	ids,ide, jds,jde, kds,kde, &
621	ims,ime, jms,jme, kms,kme, &
622	its,ite, jts,jte, kts,kte )
623	ENDIF
624
625	CASE DEFAULT
626
627	IF ( PRESENT ( CLDFRA ) .AND. &
628	PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
629	CALL cal_cldfra(CLDFRA,qc,qi,F_QC,F_QI, &
630	ids,ide, jds,jde, kds,kde, &
631	ims,ime, jms,jme, kms,kme, &
632	its,ite, jts,jte, kts,kte )
633	ENDIF
634
635	END SELECT lwrad_cldfra_select
636
637	lwrad_select: SELECT CASE(lw_physics)
638
639
640	CASE (RRTMSCHEME)
641	CALL wrf_debug (100, 'CALL rrtm')
642
643	CALL RRTMLWRAD( &
644	RTHRATEN=RTHRATEN,GLW=GLW,OLR=RLWTOA,EMISS=EMISS &
645	,QV3D=QV &
646	,QC3D=QC &
647	,QR3D=QR &
648	,QI3D=QI &
649	,QS3D=QS &
650	,QG3D=QG &
651	,P8W=p8w,P3D=p,PI3D=pi,DZ8W=dz8w,TSK=tsk,T3D=t &
652	,T8W=t8w,RHO3D=rho, CLDFRA3D=CLDFRA,R=R_d,G=G &
653	,F_QV=F_QV,F_QC=F_QC,F_QR=F_QR &
654	,F_QI=F_QI,F_QS=F_QS,F_QG=F_QG &
655	,ICLOUD=icloud,WARM_RAIN=warm_rain &
656	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
657	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
658	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
659	)
660
661	CASE (GFDLLWSCHEME)
662
663	CALL wrf_debug (100, 'CALL gfdllw')
664
665	IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
666	PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
667	PRESENT(qv) .AND. PRESENT(qc) ) THEN
668	IF ( F_QV .AND. F_QC .AND. F_QS) THEN
669	gfdl_lw = .true.
670	CALL ETARA( &
671	DT=dt,XLAND=xland &
672	,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
673	,QV=qv,QW=qc_temp,QI=qi,QS=qs &
674	,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
675	,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
676	,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
677	,HBOTR=hbotr, HTOPR=htopr &
678	,ALBEDO=albedo,CUPPT=cuppt &
679	,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
680	,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
681	,XTIME=xtime,JULIAN=julian &
682	,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
683	,JULYR=julyr,JULDAY=julday &
684	,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
685	,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
686	,ACFRST=acfrst,NCFRST=ncfrst &
687	,ACFRCV=acfrcv,NCFRCV=ncfrcv &
688	,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
689	,THRATEN=rthraten,THRATENLW=rthratenlw &
690	,THRATENSW=rthratensw &
691	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
692	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
693	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
694	)
695	ELSE
696	CALL wrf_error_fatal('Can not call ETARA (1a). Missing moisture fields.')
697	ENDIF
698	ELSE
699	CALL wrf_error_fatal('Can not call ETARA (1b). Missing moisture fields.')
700	ENDIF
701	CASE (CAMLWSCHEME)
702	CALL wrf_debug(100, 'CALL camrad lw')
703	IF(cam_abs_dim1 .ne. 4 .or. cam_abs_dim2 .ne. kde .or. &
704	paerlev .ne. 29 .or. levsiz .ne. 59 )THEN
705	WRITE( wrf_err_message , * ) &
706	'set paerlev=29, levsiz=59, cam_abs_dim1=4, and cam_abs_dim2=number of levels (e_vert) in physics namelist for CAM radiation'
707	CALL wrf_error_fatal ( wrf_err_message )
708	ENDIF
709	IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
710	PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
711	PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
712	.AND. PRESENT(AEROSOLC_2) ) THEN
713	CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
714	dolw=.true.,dosw=.false., &
715	SWUPT=SWUPT,SWUPTC=SWUPTC, &
716	SWDNT=SWDNT,SWDNTC=SWDNTC, &
717	LWUPT=LWUPT,LWUPTC=LWUPTC, &
718	LWDNT=LWDNT,LWDNTC=LWDNTC, &
719	SWUPB=SWUPB,SWUPBC=SWUPBC, &
720	SWDNB=SWDNB,SWDNBC=SWDNBC, &
721	LWUPB=LWUPB,LWUPBC=LWUPBC, &
722	LWDNB=LWDNB,LWDNBC=LWDNBC, &
723	SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
724	TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
725	GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
726	ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
727	,QV3D=qv &
728	,QC3D=qc &
729	,QR3D=qr &
730	,QI3D=qi &
731	,QS3D=qs &
732	,QG3D=qg &
733	,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
734	,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
735	,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
736	,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
737	dz8w=dz8w, &
738	CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
739	ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
740	num_months=n_ozmixm, &
741	m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
742	aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
743	paerlev=paerlev, naer_c=n_aerosolc, &
744	cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
745	GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
746	SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
747	,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
748	,doabsems=doabsems &
749	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
750	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
751	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
752	)
753	ELSE
754	CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
755	ENDIF
756
757
758	CASE (HELDSUAREZ)
759	CALL wrf_debug (100, 'CALL heldsuarez')
760
761	CALL HSRAD(RTHRATEN,p8w,p,pi,dz8w,t, &
762	t8w, rho, R_d,G,CP, dt, xlat, degrad, &
763	ids,ide, jds,jde, kds,kde, &
764	ims,ime, jms,jme, kms,kme, &
765	its,ite, jts,jte, kts,kte )
766
767	CASE DEFAULT
768
769	WRITE( wrf_err_message , * ) 'The longwave option does not exist: lw_physics = ', lw_physics
770	CALL wrf_error_fatal ( wrf_err_message )
771
772	END SELECT lwrad_select
773
774	IF (lw_physics .gt. 0 .and. .not.gfdl_lw) THEN
775	DO j=jts,jte
776	DO k=kts,kte
777	DO i=its,ite
778	RTHRATENLW(I,K,J)=RTHRATEN(I,K,J)
779	! OLR ALSO WILL CONTAIN OUTGOING LONGWAVE FOR RRTM (NMM HAS NO OLR ARRAY)
780	IF(PRESENT(OLR) .AND. K .EQ. 1)OLR(I,J)=RLWTOA(I,J)
781	ENDDO
782	ENDDO
783	ENDDO
784	ENDIF
785	!
786	swrad_cldfra_select: SELECT CASE(sw_physics)
787
788	CASE (CAMSWSCHEME)
789
790	IF ( PRESENT ( CLDFRA ) .AND. &
791	PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
792	! Call to cloud fraction routine based on Randall 1994 (Hong Pan 1998)
793
794	CALL cal_cldfra2(CLDFRA,qv,qc,qi,qs, &
795	F_QV,F_QC,F_QI,F_QS,t,p, &
796	F_ICE_PHY,F_RAIN_PHY, &
797	ids,ide, jds,jde, kds,kde, &
798	ims,ime, jms,jme, kms,kme, &
799	its,ite, jts,jte, kts,kte )
800	ENDIF
801
802	CASE DEFAULT
803
804	END SELECT swrad_cldfra_select
805
806	swrad_select: SELECT CASE(sw_physics)
807
808	CASE (SWRADSCHEME)
809	CALL wrf_debug(100, 'CALL swrad')
810	CALL SWRAD( &
811	DT=dt,RTHRATEN=rthraten,GSW=gsw &
812	,XLAT=xlat,XLONG=xlong,ALBEDO=albedo &
813	#ifdef WRF_CHEM
814	,PM2_5_DRY=pm2_5_dry,PM2_5_WATER=pm2_5_water &
815	,PM2_5_DRY_EC=pm2_5_dry_ec &
816	#endif
817	,RHO_PHY=rho,T3D=t &
818	,P3D=p,PI3D=pi,DZ8W=dz8w,GMT=gmt &
819	,R=r_d,CP=cp,G=g,JULDAY=julday &
820	,XTIME=xtime,DECLIN=declin,SOLCON=solcon &
821	! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
822	,RADFRQ=radt,ICLOUD=icloud,DEGRAD=degrad &
823	,warm_rain=warm_rain &
824	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
825	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
826	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
827	,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
828	,QV3D=qv &
829	,QC3D=qc &
830	,QR3D=qr &
831	,QI3D=qi &
832	,QS3D=qs &
833	,QG3D=qg &
834	,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
835	,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
836	,slope_rad=slope_rad,topo_shading=topo_shading &
837	,shadowmask=shadowmask &
838	,ht=ht,dx=dx,dy=dy,sina=sina,cosa=cosa )
839
840	CASE (GSFCSWSCHEME)
841	CALL wrf_debug(100, 'CALL gsfcswrad')
842	CALL GSFCSWRAD( &
843	RTHRATEN=rthraten,GSW=gsw,XLAT=xlat,XLONG=xlong &
844	,ALB=albedo,T3D=t,P3D=p,P8W3D=p8w,pi3D=pi &
845	,DZ8W=dz8w,RHO_PHY=rho &
846	,CLDFRA3D=cldfra,RSWTOA=rswtoa &
847	,GMT=gmt,CP=cp,G=g &
848	! ,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
849	,JULDAY=julday,XTIME=xtime &
850	,DECLIN=declin,SOLCON=solcon &
851	,RADFRQ=radt,DEGRAD=degrad &
852	,TAUCLDI=taucldi,TAUCLDC=taucldc &
853	,WARM_RAIN=warm_rain &
854	#ifdef WRF_CHEM
855	,TAUAER300=tauaer300,TAUAER400=tauaer400 & ! jcb
856	,TAUAER600=tauaer600,TAUAER999=tauaer999 & ! jcb
857	,GAER300=gaer300,GAER400=gaer400 & ! jcb
858	,GAER600=gaer600,GAER999=gaer999 & ! jcb
859	,WAER300=waer300,WAER400=waer400 & ! jcb
860	,WAER600=waer600,WAER999=waer999 & ! jcb
861	,aer_ra_feedback=aer_ra_feedback &
862	#endif
863	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
864	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
865	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
866	,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d & !urban
867	,QV3D=qv &
868	,QC3D=qc &
869	,QR3D=qr &
870	,QI3D=qi &
871	,QS3D=qs &
872	,QG3D=qg &
873	,QNDROP3D=qndrop &
874	,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
875	,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
876	,F_QNDROP=f_qndrop &
877	)
878	CASE (CAMSWSCHEME)
879	CALL wrf_debug(100, 'CALL camrad sw')
880	IF(cam_abs_dim1 .ne. 4 .or. cam_abs_dim2 .ne. kde .or. &
881	paerlev .ne. 29 .or. levsiz .ne. 59 )THEN
882	WRITE( wrf_err_message , * ) &
883	'set paerlev=29, levsiz=59, cam_abs_dim1=4, and cam_abs_dim2=number of levels (e_vert) in physics namelist for CAM radiation'
884	CALL wrf_error_fatal ( wrf_err_message )
885	ENDIF
886	IF ( PRESENT( OZMIXM ) .AND. PRESENT( PIN ) .AND. &
887	PRESENT(M_PS_1) .AND. PRESENT(M_PS_2) .AND. &
888	PRESENT(M_HYBI0) .AND. PRESENT(AEROSOLC_1) &
889	.AND. PRESENT(AEROSOLC_2) ) THEN
890	CALL CAMRAD(RTHRATENLW=RTHRATEN,RTHRATENSW=RTHRATENSW, &
891	dolw=.false.,dosw=.true., &
892	SWUPT=SWUPT,SWUPTC=SWUPTC, &
893	SWDNT=SWDNT,SWDNTC=SWDNTC, &
894	LWUPT=LWUPT,LWUPTC=LWUPTC, &
895	LWDNT=LWDNT,LWDNTC=LWDNTC, &
896	SWUPB=SWUPB,SWUPBC=SWUPBC, &
897	SWDNB=SWDNB,SWDNBC=SWDNBC, &
898	LWUPB=LWUPB,LWUPBC=LWUPBC, &
899	LWDNB=LWDNB,LWDNBC=LWDNBC, &
900	SWCF=SWCF,LWCF=LWCF,OLR=RLWTOA,CEMISS=CEMISS, &
901	TAUCLDC=TAUCLDC,TAUCLDI=TAUCLDI,COSZR=COSZR, &
902	GSW=GSW,GLW=GLW,XLAT=XLAT,XLONG=XLONG, &
903	ALBEDO=ALBEDO,t_phy=t,TSK=TSK,EMISS=EMISS &
904	,QV3D=qv &
905	,QC3D=qc &
906	,QR3D=qr &
907	,QI3D=qi &
908	,QS3D=qs &
909	,QG3D=qg &
910	,F_QV=f_qv,F_QC=f_qc,F_QR=f_qr &
911	,F_QI=f_qi,F_QS=f_qs,F_QG=f_qg &
912	,f_ice_phy=f_ice_phy,f_rain_phy=f_rain_phy &
913	,p_phy=p,p8w=p8w,z=z,pi_phy=pi,rho_phy=rho, &
914	dz8w=dz8w, &
915	CLDFRA=CLDFRA,XLAND=XLAND,XICE=XICE,SNOW=SNOW, &
916	ozmixm=ozmixm,pin0=pin,levsiz=levsiz, &
917	num_months=n_ozmixm, &
918	m_psp=m_ps_1,m_psn=m_ps_2,aerosolcp=aerosolc_1, &
919	aerosolcn=aerosolc_2,m_hybi0=m_hybi0, &
920	paerlev=paerlev, naer_c=n_aerosolc, &
921	cam_abs_dim1=cam_abs_dim1, cam_abs_dim2=cam_abs_dim2, &
922	GMT=GMT,JULDAY=JULDAY,JULIAN=JULIAN,DT=DT,XTIME=XTIME,DECLIN=DECLIN, &
923	SOLCON=SOLCON,RADT=RADT,DEGRAD=DEGRAD,n_cldadv=3 &
924	,abstot_3d=abstot,absnxt_3d=absnxt,emstot_3d=emstot &
925	,doabsems=doabsems &
926	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
927	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
928	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
929	)
930	ELSE
931	CALL wrf_error_fatal ( 'arguments not present for calling cam radiation' )
932	ENDIF
933	DO j=jts,jte
934	DO k=kts,kte
935	DO i=its,ite
936	RTHRATEN(I,K,J)=RTHRATEN(I,K,J)+RTHRATENSW(I,K,J)
937	ENDDO
938	ENDDO
939	ENDDO
940
941	CASE (GFDLSWSCHEME)
942
943	CALL wrf_debug (100, 'CALL gfdlsw')
944
945	IF ( PRESENT(F_QV) .AND. PRESENT(F_QC) .AND. &
946	PRESENT(F_QS) .AND. PRESENT(qs) .AND. &
947	PRESENT(qv) .AND. PRESENT(qc) ) THEN
948	IF ( F_QV .AND. F_QC .AND. F_QS ) THEN
949	gfdl_sw = .true.
950	CALL ETARA( &
951	DT=dt,XLAND=xland &
952	,P8W=p8w,DZ8W=dz8w,RHO_PHY=rho,P_PHY=p,T=t &
953	,QV=qv,QW=qc_temp,QI=qi,QS=qs &
954	,TSK2D=tsk,GLW=GLW,RSWIN=SWDOWN,GSW=GSW &
955	,RSWINC=SWDOWNC,CLDFRA=CLDFRA,PI3D=pi &
956	,GLAT=glat,GLON=glon,HTOP=htop,HBOT=hbot &
957	,HBOTR=hbotr, HTOPR=htopr &
958	,ALBEDO=albedo,CUPPT=cuppt &
959	,VEGFRA=vegfra,SNOW=snow,G=g,GMT=gmt &
960	,NSTEPRA=stepra,NPHS=nphs,ITIMESTEP=itimestep &
961	,XTIME=xtime,JULIAN=julian &
962	,COSZ_URB2D=COSZ_URB2D ,OMG_URB2D=omg_urb2d &
963	,JULYR=julyr,JULDAY=julday &
964	,GFDL_LW=gfdl_lw,GFDL_SW=gfdl_sw &
965	,CFRACL=cfracl,CFRACM=cfracm,CFRACH=cfrach &
966	,ACFRST=acfrst,NCFRST=ncfrst &
967	,ACFRCV=acfrcv,NCFRCV=ncfrcv &
968	,RSWTOA=rswtoa,RLWTOA=rlwtoa,CZMEAN=czmean &
969	,THRATEN=rthraten,THRATENLW=rthratenlw &
970	,THRATENSW=rthratensw &
971	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
972	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
973	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
974	)
975	ELSE
976	CALL wrf_error_fatal('Can not call ETARA (2a). Missing moisture fields.')
977	ENDIF
978	ELSE
979	CALL wrf_error_fatal('Can not call ETARA (2b). Missing moisture fields.')
980	ENDIF
981
982	CASE (0)
983
984	! Here in case we don't want to call a sw radiation scheme
985	! For example, the Held-Suarez idealized test case
986	IF (lw_physics /= HELDSUAREZ) THEN
987	WRITE( wrf_err_message , * ) &
988	'You have selected a longwave radiation option, but not a shortwave option (sw_physics = 0, lw_physics = ',lw_physics,')'
989	CALL wrf_error_fatal ( wrf_err_message )
990	END IF
991
992	CASE DEFAULT
993
994	WRITE( wrf_err_message , * ) 'The shortwave option does not exist: sw_physics = ', sw_physics
995	CALL wrf_error_fatal ( wrf_err_message )
996
997	END SELECT swrad_select
998
999	IF (sw_physics .gt. 0 .and. .not.gfdl_sw) THEN
1000	DO j=jts,jte
1001	DO k=kts,kte
1002	DO i=its,ite
1003	RTHRATENSW(I,K,J)=RTHRATEN(I,K,J)-RTHRATENLW(I,K,J)
1004	ENDDO
1005	ENDDO
1006	ENDDO
1007
1008	DO j=jts,jte
1009	DO i=its,ite
1010	SWDOWN(I,J)=GSW(I,J)/(1.-ALBEDO(I,J))
1011	ENDDO
1012	ENDDO
1013
1014	ENDIF
1015
1016	IF ( PRESENT( cu_rad_feedback ) ) THEN
1017	IF ( PRESENT( qc ) .AND. PRESENT( qc_adjust ) .AND. cu_rad_feedback ) THEN
1018	DO j=jts,jte
1019	DO k=kts,kte
1020	DO i=its,ite
1021	qc(i,k,j) = qc_save(i,k,j)
1022	ENDDO
1023	ENDDO
1024	ENDDO
1025	ENDIF
1026	IF ( PRESENT( qi ) .AND. PRESENT( qi_adjust ) .AND. cu_rad_feedback ) THEN
1027	DO j=jts,jte
1028	DO k=kts,kte
1029	DO i=its,ite
1030	qi(i,k,j) = qi_save(i,k,j)
1031	ENDDO
1032	ENDDO
1033	ENDDO
1034	ENDIF
1035	ENDIF
1036
1037	ENDDO
1038	!$OMP END PARALLEL DO
1039
1040	ENDIF Radiation_step
1041
1042	accumulate_lw_select: SELECT CASE(lw_physics)
1043
1044	CASE (CAMLWSCHEME)
1045	IF(PRESENT(LWUPTC))THEN
1046	!$OMP PARALLEL DO &
1047	!$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1048
1049	DO ij = 1 , num_tiles
1050	its = i_start(ij)
1051	ite = i_end(ij)
1052	jts = j_start(ij)
1053	jte = j_end(ij)
1054
1055	DO j=jts,jte
1056	DO i=its,ite
1057	ACLWUPT(I,J) = ACLWUPT(I,J) + LWUPT(I,J)*DT
1058	ACLWUPTC(I,J) = ACLWUPTC(I,J) + LWUPTC(I,J)*DT
1059	ACLWDNT(I,J) = ACLWDNT(I,J) + LWDNT(I,J)*DT
1060	ACLWDNTC(I,J) = ACLWDNTC(I,J) + LWDNTC(I,J)*DT
1061	ACLWUPB(I,J) = ACLWUPB(I,J) + LWUPB(I,J)*DT
1062	ACLWUPBC(I,J) = ACLWUPBC(I,J) + LWUPBC(I,J)*DT
1063	ACLWDNB(I,J) = ACLWDNB(I,J) + LWDNB(I,J)*DT
1064	ACLWDNBC(I,J) = ACLWDNBC(I,J) + LWDNBC(I,J)*DT
1065	ENDDO
1066	ENDDO
1067	ENDDO
1068	!$OMP END PARALLEL DO
1069	ENDIF
1070	CASE DEFAULT
1071	END SELECT accumulate_lw_select
1072
1073	accumulate_sw_select: SELECT CASE(sw_physics)
1074
1075	CASE (CAMSWSCHEME)
1076	IF(PRESENT(SWUPTC))THEN
1077	!$OMP PARALLEL DO &
1078	!$OMP PRIVATE ( ij ,i,j,k,its,ite,jts,jte)
1079
1080	DO ij = 1 , num_tiles
1081	its = i_start(ij)
1082	ite = i_end(ij)
1083	jts = j_start(ij)
1084	jte = j_end(ij)
1085
1086	DO j=jts,jte
1087	DO i=its,ite
1088	ACSWUPT(I,J) = ACSWUPT(I,J) + SWUPT(I,J)*DT
1089	ACSWUPTC(I,J) = ACSWUPTC(I,J) + SWUPTC(I,J)*DT
1090	ACSWDNT(I,J) = ACSWDNT(I,J) + SWDNT(I,J)*DT
1091	ACSWDNTC(I,J) = ACSWDNTC(I,J) + SWDNTC(I,J)*DT
1092	ACSWUPB(I,J) = ACSWUPB(I,J) + SWUPB(I,J)*DT
1093	ACSWUPBC(I,J) = ACSWUPBC(I,J) + SWUPBC(I,J)*DT
1094	ACSWDNB(I,J) = ACSWDNB(I,J) + SWDNB(I,J)*DT
1095	ACSWDNBC(I,J) = ACSWDNBC(I,J) + SWDNBC(I,J)*DT
1096	ENDDO
1097	ENDDO
1098	ENDDO
1099	!$OMP END PARALLEL DO
1100	ENDIF
1101
1102	CASE DEFAULT
1103	END SELECT accumulate_sw_select
1104
1105	END SUBROUTINE radiation_driver
1106
1107	SUBROUTINE pre_radiation_driver ( grid, config_flags &
1108	,itimestep, ra_call_offset &
1109	,XLAT, XLONG, GMT, julian, xtime, RADT, STEPRA &
1110	,ht,dx,dy,sina,cosa,shadowmask,slope_rad ,topo_shading &
1111	,shadlen,ht_shad,ht_loc &
1112	,ht_shad_bxs, ht_shad_bxe &
1113	,ht_shad_bys, ht_shad_bye &
1114	,nested, min_ptchsz &
1115	,spec_bdy_width &
1116	,ids, ide, jds, jde, kds, kde &
1117	,ims, ime, jms, jme, kms, kme &
1118	,ips, ipe, jps, jpe, kps, kpe &
1119	,i_start, i_end &
1120	,j_start, j_end &
1121	,kts, kte &
1122	,num_tiles )
1123
1124	USE module_dm
1125	USE module_domain
1126	USE module_bc
1127	USE module_model_constants
1128
1129	IMPLICIT NONE
1130
1131	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1132	ims,ime, jms,jme, kms,kme, &
1133	ips,ipe, jps,jpe, kps,kpe, &
1134	kts,kte, &
1135	num_tiles
1136
1137	TYPE(domain) , INTENT(INOUT) :: grid
1138	TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
1139
1140	INTEGER, INTENT(IN ) :: itimestep, ra_call_offset, stepra, &
1141	slope_rad, topo_shading, &
1142	spec_bdy_width
1143
1144	INTEGER, INTENT(INOUT) :: min_ptchsz
1145
1146	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
1147	i_start,i_end,j_start,j_end
1148
1149	REAL, INTENT(IN ) :: GMT, radt, julian, xtime, dx, dy, shadlen
1150
1151	REAL, DIMENSION( ims:ime, jms:jme ), &
1152	INTENT(IN ) :: XLAT, &
1153	XLONG, &
1154	HT, &
1155	SINA, &
1156	COSA
1157
1158	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
1159
1160	REAL, DIMENSION( jms:jme , kds:kde , spec_bdy_width ), &
1161	INTENT(IN ) :: ht_shad_bxs, ht_shad_bxe
1162	REAL, DIMENSION( ims:ime , kds:kde , spec_bdy_width ), &
1163	INTENT(IN ) :: ht_shad_bys, ht_shad_bye
1164
1165	INTEGER, DIMENSION( ims:ime, jms:jme ), &
1166	INTENT(INOUT) :: shadowmask
1167
1168	LOGICAL, INTENT(IN ) :: nested
1169
1170	!Local
1171	! For orographic shading
1172	INTEGER :: niter,ni,psx,psy,idum,jdum,i,j,ij
1173	REAL :: DECLIN,SOLCON
1174
1175	! Determine minimum patch size for slope-dependent radiation
1176	if (itimestep .eq. 1) then
1177	psx = ipe-ips+1
1178	psy = jpe-jps+1
1179	min_ptchsz = min(psx,psy)
1180	idum = 0
1181	jdum = 0
1182	endif
1183
1184	# ifdef DM_PARALLEL
1185	if (itimestep .eq. 1) then
1186	call wrf_dm_minval_integer (psx,idum,jdum)
1187	call wrf_dm_minval_integer (psy,idum,jdum)
1188	min_ptchsz = min(psx,psy)
1189	endif
1190	# endif
1191
1192	! Topographic shading
1193
1194	if ((topo_shading.eq.1).and.(itimestep .eq. 1 .or. &
1195	mod(itimestep,STEPRA) .eq. 1 + ra_call_offset)) then
1196
1197	!---------------
1198	! Calculate constants for short wave radiation
1199
1200	CALL radconst(XTIME,DECLIN,SOLCON,JULIAN,DEGRAD,DPD)
1201
1202	! Make a local copy of terrain height field
1203	do j=jms,jme
1204	do i=ims,ime
1205	ht_loc(i,j) = ht(i,j)
1206	enddo
1207	enddo
1208	! Determine if iterations are necessary for shadows to propagate from one patch to another
1209	if ((ids.eq.ips).and.(ide.eq.ipe).and.(jds.eq.jps).and.(jde.eq.jpe)) then
1210	niter = 1
1211	else
1212	niter = int(shadlen/(dx*min_ptchsz)+3)
1213	endif
1214
1215
1216
1217	IF( nested ) THEN
1218
1219	!$OMP PARALLEL DO &
1220	!$OMP PRIVATE ( ij )
1221
1222	DO ij = 1 , num_tiles
1223
1224	CALL spec_bdyfield(ht_shad, &
1225	ht_shad_bxs, ht_shad_bxe, &
1226	ht_shad_bys, ht_shad_bye, &
1227	'm', config_flags, spec_bdy_width, 2,&
1228	ids,ide, jds,jde, 1 ,1 , & ! domain dims
1229	ims,ime, jms,jme, 1 ,1 , & ! memory dims
1230	ips,ipe, jps,jpe, 1 ,1 , & ! patch dims
1231	i_start(ij), i_end(ij), &
1232	j_start(ij), j_end(ij), &
1233	1 , 1 )
1234	ENDDO
1235	ENDIF
1236
1237	do ni = 1, niter
1238
1239	!$OMP PARALLEL DO &
1240	!$OMP PRIVATE ( ij,i,j )
1241	do ij = 1 , num_tiles
1242
1243	call toposhad_init (ht_shad,ht_loc, &
1244	shadowmask,nested,ni, &
1245	ids,ide, jds,jde, kds,kde, &
1246	ims,ime, jms,jme, kms,kme, &
1247	ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1248	i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1249	min(j_end(ij), jde-1), kts, kte )
1250
1251	enddo
1252	!$OMP END PARALLEL DO
1253
1254
1255	!$OMP PARALLEL DO &
1256	!$OMP PRIVATE ( ij,i,j )
1257	do ij = 1 , num_tiles
1258
1259	call toposhad (xlat,xlong,sina,cosa,xtime,gmt,radt,declin, &
1260	dx,dy,ht_shad,ht_loc,ni, &
1261	shadowmask,shadlen, &
1262	ids,ide, jds,jde, kds,kde, &
1263	ims,ime, jms,jme, kms,kme, &
1264	ips,min(ipe,ide-1), jps,min(jpe,jde-1), kps,kpe, &
1265	i_start(ij),min(i_end(ij), ide-1),j_start(ij),&
1266	min(j_end(ij), jde-1), kts, kte )
1267
1268	enddo
1269	!$OMP END PARALLEL DO
1270
1271	#if defined( DM_PARALLEL ) && (EM_CORE == 1)
1272	# include "HALO_TOPOSHAD.inc"
1273	#endif
1274	enddo
1275	endif
1276
1277	END SUBROUTINE pre_radiation_driver
1278
1279	!---------------------------------------------------------------------
1280	!BOP
1281	! !IROUTINE: radconst - compute radiation terms
1282	! !INTERFAC:
1283	SUBROUTINE radconst(XTIME,DECLIN,SOLCON,JULIAN, &
1284	DEGRAD,DPD )
1285	!---------------------------------------------------------------------
1286	USE module_wrf_error
1287	IMPLICIT NONE
1288	!---------------------------------------------------------------------
1289
1290	! !ARGUMENTS:
1291	REAL, INTENT(IN ) :: DEGRAD,DPD,XTIME,JULIAN
1292	REAL, INTENT(OUT ) :: DECLIN,SOLCON
1293	REAL :: OBECL,SINOB,SXLONG,ARG, &
1294	DECDEG,DJUL,RJUL,ECCFAC
1295	!
1296	! !DESCRIPTION:
1297	! Compute terms used in radiation physics
1298	!EOP
1299
1300	! for short wave radiation
1301
1302	DECLIN=0.
1303	SOLCON=0.
1304
1305	!-----OBECL : OBLIQUITY = 23.5 DEGREE.
1306
1307	OBECL=23.5*DEGRAD
1308	SINOB=SIN(OBECL)
1309
1310	!-----CALCULATE LONGITUDE OF THE SUN FROM VERNAL EQUINOX:
1311
1312	IF(JULIAN.GE.80.)SXLONG=DPD*(JULIAN-80.)
1313	IF(JULIAN.LT.80.)SXLONG=DPD*(JULIAN+285.)
1314	SXLONG=SXLONG*DEGRAD
1315	ARG=SINOB*SIN(SXLONG)
1316	DECLIN=ASIN(ARG)
1317	DECDEG=DECLIN/DEGRAD
1318	!----SOLAR CONSTANT ECCENTRICITY FACTOR (PALTRIDGE AND PLATT 1976)
1319	DJUL=JULIAN*360./365.
1320	RJUL=DJUL*DEGRAD
1321	ECCFAC=1.000110+0.034221COS(RJUL)+0.001280SIN(RJUL)+0.000719* &
1322	COS(2RJUL)+0.000077SIN(2*RJUL)
1323	SOLCON=1370.*ECCFAC
1324
1325	END SUBROUTINE radconst
1326
1327	!---------------------------------------------------------------------
1328	!BOP
1329	! !IROUTINE: cal_cldfra - Compute cloud fraction
1330	! !INTERFACE:
1331	SUBROUTINE cal_cldfra(CLDFRA,QC,QI,F_QC,F_QI, &
1332	ids,ide, jds,jde, kds,kde, &
1333	ims,ime, jms,jme, kms,kme, &
1334	its,ite, jts,jte, kts,kte )
1335	!---------------------------------------------------------------------
1336	IMPLICIT NONE
1337	!---------------------------------------------------------------------
1338	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1339	ims,ime, jms,jme, kms,kme, &
1340	its,ite, jts,jte, kts,kte
1341
1342	!
1343	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1344	CLDFRA
1345
1346	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1347	QI, &
1348	QC
1349
1350	LOGICAL,INTENT(IN) :: F_QC,F_QI
1351
1352	REAL thresh
1353	INTEGER:: i,j,k
1354	! !DESCRIPTION:
1355	! Compute cloud fraction from input ice and cloud water fields
1356	! if provided.
1357	!
1358	! Whether QI or QC is active or not is determined from the indices of
1359	! the fields into the 4D scalar arrays in WRF. These indices are
1360	! P_QI and P_QC, respectively, and they are passed in to the routine
1361	! to enable testing to see if QI and QC represent active fields in
1362	! the moisture 4D scalar array carried by WRF.
1363	!
1364	! If a field is active its index will have a value greater than or
1365	! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1366	! this routine.
1367	!EOP
1368	!---------------------------------------------------------------------
1369	thresh=1.0e-6
1370
1371	IF ( f_qi .AND. f_qc ) THEN
1372	DO j = jts,jte
1373	DO k = kts,kte
1374	DO i = its,ite
1375	IF ( QC(i,k,j)+QI(I,k,j) .gt. thresh) THEN
1376	CLDFRA(i,k,j)=1.
1377	ELSE
1378	CLDFRA(i,k,j)=0.
1379	ENDIF
1380	ENDDO
1381	ENDDO
1382	ENDDO
1383	ELSE IF ( f_qc ) THEN
1384	DO j = jts,jte
1385	DO k = kts,kte
1386	DO i = its,ite
1387	IF ( QC(i,k,j) .gt. thresh) THEN
1388	CLDFRA(i,k,j)=1.
1389	ELSE
1390	CLDFRA(i,k,j)=0.
1391	ENDIF
1392	ENDDO
1393	ENDDO
1394	ENDDO
1395	ELSE
1396	DO j = jts,jte
1397	DO k = kts,kte
1398	DO i = its,ite
1399	CLDFRA(i,k,j)=0.
1400	ENDDO
1401	ENDDO
1402	ENDDO
1403	ENDIF
1404
1405	END SUBROUTINE cal_cldfra
1406
1407	!BOP
1408	! !IROUTINE: cal_cldfra2 - Compute cloud fraction
1409	! !INTERFACE:
1410	! cal_cldfra_xr - Compute cloud fraction.
1411	! Code adapted from that in module_ra_gfdleta.F in WRF_v2.0.3 by James Done
1412	!!
1413	!!--- Cloud fraction parameterization follows Randall, 1994
1414	!! (see Hong et al., 1998)
1415	!! (modified by Ferrier, Feb '02)
1416	!
1417	SUBROUTINE cal_cldfra2(CLDFRA, QV, QC, QI, QS, &
1418	F_QV, F_QC, F_QI, F_QS, t_phy, p_phy, &
1419	F_ICE_PHY,F_RAIN_PHY, &
1420	ids,ide, jds,jde, kds,kde, &
1421	ims,ime, jms,jme, kms,kme, &
1422	its,ite, jts,jte, kts,kte )
1423	!---------------------------------------------------------------------
1424	IMPLICIT NONE
1425	!---------------------------------------------------------------------
1426	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1427	ims,ime, jms,jme, kms,kme, &
1428	its,ite, jts,jte, kts,kte
1429
1430	!
1431	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: &
1432	CLDFRA
1433
1434	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: &
1435	QV, &
1436	QI, &
1437	QC, &
1438	QS, &
1439	t_phy, &
1440	p_phy, &
1441	F_ICE_PHY, &
1442	F_RAIN_PHY
1443
1444	LOGICAL,INTENT(IN) :: F_QC,F_QI,F_QV,F_QS
1445
1446	! REAL thresh
1447	INTEGER:: i,j,k
1448	REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight, QIMID, QWMID, QCLD, DENOM, ARG, SUBSAT
1449
1450	REAL ,PARAMETER :: ALPHA0=100., GAMMA=0.49, QCLDMIN=1.E-12, &
1451	PEXP=0.25, RHGRID=1.0
1452	REAL , PARAMETER :: SVP1=0.61078
1453	REAL , PARAMETER :: SVP2=17.2693882
1454	REAL , PARAMETER :: SVPI2=21.8745584
1455	REAL , PARAMETER :: SVP3=35.86
1456	REAL , PARAMETER :: SVPI3=7.66
1457	REAL , PARAMETER :: SVPT0=273.15
1458	REAL , PARAMETER :: r_d = 287.
1459	REAL , PARAMETER :: r_v = 461.6
1460	REAL , PARAMETER :: ep_2=r_d/r_v
1461	! !DESCRIPTION:
1462	! Compute cloud fraction from input ice and cloud water fields
1463	! if provided.
1464	!
1465	! Whether QI or QC is active or not is determined from the indices of
1466	! the fields into the 4D scalar arrays in WRF. These indices are
1467	! P_QI and P_QC, respectively, and they are passed in to the routine
1468	! to enable testing to see if QI and QC represent active fields in
1469	! the moisture 4D scalar array carried by WRF.
1470	!
1471	! If a field is active its index will have a value greater than or
1472	! equal to PARAM_FIRST_SCALAR, which is also an input argument to
1473	! this routine.
1474	!EOP
1475
1476
1477	!-----------------------------------------------------------------------
1478	!--- COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
1479	! (modified by Ferrier, Feb '02)
1480	!
1481	!--- Cloud fraction parameterization follows Randall, 1994
1482	! (see Hong et al., 1998)
1483	!-----------------------------------------------------------------------
1484	! Note: ep_2=287./461.6 Rd/Rv
1485	! Note: R_D=287.
1486
1487	! Alternative calculation for critical RH for grid saturation
1488	! RHGRID=0.90+.08((100.-DX)/95.)*.5
1489
1490	! Calculate saturation mixing ratio weighted according to the fractions of
1491	! water and ice.
1492	! Following:
1493	! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
1494	! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
1495	!
1496	! over ice over water
1497	! a = 21.8745584 17.2693882
1498	! b = 7.66 35.86
1499
1500	!---------------------------------------------------------------------
1501
1502	DO j = jts,jte
1503	DO k = kts,kte
1504	DO i = its,ite
1505	tc = t_phy(i,k,j) - SVPT0
1506	esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(i,k,j) - SVP3 ) )
1507	esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(i,k,j) - SVPI3 ) )
1508	QVSW = EP_2 * esw / ( p_phy(i,k,j) - esw )
1509	QVSI = EP_2 * esi / ( p_phy(i,k,j) - esi )
1510
1511	IF ( F_QI .and. F_QC .and. F_QS) THEN
1512	QCLD=QI(i,k,j)+QC(i,k,j)+QS(I,k,j)
1513	IF (QCLD .LT. QCLDMIN) THEN
1514	weight = 0.
1515	ELSE
1516	weight = (QI(i,k,j)+QS(I,k,j)) / QCLD
1517	ENDIF
1518	ELSE IF ( F_QC ) THEN
1519
1520	! Mixing ratios of cloud water & total ice (cloud ice + snow).
1521	! Mixing ratios of rain are not considered in this scheme.
1522	! F_ICE is fraction of ice
1523	! F_RAIN is fraction of rain
1524
1525	QIMID=QC(i,k,j)*F_ICE_PHY(i,k,j)
1526	QWMID=(QC(i,k,j)-QIMID)*(1.-F_RAIN_PHY(i,k,j))
1527
1528
1529	!
1530	!--- Total "cloud" mixing ratio, QCLD. Rain is not part of cloud,
1531	! only cloud water + cloud ice + snow
1532	!
1533	QCLD=QWMID+QIMID
1534	IF (QCLD .LT. QCLDMIN) THEN
1535	weight = 0.
1536	ELSE
1537	weight = F_ICE_PHY(i,k,j)
1538	ENDIF
1539
1540	ELSE
1541	CLDFRA(i,k,j)=0.
1542	ENDIF ! IF ( F_QI .and. F_QC )
1543
1544
1545	QVS_WEIGHT = (1-weight)QVSW + weightQVSI
1546	RHUM=QV(i,k,j)/QVS_WEIGHT !--- Relative humidity
1547	!
1548	!--- Determine cloud fraction (modified from original algorithm)
1549	!
1550	IF (QCLD .LT. QCLDMIN) THEN
1551	!
1552	!--- Assume zero cloud fraction if there is no cloud mixing ratio
1553	!
1554	CLDFRA(i,k,j)=0.
1555	ELSEIF(RHUM.GE.RHGRID)THEN
1556	!
1557	!--- Assume cloud fraction of unity if near saturation and the cloud
1558	! mixing ratio is at or above the minimum threshold
1559	!
1560	CLDFRA(i,k,j)=1.
1561	ELSE
1562	!
1563	!--- Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
1564	! modified based on assumed grid-scale saturation at RH=RHgrid.
1565	!
1566	SUBSAT=MAX(1.E-10,RHGRID*QVS_WEIGHT-QV(i,k,j))
1567	DENOM=(SUBSAT)**GAMMA
1568	ARG=MAX(-6.9, -ALPHA0*QCLD/DENOM) ! <-- EXP(-6.9)=.001
1569	! prevent negative values (new)
1570	RHUM=MAX(1.E-10, RHUM)
1571	CLDFRA(i,k,j)=(RHUM/RHGRID)*PEXP(1.-EXP(ARG))
1572	!! ARG=-1000*QCLD/(RHUM-RHGRID)
1573	!! ARG=MAX(ARG, ARGMIN)
1574	!! CLDFRA(i,k,j)=(RHUM/RHGRID)*(1.-EXP(ARG))
1575	IF (CLDFRA(i,k,j) .LT. .01) CLDFRA(i,k,j)=0.
1576	ENDIF !--- End IF (QCLD .LT. QCLDMIN) ...
1577	ENDDO !--- End DO i
1578	ENDDO !--- End DO k
1579	ENDDO !--- End DO j
1580
1581	END SUBROUTINE cal_cldfra2
1582
1583
1584	SUBROUTINE toposhad_init(ht_shad,ht_loc,shadowmask,nested,iter, &
1585	ids,ide, jds,jde, kds,kde, &
1586	ims,ime, jms,jme, kms,kme, &
1587	ips,ipe, jps,jpe, kps,kpe, &
1588	its,ite, jts,jte, kts,kte )
1589
1590	USE module_model_constants
1591
1592	implicit none
1593
1594	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
1595	ims,ime, jms,jme, kms,kme, &
1596	ips,ipe, jps,jpe, kps,kpe, &
1597	its,ite, jts,jte, kts,kte
1598
1599	LOGICAL, INTENT(IN) :: nested
1600
1601	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad, ht_loc
1602
1603	INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
1604	INTEGER, INTENT(IN) :: iter
1605
1606	! Local variables
1607
1608	INTEGER :: i, j
1609
1610	if (iter.eq.1) then
1611
1612	! Initialize shadow mask
1613	do j=jts,jte
1614	do i=its,ite
1615	shadowmask(i,j) = 0
1616	ENDDO
1617	ENDDO
1618
1619	! Initialize shading height
1620
1621	IF ( nested ) THEN ! Do not overwrite input from parent domain
1622	do j=max(jts,jds+2),min(jte,jde-3)
1623	do i=max(its,ids+2),min(ite,ide-3)
1624	ht_shad(i,j) = ht_loc(i,j)-0.001
1625	ENDDO
1626	ENDDO
1627	ELSE
1628	do j=jts,jte
1629	do i=its,ide
1630	ht_shad(i,j) = ht_loc(i,j)-0.001
1631	ENDDO
1632	ENDDO
1633	ENDIF
1634
1635	IF ( nested ) THEN ! Check if a shadow exceeding the topography height is available at the lateral domain edge from nesting
1636	if (its.eq.ids) then
1637	do j=jts,jte
1638	if (ht_shad(its,j) .gt. ht_loc(its,j)) then
1639	shadowmask(its,j) = 1
1640	ht_loc(its,j) = ht_shad(its,j)
1641	endif
1642	if (ht_shad(its+1,j) .gt. ht_loc(its+1,j)) then
1643	shadowmask(its+1,j) = 1
1644	ht_loc(its+1,j) = ht_shad(its+1,j)
1645	endif
1646	enddo
1647	endif
1648	if (ite.eq.ide-1) then
1649	do j=jts,jte
1650	if (ht_shad(ite,j) .gt. ht_loc(ite,j)) then
1651	shadowmask(ite,j) = 1
1652	ht_loc(ite,j) = ht_shad(ite,j)
1653	endif
1654	if (ht_shad(ite-1,j) .gt. ht_loc(ite-1,j)) then
1655	shadowmask(ite-1,j) = 1
1656	ht_loc(ite-1,j) = ht_shad(ite-1,j)
1657	endif
1658	enddo
1659	endif
1660	if (jts.eq.jds) then
1661	do i=its,ite
1662	if (ht_shad(i,jts) .gt. ht_loc(i,jts)) then
1663	shadowmask(i,jts) = 1
1664	ht_loc(i,jts) = ht_shad(i,jts)
1665	endif
1666	if (ht_shad(i,jts+1) .gt. ht_loc(i,jts+1)) then
1667	shadowmask(i,jts+1) = 1
1668	ht_loc(i,jts+1) = ht_shad(i,jts+1)
1669	endif
1670	enddo
1671	endif
1672	if (jte.eq.jde-1) then
1673	do i=its,ite
1674	if (ht_shad(i,jte) .gt. ht_loc(i,jte)) then
1675	shadowmask(i,jte) = 1
1676	ht_loc(i,jte) = ht_shad(i,jte)
1677	endif
1678	if (ht_shad(i,jte-1) .gt. ht_loc(i,jte-1)) then
1679	shadowmask(i,jte-1) = 1
1680	ht_loc(i,jte-1) = ht_shad(i,jte-1)
1681	endif
1682	enddo
1683	endif
1684	ENDIF
1685
1686	else
1687
1688	! Fill the local topography field at the points next to internal tile boundaries with ht_shad values
1689	! A 2-pt halo has been applied to the ht_shad before the repeated call of this subroutine
1690
1691	if ((its.ne.ids).and.(its.eq.ips)) then
1692	do j=jts-2,jte+2
1693	ht_loc(its-1,j) = max(ht_loc(its-1,j),ht_shad(its-1,j))
1694	ht_loc(its-2,j) = max(ht_loc(its-2,j),ht_shad(its-2,j))
1695	enddo
1696	endif
1697	if ((ite.ne.ide-1).and.(ite.eq.ipe)) then
1698	do j=jts-2,jte+2
1699	ht_loc(ite+1,j) = max(ht_loc(ite+1,j),ht_shad(ite+1,j))
1700	ht_loc(ite+2,j) = max(ht_loc(ite+2,j),ht_shad(ite+2,j))
1701	enddo
1702	endif
1703	if ((jts.ne.jds).and.(jts.eq.jps)) then
1704	do i=its-2,ite+2
1705	ht_loc(i,jts-1) = max(ht_loc(i,jts-1),ht_shad(i,jts-1))
1706	ht_loc(i,jts-2) = max(ht_loc(i,jts-2),ht_shad(i,jts-2))
1707	enddo
1708	endif
1709	if ((jte.ne.jde-1).and.(jte.eq.jpe)) then
1710	do i=its-2,ite+2
1711	ht_loc(i,jte+1) = max(ht_loc(i,jte+1),ht_shad(i,jte+1))
1712	ht_loc(i,jte+2) = max(ht_loc(i,jte+2),ht_shad(i,jte+2))
1713	enddo
1714	endif
1715
1716	endif
1717
1718	END SUBROUTINE toposhad_init
1719
1720
1721
1722
1723	SUBROUTINE toposhad(xlat,xlong,sina,cosa,xtime,gmt,radfrq,declin, &
1724	dx,dy,ht_shad,ht_loc,iter, &
1725	shadowmask,shadlen, &
1726	ids,ide, jds,jde, kds,kde, &
1727	ims,ime, jms,jme, kms,kme, &
1728	ips,ipe, jps,jpe, kps,kpe, &
1729	its,ite, jts,jte, kts,kte )
1730
1731
1732	USE module_model_constants
1733
1734	implicit none
1735
1736	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
1737	ims,ime, jms,jme, kms,kme, &
1738	ips,ipe, jps,jpe, kps,kpe, &
1739	its,ite, jts,jte, kts,kte
1740
1741	INTEGER, INTENT(IN) :: iter
1742
1743	REAL, INTENT(IN) :: RADFRQ,XTIME,DECLIN,dx,dy,gmt,shadlen
1744
1745	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT, XLONG, sina, cosa
1746
1747	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ht_shad,ht_loc
1748
1749	INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: shadowmask
1750
1751	! Local variables
1752
1753	REAL :: pi, xt24, wgt, ri, rj, argu, sol_azi, topoelev, dxabs, tloctm, hrang, xxlat, csza
1754	INTEGER :: gpshad, ii, jj, i1, i2, j1, j2, i, j
1755
1756
1757
1758	XT24=MOD(XTIME+RADFRQ*0.5,1440.)
1759	pi = 4.*atan(1.)
1760	gpshad = int(shadlen/dx+1.)
1761
1762	if (iter.eq.1) then
1763
1764
1765	j_loop1: DO J=jts,jte
1766	i_loop1: DO I=its,ite
1767
1768	TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
1769	HRANG=15.(TLOCTM-12.)DEGRAD
1770	XXLAT=XLAT(i,j)*DEGRAD
1771	CSZA=SIN(XXLAT)SIN(DECLIN)+COS(XXLAT)COS(DECLIN)*COS(HRANG)
1772
1773	if (csza.lt.1.e-2) then ! shadow mask does not need to be computed
1774	shadowmask(i,j) = 0
1775	ht_shad(i,j) = ht_loc(i,j)-0.001
1776	goto 120
1777	endif
1778
1779	! Solar azimuth angle
1780
1781	argu=(cszasin(XXLAT)-sin(DECLIN))/(sin(acos(csza))cos(XXLAT))
1782	if (argu.gt.1) argu = 1
1783	if (argu.lt.-1) argu = -1
1784	sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
1785	if (cosa(i,j).ge.0) then
1786	sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
1787	else
1788	sol_azi = sol_azi + pi - asin(sina(i,j))
1789	endif
1790
1791	! Scan for higher surrounding topography
1792
1793	if ((sol_azi.gt.1.75pi).or.(sol_azi.lt.0.25pi)) then ! sun is in the northern quarter
1794
1795	do jj = j+1,j+gpshad
1796	ri = i + (jj-j)*tan(sol_azi)
1797	i1 = int(ri)
1798	i2 = i1+1
1799	wgt = ri-i1
1800	dxabs = sqrt((dy(jj-j))2+(dx(ri-i))**2)
1801	if ((jj.ge.jpe+1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
1802	if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
1803	goto 120
1804	endif
1805	topoelev=atan((wgtht_loc(i2,jj)+(1.-wgt)ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
1806	if (sin(topoelev).ge.csza) then
1807	shadowmask(i,j) = 1
1808	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1809	endif
1810	enddo
1811
1812	else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
1813	do ii = i+1,i+gpshad
1814	rj = j - (ii-i)*tan(pi/2.+sol_azi)
1815	j1 = int(rj)
1816	j2 = j1+1
1817	wgt = rj-j1
1818	dxabs = sqrt((dx(ii-i))2+(dy(rj-j))**2)
1819	if ((ii.ge.ipe+1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
1820	if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
1821	goto 120
1822	endif
1823	topoelev=atan((wgtht_loc(ii,j2)+(1.-wgt)ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
1824	if (sin(topoelev).ge.csza) then
1825	shadowmask(i,j) = 1
1826	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1827	endif
1828	enddo
1829
1830	else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
1831	do jj = j-1,j-gpshad,-1
1832	ri = i + (jj-j)*tan(sol_azi)
1833	i1 = int(ri)
1834	i2 = i1+1
1835	wgt = ri-i1
1836	dxabs = sqrt((dy(jj-j))2+(dx(ri-i))**2)
1837	if ((jj.le.jps-1).or.(i1.le.ips-1).or.(i2.ge.ipe+1)) then
1838	if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
1839	goto 120
1840	endif
1841	topoelev=atan((wgtht_loc(i2,jj)+(1.-wgt)ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
1842	if (sin(topoelev).ge.csza) then
1843	shadowmask(i,j) = 1
1844	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1845	endif
1846	enddo
1847
1848	else ! sun is in the western quarter
1849	do ii = i-1,i-gpshad,-1
1850	rj = j - (ii-i)*tan(pi/2.+sol_azi)
1851	j1 = int(rj)
1852	j2 = j1+1
1853	wgt = rj-j1
1854	dxabs = sqrt((dx(ii-i))2+(dy(rj-j))**2)
1855	if ((ii.le.ips-1).or.(j1.le.jps-1).or.(j2.ge.jpe+1)) then
1856	if (shadowmask(i,j).eq.0) shadowmask(i,j) = -1
1857	goto 120
1858	endif
1859	topoelev=atan((wgtht_loc(ii,j2)+(1.-wgt)ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
1860	if (sin(topoelev).ge.csza) then
1861	shadowmask(i,j) = 1
1862	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1863	endif
1864	enddo
1865	endif
1866
1867	120 continue
1868
1869	ENDDO i_loop1
1870	ENDDO j_loop1
1871
1872	else ! iteration > 1
1873
1874
1875	j_loop2: DO J=jts,jte
1876	i_loop2: DO I=its,ite
1877
1878	! if (shadowmask(i,j).eq.-1) then ! this indicates that the search ended at a lateral boundary during iteration 1
1879
1880	TLOCTM=GMT+XT24/60.+XLONG(i,j)/15.
1881	HRANG=15.(TLOCTM-12.)DEGRAD
1882	XXLAT=XLAT(i,j)*DEGRAD
1883	CSZA=SIN(XXLAT)SIN(DECLIN)+COS(XXLAT)COS(DECLIN)*COS(HRANG)
1884
1885	! Solar azimuth angle
1886
1887	argu=(cszasin(XXLAT)-sin(DECLIN))/(sin(acos(csza))cos(XXLAT))
1888	if (argu.gt.1) argu = 1
1889	if (argu.lt.-1) argu = -1
1890	sol_azi = sign(acos(argu),sin(HRANG))+pi ! azimuth angle of the sun
1891	if (cosa(i,j).ge.0) then
1892	sol_azi = sol_azi + asin(sina(i,j)) ! rotation towards WRF grid
1893	else
1894	sol_azi = sol_azi + pi - asin(sina(i,j))
1895	endif
1896
1897	! Scan for higher surrounding topography
1898
1899	if ((sol_azi.gt.1.75pi).or.(sol_azi.lt.0.25pi)) then ! sun is in the northern quarter
1900
1901	do jj = j+1,j+gpshad
1902	ri = i + (jj-j)*tan(sol_azi)
1903	i1 = int(ri)
1904	i2 = i1+1
1905	wgt = ri-i1
1906	dxabs = sqrt((dy(jj-j))2+(dx(ri-i))**2)
1907	if ((jj.ge.min(jde,jpe+3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220
1908	topoelev=atan((wgtht_loc(i2,jj)+(1.-wgt)ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
1909	if (sin(topoelev).ge.csza) then
1910	shadowmask(i,j) = 1
1911	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1912	endif
1913	enddo
1914
1915	else if (sol_azi.lt.0.75*pi) then ! sun is in the eastern quarter
1916	do ii = i+1,i+gpshad
1917	rj = j - (ii-i)*tan(pi/2.+sol_azi)
1918	j1 = int(rj)
1919	j2 = j1+1
1920	wgt = rj-j1
1921	dxabs = sqrt((dx(ii-i))2+(dy(rj-j))**2)
1922	if ((ii.ge.min(ide,ipe+3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220
1923	topoelev=atan((wgtht_loc(ii,j2)+(1.-wgt)ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
1924	if (sin(topoelev).ge.csza) then
1925	shadowmask(i,j) = 1
1926	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1927	endif
1928	enddo
1929
1930	else if (sol_azi.lt.1.25*pi) then ! sun is in the southern quarter
1931	do jj = j-1,j-gpshad,-1
1932	ri = i + (jj-j)*tan(sol_azi)
1933	i1 = int(ri)
1934	i2 = i1+1
1935	wgt = ri-i1
1936	dxabs = sqrt((dy(jj-j))2+(dx(ri-i))**2)
1937	if ((jj.le.max(jds-1,jps-3)).or.(i1.le.max(ids-1,ips-3)).or.(i2.ge.min(ide,ipe+3))) goto 220
1938	topoelev=atan((wgtht_loc(i2,jj)+(1.-wgt)ht_loc(i1,jj)-ht_loc(i,j))/dxabs)
1939	if (sin(topoelev).ge.csza) then
1940	shadowmask(i,j) = 1
1941	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1942	endif
1943	enddo
1944
1945	else ! sun is in the western quarter
1946	do ii = i-1,i-gpshad,-1
1947	rj = j - (ii-i)*tan(pi/2.+sol_azi)
1948	j1 = int(rj)
1949	j2 = j1+1
1950	wgt = rj-j1
1951	dxabs = sqrt((dx(ii-i))2+(dy(rj-j))**2)
1952	if ((ii.le.max(ids-1,ips-3)).or.(j1.le.max(jds-1,jps-3)).or.(j2.ge.min(jde,jpe+3))) goto 220
1953	topoelev=atan((wgtht_loc(ii,j2)+(1.-wgt)ht_loc(ii,j1)-ht_loc(i,j))/dxabs)
1954	if (sin(topoelev).ge.csza) then
1955	shadowmask(i,j) = 1
1956	ht_shad(i,j) = max(ht_shad(i,j),ht_loc(i,j)+dxabs*(tan(topoelev)-tan(asin(csza))))
1957	endif
1958	enddo
1959	endif
1960
1961	220 continue
1962	! endif
1963
1964	ENDDO i_loop2
1965	ENDDO j_loop2
1966
1967	endif ! iteration
1968
1969	END SUBROUTINE toposhad
1970
1971	END MODULE module_radiation_driver

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

Download in other formats:

Original Format