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