Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_radiation_driver.F @ 1379

Last change on this file since 1379 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format