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