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