Line data Source code
1 : module vertical_diffusion
2 :
3 : !----------------------------------------------------------------------------------------------------- !
4 : ! Module to compute vertical diffusion of momentum, moisture, trace constituents !
5 : ! and static energy. Separate modules compute !
6 : ! 1. stresses associated with turbulent flow over orography !
7 : ! ( turbulent mountain stress ) !
8 : ! 2. eddy diffusivities, including nonlocal tranport terms !
9 : ! 3. molecular diffusivities !
10 : ! Lastly, a implicit diffusion solver is called, and tendencies retrieved by !
11 : ! differencing the diffused and initial states. !
12 : ! !
13 : ! Calling sequence: !
14 : ! !
15 : ! vertical_diffusion_init Initializes vertical diffustion constants and modules !
16 : ! init_molec_diff Initializes molecular diffusivity module !
17 : ! init_eddy_diff Initializes eddy diffusivity module (includes PBL) !
18 : ! init_tms Initializes turbulent mountain stress module !
19 : ! init_vdiff Initializes diffusion solver module !
20 : ! vertical_diffusion_ts_init Time step initialization (only used for upper boundary condition) !
21 : ! vertical_diffusion_tend Computes vertical diffusion tendencies !
22 : ! compute_tms Computes turbulent mountain stresses !
23 : ! compute_eddy_diff Computes eddy diffusivities and countergradient terms !
24 : ! compute_vdiff Solves vertical diffusion equations, including molecular diffusivities !
25 : ! !
26 : !----------------------------------------------------------------------------------------------------- !
27 : ! Some notes on refactoring changes made in 2015, which were not quite finished. !
28 : ! !
29 : ! - eddy_diff_tend should really only have state, pbuf, and cam_in as inputs. The process of !
30 : ! removing these arguments, and referring to pbuf fields instead, is not complete. !
31 : ! !
32 : ! - compute_vdiff was intended to be split up into three components: !
33 : ! !
34 : ! 1. Diffusion of winds and heat ("U", "V", and "S" in the fieldlist object). !
35 : ! !
36 : ! 2. Turbulent diffusion of a single constituent !
37 : ! !
38 : ! 3. Molecular diffusion of a single constituent !
39 : ! !
40 : ! This reorganization would allow the three resulting functions to each use a simpler interface !
41 : ! than the current combined version, and possibly also remove the need to use the fieldlist !
42 : ! object at all. !
43 : ! !
44 : ! - The conditionals controlled by "do_pbl_diags" are somewhat scattered. It might be better to !
45 : ! pull out these diagnostic calculations and outfld calls into separate functions. !
46 : ! !
47 : !---------------------------Code history-------------------------------------------------------------- !
48 : ! J. Rosinski : Jun. 1992 !
49 : ! J. McCaa : Sep. 2004 !
50 : ! S. Park : Aug. 2006, Dec. 2008. Jan. 2010 !
51 : !----------------------------------------------------------------------------------------------------- !
52 :
53 : use shr_kind_mod, only : r8 => shr_kind_r8, i4=> shr_kind_i4
54 : use ppgrid, only : pcols, pver, pverp
55 : use constituents, only : pcnst
56 : use diffusion_solver, only : vdiff_selector
57 : use cam_abortutils, only : endrun
58 : use error_messages, only : handle_errmsg
59 : use physconst, only : &
60 : cpair , & ! Specific heat of dry air
61 : gravit , & ! Acceleration due to gravity
62 : rair , & ! Gas constant for dry air
63 : zvir , & ! rh2o/rair - 1
64 : latvap , & ! Latent heat of vaporization
65 : latice , & ! Latent heat of fusion
66 : karman , & ! von Karman constant
67 : mwdry , & ! Molecular weight of dry air
68 : avogad ! Avogadro's number
69 : use cam_history, only : fieldname_len
70 : use perf_mod
71 : use cam_logfile, only : iulog
72 : use ref_pres, only : do_molec_diff, nbot_molec
73 : use phys_control, only : phys_getopts
74 : use time_manager, only : is_first_step
75 : implicit none
76 : private
77 : save
78 :
79 : ! ----------------- !
80 : ! Public interfaces !
81 : ! ----------------- !
82 :
83 : public vd_readnl
84 : public vd_register ! Register multi-time-level variables with physics buffer
85 : public vertical_diffusion_init ! Initialization
86 : public vertical_diffusion_ts_init ! Time step initialization (only used for upper boundary condition)
87 : public vertical_diffusion_tend ! Full vertical diffusion routine
88 :
89 : ! ------------ !
90 : ! Private data !
91 : ! ------------ !
92 :
93 : character(len=16) :: eddy_scheme ! Default set in phys_control.F90, use namelist to change
94 : ! 'HB' = Holtslag and Boville (default)
95 : ! 'HBR' = Holtslag and Boville and Rash
96 : ! 'diag_TKE' = Bretherton and Park ( UW Moist Turbulence Scheme )
97 : logical, parameter :: wstarent = .true. ! Use wstar (.true.) or TKE (.false.) entrainment closure
98 : ! ( when 'diag_TKE' scheme is selected )
99 : logical :: do_pseudocon_diff = .false. ! If .true., do pseudo-conservative variables diffusion
100 :
101 : character(len=16) :: shallow_scheme ! Shallow convection scheme
102 :
103 : type(vdiff_selector) :: fieldlist_wet ! Logical switches for moist mixing ratio diffusion
104 : type(vdiff_selector) :: fieldlist_dry ! Logical switches for dry mixing ratio diffusion
105 : type(vdiff_selector) :: fieldlist_molec ! Logical switches for molecular diffusion
106 : integer :: tke_idx, kvh_idx, kvm_idx ! TKE and eddy diffusivity indices for fields in the physics buffer
107 : integer :: kvt_idx ! Index for kinematic molecular conductivity
108 : integer :: tauresx_idx, tauresy_idx ! Redisual stress for implicit surface stress
109 :
110 : character(len=fieldname_len) :: vdiffnam(pcnst) ! Names of vertical diffusion tendencies
111 : integer :: ixcldice, ixcldliq ! Constituent indices for cloud liquid and ice water
112 : integer :: ixnumice, ixnumliq
113 :
114 : integer :: pblh_idx, tpert_idx, qpert_idx
115 :
116 : ! pbuf fields for unicon
117 : integer :: qtl_flx_idx = -1 ! for use in cloud macrophysics when UNICON is on
118 : integer :: qti_flx_idx = -1 ! for use in cloud macrophysics when UNICON is on
119 :
120 : ! pbuf fields for tms
121 : integer :: ksrftms_idx = -1
122 : integer :: tautmsx_idx = -1
123 : integer :: tautmsy_idx = -1
124 :
125 : ! pbuf fields for blj (Beljaars)
126 : integer :: dragblj_idx = -1
127 : integer :: taubljx_idx = -1
128 : integer :: taubljy_idx = -1
129 :
130 : ! pbuf field for clubb top above which HB (Holtslag Boville) scheme may be enabled
131 : integer :: clubbtop_idx = -1
132 :
133 : logical :: diff_cnsrv_mass_check ! do mass conservation check
134 : logical :: do_iss ! switch for implicit turbulent surface stress
135 : logical :: prog_modal_aero = .false. ! set true if prognostic modal aerosols are present
136 : integer :: pmam_ncnst = 0 ! number of prognostic modal aerosol constituents
137 : integer, allocatable :: pmam_cnst_idx(:) ! constituent indices of prognostic modal aerosols
138 :
139 : logical :: do_pbl_diags = .false.
140 : logical :: waccmx_mode = .false.
141 : logical :: do_hb_above_clubb = .false.
142 :
143 : real(r8),allocatable :: kvm_sponge(:)
144 :
145 : contains
146 :
147 : ! =============================================================================== !
148 : ! !
149 : ! =============================================================================== !
150 1536 : subroutine vd_readnl(nlfile)
151 :
152 : use namelist_utils, only: find_group_name
153 : use units, only: getunit, freeunit
154 : use spmd_utils, only: masterproc, masterprocid, mpi_logical, mpicom
155 : use shr_log_mod, only: errMsg => shr_log_errMsg
156 : use trb_mtn_stress_cam, only: trb_mtn_stress_readnl
157 : use beljaars_drag_cam, only: beljaars_drag_readnl
158 : use eddy_diff_cam, only: eddy_diff_readnl
159 :
160 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
161 :
162 : ! Local variables
163 : integer :: unitn, ierr
164 : character(len=*), parameter :: subname = 'vd_readnl'
165 :
166 : namelist /vert_diff_nl/ diff_cnsrv_mass_check, do_iss
167 : !-----------------------------------------------------------------------------
168 :
169 1536 : if (masterproc) then
170 2 : unitn = getunit()
171 2 : open( unitn, file=trim(nlfile), status='old' )
172 2 : call find_group_name(unitn, 'vert_diff_nl', status=ierr)
173 2 : if (ierr == 0) then
174 2 : read(unitn, vert_diff_nl, iostat=ierr)
175 2 : if (ierr /= 0) then
176 0 : call endrun(subname // ':: ERROR reading namelist')
177 : end if
178 : end if
179 2 : close(unitn)
180 2 : call freeunit(unitn)
181 : end if
182 :
183 1536 : call mpi_bcast(diff_cnsrv_mass_check, 1, mpi_logical, masterprocid, mpicom, ierr)
184 1536 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
185 1536 : call mpi_bcast(do_iss, 1, mpi_logical, masterprocid, mpicom, ierr)
186 1536 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
187 :
188 : ! Get eddy_scheme setting from phys_control.
189 : call phys_getopts( eddy_scheme_out = eddy_scheme, &
190 1536 : shallow_scheme_out = shallow_scheme )
191 :
192 : ! TMS reads its own namelist.
193 1536 : call trb_mtn_stress_readnl(nlfile)
194 :
195 : ! Beljaars reads its own namelist.
196 1536 : call beljaars_drag_readnl(nlfile)
197 :
198 1536 : if (eddy_scheme == 'diag_TKE' .or. eddy_scheme == 'SPCAM_m2005' ) call eddy_diff_readnl(nlfile)
199 :
200 1536 : end subroutine vd_readnl
201 :
202 : ! =============================================================================== !
203 : ! !
204 : ! =============================================================================== !
205 :
206 1536 : subroutine vd_register()
207 :
208 : !------------------------------------------------ !
209 : ! Register physics buffer fields and constituents !
210 : !------------------------------------------------ !
211 :
212 1536 : use physics_buffer, only : pbuf_add_field, dtype_r8, dtype_i4
213 : use trb_mtn_stress_cam, only : trb_mtn_stress_register
214 : use beljaars_drag_cam, only : beljaars_drag_register
215 : use eddy_diff_cam, only : eddy_diff_register
216 :
217 : ! Add fields to physics buffer
218 :
219 : ! kvt is used by gw_drag. only needs physpkg scope.
220 1536 : call pbuf_add_field('kvt', 'physpkg', dtype_r8, (/pcols,pverp/), kvt_idx)
221 :
222 :
223 1536 : if (eddy_scheme /= 'CLUBB_SGS') then
224 0 : call pbuf_add_field('kvh', 'global', dtype_r8, (/pcols, pverp/), kvh_idx)
225 : end if
226 :
227 1536 : call pbuf_add_field('kvm', 'global', dtype_r8, (/pcols, pverp/), kvm_idx )
228 1536 : call pbuf_add_field('pblh', 'global', dtype_r8, (/pcols/), pblh_idx)
229 1536 : call pbuf_add_field('tke', 'global', dtype_r8, (/pcols, pverp/), tke_idx)
230 :
231 1536 : call pbuf_add_field('tauresx', 'global', dtype_r8, (/pcols/), tauresx_idx)
232 1536 : call pbuf_add_field('tauresy', 'global', dtype_r8, (/pcols/), tauresy_idx)
233 :
234 1536 : call pbuf_add_field('tpert', 'global', dtype_r8, (/pcols/), tpert_idx)
235 1536 : call pbuf_add_field('qpert', 'global', dtype_r8, (/pcols,pcnst/), qpert_idx)
236 :
237 1536 : if (trim(shallow_scheme) == 'UNICON') then
238 0 : call pbuf_add_field('qtl_flx', 'global', dtype_r8, (/pcols, pverp/), qtl_flx_idx)
239 0 : call pbuf_add_field('qti_flx', 'global', dtype_r8, (/pcols, pverp/), qti_flx_idx)
240 : end if
241 :
242 : ! diag_TKE fields
243 1536 : if (eddy_scheme == 'diag_TKE' .or. eddy_scheme == 'SPCAM_m2005') then
244 0 : call eddy_diff_register()
245 : end if
246 :
247 : ! TMS fields
248 1536 : call trb_mtn_stress_register()
249 :
250 : ! Beljaars fields
251 1536 : call beljaars_drag_register()
252 :
253 1536 : end subroutine vd_register
254 :
255 : ! =============================================================================== !
256 : ! !
257 : ! =============================================================================== !
258 :
259 1536 : subroutine vertical_diffusion_init(pbuf2d)
260 :
261 : !------------------------------------------------------------------!
262 : ! Initialization of time independent fields for vertical diffusion !
263 : ! Calls initialization routines for subsidiary modules !
264 : !----------------------------------------------------------------- !
265 :
266 1536 : use cam_history, only : addfld, add_default, horiz_only
267 : use cam_history, only : register_vector_field
268 : use eddy_diff_cam, only : eddy_diff_init
269 : use hb_diff, only : init_hb_diff
270 : use molec_diff, only : init_molec_diff
271 : use diffusion_solver, only : init_vdiff, new_fieldlist_vdiff, vdiff_select
272 : use constituents, only : cnst_get_ind, cnst_get_type_byind, cnst_name, cnst_get_molec_byind
273 : use spmd_utils, only : masterproc
274 : use ref_pres, only : press_lim_idx, pref_mid
275 : use physics_buffer, only : pbuf_set_field, pbuf_get_index, physics_buffer_desc
276 : use rad_constituents, only : rad_cnst_get_info, rad_cnst_get_mode_num_idx, &
277 : rad_cnst_get_mam_mmr_idx
278 : use trb_mtn_stress_cam,only : trb_mtn_stress_init
279 : use beljaars_drag_cam, only : beljaars_drag_init
280 : use upper_bc, only : ubc_init
281 : use phys_control, only : waccmx_is, fv_am_correction
282 : use ref_pres, only : ptop_ref
283 :
284 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
285 : character(128) :: errstring ! Error status for init_vdiff
286 : integer :: ntop_eddy ! Top interface level to which eddy vertical diffusion is applied ( = 1 )
287 : integer :: nbot_eddy ! Bottom interface level to which eddy vertical diffusion is applied ( = pver )
288 : integer :: k ! Vertical loop index
289 :
290 : real(r8), parameter :: ntop_eddy_pres = 1.e-7_r8 ! Pressure below which eddy diffusion is not done in WACCM-X. (Pa)
291 :
292 : integer :: im, l, m, nmodes, nspec, ierr
293 :
294 : logical :: history_amwg ! output the variables used by the AMWG diag package
295 : logical :: history_eddy ! output the eddy variables
296 : logical :: history_budget ! Output tendencies and state variables for CAM4 T, qv, ql, qi
297 : integer :: history_budget_histfile_num ! output history file number for budget fields
298 : logical :: history_waccm ! output variables of interest for WACCM runs
299 :
300 : !
301 : ! add sponge layer vertical diffusion
302 : !
303 1536 : if (ptop_ref>1e-1_r8.and.ptop_ref<100.0_r8) then
304 : !
305 : ! CAM7 FMT (but not CAM6 top (~225 Pa) or CAM7 low top or lower)
306 : !
307 1536 : allocate(kvm_sponge(4), stat=ierr)
308 1536 : if( ierr /= 0 ) then
309 0 : write(iulog,*) 'vertical_diffusion_init: kvm_sponge allocation error = ',ierr
310 0 : call endrun('vertical_diffusion_init: failed to allocate kvm_sponge array')
311 : end if
312 1536 : kvm_sponge(1) = 2E6_r8
313 1536 : kvm_sponge(2) = 2E6_r8
314 1536 : kvm_sponge(3) = 0.5E6_r8
315 1536 : kvm_sponge(4) = 0.1E6_r8
316 0 : else if (ptop_ref>1e-4_r8) then
317 : !
318 : ! WACCM and WACCM-x
319 : !
320 0 : allocate(kvm_sponge(6), stat=ierr)
321 0 : if( ierr /= 0 ) then
322 0 : write(iulog,*) 'vertical_diffusion_init: kvm_sponge allocation error = ',ierr
323 0 : call endrun('vertical_diffusion_init: failed to allocate kvm_sponge array')
324 : end if
325 0 : kvm_sponge(1) = 2E6_r8
326 0 : kvm_sponge(2) = 2E6_r8
327 0 : kvm_sponge(3) = 1.5E6_r8
328 0 : kvm_sponge(4) = 1.0E6_r8
329 0 : kvm_sponge(5) = 0.5E6_r8
330 0 : kvm_sponge(6) = 0.1E6_r8
331 : end if
332 :
333 1536 : if (masterproc) then
334 2 : write(iulog,*)'Initializing vertical diffusion (vertical_diffusion_init)'
335 2 : if (allocated(kvm_sponge)) then
336 2 : write(iulog,*)'Artificial sponge layer vertical diffusion added:'
337 10 : do k=1,size(kvm_sponge(:),1)
338 8 : write(iulog,'(a44,i2,a17,e7.2,a8)') 'vertical diffusion coefficient at interface',k,' is increased by ', &
339 18 : kvm_sponge(k),' m2 s-2'
340 : end do
341 : end if !allocated
342 : end if
343 :
344 : ! Check to see if WACCM-X is on (currently we don't care whether the
345 : ! ionosphere is on or not, since this neutral diffusion code is the
346 : ! same either way).
347 1536 : waccmx_mode = waccmx_is('ionosphere') .or. waccmx_is('neutral')
348 :
349 : ! ----------------------------------------------------------------- !
350 : ! Get indices of cloud liquid and ice within the constituents array !
351 : ! ----------------------------------------------------------------- !
352 :
353 1536 : call cnst_get_ind( 'CLDLIQ', ixcldliq )
354 1536 : call cnst_get_ind( 'CLDICE', ixcldice )
355 : ! These are optional; with the CAM4 microphysics, there are no number
356 : ! constituents.
357 1536 : call cnst_get_ind( 'NUMLIQ', ixnumliq, abort=.false. )
358 1536 : call cnst_get_ind( 'NUMICE', ixnumice, abort=.false. )
359 :
360 : ! prog_modal_aero determines whether prognostic modal aerosols are present in the run.
361 1536 : call phys_getopts(prog_modal_aero_out=prog_modal_aero)
362 1536 : if (prog_modal_aero) then
363 :
364 : ! Get the constituent indices of the number and mass mixing ratios of the modal
365 : ! aerosols.
366 : !
367 : ! N.B. - This implementation assumes that the prognostic modal aerosols are
368 : ! impacting the climate calculation (i.e., can get info from list 0).
369 : !
370 :
371 : ! First need total number of mam constituents
372 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
373 7680 : do m = 1, nmodes
374 6144 : call rad_cnst_get_info(0, m, nspec=nspec)
375 7680 : pmam_ncnst = pmam_ncnst + 1 + nspec
376 : end do
377 :
378 4608 : allocate(pmam_cnst_idx(pmam_ncnst))
379 :
380 : ! Get the constituent indicies
381 1536 : im = 1
382 7680 : do m = 1, nmodes
383 6144 : call rad_cnst_get_mode_num_idx(m, pmam_cnst_idx(im))
384 6144 : im = im + 1
385 6144 : call rad_cnst_get_info(0, m, nspec=nspec)
386 36864 : do l = 1, nspec
387 23040 : call rad_cnst_get_mam_mmr_idx(m, l, pmam_cnst_idx(im))
388 29184 : im = im + 1
389 : end do
390 : end do
391 : end if
392 :
393 : ! Initialize upper boundary condition module
394 :
395 1536 : call ubc_init()
396 :
397 : ! ---------------------------------------------------------------------------------------- !
398 : ! Initialize molecular diffusivity module !
399 : ! Note that computing molecular diffusivities is a trivial expense, but constituent !
400 : ! diffusivities depend on their molecular weights. Decomposing the diffusion matrix !
401 : ! for each constituent is a needless expense unless the diffusivity is significant. !
402 : ! ---------------------------------------------------------------------------------------- !
403 :
404 : !----------------------------------------------------------------------------------------
405 : ! Initialize molecular diffusion and get top and bottom molecular diffusion limits
406 : !----------------------------------------------------------------------------------------
407 :
408 1536 : if( do_molec_diff ) then
409 : call init_molec_diff( r8, pcnst, mwdry, avogad, &
410 0 : errstring)
411 :
412 0 : call handle_errmsg(errstring, subname="init_molec_diff")
413 :
414 0 : call addfld( 'TTPXMLC', horiz_only, 'A', 'K/S', 'Top interf. temp. flux: molec. viscosity' )
415 0 : if( masterproc ) write(iulog,fmt='(a,i3,5x,a,i3)') 'NBOT_MOLEC =', nbot_molec
416 : end if
417 :
418 : ! ---------------------------------- !
419 : ! Initialize eddy diffusivity module !
420 : ! ---------------------------------- !
421 :
422 : ! ntop_eddy must be 1 or <= nbot_molec
423 : ! Currently, it is always 1 except for WACCM-X.
424 1536 : if ( waccmx_mode ) then
425 0 : ntop_eddy = press_lim_idx(ntop_eddy_pres, top=.true.)
426 : else
427 1536 : ntop_eddy = 1
428 : end if
429 1536 : nbot_eddy = pver
430 :
431 1536 : if (masterproc) write(iulog, fmt='(a,i3,5x,a,i3)') 'NTOP_EDDY =', ntop_eddy, 'NBOT_EDDY =', nbot_eddy
432 :
433 1536 : call phys_getopts(do_hb_above_clubb_out=do_hb_above_clubb)
434 :
435 0 : select case ( eddy_scheme )
436 : case ( 'diag_TKE', 'SPCAM_m2005' )
437 0 : if( masterproc ) write(iulog,*) &
438 0 : 'vertical_diffusion_init: eddy_diffusivity scheme: UW Moist Turbulence Scheme by Bretherton and Park'
439 0 : call eddy_diff_init(pbuf2d, ntop_eddy, nbot_eddy)
440 : case ( 'HB', 'HBR', 'SPCAM_sam1mom')
441 0 : if( masterproc ) write(iulog,*) 'vertical_diffusion_init: eddy_diffusivity scheme: Holtslag and Boville'
442 : call init_hb_diff(gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, &
443 0 : karman, eddy_scheme)
444 0 : call addfld('HB_ri', (/ 'lev' /), 'A', 'no', 'Richardson Number (HB Scheme), I' )
445 : case ( 'CLUBB_SGS' )
446 1536 : do_pbl_diags = .true.
447 1536 : call init_hb_diff(gravit, cpair, ntop_eddy, nbot_eddy, pref_mid, karman, eddy_scheme)
448 : !
449 : ! run HB scheme where CLUBB is not active when running cam7 or cam6 physics
450 : ! else init_hb_diff is called just for diagnostic purposes
451 : !
452 3072 : if (do_hb_above_clubb) then
453 1536 : if( masterproc ) then
454 2 : write(iulog,*) 'vertical_diffusion_init: '
455 2 : write(iulog,*) 'eddy_diffusivity scheme where CLUBB is not active: Holtslag and Boville'
456 : end if
457 3072 : call addfld('HB_ri', (/ 'lev' /), 'A', 'no', 'Richardson Number (HB Scheme), I' )
458 : end if
459 : end select
460 :
461 : ! ------------------------------------------- !
462 : ! Initialize turbulent mountain stress module !
463 : ! ------------------------------------------- !
464 :
465 1536 : call trb_mtn_stress_init()
466 :
467 : ! ----------------------------------- !
468 : ! Initialize Beljaars SGO drag module !
469 : ! ----------------------------------- !
470 :
471 1536 : call beljaars_drag_init()
472 :
473 : ! ---------------------------------- !
474 : ! Initialize diffusion solver module !
475 : ! ---------------------------------- !
476 :
477 1536 : call init_vdiff(r8, iulog, rair, cpair, gravit, do_iss, fv_am_correction, errstring)
478 1536 : call handle_errmsg(errstring, subname="init_vdiff")
479 :
480 : ! Use fieldlist_wet to select the fields which will be diffused using moist mixing ratios ( all by default )
481 : ! Use fieldlist_dry to select the fields which will be diffused using dry mixing ratios.
482 :
483 1536 : fieldlist_wet = new_fieldlist_vdiff( pcnst)
484 1536 : fieldlist_dry = new_fieldlist_vdiff( pcnst)
485 1536 : fieldlist_molec = new_fieldlist_vdiff( pcnst)
486 :
487 1536 : if( vdiff_select( fieldlist_wet, 'u' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'u' ) )
488 1536 : if( vdiff_select( fieldlist_wet, 'v' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'v' ) )
489 1536 : if( vdiff_select( fieldlist_wet, 's' ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 's' ) )
490 :
491 64512 : constit_loop: do k = 1, pcnst
492 :
493 62976 : if (prog_modal_aero) then
494 : ! Do not diffuse droplet number - treated in dropmixnuc
495 62976 : if (k == ixnumliq) cycle constit_loop
496 : ! Don't diffuse modal aerosol - treated in dropmixnuc
497 936960 : do m = 1, pmam_ncnst
498 936960 : if (k == pmam_cnst_idx(m)) cycle constit_loop
499 : enddo
500 : end if
501 :
502 : ! Convert all constituents to wet before doing diffusion.
503 32256 : if( vdiff_select( fieldlist_wet, 'q', k ) .ne. '' ) call endrun( vdiff_select( fieldlist_wet, 'q', k ) )
504 :
505 : ! ----------------------------------------------- !
506 : ! Select constituents for molecular diffusion !
507 : ! ----------------------------------------------- !
508 33792 : if ( cnst_get_molec_byind(k) .eq. 'minor' ) then
509 95232 : if( vdiff_select(fieldlist_molec,'q',k) .ne. '' ) call endrun( vdiff_select( fieldlist_molec,'q',k ) )
510 : endif
511 :
512 : end do constit_loop
513 :
514 : ! ------------------------ !
515 : ! Diagnostic output fields !
516 : ! ------------------------ !
517 :
518 64512 : do k = 1, pcnst
519 62976 : vdiffnam(k) = 'VD'//cnst_name(k)
520 62976 : if( k == 1 ) vdiffnam(k) = 'VD01' !**** compatibility with old code ****
521 127488 : call addfld( vdiffnam(k), (/ 'lev' /), 'A', 'kg/kg/s', 'Vertical diffusion of '//cnst_name(k) )
522 : end do
523 :
524 1536 : if (.not. do_pbl_diags) then
525 0 : call addfld( 'PBLH' , horiz_only , 'A', 'm' , 'PBL height' )
526 0 : call addfld( 'QT' , (/ 'lev' /) , 'A', 'kg/kg' , 'Total water mixing ratio' )
527 0 : call addfld( 'SL' , (/ 'lev' /) , 'A', 'J/kg' , 'Liquid water static energy' )
528 0 : call addfld( 'SLV' , (/ 'lev' /) , 'A', 'J/kg' , 'Liq wat virtual static energy' )
529 0 : call addfld( 'SLFLX' , (/ 'ilev' /) , 'A', 'W/m2' , 'Liquid static energy flux' )
530 0 : call addfld( 'QTFLX' , (/ 'ilev' /) , 'A', 'W/m2' , 'Total water flux' )
531 0 : call addfld( 'TKE' , (/ 'ilev' /) , 'A', 'm2/s2' , 'Turbulent Kinetic Energy' )
532 0 : call addfld( 'TPERT' , horiz_only , 'A', 'K' , 'Perturbation temperature (eddies in PBL)' )
533 0 : call addfld( 'QPERT' , horiz_only , 'A', 'kg/kg' , 'Perturbation specific humidity (eddies in PBL)' )
534 :
535 0 : call addfld( 'UFLX' , (/ 'ilev' /) , 'A', 'W/m2' , 'Zonal momentum flux' )
536 0 : call addfld( 'VFLX' , (/ 'ilev' /) , 'A', 'W/m2' , 'Meridional momentm flux' )
537 0 : call register_vector_field('UFLX', 'VFLX')
538 : end if
539 :
540 1536 : call addfld( 'USTAR' , horiz_only , 'A', 'm/s' , 'Surface friction velocity' )
541 3072 : call addfld( 'KVH' , (/ 'ilev' /) , 'A', 'm2/s' , 'Vertical diffusion diffusivities (heat/moisture)' )
542 3072 : call addfld( 'KVM' , (/ 'ilev' /) , 'A', 'm2/s' , 'Vertical diffusion diffusivities (momentum)' )
543 3072 : call addfld( 'KVT' , (/ 'ilev' /) , 'A', 'm2/s' , 'Vertical diffusion kinematic molecular conductivity')
544 3072 : call addfld( 'CGS' , (/ 'ilev' /) , 'A', 's/m2' , 'Counter-gradient coeff on surface kinematic fluxes' )
545 3072 : call addfld( 'DTVKE' , (/ 'lev' /) , 'A', 'K/s' , 'dT/dt vertical diffusion KE dissipation' )
546 3072 : call addfld( 'DTV' , (/ 'lev' /) , 'A', 'K/s' , 'T vertical diffusion' )
547 3072 : call addfld( 'DUV' , (/ 'lev' /) , 'A', 'm/s2' , 'U vertical diffusion' )
548 3072 : call addfld( 'DVV' , (/ 'lev' /) , 'A', 'm/s2' , 'V vertical diffusion' )
549 :
550 : ! ---------------------------------------------------------------------------- !
551 : ! Below ( with '_PBL') are for detailed analysis of UW Moist Turbulence Scheme !
552 : ! ---------------------------------------------------------------------------- !
553 :
554 1536 : if (.not. do_pbl_diags) then
555 :
556 0 : call addfld( 'qt_pre_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qt_prePBL' )
557 0 : call addfld( 'sl_pre_PBL', (/ 'lev' /) , 'A', 'J/kg' , 'sl_prePBL' )
558 0 : call addfld( 'slv_pre_PBL', (/ 'lev' /) , 'A', 'J/kg' , 'slv_prePBL' )
559 0 : call addfld( 'u_pre_PBL', (/ 'lev' /) , 'A', 'm/s' , 'u_prePBL' )
560 0 : call addfld( 'v_pre_PBL', (/ 'lev' /) , 'A', 'm/s' , 'v_prePBL' )
561 0 : call addfld( 'qv_pre_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qv_prePBL' )
562 0 : call addfld( 'ql_pre_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'ql_prePBL' )
563 0 : call addfld( 'qi_pre_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qi_prePBL' )
564 0 : call addfld( 't_pre_PBL', (/ 'lev' /) , 'A', 'K' , 't_prePBL' )
565 0 : call addfld( 'rh_pre_PBL', (/ 'lev' /) , 'A', '%' , 'rh_prePBL' )
566 :
567 0 : call addfld( 'qt_aft_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qt_afterPBL' )
568 0 : call addfld( 'sl_aft_PBL', (/ 'lev' /) , 'A', 'J/kg' , 'sl_afterPBL' )
569 0 : call addfld( 'slv_aft_PBL', (/ 'lev' /) , 'A', 'J/kg' , 'slv_afterPBL' )
570 0 : call addfld( 'u_aft_PBL', (/ 'lev' /) , 'A', 'm/s' , 'u_afterPBL' )
571 0 : call addfld( 'v_aft_PBL', (/ 'lev' /) , 'A', 'm/s' , 'v_afterPBL' )
572 0 : call addfld( 'qv_aft_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qv_afterPBL' )
573 0 : call addfld( 'ql_aft_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'ql_afterPBL' )
574 0 : call addfld( 'qi_aft_PBL', (/ 'lev' /) , 'A', 'kg/kg' , 'qi_afterPBL' )
575 0 : call addfld( 't_aft_PBL', (/ 'lev' /) , 'A', 'K' , 't_afterPBL' )
576 0 : call addfld( 'rh_aft_PBL', (/ 'lev' /) , 'A', '%' , 'rh_afterPBL' )
577 :
578 0 : call addfld( 'slflx_PBL', (/ 'ilev' /) , 'A', 'J/m2/s' , 'sl flux by PBL' )
579 0 : call addfld( 'qtflx_PBL', (/ 'ilev' /) , 'A', 'kg/m2/s', 'qt flux by PBL' )
580 0 : call addfld( 'uflx_PBL', (/ 'ilev' /) , 'A', 'kg/m/s2', 'u flux by PBL' )
581 0 : call addfld( 'vflx_PBL', (/ 'ilev' /) , 'A', 'kg/m/s2', 'v flux by PBL' )
582 :
583 0 : call addfld( 'slflx_cg_PBL', (/ 'ilev' /) , 'A', 'J/m2/s' , 'sl_cg flux by PBL' )
584 0 : call addfld( 'qtflx_cg_PBL', (/ 'ilev' /) , 'A', 'kg/m2/s', 'qt_cg flux by PBL' )
585 0 : call addfld( 'uflx_cg_PBL', (/ 'ilev' /) , 'A', 'kg/m/s2', 'u_cg flux by PBL' )
586 0 : call addfld( 'vflx_cg_PBL', (/ 'ilev' /) , 'A', 'kg/m/s2', 'v_cg flux by PBL' )
587 :
588 0 : call addfld( 'qtten_PBL', (/ 'lev' /) , 'A', 'kg/kg/s', 'qt tendency by PBL' )
589 0 : call addfld( 'slten_PBL', (/ 'lev' /) , 'A', 'J/kg/s' , 'sl tendency by PBL' )
590 0 : call addfld( 'uten_PBL', (/ 'lev' /) , 'A', 'm/s2' , 'u tendency by PBL' )
591 0 : call addfld( 'vten_PBL', (/ 'lev' /) , 'A', 'm/s2' , 'v tendency by PBL' )
592 0 : call addfld( 'qvten_PBL', (/ 'lev' /) , 'A', 'kg/kg/s', 'qv tendency by PBL' )
593 0 : call addfld( 'qlten_PBL', (/ 'lev' /) , 'A', 'kg/kg/s', 'ql tendency by PBL' )
594 0 : call addfld( 'qiten_PBL', (/ 'lev' /) , 'A', 'kg/kg/s', 'qi tendency by PBL' )
595 0 : call addfld( 'tten_PBL', (/ 'lev' /) , 'A', 'K/s' , 'T tendency by PBL' )
596 0 : call addfld( 'rhten_PBL', (/ 'lev' /) , 'A', '%/s' , 'RH tendency by PBL' )
597 :
598 : end if
599 :
600 1536 : call addfld ('ustar',horiz_only, 'A', ' ',' ')
601 1536 : call addfld ('obklen',horiz_only, 'A', ' ',' ')
602 :
603 : ! ----------------------------
604 : ! determine default variables
605 : ! ----------------------------
606 :
607 : call phys_getopts( history_amwg_out = history_amwg, &
608 : history_eddy_out = history_eddy, &
609 : history_budget_out = history_budget, &
610 : history_budget_histfile_num_out = history_budget_histfile_num, &
611 1536 : history_waccm_out = history_waccm)
612 :
613 1536 : if (history_amwg) then
614 1536 : call add_default( vdiffnam(1), 1, ' ' )
615 1536 : call add_default( 'DTV' , 1, ' ' )
616 1536 : if (.not. do_pbl_diags) then
617 0 : call add_default( 'PBLH' , 1, ' ' )
618 : end if
619 : endif
620 :
621 1536 : if (history_eddy) then
622 0 : call add_default( 'UFLX ', 1, ' ' )
623 0 : call add_default( 'VFLX ', 1, ' ' )
624 : endif
625 :
626 1536 : if( history_budget ) then
627 0 : call add_default( vdiffnam(ixcldliq), history_budget_histfile_num, ' ' )
628 0 : call add_default( vdiffnam(ixcldice), history_budget_histfile_num, ' ' )
629 0 : if( history_budget_histfile_num > 1 ) then
630 0 : call add_default( vdiffnam(1), history_budget_histfile_num, ' ' )
631 0 : call add_default( 'DTV' , history_budget_histfile_num, ' ' )
632 : end if
633 : end if
634 :
635 1536 : if ( history_waccm ) then
636 0 : if (do_molec_diff) then
637 0 : call add_default ( 'TTPXMLC', 1, ' ' )
638 : end if
639 0 : call add_default( 'DUV' , 1, ' ' )
640 0 : call add_default( 'DVV' , 1, ' ' )
641 : end if
642 : ! ----------------------------
643 :
644 :
645 1536 : ksrftms_idx = pbuf_get_index('ksrftms')
646 1536 : tautmsx_idx = pbuf_get_index('tautmsx')
647 1536 : tautmsy_idx = pbuf_get_index('tautmsy')
648 :
649 1536 : dragblj_idx = pbuf_get_index('dragblj')
650 1536 : taubljx_idx = pbuf_get_index('taubljx')
651 1536 : taubljy_idx = pbuf_get_index('taubljy')
652 :
653 1536 : if (eddy_scheme == 'CLUBB_SGS') then
654 1536 : kvh_idx = pbuf_get_index('kvh')
655 : end if
656 :
657 1536 : if (do_hb_above_clubb) then
658 : ! pbuf field denoting top of clubb
659 1536 : clubbtop_idx = pbuf_get_index('clubbtop')
660 : end if
661 :
662 : ! Initialization of some pbuf fields
663 1536 : if (is_first_step()) then
664 : ! Initialization of pbuf fields tke, kvh, kvm are done in phys_inidat
665 768 : call pbuf_set_field(pbuf2d, tauresx_idx, 0.0_r8)
666 768 : call pbuf_set_field(pbuf2d, tauresy_idx, 0.0_r8)
667 768 : if (trim(shallow_scheme) == 'UNICON') then
668 0 : call pbuf_set_field(pbuf2d, qtl_flx_idx, 0.0_r8)
669 0 : call pbuf_set_field(pbuf2d, qti_flx_idx, 0.0_r8)
670 : end if
671 : end if
672 1536 : end subroutine vertical_diffusion_init
673 :
674 : ! =============================================================================== !
675 : ! !
676 : ! =============================================================================== !
677 :
678 741888 : subroutine vertical_diffusion_ts_init( pbuf2d, state )
679 :
680 : !-------------------------------------------------------------- !
681 : ! Timestep dependent setting, !
682 : ! At present only invokes upper bc code !
683 : !-------------------------------------------------------------- !
684 1536 : use upper_bc, only : ubc_timestep_init
685 : use physics_types , only : physics_state
686 : use ppgrid , only : begchunk, endchunk
687 :
688 : use physics_buffer, only : physics_buffer_desc
689 :
690 : type(physics_state), intent(in) :: state(begchunk:endchunk)
691 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
692 :
693 370944 : call ubc_timestep_init( pbuf2d, state)
694 :
695 370944 : end subroutine vertical_diffusion_ts_init
696 :
697 : ! =============================================================================== !
698 : ! !
699 : ! =============================================================================== !
700 :
701 62545392 : subroutine vertical_diffusion_tend( &
702 : ztodt , state , cam_in, &
703 : ustar , obklen , ptend , &
704 : cldn , pbuf)
705 : !---------------------------------------------------- !
706 : ! This is an interface routine for vertical diffusion !
707 : !---------------------------------------------------- !
708 370944 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_set_field
709 : use physics_types, only : physics_state, physics_ptend, physics_ptend_init
710 : use physics_types, only : set_dry_to_wet, set_wet_to_dry
711 :
712 : use camsrfexch, only : cam_in_t
713 : use cam_history, only : outfld
714 :
715 : use trb_mtn_stress_cam, only : trb_mtn_stress_tend
716 : use beljaars_drag_cam, only : beljaars_drag_tend
717 : use eddy_diff_cam, only : eddy_diff_tend
718 : use hb_diff, only : compute_hb_diff, compute_hb_free_atm_diff
719 : use wv_saturation, only : qsat
720 : use molec_diff, only : compute_molec_diff, vd_lu_qdecomp
721 : use constituents, only : qmincg, qmin, cnst_type
722 : use diffusion_solver, only : compute_vdiff, any, operator(.not.)
723 : use air_composition, only : cpairv, rairv !Needed for calculation of upward H flux
724 : use time_manager, only : get_nstep
725 : use constituents, only : cnst_get_type_byind, cnst_name, &
726 : cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx
727 : use physconst, only : pi
728 : use pbl_utils, only : virtem, calc_obklen, calc_ustar
729 : use upper_bc, only : ubc_get_vals, ubc_fixed_temp
730 : use upper_bc, only : ubc_get_flxs
731 : use coords_1d, only : Coords1D
732 : use phys_control, only : cam_physpkg_is
733 : use ref_pres, only : ptop_ref
734 :
735 : ! --------------- !
736 : ! Input Arguments !
737 : ! --------------- !
738 :
739 : type(physics_state), intent(inout) :: state ! Physics state variables
740 : type(cam_in_t), intent(in) :: cam_in ! Surface inputs
741 :
742 : real(r8), intent(in) :: ztodt ! 2 delta-t [ s ]
743 : real(r8), intent(in) :: cldn(pcols,pver) ! New stratus fraction [ fraction ]
744 :
745 : ! ---------------------- !
746 : ! Input-Output Arguments !
747 : ! ---------------------- !
748 :
749 : type(physics_ptend), intent(out) :: ptend ! Individual parameterization tendencies
750 : type(physics_buffer_desc), pointer :: pbuf(:)
751 :
752 : ! ---------------- !
753 : ! Output Arguments !
754 : ! ---------------- !
755 :
756 : real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ]
757 : real(r8), intent(out) :: obklen(pcols) ! Obukhov length [ m ]
758 :
759 : ! --------------- !
760 : ! Local Variables !
761 : ! --------------- !
762 :
763 : character(128) :: errstring ! Error status for compute_vdiff
764 :
765 : integer :: lchnk ! Chunk identifier
766 : integer :: ncol ! Number of atmospheric columns
767 : integer :: i, k, l, m ! column, level, constituent indices
768 :
769 : real(r8) :: dtk(pcols,pver) ! T tendency from KE dissipation
770 1489176 : real(r8), pointer :: tke(:,:) ! Turbulent kinetic energy [ m2/s2 ]
771 :
772 1489176 : real(r8), pointer :: qtl_flx(:,:) ! overbar(w'qtl') where qtl = qv + ql
773 1489176 : real(r8), pointer :: qti_flx(:,:) ! overbar(w'qti') where qti = qv + qi
774 :
775 : real(r8) :: cgs(pcols,pverp) ! Counter-gradient star [ cg/flux ]
776 : real(r8) :: cgh(pcols,pverp) ! Counter-gradient term for heat
777 : real(r8) :: rztodt ! 1./ztodt [ 1/s ]
778 1489176 : real(r8), pointer :: ksrftms(:) ! Turbulent mountain stress surface drag coefficient [ kg/s/m2 ]
779 1489176 : real(r8), pointer :: tautmsx(:) ! U component of turbulent mountain stress [ N/m2 ]
780 1489176 : real(r8), pointer :: tautmsy(:) ! V component of turbulent mountain stress [ N/m2 ]
781 : real(r8) :: tautotx(pcols) ! U component of total surface stress [ N/m2 ]
782 : real(r8) :: tautoty(pcols) ! V component of total surface stress [ N/m2 ]
783 :
784 1489176 : real(r8), pointer :: dragblj(:,:) ! Beljaars SGO form drag profile [ 1/s ]
785 1489176 : real(r8), pointer :: taubljx(:) ! U component of turbulent mountain stress [ N/m2 ]
786 1489176 : real(r8), pointer :: taubljy(:) ! V component of turbulent mountain stress [ N/m2 ]
787 :
788 1489176 : real(r8), pointer :: kvh_in(:,:) ! kvh from previous timestep [ m2/s ]
789 1489176 : real(r8), pointer :: kvm_in(:,:) ! kvm from previous timestep [ m2/s ]
790 1489176 : real(r8), pointer :: kvt(:,:) ! Molecular kinematic conductivity for temperature [ ]
791 : real(r8) :: kvq(pcols,pverp) ! Eddy diffusivity for constituents [ m2/s ]
792 : real(r8) :: kvh(pcols,pverp) ! Eddy diffusivity for heat [ m2/s ]
793 : real(r8) :: kvm(pcols,pverp) ! Eddy diffusivity for momentum [ m2/s ]
794 : real(r8) :: kvm_temp(pcols,pverp) ! Dummy eddy diffusivity for momentum (unused) [ m2/s ]
795 : real(r8) :: dtk_temp(pcols,pverp) ! Unused output from second compute_vdiff call
796 : real(r8) :: tautmsx_temp(pcols) ! Unused output from second compute_vdiff call
797 : real(r8) :: tautmsy_temp(pcols) ! Unused output from second compute_vdiff call
798 : real(r8) :: topflx_temp(pcols) ! Unused output from second compute_vdiff call
799 : real(r8) :: sprod(pcols,pverp) ! Shear production of tke [ m2/s3 ]
800 : real(r8) :: sfi(pcols,pverp) ! Saturation fraction at interfaces [ fraction ]
801 : real(r8) :: sl(pcols,pver)
802 : real(r8) :: qt(pcols,pver)
803 : real(r8) :: slv(pcols,pver)
804 : real(r8) :: sl_prePBL(pcols,pver)
805 : real(r8) :: qt_prePBL(pcols,pver)
806 : real(r8) :: slv_prePBL(pcols,pver)
807 : real(r8) :: slten(pcols,pver)
808 : real(r8) :: qtten(pcols,pver)
809 : real(r8) :: slflx(pcols,pverp)
810 : real(r8) :: qtflx(pcols,pverp)
811 : real(r8) :: uflx(pcols,pverp)
812 : real(r8) :: vflx(pcols,pverp)
813 : real(r8) :: slflx_cg(pcols,pverp)
814 : real(r8) :: qtflx_cg(pcols,pverp)
815 : real(r8) :: uflx_cg(pcols,pverp)
816 : real(r8) :: vflx_cg(pcols,pverp)
817 : real(r8) :: th(pcols,pver) ! Potential temperature
818 : real(r8) :: topflx(pcols) ! Molecular heat flux at top interface
819 : real(r8) :: rhoair
820 :
821 : real(r8) :: ri(pcols,pver) ! richardson number (HB output)
822 :
823 : ! for obklen calculation outside HB
824 : real(r8) :: thvs(pcols) ! Virtual potential temperature at surface
825 : real(r8) :: rrho(pcols) ! Reciprocal of density at surface
826 : real(r8) :: khfs(pcols) ! sfc kinematic heat flux [mK/s]
827 : real(r8) :: kqfs(pcols) ! sfc kinematic water vapor flux [m/s]
828 : real(r8) :: kbfs(pcols) ! sfc kinematic buoyancy flux [m^2/s^3]
829 :
830 : real(r8) :: ftem(pcols,pver) ! Saturation vapor pressure before PBL
831 : real(r8) :: ftem_prePBL(pcols,pver) ! Saturation vapor pressure before PBL
832 : real(r8) :: ftem_aftPBL(pcols,pver) ! Saturation vapor pressure after PBL
833 : real(r8) :: tem2(pcols,pver) ! Saturation specific humidity and RH
834 : real(r8) :: t_aftPBL(pcols,pver) ! Temperature after PBL diffusion
835 : real(r8) :: tten(pcols,pver) ! Temperature tendency by PBL diffusion
836 : real(r8) :: rhten(pcols,pver) ! RH tendency by PBL diffusion
837 : real(r8) :: qv_aft_PBL(pcols,pver) ! qv after PBL diffusion
838 : real(r8) :: ql_aft_PBL(pcols,pver) ! ql after PBL diffusion
839 : real(r8) :: qi_aft_PBL(pcols,pver) ! qi after PBL diffusion
840 : real(r8) :: s_aft_PBL(pcols,pver) ! s after PBL diffusion
841 : real(r8) :: u_aft_PBL(pcols,pver) ! u after PBL diffusion
842 : real(r8) :: v_aft_PBL(pcols,pver) ! v after PBL diffusion
843 : real(r8) :: qv_pro(pcols,pver)
844 : real(r8) :: ql_pro(pcols,pver)
845 : real(r8) :: qi_pro(pcols,pver)
846 : real(r8) :: s_pro(pcols,pver)
847 : real(r8) :: t_pro(pcols,pver)
848 1489176 : real(r8), pointer :: tauresx(:) ! Residual stress to be added in vdiff to correct
849 1489176 : real(r8), pointer :: tauresy(:) ! for turb stress mismatch between sfc and atm accumulated.
850 :
851 : ! Interpolated interface values.
852 : real(r8) :: tint(pcols,pver+1) ! Temperature [ K ]
853 : real(r8) :: rairi(pcols,pver+1) ! Gas constant [ J/K/kg ]
854 : real(r8) :: rhoi(pcols,pver+1) ! Density of air [ kg/m^3 ]
855 : real(r8) :: rhoi_dry(pcols,pver+1) ! Density of air based on dry air pressure [ kg/m^3 ]
856 :
857 : ! Upper boundary conditions
858 : real(r8) :: ubc_t(pcols) ! Temperature [ K ]
859 : real(r8) :: ubc_mmr(pcols,pcnst) ! Mixing ratios [ kg/kg ]
860 : real(r8) :: ubc_flux(pcols,pcnst) ! Constituent upper boundary flux (kg/s/m^2)
861 :
862 : ! Pressure coordinates used by the solver.
863 1489176 : type(Coords1D) :: p
864 1489176 : type(Coords1D) :: p_dry
865 :
866 1489176 : real(r8), pointer :: tpert(:)
867 1489176 : real(r8), pointer :: qpert(:)
868 1489176 : real(r8), pointer :: pblh(:)
869 :
870 : real(r8) :: tmp1(pcols) ! Temporary storage
871 :
872 : integer :: nstep
873 : real(r8) :: sum1, sum2, sum3, pdelx
874 : real(r8) :: sflx
875 :
876 : ! Copy state so we can pass to intent(inout) routines that return
877 : ! new state instead of a tendency.
878 : real(r8) :: s_tmp(pcols,pver)
879 : real(r8) :: u_tmp(pcols,pver)
880 : real(r8) :: v_tmp(pcols,pver)
881 : real(r8) :: q_tmp(pcols,pver,pcnst)
882 :
883 : ! kq_fac*sqrt(T)*m_d/rho for molecular diffusivity
884 : real(r8) :: kq_scal(pcols,pver+1)
885 : ! composition dependent mw_fac on interface level
886 : real(r8) :: mw_fac(pcols,pver+1,pcnst)
887 :
888 : ! Dry static energy top boundary condition.
889 : real(r8) :: dse_top(pcols)
890 :
891 : ! Copies of flux arrays used to zero out any parts that are applied
892 : ! elsewhere (e.g. by CLUBB).
893 : real(r8) :: taux(pcols)
894 : real(r8) :: tauy(pcols)
895 : real(r8) :: shflux(pcols)
896 : real(r8) :: cflux(pcols,pcnst)
897 1489176 : integer, pointer :: clubbtop(:) ! (pcols)
898 :
899 : logical :: lq(pcnst)
900 :
901 : ! ----------------------- !
902 : ! Main Computation Begins !
903 : ! ----------------------- !
904 :
905 : ! Assume 'wet' mixing ratios in diffusion code.
906 1489176 : call set_dry_to_wet(state, convert_cnst_type='dry')
907 :
908 1489176 : rztodt = 1._r8 / ztodt
909 1489176 : lchnk = state%lchnk
910 1489176 : ncol = state%ncol
911 :
912 1489176 : call pbuf_get_field(pbuf, tauresx_idx, tauresx)
913 1489176 : call pbuf_get_field(pbuf, tauresy_idx, tauresy)
914 1489176 : call pbuf_get_field(pbuf, tpert_idx, tpert)
915 1489176 : call pbuf_get_field(pbuf, qpert_idx, qpert)
916 1489176 : call pbuf_get_field(pbuf, pblh_idx, pblh)
917 :
918 : ! Interpolate temperature to interfaces.
919 138493368 : do k = 2, pver
920 2289140568 : do i = 1, ncol
921 2287651392 : tint(i,k) = 0.5_r8 * ( state%t(i,k) + state%t(i,k-1) )
922 : end do
923 : end do
924 24865776 : tint(:ncol,pver+1) = state%t(:ncol,pver)
925 :
926 : ! Get upper boundary values
927 1489176 : call ubc_get_vals( state%lchnk, ncol, state%pint, state%zi, ubc_t, ubc_mmr )
928 :
929 1489176 : if (waccmx_mode) then
930 0 : call ubc_get_flxs( state%lchnk, ncol, state%pint, state%zi, state%t, state%q, state%phis, ubc_flux )
931 : ! For WACCM-X, set ubc temperature to extrapolate from next two lower interface level temperatures
932 0 : tint(:ncol,1) = 1.5_r8*tint(:ncol,2)-.5_r8*tint(:ncol,3)
933 1489176 : else if(ubc_fixed_temp) then
934 0 : tint(:ncol,1) = ubc_t(:ncol)
935 : else
936 24865776 : tint(:ncol,1) = state%t(:ncol,1)
937 : end if
938 :
939 : ! Set up pressure coordinates for solver calls.
940 1489176 : p = Coords1D(state%pint(:ncol,:))
941 1489176 : p_dry = Coords1D(state%pintdry(:ncol,:))
942 :
943 : !------------------------------------------------------------------------
944 : ! Check to see if constituent dependent gas constant needed (WACCM-X)
945 : !------------------------------------------------------------------------
946 1489176 : if (waccmx_mode) then
947 0 : rairi(:ncol,1) = rairv(:ncol,1,lchnk)
948 0 : do k = 2, pver
949 0 : do i = 1, ncol
950 0 : rairi(i,k) = 0.5_r8 * (rairv(i,k,lchnk)+rairv(i,k-1,lchnk))
951 : end do
952 : end do
953 0 : rairi(:ncol,pver+1) = rairv(:ncol,pver,lchnk)
954 : else
955 2338872120 : rairi(:ncol,:pver+1) = rair
956 : endif
957 :
958 : ! Compute rho at interfaces.
959 141471720 : do k = 1, pver+1
960 2338872120 : do i = 1, ncol
961 2337382944 : rhoi(i,k) = p%ifc(i,k) / (rairi(i,k)*tint(i,k))
962 : end do
963 : end do
964 :
965 : ! Compute rho_dry at interfaces.
966 141471720 : do k = 1, pver+1
967 2338872120 : do i = 1, ncol
968 2337382944 : rhoi_dry(i,k) = p_dry%ifc(i,k) / (rairi(i,k)*tint(i,k))
969 : end do
970 : end do
971 :
972 : ! ---------------------------------------- !
973 : ! Computation of turbulent mountain stress !
974 : ! ---------------------------------------- !
975 :
976 : ! Consistent with the computation of 'normal' drag coefficient, we are using
977 : ! the raw input (u,v) to compute 'ksrftms', not the provisionally-marched 'u,v'
978 : ! within the iteration loop of the PBL scheme.
979 :
980 1489176 : call trb_mtn_stress_tend(state, pbuf, cam_in)
981 :
982 1489176 : call pbuf_get_field(pbuf, ksrftms_idx, ksrftms)
983 1489176 : call pbuf_get_field(pbuf, tautmsx_idx, tautmsx)
984 1489176 : call pbuf_get_field(pbuf, tautmsy_idx, tautmsy)
985 :
986 24865776 : tautotx(:ncol) = cam_in%wsx(:ncol) + tautmsx(:ncol)
987 24865776 : tautoty(:ncol) = cam_in%wsy(:ncol) + tautmsy(:ncol)
988 :
989 : ! ------------------------------------- !
990 : ! Computation of Beljaars SGO form drag !
991 : ! ------------------------------------- !
992 :
993 1489176 : call beljaars_drag_tend(state, pbuf, cam_in)
994 :
995 1489176 : call pbuf_get_field(pbuf, dragblj_idx, dragblj)
996 1489176 : call pbuf_get_field(pbuf, taubljx_idx, taubljx)
997 1489176 : call pbuf_get_field(pbuf, taubljy_idx, taubljy)
998 :
999 : ! Add Beljaars integrated drag
1000 :
1001 24865776 : tautotx(:ncol) = tautotx(:ncol) + taubljx(:ncol)
1002 24865776 : tautoty(:ncol) = tautoty(:ncol) + taubljy(:ncol)
1003 :
1004 : !----------------------------------------------------------------------- !
1005 : ! Computation of eddy diffusivities - Select appropriate PBL scheme !
1006 : !----------------------------------------------------------------------- !
1007 1489176 : call pbuf_get_field(pbuf, kvm_idx, kvm_in)
1008 1489176 : call pbuf_get_field(pbuf, kvh_idx, kvh_in)
1009 1489176 : call pbuf_get_field(pbuf, tke_idx, tke)
1010 :
1011 : ! Get potential temperature.
1012 2314006344 : th(:ncol,:pver) = state%t(:ncol,:pver) * state%exner(:ncol,:pver)
1013 :
1014 0 : select case (eddy_scheme)
1015 : case ( 'diag_TKE', 'SPCAM_m2005' )
1016 :
1017 : call eddy_diff_tend(state, pbuf, cam_in, &
1018 : ztodt, p, tint, rhoi, cldn, wstarent, &
1019 : kvm_in, kvh_in, ksrftms, dragblj, tauresx, tauresy, &
1020 : rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, &
1021 0 : tke, sprod, sfi)
1022 :
1023 : ! The diag_TKE scheme does not calculate the Monin-Obukhov length, which is used in dry deposition calculations.
1024 : ! Use the routines from pbl_utils to accomplish this. Assumes ustar and rrho have been set.
1025 0 : call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
1026 0 : call calc_obklen(ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
1027 : cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
1028 0 : khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol))
1029 :
1030 :
1031 : case ( 'HB', 'HBR', 'SPCAM_sam1mom' )
1032 :
1033 : ! Modification : We may need to use 'taux' instead of 'tautotx' here, for
1034 : ! consistency with the previous HB scheme.
1035 :
1036 :
1037 : call compute_hb_diff(ncol , &
1038 : th , state%t , state%q , state%zm , state%zi , &
1039 : state%pmid, state%u , state%v , tautotx , tautoty , &
1040 : cam_in%shf, cam_in%cflx(:,1), obklen , ustar , pblh , &
1041 : kvm , kvh , kvq , cgh , cgs , &
1042 : tpert , qpert , cldn , cam_in%ocnfrac , tke , &
1043 : ri , &
1044 0 : eddy_scheme)
1045 :
1046 0 : call outfld( 'HB_ri', ri, pcols, lchnk )
1047 :
1048 : case ( 'CLUBB_SGS' )
1049 : !
1050 : ! run HB scheme where CLUBB is not active when running cam7
1051 : !
1052 1489176 : if (do_hb_above_clubb) then
1053 : call compute_hb_free_atm_diff( ncol , &
1054 : th , state%t , state%q , state%zm , &
1055 : state%pmid, state%u , state%v , tautotx , tautoty , &
1056 : cam_in%shf, cam_in%cflx(:,1), obklen , ustar , &
1057 : kvm , kvh , kvq , cgh , cgs , &
1058 1489176 : ri )
1059 :
1060 1489176 : call pbuf_get_field(pbuf, clubbtop_idx, clubbtop)
1061 : !
1062 : ! zero out HB where CLUBB is active
1063 : !
1064 24865776 : do i=1,ncol
1065 841988340 : do k=clubbtop(i),pverp
1066 817122564 : kvm(i,k) = 0.0_r8
1067 817122564 : kvh(i,k) = 0.0_r8
1068 817122564 : kvq(i,k) = 0.0_r8
1069 817122564 : cgs(i,k) = 0.0_r8
1070 840499164 : cgh(i,k) = 0.0_r8
1071 : end do
1072 : end do
1073 :
1074 1489176 : call outfld( 'HB_ri', ri, pcols, lchnk )
1075 : else
1076 : ! CLUBB has only a bare-bones placeholder here. If using CLUBB, the
1077 : ! PBL diffusion will happen before coupling, so vertical_diffusion
1078 : ! is only handling other things, e.g. some boundary conditions, tms,
1079 : ! and molecular diffusion.
1080 :
1081 0 : call virtem(ncol, th(:ncol,pver),state%q(:ncol,pver,1), thvs(:ncol))
1082 :
1083 0 : call calc_ustar( ncol, state%t(:ncol,pver), state%pmid(:ncol,pver), &
1084 0 : cam_in%wsx(:ncol), cam_in%wsy(:ncol), rrho(:ncol), ustar(:ncol))
1085 : ! Use actual qflux, not lhf/latvap as was done previously
1086 0 : call calc_obklen( ncol, th(:ncol,pver), thvs(:ncol), cam_in%cflx(:ncol,1), &
1087 : cam_in%shf(:ncol), rrho(:ncol), ustar(:ncol), &
1088 0 : khfs(:ncol), kqfs(:ncol), kbfs(:ncol), obklen(:ncol))
1089 : ! These tendencies all applied elsewhere.
1090 0 : kvm = 0._r8
1091 0 : kvh = 0._r8
1092 0 : kvq = 0._r8
1093 : ! Not defined since PBL is not actually running here.
1094 0 : cgh = 0._r8
1095 0 : cgs = 0._r8
1096 : end if
1097 : end select
1098 :
1099 1489176 : call outfld( 'ustar', ustar(:), pcols, lchnk )
1100 1489176 : call outfld( 'obklen', obklen(:), pcols, lchnk )
1101 : !
1102 : ! add sponge layer vertical diffusion
1103 : !
1104 1489176 : if (allocated(kvm_sponge)) then
1105 7445880 : do k=1,size(kvm_sponge(:),1)
1106 100952280 : kvm(:ncol,1) = kvm(:ncol,1)+kvm_sponge(k)
1107 : end do
1108 : end if
1109 :
1110 : ! kvh (in pbuf) is used by other physics parameterizations, and as an initial guess in compute_eddy_diff
1111 : ! on the next timestep. It is not updated by the compute_vdiff call below.
1112 1489176 : call pbuf_set_field(pbuf, kvh_idx, kvh)
1113 :
1114 : ! kvm (in pbuf) is only used as an initial guess in compute_eddy_diff on the next timestep.
1115 : ! The contributions for molecular diffusion made to kvm by the call to compute_vdiff below
1116 : ! are not included in the pbuf as these are not needed in the initial guess by compute_eddy_diff.
1117 1489176 : call pbuf_set_field(pbuf, kvm_idx, kvm)
1118 :
1119 : !------------------------------------ !
1120 : ! Application of diffusivities !
1121 : !------------------------------------ !
1122 :
1123 : ! Set arrays from input state.
1124 94875749280 : q_tmp(:ncol,:,:) = state%q(:ncol,:,:)
1125 2314006344 : s_tmp(:ncol,:) = state%s(:ncol,:)
1126 2314006344 : u_tmp(:ncol,:) = state%u(:ncol,:)
1127 2314006344 : v_tmp(:ncol,:) = state%v(:ncol,:)
1128 :
1129 : !------------------------------------------------------ !
1130 : ! Write profile output before applying diffusion scheme !
1131 : !------------------------------------------------------ !
1132 :
1133 1489176 : if (.not. do_pbl_diags) then
1134 0 : sl_prePBL(:ncol,:pver) = s_tmp(:ncol,:) - latvap * q_tmp(:ncol,:,ixcldliq) &
1135 0 : - ( latvap + latice) * q_tmp(:ncol,:,ixcldice)
1136 : qt_prePBL(:ncol,:pver) = q_tmp(:ncol,:,1) + q_tmp(:ncol,:,ixcldliq) &
1137 0 : + q_tmp(:ncol,:,ixcldice)
1138 0 : slv_prePBL(:ncol,:pver) = sl_prePBL(:ncol,:pver) * ( 1._r8 + zvir*qt_prePBL(:ncol,:pver) )
1139 :
1140 0 : do k = 1, pver
1141 0 : call qsat(state%t(1:ncol,k), state%pmid(1:ncol,k), tem2(1:ncol,k), ftem(1:ncol,k), ncol)
1142 : end do
1143 0 : ftem_prePBL(:ncol,:) = state%q(:ncol,:,1)/ftem(:ncol,:)*100._r8
1144 :
1145 0 : call outfld( 'qt_pre_PBL ', qt_prePBL, pcols, lchnk )
1146 0 : call outfld( 'sl_pre_PBL ', sl_prePBL, pcols, lchnk )
1147 0 : call outfld( 'slv_pre_PBL ', slv_prePBL, pcols, lchnk )
1148 0 : call outfld( 'u_pre_PBL ', state%u, pcols, lchnk )
1149 0 : call outfld( 'v_pre_PBL ', state%v, pcols, lchnk )
1150 0 : call outfld( 'qv_pre_PBL ', state%q(:,:,1), pcols, lchnk )
1151 0 : call outfld( 'ql_pre_PBL ', state%q(:,:,ixcldliq), pcols, lchnk )
1152 0 : call outfld( 'qi_pre_PBL ', state%q(:,:,ixcldice), pcols, lchnk )
1153 0 : call outfld( 't_pre_PBL ', state%t, pcols, lchnk )
1154 0 : call outfld( 'rh_pre_PBL ', ftem_prePBL, pcols, lchnk )
1155 :
1156 : end if
1157 :
1158 : ! --------------------------------------------------------------------------------- !
1159 : ! Call the diffusivity solver and solve diffusion equation !
1160 : ! The final two arguments are optional function references to !
1161 : ! constituent-independent and constituent-dependent moleculuar diffusivity routines !
1162 : ! --------------------------------------------------------------------------------- !
1163 :
1164 : ! Modification : We may need to output 'tautotx_im,tautoty_im' from below 'compute_vdiff' and
1165 : ! separately print out as diagnostic output, because these are different from
1166 : ! the explicit 'tautotx, tautoty' computed above.
1167 : ! Note that the output 'tauresx,tauresy' from below subroutines are fully implicit ones.
1168 :
1169 1489176 : call pbuf_get_field(pbuf, kvt_idx, kvt)
1170 :
1171 1489176 : if (do_molec_diff .and. .not. waccmx_mode) then
1172 : ! Top boundary condition for dry static energy
1173 0 : dse_top(:ncol) = cpairv(:ncol,1,lchnk) * tint(:ncol,1) + &
1174 0 : gravit * state%zi(:ncol,1)
1175 : else
1176 24865776 : dse_top(:ncol) = 0._r8
1177 : end if
1178 :
1179 1489176 : select case (eddy_scheme)
1180 : case ('CLUBB_SGS')
1181 : ! CLUBB applies some fluxes itself, but we still want constituent
1182 : ! fluxes applied here (except water vapor).
1183 1489176 : taux = 0._r8
1184 1489176 : tauy = 0._r8
1185 1489176 : shflux = 0._r8
1186 25315992 : cflux(:,1) = 0._r8
1187 1489176 : if (cam_physpkg_is("cam7")) then
1188 : ! surface fluxes applied in clubb emissions module
1189 1014128856 : cflux(:,2:) = 0._r8
1190 : else
1191 0 : cflux(:,2:) = cam_in%cflx(:,2:)
1192 : end if
1193 : case default
1194 0 : taux = cam_in%wsx
1195 0 : tauy = cam_in%wsy
1196 0 : shflux = cam_in%shf
1197 1489176 : cflux = cam_in%cflx
1198 : end select
1199 :
1200 1489176 : if( any(fieldlist_wet) ) then
1201 :
1202 1489176 : if (do_molec_diff) then
1203 : call compute_molec_diff(state%lchnk, pcols, pver, pcnst, ncol, &
1204 : kvm, kvt, tint, rhoi, kq_scal, cnst_mw, &
1205 0 : mw_fac, nbot_molec)
1206 : end if
1207 :
1208 : call compute_vdiff( state%lchnk , &
1209 : pcols , pver , pcnst , ncol , tint , &
1210 : p , state%t , rhoi, ztodt , taux , &
1211 : tauy , shflux , cflux , &
1212 : kvh , kvm , kvq , cgs , cgh , &
1213 : state%zi , ksrftms , dragblj , &
1214 : qmincg , fieldlist_wet , fieldlist_molec,&
1215 : u_tmp , v_tmp , q_tmp , s_tmp , &
1216 : tautmsx , tautmsy , dtk , topflx , errstring , &
1217 : tauresx , tauresy , 1 , cpairv(:,:,state%lchnk), dse_top, &
1218 : do_molec_diff, waccmx_mode, &
1219 : vd_lu_qdecomp, &
1220 : ubc_mmr, ubc_flux, kvt, state%pmid, &
1221 : cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
1222 1489176 : kq_scal, mw_fac)
1223 :
1224 : call handle_errmsg(errstring, subname="compute_vdiff", &
1225 1489176 : extra_msg="Error in fieldlist_wet call from vertical_diffusion.")
1226 :
1227 : end if
1228 :
1229 1489176 : if( any( fieldlist_dry ) ) then
1230 :
1231 0 : if( do_molec_diff ) then
1232 : ! kvm is unused in the output here (since it was assigned
1233 : ! above), so we use a temp kvm for the inout argument, and
1234 : ! ignore the value output by compute_molec_diff.
1235 0 : kvm_temp = kvm
1236 : call compute_molec_diff(state%lchnk, pcols, pver, pcnst, ncol, &
1237 : kvm_temp, kvt, tint, rhoi_dry, kq_scal, cnst_mw, &
1238 0 : mw_fac, nbot_molec)
1239 : end if
1240 :
1241 : call compute_vdiff( state%lchnk , &
1242 : pcols , pver , pcnst , ncol , tint , &
1243 : p_dry , state%t , rhoi_dry, ztodt , taux , &
1244 : tauy , shflux , cflux , &
1245 : kvh , kvm , kvq , cgs , cgh , &
1246 : state%zi , ksrftms , dragblj , &
1247 : qmincg , fieldlist_dry , fieldlist_molec,&
1248 : u_tmp , v_tmp , q_tmp , s_tmp , &
1249 : tautmsx_temp , tautmsy_temp , dtk_temp , topflx_temp , errstring , &
1250 : tauresx , tauresy , 1 , cpairv(:,:,state%lchnk), dse_top, &
1251 : do_molec_diff , waccmx_mode, &
1252 : vd_lu_qdecomp, &
1253 : ubc_mmr, ubc_flux, kvt, state%pmiddry, &
1254 : cnst_mw, cnst_fixed_ubc, cnst_fixed_ubflx, nbot_molec, &
1255 0 : kq_scal, mw_fac)
1256 :
1257 : call handle_errmsg(errstring, subname="compute_vdiff", &
1258 0 : extra_msg="Error in fieldlist_dry call from vertical_diffusion.")
1259 :
1260 : end if
1261 :
1262 1489176 : if (prog_modal_aero) then
1263 :
1264 : ! Modal aerosol species not diffused, so just add the explicit surface fluxes to the
1265 : ! lowest layer. **NOTE** This code assumes wet mmr.
1266 :
1267 24865776 : tmp1(:ncol) = ztodt * gravit * state%rpdel(:ncol,pver)
1268 29783520 : do m = 1, pmam_ncnst
1269 28294344 : l = pmam_cnst_idx(m)
1270 473938920 : q_tmp(:ncol,pver,l) = q_tmp(:ncol,pver,l) + tmp1(:ncol) * cflux(:ncol,l)
1271 : enddo
1272 : end if
1273 :
1274 : ! -------------------------------------------------------- !
1275 : ! Diagnostics and output writing after applying PBL scheme !
1276 : ! -------------------------------------------------------- !
1277 :
1278 1489176 : if (.not. do_pbl_diags) then
1279 :
1280 0 : sl(:ncol,:pver) = s_tmp(:ncol,:) - latvap * q_tmp(:ncol,:,ixcldliq) &
1281 0 : - ( latvap + latice) * q_tmp(:ncol,:,ixcldice)
1282 : qt(:ncol,:pver) = q_tmp(:ncol,:,1) + q_tmp(:ncol,:,ixcldliq) &
1283 0 : + q_tmp(:ncol,:,ixcldice)
1284 0 : slv(:ncol,:pver) = sl(:ncol,:pver) * ( 1._r8 + zvir*qt(:ncol,:pver) )
1285 :
1286 0 : slflx(:ncol,1) = 0._r8
1287 0 : qtflx(:ncol,1) = 0._r8
1288 0 : uflx(:ncol,1) = 0._r8
1289 0 : vflx(:ncol,1) = 0._r8
1290 :
1291 0 : slflx_cg(:ncol,1) = 0._r8
1292 0 : qtflx_cg(:ncol,1) = 0._r8
1293 0 : uflx_cg(:ncol,1) = 0._r8
1294 0 : vflx_cg(:ncol,1) = 0._r8
1295 :
1296 0 : do k = 2, pver
1297 0 : do i = 1, ncol
1298 0 : rhoair = state%pint(i,k) / ( rair * ( ( 0.5_r8*(slv(i,k)+slv(i,k-1)) - gravit*state%zi(i,k))/cpair ) )
1299 0 : slflx(i,k) = kvh(i,k) * &
1300 0 : ( - rhoair*(sl(i,k-1)-sl(i,k))/(state%zm(i,k-1)-state%zm(i,k)) &
1301 0 : + cgh(i,k) )
1302 : qtflx(i,k) = kvh(i,k) * &
1303 0 : ( - rhoair*(qt(i,k-1)-qt(i,k))/(state%zm(i,k-1)-state%zm(i,k)) &
1304 0 : + rhoair*(cam_in%cflx(i,1)+cam_in%cflx(i,ixcldliq)+cam_in%cflx(i,ixcldice))*cgs(i,k) )
1305 : uflx(i,k) = kvm(i,k) * &
1306 0 : ( - rhoair*(u_tmp(i,k-1)-u_tmp(i,k))/(state%zm(i,k-1)-state%zm(i,k)))
1307 : vflx(i,k) = kvm(i,k) * &
1308 0 : ( - rhoair*(v_tmp(i,k-1)-v_tmp(i,k))/(state%zm(i,k-1)-state%zm(i,k)))
1309 0 : slflx_cg(i,k) = kvh(i,k) * cgh(i,k)
1310 : qtflx_cg(i,k) = kvh(i,k) * rhoair * ( cam_in%cflx(i,1) + &
1311 0 : cam_in%cflx(i,ixcldliq) + cam_in%cflx(i,ixcldice) ) * cgs(i,k)
1312 0 : uflx_cg(i,k) = 0._r8
1313 0 : vflx_cg(i,k) = 0._r8
1314 : end do
1315 : end do
1316 :
1317 : ! Modification : I should check whether slflx(:ncol,pverp) is correctly computed.
1318 : ! Note also that 'tautotx' is explicit total stress, different from
1319 : ! the ones that have been actually added into the atmosphere.
1320 :
1321 0 : slflx(:ncol,pverp) = cam_in%shf(:ncol)
1322 0 : qtflx(:ncol,pverp) = cam_in%cflx(:ncol,1)
1323 0 : uflx(:ncol,pverp) = tautotx(:ncol)
1324 0 : vflx(:ncol,pverp) = tautoty(:ncol)
1325 :
1326 0 : slflx_cg(:ncol,pverp) = 0._r8
1327 0 : qtflx_cg(:ncol,pverp) = 0._r8
1328 0 : uflx_cg(:ncol,pverp) = 0._r8
1329 0 : vflx_cg(:ncol,pverp) = 0._r8
1330 :
1331 0 : if (trim(shallow_scheme) == 'UNICON') then
1332 0 : call pbuf_get_field(pbuf, qtl_flx_idx, qtl_flx)
1333 0 : call pbuf_get_field(pbuf, qti_flx_idx, qti_flx)
1334 0 : qtl_flx(:ncol,1) = 0._r8
1335 0 : qti_flx(:ncol,1) = 0._r8
1336 0 : do k = 2, pver
1337 0 : do i = 1, ncol
1338 : ! For use in the cloud macrophysics
1339 : ! Note that density is not added here. Also, only consider local transport term.
1340 0 : qtl_flx(i,k) = - kvh(i,k)*(q_tmp(i,k-1,1)-q_tmp(i,k,1)+q_tmp(i,k-1,ixcldliq)-q_tmp(i,k,ixcldliq))/&
1341 0 : (state%zm(i,k-1)-state%zm(i,k))
1342 0 : qti_flx(i,k) = - kvh(i,k)*(q_tmp(i,k-1,1)-q_tmp(i,k,1)+q_tmp(i,k-1,ixcldice)-q_tmp(i,k,ixcldice))/&
1343 0 : (state%zm(i,k-1)-state%zm(i,k))
1344 : end do
1345 : end do
1346 0 : do i = 1, ncol
1347 0 : rhoair = state%pint(i,pverp)/(rair*((slv(i,pver)-gravit*state%zi(i,pverp))/cpair))
1348 0 : qtl_flx(i,pverp) = cam_in%cflx(i,1)/rhoair
1349 0 : qti_flx(i,pverp) = cam_in%cflx(i,1)/rhoair
1350 : end do
1351 : end if
1352 :
1353 : end if
1354 :
1355 : ! --------------------------------------------------------------- !
1356 : ! Convert the new profiles into vertical diffusion tendencies. !
1357 : ! Convert KE dissipative heat change into "temperature" tendency. !
1358 : ! --------------------------------------------------------------- !
1359 : ! All variables are modified by vertical diffusion
1360 :
1361 62545392 : lq(:) = .TRUE.
1362 : call physics_ptend_init(ptend,state%psetcols, "vertical diffusion", &
1363 1489176 : ls=.true., lu=.true., lv=.true., lq=lq)
1364 :
1365 2314006344 : ptend%s(:ncol,:) = ( s_tmp(:ncol,:) - state%s(:ncol,:) ) * rztodt
1366 2314006344 : ptend%u(:ncol,:) = ( u_tmp(:ncol,:) - state%u(:ncol,:) ) * rztodt
1367 2314006344 : ptend%v(:ncol,:) = ( v_tmp(:ncol,:) - state%v(:ncol,:) ) * rztodt
1368 94875749280 : ptend%q(:ncol,:pver,:) = ( q_tmp(:ncol,:pver,:) - state%q(:ncol,:pver,:) ) * rztodt
1369 :
1370 : ! Convert tendencies of dry constituents to dry basis.
1371 62545392 : do m = 1,pcnst
1372 62545392 : if (cnst_type(m).eq.'dry') then
1373 69420190320 : ptend%q(:ncol,:pver,m) = ptend%q(:ncol,:pver,m)*state%pdel(:ncol,:pver)/state%pdeldry(:ncol,:pver)
1374 : endif
1375 : end do
1376 : ! convert wet mmr back to dry before conservation check
1377 1489176 : call set_wet_to_dry(state, convert_cnst_type='dry')
1378 :
1379 1489176 : if (.not. do_pbl_diags) then
1380 0 : slten(:ncol,:) = ( sl(:ncol,:) - sl_prePBL(:ncol,:) ) * rztodt
1381 0 : qtten(:ncol,:) = ( qt(:ncol,:) - qt_prePBL(:ncol,:) ) * rztodt
1382 : end if
1383 :
1384 : ! ------------------------------------------------------------ !
1385 : ! In order to perform 'pseudo-conservative variable diffusion' !
1386 : ! perform the following two stages: !
1387 : ! !
1388 : ! I. Re-set (1) 'qvten' by 'qtten', and 'qlten = qiten = 0' !
1389 : ! (2) 'sten' by 'slten', and !
1390 : ! (3) 'qlten = qiten = 0' !
1391 : ! !
1392 : ! II. Apply 'positive_moisture' !
1393 : ! !
1394 : ! ------------------------------------------------------------ !
1395 :
1396 1489176 : if( (eddy_scheme .eq. 'diag_TKE' .or. eddy_scheme .eq. 'SPCAM_m2005') .and. do_pseudocon_diff ) then
1397 :
1398 0 : ptend%q(:ncol,:pver,1) = qtten(:ncol,:pver)
1399 0 : ptend%s(:ncol,:pver) = slten(:ncol,:pver)
1400 0 : ptend%q(:ncol,:pver,ixcldliq) = 0._r8
1401 0 : ptend%q(:ncol,:pver,ixcldice) = 0._r8
1402 0 : if (ixnumliq > 0) ptend%q(:ncol,:pver,ixnumliq) = 0._r8
1403 0 : if (ixnumice > 0) ptend%q(:ncol,:pver,ixnumice) = 0._r8
1404 :
1405 0 : do i = 1, ncol
1406 0 : do k = 1, pver
1407 0 : qv_pro(i,k) = state%q(i,k,1) + ptend%q(i,k,1) * ztodt
1408 0 : ql_pro(i,k) = state%q(i,k,ixcldliq) + ptend%q(i,k,ixcldliq) * ztodt
1409 0 : qi_pro(i,k) = state%q(i,k,ixcldice) + ptend%q(i,k,ixcldice) * ztodt
1410 0 : s_pro(i,k) = state%s(i,k) + ptend%s(i,k) * ztodt
1411 0 : t_pro(i,k) = state%t(i,k) + (1._r8/cpair)*ptend%s(i,k) * ztodt
1412 : end do
1413 : end do
1414 :
1415 0 : call positive_moisture( cpair, latvap, latvap+latice, ncol, pver, ztodt, qmin(1), qmin(ixcldliq), qmin(ixcldice), &
1416 0 : state%pdel(:ncol,pver:1:-1), qv_pro(:ncol,pver:1:-1), ql_pro(:ncol,pver:1:-1), &
1417 0 : qi_pro(:ncol,pver:1:-1), t_pro(:ncol,pver:1:-1), s_pro(:ncol,pver:1:-1), &
1418 0 : ptend%q(:ncol,pver:1:-1,1), ptend%q(:ncol,pver:1:-1,ixcldliq), &
1419 0 : ptend%q(:ncol,pver:1:-1,ixcldice), ptend%s(:ncol,pver:1:-1) )
1420 :
1421 : end if
1422 :
1423 : ! ----------------------------------------------------------------- !
1424 : ! Re-calculate diagnostic output variables after vertical diffusion !
1425 : ! ----------------------------------------------------------------- !
1426 :
1427 1489176 : if (.not. do_pbl_diags) then
1428 :
1429 0 : qv_aft_PBL(:ncol,:pver) = state%q(:ncol,:pver,1) + ptend%q(:ncol,:pver,1) * ztodt
1430 0 : ql_aft_PBL(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq) + ptend%q(:ncol,:pver,ixcldliq) * ztodt
1431 0 : qi_aft_PBL(:ncol,:pver) = state%q(:ncol,:pver,ixcldice) + ptend%q(:ncol,:pver,ixcldice) * ztodt
1432 0 : s_aft_PBL(:ncol,:pver) = state%s(:ncol,:pver) + ptend%s(:ncol,:pver) * ztodt
1433 0 : t_aftPBL(:ncol,:pver) = ( s_aft_PBL(:ncol,:pver) - gravit*state%zm(:ncol,:pver) ) / cpair
1434 :
1435 0 : u_aft_PBL(:ncol,:pver) = state%u(:ncol,:pver) + ptend%u(:ncol,:pver) * ztodt
1436 0 : v_aft_PBL(:ncol,:pver) = state%v(:ncol,:pver) + ptend%v(:ncol,:pver) * ztodt
1437 :
1438 0 : do k = 1, pver
1439 0 : call qsat(t_aftPBL(1:ncol,k), state%pmid(1:ncol,k), tem2(1:ncol,k), ftem(1:ncol,k), ncol)
1440 : end do
1441 0 : ftem_aftPBL(:ncol,:pver) = qv_aft_PBL(:ncol,:pver) / ftem(:ncol,:pver) * 100._r8
1442 :
1443 0 : tten(:ncol,:pver) = ( t_aftPBL(:ncol,:pver) - state%t(:ncol,:pver) ) * rztodt
1444 0 : rhten(:ncol,:pver) = ( ftem_aftPBL(:ncol,:pver) - ftem_prePBL(:ncol,:pver) ) * rztodt
1445 :
1446 : end if
1447 :
1448 :
1449 : ! -------------------------------------------------------------- !
1450 : ! mass conservation check.........
1451 : ! -------------------------------------------------------------- !
1452 1489176 : if (diff_cnsrv_mass_check) then
1453 :
1454 : ! Conservation check
1455 0 : do m = 1, pcnst
1456 0 : fixed_ubc: if ((.not.cnst_fixed_ubc(m)).and.(.not.cnst_fixed_ubflx(m))) then
1457 0 : col_loop: do i = 1, ncol
1458 0 : sum1 = 0._r8
1459 0 : sum2 = 0._r8
1460 0 : sum3 = 0._r8
1461 0 : do k = 1, pver
1462 0 : if(cnst_get_type_byind(m).eq.'wet') then
1463 0 : pdelx = state%pdel(i,k)
1464 : else
1465 0 : pdelx = state%pdeldry(i,k)
1466 : endif
1467 0 : sum1 = sum1 + state%q(i,k,m)*pdelx/gravit ! total column
1468 0 : sum2 = sum2 +(state%q(i,k,m)+ptend%q(i,k,m)*ztodt)*pdelx/ gravit ! total column after tendancy is applied
1469 0 : sum3 = sum3 +( ptend%q(i,k,m)*ztodt)*pdelx/ gravit ! rate of change in column
1470 : enddo
1471 0 : sum1 = sum1 + (cam_in%cflx(i,m) * ztodt) ! add in surface flux (kg/m2)
1472 0 : sflx = (cam_in%cflx(i,m) * ztodt)
1473 0 : if (sum1>1.e-36_r8) then
1474 0 : if( abs((sum2-sum1)/sum1) .gt. 1.e-12_r8 ) then
1475 0 : nstep = get_nstep()
1476 : write(iulog,'(a,a8,a,I4,2f8.3,5e25.16)') &
1477 0 : 'MASSCHECK vert diff : nstep,lon,lat,mass1,mass2,sum3,sflx,rel-diff : ', &
1478 0 : trim(cnst_name(m)), ' : ', nstep, state%lon(i)*180._r8/pi, state%lat(i)*180._r8/pi, &
1479 0 : sum1, sum2, sum3, sflx, abs(sum2-sum1)/sum1
1480 : !xxx call endrun('vertical_diffusion_tend : mass not conserved' )
1481 : endif
1482 : endif
1483 : enddo col_loop
1484 : endif fixed_ubc
1485 : enddo
1486 : endif
1487 :
1488 : ! -------------------------------------------------------------- !
1489 : ! Writing state variables after PBL scheme for detailed analysis !
1490 : ! -------------------------------------------------------------- !
1491 :
1492 1489176 : if (.not. do_pbl_diags) then
1493 :
1494 0 : call outfld( 'sl_aft_PBL' , sl, pcols, lchnk )
1495 0 : call outfld( 'qt_aft_PBL' , qt, pcols, lchnk )
1496 0 : call outfld( 'slv_aft_PBL' , slv, pcols, lchnk )
1497 0 : call outfld( 'u_aft_PBL' , u_aft_PBL, pcols, lchnk )
1498 0 : call outfld( 'v_aft_PBL' , v_aft_PBL, pcols, lchnk )
1499 0 : call outfld( 'qv_aft_PBL' , qv_aft_PBL, pcols, lchnk )
1500 0 : call outfld( 'ql_aft_PBL' , ql_aft_PBL, pcols, lchnk )
1501 0 : call outfld( 'qi_aft_PBL' , qi_aft_PBL, pcols, lchnk )
1502 0 : call outfld( 't_aft_PBL ' , t_aftPBL, pcols, lchnk )
1503 0 : call outfld( 'rh_aft_PBL' , ftem_aftPBL, pcols, lchnk )
1504 0 : call outfld( 'slflx_PBL' , slflx, pcols, lchnk )
1505 0 : call outfld( 'qtflx_PBL' , qtflx, pcols, lchnk )
1506 0 : call outfld( 'uflx_PBL' , uflx, pcols, lchnk )
1507 0 : call outfld( 'vflx_PBL' , vflx, pcols, lchnk )
1508 0 : call outfld( 'slflx_cg_PBL' , slflx_cg, pcols, lchnk )
1509 0 : call outfld( 'qtflx_cg_PBL' , qtflx_cg, pcols, lchnk )
1510 0 : call outfld( 'uflx_cg_PBL' , uflx_cg, pcols, lchnk )
1511 0 : call outfld( 'vflx_cg_PBL' , vflx_cg, pcols, lchnk )
1512 0 : call outfld( 'slten_PBL' , slten, pcols, lchnk )
1513 0 : call outfld( 'qtten_PBL' , qtten, pcols, lchnk )
1514 0 : call outfld( 'uten_PBL' , ptend%u, pcols, lchnk )
1515 0 : call outfld( 'vten_PBL' , ptend%v, pcols, lchnk )
1516 0 : call outfld( 'qvten_PBL' , ptend%q(:,:,1), pcols, lchnk )
1517 0 : call outfld( 'qlten_PBL' , ptend%q(:,:,ixcldliq), pcols, lchnk )
1518 0 : call outfld( 'qiten_PBL' , ptend%q(:,:,ixcldice), pcols, lchnk )
1519 0 : call outfld( 'tten_PBL' , tten, pcols, lchnk )
1520 0 : call outfld( 'rhten_PBL' , rhten, pcols, lchnk )
1521 :
1522 : end if
1523 :
1524 : ! ------------------------------------------- !
1525 : ! Writing the other standard output variables !
1526 : ! ------------------------------------------- !
1527 :
1528 1489176 : if (.not. do_pbl_diags) then
1529 0 : call outfld( 'QT' , qt, pcols, lchnk )
1530 0 : call outfld( 'SL' , sl, pcols, lchnk )
1531 0 : call outfld( 'SLV' , slv, pcols, lchnk )
1532 0 : call outfld( 'SLFLX' , slflx, pcols, lchnk )
1533 0 : call outfld( 'QTFLX' , qtflx, pcols, lchnk )
1534 0 : call outfld( 'UFLX' , uflx, pcols, lchnk )
1535 0 : call outfld( 'VFLX' , vflx, pcols, lchnk )
1536 0 : call outfld( 'TKE' , tke, pcols, lchnk )
1537 :
1538 0 : call outfld( 'PBLH' , pblh, pcols, lchnk )
1539 0 : call outfld( 'TPERT' , tpert, pcols, lchnk )
1540 0 : call outfld( 'QPERT' , qpert, pcols, lchnk )
1541 : end if
1542 1489176 : call outfld( 'USTAR' , ustar, pcols, lchnk )
1543 1489176 : call outfld( 'KVH' , kvh, pcols, lchnk )
1544 1489176 : call outfld( 'KVT' , kvt, pcols, lchnk )
1545 1489176 : call outfld( 'KVM' , kvm, pcols, lchnk )
1546 1489176 : call outfld( 'CGS' , cgs, pcols, lchnk )
1547 2314006344 : dtk(:ncol,:) = dtk(:ncol,:) / cpair / ztodt ! Normalize heating for history
1548 1489176 : call outfld( 'DTVKE' , dtk, pcols, lchnk )
1549 2314006344 : dtk(:ncol,:) = ptend%s(:ncol,:) / cpair ! Normalize heating for history using dtk
1550 1489176 : call outfld( 'DTV' , dtk, pcols, lchnk )
1551 1489176 : call outfld( 'DUV' , ptend%u, pcols, lchnk )
1552 1489176 : call outfld( 'DVV' , ptend%v, pcols, lchnk )
1553 62545392 : do m = 1, pcnst
1554 62545392 : call outfld( vdiffnam(m) , ptend%q(1,1,m), pcols, lchnk )
1555 : end do
1556 1489176 : if( do_molec_diff ) then
1557 0 : call outfld( 'TTPXMLC' , topflx, pcols, lchnk )
1558 : end if
1559 :
1560 1489176 : call p%finalize()
1561 1489176 : call p_dry%finalize()
1562 :
1563 2978352 : end subroutine vertical_diffusion_tend
1564 :
1565 : ! =============================================================================== !
1566 : ! !
1567 : ! =============================================================================== !
1568 :
1569 0 : subroutine positive_moisture( cp, xlv, xls, ncol, mkx, dt, qvmin, qlmin, qimin, &
1570 0 : dp, qv, ql, qi, t, s, qvten, qlten, qiten, sten )
1571 : ! ------------------------------------------------------------------------------- !
1572 : ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, !
1573 : ! force them to be larger than minimum value by (1) condensating water vapor !
1574 : ! into liquid or ice, and (2) by transporting water vapor from the very lower !
1575 : ! layer. '2._r8' is multiplied to the minimum values for safety. !
1576 : ! Update final state variables and tendencies associated with this correction. !
1577 : ! If any condensation happens, update (s,t) too. !
1578 : ! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
1579 : ! input tendencies. !
1580 : ! Be careful the order of k : '1': near-surface layer, 'mkx' : top layer !
1581 : ! ------------------------------------------------------------------------------- !
1582 : implicit none
1583 : integer, intent(in) :: ncol, mkx
1584 : real(r8), intent(in) :: cp, xlv, xls
1585 : real(r8), intent(in) :: dt, qvmin, qlmin, qimin
1586 : real(r8), intent(in) :: dp(ncol,mkx)
1587 : real(r8), intent(inout) :: qv(ncol,mkx), ql(ncol,mkx), qi(ncol,mkx), t(ncol,mkx), s(ncol,mkx)
1588 : real(r8), intent(inout) :: qvten(ncol,mkx), qlten(ncol,mkx), qiten(ncol,mkx), sten(ncol,mkx)
1589 : integer i, k
1590 : real(r8) dql, dqi, dqv, sum, aa, dum
1591 :
1592 : ! Modification : I should check whether this is exactly same as the one used in
1593 : ! shallow convection and cloud macrophysics.
1594 :
1595 0 : do i = 1, ncol
1596 0 : do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface
1597 0 : dql = max(0._r8,1._r8*qlmin-ql(i,k))
1598 0 : dqi = max(0._r8,1._r8*qimin-qi(i,k))
1599 0 : qlten(i,k) = qlten(i,k) + dql/dt
1600 0 : qiten(i,k) = qiten(i,k) + dqi/dt
1601 0 : qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
1602 0 : sten(i,k) = sten(i,k) + xlv * (dql/dt) + xls * (dqi/dt)
1603 0 : ql(i,k) = ql(i,k) + dql
1604 0 : qi(i,k) = qi(i,k) + dqi
1605 0 : qv(i,k) = qv(i,k) - dql - dqi
1606 0 : s(i,k) = s(i,k) + xlv * dql + xls * dqi
1607 0 : t(i,k) = t(i,k) + (xlv * dql + xls * dqi)/cp
1608 0 : dqv = max(0._r8,1._r8*qvmin-qv(i,k))
1609 0 : qvten(i,k) = qvten(i,k) + dqv/dt
1610 0 : qv(i,k) = qv(i,k) + dqv
1611 0 : if( k .ne. 1 ) then
1612 0 : qv(i,k-1) = qv(i,k-1) - dqv*dp(i,k)/dp(i,k-1)
1613 0 : qvten(i,k-1) = qvten(i,k-1) - dqv*dp(i,k)/dp(i,k-1)/dt
1614 : endif
1615 0 : qv(i,k) = max(qv(i,k),qvmin)
1616 0 : ql(i,k) = max(ql(i,k),qlmin)
1617 0 : qi(i,k) = max(qi(i,k),qimin)
1618 : end do
1619 : ! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally
1620 : ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
1621 : ! preserves column moisture.
1622 0 : if( dqv .gt. 1.e-20_r8 ) then
1623 : sum = 0._r8
1624 0 : do k = 1, mkx
1625 0 : if( qv(i,k) .gt. 2._r8*qvmin ) sum = sum + qv(i,k)*dp(i,k)
1626 : enddo
1627 0 : aa = dqv*dp(i,1)/max(1.e-20_r8,sum)
1628 0 : if( aa .lt. 0.5_r8 ) then
1629 0 : do k = 1, mkx
1630 0 : if( qv(i,k) .gt. 2._r8*qvmin ) then
1631 0 : dum = aa*qv(i,k)
1632 0 : qv(i,k) = qv(i,k) - dum
1633 0 : qvten(i,k) = qvten(i,k) - dum/dt
1634 : endif
1635 : enddo
1636 : else
1637 0 : write(iulog,*) 'Full positive_moisture is impossible in vertical_diffusion'
1638 : endif
1639 : endif
1640 : end do
1641 0 : return
1642 :
1643 1489176 : end subroutine positive_moisture
1644 :
1645 : end module vertical_diffusion
|
