Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_camuwpbl_driver.F @ 1564

Last change on this file since 1564 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:

Property svn:executable set to *

File size: 62.8 KB

Line
1	MODULE module_bl_camuwpbl_driver
2	!Note: Comments starting with "!!" are directly taken from CAM's interface routine for the UW PBL
3
4	!!----------------------------------------------------------------------------------------------------- !
5	!! Module to compute vertical diffusion of momentum, moisture, trace constituents !
6	!! and static energy. Separate modules compute !
7	!! 1. stresses associated with turbulent flow over orography !
8	!! ( turbulent mountain stress ) !
9	!! 2. eddy diffusivities, including nonlocal tranport terms !
10	!! 3. molecular diffusivities !
11	!! 4. coming soon... gravity wave drag !
12	!! Lastly, a implicit diffusion solver is called, and tendencies retrieved by !
13	!! differencing the diffused and initial states. !
14	!! !
15	!! Calling sequence: !
16	!! !
17	!! vertical_diffusion_init Initializes vertical diffustion constants and modules !
18	!! init_molec_diff Initializes molecular diffusivity module !
19	!! init_eddy_diff Initializes eddy diffusivity module (includes PBL) !
20	!! init_tms Initializes turbulent mountain stress module !
21	!! init_vdiff Initializes diffusion solver module !
22	!! vertical_diffusion_ts_init Time step initialization (only used for upper boundary condition) !
23	!! vertical_diffusion_tend Computes vertical diffusion tendencies !
24	!! compute_tms Computes turbulent mountain stresses !
25	!! compute_eddy_diff Computes eddy diffusivities and countergradient terms !
26	!! compute_vdiff Solves vertical diffusion equations, including molecular diffusivities !
27	!! !
28	!!---------------------------Code history-------------------------------------------------------------- !
29	!! J. Rosinski : Jun. 1992 !
30	!! J. McCaa : Sep. 2004 !
31	!! S. Park : Aug. 2006, Dec. 2008. Jan. 2010 !
32	! B. Singh : Nov. 2010 (ported to WRF by balwinder.singh@pnl.gov)
33	!!----------------------------------------------------------------------------------------------------- !
34
35	use shr_kind_mod, only : r8 => shr_kind_r8
36	use module_cam_support, only : pcols, pver, pverp, endrun, iulog,fieldname_len,pcnst
37	use constituents, only : qmin
38	use diffusion_solver, only : vdiff_selector
39	use physconst, only : &
40	cpair , & ! Specific heat of dry air
41	gravit , & ! Acceleration due to gravity
42	rair , & ! Gas constant for dry air
43	zvir , & ! rh2o/rair - 1
44	latvap , & ! Latent heat of vaporization
45	latice , & ! Latent heat of fusion
46	karman , & ! von Karman constant
47	mwdry , & ! Molecular weight of dry air
48	avogad , & ! Avogadro's number
49	boltz ! Boltzman's constant
50
51	implicit none
52	private
53	save
54
55	!! ----------------- !
56	!! Public interfaces !
57	!! ----------------- !
58
59	public camuwpblinit ! Initialization
60	public camuwpbl ! Driver for the PBL scheme
61	public vd_register ! Init routine for constituents
62
63	!! ------------ !
64	!! Private data !
65	!! ------------ !
66
67	character(len=16) :: eddy_scheme !! Default set in phys_control.F90, use namelist to change
68	!! 'HB' = Holtslag and Boville (default)
69	!! 'HBR' = Holtslag and Boville and Rash
70	!! 'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme )
71	integer, parameter :: nturb = 5 !! Number of iterations for solution ( when 'diag_TKE' scheme is selected )
72	logical, parameter :: wstarent = .true. !! Use wstar (.true.) or TKE (.false.) entrainment closure ( when 'diag_TKE' scheme is selected )
73	logical :: do_pseudocon_diff = .false. !! If .true., do pseudo-conservative variables diffusion
74
75
76	character(len=16) :: shallow_scheme !! For checking compatibility between eddy diffusion and shallow convection schemes
77	!! 'Hack' = Hack Shallow Convection Scheme
78	!! 'UW' = Park and Bretherton ( UW Shallow Convection Scheme )
79	character(len=16) :: microp_scheme !! Microphysics scheme
80
81	logical :: do_molec_diff = .false. !! Switch for molecular diffusion
82	logical :: do_tms !! Switch for turbulent mountain stress
83	real(r8) :: tms_orocnst !! Converts from standard deviation to height
84
85	type(vdiff_selector) :: fieldlist_wet !! Logical switches for moist mixing ratio diffusion
86	type(vdiff_selector) :: fieldlist_dry !! Logical switches for dry mixing ratio diffusion
87	integer :: ntop !! Top interface level to which vertical diffusion is applied ( = 1 ).
88	integer :: nbot !! Bottom interface level to which vertical diffusion is applied ( = pver ).
89	integer :: tke_idx, kvh_idx, kvm_idx !! TKE and eddy diffusivity indices for fields in the physics buffer
90	integer :: turbtype_idx, smaw_idx !! Turbulence type and instability functions
91	integer :: tauresx_idx, tauresy_idx !! Redisual stress for implicit surface stress
92
93	integer :: ixcldice, ixcldliq !! Constituent indices for cloud liquid and ice water
94	integer :: ixnumice, ixnumliq
95	integer :: wgustd_index
96	logical :: vd_registered = .false. !! Detect if vd_register called
97
98
99	CONTAINS
100
101	subroutine camuwpbl(dt,u_phy,v_phy,th_phy,rho,qv_curr,hfx &
102	,qfx,ustar,rublten,rvblten,rthblten,rqvblten,rqcblten &
103	,tke_pbl,pblh2d,kpbl2d,p8w,p_phy,z,t_phy,qc_curr,qi_curr,z_at_w &
104	,cldfra,ht,rthratenlw,exner,itimestep,tauresx2d,tauresy2d,kvh3d &
105	,kvm3d,tpert2d,qpert2d,wpert2d &
106	,ids,ide, jds,jde, kds,kde &
107	,ims,ime, jms,jme, kms,kme &
108	,its,ite, jts,jte, kts,kte)
109
110
111	!!---------------------------------------------------- !
112	!! This is an interface routine for vertical diffusion !
113	!!---------------------------------------------------- !
114	use module_cam_support, only : pcols
115	use trb_mtn_stress, only : compute_tms
116	use eddy_diff, only : compute_eddy_diff
117	use wv_saturation, only : fqsatd, aqsat
118	use molec_diff, only : compute_molec_diff, vd_lu_qdecomp
119	use constituents, only : qmincg, qmin
120	use diffusion_solver !!, only : compute_vdiff, any, operator(.not.)
121
122	implicit none
123
124	!------------------------------------------------------------------------!
125	! Input !
126	!------------------------------------------------------------------------!
127	integer, intent(in) :: ids,ide, jds,jde, kds,kde
128	integer, intent(in) :: ims,ime, jms,jme, kms,kme
129	integer, intent(in) :: its,ite, jts,jte, kts,kte
130	integer, intent(in) :: itimestep
131
132	real, intent(in) :: dt ! Time step (s)
133
134	real, dimension( ims:ime,jms:jme ), intent(in) :: ustar ! Friction velocity (m/s)
135	real, dimension( ims:ime,jms:jme ), intent(in) :: hfx ! Sensible heat flux at surface (w/m2)
136	real, dimension( ims:ime,jms:jme ), intent(in) :: qfx ! Moisture flux at surface (kg m-2 s-1)
137	real, dimension( ims:ime,jms:jme ), intent(in) :: ht ! Terrain height (m)
138
139
140	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rho ! Air density (kg/m3)
141	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: th_phy ! Potential temperature (K)
142	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: u_phy ! X-component of wind (m/s)
143	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: v_phy ! Y-component of wind (m/s)
144	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p_phy ! Pressure at mid-level (Pa)
145	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p8w ! Pressure at level interface (Pa)
146	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z ! Height above sea level at mid-level (m)
147	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: t_phy ! temperature (K)
148	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qv_curr ! Water vapor mixing ratio - Moisture (kg/kg)
149	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qc_curr ! Cloud water mixing ratio - Cloud liq (kg/kg)
150	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qi_curr ! Ice mixing ratio -Cloud ice (kg/kg)
151	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_at_w ! Height above sea level at layer interfaces (m)
152	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra ! Cloud fraction [unitless]
153	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: exner ! Dimensionless pressure [unitless]
154	real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rthratenlw ! Tendency for LW ( K/s)
155
156
157	!------------------------------------------------------------------------!
158	! Output !
159	!------------------------------------------------------------------------!
160	integer, dimension( ims:ime,jms:jme ), intent(out) :: kpbl2d ! Layer index containing PBL top within or at the base interface
161	real, dimension( ims:ime,jms:jme ), intent(out) :: pblh2d ! Planetary boundary layer height (m)
162
163	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: rublten !Tendency for u_phy (Pa m s-2)
164	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: rvblten !Tendency for v_phy (Pa m s-2)
165	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: rthblten !Tendency for th_phy (Pa K s-1)
166	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: rqvblten !Tendency for qv_curr (Pa kg kg-1 s-1)
167	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: rqcblten !Tendency for qc_curr (Pa kg kg-1 s-1)
168	real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: tke_pbl !Turbulence kinetic energy from PBL scheme (m^2/s^2)
169
170
171	!---------------------------------------------------------------------------!
172	! Local, carried on from one timestep to the other (defined in registry)!
173	!---------------------------------------------------------------------------!
174
175	real, dimension( ims:ime, jms:jme ) , intent(inout ):: tauresx2d,tauresy2d !X AND Y-COMP OF RESIDUAL STRESSES(m^2/s^2)
176	real, dimension( ims:ime, jms:jme ) , intent(out) :: tpert2d ! Convective temperature excess (K)
177	real, dimension( ims:ime, jms:jme ) , intent(out) :: qpert2d ! Convective humidity excess (kg/kg)
178	real, dimension( ims:ime, jms:jme ) , intent(out) :: wpert2d ! Turbulent velocity excess (m/s)
179
180	real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: kvm3d,kvh3d !Eddy diffusivity for momentum and heat(m^2/s)
181
182	!---------------------------------------------------------------------------!
183	! Local !
184	!---------------------------------------------------------------------------!
185
186	character(128) :: errstring ! Error status for compute_vdiff
187	logical :: kvinit ! Tell compute_eddy_diff/ caleddy to initialize kvh, kvm (uses kvf)
188	logical :: is_first_step ! Flag to know if this a first time step
189	integer :: i,j,k,itsp1,itile_len,ktep1,kflip,ncol,ips
190	integer :: lchnk ! Chunk identifier
191	integer :: pcnstmax ! Max number of constituents considered
192	real(r8) :: tauFac, uMean, dp
193	real(r8) :: ztodt ! 2*delta-t (s)
194	real(r8) :: rztodt ! 1./ztodt (1/s)
195
196	real(r8) :: topflx( pcols) ! Molecular heat flux at top interface
197	real(r8) :: wpert( pcols) ! Turbulent velocity excess (m/s)
198	real(r8) :: tauresx( pcols) ! [Residual stress to be added in vdiff to correct...
199	real(r8) :: tauresy( pcols) ! for turb stress mismatch between sfc and atm accumulated.] (N/m2)
200	real(r8) :: ipbl( pcols) ! If 1, PBL is CL, while if 0, PBL is STL.
201	real(r8) :: kpblh( pcols) ! Layer index containing PBL top within or at the base interface
202	real(r8) :: wstarPBL(pcols) ! Convective velocity within PBL (m/s)
203	real(r8) :: sgh( pcols) ! Standard deviation of orography (m)
204	real(r8) :: landfrac(pcols) ! Land fraction [unitless]
205	real(r8) :: taux( pcols) ! x surface stress (N/m2)
206	real(r8) :: tauy( pcols) ! y surface stress (N/m2)
207	real(r8) :: tautotx( pcols) ! U component of total surface stress (N/m2)
208	real(r8) :: tautoty( pcols) ! V component of total surface stress (N/m2)
209	real(r8) :: ksrftms( pcols) ! Turbulent mountain stress surface drag coefficient (kg/s/m2)
210	real(r8) :: tautmsx( pcols) ! U component of turbulent mountain stress (N/m2)
211	real(r8) :: tautmsy( pcols) ! V component of turbulent mountain stress (N/m2)
212	real(r8) :: ustar8( pcols) ! Surface friction velocity (m/s)
213	real(r8) :: pblh( pcols) ! Planetary boundary layer height (m)
214	real(r8) :: tpert( pcols) ! Convective temperature excess (K)
215	real(r8) :: qpert( pcols) ! Convective humidity excess (kg/kg)
216	real(r8) :: shflx( pcols) ! Surface sensible heat flux (w/m2)
217	real(r8) :: phis( pcols) ! Geopotential at terrain height (m2/s2)
218
219	real(r8) :: cldn8( pcols,kte) ! New stratus fraction (fraction)
220	real(r8) :: qrl8( pcols,kte) ! LW radiative cooling rate(W/kg*Pa)
221	real(r8) :: wsedl8( pcols,kte) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
222	real(r8) :: dtk( pcols,kte) ! T tendency from KE dissipation
223	real(r8) :: qt( pcols,kte) !
224	real(r8) :: sl_prePBL( pcols,kte) !
225	real(r8) :: qt_prePBL( pcols,kte) !
226	real(r8) :: slten( pcols,kte) !
227	real(r8) :: qtten( pcols,kte) !
228	real(r8) :: sl( pcols,kte) !
229	real(r8) :: ftem( pcols,kte) ! Saturation vapor pressure before PBL
230	real(r8) :: ftem_prePBL(pcols,kte) ! Saturation vapor pressure before PBL
231	real(r8) :: ftem_aftPBL(pcols,kte) ! Saturation vapor pressure after PBL
232	real(r8) :: tem2( pcols,kte) ! Saturation specific humidity and RH
233	real(r8) :: t_aftPBL( pcols,kte) ! Temperature after PBL diffusion
234	real(r8) :: tten( pcols,kte) ! Temperature tendency by PBL diffusion
235	real(r8) :: rhten( pcols,kte) ! RH tendency by PBL diffusion
236	real(r8) :: qv_aft_PBL( pcols,kte) ! qv after PBL diffusion
237	real(r8) :: ql_aft_PBL( pcols,kte) ! ql after PBL diffusion
238	real(r8) :: qi_aft_PBL( pcols,kte) ! qi after PBL diffusion
239	real(r8) :: s_aft_PBL( pcols,kte) ! s after PBL diffusion
240	real(r8) :: u_aft_PBL( pcols,kte) ! u after PBL diffusion
241	real(r8) :: v_aft_PBL( pcols,kte) ! v after PBL diffusion
242	real(r8) :: qv_pro( pcols,kte) !
243	real(r8) :: ql_pro( pcols,kte) !
244	real(r8) :: qi_pro( pcols,kte) !
245	real(r8) :: s_pro( pcols,kte) !
246	real(r8) :: t_pro( pcols,kte) !
247	real(r8) :: u8( pcols,kte) ! x component of velocity in CAM's data structure (m/s)
248	real(r8) :: v8( pcols,kte) ! y component of velocity in CAM's data structure (m/s)
249	real(r8) :: t8( pcols,kte) !
250	real(r8) :: pmid8( pcols,kte) ! Pressure at the midpoints in CAM's data structure (Pa)
251	real(r8) :: pmiddry8( pcols,kte) ! Dry Pressure at the midpoints in CAM's data structure (Pa)
252	real(r8) :: zm8( pcols,kte) ! Height at the midpoints in CAM's data structure (m)
253	real(r8) :: exner8( pcols,kte) ! exner function in CAM's data structure
254	real(r8) :: s8( pcols,kte) ! Dry static energy (m2/s2)
255	real(r8) :: rpdel8( pcols,kte) ! Inverse of pressure difference (1/Pa)
256	real(r8) :: pdel8( pcols,kte) ! Pressure difference (Pa)
257	real(r8) :: rpdeldry8( pcols,kte) ! Inverse of dry pressure difference (1/Pa)
258	REAL(r8) :: stnd( pcols,kte) ! Heating rate (dry static energy tendency, W/kg)
259
260	real(r8) :: tke8( pcols,kte+1) ! Turbulent kinetic energy [ m2/s2 ]
261	real(r8) :: turbtype( pcols,kte+1) ! Turbulent interface types [ no unit ]
262	real(r8) :: smaw( pcols,kte+1) ! Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
263	! This is 1 when neutral condition (Ri=0), 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19.
264	real(r8) :: cgs( pcols,kte+1) ! Counter-gradient star [ cg/flux ]
265	real(r8) :: cgh( pcols,kte+1) ! Counter-gradient term for heat [ J/kg/m ]
266	real(r8) :: kvh( pcols,kte+1) ! Eddy diffusivity for heat [ m2/s ]
267	real(r8) :: kvm( pcols,kte+1) ! Eddy diffusivity for momentum [ m2/s ]
268	real(r8) :: kvq( pcols,kte+1) ! Eddy diffusivity for constituents [ m2/s ]
269	real(r8) :: kvh_in( pcols,kte+1) ! kvh from previous timestep [ m2/s ]
270	real(r8) :: kvm_in( pcols,kte+1) ! kvm from previous timestep [ m2/s ]
271	real(r8) :: bprod( pcols,kte+1) ! Buoyancy production of tke [ m2/s3 ]
272	real(r8) :: sprod( pcols,kte+1) ! Shear production of tke [ m2/s3 ]
273	real(r8) :: sfi( pcols,kte+1) ! Saturation fraction at interfaces [ fraction ]
274	real(r8) :: pint8( pcols,kte+1) ! Pressure at the interfaces in CAM's data structure (Pa)
275	real(r8) :: pintdry8( pcols,kte+1) ! Dry pressure at the interfaces in CAM's data structure (Pa)
276	real(r8) :: zi8( pcols,kte+1) ! Height at the interfacesin CAM's data structure (m)
277
278	real(r8) :: cloud( pcols,kte,3) ! Holder for cloud water and ice (q in cam)
279	real(r8) :: cloudtnd( pcols,kte,3) ! Holder for cloud tendencies (ptend_loc%q in cam)
280	real(r8) :: wind_tends(pcols,kte,2) ! Wind component tendencies (m/s2)
281
282	real(r8) :: cflx(pcols,pcnst) ! Surface constituent flux [ kg/m2/s ]
283
284
285	!! ----------------------- !
286	!! Main Computation Begins !
287	!! ----------------------- !
288
289	is_first_step = .false.
290	if(itimestep == 1) then
291	is_first_step = .true.
292	endif
293	!-------------------------------------------------------------------------------------!
294	!Declare maximum number of constituents to be considered. pcnst is 7 in WRF but we are!
295	!using 1 constituent as we have cflx for only water vapours !
296	!-------------------------------------------------------------------------------------!
297
298	pcnstmax = 1
299
300	ncol = pcols
301	ztodt = 2.0_r8 * DT
302	rztodt = 1.0_r8 / ztodt
303
304
305	!Few definitions in this subroutine require that ncol==1
306	if(ncol .NE. 1) then
307	call wrf_error_fatal('Number of CAM Columns (NCOL) in CAMUWPBL scheme must be 1')
308	endif
309
310	!Initialize all variables which will be used by CAM's modules to zero
311	!This is done to avoid any uninitialized variable
312
313	errstring = ''
314	topflx( :) = 0.0_r8
315	wpert( :) = 0.0_r8
316	tauresx( :) = 0.0_r8
317	tauresy( :) = 0.0_r8
318	ipbl( :) = 0.0_r8
319	kpblh( :) = 0.0_r8
320	wstarPBL(:) = 0.0_r8
321	sgh( :) = 0.0_r8
322	landfrac(:) = 0.0_r8
323	taux( :) = 0.0_r8
324	tauy( :) = 0.0_r8
325	tautotx( :) = 0.0_r8
326	tautoty( :) = 0.0_r8
327	ksrftms( :) = 0.0_r8
328	tautmsx( :) = 0.0_r8
329	tautmsy( :) = 0.0_r8
330	ustar8( :) = 0.0_r8
331	pblh( :) = 0.0_r8
332	tpert( :) = 0.0_r8
333	qpert( :) = 0.0_r8
334	shflx( :) = 0.0_r8
335	phis( :) = 0.0_r8
336
337
338	cldn8( :,:) = 0.0_r8
339	qrl8( :,:) = 0.0_r8
340	wsedl8( :,:) = 0.0_r8
341	dtk( :,:) = 0.0_r8
342	qt( :,:) = 0.0_r8
343	sl_prePBL( :,:) = 0.0_r8
344	qt_prePBL( :,:) = 0.0_r8
345	slten( :,:) = 0.0_r8
346	qtten( :,:) = 0.0_r8
347	sl( :,:) = 0.0_r8
348	ftem( :,:) = 0.0_r8
349	ftem_prePBL(:,:) = 0.0_r8
350	ftem_aftPBL(:,:) = 0.0_r8
351	tem2( :,:) = 0.0_r8
352	t_aftPBL( :,:) = 0.0_r8
353	tten( :,:) = 0.0_r8
354	rhten( :,:) = 0.0_r8
355	qv_aft_PBL( :,:) = 0.0_r8
356	ql_aft_PBL( :,:) = 0.0_r8
357	qi_aft_PBL( :,:) = 0.0_r8
358	s_aft_PBL( :,:) = 0.0_r8
359	u_aft_PBL( :,:) = 0.0_r8
360	v_aft_PBL( :,:) = 0.0_r8
361	qv_pro( :,:) = 0.0_r8
362	ql_pro( :,:) = 0.0_r8
363	qi_pro( :,:) = 0.0_r8
364	s_pro( :,:) = 0.0_r8
365	t_pro( :,:) = 0.0_r8
366	u8( :,:) = 0.0_r8
367	v8( :,:) = 0.0_r8
368	t8( :,:) = 0.0_r8
369	pmid8( :,:) = 0.0_r8
370	pmiddry8( :,:) = 0.0_r8
371	zm8( :,:) = 0.0_r8
372	exner8( :,:) = 0.0_r8
373	s8( :,:) = 0.0_r8
374	rpdel8( :,:) = 0.0_r8
375	pdel8( :,:) = 0.0_r8
376	rpdeldry8( :,:) = 0.0_r8
377	stnd( :,:) = 0.0_r8
378
379	tke8( :,:) = 0.0_r8
380	turbtype( :,:) = 0.0_r8
381	smaw( :,:) = 0.0_r8
382
383	cgs( :,:) = 0.0_r8
384	cgh( :,:) = 0.0_r8
385	kvh( :,:) = 0.0_r8
386	kvm( :,:) = 0.0_r8
387	kvq( :,:) = 0.0_r8
388	kvh_in( :,:) = 0.0_r8
389	kvm_in( :,:) = 0.0_r8
390	bprod( :,:) = 0.0_r8
391	sprod( :,:) = 0.0_r8
392	sfi( :,:) = 0.0_r8
393	pint8( :,:) = 0.0_r8
394	pintdry8( :,:) = 0.0_r8
395	zi8( :,:) = 0.0_r8
396
397	cloud( :,:,:) = 0.0_r8
398	cloudtnd( :,:,:) = 0.0_r8
399	wind_tends(:,:,:) = 0.0_r8
400	tke_pbl( :,:,:) = 0.0_r8
401	cflx(:,:) = 0.0_r8
402
403
404	!-------------------------------------------------------------------------------------!
405	!Map CAM variables to the corresponding WRF variables !
406	!Loop over the points in the tile and treat them each as a CAM Chunk !
407	!-------------------------------------------------------------------------------------!
408	itsp1 = its - 1
409	itile_len = ite - itsp1
410	do j = jts , jte
411	do i = its , ite
412
413	lchnk = (j - jts) * itile_len + (i - itsp1) !1-D index location from a 2-D tile
414	phis(1) = ht(i,j) * gravit !Used for computing dry static energy
415
416	!Flip vertically quantities computed at the mid points
417	ktep1 = kte + 1
418	do k = kts,kte
419	kflip = ktep1 - k
420
421	u8( 1,kflip) = u_phy(i,k,j) ! X-component of velocity at the mid-points (m/s) [state%u in CAM]
422	v8( 1,kflip) = v_phy(i,k,j) ! Y-component of velocity at the mid-points (m/s) [state%v in CAM]
423
424	pmid8( 1,kflip) = p_phy(i,k,j) ! Pressure at the mid-points (Pa) [state%pmid in CAM]
425
426	!The following pmiddry mapping is wrong. Presently, it is not being used in the computations
427	pmiddry8(1,kflip) = p_phy(i,k,j) ! Dry pressure at the mid-points (Pa) [state%pmiddry in CAM]
428	dp = p8w(i,k,j) - p8w(i,k+1,j) !Change in pressure (Pa)
429	pdel8( 1,kflip) = dp
430	rpdel8( 1,kflip) = 1.0_r8/dp ! Reciprocal of pressure difference (1/Pa) [state%rpdel in CAM]
431	zm8( 1,kflip) = z(i,k,j)-ht(i,j) ! Height above the ground at the midpoints (m) [state%zm in CAM]
432	t8( 1,kflip) = t_phy(i,k,j) ! Temprature at the mid points (K) [state%t in CAM]
433
434	s8( 1,kflip) = cpair t8(1,kflip) + gravitzm8(1,kflip) + phis(1) ! Dry static energy (m2/s2) [state%s in CAM]-Formula tested in vertical_diffusion.F90 of CAM
435	qrl8( 1,kflip) = rthratenlw(i,k,j)* cpair * dp ! Long Wave heating rate (W/kg*Pa)- Formula obtained from definition of qrlw(pcols,pver) in eddy_diff.F90 in CAM
436
437	!The following is set to zero currently. This value should be passed on from the microphysics subroutine in future
438	wsedl8( 1,kflip) = 0.0_r8 ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
439
440	!Following three formulas are obtained from ported CAM's ZM cumulus scheme
441	!Values of 0 cause a crash in entropy
442	cloud( 1,kflip,1) = max( qv_curr(i,k,j)/(1.+qv_curr(i,k,j)), 1e-30 ) !Specific humidity [state%q(:,:,1) in CAM]
443	cloud( 1,kflip,2) = qc_curr(i,k,j)/(1.+qv_curr(i,k,j)) !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
444	cloud( 1,kflip,3) = qi_curr(i,k,j)/(1.+qv_curr(i,k,j)) !cloud ice [state%q(:,:,3) in CAM]
445
446	exner8(1,kflip) = exner(i,k,j) !Exner function (no units)
447	cldn8( 1,kflip) = cldfra(i,k,j) !Cloud fraction (no unit)
448
449	enddo
450
451	!Future work: Merge the above and below do-loops to avoid extra do-loop, if possible
452
453	do k = kts,kte+1
454	kflip = kte - k + 2
455
456	pint8( 1,kflip) = p8w( i,k,j) ! Pressure at interfaces [state%pint in CAM]
457
458	!The following pintdry mapping is wrong.Presently, it is not being used in the computations
459	pintdry8(1,kflip) = p8w( i,k,j) ! Dry pressure at interfaces [state%pintdry in CAM]
460	zi8( 1,kflip) = z_at_w(i,k,j) -ht(i,j) ! Height at interfaces [state%zi in CAM]
461
462	!Initialize Variables to zero, these are outputs from the "compute_eddy_diff"
463	kvq(1,kflip) = 0.0_r8 ! Eddy diffusivity for constituents (m2/s)
464	cgh(1,kflip) = 0.0_r8 ! Counter-gradient term for heat
465	cgs(1,kflip) = 0.0_r8 ! Counter-gradient star (cg/flux)
466	if( is_first_step ) then
467	kvh3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for heat (m2/s)
468	kvm3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for momentum (m2/s)
469	endif
470	kvh(1,kflip) = kvh3d(i,k,j)
471	kvm(1,kflip) = kvm3d(i,k,j)
472	end do
473
474	shflx( ncol) = hfx(i,j) ! Surface sensible heat flux (w/m2)
475
476	!SGH and LANDFRAC are inputs for the compute_tms subroutine. Presently set to zero as do_tms is always false for now
477	sgh( ncol) = 0.0_r8 ! Standard deviation of orography (m)
478	landfrac(ncol) = 0.0_r8 ! Fraction (unitless)
479
480	uMean = sqrt( u_phy(i,kts,j) * u_phy(i,kts,j) + v_phy(i,kts,j) * v_phy(i,kts,j) ) ! Mean velocity
481	tauFac = rho(i,kts,j) * ustar(i,j) *ustar(i,j)/uMean
482
483	taux(ncol) = -tauFac * u_phy(i,kts,j) ! x surface stress (N/m2) [Formulation obtained from CAM's BareGround.F90]
484	tauy(ncol) = -tauFac * v_phy(i,kts,j) ! y surface stress (N/m2)
485
486	!! Retrieve 'tauresx, tauresy' from from the last timestep
487	if( is_first_step ) then
488	tauresx2d(i,j) = 0._r8
489	tauresy2d(i,j) = 0._r8
490	endif
491	tauresx(:ncol) = tauresx2d(i,j)
492	tauresy(:ncol) = tauresy2d(i,j)
493
494	!! All variables are modified by vertical diffusion
495
496	!!------------------------------------------!
497	!! Computation of turbulent mountain stress !
498	!!------------------------------------------!
499
500	!! Consistent with the computation of 'normal' drag coefficient, we are using
501	!! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v'
502	!! within the iteration loop of the PBL scheme.
503
504	if( do_tms ) then
505	call compute_tms( pcols , pver , ncol , &
506	u8 , v8 , t8 , pmid8 , &
507	exner8 , zm8 , sgh , ksrftms , &
508	tautmsx , tautmsy , landfrac )
509	!! Here, both 'taux, tautmsx' are explicit surface stresses.
510	!! Note that this 'tautotx, tautoty' are different from the total stress
511	!! that has been actually added into the atmosphere. This is because both
512	!! taux and tautmsx are fully implicitly treated within compute_vdiff.
513	!! However, 'tautotx, tautoty' are not used in the actual numerical
514	!! computation in this module.
515	tautotx(:ncol) = taux(:ncol) + tautmsx(:ncol)
516	tautoty(:ncol) = tauy(:ncol) + tautmsy(:ncol)
517	else
518	ksrftms(:ncol) = 0.0_r8
519	tautotx(:ncol) = taux(:ncol)
520	tautoty(:ncol) = tauy(:ncol)
521	endif
522
523	!-------------------------------------------------------------------------------------!
524	!We are currenly using just water vapour flux at the surface, therefore we will use !
525	!just one constituent for the flux !
526	!-------------------------------------------------------------------------------------!
527	cflx(:pcols,1) = qfx(i,j) ! Surface constituent flux (kg/m2/s)
528
529	!Following variables are initialized to zero, they are the output from the "compute_eddy_diff" call
530	ustar8( :pcols) = 0.0_r8 ! Surface friction velocity (m/s)
531	pblh( :pcols) = 0.0_r8 ! Planetary boundary layer height (m )
532	ipbl( :pcols) = 0.0_r8 ! If 1, PBL is CL, while if 0, PBL is STL.
533	kpblh( :pcols) = 0.0_r8 ! Layer index containing PBL top within or at the base interface
534	wstarPBL(:pcols) = 0.0_r8 ! Convective velocity within PBL (m/s)
535
536
537	!!----------------------------------------------------------------------- !
538	!! Computation of eddy diffusivities - Select appropriate PBL scheme !
539	!!----------------------------------------------------------------------- !
540
541	select case (eddy_scheme)
542	case ( 'diag_TKE' )
543
544	!! ---------------------------------------------------------------- !
545	!! At first time step, have eddy_diff.F90:caleddy() use kvh=kvm=kvf !
546	!! This has to be done in compute_eddy_diff after kvf is calculated !
547	!! ---------------------------------------------------------------- !
548
549	kvinit = .false.
550	if(is_first_step) then
551	kvinit = .true.
552	endif
553	!! ---------------------------------------------- !
554	!! Get LW radiative heating out of physics buffer !
555	!! ---------------------------------------------- !
556
557	!! Retrieve eddy diffusivities for heat and momentum from physics buffer
558	!! from last timestep ( if first timestep, has been initialized by inidat.F90 )
559
560	kvm_in(:ncol,:) = kvm(:ncol,:)
561	kvh_in(:ncol,:) = kvh(:ncol,:)
562	call compute_eddy_diff( lchnk , pcols , pver , ncol , t8 , cloud(:,:,1) , ztodt , &
563	cloud(:,:,2), cloud(:,:,3), s8 , rpdel8 , cldn8 , qrl8 , wsedl8 , zm8 , zi8 , &
564	pmid8 , pint8 , u8 , v8 , taux , tauy , shflx , cflx(:,1), wstarent, &
565	nturb , ustar8 , pblh , kvm_in , kvh_in , kvm , kvh , kvq , cgh , &
566	cgs , tpert , qpert , wpert , tke8 , bprod , sprod , sfi , fqsatd , &
567	kvinit , tauresx , tauresy, ksrftms, ipbl(:), kpblh(:), wstarPBL(:) , turbtype , smaw )
568
569	!! ----------------------------------------------- !
570	!! Store TKE in pbuf for use by shallow convection !
571	!! ----------------------------------------------- !
572
573	!! Store updated kvh, kvm in pbuf to use here on the next timestep
574	do k = kts,kte+1
575	kflip = kte - k + 2
576
577	kvh3d(i,k,j) = kvh(1,kflip)
578	kvm3d(i,k,j) = kvm(1,kflip)
579	end do
580
581
582	end select
583
584	!!------------------------------------ !
585	!! Application of diffusivities !
586	!!------------------------------------ !
587	cloudtnd( :ncol,:,:) = cloud(:ncol,:,:)
588	stnd( :ncol,: ) = s8( :ncol,: )
589	wind_tends(:ncol,:,1) = u8( :ncol,: )
590	wind_tends(:ncol,:,2) = v8( :ncol,: )
591
592	!!------------------------------------------------------ !
593	!! Write profile output before applying diffusion scheme !
594	!!------------------------------------------------------ !
595
596	sl_prePBL(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
597	- ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
598	qt_prePBL(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
599	+ cloudtnd(:ncol,:pver,ixcldice)
600
601	call aqsat( t8, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
602	ftem_prePBL(:ncol,:) = cloud(:ncol,:,1)/ftem(:ncol,:)*100._r8
603
604	!! --------------------------------------------------------------------------------- !
605	!! Call the diffusivity solver and solve diffusion equation !
606	!! The final two arguments are optional function references to !
607	!! constituent-independent and constituent-dependent moleculuar diffusivity routines !
608	!! --------------------------------------------------------------------------------- !
609
610	!! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and
611	!! separately print out as diagnostic output, because these are different from
612	!! the explicit 'tautotx, tautoty' computed above.
613	!! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones.
614
615	if( any(fieldlist_wet) ) then
616	call compute_vdiff( lchnk, pcols, pver, pcnstmax, ncol, pmid8, pint8, rpdel8, t8, ztodt, &
617	taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, qmincg, &
618	fieldlist_wet, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
619	tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
620	compute_molec_diff, vd_lu_qdecomp)
621	end if
622
623	if( errstring .ne. '' ) then
624	call wrf_error_fatal(errstring)
625	endif
626
627	if( any( fieldlist_dry ) ) then
628	if( do_molec_diff ) then
629	errstring = "Design flaw: dry vdiff not currently supported with molecular diffusion"
630	call wrf_error_fatal(errstring)
631	end if
632
633	call compute_vdiff( lchnk, pcols, pver, pcnstmax, ncol, pmiddry8, pintdry8, rpdeldry8, t8, &
634	ztodt, taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, &
635	qmincg, fieldlist_dry, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
636	tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
637	compute_molec_diff, vd_lu_qdecomp)
638
639	if( errstring .ne. '' ) call wrf_error_fatal(errstring)
640
641	end if
642
643	! Store updated tauresx, tauresy to use here on the next timestep
644	tauresx2d(i,j) = tauresx(ncol)
645	tauresy2d(i,j) = tauresy(ncol)
646
647	!! -------------------------------------------------------- !
648	!! Diagnostics and output writing after applying PBL scheme !
649	!! -------------------------------------------------------- !
650	sl(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
651	- ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
652	qt(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
653	+ cloudtnd(:ncol,:pver,ixcldice)
654
655
656
657	!! --------------------------------------------------------------- !
658	!! Convert the new profiles into vertical diffusion tendencies. !
659	!! Convert KE dissipative heat change into "temperature" tendency. !
660	!! --------------------------------------------------------------- !
661
662	slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
663	qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt
664	stnd(:ncol,:) = ( stnd(:ncol,:) - s8(:ncol,:) ) * rztodt
665	wind_tends(:ncol,:,1) = ( wind_tends(:ncol,:,1) - u8(:ncol,:) ) * rztodt
666	wind_tends(:ncol,:,2) = ( wind_tends(:ncol,:,2) - v8(:ncol,:) ) * rztodt
667	cloudtnd(:ncol,:pver,:) = ( cloudtnd(:ncol,:pver,:) - cloud(:ncol,:pver,:) ) * rztodt
668
669	!! ----------------------------------------------------------- !
670	!! In order to perform 'pseudo-conservative varible diffusion' !
671	!! perform the following two stages: !
672	!! !
673	!! I. Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0' !
674	!! (2) 'sten' by 'slten', and !
675	!! (3) 'qlten = qiten = 0' !
676	!! !
677	!! II. Apply 'positive_moisture' !
678	!! !
679	!! ----------------------------------------------------------- !
680	if( eddy_scheme .eq. 'diag_TKE' .and. do_pseudocon_diff ) then
681	cloudtnd(:ncol,:pver,1) = qtten(:ncol,:pver)
682	stnd(:ncol,:pver) = slten(:ncol,:pver)
683	cloudtnd(:ncol,:pver,ixcldliq) = 0._r8
684	cloudtnd(:ncol,:pver,ixcldice) = 0._r8
685	cloudtnd(:ncol,:pver,ixnumliq) = 0._r8
686	cloudtnd(:ncol,:pver,ixnumice) = 0._r8
687
688	do ips = 1, ncol
689	do k = 1, pver
690	qv_pro(ips,k) = cloud(ips,k,1) + cloudtnd(ips,k,1) * ztodt
691	ql_pro(ips,k) = cloud(ips,k,ixcldliq) + cloudtnd(ips,k,ixcldliq) * ztodt
692	qi_pro(ips,k) = cloud(ips,k,ixcldice) + cloudtnd(ips,k,ixcldice) * ztodt
693	s_pro(ips,k) = s8(ips,k) + stnd(ips,k) * ztodt
694	t_pro(ips,k) = t8(ips,k) + (1._r8/cpair)stnd(ips,k) ztodt
695
696	end do
697	end do
698	call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(2), qmin(3), &
699	pdel8(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), &
700	qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1), &
701	cloudtnd(:ncol,pver:1:-1,1), cloudtnd(:ncol,pver:1:-1,ixcldliq), &
702	cloudtnd(:ncol,pver:1:-1,ixcldice), stnd(:ncol,pver:1:-1) )
703
704	end if
705
706	!! ----------------------------------------------------------------- !
707	!! Re-calculate diagnostic output variables after vertical diffusion !
708	!! ----------------------------------------------------------------- !
709	qv_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,1) + cloudtnd(:ncol,:pver,1) * ztodt
710	ql_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldliq) + cloudtnd(:ncol,:pver,ixcldliq) * ztodt
711	qi_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldice) + cloudtnd(:ncol,:pver,ixcldice) * ztodt
712
713	s_aft_PBL(:ncol,:pver) = s8(:ncol,:pver) + stnd(:ncol,:pver) * ztodt
714	t_aftPBL(:ncol,:pver) = ( s_aft_PBL(:ncol,:pver) - gravit*zm8(:ncol,:pver) ) / cpair
715
716	u_aft_PBL(:ncol,:pver) = u8(:ncol,:pver) + wind_tends(:ncol,:pver,1) * ztodt
717	v_aft_PBL(:ncol,:pver) = v8(:ncol,:pver) + wind_tends(:ncol,:pver,2) * ztodt
718
719	call aqsat( t_aftPBL, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
720	ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8
721
722	tten(:ncol,:pver) = ( t_aftPBL(:ncol,:pver) - t8(:ncol,:pver) ) * rztodt
723	rhten(:ncol,:pver) = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) ) * rztodt
724
725
726	!Post processing of the output from CAM
727	do k=kts,kte
728
729	kflip = kte-k+1
730
731	rublten(i,k,j) = wind_tends(1,kflip,1)
732	rvblten(i,k,j) = wind_tends(1,kflip,2)
733	rthblten(i,k,j) = stnd(1,kflip)/cpair/exner8(1,kflip)
734	!rqvblten tendency is making the code unstable. Assigned zero for now.Still looking into this issue...
735	!rqiblten tendency computed but not output (JD) ? also above comment may be no longer valid
736	rqvblten(i,k,j) = cloudtnd(1,kflip,1)/(1. - qv_curr(i,k,j))
737	rqcblten(i,k,j) = cloudtnd(1,kflip,2)/(1. - qv_curr(i,k,j))
738	tke_pbl (i,k,j) = tke8(1,kflip)
739	!if(k==kts+10)print*,i,j,k,rublten(i,k,j),rvblten(i,k,j), rqvblten(i,k,j),rqcblten(i,k,j) ,tke_pbl (i,k,j)
740
741	end do
742
743	kpbl2d(i,j) = int(kpblh(1))
744	pblh2d(i,j) = pblh( 1)
745	tpert2d(i,j) = tpert(1)
746	qpert2d(i,j) = qpert(1)
747	wpert2d(i,j) = wpert(1)
748
749	!End Post processing of the output from CAM
750	enddo !loop of i
751	enddo !loop of j
752	return
753	end subroutine camuwpbl
754
755
756
757
758	!-----------------------------------------------------------------------------------------
759	subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, &
760	dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten )
761	!! ------------------------------------------------------------------------------- !
762	!! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
763	!! force them to be larger than minimum value by (1) condensating water vapor !
764	!! into liquid or ice, and (2) by transporting water vapor from the very lower !
765	!! layer. '2._r8' is multiplied to the minimum values for safety. !
766	!! Update final state variables and tendencies associated with this correction. !
767	!! If any condensation happens, update (s,t) too. !
768	!! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
769	!! input tendencies. !
770	!! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer !
771	!! ------------------------------------------------------------------------------- !
772	!-----------------------------------------------------------------------------------------
773	implicit none
774	integer, intent(in) :: ncol, mkx
775	real(r8), intent(in) :: cp, xlv, xls
776	real(r8), intent(in) :: dt, qvmin, qlmin, qimin
777	real(r8), intent(in) :: dp(ncol,mkx)
778	real(r8), intent(inout) :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx)
779	real(r8), intent(inout) :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx)
780	integer i, k
781	real(r8) dql, dqi, dqv, sum, aa, dum
782
783	!! Modification : I should check whether this is exactly same as the one used in
784	!! shallow convection and cloud macrophysics.
785
786	do i = 1, ncol
787	do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface
788	dql = max(0._r8,1._r8*qlmin-ql(i,k))
789	dqi = max(0._r8,1._r8*qimin-qi(i,k))
790	qlten(i,k) = qlten(i,k) + dql/dt
791	qiten(i,k) = qiten(i,k) + dqi/dt
792	qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
793	sten(i,k) = sten(i,k) + xlv * (dql/dt) + xls * (dqi/dt)
794	ql(i,k) = ql(i,k) + dql
795	qi(i,k) = qi(i,k) + dqi
796	qv(i,k) = qv(i,k) - dql - dqi
797	s(i,k) = s(i,k) + xlv * dql + xls * dqi
798	t(i,k) = t(i,k) + (xlv * dql + xls * dqi)/cp
799	dqv = max(0._r8,1._r8*qvmin-qv(i,k))
800	qvten(i,k) = qvten(i,k) + dqv/dt
801	qv(i,k) = qv(i,k) + dqv
802	if( k .ne. 1 ) then
803	qv(i,k-1) = qv(i,k-1) - dqv*dp(i,k)/dp(i,k-1)
804	qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt
805	endif
806	qv(i,k) = max(qv(i,k),qvmin)
807	ql(i,k) = max(ql(i,k),qlmin)
808	qi(i,k) = max(qi(i,k),qimin)
809	end do
810	!! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
811	!! extracted from all the layers that has 'qv > 2*qvmin'. This fully
812	!! preserves column moisture.
813	if( dqv .gt. 1.e-20_r8 ) then
814	sum = 0._r8
815	do k = 1, mkx
816	if( qv(i,k) .gt. 2._r8qvmin ) sum = sum + qv(i,k)dp(i,k)
817	enddo
818	aa = dqv*dp(i,1)/max(1.e-20_r8,sum)
819	if( aa .lt. 0.5_r8 ) then
820	do k = 1, mkx
821	if( qv(i,k) .gt. 2._r8*qvmin ) then
822	dum = aa*qv(i,k)
823	qv(i,k) = qv(i,k) - dum
824	qvten(i,k) = qvten(i,k) - dum/dt
825	endif
826	enddo
827	else
828	write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion'
829	endif
830	endif
831	end do
832	return
833
834	end subroutine positive_moisture
835
836
837
838	!-----------------------------------------------------------------------------------------
839	subroutine camuwpblinit(rublten,rvblten,rthblten,rqvblten, &
840	restart,tke_pbl,grid_id, &
841	ids,ide,jds,jde,kds,kde, &
842	ims,ime,jms,jme,kms,kme, &
843	its,ite,jts,jte,kts,kte)
844	!!------------------------------------------------------------------!
845	!! Initialization of time independent fields for vertical diffusion !
846	!! Calls initialization routines for subsidiary modules !
847	!!----------------------------------------------------------------- !
848
849	!This subroutine is based on vertical_diffusion_init of CAM. This subroutine
850	!initializes variables for vertical diffusion subroutine calls. The layout
851	!is kept similar to the original CAM subroutine to facillitate future adaptations.
852	!-----------------------------------------------------------------------------------------
853
854	use eddy_diff, only : init_eddy_diff
855	use molec_diff, only : init_molec_diff
856	use trb_mtn_stress, only : init_tms
857	use diffusion_solver, only : init_vdiff, vdiff_select
858	use constituents, only : cnst_get_ind, cnst_get_type_byind, cnst_name
859	use module_cam_support, only : masterproc
860	use module_model_constants, only : epsq2
861
862	implicit none
863
864	!-------------------------------------------------------------------------------------!
865	!Input and output variables !
866	!-------------------------------------------------------------------------------------!
867	logical,intent(in) :: restart
868	integer,intent(in) :: grid_id
869	integer,intent(in) :: ids,ide,jds,jde,kds,kde
870	integer,intent(in) :: ims,ime,jms,jme,kms,kme
871	integer,intent(in) :: its,ite,jts,jte,kts,kte
872
873	real,dimension(ims:ime,kms:kme,jms:jme),intent(out) :: &
874	rublten, rvblten, rthblten, rqvblten,TKE_PBL
875
876	!-------------------------------------------------------------------------------------!
877	!Local Variables !
878	!-------------------------------------------------------------------------------------!
879	integer :: i,j,k,itf,jtf,ktf
880	integer :: ntop_eddy !! Top interface level to which eddy vertical diffusion is applied ( = 1 )
881	integer :: nbot_eddy !! Bottom interface level to which eddy vertical diffusion is applied ( = pver )
882	integer :: ntop_molec !! Top interface level to which molecular vertical diffusion is applied ( = 1 )
883	integer :: nbot_molec !! Bottom interface level to which molecular vertical diffusion is applied
884	integer :: pcnstmax ! number of constituents
885	character(128) :: errstring !! Error status for init_vdiff
886	real(r8) :: hypm(kte) !! reference state midpoint pressures
887
888
889	!Declare maximum number of constituents to be considered. pcnst is 7 in WRF but we are using 4 constituents
890	pcnstmax = 4
891
892	jtf = min(jte,jde-1)
893	ktf = min(kte,kde-1)
894	itf = min(ite,ide-1)
895
896	!Map CAM veritcal level variables
897	pver = ktf - kts + 1
898	pverp = pver + 1
899
900	!Initialize flags and add constituents
901	call vd_register()
902
903	!! ----------------------------------------------------------------- !
904	!! Get indices of cloud liquid and ice within the constituents array !
905	!! ----------------------------------------------------------------- !
906
907	call cnst_get_ind( 'CLDLIQ', ixcldliq )
908	call cnst_get_ind( 'CLDICE', ixcldice )
909	if( microp_scheme .eq. 'MG' ) then
910	call cnst_get_ind( 'NUMLIQ', ixnumliq )
911	call cnst_get_ind( 'NUMICE', ixnumice )
912	endif
913
914	if (masterproc) then
915	write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)'
916	end if
917
918	!! ---------------------------------------------------------------------------------------- !
919	!! Initialize molecular diffusivity module !
920	!! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa). !
921	!! Note that computing molecular diffusivities is a trivial expense, but constituent !
922	!! diffusivities depend on their molecular weights. Decomposing the diffusion matric !
923	!! for each constituent is a needless expense unless the diffusivity is significant. !
924	!! ---------------------------------------------------------------------------------------- !
925
926	ntop_molec = 1 !! Should always be 1
927	nbot_molec = 0 !! Should be set below about 70 km
928
929	!! ---------------------------------- !
930	!! Initialize eddy diffusivity module !
931	!! ---------------------------------- !
932
933	ntop_eddy = 1 !! No reason not to make this 1, if > 1, must be <= nbot_molec
934	nbot_eddy = pver !! Should always be pver
935
936	if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy
937
938	select case ( eddy_scheme )
939	case ( 'diag_TKE' )
940	if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme: UW Moist Turbulence Scheme by Bretherton and Park'
941
942	!! Check compatibility of eddy and shallow scheme
943	if( shallow_scheme .ne. 'UW' ) then
944	write(iulog,*) 'ERROR: shallow convection scheme ', shallow_scheme,' is incompatible with eddy scheme ', eddy_scheme
945	call wrf_error_fatal( 'convect_shallow_init: shallow_scheme and eddy_scheme are incompatible' )
946	endif
947
948	call init_eddy_diff( r8, pver, gravit, cpair, rair, zvir, latvap, latice, &
949	ntop_eddy, nbot_eddy, hypm, karman )
950
951	if( masterproc ) write(iulog,*) 'vertical_diffusion: nturb, ntop_eddy, nbot_eddy ', nturb, ntop_eddy, nbot_eddy
952	end select
953
954	!!-------------------------------------------------------------------------------------!
955	!! The vertical diffusion solver must operate !
956	!! over the full range of molecular and eddy diffusion !
957	!!-------------------------------------------------------------------------------------!
958
959	ntop = min(ntop_molec,ntop_eddy)
960	nbot = max(nbot_molec,nbot_eddy)
961
962	!! ------------------------------------------- !
963	!! Initialize turbulent mountain stress module !
964	!! ------------------------------------------- !
965
966	if( do_tms ) then
967	call init_tms( r8, tms_orocnst, karman, gravit, rair )
968	if (masterproc) then
969	write(iulog,*)'Using turbulent mountain stress module'
970	write(iulog,*)' tms_orocnst = ',tms_orocnst
971	end if
972	endif
973
974	!! ---------------------------------- !
975	!! Initialize diffusion solver module !
976	!! ---------------------------------- !
977
978	call init_vdiff( r8, pcnstmax, rair, gravit, fieldlist_wet, fieldlist_dry, errstring )
979	if( errstring .ne. '' ) call wrf_error_fatal( errstring )
980
981	!!-------------------------------------------------------------------------------------
982	!! Use fieldlist_wet to select the fields which will be diffused using moist mixing ratios ( all by default )
983	!! Use fieldlist_dry to select the fields which will be diffused using dry mixing ratios.
984	!!-------------------------------------------------------------------------------------
985
986	if( vdiff_select( fieldlist_wet, 'u' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'u' ) )
987	if( vdiff_select( fieldlist_wet, 'v' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'v' ) )
988	if( vdiff_select( fieldlist_wet, 's' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 's' ) )
989
990	do k = 1, pcnstmax
991	if( cnst_get_type_byind(k) .eq. 'wet' ) then
992	if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'q', k ) )
993	else
994	if( vdiff_select( fieldlist_dry, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_dry, 'q', k ) )
995	endif
996	enddo
997
998	!Initialize the tendencies
999	jtf=min0(jte,jde-1)
1000	ktf=min0(kte,kde-1)
1001	itf=min0(ite,ide-1)
1002
1003	if(.not.restart)then
1004	do j = jts , jtf
1005	do k = kts , ktf
1006	do i = its , itf
1007	tke_pbl(i,k,j) = epsq2
1008	rublten(i,k,j) = 0.
1009	rvblten(i,k,j) = 0.
1010	rthblten(i,k,j) = 0.
1011	rqvblten(i,k,j) = 0.
1012	enddo
1013	enddo
1014	enddo
1015	endif
1016	end subroutine camuwpblinit
1017
1018
1019
1020	!-----------------------------------------------------------------------------------------
1021	subroutine vd_register()
1022	!
1023	!This subroutine is based on the vd_register subroutine of CAM. Some additional
1024	!initializations are included in this subroutine.
1025	!
1026	!-----------------------------------------------------------------------------------------
1027
1028	use module_cam_esinti, only : esinti
1029	use physconst, only : mwh2o, cpwv, epsilo, latvap, latice, rh2o, cpair, tmelt
1030	use constituents, only : cnst_add
1031
1032	!Local Workspace
1033	integer :: mm
1034	logical :: moist_physics
1035
1036	!Following declarations are from CAM's stratiform.F90 module
1037	integer, parameter :: ncnstmax = 4 ! Number of constituents
1038
1039	character(len=8), dimension(ncnstmax), parameter :: cnst_names = &
1040	(/'CLDLIQ', 'CLDICE','NUMLIQ','NUMICE'/) ! Constituent names
1041
1042	if(vd_registered) return
1043
1044	!Set flags for the PBL scheme (these flags are obtained using phys_getopts in CAM)
1045	microp_scheme = 'MG' !Used for adding constituents
1046	eddy_scheme = 'diag_TKE' !Used for calling eddy scheme
1047
1048	!The following flag is deliberately set to UW for now.
1049	shallow_scheme = 'UW' !Eddy scheme is ONLY compaticle with 'UW' shallow scheme;
1050
1051	do_tms = .false. !To include stresses due to orography
1052	tms_orocnst = 1._r8 !Orography constant
1053
1054	!-------------------------------------------------------------------------------------!
1055	!Calls to add constituents (these calls are imported from in initindx.F90 in CAM) !
1056	! !
1057	! Register water vapor. !
1058	! This must be the first call to cnst_add so that water vapor is constituent 1. !
1059	!-------------------------------------------------------------------------------------!
1060
1061	moist_physics = .true. !This declaration is included to mimic CAM
1062
1063	if (moist_physics) then
1064	call cnst_add('Q', mwh2o, cpwv, 1.E-12_r8, mm, &
1065	longname='Specific humidity', readiv=.true. )
1066	else
1067	call cnst_add('Q', mwh2o, cpwv, 0.0_r8 , mm, &
1068	longname='Specific humidity', readiv=.false.)
1069	end if
1070
1071	!Following add constituent calls are imported from the stratiform.F90 in CAM
1072
1073	call cnst_add(cnst_names(1), mwdry, cpair, 0._r8, ixcldliq, &
1074	longname='Grid box averaged cloud liquid amount')
1075	call cnst_add(cnst_names(2), mwdry, cpair, 0._r8, ixcldice, &
1076	longname='Grid box averaged cloud ice amount' )
1077
1078	if ( microp_scheme .eq. 'MG' ) then
1079	call cnst_add(cnst_names(3), mwdry, cpair, 0._r8, ixnumliq, &
1080	longname='Grid box averaged cloud liquid number')
1081	call cnst_add(cnst_names(4), mwdry, cpair, 0._r8, ixnumice, &
1082	longname='Grid box averaged cloud ice number' )
1083	end if
1084
1085	!-------------------------------------------------------------------------------------!
1086	! Initialize the saturation vapor pressure look-up table... !
1087	! This could be moved to a master cam-init subroutine too if it is needed !
1088	! by more than one CAM parameterization. In CAM this is called from !
1089	! phys_init. !
1090	!-------------------------------------------------------------------------------------!
1091
1092	call esinti(epsilo, latvap, latice, rh2o, cpair, tmelt)
1093
1094	vd_registered = .true.
1095
1096	end subroutine vd_register
1097
1098	end module module_bl_camuwpbl_driver
1099

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

Download in other formats:

Original Format