| blob | f70307176532c173665af1fb73aa37c00e67c4c1 |
1 !WRF:MEDIATION_LAYER:PHYSICS
2 !
4 MODULE module_pbl_driver
5 real, private, parameter :: QCLOUD_MIN = 0.0
6 real, private, parameter :: QVAPOR_MIN = 0.0
7 real, private, parameter :: QKE_MIN = 1.e-4
8 CONTAINS
10 !------------------------------------------------------------------
11 SUBROUTINE pbl_driver( &
12 itimestep,dt,u_frame,v_frame &
13 ,bldt,curr_secs,adapt_step_flag &
14 ,bldtacttime &
15 ,rublten,rvblten,rthblten &
16 ,tsk,xland,znt &
17 ,ht &
18 ,ust,pblh,hfx,qfx,grdflx,hpbl_hold &
19 ,u_phy,v_phy,w,th_phy,rho &
20 ,p_phy,pi_phy,p8w,t_phy,dz8w,z &
21 ,exch_h,exch_m,akhs,akms &
22 ,thz0,qz0,uz0,vz0,qsfc,f &
23 ,lowlyr,u10,v10,uoce,voce,t2 &
24 ,psim,psih,fm,fhh,gz1oz0, wspd,br,chklowq &
25 ,bl_pbl_physics, ra_lw_physics, dx, dy &
26 ,dx2d ,area2d &
27 ,stepbl,warm_rain &
28 ,kpbl,mixht,ct,lh,snow,xice &
29 ,znu, znw, mut, p_top &
30 ,ctopo,ctopo2,windfarm_opt,power &
31 ,windfarm_wake_model, windfarm_overlap_method &
32 ,ysu_topdown_pblmix &
33 ,shinhong_tke_diag &
34 ! OPTIONAL for TEMF scheme
35 ,te_temf,km_temf,kh_temf &
36 ,shf_temf,qf_temf,uw_temf,vw_temf &
37 ,hd_temf,lcl_temf,hct_temf &
38 ,wupd_temf,mf_temf,thup_temf,qtup_temf,qlup_temf &
39 ,exch_temf,cf3d_temf,cfm_temf &
40 ,flhc,flqc &
41 ! MYNN
42 ,qke,Sh3d,Sm3d &
43 ,qke_adv,bl_mynn_tkeadvect &
44 ,tsq,qsq,cov,rmol,ch,qcg,grav_settling &
45 ,dqke,qWT,qSHEAR,qBUOY,qDISS,tke_budget &
46 ,bl_mynn_closure,bl_mynn_cloudpdf &
47 ,bl_mynn_mixlength &
48 ,icloud_bl,qc_bl,qi_bl,cldfra_bl &
49 ,bl_mynn_edmf,bl_mynn_edmf_mom,bl_mynn_edmf_tke &
50 ,bl_mynn_mixscalars,bl_mynn_output &
51 ,bl_mynn_cloudmix,bl_mynn_mixqt &
52 ,edmf_a,edmf_w,edmf_thl &
53 ,edmf_qt,edmf_ent,edmf_qc &
54 ,sub_thl3D,sub_sqv3D &
55 ,det_thl3D,det_sqv3D &
56 ,vdfg &
57 ,maxwidth,maxMF,ztop_plume,ktop_plume &
58 ,spp_pbl,pattern_spp_pbl &
59 ! EEPS
60 ,pek,pep,pek_adv,pep_adv &
61 ,restart,cycling &
62 ,ids,ide, jds,jde, kds,kde &
63 ,ims,ime, jms,jme, kms,kme &
64 ,i_start,i_end, j_start,j_end, kts,kte, num_tiles &
65 ! Optional
66 ,hol, mol, regime &
67 ! Optional gravity-wave drag
68 ,gwd_opt,gwd_diags &
69 ,dtaux3d,dtauy3d &
70 ,dusfcg,dvsfcg,var2d,oc12d &
71 ,oa1,oa2,oa3,oa4,ol1,ol2,ol3,ol4 &
72 ,sina,cosa &
73 ,dtaux3d_ls,dtauy3d_ls &
74 ,dtaux3d_bl,dtauy3d_bl &
75 ,dtaux3d_ss,dtauy3d_ss &
76 ,dtaux3d_fd,dtauy3d_fd &
77 ,dusfcg_ls,dvsfcg_ls &
78 ,dusfcg_bl,dvsfcg_bl &
79 ,dusfcg_ss,dvsfcg_ss &
80 ,dusfcg_fd,dvsfcg_fd &
81 ,var2dls,oc12dls &
82 ,oa1ls,oa2ls,oa3ls,oa4ls,ol1ls,ol2ls,ol3ls,ol4ls &
83 ,var2dss,oc12dss &
84 ,oa1ss,oa2ss,oa3ss,oa4ss,ol1ss,ol2ss,ol3ss,ol4ss &
85 ! Optional moisture tracers
86 ,qv_curr, qc_curr, qr_curr &
87 ,qi_curr, qs_curr, qg_curr &
88 ,rqvblten,rqcblten,rqiblten &
89 ,rqrblten,rqsblten,rqgblten &
90 ! Optional aerosol
91 ,qnwfa_curr,f_qnwfa &
92 ,qnifa_curr,f_qnifa &
93 ,qnbca_curr,f_qnbca &
94 ! Optional moisture tracer flags
95 ,f_qv,f_qc,f_qr &
96 ,f_qi,f_qs,f_qg &
97 ! variables added for BEP
98 ,frc_urb2d &
99 ,a_u_bep,a_v_bep,a_t_bep,a_q_bep &
100 ,b_u_bep,b_v_bep,b_t_bep,b_q_bep &
101 ,sf_bep,vl_bep &
102 ,sf_sfclay_physics,sf_urban_physics &
103 ,tke_pbl,diss_pbl,tpe_pbl &
104 ,tke_adv,diss_adv,tpe_adv &
105 ,pr_pbl,el_pbl &
106 ,wu_tur,wv_tur,wt_tur,wq_tur &
107 ! variables added for AHE
108 , gmt, xtime, julday, julyr, ahe &
109 , distributed_ahe_opt &
110 ! variables for GBM PBL
111 ,exch_tke, rthraten &
112 ,a_e_bep,b_e_bep,dlg_bep,dl_u_bep &
113 ,mfshconv, massflux_EDKF, entr_EDKF, detr_EDKF &
114 ,thl_up, thv_up, rt_up ,rv_up, rc_up, u_up, v_up &
115 ,frac_up,rc_mf &
116 ! Wind Turbine Parameterizations
117 ,phb,xlat_u,xlong_u,xlat_v,xlong_v,id &
118 ! WRF-Solar EPS
119 ,multi_perturb, pert_mynn &
120 ,perts_qvapor &
121 ,perts_qcloud &
122 ,perts_th &
123 ,perts_tke &
124 ,pert_mynn_qv, pert_mynn_qc, pert_mynn_t &
125 ,pert_mynn_qke &
126 ! variables required for camuwpbl scheme
127 , z_at_w,cldfra_old_mp,cldfra, rthratenlw &
128 , tauresx2d,tauresy2d &
129 , tpert2d,qpert2d,wpert2d,wsedl3d &
130 , turbtype3d,smaw3d &
131 , fnm,fnp &
132 ! variables required for camuwpbl scheme (optional)
133 ,qnc_curr,f_qnc,qni_curr,f_qni,rqniblten &
134 ,IS_CAMMGMP_USED &
135 ! for grims shallow convection with ysupbl/MYNN
136 , wstar,delta &
137 ! for pbl mixing of scalars and tracers
138 & ,scalar,scalar_tend,num_scalar &
139 & ,tracer,tracer_tend,num_tracer &
140 & ,scalar_pblmix,tracer_pblmix &
141 #if (WRF_CHEM == 1)
142 ,mynn_chem_vertmx,chem,vd,nchem,kdvel &
143 & ,ndvel,num_vert_mix &
144 #endif
145 !
146 ! FASDAS
147 !
148 ,QNORM, fasdas &
149 !
150 ! END FASDAS
151 !
152 )
154 !------------------------------------------------------------------
155 #if (EM_CORE==1)
156 USE module_state_description, ONLY : &
157 YSUSCHEME,MRFSCHEME,GFSSCHEME,MYJPBLSCHEME,ACMPBLSCHEME,&
158 QNSEPBLSCHEME,MYNNPBLSCHEME,BOULACSCHEME, &
159 CAMUWPBLSCHEME,BEPSCHEME,BEP_BEMSCHEME,MYJSFCSCHEME, &
160 FITCHSCHEME,SHINHONGSCHEME, &
161 TEMFPBLSCHEME,GBMPBLSCHEME,EEPSSCHEME,KEPSSCHEME, &
162 MAVSCHEME, & ! Yulong add for WLM
163 CAMMGMPSCHEME,p_qi,p_qni,p_qnc,param_first_scalar,& !CAMMGMPSCHEME, p_qni,p_qnc is used for camuwpbl scheme
164 p_qnwfa,p_qnifa,p_qnbca
165 #if ( WRFPLUS == 1 )
166 USE module_state_description, ONLY : SURFDRAGSCHEME
167 #endif
168 #else
169 USE module_state_description, ONLY : &
170 YSUSCHEME,MRFSCHEME,GFSSCHEME,MYJPBLSCHEME,ACMPBLSCHEME &
171 , QNSEPBLSCHEME, p_qi,param_first_scalar &
172 , TEMFPBLSCHEME, GFSEDMFSCHEME &
173 , CAMUWPBLSCHEME &
174 , FITCHSCHEME, SHINHONGSCHEME &
175 , MAVSCHEME ! Yulong add for WLM
176 , GBMPBLSCHEME, MYJSFCSCHEME
177 #endif
179 USE module_model_constants
181 ! *** add new modules of schemes here
183 USE module_bl_myjpbl
184 USE module_bl_qnsepbl
185 USE module_bl_ysu
186 USE module_bl_shinhong
187 USE module_bl_mrf
188 USE module_bl_gfs
189 USE module_bl_acm
190 USE module_bl_gwdo
191 USE module_bl_gwdo_gsl
192 USE module_bl_myjurb
193 USE module_bl_boulac
194 #if ( WRFPLUS == 1 )
195 USE module_bl_surface_drag
196 #endif
197 USE module_bl_camuwpbl_driver, only:CAMUWPBL
198 USE module_bl_temf
199 USE module_bl_mfshconvpbl
200 USE module_bl_gbmpbl
201 #if (EM_CORE==1)
202 USE module_bl_mynn_wrapper
203 USE module_bl_eeps
204 USE module_bl_keps
205 USE module_bl_fogdes
206 USE module_wind_fitch
207 USE module_wind_mav ! Yulong add for WLM
208 #endif
209 use module_ra_gfdleta, only: cal_mon_day
211 ! This driver calls subroutines for the PBL parameterizations.
212 !
213 ! pbl scheme:
214 ! 1. ysupbl
215 ! 2. myjpbl
216 ! 4. qnsepbl
217 ! 5. mynnpbl2
218 ! 6. mynnpbl3
219 ! 7. acmpbl
220 ! 8. boulacpbl
221 ! 9. camuwpbl
222 ! 10. temfpbl
223 ! 11. shinhongpbl
224 ! 17. kepsstpe
225 ! 98. surface_drag
226 ! 99. mrfpbl
227 ! 12. gbmpbl
228 !
229 !------------------------------------------------------------------
230 IMPLICIT NONE
231 !======================================================================
232 ! Grid structure in physics part of WRF
233 !----------------------------------------------------------------------
234 ! The horizontal velocities used in the physics are unstaggered
235 ! relative to temperature/moisture variables. All predicted
236 ! variables are carried at half levels except w, which is at full
237 ! levels. Some arrays with names (*8w) are at w (full) levels.
238 !
239 !----------------------------------------------------------------------
240 ! In WRF, kms (smallest number) is the bottom level and kme (largest
241 ! number) is the top level. In your scheme, if 1 is at the top level,
242 ! then you have to reverse the order in the k direction.
243 !
244 ! kme - half level (no data at this level)
245 ! kme ----- full level
246 ! kme-1 - half level
247 ! kme-1 ----- full level
248 ! .
249 ! .
250 ! .
251 ! kms+2 - half level
252 ! kms+2 ----- full level
253 ! kms+1 - half level
254 ! kms+1 ----- full level
255 ! kms - half level
256 ! kms ----- full level
257 !
258 !======================================================================
259 ! Definitions
260 !-----------
261 ! Rho_d dry density (kg/m^3)
262 ! Theta_m moist potential temperature (K)
263 ! Qv water vapor mixing ratio (kg/kg)
264 ! Qc cloud water mixing ratio (kg/kg)
265 ! Qr rain water mixing ratio (kg/kg)
266 ! Qi cloud ice mixing ratio (kg/kg)
267 ! Qs snow mixing ratio (kg/kg)
268 ! QNC cloud Liq number concentration (#/kg) !For CAMUWPBL scheme
269 ! QNI cloud ice number concentration (#/kg) !For CAMUWPBL scheme
270 !-----------------------------------------------------------------
271 !-- RUBLTEN U tendency due to
272 ! PBL parameterization (m/s^2)
273 !-- RVBLTEN V tendency due to
274 ! PBL parameterization (m/s^2)
275 !-- RTHBLTEN Theta tendency due to
276 ! PBL parameterization (K/s)
277 !-- RQVBLTEN Qv tendency due to
278 ! PBL parameterization (kg/kg/s)
279 !-- RQCBLTEN Qc tendency due to
280 ! PBL parameterization (kg/kg/s)
281 !-- RQIBLTEN Qi tendency due to
282 ! PBL parameterization (kg/kg/s)
283 !-- RQNIBLTEN Qni tendency due to
284 ! PBL parameterization (#/kg/s) !For CAMUWPBL scheme
285 !-- id WRF grid id (optional, only needed by turbine drag schemes)
286 !-- itimestep number of time steps
287 !-- GLW downward long wave flux at ground surface (W/m^2)
288 !-- GSW downward short wave flux at ground surface (W/m^2)
289 !-- EMISS surface emissivity (between 0 and 1)
290 !-- TSK surface temperature (K)
291 !-- TMN soil temperature at lower boundary (K)
292 !-- XLAND land mask (1 for land, 2 for water)
293 !-- ZNT thermal roughness length (m)
294 !-- MZNT momentum roughness length (m)
295 !-- MAVAIL surface moisture availability (between 0 and 1)
296 !-- UST u* in similarity theory (m/s)
297 !-- MOL T* (similarity theory) (K)
298 !-- HOL PBL height over Monin-Obukhov length
299 !-- PBLH PBL height (m)
300 !-- CAPG heat capacity for soil (J/K/m^3)
301 !-- THC thermal inertia (Cal/cm/K/s^0.5)
302 !-- SNOWC flag indicating snow coverage (1 for snow cover)
303 !-- HFX upward heat flux at the surface (W/m^2)
304 !-- QFX upward moisture flux at the surface (kg/m^2/s)
305 !-- REGIME flag indicating PBL regime (stable, unstable, etc.)
306 !-- exch_m exchange coefficient for momentum, m^2/s
307 !-- exch_h exchange coefficient for heat, K m/s
308 !-- exch_tke exchange coeff. for TKE [enhanced], m^2/s (gbmpbl scheme)
309 !-- rthraten tendency from radiation, used in GBM PBL scheme
310 !-- akhs sfc exchange coefficient of heat/moisture from MYJ
311 !-- akms sfc exchange coefficient of momentum from MYJ
312 !-- tke_pbl turbulence kinetic energy from PBL schemes (m^2/s^2)
313 !-- el_pbl length scale from PBL schemes (m)
314 !-- wu_tur turbulent flux of momentum (x) (m^2/s^2)
315 !-- wv_tur turbulent flux of momentum (y) (m^2/s^2)
316 !-- wt_tur turbulent flux of potential temperature (K m/s)
317 !-- wq_tur turbulent flux of water vapor (- m/s)
318 !-- te_temf Total energy from TEMF BL scheme
319 !-- km_temf Exchange coefficient for momentum from TEMF BL scheme
320 !-- kh_temf Exchange coefficient for heat from TEMF BL scheme
321 !-- shf_temf Sensible heat flux from TEMF BL scheme
322 !-- qf_temf Water vapor flux from TEMF BL scheme
323 !-- uw_temf Momentum flux in U direction from TEMF BL scheme
324 !-- vw_temf Momentum flux in V direction from TEMF BL scheme
325 !-- wupd_temf Updraft velocity from TEMF BL scheme
326 !-- mf_temf Mass flux from TEMF BL scheme
327 !-- thup_temf Updraft thetal from TEMF BL scheme
328 !-- qtup_temf Updraft qt from TEMF BL scheme
329 !-- qlup_temf Updraft ql from TEMF BL scheme
330 !-- cf3d_temf 3D cloud fraction from TEMF PBL
331 !-- cfm_temf Column cloud fraction from TEMF PBL
332 !-- exch_temf Surface exchange coefficient (as for moisture) from TEMF surface layer scheme
333 !-- flhc Surface exchange coefficient for heat (for TEMF)
334 !-- flqc Surface exchange coefficient for moisture (for TEMF)
335 !-- thz0 potential temperature at roughness length (K)
336 !-- uz0 u wind component at roughness length (m/s)
337 !-- vz0 v wind component at roughness length (m/s)
338 !-- qsfc specific humidity at lower boundary (kg/kg)
339 !-- UOCE sea surface zonal currents (m s-1)
340 !-- VOCE sea surface meridional currents (m s-1)
341 !-- th2 diagnostic 2-m theta from surface layer and lsm
342 !-- t2 diagnostic 2-m temperature from surface layer and lsm
343 !-- q2 diagnostic 2-m mixing ratio from surface layer and lsm
344 !-- lowlyr index of lowest model layer above ground
345 !-- rr dry air density (kg/m^3)
346 !-- u_phy u-velocity interpolated to theta points (m/s)
347 !-- v_phy v-velocity interpolated to theta points (m/s)
348 !-- th_phy potential temperature (K)
349 !-- p_phy pressure (Pa)
350 !-- pi_phy exner function (dimensionless)
351 !-- p8w pressure at full levels (Pa)
352 !-- t_phy temperature (K)
353 !-- dz8w dz between full levels (m)
354 !-- z height above sea level (m)
355 !-- DX horizontal space interval (m)
356 !-- DT time step (second)
357 !-- n_moist number of moisture species
358 !-- PSFC pressure at the surface (Pa)
359 !-- TSLB
360 !-- ZS
361 !-- DZS
362 !-- num_soil_layers number of soil layer
363 !-- IFSNOW ifsnow=1 for snow-cover effects
364 !-- z_at_w Height above sea level at layer interfaces (m)
365 !-- cldfra Cloud fraction [unitless]
366 !-- cldfra_old_mp Cloud fraction [unitless]
367 !-- rthratenlw Tendency for LW ( K/s)
368 !-- tauresx2d X-COMP OF RESIDUAL STRESS(m^2/s^2)
369 !-- tauresy2d Y-COMP OF RESIDUAL STRESS(m^2/s^2)
370 !-- tpert2d Convective temperature excess (K)
371 !-- qpert2d Convective humidity excess (kg/kg)
372 !-- wpert2d Turbulent velocity excess (m/s)
373 !-- wsedl3d Sedimentation velocity of stratiform liquid cloud droplet (m/s)
374 !-- turbtype3d Turbulent interface types [ no unit ]
375 !-- smaw3d Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
376 !
377 !-- P_QV species index for water vapor
378 !-- P_QC species index for cloud water
379 !-- P_QR species index for rain water
380 !-- P_QI species index for cloud ice
381 !-- P_QNC species index for cloud liq number concentration !For CAMUWPBL scheme
382 !-- P_QNI species index for cloud ice number concentration !For CAMUWPBL scheme
383 !-- P_QS species index for snow
384 !-- P_QG species index for graupel
385 !-- ids start index for i in domain
386 !-- ide end index for i in domain
387 !-- jds start index for j in domain
388 !-- jde end index for j in domain
389 !-- kds start index for k in domain
390 !-- kde end index for k in domain
391 !-- ims start index for i in memory
392 !-- ime end index for i in memory
393 !-- jms start index for j in memory
394 !-- jme end index for j in memory
395 !-- kms start index for k in memory
396 !-- kme end index for k in memory
397 !-- jts start index for j in tile
398 !-- jte end index for j in tile
399 !-- kts start index for k in tile
400 !-- kte end index for k in tile
401 !
402 !******************************************************************
403 !------------------------------------------------------------------
404 !
407 INTEGER, INTENT(IN ) :: bl_pbl_physics, ra_lw_physics,sf_sfclay_physics,sf_urban_physics,windfarm_opt
408 INTEGER, INTENT(IN ) :: ysu_topdown_pblmix
409 INTEGER, OPTIONAL, INTENT(IN ) :: shinhong_tke_diag
410 INTEGER, OPTIONAL, INTENT(IN ) :: scalar_pblmix, tracer_pblmix
412 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
413 ims,ime, jms,jme, kms,kme, &
414 kts,kte, num_tiles
415 INTEGER, INTENT(IN ) :: num_scalar, num_tracer
417 INTEGER, DIMENSION(num_tiles), INTENT(IN) :: &
418 & i_start,i_end,j_start,j_end
420 INTEGER, INTENT(IN ) :: itimestep,STEPBL
421 INTEGER, DIMENSION( ims:ime , jms:jme ), &
422 INTENT(IN ) :: LOWLYR
423 !
424 LOGICAL, INTENT(IN ) :: warm_rain
425 LOGICAL, INTENT(IN ) :: is_CAMMGMP_used !BSINGH:01/31/2013: Added for CAMUWPBL
426 LOGICAL, OPTIONAL, INTENT(IN ) :: restart,cycling !used by mynn
428 REAL, DIMENSION( kms:kme ), &
429 OPTIONAL, INTENT(IN ) :: znu, &
430 znw
431 !
432 REAL, INTENT(IN ) :: DT,DX,DY
433 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) , OPTIONAL :: &
434 dx2d, area2d
435 REAL, INTENT(IN ),OPTIONAL :: bldt
436 REAL, INTENT(IN ),OPTIONAL :: curr_secs
437 LOGICAL, INTENT(IN ),OPTIONAL :: adapt_step_flag
438 REAL, INTENT(INOUT),OPTIONAL :: bldtacttime
439 ! Optional for Wind Turbine Parameterizations
440 REAL, DIMENSION( ims:ime, kms:kme ,jms:jme ), &
441 INTENT(IN), OPTIONAL :: phb
442 REAL, DIMENSION( ims:ime, jms:jme ), &
443 INTENT(IN), OPTIONAL :: xlat_u,xlong_u,xlat_v,xlong_v
445 ! Yulong add for WLM
446 INTEGER, INTENT(IN ) :: windfarm_wake_model, windfarm_overlap_method
448 REAL, DIMENSION( ims:ime, kms:kme ,jms:jme ), &
449 INTENT(IN), OPTIONAL :: w
450 !
451 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
452 INTENT(IN ) :: p_phy, &
453 pi_phy, &
454 p8w, &
455 rho, &
456 t_phy, &
457 u_phy, &
458 v_phy, &
459 dz8w, &
460 z
461 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) :: th_phy
463 ! WRF-Solar EPS
464 integer, intent (in) :: multi_perturb
465 LOGICAL, intent (in) :: pert_mynn
466 real, intent (in) :: pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke
467 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_qvapor, &
468 perts_qcloud, perts_th, perts_tke
470 !1D variables required for CAMUWPBL scheme
471 REAL , DIMENSION( kms:kme ) , &
472 INTENT(IN ) , OPTIONAL :: fnm, & !Factors for interpolation at "w" grid (interfaces)
473 fnp
474 !3D Variables for camuwpbl scheme
475 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
476 INTENT(IN ), OPTIONAL :: z_at_w, &
477 cldfra_old_mp, &
478 cldfra, &
479 rthratenlw, &
480 wsedl3d
482 !For EEPS PBL scheme
483 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
484 INTENT(INOUT ) :: pek, pep
485 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
486 INTENT(INOUT) :: pek_adv,pep_adv
488 !2D Variables required by camuwpbl scheme
489 REAL, DIMENSION( ims:ime , jms:jme ), &
490 INTENT(INOUT ), OPTIONAL :: tauresx2d, &
491 tauresy2d, &
492 tpert2d, &
493 qpert2d, &
494 wpert2d
497 !3D Variables for camuwpbl scheme - out
498 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
499 INTENT(OUT) , OPTIONAL :: turbtype3d, &
500 smaw3d
501 !
502 ! for grims shallow convection with ysupbl
503 !
504 REAL, DIMENSION( ims:ime, jms:jme ) , &
505 INTENT(INOUT ), OPTIONAL :: wstar
506 REAL, DIMENSION( ims:ime, jms:jme ) , &
507 INTENT(INOUT ), OPTIONAL :: delta
508 !
509 REAL, DIMENSION( ims:ime , jms:jme ), &
510 INTENT(IN ) :: XLAND, &
511 HT, &
512 PSIM, &
513 PSIH, &
514 FM, &
515 FHH, &
516 GZ1OZ0, &
517 BR, &
518 F, &
519 CHKLOWQ
521 !
522 REAL, DIMENSION( ims:ime, jms:jme ) , &
523 INTENT(INOUT) :: TSK, &
524 UST, &
525 PBLH, &
526 HFX, &
527 QFX, &
528 ZNT, &
529 QSFC, &
530 AKHS, &
531 AKMS, &
532 MIXHT, &
533 QZ0, &
534 THZ0, &
535 UZ0, &
536 VZ0, &
537 CT, &
538 GRDFLX, &
539 U10, &
540 V10, &
541 UOCE, &
542 VOCE, &
543 T2, &
544 POWER, &
545 WSPD
546 REAL, DIMENSION( ims:ime, jms:jme ) , &
547 INTENT(INOUT) :: HPBL_HOLD
549 !
550 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
551 INTENT(INOUT) :: RUBLTEN, &
552 RVBLTEN, &
553 RTHBLTEN, &
554 EXCH_H,EXCH_M,TKE_PBL, &
555 DISS_PBL,TPE_PBL, &
556 TKE_ADV,DISS_ADV,TPE_ADV, &
557 PR_PBL,RTHRATEN
559 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
560 INTENT(OUT) :: WU_TUR,WV_TUR,WT_TUR,WQ_TUR
561 !
563 !MYNN
564 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
565 INTENT(INOUT) :: tsq,qsq,cov, & !,k_m,k_h,k_q &
566 qke,Sh3d,Sm3d, &
567 dqke,qWT,qSHEAR,qBUOY,qDISS, &
568 qc_bl,qi_bl,cldfra_bl
570 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
571 INTENT(INOUT) :: edmf_a,edmf_w,edmf_thl, & !JOE- MYNN edmf
572 edmf_qt,edmf_ent,edmf_qc, & !JOE- MYNN edmf
573 sub_thl3D,sub_sqv3D, &
574 det_thl3D,det_sqv3D
576 INTEGER, OPTIONAL, INTENT(IN) :: grav_settling, &
577 bl_mynn_cloudpdf, &
578 bl_mynn_mixlength, &
579 bl_mynn_edmf, &
580 bl_mynn_edmf_mom, &
581 bl_mynn_edmf_tke, &
582 bl_mynn_mixscalars, &
583 bl_mynn_output, &
584 bl_mynn_cloudmix, &
585 bl_mynn_mixqt, &
586 icloud_bl, &
587 tke_budget
589 REAL, OPTIONAL, INTENT(IN) :: bl_mynn_closure
591 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
592 & INTENT(INOUT) :: qke_adv
593 LOGICAL, OPTIONAL, INTENT(IN) :: bl_mynn_tkeadvect
595 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
596 & INTENT(INOUT):: vdfg
598 INTEGER, OPTIONAL, DIMENSION( ims:ime , jms:jme ), &
599 & INTENT(OUT) :: ktop_plume
600 REAL, OPTIONAL, DIMENSION( ims:ime , jms:jme ), &
601 & INTENT(OUT) :: maxwidth,maxMF,ztop_plume
603 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
604 INTENT(INOUT) :: qnwfa_curr,qnifa_curr,qnbca_curr
606 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
607 & INTENT(INOUT) :: exch_tke ! for GBM PBL scheme
610 INTEGER, OPTIONAL :: id
612 REAL, DIMENSION( ims:ime , jms:jme ), &
613 &OPTIONAL, INTENT(IN) :: &
614 & qcg, ch
615 REAL, DIMENSION( ims:ime , jms:jme ), &
616 &OPTIONAL, INTENT(INOUT) :: rmol
620 !
621 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
622 INTENT(OUT) :: EL_PBL
624 REAL, INTENT(IN) :: gmt, xtime
625 INTEGER, INTENT(IN) :: julday, julyr
626 REAL, OPTIONAL, DIMENSION( ims:ime, 0:287, jms:jme ), INTENT(IN) :: ahe
627 INTEGER, INTENT(IN) :: distributed_ahe_opt
629 REAL , INTENT(IN ) :: u_frame, &
630 v_frame
631 !
633 INTEGER, DIMENSION( ims:ime , jms:jme ), &
634 INTENT(INOUT) :: KPBL
636 REAL, DIMENSION( ims:ime , jms:jme ), &
637 INTENT(IN) :: XICE, SNOW, LH
639 ! Stochastic parameter perturbations
640 INTEGER, INTENT(IN), OPTIONAL :: spp_pbl
641 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN),OPTIONAL ::pattern_spp_pbl
643 ! Bep changes: variable added for urban
644 real, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::FRC_URB2D ! URBAN Landuse fraction
645 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep ! Implicit component for the momemtum in X-direction
646 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep ! Implicit component for the momemtum in Y-direction
647 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep ! Implicit component for the Pot. Temp.
648 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep ! Implicit component for Moisture
650 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep ! Implicit component for the TKE
651 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep ! Explicit component for the momemtum in X-direction
652 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep ! Explicit component for the momemtum in Y-direction
653 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep ! Explicit component for the Pot. Temp.
654 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep ! Explicit component for Moisture
656 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep ! Explicit component for the TKE
658 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).
659 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep ! Length scale (lb in formula (22) ofthe BLM paper).
660 ! urban surface and volumes
661 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep ! surfaces
662 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep ! volumes
664 ! Bep changes end
665 ! New variables for TEMF scheme
666 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
667 INTENT(INOUT) :: te_temf
668 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
669 INTENT( OUT) :: km_temf, kh_temf, &
670 shf_temf,qf_temf,uw_temf,vw_temf, &
671 wupd_temf,mf_temf,thup_temf,qtup_temf, &
672 qlup_temf,cf3d_temf
673 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
674 INTENT(INOUT) :: flhc,flqc
675 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
676 INTENT( OUT) :: hd_temf, lcl_temf, hct_temf, cfm_temf
677 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), &
678 INTENT(INOUT) :: exch_temf
680 !
681 !
682 ! Optional
683 !
684 !
685 ! Flags relating to the optional tendency arrays declared above
686 ! Models that carry the optional tendencies will provdide the
687 ! optional arguments at compile time; these flags all the model
688 ! to determine at run-time whether a particular tracer is in
689 ! use or not.
690 !
691 LOGICAL, INTENT(IN), OPTIONAL :: &
692 f_qv &
693 ,f_qc &
694 ,f_qr &
695 ,f_qi &
696 ,f_qs &
697 ,f_qg &
698 ,f_qnc & !used in CAMUWPBL
699 ,f_qni & !used in CAMUWPBL
700 ,f_qnwfa &
701 ,f_qnifa &
702 ,f_qnbca
704 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
705 OPTIONAL, INTENT(INOUT) :: &
706 ! optional moisture tracers
707 ! 2 time levels; if only one then use CURR
708 qv_curr, qc_curr, qr_curr &
709 ,qi_curr, qs_curr, qg_curr &
710 ,qnc_curr,qni_curr & !used in CAMUWPBL
711 ,rqvblten,rqcblten,rqrblten &
712 ,rqiblten,rqsblten,rqgblten,rqniblten !rqniblten used in CAMUWPBL
714 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: & !local: added for MYNN
715 rqncblten,rqnwfablten,rqnifablten,rqnbcablten
717 REAL, DIMENSION( ims:ime, jms:jme ) , &
718 OPTIONAL , &
719 INTENT(INOUT) :: HOL, &
720 MOL, &
721 REGIME
722 REAL, DIMENSION( ims:ime, jms:jme ) , &
723 OPTIONAL , &
724 INTENT(IN) :: mut
725 !
726 INTEGER, OPTIONAL, INTENT(IN) :: gwd_opt,gwd_diags
727 REAL, OPTIONAL, INTENT(IN) :: p_top
728 !
729 real, dimension( ims:ime, kms:kme, jms:jme ) , &
730 optional , &
731 intent(inout ) :: dtaux3d,dtauy3d, &
732 dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, &
733 dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd
735 !
736 real, dimension( ims:ime, jms:jme ) , &
737 optional , &
738 intent(inout ) :: dusfcg,dvsfcg, &
739 dusfcg_ls,dvsfcg_ls,dusfcg_bl,dvsfcg_bl, &
740 dusfcg_ss,dvsfcg_ss,dusfcg_fd,dvsfcg_fd
742 !
743 real, dimension( ims:ime, jms:jme ) , &
744 optional , &
745 intent(in ) :: var2d, &
746 oc12d, &
747 oa1,oa2,oa3,oa4, &
748 ol1,ol2,ol3,ol4, &
749 sina,cosa
751 ! Addtional topographic statistics based on 2.5min (large-scale) and 30sec orographic
752 ! (small-scale) GSL drag schemes
753 real, dimension( ims:ime, jms:jme ) , &
754 optional , &
755 intent(in ) :: &
756 var2dls,oc12dls, &
757 oa1ls,oa2ls,oa3ls,oa4ls,ol1ls,ol2ls,ol3ls,ol4ls, &
758 var2dss,oc12dss, &
759 oa1ss,oa2ss,oa3ss,oa4ss,ol1ss,ol2ss,ol3ss,ol4ss
761 ! paj
762 REAL, OPTIONAL, DIMENSION( ims:ime , jms:jme ), & !mchen
763 INTENT(IN ) :: CTOPO, &
764 CTOPO2
765 real, dimension (:, :, :), allocatable :: qke_tmp
767 ! Variables and Diagnostic for QNSE and EDKF JP
768 INTEGER, INTENT(IN) :: mfshconv
770 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), OPTIONAL, &
771 INTENT(OUT) :: massflux_EDKF, entr_EDKF, detr_EDKF &
772 ,thl_up, thv_up, rt_up &
773 ,rv_up, rc_up, u_up, v_up &
774 ,frac_up,rc_mf
776 REAL , OPTIONAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_tracer),INTENT(INOUT) :: tracer
777 REAL , OPTIONAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_tracer),INTENT(INOUT) :: tracer_tend
778 REAL , OPTIONAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT) :: scalar
779 REAL , OPTIONAL ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT) :: scalar_tend
781 #if (WRF_CHEM == 1)
782 ! ACM Chem
783 INTEGER, INTENT(IN) :: nchem, kdvel,ndvel,num_vert_mix
784 REAL, INTENT(INOUT), DIMENSION( ims:ime, kms:kme, jms:jme,nchem ) :: CHEM
785 REAL, INTENT(IN) , DIMENSION( ims:ime, kdvel, jms:jme, ndvel ) :: VD
786 LOGICAL, INTENT(IN) :: mynn_chem_vertmx
787 #endif
789 ! LOCAL VAR
791 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) ::v_phytmp
792 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) ::u_phytmp
794 REAL, DIMENSION( ims:ime, jms:jme ) :: TSKOLD, &
795 USTOLD, &
796 ZNTOLD, &
797 ZOL, &
798 PSFC
800 ! make these allocatable depending on the setting of idiff
801 ! Typically, we try to avoide allocating and deallocating local storage like this
802 ! so as not to fragment the stack. But at this point, the idiff = 1 case is disabled
803 ! (set to 0 for all cases) and has to be set manually by users who want to work with
804 ! it. When it becomes a more standard option, this should be redone, either defining
805 ! these as state with package clauses to turn them on and off and passing them in,
806 ! or pass in an integer flag that can be used to dimension the arrays to 1:1:1 as
807 ! local variables. JM 20100316
809 REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_u ! Implicit component for the momemtum in X-direction
810 REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_v ! Implicit component for the momemtum in Y-direction
811 REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_t ! Implicit component for the Pot. Temp.
812 REAL, ALLOCATABLE, DIMENSION( :, :, : )::a_q ! Implicit component for the water vapor
814 REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_u ! Explicit component for the momemtum in X-direction
815 REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_v ! Explicit component for the momemtum in Y-direction
816 REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_t ! Explicit component for the Pot. Temp.
817 REAL, ALLOCATABLE, DIMENSION( :, :, : )::b_q ! Explicit component for the water vapor
819 REAL, ALLOCATABLE, DIMENSION( :, :, : )::sf ! surfaces
820 REAL, ALLOCATABLE, DIMENSION( :, :, : )::vl ! volumes
822 REAL :: DTMIN,DTBL
823 !
824 INTEGER :: initflag
825 !
826 INTEGER :: i,J,K,NK,jj,ij,its,ite,jts,jte
827 LOGICAL :: radiation
828 LOGICAL :: flag_bep
829 LOGICAL :: flag_myjsfc
830 LOGICAL :: flag_qv, flag_qc, flag_qr, flag_qi, flag_qs, flag_qg, &
831 flag_qnc,flag_qni, flag_qnwfa, flag_qnifa, flag_qnbca
832 CHARACTER*256 :: message
833 REAL :: next_bl_time
834 LOGICAL :: run_param , doing_adapt_dt , decided
835 LOGICAL :: do_adapt
836 integer iu_bep,iurb,idiff
837 real seamask,thsk,zzz,unew,vnew,tnew,qnew,umom,vmom
838 REAL :: z0,z1,z2,w1,w2
839 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: TKE_windfarm ! Yulong add for WLM
840 INTEGER :: ihour, jmonth, jday
841 !
842 ! FASDAS
843 !
844 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT( OUT) , OPTIONAL :: QNORM
845 INTEGER, INTENT(IN ) :: fasdas
846 REAL :: scrp1
847 !
848 ! END FASDAS
849 !
850 ! To accommodate shared physics
851 character*256 :: errmsg
852 integer :: errflg
853 !
854 !------------------------------------------------------------------
855 !
856 !!!!!!!if using BEP set flag_bep to true
857 SELECT CASE(sf_urban_physics)
858 CASE (BEPSCHEME)
859 flag_bep=.true.
860 CASE (BEP_BEMSCHEME)
861 flag_bep=.true.
862 CASE DEFAULT
863 flag_bep=.false.
864 END SELECT
866 SELECT CASE(sf_sfclay_physics)
867 CASE (MYJSFCSCHEME)
868 flag_myjsfc=.true.
869 CASE DEFAULT
870 flag_myjsfc=.false.
871 END SELECT
872 !
873 flag_qv = .FALSE. ; IF ( PRESENT( F_QV ) ) flag_qv = F_QV
874 flag_qc = .FALSE. ; IF ( PRESENT( F_QC ) ) flag_qc = F_QC
875 flag_qr = .FALSE. ; IF ( PRESENT( F_QR ) ) flag_qr = F_QR
876 flag_qi = .FALSE. ; IF ( PRESENT( F_QI ) ) flag_qi = F_QI
877 flag_qs = .FALSE. ; IF ( PRESENT( F_QS ) ) flag_qs = F_QS
878 flag_qg = .FALSE. ; IF ( PRESENT( F_QG ) ) flag_qg = F_QG
879 flag_qnc = .FALSE. ; IF ( PRESENT( F_QNC ) ) flag_qnc = F_QNC !Used in CAMUWPBL
880 flag_qni = .FALSE. ; IF ( PRESENT( F_QNI ) ) flag_qni = F_QNI !Used in CAMUWPBL
881 flag_qnwfa = .FALSE. ; IF ( PRESENT( F_QNWFA ) ) flag_qnwfa = F_QNWFA
882 flag_qnifa = .FALSE. ; IF ( PRESENT( F_QNIFA ) ) flag_qnifa = F_QNIFA
883 flag_qnbca = .FALSE. ; IF ( PRESENT( F_QNBCA ) ) flag_qnbca = F_QNBCA
885 if (bl_pbl_physics .eq. 0) return
887 ! RAINBL in mm (Accumulation between PBL calls)
888 !
889 doing_adapt_dt = .FALSE.
890 IF ( PRESENT(adapt_step_flag) ) THEN
891 IF ( adapt_step_flag ) THEN
892 doing_adapt_dt = .TRUE.
893 IF ( bldtacttime .eq. 0. ) THEN
894 bldtacttime = CURR_SECS + bldt*60.
895 END IF
896 END IF
897 END IF
899 ! Do we run through this scheme or not?
901 ! Test 1: If this is the initial model time, then yes.
902 ! ITIMESTEP=1
903 ! Test 2: If the user asked for the pbl to be run every time step, then yes.
904 ! BLDT=0 or STEPBL=1
905 ! Test 3: If not adaptive dt, and this is on the requested pbl frequency, then yes.
906 ! MOD(ITIMESTEP,STEPBL)=0
907 ! Test 4: If using adaptive dt and the current time is past the last requested activate pbl time, then yes.
908 ! CURR_SECS >= BLDTACTTIME
910 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
911 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
912 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
914 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
915 ! pbl run.
917 run_param = .FALSE.
918 decided = .FALSE.
919 IF ( ( .NOT. decided ) .AND. &
920 ( itimestep .EQ. 1 ) ) THEN
921 run_param = .TRUE.
922 decided = .TRUE.
923 END IF
925 IF ( PRESENT(bldt) )THEN
926 IF ( ( .NOT. decided ) .AND. &
927 ( ( bldt .EQ. 0. ) .OR. ( stepbl .EQ. 1 ) ) ) THEN
928 run_param = .TRUE.
929 decided = .TRUE.
930 END IF
931 ELSE
932 IF ( ( .NOT. decided ) .AND. &
933 ( stepbl .EQ. 1 ) ) THEN
934 run_param = .TRUE.
935 decided = .TRUE.
936 END IF
937 END IF
939 IF ( ( .NOT. decided ) .AND. &
940 ( .NOT. doing_adapt_dt ) .AND. &
941 ( MOD(itimestep,stepbl) .EQ. 0 ) ) THEN
942 run_param = .TRUE.
943 decided = .TRUE.
944 END IF
946 IF ( ( .NOT. decided ) .AND. &
947 ( doing_adapt_dt ) .AND. &
948 ( curr_secs .GE. bldtacttime ) ) THEN
949 run_param = .TRUE.
950 decided = .TRUE.
951 bldtacttime = curr_secs + bldt*60
952 END IF
955 IF (run_param) THEN
956 radiation = .false.
957 IF (ra_lw_physics .gt. 0) radiation = .true.
959 !----
960 ! CALCULATE CONSTANT
962 DTMIN=DT/60.
963 ! PBL schemes need PBL time step for updates
965 if (PRESENT(adapt_step_flag)) then
966 if (adapt_step_flag) then
967 do_adapt = .TRUE.
968 else
969 do_adapt = .FALSE.
970 endif
971 else
972 do_adapt = .FALSE.
973 endif
975 if (PRESENT(BLDT)) then
976 if (bldt .eq. 0) then
977 DTBL = dt
978 ELSE
979 if (do_adapt) then
980 IF ( curr_secs .LT. 2. * dt ) THEN
981 call wrf_message("WARNING: When using an adaptive time-step the boundary layer"// &
982 " time-step should be 0 (i.e., equivalent to model time-step). ")
983 call wrf_message("In order to proceed, for boundary layer calculations, the "// &
984 "boundary layer time-step "// &
985 "will be rounded to the nearest minute," )
986 call wrf_message("possibly resulting in innacurate results.")
987 END IF
988 DTBL=bldt*60
989 else
990 DTBL=DT*STEPBL
991 endif
992 endif
993 else
994 DTBL=DT*STEPBL
995 endif
997 !!Alberto. If idiff=1, then compute the tendencies at the end in the routine diff3d
998 !!Alberto. If idiff=0, then compute the tendencies in each PBL routine (standard version).
999 idiff=0
1001 IF ( idiff .EQ. 1 ) THEN
1002 ALLOCATE (a_u(ims:ime,kms:kme,jms:jme)) ! Implicit component for the momemtum in X-direction
1003 ALLOCATE (a_v(ims:ime,kms:kme,jms:jme)) ! Implicit component for the momemtum in Y-direction
1004 ALLOCATE (a_t(ims:ime,kms:kme,jms:jme)) ! Implicit component for the Pot. Temp.
1005 ALLOCATE (a_q(ims:ime,kms:kme,jms:jme)) ! Implicit component for the water vapor
1006 ALLOCATE (b_u(ims:ime,kms:kme,jms:jme)) ! Explicit component for the momemtum in X-direction
1007 ALLOCATE (b_v(ims:ime,kms:kme,jms:jme)) ! Explicit component for the momemtum in Y-direction
1008 ALLOCATE (b_t(ims:ime,kms:kme,jms:jme)) ! Explicit component for the Pot. Temp.
1009 ALLOCATE (b_q(ims:ime,kms:kme,jms:jme)) ! Explicit component for the water vapor
1010 ALLOCATE (sf(ims:ime,kms:kme,jms:jme) ) ! surfaces
1011 ALLOCATE (vl(ims:ime,kms:kme,jms:jme) ) ! volumes
1012 ENDIF
1014 ! SAVE OLD VALUES
1016 !$OMP PARALLEL DO &
1017 !$OMP PRIVATE ( ij,i,j,k )
1019 DO ij = 1 , num_tiles
1020 DO j=j_start(ij),j_end(ij)
1021 DO i=i_start(ij),i_end(ij)
1022 TSKOLD(i,j)=TSK(i,j)
1023 USTOLD(i,j)=UST(i,j)
1024 ZNTOLD(i,j)=ZNT(i,j)
1025 HPBL_HOLD(i,j)=PBLH(i,j)
1027 DO k=kts,kte
1028 v_phytmp(i,k,j)=v_phy(i,k,j)+v_frame
1029 u_phytmp(i,k,j)=u_phy(i,k,j)+u_frame
1030 ENDDO
1032 ! PSFC : in Pa
1034 PSFC(I,J)=p8w(I,kms,J)
1036 DO k=kts,min(kte+1,kde)
1037 RTHBLTEN(I,K,J)=0.
1038 RUBLTEN(I,K,J)=0.
1039 RVBLTEN(I,K,J)=0.
1040 IF ( PRESENT( RQCBLTEN )) RQCBLTEN(I,K,J)=0.
1041 IF ( PRESENT( RQVBLTEN )) RQVBLTEN(I,K,J)=0.
1042 ENDDO
1044 IF (flag_QI .AND. PRESENT(RQIBLTEN) ) THEN
1045 DO k=kts,min(kte+1,kde)
1046 RQIBLTEN(I,K,J)=0.
1047 ENDDO
1048 ENDIF
1049 !Following if condition is added for CAMUWPBL scheme
1050 IF (flag_QNI .AND. PRESENT(RQNIBLTEN) ) THEN
1051 DO k=kts,min(kte+1,kde)
1052 RQNIBLTEN(I,K,J)=0.
1053 ENDDO
1054 ENDIF
1055 ENDDO
1056 ENDDO
1058 IF (bl_pbl_physics .EQ. QNSEPBLSCHEME ) THEN
1059 ! use u_phytmp instead of wthv_mf, and v_phytmp instead of lm_bl89
1060 DO j=j_start(ij),j_end(ij)
1061 DO i=i_start(ij),i_end(ij)
1062 DO k=kms,kme
1063 u_phytmp(i,k,j)=0.
1064 v_phytmp(i,k,j)=0.
1065 ENDDO
1066 ENDDO
1067 ENDDO
1068 ENDIF
1070 IF ( idiff.eq.1 ) THEN
1071 !Alberto
1072 ! Here we should put a switch:
1073 ! if we do not use BEP, just initialize the values of the a, and b to zero, except for the lowest model level,
1074 ! where the heat and momentum fluxes are introduced
1075 ! if we use BEP, use the values computed by BEP weigthed by the urban fraction.
1076 ! for the moment I put a simple "if" but it should be changed to a more elegant way after(using the CASE, maybe?)
1077 !!!!!!!!!!!!!!
1078 ! 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.
1079 ! if we need to be as tight as possible for the memory we can think how to reduce the storage.
1080 !!!!!!!!!!!!!!!!!!
1081 ! This stuff is becoming large, probably better to put in a subroutine
1082 !!!!!!!!!!!!!!!!!!!
1083 ! preparing the a, b, sf, and vl terms
1085 IF (flag_bep) THEN
1086 do j=j_start(ij),j_end(ij)
1087 do i=i_start(ij),i_end(ij)
1088 do k=kts,kte
1089 a_u(i,k,j)=a_u_bep(i,k,j)
1090 a_v(i,k,j)=a_v_bep(i,k,j)
1091 a_t(i,k,j)=a_t_bep(i,k,j)
1092 a_q(i,k,j)=a_q_bep(i,k,j)
1093 b_u(i,k,j)=b_u_bep(i,k,j)
1094 b_v(i,k,j)=b_v_bep(i,k,j)
1095 b_t(i,k,j)=b_t_bep(i,k,j)
1096 b_q(i,k,j)=b_q_bep(i,k,j)
1097 vl(i,k,j)=vl_bep(i,k,j)
1098 sf(i,k,j)=sf_bep(i,k,j)
1099 enddo
1100 sf(i,kte+1,j)=1.
1101 enddo
1102 enddo
1103 ELSE
1104 do j=j_start(ij),j_end(ij)
1105 do i=i_start(ij),i_end(ij)
1106 do k=kts,kte
1107 a_u(i,k,j)=0.
1108 a_v(i,k,j)=0.
1109 a_t(i,k,j)=0.
1110 a_q(i,k,j)=0.
1111 b_u(i,k,j)=0.
1112 b_v(i,k,j)=0.
1113 b_t(i,k,j)=0.
1114 b_q(i,k,j)=0.
1115 vl(i,k,j)=1.
1116 sf(i,k,j)=1.
1117 enddo
1118 sf(i,kte+1,j)=1.
1119 ! fix the values for the surface fluxes (source/sink terms)
1120 ! here for momentum the resolution is done implicitely
1121 ! while for heat and moisture is done explicitely
1122 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)
1123 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)
1124 b_t(i,1,j)=hfx(i,j)/rho(i,1,j)/CP/dz8w(i,1,j)
1125 b_q(i,1,j)=qfx(i,j)/rho(i,1,j)/dz8w(i,1,j)
1126 !! if you want to solve also the heat and moisture fluxes implicitely, uncomment the following lines.
1127 !! (of course you will need the values of thz0,qz0,akhs).
1128 ! b_t(i,1,j)=akhs(i,j)*(thz0(I,J))/dz8w(i,1,j)
1129 ! b_q(i,1,j)=akhs(i,j)*(qz0(I,J))/dz8w(i,1,j)
1130 ! a_t(i,1,j)=-1.*akhs(i,j)/dz8w(i,1,j)
1131 ! a_q(i,1,j)=-1.*akhs(i,j)/dz8w(i,1,j)
1132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1133 enddo
1134 enddo
1136 ENDIF
1138 endif
1142 !Alberto. From this point if some PBL scheme has a non local term
1143 ! (not dependent on the local values of the variable)
1144 ! this can be added to b_t, b_q, or b_u, b_v.
1145 !!!!!!!!!!!!!!!!!!!
1147 ENDDO
1148 !$OMP END PARALLEL DO
1149 !
1150 !$OMP PARALLEL DO &
1151 !$OMP PRIVATE ( ij, i,j,k, its, ite, jts, jte, z0, z1, z2, w1, w2, message, initflag )
1152 DO ij = 1 , num_tiles
1154 its = i_start(ij)
1155 ite = i_end(ij)
1156 jts = j_start(ij)
1157 jte = j_end(ij)
1159 pbl_select: SELECT CASE(bl_pbl_physics)
1161 CASE (TEMFPBLSCHEME)
1162 ! WA Test
1163 ! WRITE( message , * ) 'Calling TEMF PBL scheme, timestep = ', itimestep
1164 ! CALL wrf_message ( message )
1165 ! print *,'Calling TEMF PBL scheme'
1166 CALL wrf_debug(100,'in TEMF PBL')
1167 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1168 PRESENT( qi_curr ) .AND. &
1169 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1170 PRESENT( rqiblten ) .AND. &
1171 PRESENT( te_temf ) .AND. PRESENT( cfm_temf ) .AND. &
1172 PRESENT( hol ) ) THEN
1173 CALL temfpbl( &
1174 U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
1175 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
1176 ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy,RHO=rho &
1177 ,RUBLTEN=rublten,RVBLTEN=rvblten &
1178 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1179 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1180 ,FLAG_QI=flag_qi &
1181 ,g=g,cp=cp,rcp=rcp,r_d=r_d,r_v=r_v,cpv=cpv &
1182 ,Z=z,XLV=XLV,PSFC=PSFC &
1183 ,P_TOP=p_top &
1184 ,ZNT=znt,HT=ht,UST=ust,ZOL=zol,HOL=hol,HPBL=pblh &
1185 ,PSIM=psim,PSIH=psih,XLAND=xland &
1186 ,HFX=hfx,QFX=qfx,TSK=tskold,QSFC=qsfc,GZ1OZ0=gz1oz0 &
1187 ,U10=u10,V10=v10,T2=t2 &
1188 ,WSPD=wspd,BR=br,DT=dtbl,DTMIN=dtmin,KPBL2D=kpbl &
1189 ,SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0 &
1190 ,EP1=ep_1,EP2=ep_2,KARMAN=karman,EOMEG=eomeg &
1191 ,STBOLT=stbolt &
1192 ,te_temf=te_temf,kh_temf=kh_temf,km_temf=km_temf &
1193 ,shf_temf=shf_temf,qf_temf=qf_temf &
1194 ,uw_temf=uw_temf,vw_temf=vw_temf &
1195 ,hd_temf=hd_temf,lcl_temf=lcl_temf,hct_temf=hct_temf &
1196 ,wupd_temf=wupd_temf,mf_temf=mf_temf &
1197 ,thup_temf=thup_temf,qtup_temf=qtup_temf,qlup_temf=qlup_temf &
1198 ,cf3d_temf=cf3d_temf,cfm_temf=cfm_temf &
1199 ,flhc=flhc,flqc=flqc,exch_temf=exch_temf &
1200 ,fCor=f &
1201 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1202 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1203 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1204 )
1205 ELSE
1206 CALL wrf_error_fatal('Lack arguments to call TEMF pbl')
1207 ENDIF
1209 CASE (YSUSCHEME)
1210 CALL wrf_debug(100,'in YSU PBL')
1211 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1212 PRESENT( qi_curr ) .AND. &
1213 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1214 PRESENT( rqiblten ) .AND. &
1215 PRESENT( hol ) ) THEN
1216 !
1217 CALL ysu( &
1218 U3D=u_phytmp,V3D=v_phytmp,T3D=t_phy &
1219 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
1220 ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy &
1221 ,RUBLTEN=rublten,RVBLTEN=rvblten &
1222 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1223 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1224 ,FLAG_QI=flag_qi,FLAG_QC=flag_qc &
1225 ,CP=cp,G=g,ROVCP=rcp,RD=r_D,ROVG=rovg &
1226 ,DZ8W=dz8w,XLV=XLV,RV=r_v,PSFC=PSFC &
1227 ,ZNT=znt,UST=ust,HPBL=pblh &
1228 ,PSIM=fm,PSIH=fhh,XLAND=xland &
1229 ,HFX=hfx,QFX=qfx &
1230 ,U10=u10,V10=v10 &
1231 ,UOCE=uoce,VOCE=voce &
1232 ! paj
1233 ,CTOPO=ctopo,CTOPO2=ctopo2 &
1234 ,YSU_TOPDOWN_PBLMIX=ysu_topdown_pblmix &
1235 ,WSPD=wspd,BR=br,DT=dtbl,KPBL2D=kpbl &
1236 ,EP1=ep_1,EP2=ep_2,KARMAN=karman &
1237 ,EXCH_H=exch_h,EXCH_M=exch_m &
1238 ,RTHRATEN=RTHRATEN &
1239 ! for multilayer UCM
1240 ,IDIFF=idiff,FLAG_BEP=flag_bep,FRC_URB2D=frc_urb2d &
1241 ,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1242 ,A_Q_BEP=a_q_bep &
1243 ,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1244 ,B_T_BEP=b_t_bep,B_Q_BEP=b_q_bep &
1245 ,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep &
1246 ,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep &
1247 ! for grims shallow convection with ysupbl
1248 ,WSTAR=wstar,DELTA=delta &
1249 ,errmsg=errmsg,errflg=errflg &
1250 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1251 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1252 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1253 )
1254 !
1255 ! FASDAS
1256 !
1257 IF( fasdas == 1 ) THEN
1258 DO j=j_start(ij),j_end(ij)
1259 DO i=i_start(ij),i_end(ij)
1260 scrp1 = rqvblten(i,1,j)
1261 scrp1 = scrp1*(2.0*DT)/qv_curr(i,1,j)
1262 scrp1 = scrp1*100.0
1263 scrp1 = amax1(0.0,scrp1)
1264 scrp1 = amin1(5.0,scrp1)
1265 QNORM(I,J)= scrp1/100.0
1266 ENDDO
1267 ENDDO
1268 ENDIF
1269 !
1270 ! END FASDAS
1271 !
1273 ELSE
1274 WRITE ( message , FMT = '(A,7(L1,1X))' ) &
1275 'present: '// &
1276 'qv_curr, '// &
1277 'qc_curr, '// &
1278 'qi_curr, '// &
1279 'rqvblten, '// &
1280 'rqcblten, '// &
1281 'rqiblten, '// &
1282 'hol = ' , &
1283 PRESENT( qv_curr ) , &
1284 PRESENT( qc_curr ) , &
1285 PRESENT( qi_curr ) , &
1286 PRESENT( rqvblten ) , &
1287 PRESENT( rqcblten ) , &
1288 PRESENT( rqiblten ) , &
1289 PRESENT( hol )
1290 CALL wrf_debug(0,message)
1291 CALL wrf_error_fatal('Lack arguments to call YSU pbl')
1292 ENDIF
1294 CASE (SHINHONGSCHEME)
1295 CALL wrf_debug(100,'in SHINHONG PBL')
1296 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1297 PRESENT( qi_curr ) .AND. &
1298 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1299 PRESENT( rqiblten ) .AND. &
1300 PRESENT( hol ) ) THEN
1301 !
1302 CALL shinhong( &
1303 U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
1304 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
1305 ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy &
1306 ,RUBLTEN=rublten,RVBLTEN=rvblten &
1307 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1308 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1309 ,FLAG_QI=flag_qi &
1310 ,CP=cp,G=g,ROVCP=rcp,RD=r_D,ROVG=rovg &
1311 ,DZ8W=dz8w,XLV=XLV,RV=r_v,PSFC=PSFC &
1312 ,ZNU=znu,ZNW=znw,P_TOP=p_top &
1313 ,ZNT=znt,UST=ust,HPBL=pblh &
1314 ,PSIM=fm,PSIH=fhh,XLAND=xland &
1315 ,HFX=hfx,QFX=qfx &
1316 ,U10=u10,V10=v10 &
1317 ! paj
1318 ,CTOPO=ctopo,CTOPO2=ctopo2 &
1319 ,SHINHONG_TKE_DIAG=shinhong_tke_diag &
1320 ,WSPD=wspd,BR=br,DT=dtbl,KPBL2D=kpbl &
1321 ,EP1=ep_1,EP2=ep_2,KARMAN=karman &
1322 ,EXCH_H=exch_h,REGIME=regime &
1323 ! for grims shallow convection with shinhongpbl
1324 ,WSTAR=wstar,DELTA=delta &
1325 ,TKE_PBL=tke_pbl,EL_PBL=el_pbl,CORF=f &
1326 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1327 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1328 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1329 ,DX=dx,DY=dy &
1330 )
1331 ELSE
1332 WRITE ( message , FMT = '(A,7(L1,1X))' ) &
1333 'present: '// &
1334 'qv_curr, '// &
1335 'qc_curr, '// &
1336 'qi_curr, '// &
1337 'rqvblten, '// &
1338 'rqcblten, '// &
1339 'rqiblten, '// &
1340 'hol = ' , &
1341 PRESENT( qv_curr ) , &
1342 PRESENT( qc_curr ) , &
1343 PRESENT( qi_curr ) , &
1344 PRESENT( rqvblten ) , &
1345 PRESENT( rqcblten ) , &
1346 PRESENT( rqiblten ) , &
1347 PRESENT( hol )
1348 CALL wrf_debug(0,message)
1349 CALL wrf_error_fatal('Lack arguments to call SHINHONG pbl')
1350 ENDIF
1352 CASE (MRFSCHEME)
1353 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1354 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1355 PRESENT( hol ) .AND. &
1356 .TRUE. ) THEN
1358 CALL wrf_debug(100,'in MRF')
1359 CALL mrf( &
1360 U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
1361 ,QV3D=qv_curr &
1362 ,QC3D=qc_curr &
1363 ,QI3D=qi_curr &
1364 ,P3D=p_phy,PI3D=pi_phy &
1365 ,RUBLTEN=rublten,RVBLTEN=rvblten &
1366 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1367 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1368 ,CP=cp,G=g,ROVCP=rcp,R=r_d,ROVG=rovg &
1369 ,DZ8W=dz8w,Z=z,XLV=xlv,RV=r_v,PSFC=psfc &
1370 ,P1000MB=p1000mb &
1371 ,ZNT=znt,UST=ust,ZOL=zol,HOL=hol &
1372 ,PBL=pblh,PSIM=psim,PSIH=psih &
1373 ,XLAND=xland,HFX=hfx,QFX=qfx,TSK=tskold &
1374 ,GZ1OZ0=gz1oz0,WSPD=wspd,BR=br &
1375 ,DT=dtbl,DTMIN=dtmin,KPBL2D=kpbl &
1376 ,SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0 &
1377 ,EP1=ep_1,EP2=ep_2,KARMAN=karman,EOMEG=eomeg &
1378 ,STBOLT=stbolt,REGIME=regime &
1379 ,FLAG_QI=flag_qi &
1380 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1381 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1382 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1383 )
1384 ELSE
1385 WRITE ( message , FMT = '(A,5(L1,1X))' ) &
1386 'present: '// &
1387 'qv_curr, '// &
1388 'qc_curr, '// &
1389 'rqvblten, '// &
1390 'rqcblten, '// &
1391 'hol = ' , &
1392 PRESENT( qv_curr ) , &
1393 PRESENT( qc_curr ) , &
1394 PRESENT( rqvblten ) , &
1395 PRESENT( rqcblten ) , &
1396 PRESENT( hol )
1397 CALL wrf_debug(0,message)
1398 CALL wrf_error_fatal('Lack arguments to call MRF pbl')
1399 ENDIF
1401 CASE (GFSSCHEME)
1402 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1403 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1404 .TRUE. ) THEN
1405 CALL wrf_debug(100,'in GFS')
1406 CALL bl_gfs( &
1407 U3D=u_phytmp,V3D=v_phytmp &
1408 ,TH3D=th_phy,T3D=t_phy &
1409 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
1410 ,P3D=p_phy,PI3D=pi_phy &
1411 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1412 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1413 ,RQIBLTEN=rqiblten &
1414 ,CP=cp,G=g,ROVCP=rcp,R=r_d,ROVG=rovg &
1415 ,DZ8W=dz8w,z=z,PSFC=psfc &
1416 ,UST=ust,PBL=pblh,PSIM=psim,PSIH=psih &
1417 ,HFX=hfx,QFX=qfx,TSK=tskold,GZ1OZ0=gz1oz0 &
1418 ,WSPD=wspd,BR=br &
1419 ,DT=dtbl,KPBL2D=kpbl,EP1=ep_1,KARMAN=karman &
1420 ,P_QI=p_qi,P_FIRST_SCALAR=param_first_scalar &
1421 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1422 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1423 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1424 )
1425 ELSE
1426 WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1427 'present: '// &
1428 'qv_curr, '// &
1429 'qc_curr, '// &
1430 'rqvblten, '// &
1431 'rqcblten = ', &
1432 PRESENT( qv_curr ) , &
1433 PRESENT( qc_curr ) , &
1434 PRESENT( rqvblten ) , &
1435 PRESENT( rqcblten )
1436 CALL wrf_debug(0,message)
1437 CALL wrf_error_fatal('Lack arguments to call GFS pbl')
1438 ENDIF
1440 CASE (MYJPBLSCHEME)
1441 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1442 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1443 .TRUE. ) THEN
1445 CALL wrf_debug(100,'in MYJPBL')
1446 IF ( .not.flag_bep .and. idiff.ne.1) THEN
1447 CALL myjpbl(DT=dt,STEPBL=stepbl,HT=ht,DZ=dz8w &
1448 ,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1449 ,QV=qv_curr,QCW=qc_curr,QCI=qi_curr,QCS=qs_curr & !BSF
1450 ,QCR=qr_curr,QCG=qg_curr & !BSF
1451 ,U=u_phy,V=v_phy,RHO=rho &
1452 ,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1453 ,QZ0=qz0,UZ0=uz0,VZ0=vz0 &
1454 ,LOWLYR=lowlyr &
1455 ,XLAND=xland,SICE=xice,SNOW=snow &
1456 ,TKE_MYJ=tke_pbl,EXCH_H=exch_h,USTAR=ust,ZNT=znt &
1457 ,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1458 ,AKHS=akhs,AKMS=akms,ELFLX=lh,MIXHT=mixht &
1459 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1460 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1461 ,RQIBLTEN=rqiblten,RQSBLTEN=rqsblten & !BSF
1462 ,RQRBLTEN=rqrblten,RQGBLTEN=rqgblten & !BSF
1463 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1464 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1465 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1466 )
1467 ELSE !(SF_URBAN_PHYSICS.EQ.2)
1468 ! Bep changes begin
1470 CALL myjurb(IDIFF=idiff,FLAG_BEP=flag_bep &
1471 ,DT=dtbl,STEPBL=stepbl,HT=ht,DZ=dz8w &
1472 ,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1473 ,QV=qv_curr, CWM=qc_curr &
1474 ,U=u_phy,V=v_phy,RHO=rho &
1475 ,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1476 ,QZ0=qz0,UZ0=uz0,VZ0=vz0 &
1477 ,LOWLYR=lowlyr &
1478 ,XLAND=xland,SICE=xice,SNOW=snow &
1479 ,TKE_MYJ=tke_pbl,EXCH_H=exch_h,EXCH_M=exch_m &
1480 ,USTAR=ust,ZNT=znt &
1481 ,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1482 ,AKHS=akhs,AKMS=akms,ELFLX=lh &
1483 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1484 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1485 ! Multi-layer UCM
1486 ,FRC_URB2D=frc_urb2d &
1487 ,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1488 ,A_Q_BEP=a_q_bep &
1489 ,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1490 ,B_T_BEP=b_t_bep,B_Q_BEP=b_q_bep &
1491 ,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep &
1492 ,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep &
1493 ! UCM end
1494 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1495 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1496 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1497 )
1498 ENDIF
1500 ELSE
1501 WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1502 'present: '// &
1503 'qv_curr, '// &
1504 'qc_curr, '// &
1505 'rqvblten, '// &
1506 'rqcblten = ' , &
1507 PRESENT( qv_curr ) , &
1508 PRESENT( qc_curr ) , &
1509 PRESENT( rqvblten ) , &
1510 PRESENT( rqcblten )
1511 CALL wrf_debug(0,message)
1512 CALL wrf_error_fatal('Lack arguments to call MYJ pbl')
1513 ENDIF
1515 CASE (QNSEPBLSCHEME)
1516 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1517 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1518 .TRUE. ) THEN
1519 IF ( MFSHCONV.EQ.1 )THEN
1520 CALL wrf_debug(100,'in MFSHCONVPBL')
1521 CALL mfshconvpbl(DT=dt,STEPBL=stepbl,HT=ht,DZ=dz8w &
1522 ,RHO=rho,PMID=p_phy,PINT=p8w,TH=th_phy,EXNER=pi_phy &
1523 ,QV=qv_curr, QC=qc_curr &
1524 ,U=u_phy,V=v_phy &
1525 ,HFX=hfx, QFX=qfx,TKE=tke_pbl &
1526 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1527 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1528 ,IDS=IDS,IDE=IDE,JDS=JDS,JDE=JDE,KDS=KDS,KDE=KDE &
1529 ,IMS=IMS,IME=IME,JMS=JMS,JME=JME,KMS=KMS,KME=KME &
1530 ,ITS=ITS,ITE=ITE,JTS=JTS,JTE=JTE,KTS=KTS,KTE=KTE &
1531 ,KRR=2,MASSFLUX_EDKF=massflux_EDKF &
1532 ,ENTR_EDKF=entr_EDKF, DETR_EDKF=detr_EDKF &
1533 ,THL_UP=thl_up, THV_UP=thv_up, RT_UP=rt_up &
1534 ,RV_UP=rv_up,RC_UP=rc_up, U_UP=u_up, V_UP=v_up &
1535 ,FRAC_UP=frac_up, RC_MF=rc_mf,WTHV=u_phytmp &
1536 ,PLM_BL89=v_phytmp )
1537 ENDIF
1539 CALL wrf_debug(100,'in QNSEPBL')
1540 CALL qnsepbl( &
1541 DT=dt,STEPBL=stepbl,HT=ht,DZ=dz8w &
1542 ,PMID=p_phy,PINT=p8w,TH=th_phy,T=t_phy,EXNER=pi_phy &
1543 ,QV=qv_curr, CWM=qc_curr &
1544 ,U=u_phy,V=v_phy,RHO=rho &
1545 ,TSK=tsk,QSFC=qsfc,CHKLOWQ=chklowq,THZ0=thz0 &
1546 ,QZ0=qz0,UZ0=uz0,VZ0=vz0,CORF=f &
1547 ,LOWLYR=lowlyr &
1548 ,XLAND=xland,SICE=xice,SNOW=snow &
1549 ,TKE=tke_pbl,EXCH_H=exch_h,EXCH_M=exch_m,USTAR=ust,ZNT=znt &
1550 ,EL_MYJ=el_pbl,PBLH=pblh,KPBL=kpbl,CT=ct &
1551 ,AKHS=akhs,AKMS=akms,ELFLX=lh &
1552 ,HFX=hfx,WTHV_MF=u_phytmp,LM_BL89=v_phytmp &
1553 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1554 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1555 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1556 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1557 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1558 )
1559 ! note the tendencies have been added inside qnse
1560 ! IF ( PRESENT (MFSHCONV) )THEN
1561 ! IF (MFSHCONV.EQ.1)THEN
1562 ! rublten(:,:,:)=rublten(:,:,:)+rumften(:,:,:)
1563 ! rvblten(:,:,:)=rvblten(:,:,:)+rvmften(:,:,:)
1564 ! rthblten(:,:,:)=rthblten(:,:,:)+rthmften(:,:,:)
1565 ! rqvblten(:,:,:)=rqvblten(:,:,:)+rqvmften(:,:,:)
1566 ! rqcblten(:,:,:)=rqcblten(:,:,:)+rqcmften(:,:,:)
1567 ! ENDIF
1568 ! ENDIF
1570 ELSE
1571 WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1572 'present: '// &
1573 'qv_curr, '// &
1574 'qc_curr, '// &
1575 'rqvblten, '// &
1576 'rqcblten = ' , &
1577 PRESENT( qv_curr ) , &
1578 PRESENT( qc_curr ) , &
1579 PRESENT( rqvblten ) , &
1580 PRESENT( rqcblten )
1581 CALL wrf_debug(0,message)
1582 CALL wrf_error_fatal('Lack arguments to call QNSE pbl')
1583 ENDIF
1585 CASE (ACMPBLSCHEME)
1587 !! These are values that are not supplied to pbl driver, but are required by ACM
1588 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1589 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1590 .TRUE. ) THEN
1591 CALL wrf_debug(100,'in ACM PBL')
1593 CALL ACMPBL( &
1594 XTIME=itimestep, DTPBL=dtbl,U3D=u_phytmp, V3D=v_phytmp &
1595 ,PP3D=p_phy, DZ8W=dz8w, TH3D=th_phy, T3D=t_phy &
1596 ,QV3D=qv_curr, QC3D=qc_curr, QI3D=qi_curr, RR3D=rho &
1597 #if (WRF_CHEM == 1)
1598 ,CHEM3D=CHEM, VD3D=VD, NCHEM=nchem, KDVEL=kdvel &
1599 ,NDVEL=ndvel, NUM_VERT_MIX=num_vert_mix &
1600 #endif
1601 ,UST=UST, HFX=HFX, QFX=QFX, TSK=tsk &
1602 ,PSFC=PSFC, EP1=EP_1, G=g, ROVCP=rcp,RD=r_D,CPD=cp &
1603 ,PBLH=pblh, KPBL2D=kpbl, EXCH_H=exch_h, EXCH_M=exch_m &
1604 ,REGIME=regime &
1605 ,GZ1OZ0=gz1oz0,WSPD=wspd,PSIM=psim, MUT=mut, RMOL=rmol &
1606 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1607 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
1608 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1609 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1610 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1611 )
1612 ELSE
1613 WRITE ( message , FMT = '(A,4(L1,1X))' ) &
1614 'present: '// &
1615 'qv_curr, '// &
1616 'qc_curr, '// &
1617 'rqvblten, '// &
1618 'rqcblten = ' , &
1619 PRESENT( qv_curr ) , &
1620 PRESENT( qc_curr ) , &
1621 PRESENT( rqvblten ) , &
1622 PRESENT( rqcblten )
1623 CALL wrf_debug(0,message)
1624 CALL wrf_error_fatal('Lack arguments to call ACM2 pbl')
1625 ENDIF
1627 #if (EM_CORE==1)
1629 CASE (MYNNPBLSCHEME)
1631 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1632 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1633 PRESENT( qi_curr ) .AND. PRESENT( rqiblten ) .AND. &
1634 PRESENT( rqniblten ) .AND. PRESENT( qni_curr ) .AND.&
1635 PRESENT(qke) .AND. PRESENT(tsq) .AND. &
1636 PRESENT(qsq) .AND. PRESENT(cov) .AND. &
1637 PRESENT(rmol) .AND. PRESENT(ch) .AND. &
1638 PRESENT(tke_budget) .AND. PRESENT(qke_adv) .AND. &
1639 PRESENT(bl_mynn_tkeadvect) ) THEN
1641 CALL wrf_debug(100,'in MYNNPBL')
1643 IF (itimestep==1) THEN
1644 initflag=1
1645 ELSE
1646 initflag=0
1647 ENDIF
1649 if (pert_mynn .and. multi_perturb == 1) then
1650 allocate (qke_tmp(its:ite, kts:kte, jts:jte))
1652 call Add_multi_perturb_pbl_perturbations ( &
1653 perts_qvapor, perts_qcloud, perts_th, perts_tke, &
1654 pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke,&
1655 th_phy, qke, qke_adv, qv_curr, qc_curr, &
1656 bl_mynn_tkeadvect, qke_tmp, its, ite, jts, jte, &
1657 ims, ime, jms, jme, kms, kme, kts, kte)
1658 end if
1660 CALL mynnedmf_wrapper_run( &
1661 &initflag=initflag,restart=restart,cycling=cycling, &
1662 &delt=dtbl,dz=dz8w,dxc=dx,znt=znt, &
1663 &u=u_phy,v=v_phy,w=w,th=th_phy,qv=qv_curr, &
1664 &qc=qc_curr,qi=qi_curr,qs=qs_curr, &
1665 &qnc=qnc_curr,qni=qni_curr, &
1666 &QNWFA=qnwfa_curr,QNIFA=qnifa_curr,QNBCA=qnbca_curr, &
1667 ! &ozone=ozone, &
1668 &p=p_phy,exner=pi_phy,rho=rho,T3D=t_phy, &
1669 &xland=xland,ts=tsk,qsfc=qsfc,ps=psfc, &
1670 &ust=ust,ch=ch,hfx=hfx,qfx=qfx,rmol=rmol,wspd=wspd, &
1671 &uoce=uoce,voce=voce, & !Ocean currents
1672 &Qke=qke,qke_adv=qke_adv,Sh3d=Sh3d,Sm3d=Sm3d, &
1673 #if (WRF_CHEM == 1)
1674 &mix_chem=mynn_chem_vertmx, &
1675 &chem3d=chem,vd3d=vd,nchem=nchem,kdvel=kdvel, &
1676 &ndvel=ndvel,num_vert_mix=num_vert_mix, &
1677 ! &FRP_MEAN=FRP_MEAN,EMIS_ANT_NO=EMIS_ANT_NO, & !to be included soon
1678 ! &fire_turb=enh_vermix, & !to be included soon
1679 #endif
1680 &Tsq=tsq,Qsq=qsq,Cov=cov, &
1681 &RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten, &
1682 &RQVBLTEN=rqvblten,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten,&
1683 &RQNCBLTEN=rqncblten,RQNIBLTEN=rqniblten, &
1684 &RQSBLTEN=rqsblten, &
1685 &RQNWFABLTEN=rqnwfablten,RQNIFABLTEN=rqnifablten, &
1686 &RQNBCABLTEN=rqnbcablten, &
1687 ! &Ro3BLTEN=ro3blten, &
1688 &EXCH_H=exch_h,EXCH_M=exch_m, &
1689 &pblh=pblh,KPBL=KPBL, &
1690 &el_pbl=el_pbl, &
1691 &dqke=dqke, &
1692 &qWT=qWT,qSHEAR=qSHEAR,qBUOY=qBUOY,qDISS=qDISS, &
1693 &qc_bl=qc_bl,qi_bl=qi_bl,cldfra_bl=cldfra_bl, &
1694 &edmf_a=edmf_a,edmf_w=edmf_w,edmf_qt=edmf_qt, &
1695 &edmf_thl=edmf_thl,edmf_ent=edmf_ent,edmf_qc=edmf_qc, &
1696 &sub_thl3D=sub_thl3D,sub_sqv3D=sub_sqv3D, &
1697 &det_thl3D=det_thl3D,det_sqv3D=det_sqv3D, &
1698 &maxwidth=maxwidth,maxMF=maxMF, &
1699 &ztop_plume=ztop_plume,ktop_plume=ktop_plume, &
1700 &RTHRATEN=RTHRATEN, &
1701 &bl_mynn_tkeadvect=bl_mynn_tkeadvect, &
1702 &tke_budget=tke_budget, &
1703 &bl_mynn_cloudpdf=bl_mynn_cloudpdf, &
1704 &bl_mynn_mixlength=bl_mynn_mixlength, &
1705 &icloud_bl=icloud_bl, &
1706 &bl_mynn_edmf=bl_mynn_edmf, &
1707 &bl_mynn_edmf_mom=bl_mynn_edmf_mom, &
1708 &bl_mynn_edmf_tke=bl_mynn_edmf_tke, &
1709 &bl_mynn_mixscalars=bl_mynn_mixscalars, &
1710 &bl_mynn_output=bl_mynn_output, &
1711 &bl_mynn_cloudmix=bl_mynn_cloudmix, &
1712 &bl_mynn_mixqt=bl_mynn_mixqt, &
1713 &bl_mynn_closure=bl_mynn_closure, &
1714 &spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl, &
1715 &FLAG_QC=flag_qc,FLAG_QI=flag_qi,FLAG_QS=flag_qs, &
1716 &FLAG_QNC=flag_qnc,FLAG_QNI=flag_qni, &
1717 &FLAG_QNWFA=flag_qnwfa,FLAG_QNIFA=flag_qnifa, &
1718 &FLAG_QNBCA=flag_qnbca, &
1719 &IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde, &
1720 &IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme, &
1721 &ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
1723 if (pert_mynn .and. multi_perturb == 1) then
1724 call Remove_multi_perturb_pbl_perturbations ( &
1725 perts_qvapor, perts_qcloud, perts_th, perts_tke, &
1726 pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke,&
1727 th_phy, qke, qke_adv, qv_curr, qc_curr, &
1728 bl_mynn_tkeadvect, qke_tmp, its, ite, jts, jte, &
1729 ims, ime, jms, jme, kms, kme, kts, kte)
1731 deallocate (qke_tmp)
1732 end if
1733 ELSE
1734 WRITE ( message , FMT = '(A,17(L1,1X))' ) &
1735 'present: '// &
1736 'qv_curr, '// &
1737 'qc_curr, '// &
1738 'qi_curr, '// &
1739 'qni_curr, '// &
1740 'rqvblten, '// &
1741 'rqcblten, '// &
1742 'rqiblten, '// &
1743 'rqniblten, '// &
1744 'qke, '// &
1745 'tsq, '// &
1746 'qsq, '// &
1747 'cov, '// &
1748 'rmol, '// &
1749 'ch, '// &
1750 'tke_budget, '// &
1751 'qke_adv, '// &
1752 'bl_mynn_tkeadvect ' , &
1753 PRESENT( qv_curr ) , &
1754 PRESENT( qc_curr ) , &
1755 PRESENT( qi_curr ) , &
1756 PRESENT( qni_curr ) , &
1757 PRESENT( rqvblten ) , &
1758 PRESENT( rqcblten ) , &
1759 PRESENT( rqiblten ) , &
1760 PRESENT( rqniblten ) , &
1761 PRESENT( qke ) , &
1762 PRESENT( tsq ) , &
1763 PRESENT( qsq ) , &
1764 PRESENT( cov ) , &
1765 PRESENT( rmol ) , &
1766 PRESENT( ch ) , &
1767 PRESENT( tke_budget) , &
1768 PRESENT( qke_adv ) , &
1769 PRESENT( bl_mynn_tkeadvect)
1770 CALL wrf_debug(0,message)
1771 CALL wrf_error_fatal('Lack arguments to call MYNN pbl')
1772 ENDIF
1774 !fill scalar_tend array. The RQN*BLTEN arrays are no longer
1775 !passed to module_physics_addtendc.F (for MYNN)
1776 IF (bl_mynn_mixscalars .eq. 1) THEN
1777 IF (PRESENT( qni_curr ) .AND. Flag_qni) THEN
1778 !print*,"Updating qni after mynn-edmf",P_QNI
1779 DO j=j_start(ij),j_end(ij)
1780 DO k=kts,kte
1781 DO i=i_start(ij),i_end(ij)
1782 scalar_tend(i,k,j,P_QNI)=RQNIBLTEN(i,k,j)
1783 enddo
1784 enddo
1785 enddo
1786 ENDIF
1787 IF (PRESENT( qnc_curr ) .AND. Flag_qnc) THEN
1788 !print*,"Updating qnc after mynn-edmf",P_QNC
1789 DO j=j_start(ij),j_end(ij)
1790 DO k=kts,kte
1791 DO i=i_start(ij),i_end(ij)
1792 scalar_tend(i,k,j,P_QNC)=RQNCBLTEN(i,k,j)
1793 enddo
1794 enddo
1795 enddo
1796 ENDIF
1797 IF (PRESENT( qnwfa_curr ) .AND. Flag_qnwfa) THEN
1798 !print*,"Updating qnwfa after mynn-edmf",P_QNWFA
1799 DO j=j_start(ij),j_end(ij)
1800 DO k=kts,kte
1801 DO i=i_start(ij),i_end(ij)
1802 scalar_tend(i,k,j,P_QNWFA)=RQNWFABLTEN(i,k,j)
1803 enddo
1804 enddo
1805 enddo
1806 ENDIF
1807 IF (PRESENT( qnifa_curr ) .AND. Flag_qnifa) THEN
1808 !print*,"Updating qnifa after mynn-edmf",P_QNIFA
1809 DO j=j_start(ij),j_end(ij)
1810 DO k=kts,kte
1811 DO i=i_start(ij),i_end(ij)
1812 scalar_tend(i,k,j,P_QNIFA)=RQNIFABLTEN(i,k,j)
1813 enddo
1814 enddo
1815 enddo
1816 ENDIF
1817 IF (PRESENT( qnbca_curr ) .AND. Flag_qnbca) THEN
1818 !print*,"Updating qnbca after mynn-edmf",P_QNBCA
1819 DO j=j_start(ij),j_end(ij)
1820 DO k=kts,kte
1821 DO i=i_start(ij),i_end(ij)
1822 scalar_tend(i,k,j,P_QNBCA)=RQNBCABLTEN(i,k,j)
1823 enddo
1824 enddo
1825 enddo
1826 ENDIF
1827 ENDIF
1829 CASE (EEPSSCHEME)
1831 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1832 PRESENT( qi_curr ) .AND. &
1833 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1834 PRESENT( rqiblten ) .AND. &
1835 PRESENT( pek ) .AND. PRESENT( pep) .AND. &
1836 PRESENT( pek_adv ) .AND. PRESENT( pep_adv ) ) THEN
1838 CALL wrf_debug(100,'in EEPSPBL')
1840 CALL eeps(&
1841 U3D=u_phytmp,V3D=v_phytmp,T3D=t_phy &
1842 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
1843 ,QR3D=qr_curr,QS3D=qs_curr,QG3D=qg_curr &
1844 ,P3D=p_phy,PI3D=pi_phy,PEK=pek,PEP=pep &
1845 ,RHO3D=rho,RTHRATEN=rthraten &
1846 ,RUBLTEN=rublten,RVBLTEN=rvblten &
1847 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
1848 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten,KPBL=kpbl &
1849 ,DZ8W=dz8w,PSFC=psfc,TSK=tsk,QSFC=qsfc,WSPD=wspd &
1850 ,UST=ust,XLAND=xland,HFX=hfx,QFX=qfx,RMOL=rmol &
1851 ,DT=dtbl,DX=dx,EXCH_H=exch_h,EXCH_M=exch_m,PBLH=pblh &
1852 ,ITIMESTEP=itimestep,PEK_ADV=pek_adv,PEP_ADV=pep_adv &
1853 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1854 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1855 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
1856 )
1857 ELSE
1858 WRITE ( message , FMT = '(A,6(L1,1X))' ) &
1859 'present: '// &
1860 'qv_curr, '// &
1861 'qc_curr, '// &
1862 'qi_curr, '// &
1863 'rqvblten, '// &
1864 'rqcblten, '// &
1865 'rqiblten, ' , &
1866 PRESENT( qv_curr ) , &
1867 PRESENT( qc_curr ) , &
1868 PRESENT( qi_curr ) , &
1869 PRESENT( rqvblten ) , &
1870 PRESENT( rqcblten ) , &
1871 PRESENT( rqiblten )
1872 CALL wrf_debug(0,message)
1873 CALL wrf_error_fatal('Lack arguments to call EEPS pbl')
1874 ENDIF
1877 CASE (KEPSSCHEME)
1879 CALL wrf_debug(100,'in keps')
1880 CALL KEPS(TSK=tsk,MOL=mol,XTIME=itimestep,FRC_URB2D=frc_urb2d &
1881 , FLAG_BEP=flag_bep &
1882 ,DZ8W=dz8w,DT=dt,U_PHY=u_phytmp &
1883 ,V_PHY=v_phytmp,TH_PHY=th_phy &
1884 ,PI_PHY=pi_phy,RTHRATEN=RTHRATEN,P8W=p8w &
1885 ,RHO=rho,QV_CURR=qv_curr,QC_CURR=qc_curr,HFX=hfx &
1886 ,QFX=qfx,USTAR=ust,CP=cp,G=g &
1887 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1888 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1889 ,TKE_PBL=tke_pbl,DISS_PBL=diss_pbl,TPE_PBL=tpe_pbl &
1890 ,TKE_ADV=tke_adv,DISS_ADV=diss_adv,TPE_ADV=tpe_adv &
1891 ,PR_PBL=pr_pbl &
1892 ,WU=wu_tur,WV=wv_tur,WT=wt_tur,WQ=wq_tur &
1893 ,EXCH_H=exch_h,EXCH_M=exch_m,PBLH=pblh &
1894 ,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1895 ,A_Q_BEP=a_q_bep &
1896 ,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1897 ,B_T_BEP=b_t_bep &
1898 ,B_Q_BEP=b_q_bep &
1899 ,B_E_BEP=b_e_bep &
1900 ,SF_BEP=sf_bep,VL_BEP=vl_bep &
1901 ,BR=br,ZNT=znt &
1902 ,PSIM=psim,PSIH=psih &
1903 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1904 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1905 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
1909 CASE (BOULACSCHEME)
1911 CALL wrf_debug(100,'in boulac')
1913 CALL BOULAC(FRC_URB2D=frc_urb2d,IDIFF=idiff,FLAG_BEP=flag_bep &
1914 ,DZ8W=dz8w,DT=dt,U_PHY=u_phytmp &
1915 ,V_PHY=v_phytmp,TH_PHY=th_phy &
1916 ,RHO=rho,QV_CURR=qv_curr,QC_CURR=qc_curr,HFX=hfx &
1917 ,QFX=qfx,USTAR=ust,CP=cp,G=g &
1918 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1919 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
1920 ,TKE=tke_pbl,DLK=el_pbl,WU=wu_tur,WV=wv_tur,WT=wt_tur,WQ=wq_tur &
1921 ,EXCH_H=exch_h,EXCH_M=exch_m,PBLH=pblh &
1922 ,A_U_BEP=a_u_bep,A_V_BEP=a_v_bep,A_T_BEP=a_t_bep &
1923 ,A_Q_BEP=a_q_bep &
1924 ,A_E_BEP=a_e_bep,B_U_BEP=b_u_bep,B_V_BEP=b_v_bep &
1925 ,B_T_BEP=b_t_bep,B_E_BEP=b_e_bep,DLG_BEP=dlg_bep &
1926 ,B_Q_BEP=b_q_bep &
1927 ,DL_U_BEP=dl_u_bep,SF_BEP=sf_bep,VL_BEP=vl_bep &
1928 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
1929 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
1930 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
1932 CASE (CAMUWPBLSCHEME)
1934 CALL wrf_debug(100,'in camuwpbl')
1935 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND. &
1936 PRESENT( qi_curr ) .AND. PRESENT( qnc_curr ) .AND. &
1937 PRESENT( qni_curr ) .AND. &
1938 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1939 PRESENT( rqiblten ) .AND. PRESENT( rqniblten ) &
1940 )THEN
1941 CALL CAMUWPBL(DT=dt,U_PHY=u_phy,V_PHY=v_phy,TH_PHY=th_phy,RHO=rho &
1942 ,QV_CURR=qv_curr,HFX=hfx,QFX=qfx,USTAR=ust,P8W=p8w,P_PHY=p_phy &
1943 ,Z=z,T_PHY=t_phy,QC_CURR=qc_curr,QI_CURR=qi_curr,Z_AT_W=z_at_w &
1944 ,CLDFRA_OLD_mp=cldfra_old_mp,CLDFRA=cldfra,HT=ht &
1945 ,RTHRATENLW=rthratenlw,EXNER=pi_phy &
1946 ,is_CAMMGMP_used=is_CAMMGMP_used &
1947 ,ITIMESTEP=itimestep,QNC_CURR=qnc_curr,QNI_CURR=qni_curr &
1948 ,WSEDL3D=wsedl3d &
1950 ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde &
1951 ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme &
1952 ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte &
1953 !Intent-inout
1954 ,TAURESX2D=tauresx2d,TAURESY2D=tauresy2d &
1955 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
1956 ,RQIBLTEN=rqiblten,RQNIBLTEN=rqniblten,RQVBLTEN=rqvblten &
1957 ,RQCBLTEN=rqcblten,KVM3D=exch_m,KVH3D=exch_h &
1958 !Intent-out
1959 ,TPERT2D=tpert2d,QPERT2D=qpert2d,WPERT2D=wpert2d,SMAW3D=smaw3d &
1960 ,TURBTYPE3D=turbtype3d &
1961 ,TKE_pbl=tke_pbl,PBLH2D=pblh,KPBL2D=kpbl )
1962 ELSE
1963 WRITE ( message , FMT = '(A,8(L1,1X))' ) &
1964 'present: '// &
1965 'qv_curr, '// &
1966 'qc_curr, '// &
1967 'qi_curr, '// &
1968 'qnc_curr, '// &
1969 'qni_curr, '// &
1970 'rqvblten, '// &
1971 'rqcblten, '// &
1972 'rqiblten, '// &
1973 'rqniblten= ', &
1974 PRESENT( qv_curr ) , &
1975 PRESENT( qc_curr ) , &
1976 PRESENT( qi_curr ) , &
1977 PRESENT( qnc_curr ) , &
1978 PRESENT( qni_curr ) , &
1979 PRESENT( rqvblten ) , &
1980 PRESENT( rqcblten ) , &
1981 PRESENT( rqiblten ) , &
1982 PRESENT( rqniblten )
1984 CALL wrf_debug(0,message)
1985 CALL wrf_error_fatal('Lack arguments to call CAMUWPBL pbl')
1986 ENDIF
1988 #endif
1990 CASE (GBMPBLSCHEME)
1991 CALL wrf_debug(100,'in gbmpbl')
1992 IF ( PRESENT( qv_curr ) .AND. PRESENT( qc_curr ) .AND.&
1993 PRESENT( qi_curr ) .AND. &
1994 PRESENT( rqvblten ) .AND. PRESENT( rqcblten ) .AND. &
1995 PRESENT( rqiblten ) .AND. &
1996 PRESENT( hol ) .AND. &
1997 .TRUE. ) THEN
1998 CALL gbmpbl( &
1999 U3D=u_phytmp,V3D=v_phytmp,TH3D=th_phy,T3D=t_phy &
2000 ,QV3D=qv_curr,QC3D=qc_curr,QI3D=qi_curr &
2001 ,P3D=p_phy,PI3D=pi_phy &
2002 ,RUBLTEN=rublten,RVBLTEN=rvblten &
2003 ,RTHBLTEN=rthblten,RQVBLTEN=rqvblten &
2004 ,RQCBLTEN=rqcblten,RQIBLTEN=rqiblten &
2005 ,KZM_GBM=exch_m,KTH_GBM=exch_h,KETHL_GBM=exch_tke &
2006 ,EL_GBM=el_pbl,DZ8W=dz8w,Z=z,PSFC=PSFC &
2007 ,TKE_PBL=tke_pbl,RTHRATEN=rthraten &
2008 ,ZNT=znt,UST=ust,ZOL=zol,HOL=hol,PBL=pblh &
2009 ,KPBL2D=kpbl,REGIME=regime &
2010 ,PSIM=psim,PSIH=psih,XLAND=xland &
2011 ,HFX=hfx,QFX=qfx,TSK=tskold,GZ1OZ0=gz1oz0 &
2012 ,WSPD=wspd,BR=br,DT=dtbl,DTMIN=dtmin &
2013 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2014 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2015 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
2016 )
2017 ELSE
2018 CALL wrf_error_fatal('Lack arguments to call GBM pbl')
2019 ENDIF
2021 #if ( WRFPLUS == 1 )
2022 ! this scheme is for WRFPlus only
2023 !-----------------------------------
2024 CASE (SURFDRAGSCHEME)
2026 CALL wrf_debug(100,'in SURFDRAG scheme')
2028 CALL surface_drag( &
2029 RUBLTEN=rublten,RVBLTEN=rvblten & !
2030 ,U_PHY=u_phy, V_PHY=v_phy, Z=z & !
2031 ,XLAND=xland, HT=ht, KPBL2D=kpbl & !
2032 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2033 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2034 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
2035 )
2036 #endif
2038 CASE DEFAULT
2040 WRITE( message , * ) 'The pbl option does not exist: bl_pbl_physics = ', bl_pbl_physics
2041 CALL wrf_error_fatal ( message )
2043 END SELECT pbl_select
2045 #if (EM_CORE==1)
2046 !
2047 ! ... paj: wind farm ...
2048 !
2049 windfarm_select: SELECT CASE(windfarm_opt)
2051 CASE (fitchscheme)
2053 IF (PRESENT(id) .AND. &
2054 PRESENT(z_at_w) ) THEN
2055 !
2056 CALL wrf_debug(100,'in phys/module_wind_fitch.F')
2057 CALL dragforce( &
2058 & ID=id &
2059 &,Z_AT_W=z_at_w,u=u_phy,v=v_phy &
2060 &,DX=dx,DZ=dz8w,DT=dt &
2061 &,QKE=qke &
2062 &,DU=rublten,DV=rvblten &
2063 &,WINDFARM_OPT=windfarm_opt,POWER=power &
2064 &,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2065 &,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2066 &,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
2067 &)
2069 IF (bl_mynn_tkeadvect) THEN
2070 qke_adv=qke
2071 ENDIF
2073 ELSE
2074 WRITE ( message , FMT = '(A,6(L1,1X))' ) &
2075 'present: '// &
2076 'ID, '// &
2077 'z_at_w, '// &
2078 'xlat_u, '// &
2079 'xlong_u, '// &
2080 'xlat_v, '// &
2081 'xlong_v = ' , &
2082 PRESENT( id ) , &
2083 PRESENT( z_at_w )
2084 CALL wrf_debug(0,message)
2085 CALL wrf_error_fatal('Lack arguments to call turbine_drag')
2086 ENDIF
2088 ! Yulong add new wind farm schemes with wind turbine loss effect
2089 CASE (mavscheme)
2090 IF (PRESENT(id) .AND. &
2091 PRESENT(z_at_w) ) THEN
2092 CALL wrf_debug(100,'in phys/module_wind_mav.F')
2093 CALL dragforce_mav(itimestep &
2094 &,ID=id &
2095 &,Z_AT_W=z_at_w,z_at_m=z,u=u_phy,v=v_phy &
2096 &,DX=dx,DZ=dz8w,DT=dt &
2097 &,TKE=TKE_windfarm &
2098 &,DU=rublten,DV=rvblten &
2099 &,WINDFARM_OPT=windfarm_opt,POWER=power &
2100 &,windfarm_wake_model=windfarm_wake_model &
2101 &,windfarm_overlap_method=windfarm_overlap_method &
2102 &,xland=xland &
2103 &,cosa=cosa,sina=sina &
2104 &,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2105 &,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2106 &,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte &
2107 &)
2109 IF (bl_mynn_tkeadvect) THEN
2110 QKE = QKE + 2.*TKE_windfarm
2111 qke_adv=qke
2112 ENDIF
2114 ELSE
2115 WRITE ( message , FMT = '(A,6(L1,1X))' ) &
2116 'present: '// &
2117 'ID, '// &
2118 'z_at_w, '// &
2119 'xlat_u, '// &
2120 'xlong_u, '// &
2121 'xlat_v, '// &
2122 'xlong_v = ' , &
2123 PRESENT( id ) , &
2124 PRESENT( z_at_w )
2125 CALL wrf_debug(0,message)
2126 CALL wrf_error_fatal('Lack arguments to call dragforce_mav')
2127 ENDIF
2129 END SELECT windfarm_select
2130 #endif
2133 IF (PRESENT(GWD_opt)) THEN
2134 IF (GWD_opt .EQ. 1) THEN
2135 IF (PRESENT(dtaux3d)) THEN
2136 CALL gwdo( &
2137 U3D=u_phy,V3D=v_phy,T3D=t_phy &
2138 ,QV3D=qv_curr &
2139 ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy,Z=z &
2140 ,RUBLTEN=rublten,RVBLTEN=rvblten &
2141 ,DTAUX3D=dtaux3d,DTAUY3D=dtauy3d &
2142 ,DUSFCG=dusfcg,DVSFCG=dvsfcg &
2143 ,VAR2D=var2d,OC12D=oc12d &
2144 ,OA2D1=oa1,OA2D2=oa2,OA2D3=oa3,OA2D4=oa4 &
2145 ,OL2D1=ol1,OL2D2=ol2,OL2D3=ol3,OL2D4=ol4 &
2146 ,SINA=sina,COSA=cosa &
2147 ,ZNU=znu,ZNW=znw,P_TOP=p_top &
2148 ,CP=cp,G=g,RD=r_d &
2149 ,RV=r_v,EP1=ep_1,PI=3.141592653 &
2150 ,DT=dtbl,DX=dx2d,KPBL2D=kpbl,ITIMESTEP=itimestep &
2151 ,errmsg=errmsg,errflg=errflg &
2152 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2153 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2154 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
2155 ENDIF
2156 ELSEIF(gwd_opt .EQ. 3) THEN
2157 CALL gwdo_gsl( &
2158 U3D=u_phy,V3D=v_phy,T3D=t_phy &
2159 ,QV3D=qv_curr &
2160 ,P3D=p_phy,P3DI=p8w,PI3D=pi_phy,Z=z &
2161 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
2162 ,DTAUX3D_ls=dtaux3d_ls,DTAUY3D_ls=dtauy3d_ls &
2163 ,DTAUX3D_bl=dtaux3d_bl,DTAUY3D_bl=dtauy3d_bl &
2164 ,DTAUX3D_ss=dtaux3d_ss,DTAUY3D_ss=dtauy3d_ss &
2165 ,DTAUX3D_fd=dtaux3d_fd,DTAUY3D_fd=dtauy3d_fd &
2166 ,DUSFCG_ls=dusfcg_ls,DVSFCG_ls=dvsfcg_ls &
2167 ,DUSFCG_bl=dusfcg_bl,DVSFCG_bl=dvsfcg_bl &
2168 ,DUSFCG_ss=dusfcg_ss,DVSFCG_ss=dvsfcg_ss &
2169 ,DUSFCG_fd=dusfcg_fd,DVSFCG_fd=dvsfcg_fd &
2170 ,xland=xland,br=br &
2171 ,VAR2D=var2dls,OC12D=oc12dls &
2172 ,OA2D1=oa1ls,OA2D2=oa2ls,OA2D3=oa3ls,OA2D4=oa4ls &
2173 ,OL2D1=ol1ls,OL2D2=ol2ls,OL2D3=ol3ls,OL2D4=ol4ls &
2174 ,VAR2Dss=var2dss,OC12Dss=oc12dss &
2175 ,OA2D1ss=oa1ss,OA2D2ss=oa2ss,OA2D3ss=oa3ss &
2176 ,OA2D4ss=oa4ss,OL2D1ss=ol1ss,OL2D2ss=ol2ss &
2177 ,OL2D3ss=ol3ss,OL2D4ss=ol4ss &
2178 ,SINA=sina,COSA=cosa &
2179 ,ZNU=znu,ZNW=znw,P_TOP=p_top &
2180 ,DZ=dz8w,PBLH=PBLH &
2181 ,CP=cp,G=g,RD=r_d &
2182 ,RV=r_v,EP1=ep_1,PI=3.141592653 &
2183 ,DT=dtbl,DX=dx,KPBL2D=kpbl &
2184 ,ITIMESTEP=itimestep,gwd_opt=gwd_opt &
2185 ,gwd_diags=gwd_diags &
2186 ,spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl &
2187 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2188 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2189 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte )
2190 ENDIF
2191 ENDIF
2193 #if (EM_CORE==1)
2194 !JOE-added - fog (cloud) water deposition calculation
2195 IF (grav_settling .GE. 1) THEN
2196 IF ( PRESENT(vdfg) .AND. PRESENT(qc_curr)) THEN
2197 !NOTE 1: vdfg is calculated in the surface driver.
2198 !NOTE 2: as of now, only qc_curr settles, not # conc.
2200 CALL bl_fogdes( &
2201 vdfg,qc_curr,dtbl,rho,dz8w,grav_settling,RQCBLTEN,&
2202 ids,ide, jds,jde, kds,kde, &
2203 ims,ime, jms,jme, kms,kme, &
2204 i_start(ij),i_end(ij), &
2205 j_start(ij),j_end(ij),kts,kte )
2206 ELSE
2207 CALL wrf_error_fatal('Missing args for bl_fogdes in pbl driver')
2208 ENDIF
2209 ENDIF
2210 !JOE-END
2211 #endif
2213 IF (idiff.eq.1) THEN
2215 !Alberto: here we call the general routine to solve the diffusion
2216 ! + all the source/sink terms.
2217 ! the only thing that should be passed from the PBL schemes is the value of the exch_h, and exch_m
2218 ! (and the non local terms, if present, through the b).
2219 ! So, this routine can be used by any of the previous schemes. We only need to pass the right variables
2220 ! (and eliminate the computation of the tendencies from the previous routines, to avoid doing twice the job).
2221 ! As I explain below, in the routine, here we could extract the vertical turbulent
2222 ! fluxes (something that will be general for all the routines).
2223 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2226 CALL diff3d (DT=dtbl,CP=cp,DZ=dz8w,TH=th_phy,QV=qv_curr,QC=qc_curr,T=t_phy &
2227 ,U=u_phy,V=v_phy,RHO=rho,EXCH_H=exch_h &
2228 ,EXCH_M=exch_m &
2229 ,RUBLTEN=rublten,RVBLTEN=rvblten,RTHBLTEN=rthblten &
2230 ,WU=wu_tur,WV=wv_tur,WT=wt_tur,WQ=wq_tur &
2231 ,RQVBLTEN=rqvblten,RQCBLTEN=rqcblten &
2232 ,A_U=a_u,A_V=a_v,A_T=a_t,A_Q=a_q &
2233 ,B_U=b_u,B_V=b_v,B_T=b_t,B_Q=b_q &
2234 ,SF=sf,VL=vl &
2235 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2236 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2237 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte)
2239 DEALLOCATE (a_u) ! Implicit component for the momemtum in X-direction
2240 DEALLOCATE (a_v) ! Implicit component for the momemtum in Y-direction
2241 DEALLOCATE (a_t) ! Implicit component for the Pot. Temp.
2242 DEALLOCATE (a_q) ! Implicit component for the water vapor
2243 DEALLOCATE (b_u) ! Explicit component for the momemtum in X-direction
2244 DEALLOCATE (b_v) ! Explicit component for the momemtum in Y-direction
2245 DEALLOCATE (b_t) ! Explicit component for the Pot. Temp.
2246 DEALLOCATE (b_q) ! Explicit component for the water vapor
2247 DEALLOCATE (sf ) ! surfaces
2248 DEALLOCATE (vl ) ! volumes
2249 ENDIF !idiff
2251 IF(scalar_pblmix .GT. 0)THEN
2252 CALL diff4d (DT=dtbl,DZ=dz8w, SCALAR=scalar, is_scalar=.true. &
2253 ,RHO=rho,EXCH_H=exch_h &
2254 ,EXCH_M=exch_m &
2255 ,SCALAR_TEND=scalar_tend &
2256 ,NUM_SCALAR=num_scalar, PARAM_FIRST_SCALAR=param_first_scalar &
2257 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2258 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2259 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte)
2260 ENDIF
2261 IF(tracer_pblmix .GT. 0)THEN
2262 CALL diff4d (DT=dtbl,DZ=dz8w, SCALAR=tracer, is_scalar=.false. &
2263 ,RHO=rho,EXCH_H=exch_h &
2264 ,EXCH_M=exch_m &
2265 ,SCALAR_TEND=tracer_tend &
2266 ,NUM_SCALAR=num_tracer, PARAM_FIRST_SCALAR=param_first_scalar &
2267 ,IDS=ids,IDE=ide,JDS=jds,JDE=jde,KDS=kds,KDE=kde &
2268 ,IMS=ims,IME=ime,JMS=jms,JME=jme,KMS=kms,KME=kme &
2269 ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte)
2270 ENDIF
2272 IF (distributed_ahe_opt == 1) THEN
2273 call cal_mon_day(julday, julyr, jmonth, jday)
2274 ihour = (jmonth - 1) * 24 + MOD(INT(gmt + xtime / 60.0), 24)
2275 DO j = jts, jte
2276 DO i = its, ite
2277 ! Volumetric heat capacity of air = 1200 J/(K m3)
2278 RTHBLTEN(i, 1, j) = RTHBLTEN(i, 1, j) + ahe(i, ihour, j) / 1200 / DZ8W(i, 1, j)
2279 END DO
2280 END DO
2281 END IF
2283 ENDDO
2284 !$OMP END PARALLEL DO
2286 ENDIF
2288 END SUBROUTINE pbl_driver
2290 subroutine Add_multi_perturb_pbl_perturbations (perts_qvapor, perts_qcloud, &
2291 perts_th, perts_tke, pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke, &
2292 th_phy, qke, qke_adv, qv_curr, qc_curr, bl_mynn_tkeadvect, qke_tmp, its, ite, jts, jte, &
2293 ims, ime, jms, jme, kms, kme, kts, kte)
2296 implicit none
2298 integer, intent(in) :: its, ite, jts, jte, ims, ime, jms, jme, kms, kme, kts, kte
2299 logical, optional, intent(in) :: bl_mynn_tkeadvect
2300 real, intent(in) :: pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke
2301 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_qvapor, &
2302 perts_qcloud, perts_th, perts_tke
2303 real, dimension(ims:ime, kms:kme, jms:jme), intent (inout) :: th_phy, qke, qv_curr, qc_curr
2304 real, optional, dimension (ims:ime, kms:kme, jms:jme), intent(inout) :: qke_adv
2305 real, dimension (its:ite, kts:kte, jts:jte), intent(out) :: qke_tmp
2307 integer :: i, j, k
2310 do j = jts, jte
2311 do k = kts, kte
2312 do i = its, ite
2313 qc_curr(i, k, j) = max (QCLOUD_MIN, (1.0 + perts_qcloud(i, k, j) * pert_mynn_qc) * qc_curr(i, k, j))
2314 qv_curr(i, k, j) = max (QVAPOR_MIN, (1.0 + perts_qvapor(i, k, j) * pert_mynn_qv) * qv_curr(i, k, j))
2315 th_phy(i, k, j) = (1.0 + perts_th(i, k, j) * pert_mynn_t) * th_phy(i, k, j)
2316 end do
2317 end do
2318 end do
2320 if (bl_mynn_tkeadvect) then
2321 do j = jts, jte
2322 do k = kts, kte
2323 do i = its, ite
2324 qke_tmp(i, k, j) = qke_adv(i, k, j)
2325 qke_adv(i, k, j) = max (QKE_MIN, (1.0 + perts_tke(i, k, j) * pert_mynn_qke) * qke_adv(i, k, j))
2326 end do
2327 end do
2328 end do
2329 else
2330 do j = jts, jte
2331 do k = kts, kte
2332 do i = its, ite
2333 qke_tmp(i, k, j) = qke(i, k, j)
2334 qke(i, k, j) = max (QKE_MIN, (1.0 + perts_tke(i, k, j) * pert_mynn_qke) * qke(i, k, j))
2335 end do
2336 end do
2337 end do
2338 end if
2340 end subroutine Add_multi_perturb_pbl_perturbations
2342 subroutine Remove_multi_perturb_pbl_perturbations (perts_qvapor, perts_qcloud, &
2343 perts_th, perts_tke, pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke, &
2344 th_phy, qke, qke_adv, qv_curr, qc_curr, bl_mynn_tkeadvect, qke_tmp, its, ite, jts, jte, &
2345 ims, ime, jms, jme, kms, kme, kts, kte)
2348 implicit none
2350 integer, intent(in) :: its, ite, jts, jte, ims, ime, jms, jme, kms, kme, kts, kte
2351 logical, optional, intent(in) :: bl_mynn_tkeadvect
2352 real, intent(in) :: pert_mynn_qv, pert_mynn_qc, pert_mynn_t, pert_mynn_qke
2353 real, dimension(ims:ime, kms:kme, jms:jme), optional, intent (in) :: perts_qvapor, &
2354 perts_qcloud, perts_th, perts_tke
2355 real, dimension(ims:ime, kms:kme, jms:jme), intent (inout) :: th_phy, qke, qv_curr, qc_curr
2356 real, optional, dimension (ims:ime, kms:kme, jms:jme), intent(inout) :: qke_adv
2357 real, dimension (its:ite, kts:kte, jts:jte), intent(in) :: qke_tmp
2359 integer :: i, j, k
2362 do j = jts, jte
2363 do k = kts, kte
2364 do i = its, ite
2365 qc_curr(i, k, j) = max (QCLOUD_MIN, qc_curr(i, k, j) / (1.0 + perts_qcloud(i, k, j) * pert_mynn_qc))
2366 qv_curr(i, k, j) = max (QVAPOR_MIN, qv_curr(i, k, j) / (1.0 + perts_qvapor(i, k, j) * pert_mynn_qv))
2367 th_phy(i, k, j) = th_phy(i, k, j) / (1.0 + perts_th(i, k, j) * pert_mynn_t)
2368 end do
2369 end do
2370 end do
2372 if (bl_mynn_tkeadvect) then
2373 do j = jts, jte
2374 do k = kts, kte
2375 do i = its, ite
2376 qke_adv(i, k, j) = max (QKE_MIN, qke_adv(i, k, j) - perts_tke(i, k, j) * pert_mynn_qke * qke_tmp(i, k, j))
2377 qke(i, k, j) = qke_adv(i, k, j)
2378 end do
2379 end do
2380 end do
2381 else
2382 do j = jts, jte
2383 do k = kts, kte
2384 do i = its, ite
2385 qke(i, k, j) = max (QKE_MIN, qke(i, k, j) - perts_tke(i, k, j) * pert_mynn_qke * qke_tmp(i, k, j))
2386 end do
2387 end do
2388 end do
2389 end if
2391 end subroutine Remove_multi_perturb_pbl_perturbations
2393 !=============================================================================
2394 SUBROUTINE diff3d(DT,CP,DZ,TH ,QV,QC,T,U,V,RHO &
2395 ,EXCH_H,EXCH_M &
2396 ,RUBLTEN,RVBLTEN,RTHBLTEN &
2397 ,RQVBLTEN,RQCBLTEN &
2398 ,WU,WV,WT,WQ &
2399 ,A_U,A_V,A_T,A_Q &
2400 ,B_U,B_V,B_T,B_Q &
2401 ,SF,VL &
2402 ,IDS,IDE,JDS,JDE,KDS,KDE &
2403 ,IMS,IME,JMS,JME,KMS,KME &
2404 ,ITS,ITE,JTS,JTE,KTS,KTE &
2405 )
2406 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2407 ! Subroutine written by Alberto Martilli, CIEMAT, Spain,
2408 ! e-mail:alberto_martilli@ciemat.es
2409 ! August 2008.
2410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2411 ! ALberto
2412 ! This routine solves the vertical diffusion
2413 ! + other source/sink terms (surface fluxes, building drag, non local terms, etc.)
2414 ! for U,V, potential temperature and moisture. The equation should be written
2415 ! as follow, for a generic variable M:
2416 !
2417 ! dM 1 d dM
2418 ! ---- =---- ---(rho*K----)+AM+B
2419 ! dt rho dz dz
2420 ! where A and B are the implict and explcit coefficients of the source/sink terms
2421 ! (at the lowest level the surface fluxes should go in these terms)
2422 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2423 !-----------------------------------------------------------------------
2425 IMPLICIT NONE
2426 !
2427 !----------------------------------------------------------------------
2428 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2429 & ,IMS,IME,JMS,JME,KMS,KME &
2430 & ,ITS,ITE,JTS,JTE,KTS,KTE
2432 ! INputs
2433 real DT,CP
2434 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ ! Depth of vertical levels
2435 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: TH ! Potential Temperature
2436 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QV ! water vapor mixing ratio (kg/kg)
2437 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QC ! cloud water mixing ratio (kg/kg)
2438 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: T ! temperature
2439 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: U ! X-compnent of the wind
2440 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: V ! Y-compnent of the wind
2441 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: RHO ! Air Density
2442 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_H ! Turbulent coefficient for heat
2443 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_M ! Turbulent coefficient for momentum
2445 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_U ! Implicit coefficient for DRAG (x -component)
2446 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_U ! Explicit coefficient for DRAG (x -component)
2447 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_V ! Implicit coefficient for DRAG (y -component)
2448 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_V ! Explicit coefficient for DRAG (y -component)
2449 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_T ! Implicit coefficient for Heat flux from buildings
2450 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_T ! Explicit coefficient for Heat flux from buildings
2451 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_Q ! Implicit coefficient for moisture flux
2452 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: B_Q ! Explicit coefficient for moisture flux
2454 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: VL ! fraction of air volume in the grid cell
2455 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: SF ! fraction of air at the face of the grid cell
2456 !outputs
2457 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RUBLTEN ! tendency for x-wind component
2458 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RVBLTEN ! tendency for y-wind component
2459 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RTHBLTEN ! tendency for Potential Temperature
2460 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RQVBLTEN ! tendency for water vapor mixing ratio (kg/kg)
2461 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: RQCBLTEN ! tendency for cloud water mixing ratio (kg/kg)
2462 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WU ! turbulent flux of momentum (x)
2463 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WV ! turbulent flux of momentum (y)
2464 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WT ! turbulent flux of potential temperature
2465 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: WQ ! turbulent flux of water vapor
2466 ! Parameters
2467 REAL ELOCP
2468 ! locals (same meaning as above, but 1D)
2470 real u1d(kms:kme),v1d(kms:kme),exch_h1d(kms:kme)
2471 real the1d(kms:kme) ! Equivalent potential temperature
2472 real exch_m1d(kms:kme),qv1d(kms:kme),qc1d(kms:kme)
2473 real dz1d(kms:kme),rho1d(kms:kme),rhoz1d(kms:kme)
2474 real sf1d(kms:kme),vl1d(kms:kme)
2475 real a_u1d(kms:kme),b_u1d(kms:kme)
2476 real a_v1d(kms:kme),b_v1d(kms:kme)
2477 real a_t1d(kms:kme),b_t1d(kms:kme)
2478 real a_q1d(kms:kme),b_q1d(kms:kme)
2479 real a_qc1d(kms:kme),b_qc1d(kms:kme)
2480 real wu1d(kms:kme),wv1d(kms:kme),wt1d(kms:kme),wq1d(kms:kme),wqc1d(kms:kme)
2481 real thnew
2482 !
2483 integer i,k,j
2484 !----------------------------------------------------------------------
2485 ELOCP=2.72E6/CP
2486 u1d=0.
2487 v1d=0.
2488 exch_h1d=0.
2489 exch_m1d=0.
2490 qv1d=0.
2491 qc1d=0.
2492 dz1d=0.
2493 rho1d=0.
2494 rhoz1d=0.
2495 sf1d=0.
2496 vl1d=0.
2497 a_u1d=0.
2498 a_v1d=0.
2499 a_t1d=0.
2500 a_q1d=0.
2501 a_qc1d=0.
2502 b_u1d=0.
2503 b_v1d=0.
2504 b_t1d=0.
2505 b_q1d=0.
2506 b_qc1d=0.
2508 do j=jts,jte
2509 do i=its,ite
2511 ! put three D variables in temporary 1D variables
2513 do k=kts,kte
2514 u1d(k)=U(i,k,j)
2515 v1d(k)=V(i,k,j)
2516 the1d(k)=TH(i,k,j)*(QC(i,k,j)*(-ELOCP/T(i,k,j))+1)
2517 qv1d(k)=qv(i,k,j)
2518 dz1d(k)=dz(i,k,j)
2519 rho1d(k)=rho(i,k,j)
2520 a_u1d(k)=a_u(i,k,j)
2521 b_u1d(k)=b_u(i,k,j)
2522 a_v1d(k)=a_v(i,k,j)
2523 b_v1d(k)=b_v(i,k,j)
2524 a_t1d(k)=a_t(i,k,j)
2525 b_t1d(k)=b_t(i,k,j)
2526 a_q1d(k)=a_q(i,k,j)
2527 b_q1d(k)=b_q(i,k,j)
2528 a_qc1d(k)=0.
2529 b_qc1d(k)=0.
2530 vl1d(k)=vl(i,k,j)
2531 sf1d(k)=sf(i,k,j)
2532 enddo
2533 sf1d(kte+1)=1.
2534 do k=kts,kte
2535 exch_h1d(k)=exch_h(i,k,j)
2536 exch_m1d(k)=exch_m(i,k,j)
2537 enddo
2538 exch_h1d(kts)=0.
2539 ! exch_h1d(kte+1)=0
2540 exch_m1d(kts)=0.
2541 ! exch_m1d(kte+1)=0
2542 rhoz1d(kts)=rho1d(kts)
2543 do k=kts+1,kte
2544 rhoz1d(k)=(rho1d(k)*dz1d(k-1)+rho1d(k-1)*dz1d(k))/ &
2545 & (dz1d(k-1)+dz1d(k))
2546 enddo
2547 rhoz1d(kte+1)=rho1d(kte)
2550 ! solve the diffusion for x-component of the wind
2552 call diff(kms,kme,kts,kte,dt,u1d,rho1d,rhoz1d,exch_m1d,a_u1d,b_u1d,sf1d, &
2553 & vl1d,dz1d,wu1d)
2555 ! solve the diffusion for y-component of the wind
2557 call diff(kms,kme,kts,kte,dt,v1d,rho1d,rhoz1d,exch_m1d,a_v1d,b_v1d,sf1d, &
2558 & vl1d,dz1d,wv1d)
2560 ! solve the diffusion for equivalent potential temperature
2562 call diff(kms,kme,kts,kte,dt,the1d,rho1d,rhoz1d,exch_h1d,a_t1d,b_t1d,sf1d, &
2563 & vl1d,dz1d,wt1d)
2565 ! solve the diffusion for the water vapor mixing ratio
2567 call diff(kms,kme,kts,kte,dt,qv1d,rho1d,rhoz1d,exch_h1d,a_q1d,b_q1d,sf1d, &
2568 & vl1d,dz1d,wq1d)
2570 ! solve the diffusion for the cloud water mixing ratio
2572 call diff(kms,kme,kts,kte,dt,qc1d,rho1d,rhoz1d,exch_h1d,a_qc1d,b_qc1d,sf1d, &
2573 & vl1d,dz1d,wqc1d)
2575 !
2576 ! compute the tendencies
2578 do k=kts,kte
2579 rublten(i,k,j)=(u1d(k)-u(i,k,j))/dt
2580 rvblten(i,k,j)=(v1d(k)-v(i,k,j))/dt
2581 thnew=the1d(k)/(QC(i,k,j)*(-ELOCP/T(i,k,j))+1)
2582 rthblten(i,k,j)=(thnew-th(i,k,j))/dt
2583 rqvblten(i,k,j)=(qv1d(k)-qv(i,k,j))/dt
2584 rqcblten(i,k,j)=(qc1d(k)-qc(i,k,j))/dt
2585 wu(i,k,j)=wu1d(k)
2586 wv(i,k,j)=wv1d(k)
2587 wt(i,k,j)=wt1d(k)
2588 wq(i,k,j)=wq1d(k)
2589 enddo
2590 enddo
2591 enddo
2592 !!!!!!!!!!!!
2595 END SUBROUTINE diff3d
2596 #if (EM_CORE==1)
2597 ! ===6=8===============================================================72
2598 SUBROUTINE diff4d(DT,DZ,SCALAR,IS_SCALAR,RHO &
2599 ,EXCH_H,EXCH_M &
2600 ,SCALAR_TEND &
2601 ,NUM_SCALAR, PARAM_FIRST_SCALAR &
2602 ,IDS,IDE,JDS,JDE,KDS,KDE &
2603 ,IMS,IME,JMS,JME,KMS,KME &
2604 ,ITS,ITE,JTS,JTE,KTS,KTE &
2605 )
2606 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2607 ! Based on subroutine written by Alberto Martilli, CIEMAT, Spain,
2608 ! e-mail:alberto_martilli@ciemat.es
2609 ! August 2008.
2610 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2611 ! ALberto
2612 ! This routine solves the vertical diffusion
2613 ! + other source/sink terms (surface fluxes, building drag, non local terms, etc.)
2614 ! for U,V, potential temperature and moisture. The equation should be written
2615 ! as follow, for a generic variable M:
2616 !
2617 ! dM 1 d dM
2618 ! ---- =---- ---(rho*K----)+AM+B
2619 ! dt rho dz dz
2620 ! where A and B are the implict and explcit coefficients of the source/sink terms
2621 ! (at the lowest level the surface fluxes should go in these terms)
2622 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2623 !-----------------------------------------------------------------------
2624 USE module_state_description, ONLY: &
2625 P_QNS, P_QNR, P_QNG, P_QT, P_QNH, P_QVOLG, P_QKE_ADV
2628 IMPLICIT NONE
2629 !
2630 !----------------------------------------------------------------------
2631 INTEGER,INTENT(IN) :: NUM_SCALAR, PARAM_FIRST_SCALAR
2632 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2633 & ,IMS,IME,JMS,JME,KMS,KME &
2634 & ,ITS,ITE,JTS,JTE,KTS,KTE
2636 ! inputs
2637 REAL, INTENT(IN) :: DT
2638 LOGICAL,INTENT(IN) :: IS_SCALAR
2639 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ ! Depth of vertical levels
2640 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME,NUM_SCALAR),INTENT(IN) :: SCALAR ! 4d scalar mixing ratio
2641 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: RHO ! Air Density
2642 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_H ! Turbulent coefficient for heat
2643 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: EXCH_M ! Turbulent coefficient for momentum
2645 ! outputs
2646 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME,NUM_SCALAR),INTENT(INOUT) :: SCALAR_TEND ! 4d scalar mixing ratio tendency
2647 ! locals (same meaning as above, but 1D)
2649 real s1d(kms:kme),exch_h1d(kms:kme)
2650 real exch_m1d(kms:kme)
2651 real dz1d(kms:kme),rho1d(kms:kme),rhoz1d(kms:kme)
2652 real sf1d(kms:kme),vl1d(kms:kme)
2653 real a_s1d(kms:kme),b_s1d(kms:kme)
2654 real ws1d(kms:kme)
2655 !
2656 integer i,k,j,im
2657 !----------------------------------------------------------------------
2658 s1d=0.
2659 exch_h1d=0.
2660 exch_m1d=0.
2661 rho1d=0.
2662 rhoz1d=0.
2663 sf1d=0.
2664 vl1d=0.
2665 a_s1d=0.
2666 b_s1d=0.
2668 DO im=PARAM_FIRST_SCALAR,NUM_SCALAR
2669 ! need to avoid mixing scalars associated with precipitating species (e.g. Nr)
2670 IF((IS_SCALAR .AND. im.NE.P_QNS .AND. im.NE.P_QNR .AND. im.NE.P_QNG &
2671 .AND. im.NE.P_QNH .AND. im.NE.P_QT .AND. im.NE.P_QVOLG .AND. im.NE.P_QKE_ADV) &
2672 .OR. (.not. IS_SCALAR)) THEN
2673 do j=jts,jte
2674 do i=its,ite
2676 ! put three D variables in temporary 1D variables
2678 do k=kts,kte
2679 s1d(k)=SCALAR(i,k,j,im)
2680 dz1d(k)=dz(i,k,j)
2681 rho1d(k)=rho(i,k,j)
2682 ! a_s1d(k)=a_s(i,k,j) ! implicit source
2683 ! b_s1d(k)=b_s(i,k,j) ! explicit source
2684 vl1d(k)=1.
2685 sf1d(k)=1.
2686 enddo
2687 sf1d(kte+1)=1.
2688 do k=kts,kte
2689 exch_h1d(k)=exch_h(i,k,j)
2690 exch_m1d(k)=exch_m(i,k,j)
2691 enddo
2692 exch_h1d(kts)=0.
2693 ! exch_h1d(kte+1)=0
2694 exch_m1d(kts)=0.
2695 ! exch_m1d(kte+1)=0
2696 rhoz1d(kts)=rho1d(kts)
2697 do k=kts+1,kte
2698 rhoz1d(k)=(rho1d(k)*dz1d(k-1)+rho1d(k-1)*dz1d(k))/ &
2699 & (dz1d(k-1)+dz1d(k))
2700 enddo
2701 rhoz1d(kte+1)=rho1d(kte)
2704 ! solve the diffusion for scalar
2706 call diff(kms,kme,kts,kte,dt,s1d,rho1d,rhoz1d,exch_h1d,a_s1d,b_s1d,sf1d, &
2707 & vl1d,dz1d,ws1d)
2709 !
2710 ! compute the tendencies
2712 do k=kts,kte
2713 scalar_tend(i,k,j,im)=(s1d(k)-scalar(i,k,j,im))/dt
2714 ! ws(i,k,j)=ws1d(k) ! output fluxes
2715 enddo
2716 enddo
2717 enddo
2718 ELSE
2719 ! print *,'scalar skipped im=',im, is_scalar
2720 ENDIF
2722 ENDDO ! im loop
2723 !!!!!!!!!!!!
2726 END SUBROUTINE diff4d
2727 ! ===6=8===============================================================72
2728 #endif
2730 subroutine diff(kms,kme,kts,kte,dt,co,da,daz,cd,aa,bb,sf,vl,dz,fc)
2732 !------------------------------------------------------------------------
2733 ! Calculation of the diffusion in 1D
2734 !------------------------------------------------------------------------
2735 ! - Input:
2736 ! nz : number of points
2737 ! iz1 : first calculated point
2738 ! co : concentration of the variable of interest
2739 ! dz : vertical levels
2740 ! cd : diffusion coefficients
2741 ! da : density of the air at the center
2742 ! daz : density of the air at the face
2743 ! dt : diffusion time step
2744 ! - Output:
2745 ! co :concentration of the variable of interest
2747 ! - Internal:
2748 ! cddz : constant terms in the equations
2749 !---------------------------------------------------------------------
2751 implicit none
2752 integer iz,iz1,izf
2753 integer kms,kme,kts,kte
2754 real dt,dzv
2755 real co(kms:kme),cd(kms:kme),dz(kms:kme)
2756 real da(kms:kme),daz(kms:kme)
2757 real cddz(kms:kme),fc(kms:kme),df(kms:kme)
2758 real a(kms:kme,3),c(kms:kme)
2759 real sf(kms:kme),vl(kms:kme)
2760 real aa(kms:kme),bb(kms:kme)
2763 ! Compute cddz=2*cd/dz
2765 cddz(kts)=sf(kts)*daz(kts)*cd(kts)/dz(kts)
2766 do iz=kts+1,kte
2767 cddz(iz)=2.*sf(iz)*daz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
2768 enddo
2769 cddz(kte+1)=sf(kte+1)*daz(kte+1)*cd(kte+1)/dz(kte)
2771 iz1=1
2772 izf=1
2774 do iz=iz1,kte-1
2776 dzv=vl(iz)*dz(iz)
2777 a(iz,1)=-cddz(iz)*dt/dzv/da(iz)
2778 a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/da(iz)-aa(iz)*dt
2779 a(iz,3)=-cddz(iz+1)*dt/dzv/da(iz)
2780 c(iz)=co(iz)+bb(iz)*dt
2781 enddo
2783 do iz=kte-(izf-1),kte
2784 a(iz,1)=0.
2785 a(iz,2)=1
2786 a(iz,3)=0.
2787 c(iz)=co(iz)
2788 enddo
2789 call invert (kms,kme,kts,kte,a,c,co)
2791 !let compute the fluxes, as diagnostic variable
2793 do iz=kts,iz1
2794 fc(iz)=0.
2795 enddo
2797 do iz=iz1+1,kte
2798 fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/da(iz)
2799 enddo
2803 return
2804 end subroutine diff
2805 ! ===6=8===============================================================72
2807 subroutine invert(kms,kme,kts,kte,a,c,x)
2809 !ccccccccccccccccccccccccccccccc
2810 ! Aim: Inversion and resolution of a tridiagonal matrix
2811 ! A X = C
2812 ! Input:
2813 ! a(*,1) lower diagonal (Ai,i-1)
2814 ! a(*,2) principal diagonal (Ai,i)
2815 ! a(*,3) upper diagonal (Ai,i+1)
2816 ! c
2817 ! Output
2818 ! x results
2819 !ccccccccccccccccccccccccccccccc
2821 implicit none
2822 integer kms,kme,kts,kte,in
2823 real a(kms:kme,3),c(kms:kme),x(kms:kme)
2825 do in=kte-1,kts,-1
2826 c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
2827 a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
2828 enddo
2830 do in=kts+1,kte
2831 c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
2832 enddo
2834 do in=kts,kte
2835 x(in)=c(in)/a(in,2)
2836 enddo
2838 return
2839 end subroutine invert
2841 ! ===6=8===============================================================72
2843 END MODULE module_pbl_driver