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