| blob | 4bcabf09c9864bde7ac1f2906b080a3db028ab22 |
1 #define MODAL_AERO
2 MODULE module_bl_camuwpbl_driver
3 !Note: Comments starting with "!!" are directly taken from CAM's interface routine for the UW PBL (vertical_diffusion.F90)
4 !This modules is based on vertical_diffusion.F90 of CAM5
6 !List of possible modifications in the future:
7 !_____________________________________________
8 !1. Constituents(per Kg of DRY air or "dry" constituents) are not diffused currently. The call which
9 ! makes DRY constituents diffusion decision flag 'true' is commented out in camuwpblinit subroutine.
10 ! Therefore the flag is ALWAYS false for DRY constituent diffusion.
11 !
12 ! REASON: DRY constituents are Aerosols and other species (excepts for cloud mass mixing ratios and number concentrations).
13 ! ------ In WRF, DRY constituents appear only when WRF_CHEM simulations are run. DRY constituents are vertically mixed(diffused)
14 ! in dry_dep_driver of WRF_CHEM. Therefore, DRY constituents are not diffused in this PBL subroutine
15 !
16 !2. The liquid number concentrations is NOT diffused(or mixed) in PBL to mimic CAM5, which diffuses liquid number
17 ! concentrations in dropmixnuc.
18 !3. Surface fluxes for ALL the constituents are ZERO except for the water vapours (1st constituent)
19 !4. Molecular diffusion is turned off for now
20 !5. Mountain stresses are not computed for now (all calls to Mountain stresses currently may hold undefined variables)
21 !6. *Formulas used for computing surface stresses is based on the formula given in CAM5's BareGround.F90.
22 ! This formula may not work well for ocean surfaces (possible future modification)
23 !_____________________________________________
26 !!----------------------------------------------------------------------------------------------------- !
27 !! Module to compute vertical diffusion of momentum, moisture, trace constituents !
28 !! and static energy. Separate modules compute !
29 !! 1. stresses associated with turbulent flow over orography !
30 !! ( turbulent mountain stress ) !
31 !! 2. eddy diffusivities, including nonlocal tranport terms !
32 !! 3. molecular diffusivities !
33 !! 4. coming soon... gravity wave drag !
34 !! Lastly, a implicit diffusion solver is called, and tendencies retrieved by !
35 !! differencing the diffused and initial states. !
36 !! !
37 !! Calling sequence: !
38 !! !
39 !! vertical_diffusion_init Initializes vertical diffustion constants and modules !
40 !! init_molec_diff Initializes molecular diffusivity module !
41 !! init_eddy_diff Initializes eddy diffusivity module (includes PBL) !
42 !! init_tms Initializes turbulent mountain stress module !
43 !! init_vdiff Initializes diffusion solver module !
44 !! vertical_diffusion_ts_init Time step initialization (only used for upper boundary condition) !
45 !! vertical_diffusion_tend Computes vertical diffusion tendencies !
46 !! compute_tms Computes turbulent mountain stresses !
47 !! compute_eddy_diff Computes eddy diffusivities and countergradient terms !
48 !! compute_vdiff Solves vertical diffusion equations, including molecular diffusivities !
49 !! !
50 !!---------------------------Code history-------------------------------------------------------------- !
51 !! J. Rosinski : Jun. 1992 !
52 !! J. McCaa : Sep. 2004 !
53 !! S. Park : Aug. 2006, Dec. 2008. Jan. 2010 !
54 ! B. Singh : Nov. 2010 (ported to WRF by balwinder.singh@pnl.gov)
55 !!----------------------------------------------------------------------------------------------------- !
57 use shr_kind_mod, only : r8 => shr_kind_r8
58 use module_cam_support, only : pcols, pver, pverp, endrun, iulog,fieldname_len,pcnst =>pcnst_runtime
59 use constituents, only : qmin
60 use diffusion_solver, only : vdiff_selector
61 use physconst, only : &
62 cpair , & ! Specific heat of dry air
63 gravit , & ! Acceleration due to gravity
64 rair , & ! Gas constant for dry air
65 zvir , & ! rh2o/rair - 1
66 latvap , & ! Latent heat of vaporization
67 latice , & ! Latent heat of fusion
68 karman , & ! von Karman constant
69 mwdry , & ! Molecular weight of dry air
70 avogad , & ! Avogadro's number
71 boltz ! Boltzman's constant
73 implicit none
74 private
75 save
77 !! ----------------- !
78 !! Public interfaces !
79 !! ----------------- !
81 public camuwpblinit ! Initialization
82 public camuwpbl ! Driver for the PBL scheme
83 public vd_register ! Init routine for constituents
85 !! ------------ !
86 !! Private data !
87 !! ------------ !
89 character(len=16) :: eddy_scheme !! Default set in phys_control.F90, use namelist to change
90 !! 'HB' = Holtslag and Boville (default)
91 !! 'HBR' = Holtslag and Boville and Rash
92 !! 'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme )
93 integer, parameter :: nturb = 5 !! Number of iterations for solution ( when 'diag_TKE' scheme is selected )
94 logical, parameter :: wstarent = .true. !! Use wstar (.true.) or TKE (.false.) entrainment closure ( when 'diag_TKE' scheme is selected )
95 logical :: do_pseudocon_diff = .false. !! If .true., do pseudo-conservative variables diffusion
98 character(len=16) :: shallow_scheme !! For checking compatibility between eddy diffusion and shallow convection schemes
99 !! 'Hack' = Hack Shallow Convection Scheme
100 !! 'UW' = Park and Bretherton ( UW Shallow Convection Scheme )
101 character(len=16) :: microp_scheme !! Microphysics scheme
103 logical :: do_molec_diff = .false. !! Switch for molecular diffusion
104 logical :: do_tms !! Switch for turbulent mountain stress
105 real(r8) :: tms_orocnst !! Converts from standard deviation to height
106 real(r8) :: tms_z0fac ! Converts from standard deviation to height
108 type(vdiff_selector) :: fieldlist_wet !! Logical switches for moist mixing ratio diffusion
109 type(vdiff_selector) :: fieldlist_dry !! Logical switches for dry mixing ratio diffusion
110 integer :: ntop !! Top interface level to which vertical diffusion is applied ( = 1 ).
111 integer :: nbot !! Bottom interface level to which vertical diffusion is applied ( = pver ).
112 integer :: tke_idx, kvh_idx, kvm_idx !! TKE and eddy diffusivity indices for fields in the physics buffer
113 integer :: turbtype_idx, smaw_idx !! Turbulence type and instability functions
114 integer :: tauresx_idx, tauresy_idx !! Redisual stress for implicit surface stress
116 integer :: ixcldice, ixcldliq !! Constituent indices for cloud liquid and ice water
117 integer :: ixnumice, ixnumliq
118 #ifdef MODAL_AERO
119 integer :: ixndrop
120 #endif
121 integer :: wgustd_index
122 CONTAINS
124 subroutine camuwpbl(dt,u_phy,v_phy,th_phy,rho,qv_curr,hfx,qfx,ustar,p8w &
125 ,p_phy,z,t_phy,qc_curr,qi_curr,z_at_w,cldfra_old_mp,cldfra,ht &
126 ,rthratenlw,exner,is_CAMMGMP_used &
127 ,itimestep,qnc_curr,qni_curr,wsedl3d &
128 ,ids,ide, jds,jde, kds,kde &
129 ,ims,ime, jms,jme, kms,kme &
130 ,its,ite, jts,jte, kts,kte &
131 !Intent-inout
132 ,tauresx2d,tauresy2d &
133 ,rublten,rvblten,rthblten,rqiblten,rqniblten,rqvblten,rqcblten &
134 ,kvm3d,kvh3d &
135 !Intent-out
136 ,tpert2d,qpert2d,wpert2d,smaw3d,turbtype3d &
137 ,tke_pbl,pblh2d,kpbl2d )
139 !!---------------------------------------------------- !
140 !! This is an interface routine for vertical diffusion !
141 ! This subroutine is called by module_pbl_driver and !
142 ! it calls: wrf_error_fatal,compute_tms, !
143 ! compute_eddy_diff,aqsat,compute_vdiff and !
144 ! positive_moisture.
145 !!---------------------------------------------------- !
146 use module_cam_support, only : pcols
147 use trb_mtn_stress, only : compute_tms
148 use eddy_diff, only : compute_eddy_diff
149 use wv_saturation, only : fqsatd, aqsat
150 use molec_diff, only : compute_molec_diff, vd_lu_qdecomp
151 use constituents, only : qmincg, qmin
152 use diffusion_solver !!, only : compute_vdiff, any, operator(.not.)
153 #ifdef MODAL_AERO
154 use modal_aero_data
155 #endif
157 implicit none
159 !------------------------------------------------------------------------!
160 ! Input !
161 !------------------------------------------------------------------------!
162 logical, intent(in) :: is_CAMMGMP_used
163 integer, intent(in) :: ids,ide, jds,jde, kds,kde
164 integer, intent(in) :: ims,ime, jms,jme, kms,kme
165 integer, intent(in) :: its,ite, jts,jte, kts,kte
166 integer, intent(in) :: itimestep
168 real, intent(in) :: dt ! Time step (s)
170 real, dimension( ims:ime,jms:jme ), intent(in) :: hfx !Sensible heat flux at surface (w/m2)
171 real, dimension( ims:ime,jms:jme ), intent(in) :: qfx !Moisture flux at surface (kg m-2 s-1)
172 real, dimension( ims:ime,jms:jme ), intent(in) :: ustar !Friction velocity (m/s)
173 real, dimension( ims:ime,jms:jme ), intent(in) :: ht !Terrain height (m)
175 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: v_phy !Y-component of wind (m/s)
176 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: u_phy !X-component of wind (m/s)
177 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: th_phy !Potential temperature (K)
178 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rho !Air density (kg/m3)
179 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qv_curr !Water vapor mixing ratio - Moisture (kg/kg)
180 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qc_curr !Cloud water mixing ratio - Cloud liq (kg/kg)
181 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qi_curr !Ice mixing ratio - Cloud ice (kg/kg)
182 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qnc_curr !Liq # mixing ratio - Cloud liq # (#/kg)
183 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: qni_curr !Ice # mixing ratio - Cloud ice # (#/kg)
184 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p_phy !Pressure at mid-level (Pa)
185 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p8w !Pressure at level interface (Pa)
186 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z !Height above sea level at mid-level (m)
187 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: t_phy !Temperature (K)
188 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_at_w !Height above sea level at layer interfaces (m)
189 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra_old_mp !Cloud fraction [unitless]
190 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cldfra !Cloud fraction [unitless]
191 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: exner !Dimensionless pressure [unitless]
192 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: rthratenlw !Tendency for LW ( K/s)
193 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: wsedl3d !Sedimentation velocity of stratiform liquid cloud droplet (m/s)
197 !------------------------------------------------------------------------!
198 ! Output !
199 !------------------------------------------------------------------------!
200 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rublten !Tendency for u_phy (Pa m s-2)
201 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rvblten !Tendency for v_phy (Pa m s-2)
202 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rthblten !Tendency for th_phy (Pa K s-1)
203 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqvblten !Tendency for qv_curr (Pa kg kg-1 s-1)
204 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqcblten !Tendency for qc_curr (Pa kg kg-1 s-1)
205 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqiblten !Tendency for qi_curr (Pa kg kg-1 s-1)
206 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rqniblten !Tendency for qni_curr(Pa # kg-1 s-1)
208 integer, dimension( ims:ime,jms:jme ), intent(out) :: kpbl2d !Layer index containing PBL top within or at the base interface
209 real, dimension( ims:ime,jms:jme ), intent(out) :: pblh2d !Planetary boundary layer height (m)
211 real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: tke_pbl !Turbulence kinetic energy at midpoints (m^2/s^2)
212 real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: turbtype3d !Turbulent interface types [ no unit ]
213 real, dimension( ims:ime, kms:kme, jms:jme ), intent(out) :: smaw3d !Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
217 !---------------------------------------------------------------------------!
218 ! Local, carried on from one timestep to the other (defined in registry)!
219 !---------------------------------------------------------------------------!
221 real, dimension( ims:ime, jms:jme ) , intent(inout ):: tauresx2d,tauresy2d !X AND Y-COMP OF RESIDUAL STRESSES(m^2/s^2)
222 real, dimension( ims:ime, jms:jme ) , intent(out) :: tpert2d ! Convective temperature excess (K)
223 real, dimension( ims:ime, jms:jme ) , intent(out) :: qpert2d ! Convective humidity excess (kg/kg)
224 real, dimension( ims:ime, jms:jme ) , intent(out) :: wpert2d ! Turbulent velocity excess (m/s)
226 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: kvm3d,kvh3d !Eddy diffusivity for momentum and heat(m^2/s)
229 !---------------------------------------------------------------------------!
230 ! Local !
231 !---------------------------------------------------------------------------!
233 character(128) :: errstring ! Error status for compute_vdiff
234 logical :: kvinit ! Tell compute_eddy_diff/ caleddy to initialize kvh, kvm (uses kvf)
235 logical :: is_first_step ! Flag to know if this a first time step
236 integer :: i,j,k,itsp1,itile_len,ktep1,kflip,ncol,ips,icnst,m,kp1
237 integer :: lchnk ! Chunk identifier
238 real(r8) :: tauFac, uMean, dp, multFrc
239 real(r8) :: ztodt ! 2*delta-t (s)
240 real(r8) :: rztodt ! 1./ztodt (1/s)
241 real(r8),parameter :: invalidVal = -999888777.0_r8
243 real(r8) :: topflx( pcols) ! Molecular heat flux at top interface
244 real(r8) :: wpert( pcols) ! Turbulent velocity excess (m/s)
245 real(r8) :: tauresx( pcols) ! [Residual stress to be added in vdiff to correct...
246 real(r8) :: tauresy( pcols) ! for turb stress mismatch between sfc and atm accumulated.] (N/m2)
247 real(r8) :: ipbl( pcols) ! If 1, PBL is CL, while if 0, PBL is STL.
248 real(r8) :: kpblh( pcols) ! Layer index containing PBL top within or at the base interface
249 real(r8) :: wstarPBL(pcols) ! Convective velocity within PBL (m/s)
250 real(r8) :: sgh( pcols) ! Standard deviation of orography (m)
251 real(r8) :: landfrac(pcols) ! Land fraction [unitless]
252 real(r8) :: taux( pcols) ! x surface stress (N/m2)
253 real(r8) :: tauy( pcols) ! y surface stress (N/m2)
254 real(r8) :: tautotx( pcols) ! U component of total surface stress (N/m2)
255 real(r8) :: tautoty( pcols) ! V component of total surface stress (N/m2)
256 real(r8) :: ksrftms( pcols) ! Turbulent mountain stress surface drag coefficient (kg/s/m2)
257 real(r8) :: tautmsx( pcols) ! U component of turbulent mountain stress (N/m2)
258 real(r8) :: tautmsy( pcols) ! V component of turbulent mountain stress (N/m2)
259 real(r8) :: ustar8( pcols) ! Surface friction velocity (m/s)
260 real(r8) :: pblh( pcols) ! Planetary boundary layer height (m)
261 real(r8) :: tpert( pcols) ! Convective temperature excess (K)
262 real(r8) :: qpert( pcols) ! Convective humidity excess (kg/kg)
263 real(r8) :: shflx( pcols) ! Surface sensible heat flux (w/m2)
264 real(r8) :: phis( pcols) ! Geopotential at terrain height (m2/s2)
266 real(r8) :: cldn8( pcols,kte) ! New stratus fraction (fraction)
267 real(r8) :: qrl8( pcols,kte) ! LW radiative cooling rate(W/kg*Pa)
268 real(r8) :: wsedl8( pcols,kte) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
269 real(r8) :: dtk( pcols,kte) ! T tendency from KE dissipation
270 real(r8) :: qt( pcols,kte) !
271 real(r8) :: sl_prePBL( pcols,kte) !
272 real(r8) :: qt_prePBL( pcols,kte) !
273 real(r8) :: slten( pcols,kte) !
274 real(r8) :: qtten( pcols,kte) !
275 real(r8) :: sl( pcols,kte) !
276 real(r8) :: ftem( pcols,kte) ! Saturation vapor pressure before PBL
277 real(r8) :: ftem_prePBL(pcols,kte) ! Saturation vapor pressure before PBL
278 real(r8) :: ftem_aftPBL(pcols,kte) ! Saturation vapor pressure after PBL
279 real(r8) :: tem2( pcols,kte) ! Saturation specific humidity and RH
280 real(r8) :: t_aftPBL( pcols,kte) ! Temperature after PBL diffusion
281 real(r8) :: tten( pcols,kte) ! Temperature tendency by PBL diffusion
282 real(r8) :: rhten( pcols,kte) ! RH tendency by PBL diffusion
283 real(r8) :: qv_aft_PBL( pcols,kte) ! qv after PBL diffusion
284 real(r8) :: ql_aft_PBL( pcols,kte) ! ql after PBL diffusion
285 real(r8) :: qi_aft_PBL( pcols,kte) ! qi after PBL diffusion
286 real(r8) :: s_aft_PBL( pcols,kte) ! s after PBL diffusion
287 real(r8) :: u_aft_PBL( pcols,kte) ! u after PBL diffusion
288 real(r8) :: v_aft_PBL( pcols,kte) ! v after PBL diffusion
289 real(r8) :: qv_pro( pcols,kte) !
290 real(r8) :: ql_pro( pcols,kte) !
291 real(r8) :: qi_pro( pcols,kte) !
292 real(r8) :: s_pro( pcols,kte) !
293 real(r8) :: t_pro( pcols,kte) !
294 real(r8) :: u8( pcols,kte) ! x component of velocity in CAM's data structure (m/s)
295 real(r8) :: v8( pcols,kte) ! y component of velocity in CAM's data structure (m/s)
296 real(r8) :: t8( pcols,kte) !
297 real(r8) :: pmid8( pcols,kte) ! Pressure at the midpoints in CAM's data structure (Pa)
298 real(r8) :: pmiddry8( pcols,kte) ! Dry Pressure at the midpoints in CAM's data structure (Pa)
299 real(r8) :: zm8( pcols,kte) ! Height at the midpoints in CAM's data structure (m)
300 real(r8) :: exner8( pcols,kte) ! exner function in CAM's data structure
301 real(r8) :: s8( pcols,kte) ! Dry static energy (m2/s2)
302 real(r8) :: rpdel8( pcols,kte) ! Inverse of pressure difference (1/Pa)
303 real(r8) :: pdel8( pcols,kte) ! Pressure difference (Pa)
304 real(r8) :: rpdeldry8( pcols,kte) ! Inverse of dry pressure difference (1/Pa)
305 real(r8) :: pdeldry8( pcols,kte) ! dry pressure difference (1/Pa)
306 REAL(r8) :: stnd( pcols,kte) ! Heating rate (dry static energy tendency, W/kg)
308 real(r8) :: tke8( pcols,kte+1) ! Turbulent kinetic energy [ m2/s2 ]
309 real(r8) :: turbtype( pcols,kte+1) ! Turbulent interface types [ no unit ]
310 real(r8) :: smaw( pcols,kte+1) ! Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
311 ! This is 1 when neutral condition (Ri=0), 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19.
312 real(r8) :: cgs( pcols,kte+1) ! Counter-gradient star [ cg/flux ]
313 real(r8) :: cgh( pcols,kte+1) ! Counter-gradient term for heat [ J/kg/m ]
314 real(r8) :: kvh( pcols,kte+1) ! Eddy diffusivity for heat [ m2/s ]
315 real(r8) :: kvm( pcols,kte+1) ! Eddy diffusivity for momentum [ m2/s ]
316 real(r8) :: kvq( pcols,kte+1) ! Eddy diffusivity for constituents [ m2/s ]
317 real(r8) :: kvh_in( pcols,kte+1) ! kvh from previous timestep [ m2/s ]
318 real(r8) :: kvm_in( pcols,kte+1) ! kvm from previous timestep [ m2/s ]
319 real(r8) :: bprod( pcols,kte+1) ! Buoyancy production of tke [ m2/s3 ]
320 real(r8) :: sprod( pcols,kte+1) ! Shear production of tke [ m2/s3 ]
321 real(r8) :: sfi( pcols,kte+1) ! Saturation fraction at interfaces [ fraction ]
322 real(r8) :: pint8( pcols,kte+1) ! Pressure at the interfaces in CAM's data structure (Pa)
323 real(r8) :: pintdry8( pcols,kte+1) ! Dry pressure at the interfaces in CAM's data structure (Pa)
324 real(r8) :: zi8( pcols,kte+1) ! Height at the interfacesin CAM's data structure (m)
326 real(r8) :: cloud( pcols,kte,pcnst) ! Holder for cloud water and ice (q in cam)
327 real(r8) :: cloudtnd( pcols,kte,pcnst) ! Holder for cloud tendencies (ptend_loc%q in cam)
328 real(r8) :: wind_tends(pcols,kte,2) ! Wind component tendencies (m/s2)
330 real(r8) :: cflx(pcols,pcnst) ! Surface constituent flux [ kg/m2/s ]
332 #ifdef MODAL_AERO
333 real(r8) :: tmp1(pcols) ! Temporary storage
334 integer :: l, lspec
335 #endif
336 !! ----------------------- !
337 !! Main Computation Begins !
338 !! ----------------------- !
340 do icnst = 1 , pcnst
341 !Setting curface fluxes for all constituents to be zero. Later, cflx(:,1) is populated by the water vapour surface flux
342 cflx(:pcols,icnst) = 0.0_r8
343 enddo
345 is_first_step = .false.
346 if(itimestep == 1) then
347 is_first_step = .true.
348 endif
349 ncol = pcols
350 ztodt = DT
351 rztodt = 1.0_r8 / ztodt
353 !Few definitions in this subroutine require that ncol==1
354 if(ncol .NE. 1) then
355 call wrf_error_fatal('Number of CAM Columns (NCOL) in CAMUWPBL scheme must be 1')
356 endif
358 !Initialize errstring
359 errstring = ''
362 !-------------------------------------------------------------------------------------!
363 !Map CAM variables to the corresponding WRF variables !
364 !Loop over the points in the tile and treat them each as a CAM Chunk !
365 !-------------------------------------------------------------------------------------!
366 itsp1 = its - 1
367 itile_len = ite - itsp1
368 do j = jts , jte
369 do i = its , ite
371 lchnk = (j - jts) * itile_len + (i - itsp1) !1-D index location from a 2-D tile
372 phis(1) = ht(i,j) * gravit !Used for computing dry static energy
374 !Flip vertically quantities computed at the mid points
375 ktep1 = kte + 1
376 do k = kts,kte
377 kflip = ktep1 - k
379 !Loop is used as vector assignment may take more computational time
380 do icnst = 1 , pcnst
381 !Setting cloud and cloudtnd to values (obtained from module_cam_support) which should produce errors if used
382 !1st,2nd,3rd and 5th constituents are used (see the assignments below) and diffused presently in this scheme
383 cloud(1,kflip,icnst) = invalidVal
384 cloudtnd(1,kflip,icnst) = invalidVal
385 enddo
387 u8( 1,kflip) = u_phy(i,k,j) ! X-component of velocity at the mid-points (m/s) [state%u in CAM]
388 v8( 1,kflip) = v_phy(i,k,j) ! Y-component of velocity at the mid-points (m/s) [state%v in CAM]
390 pmid8( 1,kflip) = p_phy(i,k,j) ! Pressure at the mid-points (Pa) [state%pmid in CAM]
392 dp = p8w(i,k,j) - p8w(i,k+1,j) !Change in pressure (Pa)
393 pdel8( 1,kflip) = dp
394 rpdel8( 1,kflip) = 1.0_r8/dp ! Reciprocal of pressure difference (1/Pa) [state%rpdel in CAM]
395 zm8( 1,kflip) = z(i,k,j)-ht(i,j) ! Height above the ground at the midpoints (m) [state%zm in CAM]
396 t8( 1,kflip) = t_phy(i,k,j) ! Temprature at the mid points (K) [state%t in CAM]
398 s8( 1,kflip) = cpair *t8(1,kflip) + gravit*zm8(1,kflip) + phis(1) ! Dry static energy (m2/s2) [state%s in CAM]-Formula tested in vertical_diffusion.F90 of CAM
399 qrl8( 1,kflip) = rthratenlw(i,k,j) * exner(i,k,j) * cpair * dp ! Long Wave heating rate (W/kg*Pa)- Formula obtained from definition of qrlw(pcols,pver) in eddy_diff.F90 in CAM
401 wsedl8( 1,kflip) = wsedl3d(i,k,j) ! Sedimentation velocity of stratiform liquid cloud droplet (m/s)
403 !Following three formulas are obtained from ported CAM's ZM cumulus scheme
404 !Values of 0 cause a crash in entropy
405 multFrc = 1._r8/(1._r8 + qv_curr(i,k,j))
406 cloud( 1,kflip,1) = max( qv_curr(i,k,j) * multFrc, 1.0e-30_r8 ) !Specific humidity [state%q(:,:,1) in CAM]
407 cloud( 1,kflip,ixcldliq) = qc_curr(i,k,j) * multFrc !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
408 cloud( 1,kflip,ixcldice) = qi_curr(i,k,j) * multFrc !cloud ice [state%q(:,:,3) in CAM]
409 cloud( 1,kflip,ixnumliq) = qnc_curr(i,k,j) * multFrc
410 cloud( 1,kflip,ixnumice) = qni_curr(i,k,j) * multFrc
412 exner8(1,kflip) = exner(i,k,j) !Exner function (no units)
413 if(is_CAMMGMP_used) then
414 cldn8( 1,kflip) = cldfra_old_mp(i,k,j) !Cloud fraction (no unit)
415 else
416 cldn8( 1,kflip) = cldfra(i,k,j) !Cloud fraction (no unit)
417 endif
419 !Following formula is obtained from physics_types.F90 of CAM (CESM101)
420 pdeldry8(1,kflip) = pdel8(1,kflip)*(1._r8-cloud(1,kflip,1))
421 rpdeldry8(1,kflip) = 1._r8/pdeldry8(1,kflip)
423 enddo
425 do k = kts,kte+1
426 kflip = kte - k + 2
428 pint8( 1,kflip) = p8w( i,k,j) ! Pressure at interfaces [state%pint in CAM]
430 zi8( 1,kflip) = z_at_w(i,k,j) -ht(i,j) ! Height at interfaces [state%zi in CAM]
432 !Initialize Variables to zero, these are outputs from the "compute_eddy_diff"
433 kvq(1,kflip) = 0.0_r8 ! Eddy diffusivity for constituents (m2/s)
434 cgh(1,kflip) = 0.0_r8 ! Counter-gradient term for heat
435 cgs(1,kflip) = 0.0_r8 ! Counter-gradient star (cg/flux)
436 if( is_first_step ) then
437 kvm3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for heat (m2/s)
438 kvh3d(i,k,j) = 0.0_r8 ! Eddy diffusivity for momentum (m2/s)
439 endif
440 kvh(1,kflip) = kvh3d(i,k,j)
441 kvm(1,kflip) = kvm3d(i,k,j)
442 end do
444 !Compute pintdry8 and pmiddry8
445 !Following formula is obtained from physics_types.F90 of CAM (CESM101)
446 pintdry8(1,1) = pint8(1,1)
447 do k = 1, pver
448 pintdry8(1,k+1) = pintdry8(1,k) + pdeldry8(1,k)
449 pmiddry8(1,k) = (pintdry8(1,k+1)+ pintdry8(1,k))*0.5_r8
450 end do
452 shflx( ncol) = hfx(i,j) ! Surface sensible heat flux (w/m2)
454 !SGH and LANDFRAC are inputs for the compute_tms subroutine. Presently set to zero as do_tms is always false for now
455 sgh( ncol) = 0.0_r8 ! Standard deviation of orography (m)
456 landfrac(ncol) = 0.0_r8 ! Fraction (unitless)
458 uMean = sqrt( u_phy(i,kts,j) * u_phy(i,kts,j) + v_phy(i,kts,j) * v_phy(i,kts,j) ) ! Mean velocity
459 tauFac = rho(i,kts,j) * ustar(i,j) *ustar(i,j)/uMean
461 taux(ncol) = -tauFac * u_phy(i,kts,j) ! x surface stress (N/m2) [Formulation obtained from CAM's BareGround.F90]
462 tauy(ncol) = -tauFac * v_phy(i,kts,j) ! y surface stress (N/m2)
464 !! Retrieve 'tauresx, tauresy' from from the last timestep
465 if( is_first_step ) then
466 tauresx2d(i,j) = 0._r8
467 tauresy2d(i,j) = 0._r8
468 endif
469 tauresx(:ncol) = tauresx2d(i,j)
470 tauresy(:ncol) = tauresy2d(i,j)
472 !! All variables are modified by vertical diffusion
474 !!------------------------------------------!
475 !! Computation of turbulent mountain stress !
476 !!------------------------------------------!
478 !! Consistent with the computation of 'normal' drag coefficient, we are using
479 !! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v'
480 !! within the iteration loop of the PBL scheme.
482 if( do_tms ) then
483 call compute_tms( pcols , pver , ncol , &
484 u8 , v8 , t8 , pmid8 , &
485 exner8 , zm8 , sgh , ksrftms , &
486 tautmsx , tautmsy , landfrac )
487 !! Here, both 'taux, tautmsx' are explicit surface stresses.
488 !! Note that this 'tautotx, tautoty' are different from the total stress
489 !! that has been actually added into the atmosphere. This is because both
490 !! taux and tautmsx are fully implicitly treated within compute_vdiff.
491 !! However, 'tautotx, tautoty' are not used in the actual numerical
492 !! computation in this module.
493 tautotx(:ncol) = taux(:ncol) + tautmsx(:ncol)
494 tautoty(:ncol) = tauy(:ncol) + tautmsy(:ncol)
495 else
496 ksrftms(:ncol) = 0.0_r8
497 tautotx(:ncol) = taux(:ncol)
498 tautoty(:ncol) = tauy(:ncol)
499 endif
501 !-------------------------------------------------------------------------------------!
502 !We are currenly using just water vapour flux at the surface, rest are set to zero
503 !-------------------------------------------------------------------------------------!
504 cflx(:pcols,1) = qfx(i,j) ! Surface constituent flux (kg/m2/s)
506 !Following variables are initialized to zero, they are the output from the "compute_eddy_diff" call
507 ustar8( :pcols) = 0.0_r8 ! Surface friction velocity (m/s)
508 pblh( :pcols) = 0.0_r8 ! Planetary boundary layer height (m )
509 ipbl( :pcols) = 0.0_r8 ! If 1, PBL is CL, while if 0, PBL is STL.
510 kpblh( :pcols) = 0.0_r8 ! Layer index containing PBL top within or at the base interface
511 wstarPBL(:pcols) = 0.0_r8 ! Convective velocity within PBL (m/s)
514 !!----------------------------------------------------------------------- !
515 !! Computation of eddy diffusivities - Select appropriate PBL scheme !
516 !!----------------------------------------------------------------------- !
518 select case (eddy_scheme)
519 case ( 'diag_TKE' )
521 !! ---------------------------------------------------------------- !
522 !! At first time step, have eddy_diff.F90:caleddy() use kvh=kvm=kvf !
523 !! This has to be done in compute_eddy_diff after kvf is calculated !
524 !! ---------------------------------------------------------------- !
526 kvinit = .false.
527 if(is_first_step) then
528 kvinit = .true.
529 endif
530 !! ---------------------------------------------- !
531 !! Get LW radiative heating out of physics buffer !
532 !! ---------------------------------------------- !
534 !! Retrieve eddy diffusivities for heat and momentum from physics buffer
535 !! from last timestep ( if first timestep, has been initialized by inidat.F90 )
537 kvm_in(:ncol,:) = kvm(:ncol,:)
538 kvh_in(:ncol,:) = kvh(:ncol,:)
539 call compute_eddy_diff( lchnk , pcols , pver , ncol , t8 , cloud(:,:,1) , ztodt , &
540 cloud(:,:,2), cloud(:,:,3), s8 , rpdel8 , cldn8 , qrl8 , wsedl8 , zm8 , zi8 , &
541 pmid8 , pint8 , u8 , v8 , taux , tauy , shflx , cflx(:,1), wstarent, &
542 nturb , ustar8 , pblh , kvm_in , kvh_in , kvm , kvh , kvq , cgh , &
543 cgs , tpert , qpert , wpert , tke8 , bprod , sprod , sfi , fqsatd , &
544 kvinit , tauresx , tauresy, ksrftms, ipbl(:), kpblh(:), wstarPBL(:) , turbtype , smaw )
546 !! ----------------------------------------------- !
547 !! Store TKE in pbuf for use by shallow convection !
548 !! ----------------------------------------------- !
550 !! Store updated kvh, kvm in pbuf to use here on the next timestep
551 do k = kts,kte+1
552 kflip = kte - k + 2
554 kvh3d(i,k,j) = kvh(1,kflip)
555 kvm3d(i,k,j) = kvm(1,kflip)
556 end do
559 end select
561 !!------------------------------------ !
562 !! Application of diffusivities !
563 !!------------------------------------ !
564 cloudtnd( :ncol,:,:) = cloud(:ncol,:,:)
565 stnd( :ncol,: ) = s8( :ncol,: )
566 wind_tends(:ncol,:,1) = u8( :ncol,: )
567 wind_tends(:ncol,:,2) = v8( :ncol,: )
569 !!------------------------------------------------------ !
570 !! Write profile output before applying diffusion scheme !
571 !!------------------------------------------------------ !
573 sl_prePBL(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
574 - ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
575 qt_prePBL(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
576 + cloudtnd(:ncol,:pver,ixcldice)
578 call aqsat( t8, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
579 ftem_prePBL(:ncol,:) = cloud(:ncol,:,1)/ftem(:ncol,:)*100._r8
581 !! --------------------------------------------------------------------------------- !
582 !! Call the diffusivity solver and solve diffusion equation !
583 !! The final two arguments are optional function references to !
584 !! constituent-independent and constituent-dependent moleculuar diffusivity routines !
585 !! --------------------------------------------------------------------------------- !
587 !! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and
588 !! separately print out as diagnostic output, because these are different from
589 !! the explicit 'tautotx, tautoty' computed above.
590 !! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones.
592 if( any(fieldlist_wet) ) then
593 call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmid8, pint8, rpdel8, t8, ztodt, &
594 taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, qmincg, &
595 fieldlist_wet, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
596 tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
597 compute_molec_diff, vd_lu_qdecomp)
598 end if
600 if( errstring .ne. '' ) then
601 call wrf_error_fatal(errstring)
602 endif
604 if( any( fieldlist_dry ) ) then
605 if( do_molec_diff ) then
606 errstring = "Design flaw: dry vdiff not currently supported with molecular diffusion"
607 call wrf_error_fatal(errstring)
608 end if
610 call compute_vdiff( lchnk, pcols, pver, pcnst, ncol, pmiddry8, pintdry8, rpdeldry8, t8, &
611 ztodt, taux, tauy, shflx, cflx, ntop, nbot, kvh, kvm, kvq, cgs, cgh, zi8, ksrftms, &
612 qmincg, fieldlist_dry, wind_tends(:,:,1), wind_tends(:,:,2), cloudtnd, stnd, tautmsx, &
613 tautmsy, dtk, topflx, errstring, tauresx, tauresy, 1, do_molec_diff, &
614 compute_molec_diff, vd_lu_qdecomp)
616 if( errstring .ne. '' ) call wrf_error_fatal(errstring)
618 end if
620 !! Store updated tauresx, tauresy to use here on the next timestep
621 tauresx2d(i,j) = tauresx(ncol)
622 tauresy2d(i,j) = tauresy(ncol)
624 #ifdef MODAL_AERO
626 !! Add the explicit surface fluxes to the lowest layer
627 !! Modification : I should check whether this explicit adding is consistent with
628 !! the treatment of other tracers.
630 !The following code is commented out as the diffusion for the Aerosols and other species is handled by CAMMGMP and WRF_CHEM's dry_dep_driver subroutines
631 !tmp1(:ncol) = ztodt * gravit * rpdel8(:ncol,pver)
632 !do m = 1, ntot_amode
633 ! l = numptr_amode(m)
634 ! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l)
635 ! do lspec = 1, nspec_amode(m)
636 ! l = lmassptr_amode(lspec,m)
637 ! cloudtnd(:ncol,pver,l) = cloudtnd(:ncol,pver,l) + tmp1(:ncol) * cflx(:ncol,l)
638 ! enddo
639 !enddo
641 #endif
642 !! -------------------------------------------------------- !
643 !! Diagnostics and output writing after applying PBL scheme !
644 !! -------------------------------------------------------- !
645 sl(:ncol,:pver) = stnd(:ncol,:pver) - latvap * cloudtnd(:ncol,:pver,ixcldliq) &
646 - ( latvap + latice) * cloudtnd(:ncol,:pver,ixcldice)
647 qt(:ncol,:pver) = cloudtnd(:ncol,:pver,1) + cloudtnd(:ncol,:pver,ixcldliq) &
648 + cloudtnd(:ncol,:pver,ixcldice)
652 !! --------------------------------------------------------------- !
653 !! Convert the new profiles into vertical diffusion tendencies. !
654 !! Convert KE dissipative heat change into "temperature" tendency. !
655 !! --------------------------------------------------------------- !
657 slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
658 qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt
659 stnd(:ncol,:) = ( stnd(:ncol,:) - s8(:ncol,:) ) * rztodt
660 wind_tends(:ncol,:,1) = ( wind_tends(:ncol,:,1) - u8(:ncol,:) ) * rztodt
661 wind_tends(:ncol,:,2) = ( wind_tends(:ncol,:,2) - v8(:ncol,:) ) * rztodt
662 cloudtnd(:ncol,:pver,:) = ( cloudtnd(:ncol,:pver,:) - cloud(:ncol,:pver,:) ) * rztodt
664 !! ----------------------------------------------------------- !
665 !! In order to perform 'pseudo-conservative varible diffusion' !
666 !! perform the following two stages: !
667 !! !
668 !! I. Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0' !
669 !! (2) 'sten' by 'slten', and !
670 !! (3) 'qlten = qiten = 0' !
671 !! !
672 !! II. Apply 'positive_moisture' !
673 !! !
674 !! ----------------------------------------------------------- !
675 if( eddy_scheme .eq. 'diag_TKE' .and. do_pseudocon_diff ) then
676 cloudtnd(:ncol,:pver,1) = qtten(:ncol,:pver)
677 stnd(:ncol,:pver) = slten(:ncol,:pver)
678 cloudtnd(:ncol,:pver,ixcldliq) = 0._r8
679 cloudtnd(:ncol,:pver,ixcldice) = 0._r8
680 cloudtnd(:ncol,:pver,ixnumliq) = 0._r8
681 cloudtnd(:ncol,:pver,ixnumice) = 0._r8
683 do ips = 1, ncol
684 do k = 1, pver
685 qv_pro(ips,k) = cloud(ips,k,1) + cloudtnd(ips,k,1) * ztodt
686 ql_pro(ips,k) = cloud(ips,k,ixcldliq) + cloudtnd(ips,k,ixcldliq) * ztodt
687 qi_pro(ips,k) = cloud(ips,k,ixcldice) + cloudtnd(ips,k,ixcldice) * ztodt
688 s_pro(ips,k) = s8(ips,k) + stnd(ips,k) * ztodt
689 t_pro(ips,k) = t8(ips,k) + (1._r8/cpair)*stnd(ips,k) * ztodt
691 end do
692 end do
693 call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(2), qmin(3), &
694 pdel8(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), &
695 qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1), &
696 cloudtnd(:ncol,pver:1:-1,1), cloudtnd(:ncol,pver:1:-1,ixcldliq), &
697 cloudtnd(:ncol,pver:1:-1,ixcldice), stnd(:ncol,pver:1:-1) )
699 end if
701 !! ----------------------------------------------------------------- !
702 !! Re-calculate diagnostic output variables after vertical diffusion !
703 !! ----------------------------------------------------------------- !
704 qv_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,1) + cloudtnd(:ncol,:pver,1) * ztodt
705 ql_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldliq) + cloudtnd(:ncol,:pver,ixcldliq) * ztodt
706 qi_aft_PBL(:ncol,:pver) = cloud(:ncol,:pver,ixcldice) + cloudtnd(:ncol,:pver,ixcldice) * ztodt
708 s_aft_PBL(:ncol,:pver) = s8(:ncol,:pver) + stnd(:ncol,:pver) * ztodt
709 t_aftPBL(:ncol,:pver) = ( s_aft_PBL(:ncol,:pver) - gravit*zm8(:ncol,:pver) ) / cpair
711 u_aft_PBL(:ncol,:pver) = u8(:ncol,:pver) + wind_tends(:ncol,:pver,1) * ztodt
712 v_aft_PBL(:ncol,:pver) = v8(:ncol,:pver) + wind_tends(:ncol,:pver,2) * ztodt
714 call aqsat( t_aftPBL, pmid8, tem2, ftem, pcols, ncol, pver, 1, pver )
715 ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8
717 tten(:ncol,:pver) = ( t_aftPBL(:ncol,:pver) - t8(:ncol,:pver) ) * rztodt
718 rhten(:ncol,:pver) = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) ) * rztodt
721 !Post processing of the output from CAM's parameterization
722 do k=kts,kte
724 kflip = kte-k+1
726 rublten(i,k,j) = wind_tends(1,kflip,1)
727 rvblten(i,k,j) = wind_tends(1,kflip,2)
728 rthblten(i,k,j) = stnd(1,kflip)/cpair/exner8(1,kflip)
730 multFrc = 1._r8 + qv_curr(i,k,j)
732 rqvblten(i,k,j) = cloudtnd(1,kflip,1 ) * multFrc * multFrc
733 rqcblten(i,k,j) = cloudtnd(1,kflip,ixcldliq) * multFrc
734 rqiblten(i,k,j) = cloudtnd(1,kflip,ixcldice) * multFrc
735 !*Important* : ixnumliq is mixed in the dropmixnuc, therefore ixnumliq is NOT mixed here (ONLY if CAMMGMP is used for mp_physics)
736 rqniblten(i,k,j) = cloudtnd(1,kflip,ixnumice) * multFrc
738 !Diffusivity coefficients at the midpoints
739 kp1 = k + 1
740 end do
742 do k = kts,kte+1
743 kflip = kte - k + 2
745 tke_pbl(i,k,j) = tke8(1,kflip) !TKE at the interfaces
746 turbtype3d(i,k,j) = turbtype(1,kflip)
747 smaw3d(i,k,j) = smaw(1,kflip)
748 end do
752 kpbl2d(i,j) = kte - int(kpblh(1)) + 1 !int(kpblh(1))
753 pblh2d(i,j) = pblh( 1)
754 tpert2d(i,j) = tpert(1)
755 qpert2d(i,j) = qpert(1)
756 wpert2d(i,j) = wpert(1)
758 !End Post processing of the output from CAM
759 enddo !loop of i
760 enddo !loop of j
761 return
762 end subroutine camuwpbl
767 !-----------------------------------------------------------------------------------------
768 subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, &
769 dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten )
770 !! ------------------------------------------------------------------------------- !
771 !! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
772 !! force them to be larger than minimum value by (1) condensating water vapor !
773 !! into liquid or ice, and (2) by transporting water vapor from the very lower !
774 !! layer. '2._r8' is multiplied to the minimum values for safety. !
775 !! Update final state variables and tendencies associated with this correction. !
776 !! If any condensation happens, update (s,t) too. !
777 !! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
778 !! input tendencies. !
779 !! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer !
780 !! ------------------------------------------------------------------------------- !
781 !-----------------------------------------------------------------------------------------
782 implicit none
783 integer, intent(in) :: ncol, mkx
784 real(r8), intent(in) :: cp, xlv, xls
785 real(r8), intent(in) :: dt, qvmin, qlmin, qimin
786 real(r8), intent(in) :: dp(ncol,mkx)
787 real(r8), intent(inout) :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx)
788 real(r8), intent(inout) :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx)
789 integer i, k
790 real(r8) dql, dqi, dqv, sum, aa, dum
792 !! Modification : I should check whether this is exactly same as the one used in
793 !! shallow convection and cloud macrophysics.
795 do i = 1, ncol
796 do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface
797 dql = max(0._r8,1._r8*qlmin-ql(i,k))
798 dqi = max(0._r8,1._r8*qimin-qi(i,k))
799 qlten(i,k) = qlten(i,k) + dql/dt
800 qiten(i,k) = qiten(i,k) + dqi/dt
801 qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
802 sten(i,k) = sten(i,k) + xlv * (dql/dt) + xls * (dqi/dt)
803 ql(i,k) = ql(i,k) + dql
804 qi(i,k) = qi(i,k) + dqi
805 qv(i,k) = qv(i,k) - dql - dqi
806 s(i,k) = s(i,k) + xlv * dql + xls * dqi
807 t(i,k) = t(i,k) + (xlv * dql + xls * dqi)/cp
808 dqv = max(0._r8,1._r8*qvmin-qv(i,k))
809 qvten(i,k) = qvten(i,k) + dqv/dt
810 qv(i,k) = qv(i,k) + dqv
811 if( k .ne. 1 ) then
812 qv(i,k-1) = qv(i,k-1) - dqv*dp(i,k)/dp(i,k-1)
813 qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt
814 endif
815 qv(i,k) = max(qv(i,k),qvmin)
816 ql(i,k) = max(ql(i,k),qlmin)
817 qi(i,k) = max(qi(i,k),qimin)
818 end do
819 !! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
820 !! extracted from all the layers that has 'qv > 2*qvmin'. This fully
821 !! preserves column moisture.
822 if( dqv .gt. 1.e-20_r8 ) then
823 sum = 0._r8
824 do k = 1, mkx
825 if( qv(i,k) .gt. 2._r8*qvmin ) sum = sum + qv(i,k)*dp(i,k)
826 enddo
827 aa = dqv*dp(i,1)/max(1.e-20_r8,sum)
828 if( aa .lt. 0.5_r8 ) then
829 do k = 1, mkx
830 if( qv(i,k) .gt. 2._r8*qvmin ) then
831 dum = aa*qv(i,k)
832 qv(i,k) = qv(i,k) - dum
833 qvten(i,k) = qvten(i,k) - dum/dt
834 endif
835 enddo
836 else
837 write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion'
838 call wrf_message(iulog)
839 endif
840 endif
841 end do
842 return
844 end subroutine positive_moisture
848 !-----------------------------------------------------------------------------------------
849 subroutine camuwpblinit(rublten,rvblten,rthblten,rqvblten, &
850 restart,tke_pbl,is_CAMMGMP_used, &
851 ids,ide,jds,jde,kds,kde, &
852 ims,ime,jms,jme,kms,kme, &
853 its,ite,jts,jte,kts,kte)
854 !!------------------------------------------------------------------!
855 !! Initialization of time independent fields for vertical diffusion !
856 !! Calls initialization routines for subsidiary modules !
857 !!----------------------------------------------------------------- !
859 !This subroutine is based on vertical_diffusion_init of CAM. This subroutine
860 !initializes variables for vertical diffusion subroutine calls. The layout
861 !is kept similar to the original CAM subroutine to facillitate future adaptations.
862 !Called by : physics_init.F
863 !Calls: vd_register, cnst_get_ind, wrf_error_fatal,init_eddy_diff,init_tms,init_vdiff
864 !-----------------------------------------------------------------------------------------
866 use eddy_diff, only : init_eddy_diff
867 use molec_diff, only : init_molec_diff
868 use trb_mtn_stress, only : init_tms
869 use diffusion_solver, only : init_vdiff, vdiff_select
870 use constituents, only : cnst_get_ind, cnst_get_type_byind, cnst_name
871 use module_cam_support, only : masterproc
872 use module_model_constants, only : epsq2
873 #ifdef MODAL_AERO
874 use modal_aero_data
875 #endif
877 implicit none
879 !-------------------------------------------------------------------------------------!
880 !Input and output variables !
881 !-------------------------------------------------------------------------------------!
882 logical,intent(in) :: restart,is_CAMMGMP_used
883 integer,intent(in) :: ids,ide,jds,jde,kds,kde
884 integer,intent(in) :: ims,ime,jms,jme,kms,kme
885 integer,intent(in) :: its,ite,jts,jte,kts,kte
887 real,dimension(ims:ime,kms:kme,jms:jme),intent(inout) :: &
888 rublten, rvblten, rthblten, rqvblten
889 real,dimension(ims:ime,kms:kme,jms:jme),intent(out) :: TKE_PBL
891 !-------------------------------------------------------------------------------------!
892 !Local Variables !
893 !-------------------------------------------------------------------------------------!
894 integer :: i,j,k,itf,jtf,ktf
895 integer :: ntop_eddy !! Top interface level to which eddy vertical diffusion is applied ( = 1 )
896 integer :: nbot_eddy !! Bottom interface level to which eddy vertical diffusion is applied ( = pver )
897 integer :: ntop_molec !! Top interface level to which molecular vertical diffusion is applied ( = 1 )
898 integer :: nbot_molec !! Bottom interface level to which molecular vertical diffusion is applied
899 character(128) :: errstring !! Error status for init_vdiff
900 real(r8) :: hypm(kte) !! reference state midpoint pressures
901 #ifdef MODAL_AERO
902 integer :: m, l
903 #endif
906 jtf = min(jte,jde-1)
907 ktf = min(kte,kde-1)
908 itf = min(ite,ide-1)
910 !Map CAM veritcal level variables
911 pver = ktf - kts + 1
912 pverp = pver + 1
914 !Initialize flags and add constituents
915 call vd_register()
917 !! ----------------------------------------------------------------- !
918 !! Get indices of cloud liquid and ice within the constituents array !
919 !! ----------------------------------------------------------------- !
921 call cnst_get_ind( 'CLDLIQ', ixcldliq )
922 call cnst_get_ind( 'CLDICE', ixcldice )
923 if( microp_scheme .eq. 'MG' ) then
924 call cnst_get_ind( 'NUMLIQ', ixnumliq )
925 call cnst_get_ind( 'NUMICE', ixnumice )
926 endif
928 if (masterproc) then
929 write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)'
930 call wrf_debug(1,iulog)
931 end if
933 !! ---------------------------------------------------------------------------------------- !
934 !! Initialize molecular diffusivity module !
935 !! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km (.1 Pa). !
936 !! Note that computing molecular diffusivities is a trivial expense, but constituent !
937 !! diffusivities depend on their molecular weights. Decomposing the diffusion matric !
938 !! for each constituent is a needless expense unless the diffusivity is significant. !
939 !! ---------------------------------------------------------------------------------------- !
941 ntop_molec = 1 !! Should always be 1
942 nbot_molec = 0 !! Should be set below about 70 km
944 !! ---------------------------------- !
945 !! Initialize eddy diffusivity module !
946 !! ---------------------------------- !
948 ntop_eddy = 1 !! No reason not to make this 1, if > 1, must be <= nbot_molec
949 nbot_eddy = pver !! Should always be pver
951 if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy
952 call wrf_debug(1,iulog)
954 select case ( eddy_scheme )
955 case ( 'diag_TKE' )
956 if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme'
957 call wrf_debug(1,iulog)
958 if( masterproc ) write(iulog,*) 'UW Moist Turbulence Scheme by Bretherton and Park'
959 call wrf_debug(1,iulog)
960 !! Check compatibility of eddy and shallow scheme
961 if( shallow_scheme .ne. 'UW' ) then
962 write(iulog,*) 'ERROR: shallow convection scheme ', shallow_scheme,' is incompatible with eddy scheme ', eddy_scheme
963 call wrf_message(iulog)
964 call wrf_error_fatal( 'convect_shallow_init: shallow_scheme and eddy_scheme are incompatible' )
965 endif
967 call init_eddy_diff( r8, pver, gravit, cpair, rair, zvir, latvap, latice, &
968 ntop_eddy, nbot_eddy, hypm, karman )
970 if( masterproc ) write(iulog,*) 'vertical_diffusion: nturb, ntop_eddy, nbot_eddy ', nturb, ntop_eddy, nbot_eddy
971 call wrf_debug(1,iulog)
972 end select
974 !!-------------------------------------------------------------------------------------!
975 !! The vertical diffusion solver must operate !
976 !! over the full range of molecular and eddy diffusion !
977 !!-------------------------------------------------------------------------------------!
979 ntop = min(ntop_molec,ntop_eddy)
980 nbot = max(nbot_molec,nbot_eddy)
982 !! ------------------------------------------- !
983 !! Initialize turbulent mountain stress module !
984 !! ------------------------------------------- !
986 if( do_tms ) then
987 call init_tms( r8, tms_orocnst, tms_z0fac, karman, gravit, rair )
988 if (masterproc) then
989 write(iulog,*)'Using turbulent mountain stress module'
990 call wrf_message(iulog)
991 write(iulog,*)' tms_orocnst = ',tms_orocnst
992 call wrf_message(iulog)
993 end if
994 endif
996 !! ---------------------------------- !
997 !! Initialize diffusion solver module !
998 !! ---------------------------------- !
1000 call init_vdiff( r8, pcnst, rair, gravit, fieldlist_wet, fieldlist_dry, errstring )
1001 if( errstring .ne. '' ) call wrf_error_fatal( errstring )
1003 !!-------------------------------------------------------------------------------------
1004 !! Use fieldlist_wet to select the fields which will be diffused using moist mixing ratios ( all by default )
1005 !! Use fieldlist_dry to select the fields which will be diffused using dry mixing ratios.
1006 !!-------------------------------------------------------------------------------------
1008 if( vdiff_select( fieldlist_wet, 'u' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'u' ) )
1009 if( vdiff_select( fieldlist_wet, 'v' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'v' ) )
1010 if( vdiff_select( fieldlist_wet, 's' ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 's' ) )
1012 #ifdef MODAL_AERO
1013 !Get index for the liquid number concentration (#/kg)
1014 call cnst_get_ind( 'NUMLIQ', ixndrop )
1015 #endif
1017 do k = 1, pcnst
1019 #ifdef MODAL_AERO
1020 !! Do not diffuse droplet number - treated in dropmixnuc
1022 !BSINGH:01/31/2013: The following if condition is modified from its original form[which is:if( k == ixndrop ) cycle]
1023 !to allow the PBL scheme to mix liquid number when CAMMGMP scheme is not being used
1025 if(is_CAMMGMP_used) then!BSINGH:01/31/2013: Modified this if cond. see comment above
1026 if( k == ixndrop ) cycle
1027 endif
1028 !BSINGH-08/31/2012: In WRF, Physics init is called before chem init therefore
1029 !aerosol species are not yet added. Following condition, skip all other constituents
1030 !(which will be added in chem init) except first five (vapor, cld liq, cld ice, liq num and ice num)
1031 if( k > 5 ) cycle
1034 !The Following do-loop is commented out as Aerosols and other species are diffused in CAMMGMP and WRF_CHEM's dry_dep_driver subroutines
1036 !! Don't diffuse aerosol - treated in dropmixnuc
1037 !do m = 1, ntot_amode
1038 ! if( k == numptr_amode(m) ) cycle
1039 !! if( k == numptrcw_amode(m) ) go to 20
1040 ! do l = 1, nspec_amode(m)
1041 ! if( k == lmassptr_amode(l,m) )cycle
1042 !! if( k == lmassptrcw_amode(l,m) ) go to 20
1043 ! enddo
1044 !! if( k == lwaterptr_amode(m) ) go to 20
1045 !enddo
1046 #endif
1047 !endif
1048 if( cnst_get_type_byind(k) .eq. 'wet' ) then
1049 if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_wet, 'q', k ) )
1050 else
1051 !The DRY constituents are not diffused currently in PBL, they are diffused in WRF_CHEM's dry_dep_driver subroutine. Therefore the following line is commented out
1052 !if( vdiff_select( fieldlist_dry, 'q', k ) .ne. '' ) call wrf_error_fatal( vdiff_select( fieldlist_dry, 'q', k ) )
1053 endif
1054 enddo
1056 !Initialize the tendencies
1057 jtf=min0(jte,jde-1)
1058 ktf=min0(kte,kde-1)
1059 itf=min0(ite,ide-1)
1061 if(.not.restart)then
1062 do j = jts , jtf
1063 do k = kts , ktf
1064 do i = its , itf
1065 tke_pbl(i,k,j) = epsq2
1066 rublten(i,k,j) = 0.
1067 rvblten(i,k,j) = 0.
1068 rthblten(i,k,j) = 0.
1069 rqvblten(i,k,j) = 0.
1070 enddo
1071 enddo
1072 enddo
1073 endif
1074 end subroutine camuwpblinit
1078 !-----------------------------------------------------------------------------------------
1079 subroutine vd_register()
1080 !
1081 !This subroutine is based on the vd_register subroutine of CAM.
1082 !-----------------------------------------------------------------------------------------
1084 !Set flags for the PBL scheme (these flags are obtained using phys_getopts in CAM)
1085 microp_scheme = 'MG' !Used for adding constituents
1086 eddy_scheme = 'diag_TKE' !Used for calling eddy scheme
1088 !The following flag is deliberately set to UW for now.
1089 shallow_scheme = 'UW' !Eddy scheme is ONLY compaticle with 'UW' shallow scheme;
1091 do_tms = .false. !To include stresses due to orography
1092 tms_orocnst = 1._r8 !Orography constant
1094 end subroutine vd_register
1096 end module module_bl_camuwpbl_driver