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