Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_pbl_driver.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 78.0 KB

Line
1	!WRF:MEDIATION_LAYER:PHYSICS
2	!
3
4	MODULE module_pbl_driver
5	CONTAINS
6
7	!------------------------------------------------------------------
8	SUBROUTINE pbl_driver( &
9	itimestep,dt,u_frame,v_frame &
10	,bldt,curr_secs,adapt_step_flag &
11	,rublten,rvblten,rthblten &
12	,tsk,xland,znt,ht &
13	,ust,pblh,hfx,qfx,grdflx &
14	,u_phy,v_phy,th_phy,rho &
15	,p_phy,pi_phy,p8w,t_phy,dz8w,z &
16	,exch_h,exch_m,akhs,akms &
17	,thz0,qz0,uz0,vz0,qsfc,f &
18	,lowlyr,u10,v10,t2 &
19	,psim,psih,gz1oz0, wspd,br,chklowq &
20	,bl_pbl_physics, ra_lw_physics, dx &
21	,stepbl,warm_rain &
22	,kpbl,mixht,ct,lh,snow,xice &
23	,znu, znw, mut, p_top &
24	! OPTIONAL for TEMF scheme
25	,te_temf,km_temf,kh_temf &
26	,shf_temf,qf_temf,uw_temf,vw_temf &
27	,hd_temf,lcl_temf,hct_temf &
28	,wupd_temf,mf_temf,thup_temf,qtup_temf,qlup_temf &
29	,exch_temf,cf3d_temf,cfm_temf &
30	,flhc,flqc &
31	!
32	,qke,tsq,qsq,cov,rmol,ch,qcg,grav_settling &
33	#if (NMM_CORE==1)
34	,DISHEAT &
35	#endif
36	,ids,ide, jds,jde, kds,kde &
37	,ims,ime, jms,jme, kms,kme &
38	,i_start,i_end, j_start,j_end, kts,kte, num_tiles &
39	! Optional
40	,hol, mol, regime &
41	! Optional gravity-wave drag
42	,gwd_opt &
43	,dtaux3d,dtauy3d &
44	,dusfcg,dvsfcg,var2d,oc12d &
45	,oa1,oa2,oa3,oa4,ol1,ol2,ol3,ol4 &
46	! Optional moisture tracers
47	,qv_curr, qc_curr, qr_curr &
48	,qi_curr, qs_curr, qg_curr &
49	,rqvblten,rqcblten,rqiblten &
50	,rqrblten,rqsblten,rqgblten &
51	! Optional moisture tracer flags
52	,f_qv,f_qc,f_qr &
53	,f_qi,f_qs,f_qg &
54	! variables added for BEP
55	,frc_urb2d &
56	,a_u_bep,a_v_bep,a_t_bep,a_q_bep &
57	,b_u_bep,b_v_bep,b_t_bep,b_q_bep &
58	,sf_bep,vl_bep &
59	,sf_sfclay_physics,sf_urban_physics &
60	,tke_pbl,el_pbl &
61	#if (NMM_CORE != 1)
62	,wu_tur,wv_tur,wt_tur,wq_tur &
63	#endif
64	,a_e_bep,b_e_bep,dlg_bep &
65	,dl_u_bep &
66	! Wind Turbine Parameterizations
67	,phb,xlat_u,xlong_u,xlat_v,xlong_v,id &
68	! variables required for camuwpbl scheme
69	, z_at_w,cldfra,rthratenlw,tauresx2d,tauresy2d &
70	, tpert2d,qpert2d,wpert2d &
71	)
72	!------------------------------------------------------------------
73	#if (EM_CORE==1)
74	USE module_state_description, ONLY : &
75	YSUSCHEME,MRFSCHEME,GFSSCHEME,MYJPBLSCHEME,ACMPBLSCHEME,&
76	QNSEPBLSCHEME,MYNNPBLSCHEME2,MYNNPBLSCHEME3,BOULACSCHEME,&
77	CAMUWPBLSCHEME,BEPSCHEME,BEP_BEMSCHEME,MYJSFCSCHEME, &
78	TEMFPBLSCHEME, &
79	p_qi,param_first_scalar
80	USE module_wind_generic, ONLY : windspec
81	#else
82	USE module_state_description, ONLY : &
83	YSUSCHEME,MRFSCHEME,GFSSCHEME,MYJPBLSCHEME,ACMPBLSCHEME &
84	, QNSEPBLSCHEME, p_qi,param_first_scalar &
85	, TEMFPBLSCHEME &
86	, MYJSFCSCHEME
87	#endif
88
89	USE module_model_constants
90
91	! *** add new modules of schemes here
92
93	USE module_bl_myjpbl
94	USE module_bl_qnsepbl
95	USE module_bl_ysu
96	USE module_bl_mrf
97	USE module_bl_gfs
98	USE module_bl_acm
99	USE module_bl_gwdo
100	USE module_bl_myjurb
101	USE module_bl_boulac
102	USE module_bl_camuwpbl_driver, only:CAMUWPBL
103	USE module_bl_temf
104	#if (EM_CORE==1)
105	USE module_bl_mynn
106	USE module_wind_fitch
107	#endif
108
109	! This driver calls subroutines for the PBL parameterizations.
110	!
111	! pbl scheme:
112	! 1. ysupbl
113	! 2. myjpbl
114	! 4. qnsepbl
115	! 5. mynnpbl2
116	! 6. mynnpbl3
117	! 7. acmpbl
118	! 8. boulacpbl
119	! 9. camuwpbl
120	! 10. temfpbl
121	! 99. mrfpbl
122	!
123	!------------------------------------------------------------------
124	IMPLICIT NONE
125	!======================================================================
126	! Grid structure in physics part of WRF
127	!----------------------------------------------------------------------
128	! The horizontal velocities used in the physics are unstaggered
129	! relative to temperature/moisture variables. All predicted
130	! variables are carried at half levels except w, which is at full
131	! levels. Some arrays with names (*8w) are at w (full) levels.
132	!
133	!----------------------------------------------------------------------
134	! In WRF, kms (smallest number) is the bottom level and kme (largest
135	! number) is the top level. In your scheme, if 1 is at the top level,
136	! then you have to reverse the order in the k direction.
137	!
138	! kme - half level (no data at this level)
139	! kme ----- full level
140	! kme-1 - half level
141	! kme-1 ----- full level
142	! .
143	! .
144	! .
145	! kms+2 - half level
146	! kms+2 ----- full level
147	! kms+1 - half level
148	! kms+1 ----- full level
149	! kms - half level
150	! kms ----- full level
151	!
152	!======================================================================
153	! Definitions
154	!-----------
155	! Rho_d dry density (kg/m^3)
156	! Theta_m moist potential temperature (K)
157	! Qv water vapor mixing ratio (kg/kg)
158	! Qc cloud water mixing ratio (kg/kg)
159	! Qr rain water mixing ratio (kg/kg)
160	! Qi cloud ice mixing ratio (kg/kg)
161	! Qs snow mixing ratio (kg/kg)
162	!-----------------------------------------------------------------
163	!-- RUBLTEN U tendency due to
164	! PBL parameterization (m/s^2)
165	!-- RVBLTEN V tendency due to
166	! PBL parameterization (m/s^2)
167	!-- RTHBLTEN Theta tendency due to
168	! PBL parameterization (K/s)
169	!-- RQVBLTEN Qv tendency due to
170	! PBL parameterization (kg/kg/s)
171	!-- RQCBLTEN Qc tendency due to
172	! PBL parameterization (kg/kg/s)
173	!-- RQIBLTEN Qi tendency due to
174	! PBL parameterization (kg/kg/s)
175	!-- id WRF grid id (optional, only needed by turbine drag schemes)
176	!-- itimestep number of time steps
177	!-- GLW downward long wave flux at ground surface (W/m^2)
178	!-- GSW downward short wave flux at ground surface (W/m^2)
179	!-- EMISS surface emissivity (between 0 and 1)
180	!-- TSK surface temperature (K)
181	!-- TMN soil temperature at lower boundary (K)
182	!-- XLAND land mask (1 for land, 2 for water)
183	!-- ZNT roughness length (m)
184	!-- MAVAIL surface moisture availability (between 0 and 1)
185	!-- UST u* in similarity theory (m/s)
186	!-- MOL T* (similarity theory) (K)
187	!-- HOL PBL height over Monin-Obukhov length
188	!-- PBLH PBL height (m)
189	!-- CAPG heat capacity for soil (J/K/m^3)
190	!-- THC thermal inertia (Cal/cm/K/s^0.5)
191	!-- SNOWC flag indicating snow coverage (1 for snow cover)
192	!-- HFX upward heat flux at the surface (W/m^2)
193	!-- QFX upward moisture flux at the surface (kg/m^2/s)
194	!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
195	!-- akhs sfc exchange coefficient of heat/moisture from MYJ
196	!-- akms sfc exchange coefficient of momentum from MYJ
197	!-- tke_pbl turbulence kinetic energy from PBL schemes (m^2/s^2)
198	!-- el_pbl length scale from PBL schemes (m)
199	!-- wu_tur turbulent flux of momentum (x) (m^2/s^2)
200	!-- wv_tur turbulent flux of momentum (y) (m^2/s^2)
201	!-- wt_tur turbulent flux of potential temperature (K m/s)
202	!-- wq_tur turbulent flux of water vapor (- m/s)
203	!-- te_temf Total energy from TEMF BL scheme
204	!-- km_temf Exchange coefficient for momentum from TEMF BL scheme
205	!-- kh_temf Exchange coefficient for heat from TEMF BL scheme
206	!-- shf_temf Sensible heat flux from TEMF BL scheme
207	!-- qf_temf Water vapor flux from TEMF BL scheme
208	!-- uw_temf Momentum flux in U direction from TEMF BL scheme
209	!-- vw_temf Momentum flux in V direction from TEMF BL scheme
210	!-- wupd_temf Updraft velocity from TEMF BL scheme
211	!-- mf_temf Mass flux from TEMF BL scheme
212	!-- thup_temf Updraft thetal from TEMF BL scheme
213	!-- qtup_temf Updraft qt from TEMF BL scheme
214	!-- qlup_temf Updraft ql from TEMF BL scheme
215	!-- cf3d_temf 3D cloud fraction from TEMF PBL
216	!-- cfm_temf Column cloud fraction from TEMF PBL
217	!-- exch_temf Surface exchange coefficient (as for moisture) from TEMF surface layer scheme
218	!-- flhc Surface exchange coefficient for heat (for TEMF)
219	!-- flqc Surface exchange coefficient for moisture (for TEMF)
220	!-- thz0 potential temperature at roughness length (K)
221	!-- uz0 u wind component at roughness length (m/s)
222	!-- vz0 v wind component at roughness length (m/s)
223	!-- qsfc specific humidity at lower boundary (kg/kg)
224	!-- th2 diagnostic 2-m theta from surface layer and lsm
225	!-- t2 diagnostic 2-m temperature from surface layer and lsm
226	!-- q2 diagnostic 2-m mixing ratio from surface layer and lsm
227	!-- lowlyr index of lowest model layer above ground
228	!-- rr dry air density (kg/m^3)
229	!-- u_phy u-velocity interpolated to theta points (m/s)
230	!-- v_phy v-velocity interpolated to theta points (m/s)
231	!-- th_phy potential temperature (K)
232	!-- p_phy pressure (Pa)
233	!-- pi_phy exner function (dimensionless)
234	!-- p8w pressure at full levels (Pa)
235	!-- t_phy temperature (K)
236	!-- dz8w dz between full levels (m)
237	!-- z height above sea level (m)
238	!-- DX horizontal space interval (m)
239	!-- DT time step (second)
240	!-- n_moist number of moisture species
241	!-- PSFC pressure at the surface (Pa)
242	!-- TSLB
243	!-- ZS
244	!-- DZS
245	!-- num_soil_layers number of soil layer
246	!-- IFSNOW ifsnow=1 for snow-cover effects
247	!
248	!-- P_QV species index for water vapor
249	!-- P_QC species index for cloud water
250	!-- P_QR species index for rain water
251	!-- P_QI species index for cloud ice
252	!-- P_QS species index for snow
253	!-- P_QG species index for graupel
254	!-- ids start index for i in domain
255	!-- ide end index for i in domain
256	!-- jds start index for j in domain
257	!-- jde end index for j in domain
258	!-- kds start index for k in domain
259	!-- kde end index for k in domain
260	!-- ims start index for i in memory
261	!-- ime end index for i in memory
262	!-- jms start index for j in memory
263	!-- jme end index for j in memory
264	!-- kms start index for k in memory
265	!-- kme end index for k in memory
266	!-- jts start index for j in tile
267	!-- jte end index for j in tile
268	!-- kts start index for k in tile
269	!-- kte end index for k in tile
270	!
271	!******************************************************************
272	!------------------------------------------------------------------
273	!
274
275
276	INTEGER, INTENT(IN ) :: bl_pbl_physics, ra_lw_physics,sf_sfclay_physics,sf_urban_physics
277
278	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
279	ims,ime, jms,jme, kms,kme, &
280	kts,kte, num_tiles
281
282	INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
283	& i_start,i_end,j_start,j_end
284
285	INTEGER, INTENT(IN ) :: itimestep,STEPBL
286	INTEGER, DIMENSION( ims:ime , jms:jme ), &
287	INTENT(IN ) :: LOWLYR
288	!
289	LOGICAL, INTENT(IN ) :: warm_rain
290	#if (NMM_CORE==1)
291	LOGICAL, INTENT(IN ) :: DISHEAT ! (for HWRF)
292	#endif
293
294	REAL, DIMENSION( kms:kme ), &
295	OPTIONAL, INTENT(IN ) :: znu, &
296	znw
297	!
298	REAL, INTENT(IN ) :: DT,DX
299	REAL, INTENT(IN ),OPTIONAL :: bldt
300	REAL, INTENT(IN ),OPTIONAL :: curr_secs
301	LOGICAL, INTENT(IN ),OPTIONAL :: adapt_step_flag
302	! Optional for Wind Turbine Parameterizations
303	REAL, DIMENSION( ims:ime, kms:kme ,jms:jme ), &
304	INTENT(IN), OPTIONAL :: phb
305	REAL, DIMENSION( ims:ime, jms:jme ), &
306	INTENT(IN), OPTIONAL :: xlat_u,xlong_u,xlat_v,xlong_v
307
308	!
309	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
310	INTENT(IN ) :: p_phy, &
311	pi_phy, &
312	p8w, &
313	rho, &
314	t_phy, &
315	u_phy, &
316	v_phy, &
317	dz8w, &
318	z, &
319	th_phy
320	!3D Variables for camuwpbl scheme
321	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
322	INTENT(IN ), OPTIONAL :: z_at_w, &
323	cldfra, &
324	rthratenlw
325
326	!2D Variables required by camuwpbl scheme
327	REAL, DIMENSION( ims:ime , jms:jme ), &
328	INTENT(INOUT ), OPTIONAL :: tauresx2d, &
329	tauresy2d, &
330	tpert2d, &
331	qpert2d, &
332	wpert2d
333	!
334	!
335	REAL, DIMENSION( ims:ime , jms:jme ), &
336	INTENT(IN ) :: XLAND, &
337	HT, &
338	PSIM, &
339	PSIH, &
340	GZ1OZ0, &
341	BR, &
342	F, &
343	CHKLOWQ
344	!
345	REAL, DIMENSION( ims:ime, jms:jme ) , &
346	INTENT(INOUT) :: TSK, &
347	UST, &
348	PBLH, &
349	HFX, &
350	QFX, &
351	ZNT, &
352	QSFC, &
353	AKHS, &
354	AKMS, &
355	MIXHT, &
356	QZ0, &
357	THZ0, &
358	UZ0, &
359	VZ0, &
360	CT, &
361	GRDFLX, &
362	U10, &
363	V10, &
364	T2, &
365	WSPD
366
367	!
368	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
369	INTENT(INOUT) :: RUBLTEN, &
370	RVBLTEN, &
371	RTHBLTEN, &
372	EXCH_H,EXCH_M,TKE_PBL
373	#if (NMM_CORE != 1)
374	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
375	INTENT(OUT) :: WU_TUR,WV_TUR,WT_TUR,WQ_TUR
376	!
377	#endif
378
379	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
380	&OPTIONAL, INTENT(INOUT) :: &
381	& qke,tsq,qsq,cov!,k_m,k_h,k_q
382	INTEGER, OPTIONAL :: id
383
384	REAL, DIMENSION( ims:ime , jms:jme ), &
385	&OPTIONAL, INTENT(IN) :: &
386	& qcg, rmol, ch
387
388
389	INTEGER, OPTIONAL, INTENT(IN) :: grav_settling
390
391
392
393
394	!
395	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
396	INTENT(OUT) :: EL_PBL
397
398	REAL , INTENT(IN ) :: u_frame, &
399	v_frame
400	!
401
402	INTEGER, DIMENSION( ims:ime , jms:jme ), &
403	INTENT(INOUT) :: KPBL
404
405	REAL, DIMENSION( ims:ime , jms:jme ), &
406	INTENT(IN) :: XICE, SNOW, LH
407
408	! Bep changes: variable added for urban
409	real, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::FRC_URB2D ! URBAN Landuse fraction
410	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep ! Implicit component for the momemtum in X-direction
411	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep ! Implicit component for the momemtum in Y-direction
412	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep ! Implicit component for the Pot. Temp.
413	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep ! Implicit component for Moisture
414
415	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep ! Implicit component for the TKE
416	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep ! Explicit component for the momemtum in X-direction
417	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep ! Explicit component for the momemtum in Y-direction
418	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep ! Explicit component for the Pot. Temp.
419	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep ! Explicit component for Moisture
420
421	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep ! Explicit component for the TKE
422
423	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep ! Height above ground (L_ground in formula (24) of the BLM paper).
424	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep ! Length scale (lb in formula (22) ofthe BLM paper).
425	! urban surface and volumes
426	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep ! surfaces
427	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep ! volumes
428
429	! Bep changes end
430	! New variables for TEMF scheme
431	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
432	INTENT(INOUT) :: te_temf
433	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
434	INTENT( OUT) :: km_temf, kh_temf, &
435	shf_temf,qf_temf,uw_temf,vw_temf, &
436	wupd_temf,mf_temf,thup_temf,qtup_temf, &
437	qlup_temf,cf3d_temf
438	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
439	INTENT(INOUT) :: flhc,flqc
440	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
441	INTENT( OUT) :: hd_temf, lcl_temf, hct_temf, cfm_temf
442	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
443	INTENT(INOUT) :: exch_temf
444
445	!
446	!
447	! Optional
448	!
449	!
450	! Flags relating to the optional tendency arrays declared above
451	! Models that carry the optional tendencies will provdide the
452	! optional arguments at compile time; these flags all the model
453	! to determine at run-time whether a particular tracer is in
454	! use or not.
455	!
456	LOGICAL, INTENT(IN), OPTIONAL :: &
457	f_qv &
458	,f_qc &
459	,f_qr &
460	,f_qi &
461	,f_qs &
462	,f_qg
463
464	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
465	OPTIONAL, INTENT(INOUT) :: &
466	! optional moisture tracers
467	! 2 time levels; if only one then use CURR
468	qv_curr, qc_curr, qr_curr &
469	,qi_curr, qs_curr, qg_curr &
470	,rqvblten,rqcblten,rqrblten &
471	,rqiblten,rqsblten,rqgblten
472
473	REAL, DIMENSION( ims:ime, jms:jme ) , &
474	OPTIONAL , &
475	INTENT(INOUT) :: HOL, &
476	MOL, &
477	REGIME
478	REAL, DIMENSION( ims:ime, jms:jme ) , &
479	OPTIONAL , &
480	INTENT(IN) :: mut
481	!
482	INTEGER, OPTIONAL, INTENT(IN) :: gwd_opt
483	REAL, OPTIONAL, INTENT(IN) :: p_top
484	!
485	real, dimension( ims:ime, kms:kme, jms:jme ) , &
486	optional , &
487	intent(inout ) :: dtaux3d, &
488	dtauy3d
489	!
490	real, dimension( ims:ime, jms:jme ) , &
491	optional , &
492	intent(inout ) :: dusfcg, &
493	dvsfcg
494	!
495	real, dimension( ims:ime, jms:jme ) , &
496	optional , &
497	intent(in ) :: var2d, &
498	oc12d, &
499	oa1,oa2,oa3,oa4, &
500	ol1,ol2,ol3,ol4
501
502	! LOCAL VAR
503
504	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) ::v_phytmp
505	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) ::u_phytmp
506
507	REAL, DIMENSION( ims:ime, jms:jme ) :: TSKOLD, &
508	USTOLD, &
509	ZNTOLD, &
510	ZOL, &
511	PSFC
512	! make these allocatable depending on the setting of idiff
513	! Typically, we try to avoide allocating and deallocating local storage like this
514	! so as not to fragment the stack. But at this point, the idiff = 1 case is disabled
515	! (set to 0 for all cases) and has to be set manually by users who want to work with
516	! it. When it becomes a more standard option, this should be redone, either defining
517	! these as state with package clauses to turn them on and off and passing them in,
518	! or pass in an integer flag that can be used to dimension the arrays to 1:1:1 as
519	! local variables. JM 20100316
520
521	REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_u ! Implicit component for the momemtum in X-direction
522	REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_v ! Implicit component for the momemtum in Y-direction
523	REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_t ! Implicit component for the Pot. Temp.
524	REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_q ! Implicit component for the water vapor
525
526	REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_u ! Explicit component for the momemtum in X-direction
527	REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_v ! Explicit component for the momemtum in Y-direction
528	REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_t ! Explicit component for the Pot. Temp.
529	REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_q ! Explicit component for the water vapor
530
531	REAL, ALLOCATABLE, DIMENSION( :, :, : )::sf ! surfaces
532	REAL, ALLOCATABLE, DIMENSION( :, :, : )::vl ! volumes
533
534	REAL :: DTMIN,DTBL
535	!
536
537	INTEGER :: initflag
538	!
539
540
541
542	INTEGER :: i,J,K,NK,jj,ij,its,ite,jts,jte
543	LOGICAL :: radiation
544	LOGICAL :: flag_bep
545	LOGICAL :: flag_myjsfc
546	LOGICAL :: flag_qv, flag_qc, flag_qr, flag_qi, flag_qs, flag_qg
547	CHARACTER*256 :: message
548	REAL :: next_bl_time
549	LOGICAL :: run_param
550	LOGICAL :: do_adapt
551	integer iu_bep,iurb,idiff
552	real seamask,thsk,zzz,unew,vnew,tnew,qnew,umom,vmom
553
554	!------------------------------------------------------------------
555	!
556	!!!!!!!if using BEP set flag_bep to true
557	SELECT CASE(sf_urban_physics)
558	#if (NMM_CORE != 1)
559	CASE (BEPSCHEME)
560	flag_bep=.true.
561	CASE (BEP_BEMSCHEME)
562	flag_bep=.true.
563	#endif
564	CASE DEFAULT
565	flag_bep=.false.
566	END SELECT
567
568	SELECT CASE(sf_sfclay_physics)
569	CASE (MYJSFCSCHEME)
570	flag_myjsfc=.true.
571	CASE DEFAULT
572	flag_myjsfc=.false.
573	END SELECT
574	!
575	flag_qv = .FALSE. ; IF ( PRESENT( F_QV ) ) flag_qv = F_QV
576	flag_qc = .FALSE. ; IF ( PRESENT( F_QC ) ) flag_qc = F_QC
577	flag_qr = .FALSE. ; IF ( PRESENT( F_QR ) ) flag_qr = F_QR
578	flag_qi = .FALSE. ; IF ( PRESENT( F_QI ) ) flag_qi = F_QI
579	flag_qs = .FALSE. ; IF ( PRESENT( F_QS ) ) flag_qs = F_QS
580	flag_qg = .FALSE. ; IF ( PRESENT( F_QG ) ) flag_qg = F_QG
581
582	!print *,flag_qv, flag_qc, flag_qr, flag_qi, flag_qs, flag_qg,' flag_qv, flag_qc, flag_qr, flag_qi, flag_qs, flag_qg'
583	!print *,f_qv, f_qc, f_qr, f_qi, f_qs, f_qg,' f_qv, f_qc, f_qr, f_qi, f_qs, f_qg'
584
585	if (bl_pbl_physics .eq. 0) return
586	! RAINBL in mm (Accumulation between PBL calls)
587	!
588	! Modified for adaptive time step
589	!
590	IF ( (itimestep .EQ. 1) .OR. (MOD(itimestep,STEPBL) .EQ. 0) ) THEN
591	run_param = .TRUE.
592	ELSE
593	run_param = .FALSE.
594	ENDIF
595
596	IF (PRESENT(adapt_step_flag)) THEN
597	IF ((adapt_step_flag)) THEN
598	IF ( (itimestep .EQ. 1) .OR. (bldt .EQ. 0) .OR. &
599	( CURR_SECS + dt >= ( INT( CURR_SECS / ( bldt * 60 ) + 1 ) * bldt * 60) ) ) THEN
600	run_param = .TRUE.
601	ELSE
602	run_param = .FALSE.
603	ENDIF
604	ENDIF
605	ENDIF
606
607	IF (run_param) THEN
608	radiation = .false.
609	IF (ra_lw_physics .gt. 0) radiation = .true.
610
611	!----
612	! CALCULATE CONSTANT
613
614	DTMIN=DT/60.
615	! PBL schemes need PBL time step for updates
616
617	if (PRESENT(adapt_step_flag)) then
618	if (adapt_step_flag) then
619	do_adapt = .TRUE.
620	else
621	do_adapt = .FALSE.
622	endif
623	else
624	do_adapt = .FALSE.
625	endif
626
627	if (PRESENT(BLDT)) then
628	if (bldt .eq. 0) then
629	DTBL = dt
630	ELSE
631	if (do_adapt) then
632	call wrf_message("WARNING: When using an adaptive time-step the boundary layer"// &
633	" time-step should be 0 (i.e., equivalent to model time-step). "// &
634	"In order to proceed, for boundary layer calculations, the "// &
635	"boundary layer time-step"// &
636	" will be rounded to the nearest minute, possibly resulting in"// &
637	" innacurate results.")
638	DTBL=bldt*60
639	else
640	DTBL=DT*STEPBL
641	endif
642	endif
643	else
644	DTBL=DT*STEPBL
645	endif
646
647	#if (NMM_CORE != 1)
648	!!Alberto. If idiff=1, then compute the tendencies at the end in the routine diff3d
649	!!Alberto. If idiff=0, then compute the tendencies in each PBL routine (standard version).
650	idiff=0
651	#else
652	! Added this else clause so that idiff is always initialized regardless of which core we are using
653	idiff=0
654	#endif
655
656	IF ( idiff .EQ. 1 ) THEN
657	ALLOCATE (a_u(ims:ime,kms:kme,jms:jme)) ! Implicit component for the momemtum in X-direction
658	ALLOCATE (a_v(ims:ime,kms:kme,jms:jme)) ! Implicit component for the momemtum in Y-direction
659	ALLOCATE (a_t(ims:ime,kms:kme,jms:jme)) ! Implicit component for the Pot. Temp.
660	ALLOCATE (a_q(ims:ime,kms:kme,jms:jme)) ! Implicit component for the water vapor
661	ALLOCATE (b_u(ims:ime,kms:kme,jms:jme)) ! Explicit component for the momemtum in X-direction
662	ALLOCATE (b_v(ims:ime,kms:kme,jms:jme)) ! Explicit component for the momemtum in Y-direction
663	ALLOCATE (b_t(ims:ime,kms:kme,jms:jme)) ! Explicit component for the Pot. Temp.
664	ALLOCATE (b_q(ims:ime,kms:kme,jms:jme)) ! Explicit component for the water vapor
665	ALLOCATE (sf(ims:ime,kms:kme,jms:jme) ) ! surfaces
666	ALLOCATE (vl(ims:ime,kms:kme,jms:jme) ) ! volumes
667	ENDIF
668
669
670	! SAVE OLD VALUES
671
672	!$OMP PARALLEL DO &
673	!$OMP PRIVATE ( ij,i,j,k )
674
675	DO ij = 1 , num_tiles
676	DO j=j_start(ij),j_end(ij)
677	DO i=i_start(ij),i_end(ij)
678	TSKOLD(i,j)=TSK(i,j)
679	USTOLD(i,j)=UST(i,j)
680	ZNTOLD(i,j)=ZNT(i,j)
681
682	! REVERSE ORDER IN THE VERTICAL DIRECTION
683
684	! testing change later
685
686	DO k=kts,kte
687	v_phytmp(i,k,j)=v_phy(i,k,j)+v_frame
688	u_phytmp(i,k,j)=u_phy(i,k,j)+u_frame
689	ENDDO
690
691	! PSFC : in Pa
692
693	PSFC(I,J)=p8w(I,kms,J)
694
695	DO k=kts,min(kte+1,kde)
696	RTHBLTEN(I,K,J)=0.
697	RUBLTEN(I,K,J)=0.
698	RVBLTEN(I,K,J)=0.
699	IF ( PRESENT( RQCBLTEN )) RQCBLTEN(I,K,J)=0.
700	IF ( PRESENT( RQVBLTEN )) RQVBLTEN(I,K,J)=0.
701	ENDDO
702
703	IF (flag_QI .AND. PRESENT(RQIBLTEN) ) THEN
704	DO k=kts,min(kte+1,kde)
705	RQIBLTEN(I,K,J)=0.
706	ENDDO
707	ENDIF
708	ENDDO
709	ENDDO
710
711	#if (NMM_CORE != 1)
712	IF ( idiff.eq.1 ) THEN
713	!Alberto
714	! Here we should put a switch:
715	! if we do not use BEP, just initialize the values of the a, and b to zero, except for the lowest model level,
716	! where the heat and momentum fluxes are introduced
717	! if we use BEP, use the values computed by BEP weigthed by the urban fraction.
718	! for the moment I put a simple "if" but it should be changed to a more elegant way after(using the CASE, maybe?)
719	!!!!!!!!!!!!!!
720	! This is the more general way to go, however some of these variables (e. g. a_t, a_q) may be zero nearly all time.
721	! if we need to be as tight as possible for the memory we can think how to reduce the storage.
722	!!!!!!!!!!!!!!!!!!
723	! This stuff is becoming large, probably better to put in a subroutine
724	!!!!!!!!!!!!!!!!!!!
725	! preparing the a, b, sf, and vl terms
726
727	IF (flag_bep) THEN
728	do j=j_start(ij),j_end(ij)
729	do i=i_start(ij),i_end(ij)
730	do k=kts,kte
731	a_u(i,k,j)=a_u_bep(i,k,j)
732	a_v(i,k,j)=a_v_bep(i,k,j)
733	a_t(i,k,j)=a_t_bep(i,k,j)
734	a_q(i,k,j)=a_q_bep(i,k,j)
735	b_u(i,k,j)=b_u_bep(i,k,j)
736	b_v(i,k,j)=b_v_bep(i,k,j)
737	b_t(i,k,j)=b_t_bep(i,k,j)
738	b_q(i,k,j)=b_q_bep(i,k,j)
739	vl(i,k,j)=vl_bep(i,k,j)
740	sf(i,k,j)=sf_bep(i,k,j)
741	enddo
742	sf(i,kte+1,j)=1.
743	enddo
744	enddo
745	ELSE
746	do j=j_start(ij),j_end(ij)
747	do i=i_start(ij),i_end(ij)
748	do k=kts,kte
749	a_u(i,k,j)=0.
750	a_v(i,k,j)=0.
751	a_t(i,k,j)=0.
752	a_q(i,k,j)=0.
753	b_u(i,k,j)=0.
754	b_v(i,k,j)=0.
755	b_t(i,k,j)=0.
756	b_q(i,k,j)=0.
757	vl(i,k,j)=1.
758	sf(i,k,j)=1.
759	enddo
760	sf(i,kte+1,j)=1.
761	! fix the values for the surface fluxes (source/sink terms)
762	! here for momentum the resolution is done implicitely
763	! while for heat and moisture is done explicitely
764	a_u(i,1,j)=(-ust(I,J)ust(I,J))/dz8w(i,1,j)/((u_phy(i,1,j)2+v_phy(i,1,j)2.)*.5)
765	a_v(i,1,j)=(-ust(I,J)ust(I,J))/dz8w(i,1,j)/((u_phy(i,1,j)2+v_phy(i,1,j)2.)*.5)
766	b_t(i,1,j)=hfx(i,j)/rho(i,1,j)/CP/dz8w(i,1,j)
767	b_q(i,1,j)=qfx(i,j)/rho(i,1,j)/dz8w(i,1,j)
768	!! if you want to solve also the heat and moisture fluxes implicitely, uncomment the following lines.
769	!! (of course you will need the values of thz0,qz0,akhs).
770	! b_t(i,1,j)=akhs(i,j)*(thz0(I,J))/dz8w(i,1,j)
771	! b_q(i,1,j)=akhs(i,j)*(qz0(I,J))/dz8w(i,1,j)
772	! a_t(i,1,j)=-1.*akhs(i,j)/dz8w(i,1,j)
773	! a_q(i,1,j)=-1.*akhs(i,j)/dz8w(i,1,j)
774	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
775	enddo
776	enddo
777
778	ENDIF
779
780	endif
781
782
783
784	!Alberto. From this point if some PBL scheme has a non local term
785	! (not dependent on the local values of the variable)
786	! this can be added to b_t, b_q, or b_u, b_v.
787	!!!!!!!!!!!!!!!!!!!
788
789	#endif
790
791	ENDDO
792	!$OMP END PARALLEL DO
793	!
794	!$OMP PARALLEL DO &
795	!$OMP PRIVATE ( ij, i,j,k, its, ite, jts, jte )
796	DO ij = 1 , num_tiles
797
798	its = i_start(ij)
799	ite = i_end(ij)
800	jts = j_start(ij)
801	jte = j_end(ij)
802
803	pbl_select: SELECT CASE(bl_pbl_physics)
804
805	CASE (TEMFPBLSCHEME)
806	! WA Test
807	! WRITE( message , * ) 'Calling TEMF PBL scheme, timestep = ', itimestep
808	! CALL wrf_message ( message )
809	! print *,'Calling TEMF PBL scheme'
810	CALL wrf_debug(100,'in TEMF PBL')
811	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
812	PRESENT( qi_curr ) .AND. &
813	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
814	PRESENT( rqiblten ) .AND. &
815	PRESENT( te_temf ) .AND. PRESENT( cfm_temf ) .AND. &
816	PRESENT( hol ) ) THEN
817	CALL temfpbl( &
818	U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
819	,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
820	,P3D=p_phy,P3DI=p8w,PI3D=pi_phy,RHO=rho &
821	,RUBLTEN=rublten,RVBLTEN=rvblten &
822	,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
823	,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
824	,FLAG_QI=flag_qi &
825	,g=g,cp=cp,rcp=rcp,r_d=r_d,r_v=r_v,cpv=cpv &
826	,Z=z,XLV=XLV,PSFC=PSFC &
827	,MUT=mut,P_TOP=p_top &
828	,ZNT=znt,HT=ht,UST=ust,ZOL=zol,HOL=hol,HPBL=pblh &
829	,PSIM=psim,PSIH=psih,XLAND=xland &
830	,HFX=hfx,QFX=qfx,TSK=tskold,QSFC=qsfc,GZ1OZ0=gz1oz0 &
831	,U10=u10,V10=v10,T2=t2 &
832	,WSPD=wspd,BR=br,DT=dtbl,DTMIN=dtmin,KPBL2D=kpbl &
833	,SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0 &
834	,EP1=ep_1,EP2=ep_2,KARMAN=karman,EOMEG=eomeg &
835	,STBOLT=stbolt &
836	,te_temf=te_temf,kh_temf=kh_temf,km_temf=km_temf &
837	,shf_temf=shf_temf,qf_temf=qf_temf &
838	,uw_temf=uw_temf,vw_temf=vw_temf &
839	,hd_temf=hd_temf,lcl_temf=lcl_temf,hct_temf=hct_temf &
840	,wupd_temf=wupd_temf,mf_temf=mf_temf &
841	,thup_temf=thup_temf,qtup_temf=qtup_temf,qlup_temf=qlup_temf &
842	,cf3d_temf=cf3d_temf,cfm_temf=cfm_temf &
843	,flhc=flhc,flqc=flqc,exch_temf=exch_temf &
844	,fCor=f &
845	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
846	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
847	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
848	)
849	ELSE
850	CALL wrf_error_fatal('Lack arguments to call TEMF pbl')
851	ENDIF
852
853	CASE (YSUSCHEME)
854	CALL wrf_debug(100,'in YSU PBL')
855	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
856	PRESENT( qi_curr ) .AND. &
857	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
858	PRESENT( rqiblten ) .AND. &
859	PRESENT( hol ) ) THEN
860	CALL ysu( &
861	U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
862	,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
863	,P3D=p_phy,P3DI=p8w,PI3D=pi_phy &
864	,RUBLTEN=rublten,RVBLTEN=rvblten &
865	,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
866	,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
867	,FLAG_QI=flag_qi &
868	,CP=cp,G=g,ROVCP=rcp,RD=r_D,ROVG=rovg &
869	,DZ8W=dz8w,XLV=XLV,RV=r_v,PSFC=PSFC &
870	,ZNU=znu,ZNW=znw,MUT=mut,P_TOP=p_top &
871	,ZNT=znt,UST=ust,HPBL=pblh &
872	,PSIM=psim,PSIH=psih,XLAND=xland &
873	,HFX=hfx,QFX=qfx,GZ1OZ0=gz1oz0 &
874	,U10=u10,V10=v10 &
875	,WSPD=wspd,BR=br,DT=dtbl,KPBL2D=kpbl &
876	,EP1=ep_1,EP2=ep_2,KARMAN=karman &
877	,EXCH_H=exch_h,REGIME=regime &
878	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
879	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
880	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
881	)
882	ELSE
883	WRITE ( message , FMT = '(A,7(L1,1X))' ) &
884	'present: '// &
885	'qv_curr, '// &
886	'qc_curr, '// &
887	'qi_curr, '// &
888	'rqvblten, '// &
889	'rqcblten, '// &
890	'rqiblten, '// &
891	'hol = ' , &
892	PRESENT( qv_curr ) , &
893	PRESENT( qc_curr ) , &
894	PRESENT( qi_curr ) , &
895	PRESENT( rqvblten ) , &
896	PRESENT( rqcblten ) , &
897	PRESENT( rqiblten ) , &
898	PRESENT( hol )
899	CALL wrf_debug(0,message)
900	CALL wrf_error_fatal('Lack arguments to call YSU pbl')
901	ENDIF
902
903	CASE (MRFSCHEME)
904	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
905	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
906	PRESENT( hol ) .AND. &
907	.TRUE. ) THEN
908
909	CALL wrf_debug(100,'in MRF')
910	CALL mrf( &
911	U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
912	,QV3D=qv_curr &
913	,QC3D=qc_curr &
914	,QI3D=qi_curr &
915	,P3D=p_phy,PI3D=pi_phy &
916	,RUBLTEN=rublten,RVBLTEN=rvblten &
917	,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
918	,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
919	,CP=cp,G=g,ROVCP=rcp,R=r_d,ROVG=rovg &
920	,DZ8W=dz8w,Z=z,XLV=xlv,RV=r_v,PSFC=psfc &
921	,P1000MB=p1000mb &
922	,ZNT=znt,UST=ust,ZOL=zol,HOL=hol &
923	,PBL=pblh,PSIM=psim,PSIH=psih &
924	,XLAND=xland,HFX=hfx,QFX=qfx,TSK=tskold &
925	,GZ1OZ0=gz1oz0,WSPD=wspd,BR=br &
926	,DT=dtbl,DTMIN=dtmin,KPBL2D=kpbl &
927	,SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0 &
928	,EP1=ep_1,EP2=ep_2,KARMAN=karman,EOMEG=eomeg &
929	,STBOLT=stbolt,REGIME=regime &
930	,FLAG_QI=flag_qi &
931	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
932	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
933	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
934	)
935	ELSE
936	WRITE ( message , FMT = '(A,5(L1,1X))' ) &
937	'present: '// &
938	'qv_curr, '// &
939	'qc_curr, '// &
940	'rqvblten, '// &
941	'rqcblten, '// &
942	'hol = ' , &
943	PRESENT( qv_curr ) , &
944	PRESENT( qc_curr ) , &
945	PRESENT( rqvblten ) , &
946	PRESENT( rqcblten ) , &
947	PRESENT( hol )
948	CALL wrf_debug(0,message)
949	CALL wrf_error_fatal('Lack arguments to call MRF pbl')
950	ENDIF
951
952	CASE (GFSSCHEME)
953	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
954	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
955	.TRUE. ) THEN
956	CALL wrf_debug(100,'in GFS')
957	CALL bl_gfs( &
958	U3D=u_phytmp,V3D=v_phytmp &
959	,TH3D=th_phy,T3D=t_phy &
960	,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
961	,P3D=p_phy,PI3D=pi_phy &
962	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
963	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
964	,RQIBLTEN=rqiblten &
965	,CP=cp,G=g,ROVCP=rcp,R=r_d,ROVG=rovg &
966	,DZ8W=dz8w,z=z,PSFC=psfc &
967	,UST=ust,PBL=pblh,PSIM=psim,PSIH=psih &
968	,HFX=hfx,QFX=qfx,TSK=tskold,GZ1OZ0=gz1oz0 &
969	,WSPD=wspd,BR=br &
970	,DT=dtbl,KPBL2D=kpbl,EP1=ep_1,KARMAN=karman &
971	#if (NMM_CORE==1)
972	,DISHEAT=DISHEAT &
973	#endif
974	,P_QI=p_qi,P_FIRST_SCALAR=param_first_scalar &
975	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
976	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
977	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
978	)
979	ELSE
980	WRITE ( message , FMT = '(A,4(L1,1X))' ) &
981	'present: '// &
982	'qv_curr, '// &
983	'qc_curr, '// &
984	'rqvblten, '// &
985	'rqcblten = ', &
986	PRESENT( qv_curr ) , &
987	PRESENT( qc_curr ) , &
988	PRESENT( rqvblten ) , &
989	PRESENT( rqcblten )
990	CALL wrf_debug(0,message)
991	CALL wrf_error_fatal('Lack arguments to call GFS pbl')
992	ENDIF
993
994	CASE (MYJPBLSCHEME)
995	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
996	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
997	.TRUE. ) THEN
998
999	CALL wrf_debug(100,'in MYJPBL')
1000	IF ( .not.flag_bep .and. idiff.ne.1) THEN
1001	CALL myjpbl(DT=dt,STEPBL=stepbl,HT=ht,DZ=dz8w &
1002	,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1003	,QV=qv_curr,QCW=qc_curr,QCI=qi_curr,QCS=qs_curr & !BSF
1004	,QCR=qr_curr,QCG=qg_curr & !BSF
1005	,U=u_phy,V=v_phy,RHO=rho &
1006	,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1007	,QZ0=qz0,UZ0=uz0,VZ0=vz0 &
1008	,LOWLYR=lowlyr &
1009	,XLAND=xland,SICE=xice,SNOW=snow &
1010	,TKE_MYJ=tke_pbl,EXCH_H=exch_h,USTAR=ust,ZNT=znt &
1011	,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1012	,AKHS=akhs,AKMS=akms,ELFLX=lh,MIXHT=mixht &
1013	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1014	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1015	,RQIBLTEN=rqiblten,RQSBLTEN=rqsblten & !BSF
1016	,RQRBLTEN=rqrblten,RQGBLTEN=rqgblten & !BSF
1017	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1018	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1019	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1020	)
1021	ELSE !(SF_URBAN_PHYSICS.EQ.2)
1022	! Bep changes begin
1023
1024	CALL myjurb(IDIFF=idiff,FLAG_BEP=flag_bep &
1025	,DT=dtbl,STEPBL=stepbl,HT=ht,DZ=dz8w &
1026	,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1027	,QV=qv_curr, CWM=qc_curr &
1028	,U=u_phy,V=v_phy,RHO=rho &
1029	,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1030	,QZ0=qz0,UZ0=uz0,VZ0=vz0 &
1031	,LOWLYR=lowlyr &
1032	,XLAND=xland,SICE=xice,SNOW=snow &
1033	,TKE_MYJ=tke_pbl,EXCH_H=exch_h,EXCH_M=exch_m &
1034	,USTAR=ust,ZNT=znt &
1035	,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1036	,AKHS=akhs,AKMS=akms,ELFLX=lh &
1037	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1038	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1039	! Multi-layer UCM
1040	,FRC_URB2D=frc_urb2d &
1041	,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1042	,A_Q_BEP=a_q_bep &
1043	,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1044	,B_T_BEP=b_t_bep,B_Q_BEP=b_q_bep &
1045	,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep &
1046	,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep &
1047	! UCM end
1048	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1049	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1050	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1051	)
1052	ENDIF
1053
1054	ELSE
1055	WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1056	'present: '// &
1057	'qv_curr, '// &
1058	'qc_curr, '// &
1059	'rqvblten, '// &
1060	'rqcblten = ' , &
1061	PRESENT( qv_curr ) , &
1062	PRESENT( qc_curr ) , &
1063	PRESENT( rqvblten ) , &
1064	PRESENT( rqcblten )
1065	CALL wrf_debug(0,message)
1066	CALL wrf_error_fatal('Lack arguments to call MYJ pbl')
1067	ENDIF
1068
1069	CASE (QNSEPBLSCHEME)
1070	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1071	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1072	.TRUE. ) THEN
1073	CALL wrf_debug(100,'in QNSEPBL')
1074	CALL qnsepbl( &
1075	DT=dt,STEPBL=stepbl,HT=ht,DZ=dz8w &
1076	,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1077	,QV=qv_curr, CWM=qc_curr &
1078	,U=u_phy,V=v_phy,RHO=rho &
1079	,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1080	,QZ0=qz0,UZ0=uz0,VZ0=vz0,CORF=f &
1081	,LOWLYR=lowlyr &
1082	,XLAND=xland,SICE=xice,SNOW=snow &
1083	,TKE=tke_pbl,EXCH_H=exch_h,EXCH_M=exch_m,USTAR=ust,ZNT=znt &
1084	,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1085	,AKHS=akhs,AKMS=akms,ELFLX=lh &
1086	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1087	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1088	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1089	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1090	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1091	)
1092	ELSE
1093	WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1094	'present: '// &
1095	'qv_curr, '// &
1096	'qc_curr, '// &
1097	'rqvblten, '// &
1098	'rqcblten = ' , &
1099	PRESENT( qv_curr ) , &
1100	PRESENT( qc_curr ) , &
1101	PRESENT( rqvblten ) , &
1102	PRESENT( rqcblten )
1103	CALL wrf_debug(0,message)
1104	CALL wrf_error_fatal('Lack arguments to call QNSE pbl')
1105	ENDIF
1106
1107	CASE (ACMPBLSCHEME)
1108
1109	!! These are values that are not supplied to pbl driver, but are required by ACM
1110	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1111	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1112	.TRUE. ) THEN
1113	CALL wrf_debug(100,'in ACM PBL')
1114
1115	CALL ACMPBL( &
1116	XTIME=itimestep, DTPBL=dtbl, ZNW=znw, SIGMAH=znu &
1117	,U3D=u_phytmp, V3D=v_phytmp, PP3D=p_phy, DZ8W=dz8w, TH3D=th_phy, T3D=t_phy &
1118	,QV3D=qv_curr, QC3D=qc_curr, QI3D=qi_curr, RR3D=rho &
1119	,UST=UST, HFX=HFX, QFX=QFX, TSK=tsk &
1120	,PSFC=PSFC, EP1=EP_1, G=g, ROVCP=rcp,RD=r_D,CPD=cp &
1121	,PBLH=pblh, KPBL2D=kpbl, REGIME=regime &
1122	,GZ1OZ0=gz1oz0,WSPD=wspd,PSIM=psim, MUT=mut &
1123	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1124	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1125	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1126	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1127	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1128	)
1129	ELSE
1130	WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1131	'present: '// &
1132	'qv_curr, '// &
1133	'qc_curr, '// &
1134	'rqvblten, '// &
1135	'rqcblten = ' , &
1136	PRESENT( qv_curr ) , &
1137	PRESENT( qc_curr ) , &
1138	PRESENT( rqvblten ) , &
1139	PRESENT( rqcblten )
1140	CALL wrf_debug(0,message)
1141	CALL wrf_error_fatal('Lack arguments to call ACM2 pbl')
1142	ENDIF
1143
1144	#if (EM_CORE==1)
1145
1146	CASE (MYNNPBLSCHEME2, MYNNPBLSCHEME3)
1147
1148	IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1149	PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1150	& PRESENT(qke) .AND. PRESENT(tsq) .AND. &
1151	& PRESENT(qsq) .AND. PRESENT(cov) .AND. &
1152	& PRESENT(rmol) .AND. &
1153	& PRESENT(qcg) .AND. PRESENT(ch) .AND. &
1154	& PRESENT(grav_settling) ) THEN
1155
1156	CALL wrf_debug(100,'in MYNNPBL')
1157
1158	IF (itimestep==1) THEN
1159	initflag=1
1160	ELSE
1161	initflag=0
1162	ENDIF
1163
1164	CALL mynn_bl_driver(&
1165	&initflag=initflag,&
1166	&grav_settling=grav_settling,&
1167	&delt=dtbl,&
1168	&dz=dz8w,&
1169	&u=u_phy,v=v_phy,th=th_phy,qv=qv_curr,qc=qc_curr,&
1170	&p=p_phy,exner=pi_phy,rho=rho,&
1171	&xland=xland,ts=tsk,qsfc=qsfc,qcg=qcg,ps=psfc,&
1172	&ust=ust,ch=ch,hfx=hfx,qfx=qfx,rmol=rmol,wspd=wspd,&
1173	&Qke=qke,Tsq=tsq,Qsq=qsq,Cov=cov,&
1174	&Du=rublten,Dv=rvblten,Dth=rthblten,&
1175	&Dqv=rqvblten,Dqc=rqcblten,&
1176	&k_h=exch_h,k_m=exch_m,&
1177	&pblh=pblh&
1178	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1179	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1180	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1181	)
1182
1183	IF (windspec.GT.0 &
1184	.AND.PRESENT(id) &
1185	.AND.PRESENT(phb) &
1186	.AND.PRESENT(xlat_u).AND.PRESENT(xlong_u) &
1187	.AND.PRESENT(xlat_v).AND.PRESENT(xlong_v)) THEN
1188	CALL turbine_drag( &
1189	& ID=id &
1190	&,PHB=phb,u=u_phy,v=v_phy &
1191	&,XLAT_U=xlat_u,XLONG_U=xlong_u &
1192	&,XLAT_V=xlat_v,XLONG_V=xlong_v &
1193	&,DX=dx,DZ=dz8w,DT=dt &
1194	&,QKE=qke,DU=rublten,DV=rvblten &
1195	&,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1196	&,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1197	&,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1198	&)
1199	ENDIF
1200
1201	ELSE
1202	WRITE ( message , FMT = '(A,12(L1,1X))' ) &
1203	'present: '// &
1204	'qv_curr, '// &
1205	'qc_curr, '// &
1206	'rqvblten, '// &
1207	'rqcblten, '// &
1208	'qke, '// &
1209	'tsq = '// &
1210	'qsq = '// &
1211	'cov = '// &
1212	'rmol = '// &
1213	'qcg = '// &
1214	'ch = '// &
1215	'grav_settling = ', &
1216	PRESENT( qv_curr ) , &
1217	PRESENT( qc_curr ) , &
1218	PRESENT( rqvblten ) , &
1219	PRESENT( rqcblten ) , &
1220	PRESENT( qke ) , &
1221	PRESENT( tsq ) , &
1222	PRESENT( qsq ) , &
1223	PRESENT( cov ) , &
1224	PRESENT( rmol ) , &
1225	PRESENT( qcg ) , &
1226	PRESENT( ch ) , &
1227	PRESENT( grav_settling)
1228	CALL wrf_debug(0,message)
1229	CALL wrf_error_fatal('Lack arguments to call MYNN pbl')
1230	ENDIF
1231
1232	CASE (BOULACSCHEME)
1233
1234	CALL wrf_debug(100,'in boulac')
1235
1236	CALL BOULAC(FRC_URB2D=frc_urb2d,IDIFF=idiff,FLAG_BEP=flag_bep &
1237	,DZ8W=dz8w,DT=dt,U_PHY=u_phytmp &
1238	,V_PHY=v_phytmp,TH_PHY=th_phy &
1239	,RHO=rho,QV_CURR=qv_curr,HFX=hfx &
1240	,QFX=qfx,USTAR=ust,CP=cp,G=g &
1241	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1242	,RQVBLTEN=rqvblten,RQCBLTEN=rqvblten &
1243	,TKE=tke_pbl,DLK=el_pbl,WU=wu_tur,WV=wv_tur,WT=wt_tur,WQ=wq_tur &
1244	,EXCH_H=exch_h,EXCH_M=exch_m,PBLH=pblh &
1245	,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1246	,A_Q_BEP=a_q_bep &
1247	,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1248	,B_T_BEP=b_t_bep,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep &
1249	,B_Q_BEP=b_q_bep &
1250	,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep &
1251	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1252	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1253	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
1254
1255	CASE (CAMUWPBLSCHEME)
1256
1257	IF ( PRESENT( z_at_w ) .AND. PRESENT( tauresx2d )) THEN
1258	CALL wrf_debug(100,'in camuwpbl')
1259
1260	CALL CAMUWPBL(DT=dt,U_PHY=u_phy,V_PHY=v_phy,TH_PHY=th_phy,RHO=rho &
1261	,QV_CURR=qv_curr,HFX=hfx,QFX=qfx,USTAR=ust,RUBLTEN=rublten &
1262	,RVBLTEN=rvblten,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1263	,RQCBLTEN=rqcblten,TKE_PBL=tke_pbl,PBLH2D=pblh,KPBL2D=kpbl &
1264	,P8W=p8w,P_PHY=p_phy,Z=z,T_PHY=t_phy,QC_CURR=qc_curr &
1265	,QI_CURR=qi_curr,Z_AT_W=z_at_w,CLDFRA=cldfra,HT=ht &
1266	,RTHRATENLW=rthratenlw,EXNER=pi_phy,ITIMESTEP=itimestep &
1267	,TAURESX2D=tauresx2d,TAURESY2D=tauresy2d,KVH3D=exch_h,KVM3D=exch_m &
1268	,TPERT2D=tpert2d,QPERT2D=qpert2d,WPERT2D=wpert2d &
1269	,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1270	,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1271	,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte )
1272
1273	ELSE
1274	CALL wrf_debug(0,message)
1275	CALL wrf_error_fatal('Lack arguments to call CAMUWPBL pbl')
1276	ENDIF
1277	#endif
1278
1279
1280	CASE DEFAULT
1281
1282	WRITE( message , * ) 'The pbl option does not exist: bl_pbl_physics = ', bl_pbl_physics
1283	CALL wrf_error_fatal ( message )
1284
1285	END SELECT pbl_select
1286
1287	IF (PRESENT(dtaux3d)) THEN
1288	IF(gwd_opt .EQ. 1)THEN
1289	CALL gwdo( &
1290	U3D=u_phytmp,V3D=v_phytmp,T3D=t_phy &
1291	,QV3D=qv_curr &
1292	,P3D=p_phy,P3DI=p8w,PI3D=pi_phy,Z=z &
1293	,RUBLTEN=rublten,RVBLTEN=rvblten &
1294	,DTAUX3D=dtaux3d,DTAUY3D=dtauy3d &
1295	,DUSFCG=dusfcg,DVSFCG=dvsfcg &
1296	,VAR2D=var2d,OC12D=oc12d &
1297	,OA2D1=oa1,OA2D2=oa2,OA2D3=oa3,OA2D4=oa4 &
1298	,OL2D1=ol1,OL2D2=ol2,OL2D3=ol3,OL2D4=ol4 &
1299	,ZNU=znu,ZNW=znw,MUT=mut,P_TOP=p_top &
1300	,CP=cp,G=g,RD=r_d &
1301	,RV=r_v,EP1=ep_1,PI=3.141592653 &
1302	,DT=dtbl,DX=dx,KPBL2D=kpbl,ITIMESTEP=itimestep &
1303	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1304	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1305	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
1306	ENDIF
1307	ENDIF
1308
1309
1310	#if (NMM_CORE != 1)
1311
1312	IF (idiff.eq.1) THEN
1313
1314	!Alberto: here we call the general routine to solve the diffusion
1315	! + all the source/sink terms.
1316	! the only thing that should be passed from the PBL schemes is the value of the exch_h, and exch_m
1317	! (and the non local terms, if present, through the b).
1318	! So, this routine can be used by any of the previous schemes. We only need to pass the right variables
1319	! (and eliminate the computation of the tendencies from the previous routines, to avoid doing twice the job).
1320	! As I explain below, in the routine, here we could extract the vertical turbulent
1321	! fluxes (something that will be general for all the routines).
1322	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1323
1324
1325	CALL diff3d (DT=dtbl,CP=cp,DZ=dz8w,TH=th_phy,QV=qv_curr,QC=qc_curr,T=t_phy &
1326	,U=u_phy,V=v_phy,RHO=rho,EXCH_H=exch_h &
1327	,EXCH_M=exch_m &
1328	,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1329	,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1330	,WU=wu_tur,WV=wv_tur,WT=wt_tur,WQ=wq_tur &
1331	,A_U=a_u,A_V=a_v,A_T=a_t,A_Q=a_q &
1332	,B_U=b_u,B_V=b_v,B_T=b_t,B_Q=b_q &
1333	,SF=sf,VL=vl &
1334	,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1335	,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1336	,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte)
1337
1338	DEALLOCATE (a_u) ! Implicit component for the momemtum in X-direction
1339	DEALLOCATE (a_v) ! Implicit component for the momemtum in Y-direction
1340	DEALLOCATE (a_t) ! Implicit component for the Pot. Temp.
1341	DEALLOCATE (a_q) ! Implicit component for the water vapor
1342	DEALLOCATE (b_u) ! Explicit component for the momemtum in X-direction
1343	DEALLOCATE (b_v) ! Explicit component for the momemtum in Y-direction
1344	DEALLOCATE (b_t) ! Explicit component for the Pot. Temp.
1345	DEALLOCATE (b_q) ! Explicit component for the water vapor
1346	DEALLOCATE (sf ) ! surfaces
1347	DEALLOCATE (vl ) ! volumes
1348
1349	ENDIF !idiff
1350	#endif
1351
1352	ENDDO
1353	!$OMP END PARALLEL DO
1354
1355	ENDIF
1356	!
1357
1358
1359
1360	END SUBROUTINE pbl_driver
1361
1362	!=============================================================================
1363	SUBROUTINE diff3d(DT,CP,DZ,TH ,QV,QC,T,U,V,RHO &
1364	,EXCH_H,EXCH_M &
1365	,RUBLTEN,RVBLTEN,RTHBLTEN &
1366	,RQVBLTEN,RQCBLTEN &
1367	,WU,WV,WT,WQ &
1368	,A_U,A_V,A_T,A_Q &
1369	,B_U,B_V,B_T,B_Q &
1370	,SF,VL &
1371	,IDS,IDE,JDS,JDE,KDS,KDE &
1372	,IMS,IME,JMS,JME,KMS,KME &
1373	,ITS,ITE,JTS,JTE,KTS,KTE &
1374	)
1375	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1376	! Subroutine written by Alberto Martilli, CIEMAT, Spain,
1377	! e-mail:alberto_martilli@ciemat.es
1378	! August 2008.
1379	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1380	! ALberto
1381	! This routine solves the vertical diffusion
1382	! + other source/sink terms (surface fluxes, building drag, non local terms, etc.)
1383	! for U,V, potential temperature and moisture. The equation should be written
1384	! as follow, for a generic variable M:
1385	!
1386	! dM 1 d dM
1387	! ---- =---- ---(rho*K----)+AM+B
1388	! dt rho dz dz
1389	! where A and B are the implict and explcit coefficients of the source/sink terms
1390	! (at the lowest level the surface fluxes should go in these terms)
1391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392	!-----------------------------------------------------------------------
1393
1394	IMPLICIT NONE
1395	!
1396	!----------------------------------------------------------------------
1397	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1398	& ,IMS,IME,JMS,JME,KMS,KME &
1399	& ,ITS,ITE,JTS,JTE,KTS,KTE
1400
1401	! INputs
1402	real DT,CP
1403	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ ! Depth of vertical levels
1404	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: TH ! Potential Temperature
1405	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QV ! water vapor mixing ratio (kg/kg)
1406	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QC ! cloud water mixing ratio (kg/kg)
1407	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: T ! temperature
1408	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: U ! X-compnent of the wind
1409	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: V ! Y-compnent of the wind
1410	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: RHO ! Air Density
1411	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_H ! Turbulent coefficient for heat
1412	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_M ! Turbulent coefficient for momentum
1413
1414	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_U ! Implicit coefficient for DRAG (x -component)
1415	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_U ! Explicit coefficient for DRAG (x -component)
1416	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_V ! Implicit coefficient for DRAG (y -component)
1417	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_V ! Explicit coefficient for DRAG (y -component)
1418	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_T ! Implicit coefficient for Heat flux from buildings
1419	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_T ! Explicit coefficient for Heat flux from buildings
1420	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_Q ! Implicit coefficient for moisture flux
1421	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_Q ! Explicit coefficient for moisture flux
1422
1423	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: VL ! fraction of air volume in the grid cell
1424	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: SF ! fraction of air at the face of the grid cell
1425	!outputs
1426	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RUBLTEN ! tendency for x-wind component
1427	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RVBLTEN ! tendency for y-wind component
1428	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RTHBLTEN ! tendency for Potential Temperature
1429	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RQVBLTEN ! tendency for water vapor mixing ratio (kg/kg)
1430	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RQCBLTEN ! tendency for cloud water mixing ratio (kg/kg)
1431	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WU ! turbulent flux of momentum (x)
1432	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WV ! turbulent flux of momentum (y)
1433	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WT ! turbulent flux of potential temperature
1434	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WQ ! turbulent flux of water vapor
1435	! Parameters
1436	REAL ELOCP
1437	! locals (same meaning as above, but 1D)
1438
1439	real u1d(kms:kme),v1d(kms:kme),exch_h1d(kms:kme)
1440	real the1d(kms:kme) ! Equivalent potential temperature
1441	real exch_m1d(kms:kme),qv1d(kms:kme),qc1d(kms:kme)
1442	real dz1d(kms:kme),rho1d(kms:kme),rhoz1d(kms:kme)
1443	real sf1d(kms:kme),vl1d(kms:kme)
1444	real a_u1d(kms:kme),b_u1d(kms:kme)
1445	real a_v1d(kms:kme),b_v1d(kms:kme)
1446	real a_t1d(kms:kme),b_t1d(kms:kme)
1447	real a_q1d(kms:kme),b_q1d(kms:kme)
1448	real a_qc1d(kms:kme),b_qc1d(kms:kme)
1449	real wu1d(kms:kme),wv1d(kms:kme),wt1d(kms:kme),wq1d(kms:kme),wqc1d(kms:kme)
1450	real thnew
1451	!
1452	integer i,k,j
1453	!----------------------------------------------------------------------
1454	ELOCP=2.72E6/CP
1455	u1d=0.
1456	v1d=0.
1457	exch_h1d=0.
1458	exch_m1d=0.
1459	qv1d=0.
1460	qc1d=0.
1461	dz1d=0.
1462	rho1d=0.
1463	rhoz1d=0.
1464	sf1d=0.
1465	vl1d=0.
1466	a_u1d=0.
1467	a_v1d=0.
1468	a_t1d=0.
1469	a_q1d=0.
1470	a_qc1d=0.
1471	b_u1d=0.
1472	b_v1d=0.
1473	b_t1d=0.
1474	b_q1d=0.
1475	b_qc1d=0.
1476
1477	do j=jts,jte
1478	do i=its,ite
1479
1480	! put three D variables in temporary 1D variables
1481
1482	do k=kts,kte
1483	u1d(k)=U(i,k,j)
1484	v1d(k)=V(i,k,j)
1485	the1d(k)=TH(i,k,j)(QC(i,k,j)(-ELOCP/T(i,k,j))+1)
1486	qv1d(k)=qv(i,k,j)
1487	dz1d(k)=dz(i,k,j)
1488	rho1d(k)=rho(i,k,j)
1489	a_u1d(k)=a_u(i,k,j)
1490	b_u1d(k)=b_u(i,k,j)
1491	a_v1d(k)=a_v(i,k,j)
1492	b_v1d(k)=b_v(i,k,j)
1493	a_t1d(k)=a_t(i,k,j)
1494	b_t1d(k)=b_t(i,k,j)
1495	a_q1d(k)=a_q(i,k,j)
1496	b_q1d(k)=b_q(i,k,j)
1497	a_qc1d(k)=0.
1498	b_qc1d(k)=0.
1499	vl1d(k)=vl(i,k,j)
1500	sf1d(k)=sf(i,k,j)
1501	enddo
1502	sf1d(kte+1)=1.
1503	do k=kts,kte
1504	exch_h1d(k)=exch_h(i,k,j)
1505	exch_m1d(k)=exch_m(i,k,j)
1506	enddo
1507	exch_h1d(kts)=0.
1508	! exch_h1d(kte+1)=0
1509	exch_m1d(kts)=0.
1510	! exch_m1d(kte+1)=0
1511	rhoz1d(kts)=rho1d(kts)
1512	do k=kts+1,kte
1513	rhoz1d(k)=(rho1d(k)dz1d(k-1)+rho1d(k-1)dz1d(k))/ &
1514	& (dz1d(k-1)+dz1d(k))
1515	enddo
1516	rhoz1d(kte+1)=rho1d(kte)
1517
1518
1519	! solve the diffusion for x-component of the wind
1520
1521	call diff(kms,kme,kts,kte,dt,u1d,rho1d,rhoz1d,exch_m1d,a_u1d,b_u1d,sf1d, &
1522	& vl1d,dz1d,wu1d)
1523
1524	! solve the diffusion for y-component of the wind
1525
1526	call diff(kms,kme,kts,kte,dt,v1d,rho1d,rhoz1d,exch_m1d,a_v1d,b_v1d,sf1d, &
1527	& vl1d,dz1d,wv1d)
1528
1529	! solve the diffusion for equivalent potential temperature
1530
1531	call diff(kms,kme,kts,kte,dt,the1d,rho1d,rhoz1d,exch_h1d,a_t1d,b_t1d,sf1d, &
1532	& vl1d,dz1d,wt1d)
1533
1534	! solve the diffusion for the water vapor mixing ratio
1535
1536	call diff(kms,kme,kts,kte,dt,qv1d,rho1d,rhoz1d,exch_h1d,a_q1d,b_q1d,sf1d, &
1537	& vl1d,dz1d,wq1d)
1538
1539	! solve the diffusion for the cloud water mixing ratio
1540
1541	call diff(kms,kme,kts,kte,dt,qc1d,rho1d,rhoz1d,exch_h1d,a_qc1d,b_qc1d,sf1d, &
1542	& vl1d,dz1d,wqc1d)
1543
1544	!
1545	! compute the tendencies
1546
1547	do k=kts,kte
1548	rublten(i,k,j)=(u1d(k)-u(i,k,j))/dt
1549	rvblten(i,k,j)=(v1d(k)-v(i,k,j))/dt
1550	thnew=the1d(k)/(QC(i,k,j)*(-ELOCP/T(i,k,j))+1)
1551	rthblten(i,k,j)=(thnew-th(i,k,j))/dt
1552	rqvblten(i,k,j)=(qv1d(k)-qv(i,k,j))/dt
1553	rqcblten(i,k,j)=(qc1d(k)-qc(i,k,j))/dt
1554	wu(i,k,j)=wu1d(k)
1555	wv(i,k,j)=wv1d(k)
1556	wt(i,k,j)=wt1d(k)
1557	wq(i,k,j)=wq1d(k)
1558	enddo
1559	enddo
1560	enddo
1561	!!!!!!!!!!!!
1562
1563
1564	END SUBROUTINE diff3d
1565	! ===6=8===============================================================72
1566	! ===6=8===============================================================72
1567
1568	subroutine diff(kms,kme,kts,kte,dt,co,da,daz,cd,aa,bb,sf,vl,dz,fc)
1569
1570	!------------------------------------------------------------------------
1571	! Calculation of the diffusion in 1D
1572	!------------------------------------------------------------------------
1573	! - Input:
1574	! nz : number of points
1575	! iz1 : first calculated point
1576	! co : concentration of the variable of interest
1577	! dz : vertical levels
1578	! cd : diffusion coefficients
1579	! da : density of the air at the center
1580	! daz : density of the air at the face
1581	! dt : diffusion time step
1582	! - Output:
1583	! co :concentration of the variable of interest
1584
1585	! - Internal:
1586	! cddz : constant terms in the equations
1587	!---------------------------------------------------------------------
1588
1589	implicit none
1590	integer iz,iz1,izf
1591	integer kms,kme,kts,kte
1592	real dt,dzv
1593	real co(kms:kme),cd(kms:kme),dz(kms:kme)
1594	real da(kms:kme),daz(kms:kme)
1595	real cddz(kms:kme),fc(kms:kme),df(kms:kme)
1596	real a(kms:kme,3),c(kms:kme)
1597	real sf(kms:kme),vl(kms:kme)
1598	real aa(kms:kme),bb(kms:kme)
1599
1600
1601	! Compute cddz=2*cd/dz
1602
1603	cddz(kts)=sf(kts)daz(kts)cd(kts)/dz(kts)
1604	do iz=kts+1,kte
1605	cddz(iz)=2.sf(iz)daz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
1606	enddo
1607	cddz(kte+1)=sf(kte+1)daz(kte+1)cd(kte+1)/dz(kte)
1608
1609	iz1=1
1610	izf=1
1611
1612	do iz=iz1,kte-1
1613
1614	dzv=vl(iz)*dz(iz)
1615	a(iz,1)=-cddz(iz)*dt/dzv/da(iz)
1616	a(iz,2)=1+dt(cddz(iz)+cddz(iz+1))/dzv/da(iz)-aa(iz)dt
1617	a(iz,3)=-cddz(iz+1)*dt/dzv/da(iz)
1618	c(iz)=co(iz)+bb(iz)*dt
1619	enddo
1620
1621	do iz=kte-(izf-1),kte
1622	a(iz,1)=0.
1623	a(iz,2)=1
1624	a(iz,3)=0.
1625	c(iz)=co(iz)
1626	enddo
1627	call invert (kms,kme,kts,kte,a,c,co)
1628
1629	!let compute the fluxes, as diagnostic variable
1630
1631	do iz=kts,iz1
1632	fc(iz)=0.
1633	enddo
1634
1635	do iz=iz1+1,kte
1636	fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/da(iz)
1637	enddo
1638
1639
1640
1641	return
1642	end subroutine diff
1643	! ===6=8===============================================================72
1644
1645	subroutine invert(kms,kme,kts,kte,a,c,x)
1646
1647	!ccccccccccccccccccccccccccccccc
1648	! Aim: Inversion and resolution of a tridiagonal matrix
1649	! A X = C
1650	! Input:
1651	! a(*,1) lower diagonal (Ai,i-1)
1652	! a(*,2) principal diagonal (Ai,i)
1653	! a(*,3) upper diagonal (Ai,i+1)
1654	! c
1655	! Output
1656	! x results
1657	!ccccccccccccccccccccccccccccccc
1658
1659	implicit none
1660	integer kms,kme,kts,kte,in
1661	real a(kms:kme,3),c(kms:kme),x(kms:kme)
1662
1663	do in=kte-1,kts,-1
1664	c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
1665	a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
1666	enddo
1667
1668	do in=kts+1,kte
1669	c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
1670	enddo
1671
1672	do in=kts,kte
1673	x(in)=c(in)/a(in,2)
1674	enddo
1675
1676	return
1677	end subroutine invert
1678
1679	! ===6=8===============================================================72
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693	END MODULE module_pbl_driver

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

Download in other formats:

Original Format