Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
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

Rev	Line
[2759]	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