Context Navigation

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

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