Line data Source code
1 : module eddy_diff_cam
2 :
3 : use shr_kind_mod, only: i4 => shr_kind_i4, r8 => shr_kind_r8
4 : use ppgrid, only: pcols, pver, pverp
5 : use cam_logfile, only: iulog
6 : use cam_abortutils, only: endrun
7 : use physconst, only: gravit, cpair, rair, zvir, latvap, latice, karman
8 : use diffusion_solver, only: vdiff_selector
9 : use eddy_diff, only: ncvmax
10 : use time_manager, only: is_first_step
11 : use physics_buffer, only: physics_buffer_desc
12 : use spmd_utils, only: masterproc
13 : use phys_control, only: phys_getopts
14 :
15 : implicit none
16 : private
17 :
18 : public :: eddy_diff_readnl
19 : public :: eddy_diff_register
20 : public :: eddy_diff_init
21 : public :: eddy_diff_tend
22 :
23 : ! Is UNICON switched on (and thus interacting with eddy_diff via pbuf)?
24 : logical :: unicon_is_on
25 :
26 : ! Number of iterations for solution
27 : integer, parameter :: nturb = 5
28 :
29 : ! Logical switches for moist mixing ratio diffusion
30 : type(vdiff_selector) :: fieldlist_wet
31 : ! Logical switches for molecular diffusion
32 : ! (Molecular diffusion is not done here.)
33 : type(vdiff_selector) :: fieldlist_molec
34 :
35 : integer :: ntop_eddy, nbot_eddy
36 :
37 : ! Cloud mass constituent indices
38 : integer :: ixcldliq, ixcldice
39 :
40 : ! input pbuf field indices
41 : integer :: qrl_idx = -1
42 : integer :: wsedl_idx = -1
43 :
44 : ! output pbuf field indices for UNICON
45 : integer :: bprod_idx = -1
46 : integer :: ipbl_idx = -1
47 : integer :: kpblh_idx = -1
48 : integer :: wstarPBL_idx = -1
49 : integer :: tkes_idx = -1
50 : integer :: went_idx = -1
51 :
52 : ! Mixing lengths squared.
53 : ! Used for computing free air diffusivity.
54 : real(r8) :: ml2(pver+1)
55 :
56 : ! Various namelist options to limit or tweak the effects of eddy diffusion.
57 :
58 : ! Pressure defining the bottom of the upper atmosphere for kvh scaling (Pa)
59 : real(r8) :: kv_top_pressure = 0._r8
60 : ! Eddy diffusivity scale factor for upper atmosphere
61 : real(r8) :: kv_top_scale = 1._r8
62 : ! Eddy diffusivity scale factor for the free troposphere
63 : real(r8) :: kv_freetrop_scale = 1._r8
64 :
65 : ! The following all have to be set in all cases.
66 : real(r8), parameter :: unset_r8 = huge(1._r8)
67 : ! Maximum master length for diag_TKE
68 : real(r8) :: eddy_lbulk_max = unset_r8
69 : ! Maximum dissipation length for diag_TKE
70 : real(r8) :: eddy_leng_max = unset_r8
71 : ! Bottom pressure level (hPa) for eddy_leng_max
72 : real(r8) :: eddy_max_bot_pressure = unset_r8
73 : ! Moist entrainment enhancement param
74 : real(r8) :: eddy_moist_entrain_a2l = unset_r8
75 :
76 : contains
77 :
78 0 : subroutine eddy_diff_readnl(nlfile)
79 : use namelist_utils, only: find_group_name
80 : use units, only: getunit, freeunit
81 : use spmd_utils, only: masterprocid, mpi_real8, mpicom
82 : use shr_log_mod, only: errMsg => shr_log_errMsg
83 :
84 : ! filepath for file containing namelist input
85 : character(len=*), intent(in) :: nlfile
86 :
87 : ! file unit and error code
88 : integer :: unitn, ierr
89 :
90 : character(len=*), parameter :: subname = 'eddy_diff_readnl'
91 :
92 : namelist /eddy_diff_nl/ kv_top_pressure, kv_top_scale, &
93 : kv_freetrop_scale, eddy_lbulk_max, eddy_leng_max, &
94 : eddy_max_bot_pressure, eddy_moist_entrain_a2l
95 :
96 0 : if (masterproc) then
97 0 : unitn = getunit()
98 0 : open( unitn, file=trim(nlfile), status='old' )
99 0 : call find_group_name(unitn, 'eddy_diff_nl', status=ierr)
100 0 : if (ierr == 0) then
101 0 : read(unitn, eddy_diff_nl, iostat=ierr)
102 : end if
103 0 : if (ierr /= 0) then
104 0 : call endrun(subname // ':: ERROR reading namelist')
105 : end if
106 0 : close(unitn)
107 0 : call freeunit(unitn)
108 : end if
109 :
110 0 : call mpi_bcast(kv_top_pressure, 1, mpi_real8, masterprocid, mpicom, ierr)
111 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
112 0 : call mpi_bcast(kv_top_scale, 1, mpi_real8, masterprocid, mpicom, ierr)
113 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
114 0 : call mpi_bcast(kv_freetrop_scale, 1, mpi_real8, masterprocid, mpicom, ierr)
115 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
116 :
117 0 : call mpi_bcast(eddy_lbulk_max, 1, mpi_real8, masterprocid, mpicom, ierr)
118 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
119 0 : call mpi_bcast(eddy_leng_max, 1, mpi_real8, masterprocid, mpicom, ierr)
120 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
121 0 : call mpi_bcast(eddy_max_bot_pressure, 1, mpi_real8, masterprocid, mpicom, ierr)
122 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
123 0 : call mpi_bcast(eddy_moist_entrain_a2l, 1, mpi_real8, masterprocid, mpicom, ierr)
124 0 : if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
125 :
126 0 : end subroutine eddy_diff_readnl
127 :
128 0 : subroutine eddy_diff_register()
129 : use physics_buffer, only: pbuf_add_field, dtype_r8, dtype_i4
130 :
131 : character(len=16) :: shallow_scheme
132 :
133 : ! Check for UNICON and add relevant pbuf entries.
134 0 : call phys_getopts(shallow_scheme_out=shallow_scheme)
135 :
136 0 : unicon_is_on = (shallow_scheme == "UNICON")
137 :
138 0 : if (unicon_is_on) then
139 0 : call pbuf_add_field('bprod', 'global', dtype_r8, (/pcols,pverp/), bprod_idx)
140 0 : call pbuf_add_field('ipbl', 'global', dtype_i4, (/pcols/), ipbl_idx)
141 0 : call pbuf_add_field('kpblh', 'global', dtype_i4, (/pcols/), kpblh_idx)
142 0 : call pbuf_add_field('wstarPBL', 'global', dtype_r8, (/pcols/), wstarPBL_idx)
143 0 : call pbuf_add_field('tkes', 'global', dtype_r8, (/pcols/), tkes_idx)
144 0 : call pbuf_add_field('went', 'global', dtype_r8, (/pcols/), went_idx)
145 : end if
146 :
147 0 : end subroutine eddy_diff_register
148 :
149 0 : subroutine eddy_diff_init(pbuf2d, ntop_eddy_in, nbot_eddy_in)
150 :
151 0 : use error_messages, only: handle_errmsg
152 : use cam_history, only: addfld, add_default, horiz_only
153 : use constituents, only: cnst_get_ind
154 : use ref_pres, only: pref_mid
155 : use diffusion_solver, only: new_fieldlist_vdiff, vdiff_select
156 : use eddy_diff, only: init_eddy_diff
157 : use physics_buffer, only: pbuf_set_field, pbuf_get_index
158 :
159 : type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer
160 : integer, intent(in) :: ntop_eddy_in ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
161 : integer, intent(in) :: nbot_eddy_in ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )
162 :
163 : character(len=128) :: errstring
164 :
165 : real(r8) :: leng_max(pver)
166 : integer :: k
167 :
168 : logical :: history_amwg
169 :
170 0 : ntop_eddy = ntop_eddy_in
171 0 : nbot_eddy = nbot_eddy_in
172 :
173 0 : do k = 1, pver
174 0 : if (pref_mid(k) <= eddy_max_bot_pressure*1.e2_r8) then
175 0 : leng_max(k) = eddy_leng_max
176 : else
177 0 : leng_max(k) = 40.e3_r8
178 : end if
179 : end do
180 :
181 0 : if (masterproc) then
182 0 : write(iulog,*)'init_eddy_diff: nturb=',nturb
183 0 : write(iulog,*)'init_eddy_diff: eddy_leng_max=',eddy_leng_max,' lbulk_max=',eddy_lbulk_max
184 0 : do k = 1,pver
185 0 : write(iulog,*)'init_eddy_diff:',k,pref_mid(k),'leng_max=',leng_max(k)
186 : end do
187 : end if
188 :
189 : call init_eddy_diff(pver, gravit, cpair, rair, zvir, &
190 : latvap, latice, ntop_eddy, nbot_eddy, karman, &
191 : eddy_lbulk_max, leng_max, &
192 0 : eddy_moist_entrain_a2l, errstring)
193 :
194 0 : call handle_errmsg(errstring, subname="init_eddy_diff")
195 :
196 : ! Set the square of the mixing lengths.
197 0 : ml2(1:ntop_eddy) = 0._r8
198 0 : do k = ntop_eddy + 1, nbot_eddy
199 0 : ml2(k) = 30.0_r8**2
200 : end do
201 0 : ml2(nbot_eddy+1:pver+1) = 0._r8
202 :
203 : ! Get fieldlists to pass to diffusion solver.
204 0 : fieldlist_wet = new_fieldlist_vdiff(1)
205 0 : fieldlist_molec = new_fieldlist_vdiff(1)
206 :
207 : call handle_errmsg(vdiff_select(fieldlist_wet,'s'), &
208 0 : subname="vdiff_select")
209 : call handle_errmsg(vdiff_select(fieldlist_wet,'q',1), &
210 0 : subname="vdiff_select")
211 : call handle_errmsg(vdiff_select(fieldlist_wet,'u'), &
212 0 : subname="vdiff_select")
213 : call handle_errmsg(vdiff_select(fieldlist_wet,'v'), &
214 0 : subname="vdiff_select")
215 :
216 : ! Cloud mass constituents
217 0 : call cnst_get_ind('CLDLIQ', ixcldliq)
218 0 : call cnst_get_ind('CLDICE', ixcldice)
219 :
220 : ! Input pbuf fields
221 0 : qrl_idx = pbuf_get_index('QRL')
222 0 : wsedl_idx = pbuf_get_index('WSEDL')
223 :
224 : ! Initialize output pbuf fields
225 0 : if (is_first_step() .and. unicon_is_on) then
226 0 : call pbuf_set_field(pbuf2d, bprod_idx, 1.0e-5_r8)
227 0 : call pbuf_set_field(pbuf2d, ipbl_idx, 0 )
228 0 : call pbuf_set_field(pbuf2d, kpblh_idx, 1 )
229 0 : call pbuf_set_field(pbuf2d, wstarPBL_idx, 0.0_r8)
230 0 : call pbuf_set_field(pbuf2d, tkes_idx, 0.0_r8)
231 0 : call pbuf_set_field(pbuf2d, went_idx, 0.0_r8)
232 : end if
233 :
234 : ! Scheme-specific default output.
235 0 : call phys_getopts(history_amwg_out=history_amwg)
236 :
237 0 : call addfld('WGUSTD', horiz_only, 'A', 'm/s', 'wind gusts from turbulence' )
238 0 : if (history_amwg) then
239 0 : call add_default( 'WGUSTD ', 1, ' ' )
240 : end if
241 :
242 : ! ------------------------------------------------------------------- !
243 : ! Writing outputs for detailed analysis of UW moist turbulence scheme !
244 : ! ------------------------------------------------------------------- !
245 :
246 0 : call addfld( 'BPROD', ['ilev'], 'A', 'm2/s3', 'Buoyancy Production' )
247 0 : call addfld( 'SFI', ['ilev'], 'A', '1', 'Interface-layer sat frac' )
248 0 : call addfld( 'SPROD', ['ilev'], 'A', 'm2/s3', 'Shear Production' )
249 :
250 :
251 0 : call addfld('UW_errorPBL',horiz_only,'A', 'm2/s', 'Error function of UW PBL')
252 0 : call addfld('UW_n2', ['lev'], 'A', 's-2', 'Buoyancy Frequency, LI')
253 0 : call addfld('UW_s2', ['lev'], 'A', 's-2', 'Shear Frequency, LI')
254 0 : call addfld('UW_ri', ['lev'], 'A', '1', 'Interface Richardson Number, I')
255 0 : call addfld('UW_sfuh', ['lev'], 'A', '1', 'Upper-Half Saturation Fraction, L')
256 0 : call addfld('UW_sflh', ['lev'], 'A', '1', 'Lower-Half Saturation Fraction, L')
257 0 : call addfld('UW_sfi', ['ilev'], 'A', '1', 'Interface Saturation Fraction, I')
258 0 : call addfld('UW_cldn', ['lev'], 'A', '1', 'Cloud Fraction, L')
259 0 : call addfld('UW_qrl', ['lev'], 'A', 'gravity W/m2', 'LW cooling rate, L')
260 0 : call addfld('UW_ql', ['lev'], 'A', 'kg/kg', 'ql(LWC), L')
261 0 : call addfld('UW_chu', ['ilev'], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chu, I')
262 0 : call addfld('UW_chs', ['ilev'], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chs, I')
263 0 : call addfld('UW_cmu', ['ilev'], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cmu, I')
264 0 : call addfld('UW_cms', ['ilev'], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cms, I')
265 0 : call addfld('UW_tke', ['ilev'], 'A', 'm2/s2', 'TKE, I')
266 0 : call addfld('UW_wcap', ['ilev'], 'A', 'm2/s2', 'Wcap, I')
267 0 : call addfld('UW_bprod', ['ilev'], 'A', 'm2/s3', 'Buoyancy production, I')
268 0 : call addfld('UW_sprod', ['ilev'], 'A', 'm2/s3', 'Shear production, I')
269 0 : call addfld('UW_kvh', ['ilev'], 'A', 'm2/s', 'Eddy diffusivity of heat, I')
270 0 : call addfld('UW_kvm', ['ilev'], 'A', 'm2/s', 'Eddy diffusivity of uv, I')
271 0 : call addfld('UW_pblh', horiz_only, 'A', 'm', 'PBLH, 1')
272 0 : call addfld('UW_pblhp', horiz_only, 'A', 'Pa', 'PBLH pressure, 1')
273 0 : call addfld('UW_tpert', horiz_only, 'A', 'K', 'Convective T excess, 1')
274 0 : call addfld('UW_qpert', horiz_only, 'A', 'kg/kg', 'Convective qt excess, I')
275 0 : call addfld('UW_wpert', horiz_only, 'A', 'm/s', 'Convective W excess, I')
276 0 : call addfld('UW_ustar', horiz_only, 'A', 'm/s', 'Surface Frictional Velocity, 1')
277 0 : call addfld('UW_tkes', horiz_only, 'A', 'm2/s2', 'Surface TKE, 1')
278 0 : call addfld('UW_minpblh',horiz_only, 'A', 'm', 'Minimum PBLH, 1')
279 0 : call addfld('UW_turbtype', ['ilev'], 'A', '1', 'Interface Turbulence Type, I')
280 0 : call addfld('UW_kbase_o', ['lev'], 'A', '1', 'Initial CL Base Exterbal Interface Index, CL')
281 0 : call addfld('UW_ktop_o', ['lev'], 'A', '1', 'Initial Top Exterbal Interface Index, CL')
282 0 : call addfld('UW_ncvfin_o',horiz_only,'A', '1', 'Initial Total Number of CL regimes, CL')
283 0 : call addfld('UW_kbase_mg', ['lev'], 'A', '1', 'kbase after merging, CL')
284 0 : call addfld('UW_ktop_mg', ['lev'], 'A', '1', 'ktop after merging, CL')
285 0 : call addfld('UW_ncvfin_mg',horiz_only,'A', '1', 'ncvfin after merging, CL')
286 0 : call addfld('UW_kbase_f', ['lev'], 'A', '1', 'Final kbase with SRCL, CL')
287 0 : call addfld('UW_ktop_f', ['lev'], 'A', '1', 'Final ktop with SRCL, CL')
288 0 : call addfld('UW_ncvfin_f',horiz_only,'A', '1', 'Final ncvfin with SRCL, CL')
289 0 : call addfld('UW_wet', ['lev'], 'A', 'm/s', 'Entrainment rate at CL top, CL')
290 0 : call addfld('UW_web', ['lev'], 'A', 'm/s', 'Entrainment rate at CL base, CL')
291 0 : call addfld('UW_jtbu', ['lev'], 'A', 'm/s2', 'Buoyancy jump across CL top, CL')
292 0 : call addfld('UW_jbbu', ['lev'], 'A', 'm/s2', 'Buoyancy jump across CL base, CL')
293 0 : call addfld('UW_evhc', ['lev'], 'A', '1', 'Evaporative enhancement factor, CL')
294 0 : call addfld('UW_jt2slv', ['lev'], 'A', 'J/kg', 'slv jump for evhc, CL')
295 0 : call addfld('UW_n2ht', ['lev'], 'A', 's-2', 'n2 at just below CL top interface, CL')
296 0 : call addfld('UW_n2hb', ['lev'], 'A', 's-2', 'n2 at just above CL base interface')
297 0 : call addfld('UW_lwp', ['lev'], 'A', 'kg/m2', 'LWP in the CL top layer, CL')
298 0 : call addfld('UW_optdepth', ['lev'], 'A', '1', 'Optical depth of the CL top layer, CL')
299 0 : call addfld('UW_radfrac', ['lev'], 'A', '1', 'Fraction of radiative cooling confined in the CL top')
300 0 : call addfld('UW_radf', ['lev'], 'A', 'm2/s3', 'Buoyancy production at the CL top by radf, I')
301 0 : call addfld('UW_wstar', ['lev'], 'A', 'm/s', 'Convective velocity, Wstar, CL')
302 0 : call addfld('UW_wstar3fact',['lev'], 'A', '1', 'Enhancement of wstar3 due to entrainment, CL')
303 0 : call addfld('UW_ebrk', ['lev'], 'A', 'm2/s2', 'CL-averaged TKE, CL')
304 0 : call addfld('UW_wbrk', ['lev'], 'A', 'm2/s2', 'CL-averaged W, CL')
305 0 : call addfld('UW_lbrk', ['lev'], 'A', 'm', 'CL internal thickness, CL')
306 0 : call addfld('UW_ricl', ['lev'], 'A', '1', 'CL-averaged Ri, CL')
307 0 : call addfld('UW_ghcl', ['lev'], 'A', '1', 'CL-averaged gh, CL')
308 0 : call addfld('UW_shcl', ['lev'], 'A', '1', 'CL-averaged sh, CL')
309 0 : call addfld('UW_smcl', ['lev'], 'A', '1', 'CL-averaged sm, CL')
310 0 : call addfld('UW_gh', ['ilev'], 'A', '1', 'gh at all interfaces, I')
311 0 : call addfld('UW_sh', ['ilev'], 'A', '1', 'sh at all interfaces, I')
312 0 : call addfld('UW_sm', ['ilev'], 'A', '1', 'sm at all interfaces, I')
313 0 : call addfld('UW_ria', ['ilev'], 'A', '1', 'ri at all interfaces, I')
314 0 : call addfld('UW_leng', ['ilev'], 'A', 'm/s', 'Turbulence length scale, I')
315 : ! For sedimentation-entrainment feedback analysis
316 0 : call addfld('UW_wsed', ['lev'], 'A', 'm/s', 'Sedimentation velocity at CL top, CL')
317 :
318 0 : end subroutine eddy_diff_init
319 :
320 0 : subroutine eddy_diff_tend(state, pbuf, cam_in, &
321 : ztodt, p, tint, rhoi, cldn, wstarent, &
322 : kvm_in, kvh_in, ksrftms, dragblj,tauresx, tauresy, &
323 : rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, &
324 : tke, sprod, sfi)
325 :
326 0 : use physics_types, only: physics_state
327 : use camsrfexch, only: cam_in_t
328 : use coords_1d, only: Coords1D
329 :
330 : type(physics_state), intent(in) :: state
331 : type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
332 : type(cam_in_t), intent(in) :: cam_in
333 : real(r8), intent(in) :: ztodt
334 : type(Coords1D), intent(in) :: p
335 : real(r8), intent(in) :: tint(pcols,pver+1)
336 : real(r8), intent(in) :: rhoi(pcols,pver+1)
337 : real(r8), intent(in) :: cldn(pcols,pver)
338 : logical, intent(in) :: wstarent
339 : real(r8), intent(in) :: kvm_in(pcols,pver+1)
340 : real(r8), intent(in) :: kvh_in(pcols,pver+1)
341 : real(r8), intent(in) :: ksrftms(pcols)
342 : real(r8), intent(in) :: dragblj(pcols,pver) ! Drag profile from Beljaars SGO form drag [ 1/s ]
343 : real(r8), intent(inout) :: tauresx(pcols)
344 : real(r8), intent(inout) :: tauresy(pcols)
345 : real(r8), intent(out) :: rrho(pcols)
346 : real(r8), intent(out) :: ustar(pcols)
347 : real(r8), intent(out) :: pblh(pcols)
348 : real(r8), intent(out) :: kvm(pcols,pver+1)
349 : real(r8), intent(out) :: kvh(pcols,pver+1)
350 : real(r8), intent(out) :: kvq(pcols,pver+1)
351 : real(r8), intent(out) :: cgh(pcols,pver+1)
352 : real(r8), intent(out) :: cgs(pcols,pver+1)
353 : real(r8), intent(out) :: tpert(pcols)
354 : real(r8), intent(out) :: qpert(pcols)
355 : real(r8), intent(out) :: tke(pcols,pver+1)
356 : real(r8), intent(out) :: sprod(pcols,pver+1)
357 : real(r8), intent(out) :: sfi(pcols,pver+1)
358 :
359 : integer :: i, k
360 :
361 : call compute_eddy_diff( pbuf, state%lchnk , &
362 : pcols , pver , state%ncol , state%t , tint, state%q(:,:,1) , ztodt , &
363 : state%q(:,:,ixcldliq) , state%q(:,:,ixcldice) , &
364 : state%s , p , rhoi, cldn , &
365 : state%zm , state%zi , state%pmid , state%pint , state%u , state%v , &
366 : cam_in%wsx, cam_in%wsy , cam_in%shf , cam_in%cflx(:,1) , wstarent , &
367 : rrho , ustar , pblh , kvm_in , kvh_in , kvm , &
368 : kvh , kvq , cgh , &
369 : cgs , tpert , qpert , tke , &
370 : sprod , sfi , &
371 0 : tauresx , tauresy , ksrftms , dragblj )
372 :
373 : ! The diffusivities from diag_TKE can be much larger than from HB in the free
374 : ! troposphere and upper atmosphere. These seem to be larger than observations,
375 : ! and in WACCM the gw_drag code is already applying an eddy diffusivity in the
376 : ! upper atmosphere. Optionally, adjust the diffusivities in the free troposphere
377 : ! or the upper atmosphere.
378 : !
379 : ! NOTE: Further investigation should be done as to why the diffusivities are
380 : ! larger in diag_TKE.
381 0 : if ((kv_freetrop_scale /= 1._r8) .or. ((kv_top_scale /= 1._r8) .and. (kv_top_pressure > 0._r8))) then
382 0 : do i = 1, state%ncol
383 0 : do k = 1, pverp
384 : ! Outside of the boundary layer?
385 0 : if (state%zi(i,k) > pblh(i)) then
386 : ! In the upper atmosphere?
387 0 : if (state%pint(i,k) <= kv_top_pressure) then
388 0 : kvh(i,k) = kvh(i,k) * kv_top_scale
389 0 : kvm(i,k) = kvm(i,k) * kv_top_scale
390 0 : kvq(i,k) = kvq(i,k) * kv_top_scale
391 : else
392 0 : kvh(i,k) = kvh(i,k) * kv_freetrop_scale
393 0 : kvm(i,k) = kvm(i,k) * kv_freetrop_scale
394 0 : kvq(i,k) = kvq(i,k) * kv_freetrop_scale
395 : end if
396 : else
397 : exit
398 : end if
399 : end do
400 : end do
401 : end if
402 :
403 0 : end subroutine eddy_diff_tend
404 :
405 : !=============================================================================== !
406 : ! !
407 : !=============================================================================== !
408 :
409 0 : subroutine compute_eddy_diff( pbuf, lchnk , &
410 0 : pcols , pver , ncol , t , tint, qv , ztodt , &
411 0 : ql , qi , s , p , rhoi, cldn , &
412 0 : z , zi , pmid , pi , u , v , &
413 0 : taux , tauy , shflx , qflx , wstarent , rrho , &
414 0 : ustar , pblh , kvm_in , kvh_in , kvm_out , kvh_out , kvq , &
415 0 : cgh , cgs , tpert , qpert , tke , &
416 0 : sprod , sfi , &
417 0 : tauresx, tauresy, ksrftms, dragblj )
418 :
419 : !-------------------------------------------------------------------- !
420 : ! Purpose: Interface to compute eddy diffusivities. !
421 : ! Eddy diffusivities are calculated in a fully implicit way !
422 : ! through iteration process. !
423 : ! Author: Sungsu Park. August. 2006. !
424 : ! May. 2008. !
425 : !-------------------------------------------------------------------- !
426 :
427 0 : use diffusion_solver, only: compute_vdiff
428 : use cam_history, only: outfld
429 : use phys_debug_util, only: phys_debug_col
430 : use air_composition, only: cpairv
431 : use pbl_utils, only: calc_ustar, austausch_atm
432 : use error_messages, only: handle_errmsg
433 : use coords_1d, only: Coords1D
434 : use wv_saturation, only: qsat
435 : use eddy_diff, only: trbintd, caleddy
436 : use physics_buffer, only: pbuf_get_field
437 :
438 : ! --------------- !
439 : ! Input Variables !
440 : ! --------------- !
441 :
442 : type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
443 : integer, intent(in) :: lchnk
444 : integer, intent(in) :: pcols ! Number of atmospheric columns [ # ]
445 : integer, intent(in) :: pver ! Number of atmospheric layers [ # ]
446 : integer, intent(in) :: ncol ! Number of atmospheric columns [ # ]
447 : logical, intent(in) :: wstarent ! .true. means use the 'wstar' entrainment closure.
448 : real(r8), intent(in) :: ztodt ! Physics integration time step 2 delta-t [ s ]
449 : real(r8), intent(in) :: t(pcols,pver) ! Temperature [ K ]
450 : real(r8), intent(in) :: tint(pcols,pver+1) ! Temperature defined on interfaces [ K ]
451 : real(r8), intent(in) :: qv(pcols,pver) ! Water vapor specific humidity [ kg/kg ]
452 : real(r8), intent(in) :: ql(pcols,pver) ! Liquid water specific humidity [ kg/kg ]
453 : real(r8), intent(in) :: qi(pcols,pver) ! Ice specific humidity [ kg/kg ]
454 : real(r8), intent(in) :: s(pcols,pver) ! Dry static energy [ J/kg ]
455 : type(Coords1D), intent(in) :: p ! Pressure coordinates for solver [ Pa ]
456 : real(r8), intent(in) :: rhoi(pcols,pver+1) ! Density at interfaces [ kg/m3 ]
457 : real(r8), intent(in) :: cldn(pcols,pver) ! Stratiform cloud fraction [ fraction ]
458 : real(r8), intent(in) :: z(pcols,pver) ! Layer mid-point height above surface [ m ]
459 : real(r8), intent(in) :: zi(pcols,pver+1) ! Interface height above surface [ m ]
460 : real(r8), intent(in) :: pmid(pcols,pver) ! Layer mid-point pressure [ Pa ]
461 : real(r8), intent(in) :: pi(pcols,pver+1) ! Interface pressure [ Pa ]
462 : real(r8), intent(in) :: u(pcols,pver) ! Zonal velocity [ m/s ]
463 : real(r8), intent(in) :: v(pcols,pver) ! Meridional velocity [ m/s ]
464 : real(r8), intent(in) :: taux(pcols) ! Zonal wind stress at surface [ N/m2 ]
465 : real(r8), intent(in) :: tauy(pcols) ! Meridional wind stress at surface [ N/m2 ]
466 : real(r8), intent(in) :: shflx(pcols) ! Sensible heat flux at surface [ unit ? ]
467 : real(r8), intent(in) :: qflx(pcols) ! Water vapor flux at surface [ unit ? ]
468 : real(r8), intent(in) :: kvm_in(pcols,pver+1) ! kvm saved from last timestep [ m2/s ]
469 : real(r8), intent(in) :: kvh_in(pcols,pver+1) ! kvh saved from last timestep [ m2/s ]
470 : real(r8), intent(in) :: ksrftms(pcols) ! Surface drag coefficient of turbulent mountain stress [ unit ? ]
471 : real(r8), intent(in) :: dragblj(pcols,pver) ! Drag profile from Beljaars SGO form drag [ 1/s ]
472 :
473 : ! ---------------- !
474 : ! Output Variables !
475 : ! ---------------- !
476 :
477 : real(r8), intent(out) :: kvm_out(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
478 : real(r8), intent(out) :: kvh_out(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
479 : real(r8), intent(out) :: kvq(pcols,pver+1) ! Eddy diffusivity for constituents, moisture and tracers [ m2/s ]
480 : ! (note not having '_out')
481 : real(r8), intent(out) :: rrho(pcols) ! Reciprocal of density at the lowest layer
482 : real(r8), intent(out) :: ustar(pcols) ! Surface friction velocity [ m/s ]
483 : real(r8), intent(out) :: pblh(pcols) ! PBL top height [ m ]
484 : real(r8), intent(out) :: cgh(pcols,pver+1) ! Counter-gradient term for heat [ J/kg/m ]
485 : real(r8), intent(out) :: cgs(pcols,pver+1) ! Counter-gradient star [ cg/flux ]
486 : real(r8), intent(out) :: tpert(pcols) ! Convective temperature excess [ K ]
487 : real(r8), intent(out) :: qpert(pcols) ! Convective humidity excess [ kg/kg ]
488 : real(r8), intent(out) :: tke(pcols,pver+1) ! Turbulent kinetic energy [ m2/s2 ]
489 : real(r8), intent(out) :: sprod(pcols,pver+1) ! Shear production [ m2/s3 ]
490 : real(r8), intent(out) :: sfi(pcols,pver+1) ! Interfacial layer saturation fraction [ fraction ]
491 :
492 : ! ---------------------- !
493 : ! Input-Output Variables !
494 : ! ---------------------- !
495 :
496 : real(r8), intent(inout) :: tauresx(pcols) ! Residual stress to be added in vdiff to correct for turb
497 : real(r8), intent(inout) :: tauresy(pcols) ! Stress mismatch between sfc and atm accumulated in prior timesteps
498 :
499 : ! -------------- !
500 : ! pbuf Variables !
501 : ! -------------- !
502 :
503 0 : real(r8), pointer :: qrl(:,:) ! LW radiative cooling rate
504 0 : real(r8), pointer :: wsedl(:,:) ! Sedimentation velocity
505 : ! of stratiform liquid cloud droplet [ m/s ]
506 :
507 0 : real(r8), pointer :: bprod(:,:) ! Buoyancy production of tke [ m2/s3 ]
508 0 : integer(i4), pointer :: ipbl(:) ! If 1, PBL is CL, while if 0, PBL is STL.
509 0 : integer(i4), pointer :: kpblh(:) ! Layer index containing PBL top within or at the base interface
510 0 : real(r8), pointer :: wstarPBL(:) ! Convective velocity within PBL [ m/s ]
511 0 : real(r8), pointer :: tkes(:) ! TKE at surface interface [ m2/s2 ]
512 0 : real(r8), pointer :: went(:) ! Entrainment rate at the PBL top interface [ m/s ]
513 :
514 : ! --------------- !
515 : ! Local Variables !
516 : ! --------------- !
517 :
518 : integer icol
519 : integer i, k, iturb, status
520 :
521 : character(2048) :: warnstring ! Warning(s) to print
522 : character(128) :: errstring ! Error message
523 :
524 0 : real(r8) :: kvf(pcols,pver+1) ! Free atmospheric eddy diffusivity [ m2/s ]
525 0 : real(r8) :: kvm(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
526 0 : real(r8) :: kvh(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
527 : real(r8) :: kvm_preo(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
528 : real(r8) :: kvh_preo(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
529 : real(r8) :: kvm_pre(pcols,pver+1) ! Eddy diffusivity for momentum [ m2/s ]
530 : real(r8) :: kvh_pre(pcols,pver+1) ! Eddy diffusivity for heat [ m2/s ]
531 0 : real(r8) :: errorPBL(pcols) ! Error function showing whether PBL produced convergent solution or not.
532 : ! [ unit ? ]
533 0 : real(r8) :: s2(pcols,pver) ! Shear squared, defined at interfaces except surface [ s-2 ]
534 0 : real(r8) :: n2(pcols,pver) ! Buoyancy frequency, defined at interfaces except surface [ s-2 ]
535 0 : real(r8) :: ri(pcols,pver) ! Richardson number, 'n2/s2', defined at interfaces except surface [ s-2 ]
536 0 : real(r8) :: pblhp(pcols) ! PBL top pressure [ Pa ]
537 0 : real(r8) :: minpblh(pcols) ! Minimum PBL height based on surface stress
538 :
539 0 : real(r8) :: qt(pcols,pver) ! Total specific humidity [ kg/kg ]
540 0 : real(r8) :: sfuh(pcols,pver) ! Saturation fraction in upper half-layer [ fraction ]
541 0 : real(r8) :: sflh(pcols,pver) ! Saturation fraction in lower half-layer [ fraction ]
542 0 : real(r8) :: sl(pcols,pver) ! Liquid water static energy [ J/kg ]
543 0 : real(r8) :: slv(pcols,pver) ! Liquid water virtual static energy [ J/kg ]
544 0 : real(r8) :: slslope(pcols,pver) ! Slope of 'sl' in each layer
545 0 : real(r8) :: qtslope(pcols,pver) ! Slope of 'qt' in each layer
546 0 : real(r8) :: qvfd(pcols,pver) ! Specific humidity for diffusion [ kg/kg ]
547 0 : real(r8) :: tfd(pcols,pver) ! Temperature for diffusion [ K ]
548 0 : real(r8) :: slfd(pcols,pver) ! Liquid static energy [ J/kg ]
549 0 : real(r8) :: qtfd(pcols,pver) ! Total specific humidity [ kg/kg ]
550 0 : real(r8) :: qlfd(pcols,pver) ! Liquid water specific humidity for diffusion [ kg/kg ]
551 0 : real(r8) :: ufd(pcols,pver) ! U-wind for diffusion [ m/s ]
552 0 : real(r8) :: vfd(pcols,pver) ! V-wind for diffusion [ m/s ]
553 :
554 : ! Buoyancy coefficients : w'b' = ch * w'sl' + cm * w'qt'
555 :
556 0 : real(r8) :: chu(pcols,pver+1) ! Heat buoyancy coef for dry states, defined at each interface, finally.
557 0 : real(r8) :: chs(pcols,pver+1) ! Heat buoyancy coef for sat states, defined at each interface, finally.
558 0 : real(r8) :: cmu(pcols,pver+1) ! Moisture buoyancy coef for dry states,
559 : ! defined at each interface, finally.
560 0 : real(r8) :: cms(pcols,pver+1) ! Moisture buoyancy coef for sat states,
561 : ! defined at each interface, finally.
562 :
563 0 : real(r8) :: jnk1d(pcols)
564 0 : real(r8) :: jnk2d(pcols,pver+1)
565 0 : real(r8) :: zero(pcols)
566 0 : real(r8) :: zero2d(pcols,pver+1)
567 : real(r8) :: es ! Saturation vapor pressure
568 : real(r8) :: qs ! Saturation specific humidity
569 : real(r8) :: ep2, templ, temps
570 :
571 : ! ------------------------------- !
572 : ! Variables for diagnostic output !
573 : ! ------------------------------- !
574 :
575 0 : real(r8) :: wpert(pcols) ! Turbulent velocity excess [ m/s ]
576 :
577 0 : real(r8) :: kbase_o(pcols,ncvmax) ! Original external base interface index of CL from 'exacol'
578 0 : real(r8) :: ktop_o(pcols,ncvmax) ! Original external top interface index of CL from 'exacol'
579 0 : real(r8) :: ncvfin_o(pcols) ! Original number of CLs from 'exacol'
580 0 : real(r8) :: kbase_mg(pcols,ncvmax) ! 'kbase' after extending-merging from 'zisocl'
581 0 : real(r8) :: ktop_mg(pcols,ncvmax) ! 'ktop' after extending-merging from 'zisocl'
582 0 : real(r8) :: ncvfin_mg(pcols) ! 'ncvfin' after extending-merging from 'zisocl'
583 0 : real(r8) :: kbase_f(pcols,ncvmax) ! Final 'kbase' after extending-merging & including SRCL
584 0 : real(r8) :: ktop_f(pcols,ncvmax) ! Final 'ktop' after extending-merging & including SRCL
585 0 : real(r8) :: ncvfin_f(pcols) ! Final 'ncvfin' after extending-merging & including SRCL
586 0 : real(r8) :: wet(pcols,ncvmax) ! Entrainment rate at the CL top [ m/s ]
587 0 : real(r8) :: web(pcols,ncvmax) ! Entrainment rate at the CL base [ m/s ].
588 : ! Set to zero if CL is based at surface.
589 0 : real(r8) :: jtbu(pcols,ncvmax) ! Buoyancy jump across the CL top [ m/s2 ]
590 0 : real(r8) :: jbbu(pcols,ncvmax) ! Buoyancy jump across the CL base [ m/s2 ]
591 0 : real(r8) :: evhc(pcols,ncvmax) ! Evaporative enhancement factor at the CL top
592 0 : real(r8) :: jt2slv(pcols,ncvmax) ! Jump of slv ( across two layers ) at CL top used only for evhc [ J/kg ]
593 0 : real(r8) :: n2ht(pcols,ncvmax) ! n2 defined at the CL top interface but using
594 : ! sfuh(kt) instead of sfi(kt) [ s-2 ]
595 0 : real(r8) :: n2hb(pcols,ncvmax) ! n2 defined at the CL base interface but using
596 : ! sflh(kb-1) instead of sfi(kb) [ s-2 ]
597 0 : real(r8) :: lwp(pcols,ncvmax) ! LWP in the CL top layer [ kg/m2 ]
598 0 : real(r8) :: opt_depth(pcols,ncvmax) ! Optical depth of the CL top layer
599 0 : real(r8) :: radinvfrac(pcols,ncvmax) ! Fraction of radiative cooling confined in the top portion of CL top layer
600 0 : real(r8) :: radf(pcols,ncvmax) ! Buoyancy production at the CL top due to LW radiative cooling [ m2/s3 ]
601 0 : real(r8) :: wstar(pcols,ncvmax) ! Convective velocity in each CL [ m/s ]
602 0 : real(r8) :: wstar3fact(pcols,ncvmax) ! Enhancement of 'wstar3' due to entrainment (inverse) [ no unit ]
603 0 : real(r8) :: ebrk(pcols,ncvmax) ! Net mean TKE of CL including entrainment effect [ m2/s2 ]
604 0 : real(r8) :: wbrk(pcols,ncvmax) ! Net mean normalized TKE (W) of CL,
605 : ! 'ebrk/b1' including entrainment effect [ m2/s2 ]
606 0 : real(r8) :: lbrk(pcols,ncvmax) ! Energetic internal thickness of CL [m]
607 0 : real(r8) :: ricl(pcols,ncvmax) ! CL internal mean Richardson number
608 0 : real(r8) :: ghcl(pcols,ncvmax) ! Half of normalized buoyancy production of CL
609 0 : real(r8) :: shcl(pcols,ncvmax) ! Galperin instability function of heat-moisture of CL
610 0 : real(r8) :: smcl(pcols,ncvmax) ! Galperin instability function of mementum of CL
611 0 : real(r8) :: ghi(pcols,pver+1) ! Half of normalized buoyancy production at all interfaces
612 0 : real(r8) :: shi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces
613 0 : real(r8) :: smi(pcols,pver+1) ! Galperin instability function of heat-moisture at all interfaces
614 0 : real(r8) :: rii(pcols,pver+1) ! Interfacial Richardson number defined at all interfaces
615 0 : real(r8) :: lengi(pcols,pver+1) ! Turbulence length scale at all interfaces [ m ]
616 0 : real(r8) :: wcap(pcols,pver+1) ! Normalized TKE at all interfaces [ m2/s2 ]
617 : ! For sedimentation-entrainment feedback
618 0 : real(r8) :: wsed(pcols,ncvmax) ! Sedimentation velocity at the top of each CL [ m/s ]
619 :
620 0 : integer(i4) :: turbtype(pcols,pver+1) ! Turbulence type identifier at all interfaces [ no unit ]
621 :
622 : ! ---------- !
623 : ! Parameters !
624 : ! ---------- !
625 :
626 : logical, parameter :: use_kvf = .false. ! .true. (.false.) : initialize kvh/kvm = kvf ( 0. )
627 : real(r8), parameter :: lambda = 0.5_r8 ! Under-relaxation factor ( 0 < lambda =< 1 )
628 :
629 : ! ---------- !
630 : ! Initialize !
631 : ! ---------- !
632 :
633 0 : zero(:) = 0._r8
634 0 : zero2d(:,:) = 0._r8
635 :
636 : ! ---------------------------------------------- !
637 : ! Get LW radiative heating out of physics buffer !
638 : ! ---------------------------------------------- !
639 0 : call pbuf_get_field(pbuf, qrl_idx, qrl)
640 0 : call pbuf_get_field(pbuf, wsedl_idx, wsedl)
641 :
642 : ! These fields are put into the pbuf for UNICON only.
643 0 : if (unicon_is_on) then
644 0 : call pbuf_get_field(pbuf, bprod_idx, bprod)
645 0 : call pbuf_get_field(pbuf, ipbl_idx, ipbl)
646 0 : call pbuf_get_field(pbuf, kpblh_idx, kpblh)
647 0 : call pbuf_get_field(pbuf, wstarPBL_idx, wstarPBL)
648 0 : call pbuf_get_field(pbuf, tkes_idx, tkes)
649 0 : call pbuf_get_field(pbuf, went_idx, went)
650 : else
651 0 : allocate(bprod(pcols,pverp), ipbl(pcols), kpblh(pcols), wstarPBL(pcols), tkes(pcols), went(pcols))
652 : end if
653 :
654 : ! ----------------------- !
655 : ! Main Computation Begins !
656 : ! ----------------------- !
657 :
658 0 : ufd(:ncol,:) = u(:ncol,:)
659 0 : vfd(:ncol,:) = v(:ncol,:)
660 0 : tfd(:ncol,:) = t(:ncol,:)
661 0 : qvfd(:ncol,:) = qv(:ncol,:)
662 0 : qlfd(:ncol,:) = ql(:ncol,:)
663 :
664 0 : do iturb = 1, nturb
665 :
666 : ! Total stress includes 'tms'.
667 : ! Here, in computing 'tms', we can use either iteratively changed 'ufd,vfd' or the
668 : ! initially given 'u,v' to the PBL scheme. Note that normal stress, 'taux, tauy'
669 : ! are not changed by iteration. In order to treat 'tms' in a fully implicit way,
670 : ! I am using updated wind, here.
671 :
672 : ! Compute ustar
673 : call calc_ustar( ncol, tfd(:ncol,pver), pmid(:ncol,pver), &
674 0 : taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
675 : tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
676 0 : rrho(:ncol), ustar(:ncol))
677 0 : minpblh(:ncol) = 100.0_r8 * ustar(:ncol) ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.
678 :
679 : ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)
680 :
681 : call trbintd( &
682 : pcols , pver , ncol , z , ufd , vfd , tfd , pmid , &
683 : s2 , n2 , ri , zi , pi , cldn , qtfd , qvfd , &
684 : qlfd , qi , sfi , sfuh , sflh , slfd , slv , slslope , &
685 0 : qtslope , chs , chu , cms , cmu )
686 :
687 : ! Save initial (i.e., before iterative diffusion) profile of (qt,sl) at each iteration.
688 : ! Only necessary for (qt,sl) not (u,v) because (qt,sl) are newly calculated variables.
689 :
690 0 : if( iturb == 1 ) then
691 0 : qt(:ncol,:) = qtfd(:ncol,:)
692 0 : sl(:ncol,:) = slfd(:ncol,:)
693 : endif
694 :
695 : ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
696 : if (use_kvf) then
697 : call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, &
698 : ml2, ri, s2, kvf )
699 : else
700 0 : kvf = 0._r8
701 : end if
702 :
703 : ! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number
704 : ! This is necessary for 'wstar-based' entrainment closure.
705 :
706 0 : if( iturb == 1 ) then
707 0 : if( is_first_step() ) then
708 : ! First iteration of first model timestep : Use free tropospheric value or zero.
709 0 : kvh(:ncol,:) = kvf(:ncol,:)
710 0 : kvm(:ncol,:) = kvf(:ncol,:)
711 : else
712 : ! First iteration on any model timestep except the first : Use value from previous timestep
713 0 : kvh(:ncol,:) = kvh_in(:ncol,:)
714 0 : kvm(:ncol,:) = kvm_in(:ncol,:)
715 : endif
716 : else
717 : ! Not the first iteration : Use from previous iteration
718 0 : kvh(:ncol,:) = kvh_out(:ncol,:)
719 0 : kvm(:ncol,:) = kvm_out(:ncol,:)
720 : endif
721 :
722 : ! Calculate eddy diffusivity (kvh_out,kvm_out) and (tke,bprod,sprod) using
723 : ! a given (kvh,kvm) which are used only for initializing (bprod,sprod) at
724 : ! the first part of caleddy. (bprod,sprod) are fully updated at the end of
725 : ! caleddy after calculating (kvh_out,kvm_out)
726 :
727 : call caleddy( pcols , pver , ncol , &
728 : slfd , qtfd , qlfd , slv ,ufd , &
729 : vfd , pi , z , zi , &
730 : qflx , shflx , slslope , qtslope , &
731 : chu , chs , cmu , cms ,sfuh , &
732 : sflh , n2 , s2 , ri ,rrho , &
733 : pblh , ustar , &
734 : kvh , kvm , kvh_out , kvm_out , &
735 : tpert , qpert , qrl , kvf , tke , &
736 : wstarent , bprod , sprod , minpblh , wpert , &
737 : tkes , went , turbtype , &
738 : kbase_o , ktop_o , ncvfin_o , &
739 : kbase_mg , ktop_mg , ncvfin_mg , &
740 : kbase_f , ktop_f , ncvfin_f , &
741 : wet , web , jtbu , jbbu , &
742 : evhc , jt2slv , n2ht , n2hb , &
743 : lwp , opt_depth , radinvfrac, radf , &
744 : wstar , wstar3fact, &
745 : ebrk , wbrk , lbrk , ricl , ghcl , &
746 : shcl , smcl , ghi , shi , smi , &
747 : rii , lengi , wcap , pblhp , cldn , &
748 : ipbl , kpblh , wsedl , wsed, &
749 0 : warnstring, errstring)
750 :
751 0 : if (trim(warnstring) /= "") then
752 0 : write(iulog,*) "eddy_diff_cam: Messages from caleddy follow."
753 0 : write(iulog,*) warnstring
754 : end if
755 :
756 0 : call handle_errmsg(errstring, subname="caleddy")
757 :
758 : ! Calculate errorPBL to check whether PBL produced convergent solutions or not.
759 :
760 0 : if( iturb == nturb ) then
761 0 : do i = 1, ncol
762 0 : errorPBL(i) = 0._r8
763 0 : do k = 1, pver
764 0 : errorPBL(i) = errorPBL(i) + ( kvh(i,k) - kvh_out(i,k) )**2
765 : end do
766 0 : errorPBL(i) = sqrt(errorPBL(i)/pver)
767 : end do
768 : end if
769 :
770 : ! Eddy diffusivities which will be used for the initialization of (bprod,
771 : ! sprod) in 'caleddy' at the next iteration step.
772 :
773 0 : if( iturb > 1 .and. iturb < nturb ) then
774 0 : kvm_out(:ncol,:) = lambda * kvm_out(:ncol,:) + ( 1._r8 - lambda ) * kvm(:ncol,:)
775 0 : kvh_out(:ncol,:) = lambda * kvh_out(:ncol,:) + ( 1._r8 - lambda ) * kvh(:ncol,:)
776 : endif
777 :
778 : ! Set nonlocal terms to zero for flux diagnostics, since not used by caleddy.
779 :
780 0 : cgh(:ncol,:) = 0._r8
781 0 : cgs(:ncol,:) = 0._r8
782 :
783 0 : if( iturb < nturb ) then
784 :
785 : ! Each time we diffuse the original state
786 :
787 0 : slfd(:ncol,:) = sl(:ncol,:)
788 0 : qtfd(:ncol,:) = qt(:ncol,:)
789 0 : ufd(:ncol,:) = u(:ncol,:)
790 0 : vfd(:ncol,:) = v(:ncol,:)
791 :
792 : ! Diffuse initial profile of each time step using a given (kvh_out,kvm_out)
793 : ! In the below 'compute_vdiff', (slfd,qtfd,ufd,vfd) are 'inout' variables.
794 :
795 : call compute_vdiff( lchnk , &
796 : pcols , pver , 1 , ncol , tint, &
797 : p , t , rhoi, ztodt , taux , &
798 : tauy , shflx , qflx , &
799 : kvh_out , kvm_out , kvh_out , cgs , cgh , &
800 : zi , ksrftms , dragblj , &
801 : zero , fieldlist_wet, fieldlist_molec, &
802 : ufd , vfd , qtfd , slfd , &
803 : jnk1d , jnk1d , jnk2d , jnk1d , errstring , &
804 : tauresx , tauresy , 0 , cpairv(:,:,lchnk), zero, &
805 0 : .false., .false. )
806 :
807 : call handle_errmsg(errstring, subname="compute_vdiff", &
808 0 : extra_msg="compute_vdiff called from eddy_diff_cam")
809 :
810 : ! Retrieve (tfd,qvfd,qlfd) from (slfd,qtfd) in order to
811 : ! use 'trbintd' at the next iteration.
812 :
813 0 : do k = 1, pver
814 0 : do i = 1, ncol
815 : ! ----------------------------------------------------- !
816 : ! Compute the condensate 'qlfd' in the updated profiles !
817 : ! ----------------------------------------------------- !
818 : ! Option.1 : Assume grid-mean condensate is homogeneously diffused by the moist turbulence scheme.
819 : ! This should be used if 'pseudodiff = .false.' in vertical_diffusion.F90.
820 : ! Modification : Need to be check whether below is correct in the presence of ice, qi.
821 : ! I should understand why the variation of ice, qi is neglected during diffusion.
822 0 : templ = ( slfd(i,k) - gravit*z(i,k) ) / cpair
823 0 : call qsat( templ, pmid(i,k), es, qs)
824 0 : ep2 = .622_r8
825 0 : temps = templ + ( qtfd(i,k) - qs ) / ( cpair / latvap + latvap * qs / ( rair * templ**2 ) )
826 0 : call qsat( temps, pmid(i,k), es, qs)
827 0 : qlfd(i,k) = max( qtfd(i,k) - qi(i,k) - qs ,0._r8 )
828 : ! Option.2 : Assume condensate is not diffused by the moist turbulence scheme.
829 : ! This should bs used if 'pseudodiff = .true.' in vertical_diffusion.F90.
830 : ! qlfd(i,k) = ql(i,k)
831 : ! ----------------------------- !
832 : ! Compute the other 'qvfd, tfd' !
833 : ! ----------------------------- !
834 0 : qvfd(i,k) = max( 0._r8, qtfd(i,k) - qi(i,k) - qlfd(i,k) )
835 0 : tfd(i,k) = ( slfd(i,k) + latvap * qlfd(i,k) + (latvap+latice) * qi(i,k) - gravit*z(i,k)) / cpair
836 : end do
837 : end do
838 : endif
839 :
840 : end do ! End of 'iturb' iteration
841 :
842 0 : kvq(:ncol,:) = kvh_out(:ncol,:)
843 :
844 : ! Compute 'wstar' within the PBL for use in the future convection scheme.
845 :
846 0 : do i = 1, ncol
847 0 : if(ipbl(i) == 1) then
848 0 : wstarPBL(i) = max( 0._r8, wstar(i,1) )
849 : else
850 0 : wstarPBL(i) = 0._r8
851 : endif
852 : end do
853 :
854 : ! --------------------------------------------------------------- !
855 : ! Writing for detailed diagnostic analysis of UW moist PBL scheme !
856 : ! --------------------------------------------------------------- !
857 :
858 0 : call outfld( 'WGUSTD' , wpert, pcols, lchnk )
859 :
860 0 : call outfld( 'BPROD ', bprod, pcols, lchnk )
861 0 : call outfld( 'SPROD ', sprod, pcols, lchnk )
862 0 : call outfld( 'SFI ', sfi, pcols, lchnk )
863 :
864 0 : call outfld( 'UW_errorPBL', errorPBL, pcols, lchnk )
865 :
866 0 : call outfld( 'UW_n2', n2, pcols, lchnk )
867 0 : call outfld( 'UW_s2', s2, pcols, lchnk )
868 0 : call outfld( 'UW_ri', ri, pcols, lchnk )
869 :
870 0 : call outfld( 'UW_sfuh', sfuh, pcols, lchnk )
871 0 : call outfld( 'UW_sflh', sflh, pcols, lchnk )
872 0 : call outfld( 'UW_sfi', sfi, pcols, lchnk )
873 :
874 0 : call outfld( 'UW_cldn', cldn, pcols, lchnk )
875 0 : call outfld( 'UW_qrl', qrl, pcols, lchnk )
876 0 : call outfld( 'UW_ql', qlfd, pcols, lchnk )
877 :
878 0 : call outfld( 'UW_chu', chu, pcols, lchnk )
879 0 : call outfld( 'UW_chs', chs, pcols, lchnk )
880 0 : call outfld( 'UW_cmu', cmu, pcols, lchnk )
881 0 : call outfld( 'UW_cms', cms, pcols, lchnk )
882 :
883 0 : call outfld( 'UW_tke', tke, pcols, lchnk )
884 0 : call outfld( 'UW_wcap', wcap, pcols, lchnk )
885 0 : call outfld( 'UW_bprod', bprod, pcols, lchnk )
886 0 : call outfld( 'UW_sprod', sprod, pcols, lchnk )
887 :
888 0 : call outfld( 'UW_kvh', kvh_out, pcols, lchnk )
889 0 : call outfld( 'UW_kvm', kvm_out, pcols, lchnk )
890 :
891 0 : call outfld( 'UW_pblh', pblh, pcols, lchnk )
892 0 : call outfld( 'UW_pblhp', pblhp, pcols, lchnk )
893 0 : call outfld( 'UW_tpert', tpert, pcols, lchnk )
894 0 : call outfld( 'UW_qpert', qpert, pcols, lchnk )
895 0 : call outfld( 'UW_wpert', wpert, pcols, lchnk )
896 :
897 0 : call outfld( 'UW_ustar', ustar, pcols, lchnk )
898 0 : call outfld( 'UW_tkes', tkes, pcols, lchnk )
899 0 : call outfld( 'UW_minpblh', minpblh, pcols, lchnk )
900 :
901 0 : call outfld( 'UW_turbtype', real(turbtype,r8), pcols, lchnk )
902 :
903 0 : call outfld( 'UW_kbase_o', kbase_o, pcols, lchnk )
904 0 : call outfld( 'UW_ktop_o', ktop_o, pcols, lchnk )
905 0 : call outfld( 'UW_ncvfin_o', ncvfin_o, pcols, lchnk )
906 :
907 0 : call outfld( 'UW_kbase_mg', kbase_mg, pcols, lchnk )
908 0 : call outfld( 'UW_ktop_mg', ktop_mg, pcols, lchnk )
909 0 : call outfld( 'UW_ncvfin_mg', ncvfin_mg, pcols, lchnk )
910 :
911 0 : call outfld( 'UW_kbase_f', kbase_f, pcols, lchnk )
912 0 : call outfld( 'UW_ktop_f', ktop_f, pcols, lchnk )
913 0 : call outfld( 'UW_ncvfin_f', ncvfin_f, pcols, lchnk )
914 :
915 0 : call outfld( 'UW_wet', wet, pcols, lchnk )
916 0 : call outfld( 'UW_web', web, pcols, lchnk )
917 0 : call outfld( 'UW_jtbu', jtbu, pcols, lchnk )
918 0 : call outfld( 'UW_jbbu', jbbu, pcols, lchnk )
919 0 : call outfld( 'UW_evhc', evhc, pcols, lchnk )
920 0 : call outfld( 'UW_jt2slv', jt2slv, pcols, lchnk )
921 0 : call outfld( 'UW_n2ht', n2ht, pcols, lchnk )
922 0 : call outfld( 'UW_n2hb', n2hb, pcols, lchnk )
923 0 : call outfld( 'UW_lwp', lwp, pcols, lchnk )
924 0 : call outfld( 'UW_optdepth', opt_depth, pcols, lchnk )
925 0 : call outfld( 'UW_radfrac', radinvfrac, pcols, lchnk )
926 0 : call outfld( 'UW_radf', radf, pcols, lchnk )
927 0 : call outfld( 'UW_wstar', wstar, pcols, lchnk )
928 0 : call outfld( 'UW_wstar3fact', wstar3fact, pcols, lchnk )
929 0 : call outfld( 'UW_ebrk', ebrk, pcols, lchnk )
930 0 : call outfld( 'UW_wbrk', wbrk, pcols, lchnk )
931 0 : call outfld( 'UW_lbrk', lbrk, pcols, lchnk )
932 0 : call outfld( 'UW_ricl', ricl, pcols, lchnk )
933 0 : call outfld( 'UW_ghcl', ghcl, pcols, lchnk )
934 0 : call outfld( 'UW_shcl', shcl, pcols, lchnk )
935 0 : call outfld( 'UW_smcl', smcl, pcols, lchnk )
936 :
937 0 : call outfld( 'UW_gh', ghi, pcols, lchnk )
938 0 : call outfld( 'UW_sh', shi, pcols, lchnk )
939 0 : call outfld( 'UW_sm', smi, pcols, lchnk )
940 0 : call outfld( 'UW_ria', rii, pcols, lchnk )
941 0 : call outfld( 'UW_leng', lengi, pcols, lchnk )
942 :
943 0 : call outfld( 'UW_wsed', wsed, pcols, lchnk )
944 :
945 0 : if (.not. unicon_is_on) then
946 0 : deallocate(bprod, ipbl, kpblh, wstarPBL, tkes, went)
947 : end if
948 :
949 0 : end subroutine compute_eddy_diff
950 :
951 : end module eddy_diff_cam
|
