FV3DYCORE  Version1.0.0
fv_sg.F90
Go to the documentation of this file.
1 
2 !***********************************************************************
3 !* GNU Lesser General Public License
4 !*
5 !* This file is part of the FV3 dynamical core.
6 !*
7 !* The FV3 dynamical core is free software: you can redistribute it
8 !* and/or modify it under the terms of the
9 !* GNU Lesser General Public License as published by the
10 !* Free Software Foundation, either version 3 of the License, or
11 !* (at your option) any later version.
12 !*
13 !* The FV3 dynamical core is distributed in the hope that it will be
14 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
15 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16 !* See the GNU General Public License for more details.
17 !*
18 !* You should have received a copy of the GNU Lesser General Public
19 !* License along with the FV3 dynamical core.
20 !* If not, see <http://www.gnu.org/licenses/>.
21 !***********************************************************************
22 
24 
25 ! Modules Included:
26 ! <table>
27 ! <tr>
28 ! <th>Module Name</th>
29 ! <th>Functions Included</th>
30 ! </tr>
31 ! <tr>
32 ! <td>constants_mod</td>
33 ! <td>rdgas, rvgas, cp_air, cp_vapor, hlv, hlf, kappa, grav</td>
34 ! </tr>
35 ! <tr>
36 ! <td>field_manager_mod</td>
37 ! <td>MODEL_ATMOS</td>
38 ! </tr>
39 ! <tr>
40 ! <tr>
41 ! <td>fv_mp_mod</td>
42 ! <td>mp_reduce_min, is_master</td>
43 ! </tr>
44 ! <tr>
45 ! <td>gfdl_cloud_microphys_mod</td>
46 ! <td>wqs1, wqs2, wqsat2_moist</td>
47 ! </tr>
48 ! <tr>
49 ! <td>tracer_manager_mod</td>
50 ! <td>get_tracer_index</td>
51 ! </tr>
52 ! </table>
53 
54 module fv_sg_mod
55 
56 !-----------------------------------------------------------------------
57 ! FV sub-grid mixing
58 !-----------------------------------------------------------------------
59  use constants_mod, only: rdgas, rvgas, cp_air, cp_vapor, hlv, hlf, kappa, grav
60  use tracer_manager_mod, only: get_tracer_index
61  use field_manager_mod, only: model_atmos
62 #ifndef GFS_PHYS
63  use gfdl_cloud_microphys_mod, only: wqs1, wqs2, wqsat2_moist
64 #endif
65  use fv_mp_mod, only: mp_reduce_min, is_master
66 #ifdef MULTI_GASES
67  use multi_gases_mod, only: virq, virqd, virq_qpz, vicpqd_qpz, vicvqd_qpz, vicpqd, vicvqd, num_gas
68 #endif
69 implicit none
70 private
71 
72 public fv_subgrid_z, qsmith, neg_adj3, neg_adj2
73 
74  real, parameter:: esl = 0.621971831
75  real, parameter:: tice = 273.16
76 ! real, parameter:: c_ice = 2106. !< \cite emanuel1994atmospheric table, page 566
77  real, parameter:: c_ice = 1972.
78  real, parameter:: c_liq = 4.1855e+3
79 ! real, parameter:: c_liq = 4218. !< ECMWF-IFS
80  real, parameter:: cv_vap = cp_vapor - rvgas
81  real, parameter:: c_con = c_ice
82 
83 ! real, parameter:: dc_vap = cp_vapor - c_liq ! = -2368.
84  real, parameter:: dc_vap = cv_vap - c_liq
85  real, parameter:: dc_ice = c_liq - c_ice
86 ! Values at 0 Deg C
87  real, parameter:: hlv0 = 2.5e6
88  real, parameter:: hlf0 = 3.3358e5
89 ! real, parameter:: hlv0 = 2.501e6 ! \cite emanuel1994atmospheri Appendix-2
90 ! real, parameter:: hlf0 = 3.337e5 ! \cite emanuel1994atmospheri
91  real, parameter:: t_ice = 273.16
92  real, parameter:: ri_max = 1.
93  real, parameter:: ri_min = 0.25
94  real, parameter:: t1_min = 160.
95  real, parameter:: t2_min = 165.
96  real, parameter:: t2_max = 315.
97  real, parameter:: t3_max = 325.
98  real, parameter:: lv0 = hlv0 - dc_vap*t_ice
99  real, parameter:: li0 = hlf0 - dc_ice*t_ice
100 
101  real, parameter:: zvir = rvgas/rdgas - 1.
102  real, allocatable:: table(:),des(:)
103  real:: lv00, d0_vap
104 
105 contains
106 
107 
108 #ifdef GFS_PHYS
109 
113  subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
114  tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, &
115  hydrostatic, w, delz, u_dt, v_dt, t_dt, k_bot )
116 ! Dry convective adjustment-mixing
117 !-------------------------------------------
118  integer, intent(in):: is, ie, js, je, km, nq, nwat
119  integer, intent(in):: isd, ied, jsd, jed
120  integer, intent(in):: tau
121  real, intent(in):: dt
122  real, intent(in):: pe(is-1:ie+1,km+1,js-1:je+1)
123  real, intent(in):: peln(is :ie, km+1,js :je)
124  real, intent(in):: delp(isd:ied,jsd:jed,km)
125  real, intent(in):: delz(isd:,jsd:,1:)
126  real, intent(in):: pkz(is:ie,js:je,km)
127  logical, intent(in):: hydrostatic
128  integer, intent(in), optional:: k_bot
129 !
130  real, intent(inout):: ua(isd:ied,jsd:jed,km)
131  real, intent(inout):: va(isd:ied,jsd:jed,km)
132  real, intent(inout):: w(isd:,jsd:,1:)
133  real, intent(inout):: ta(isd:ied,jsd:jed,km)
134  real, intent(inout):: qa(isd:ied,jsd:jed,km,nq)
135  real, intent(inout):: u_dt(isd:ied,jsd:jed,km)
136  real, intent(inout):: v_dt(isd:ied,jsd:jed,km)
137  real, intent(inout):: t_dt(is:ie,js:je,km)
138 !---------------------------Local variables-----------------------------
139  real, dimension(is:ie,km):: u0, v0, w0, t0, hd, te, gz, tvm, pm, den
140  real q0(is:ie,km,nq), qcon(is:ie,km)
141  real, dimension(is:ie):: gzh, lcp2, icp2, cvm, cpm, qs
142 #ifdef MULTI_GASES
143  real :: rkx, rdx, rzx, c_air
144 #endif
145  real ri_ref, ri, pt1, pt2, ratio, tv, cv, tmp, q_liq, q_sol
146  real tv1, tv2, g2, h0, mc, fra, rk, rz, rdt, tvd, tv_surf
147  real dh, dq, qsw, dqsdt, tcp3, t_max, t_min
148  integer i, j, k, kk, n, m, iq, km1, im, kbot, l
149  real, parameter:: ustar2 = 1.e-4
150  real:: cv_air, xvir
151  integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
152 
153  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
154  rk = cp_air/rdgas + 1.
155  cv = cp_air - rdgas
156 
157  g2 = 0.5*grav
158 
159  rdt = 1./ dt
160  im = ie-is+1
161 
162  if ( present(k_bot) ) then
163  if ( k_bot < 3 ) return
164  kbot = k_bot
165  else
166  kbot = km
167  endif
168  if ( pe(is,1,js) < 2. ) then
169  t_min = t1_min
170  else
171  t_min = t2_min
172  endif
173 
174  if ( k_bot < min(km,24) ) then
175  t_max = t2_max
176  else
177  t_max = t3_max
178  endif
179 
180  sphum = get_tracer_index(model_atmos, 'sphum')
181  if ( nwat == 0 ) then
182  xvir = 0.
183  rz = 0.
184  else
185  xvir = zvir
186  rz = rvgas - rdgas ! rz = zvir * rdgas
187  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
188  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
189  rainwat = get_tracer_index(model_atmos, 'rainwat')
190  snowwat = get_tracer_index(model_atmos, 'snowwat')
191  graupel = get_tracer_index(model_atmos, 'graupel')
192  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
193  endif
194 
195 !------------------------------------------------------------------------
196 ! The nonhydrostatic pressure changes if there is heating (under constant
197 ! volume and mass is locally conserved).
198 !------------------------------------------------------------------------
199  m = 3
200  fra = dt/real(tau)
201 
202 !$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
203 !$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
204 !$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra, t_max, t_min, &
205 #ifdef MULTI_GASES
206 !$OMP u_dt,rdt,v_dt,xvir,nwat,km) &
207 #else
208 !$OMP u_dt,rdt,v_dt,xvir,nwat) &
209 #endif
210 !$OMP private(kk,lcp2,icp2,tcp3,dh,dq,den,qs,qsw,dqsdt,qcon,q0, &
211 !$OMP t0,u0,v0,w0,h0,pm,gzh,tvm,tmp,cpm,cvm,q_liq,q_sol, &
212 #ifdef MULTI_GASES
213 !$OMP rkx,rdx,rzx,c_air, &
214 #endif
215 !$OMP tv,gz,hd,te,ratio,pt1,pt2,tv1,tv2,ri_ref, ri,mc,km1)
216  do 1000 j=js,je
217 
218  do iq=1, nq
219  do k=1,kbot
220  do i=is,ie
221  q0(i,k,iq) = qa(i,j,k,iq)
222  enddo
223  enddo
224  enddo
225 
226  do k=1,kbot
227  do i=is,ie
228  t0(i,k) = ta(i,j,k)
229 #ifdef MULTI_GASES
230  tvm(i,k) = t0(i,k)*virq(q0(i,k,:))
231 #else
232  tvm(i,k) = t0(i,k)*(1.+xvir*q0(i,k,sphum))
233 #endif
234  u0(i,k) = ua(i,j,k)
235  v0(i,k) = va(i,j,k)
236  pm(i,k) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
237  enddo
238  enddo
239 
240  do i=is,ie
241  gzh(i) = 0.
242  enddo
243 
244  if( hydrostatic ) then
245  do k=kbot, 1,-1
246  do i=is,ie
247  tv = rdgas*tvm(i,k)
248  den(i,k) = pm(i,k)/tv
249  gz(i,k) = gzh(i) + tv*(1.-pe(i,k,j)/pm(i,k))
250  hd(i,k) = cp_air*tvm(i,k)+gz(i,k)+0.5*(u0(i,k)**2+v0(i,k)**2)
251  gzh(i) = gzh(i) + tv*(peln(i,k+1,j)-peln(i,k,j))
252  enddo
253  enddo
254  else
255  do k=kbot, 1, -1
256  if ( nwat == 0 ) then
257  do i=is,ie
258 #ifdef MULTI_GASES
259  cpm(i) = cp_air*vicpqd(q0(i,k,:))
260  cvm(i) = cv_air*vicvqd(q0(i,k,:))
261 #else
262  cpm(i) = cp_air
263  cvm(i) = cv_air
264 #endif
265  enddo
266  elseif ( nwat==1 ) then
267  do i=is,ie
268 #ifdef MULTI_GASES
269  cpm(i) = (1.-q0(i,k,sphum))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor
270  cvm(i) = (1.-q0(i,k,sphum))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap
271 #else
272  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
273  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
274 #endif
275  enddo
276  elseif ( nwat==2 ) then ! GFS
277  do i=is,ie
278 #ifdef MULTI_GASES
279  c_air = cp_air*vicpqd(q0(i,k,:))
280  cpm(i) = c_air + (cp_vapor-c_air)*q0(i,k,sphum)/(1.0-q0(i,k,nwat))
281  c_air = cv_air*vicvqd(q0(i,k,:))
282  cvm(i) = c_air + (cv_vap-c_air)*q0(i,k,sphum)/(1.0-q0(i,k,nwat))
283 #else
284  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
285  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
286 #endif
287  enddo
288  elseif ( nwat==3 ) then
289  do i=is,ie
290  q_liq = q0(i,k,liq_wat)
291  q_sol = q0(i,k,ice_wat)
292 #ifdef MULTI_GASES
293  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
294  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
295 #else
296  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
297  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
298 #endif
299  enddo
300  elseif ( nwat==4 ) then
301  do i=is,ie
302  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
303 #ifdef MULTI_GASES
304  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
305  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq
306 #else
307  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
308  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq
309 #endif
310  enddo
311  elseif ( nwat==5 ) then
312  do i=is,ie
313  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
314  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat)
315 #ifdef MULTI_GASES
316  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
317  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
318 #else
319  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
320  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
321 #endif
322  enddo
323  else
324  do i=is,ie
325  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
326  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
327 #ifdef MULTI_GASES
328  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
329  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
330 #else
331  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
332  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
333 #endif
334  enddo
335  endif
336 
337  do i=is,ie
338  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
339  w0(i,k) = w(i,j,k)
340  gz(i,k) = gzh(i) - g2*delz(i,j,k)
341  tmp = gz(i,k) + 0.5*(u0(i,k)**2+v0(i,k)**2+w0(i,k)**2)
342  hd(i,k) = cpm(i)*t0(i,k) + tmp
343  te(i,k) = cvm(i)*t0(i,k) + tmp
344  gzh(i) = gzh(i) - grav*delz(i,j,k)
345  enddo
346  enddo
347  endif
348 
349  do n=1,m
350 
351  if ( m==3 ) then
352  if ( n==1) ratio = 0.25
353  if ( n==2) ratio = 0.5
354  if ( n==3) ratio = 0.999
355  else
356  ratio = real(n)/real(m)
357  endif
358 
359  do i=is,ie
360  gzh(i) = 0.
361  enddo
362 
363 ! Compute total condensate
364  if ( nwat<2 ) then
365  do k=1,kbot
366  do i=is,ie
367  qcon(i,k) = 0.
368  enddo
369  enddo
370  elseif ( nwat==2 ) then ! GFS_2015
371  do k=1,kbot
372  do i=is,ie
373  qcon(i,k) = q0(i,k,liq_wat)
374  enddo
375  enddo
376  elseif ( nwat==3 ) then
377  do k=1,kbot
378  do i=is,ie
379  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,ice_wat)
380  enddo
381  enddo
382  elseif ( nwat==4 ) then
383  do k=1,kbot
384  do i=is,ie
385  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat)
386  enddo
387  enddo
388  elseif ( nwat==5 ) then
389  do k=1,kbot
390  do i=is,ie
391  qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)
392  enddo
393  enddo
394  else
395  do k=1,kbot
396  do i=is,ie
397  qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)
398  enddo
399  enddo
400  endif
401 
402  do k=kbot, 2, -1
403  km1 = k-1
404  do i=is,ie
405 ! Richardson number = g*delz * del_theta/theta / (del_u**2 + del_v**2)
406 ! Use exact form for "density temperature"
407 #ifdef MULTI_GASES
408  tv1 = t0(i,km1)*virq_qpz(q0(i,km1,:),qcon(i,km1))
409  tv2 = t0(i,k )*virq_qpz(q0(i,k, :),qcon(i,k ))
410 #else
411  tv1 = t0(i,km1)*(1.+xvir*q0(i,km1,sphum)-qcon(i,km1))
412  tv2 = t0(i,k )*(1.+xvir*q0(i,k ,sphum)-qcon(i,k))
413 #endif
414  pt1 = tv1 / pkz(i,j,km1)
415  pt2 = tv2 / pkz(i,j,k )
416 !
417  ri = (gz(i,km1)-gz(i,k))*(pt1-pt2)/( 0.5*(pt1+pt2)* &
418  ((u0(i,km1)-u0(i,k))**2+(v0(i,km1)-v0(i,k))**2+ustar2) )
419  if ( tv1>t_max .and. tv1>tv2 ) then
420 ! top layer unphysically warm
421  ri = 0.
422  elseif ( tv2<t_min ) then
423  ri = min(ri, 0.1)
424  endif
425 ! Adjustment for K-H instability:
426 ! Compute equivalent mass flux: mc
427 ! Add moist 2-dz instability consideration:
428 !!! ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(500.e2,pm(i,k))/250.e2 )
429  ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(400.e2,pm(i,k))/200.e2 )
430 ! Enhancing mixing at the model top
431  if ( k==2 ) then
432  ri_ref = 4.*ri_ref
433  elseif ( k==3 ) then
434  ri_ref = 2.*ri_ref
435  elseif ( k==4 ) then
436  ri_ref = 1.5*ri_ref
437  endif
438 
439  if ( ri < ri_ref ) then
440  mc = ratio*delp(i,j,km1)*delp(i,j,k)/(delp(i,j,km1)+delp(i,j,k))*(1.-max(0.0,ri/ri_ref))**2
441  do iq=1,nq
442  h0 = mc*(q0(i,k,iq)-q0(i,km1,iq))
443  q0(i,km1,iq) = q0(i,km1,iq) + h0/delp(i,j,km1)
444  q0(i,k ,iq) = q0(i,k ,iq) - h0/delp(i,j,k )
445  enddo
446 ! Recompute qcon
447  if ( nwat<2 ) then
448  qcon(i,km1) = 0.
449  elseif ( nwat==2 ) then ! GFS_2015
450  qcon(i,km1) = q0(i,km1,liq_wat)
451  elseif ( nwat==3 ) then ! AM3/AM4
452  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat)
453  elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice
454  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat)
455  elseif ( nwat==5 ) then ! K_warm_rain scheme with fake ice
456  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
457  q0(i,km1,snowwat) + q0(i,km1,rainwat)
458  else
459  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
460  q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
461  endif
462 ! u:
463  h0 = mc*(u0(i,k)-u0(i,k-1))
464  u0(i,k-1) = u0(i,k-1) + h0/delp(i,j,k-1)
465  u0(i,k ) = u0(i,k ) - h0/delp(i,j,k )
466 ! v:
467  h0 = mc*(v0(i,k)-v0(i,k-1))
468  v0(i,k-1) = v0(i,k-1) + h0/delp(i,j,k-1)
469  v0(i,k ) = v0(i,k ) - h0/delp(i,j,k )
470 
471  if ( hydrostatic ) then
472 ! Static energy
473  h0 = mc*(hd(i,k)-hd(i,k-1))
474  hd(i,k-1) = hd(i,k-1) + h0/delp(i,j,k-1)
475  hd(i,k ) = hd(i,k ) - h0/delp(i,j,k )
476  else
477 ! Total energy
478  h0 = mc*(hd(i,k)-hd(i,k-1))
479  te(i,k-1) = te(i,k-1) + h0/delp(i,j,k-1)
480  te(i,k ) = te(i,k ) - h0/delp(i,j,k )
481 ! w:
482  h0 = mc*(w0(i,k)-w0(i,k-1))
483  w0(i,k-1) = w0(i,k-1) + h0/delp(i,j,k-1)
484  w0(i,k ) = w0(i,k ) - h0/delp(i,j,k )
485  endif
486  endif
487  enddo
488 
489 !--------------
490 ! Retrive Temp:
491 !--------------
492  if ( hydrostatic ) then
493  kk = k
494  do i=is,ie
495 #ifdef MULTI_GASES
496  rkx = cp_air/rdgas * (vicpqd(q0(i,kk,:))/virqd(q0(i,kk,:))) + 1
497  rdx = rdgas*virqd(q0(i,kk,:))
498  rzx = rvgas - rdx
499  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
500  / ( rkx - pe(i,kk,j)/pm(i,kk) )
501  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
502  rdx = rdx + rzx*q0(i,kk,sphum)/(1.0-sum(q0(i,kk,sphum+1:sphum+nwat-1)))
503  t0(i,kk) = t0(i,kk) / rdx
504 #else
505  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
506  / ( rk - pe(i,kk,j)/pm(i,kk) )
507  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
508  t0(i,kk) = t0(i,kk) / ( rdgas + rz*q0(i,kk,sphum) )
509 #endif
510  enddo
511  kk = k-1
512  do i=is,ie
513 #ifdef MULTI_GASES
514  rkx = cp_air/rdgas * (vicpqd(q0(i,kk,:))/virqd(q0(i,kk,:))) + 1
515  rdx = rdgas*virqd(q0(i,kk,:))
516  rzx = rvgas - rdx
517  rdx = rdx + rzx*q0(i,kk,sphum)/(1.0-sum(q0(i,kk,sphum+1:sphum+nwat-1)))
518  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
519  / ((rkx-pe(i,kk,j)/pm(i,kk))*rdx)
520 #else
521  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
522  / ((rk-pe(i,kk,j)/pm(i,kk))*(rdgas+rz*q0(i,kk,sphum)))
523 #endif
524  enddo
525  else
526 ! Non-hydrostatic under constant volume heating/cooling
527  do kk=k-1,k
528  if ( nwat == 0 ) then
529  do i=is,ie
530 #ifdef MULTI_GASES
531  cpm(i) = cp_air*vicpqd(q0(i,kk,:))
532  cvm(i) = cv_air*vicvqd(q0(i,kk,:))
533 #else
534  cpm(i) = cp_air
535  cvm(i) = cv_air
536 #endif
537  enddo
538  elseif ( nwat == 1 ) then
539  do i=is,ie
540 #ifdef MULTI_GASES
541  cpm(i) = (1.-q0(i,kk,sphum))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor
542  cvm(i) = (1.-q0(i,kk,sphum))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap
543 #else
544  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
545  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
546 #endif
547  enddo
548  elseif ( nwat == 2 ) then
549  do i=is,ie
550 #ifdef MULTI_GASES
551  q_liq = q0(i,kk,liq_wat)
552  c_air = cp_air*vicpqd(q0(i,kk,:))
553  cpm(i) = c_air + (cp_vapor-c_air)*q0(i,kk,sphum)/(1.0-q_liq)
554  c_air = cv_air*vicvqd(q0(i,kk,:))
555  cvm(i) = c_air + (cv_vap-c_air)*q0(i,kk,sphum)/(1.0-q_liq)
556 #else
557  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
558  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
559 #endif
560  enddo
561  elseif ( nwat == 3 ) then
562  do i=is,ie
563  q_liq = q0(i,kk,liq_wat)
564  q_sol = q0(i,kk,ice_wat)
565 #ifdef MULTI_GASES
566  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
567  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
568 #else
569  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
570  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
571 #endif
572  enddo
573  elseif ( nwat == 4 ) then
574  do i=is,ie
575  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
576 #ifdef MULTI_GASES
577  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
578  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
579 #else
580  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
581  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
582 #endif
583  enddo
584  elseif ( nwat == 5 ) then
585  do i=is,ie
586  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
587  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat)
588 #ifdef MULTI_GASES
589  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
590  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
591 #else
592  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
593  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
594 #endif
595  enddo
596  else
597  do i=is,ie
598  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
599  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
600 #ifdef MULTI_GASES
601  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
602  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
603 #else
604  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
605  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
606 #endif
607  enddo
608  endif
609 
610  do i=is,ie
611  tv = gz(i,kk) + 0.5*(u0(i,kk)**2+v0(i,kk)**2+w0(i,kk)**2)
612  t0(i,kk) = (te(i,kk)- tv) / cvm(i)
613  hd(i,kk) = cpm(i)*t0(i,kk) + tv
614  enddo
615  enddo
616  endif
617  enddo ! k-loop
618  enddo ! n-loop
619 
620 !--------------------
621  if ( fra < 1. ) then
622  do k=1, kbot
623  do i=is,ie
624  t0(i,k) = ta(i,j,k) + (t0(i,k) - ta(i,j,k))*fra
625  u0(i,k) = ua(i,j,k) + (u0(i,k) - ua(i,j,k))*fra
626  v0(i,k) = va(i,j,k) + (v0(i,k) - va(i,j,k))*fra
627  enddo
628  enddo
629 
630  if ( .not. hydrostatic ) then
631  do k=1,kbot
632  do i=is,ie
633  w0(i,k) = w(i,j,k) + (w0(i,k) - w(i,j,k))*fra
634  enddo
635  enddo
636  endif
637 
638  do iq=1,nq
639  do k=1,kbot
640  do i=is,ie
641  q0(i,k,iq) = qa(i,j,k,iq) + (q0(i,k,iq) - qa(i,j,k,iq))*fra
642  enddo
643  enddo
644  enddo
645  endif
646 
647  do k=1,kbot
648  do i=is,ie
649  u_dt(i,j,k) = rdt*(u0(i,k) - ua(i,j,k))
650  v_dt(i,j,k) = rdt*(v0(i,k) - va(i,j,k))
651  ta(i,j,k) = t0(i,k) ! *** temperature updated ***
652 #ifdef GFS_PHYS
653  ua(i,j,k) = u0(i,k)
654  va(i,j,k) = v0(i,k)
655 #endif
656  enddo
657  do iq=1,nq
658  do i=is,ie
659  qa(i,j,k,iq) = q0(i,k,iq)
660  enddo
661  enddo
662  enddo
663 
664  if ( .not. hydrostatic ) then
665  do k=1,kbot
666  do i=is,ie
667  w(i,j,k) = w0(i,k) ! w updated
668  enddo
669  enddo
670  endif
671 
672 1000 continue
673 
674  end subroutine fv_subgrid_z
675 
676 #else
677  subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
678  tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, &
679  hydrostatic, w, delz, u_dt, v_dt, t_dt, q_dt, k_bot )
680 ! Dry convective adjustment-mixing
681 !-------------------------------------------
682  integer, intent(in):: is, ie, js, je, km, nq, nwat
683  integer, intent(in):: isd, ied, jsd, jed
684  integer, intent(in):: tau
685  real, intent(in):: dt
686  real, intent(in):: pe(is-1:ie+1,km+1,js-1:je+1)
687  real, intent(in):: peln(is :ie, km+1,js :je)
688  real, intent(in):: delp(isd:ied,jsd:jed,km)
689  real, intent(in):: delz(isd:,jsd:,1:)
690  real, intent(in):: pkz(is:ie,js:je,km)
691  logical, intent(in):: hydrostatic
692  integer, intent(in), optional:: k_bot
693 !
694  real, intent(inout):: ua(isd:ied,jsd:jed,km)
695  real, intent(inout):: va(isd:ied,jsd:jed,km)
696  real, intent(inout):: w(isd:,jsd:,1:)
697  real, intent(inout):: ta(isd:ied,jsd:jed,km)
698  real, intent(inout):: qa(isd:ied,jsd:jed,km,nq)
699  real, intent(inout):: u_dt(isd:ied,jsd:jed,km)
700  real, intent(inout):: v_dt(isd:ied,jsd:jed,km)
701  real, intent(inout):: t_dt(is:ie,js:je,km)
702  real, intent(inout):: q_dt(is:ie,js:je,km,nq)
703 !---------------------------Local variables-----------------------------
704  real, dimension(is:ie,km):: u0, v0, w0, t0, hd, te, gz, tvm, pm, den
705  real q0(is:ie,km,nq), qcon(is:ie,km)
706  real, dimension(is:ie):: gzh, lcp2, icp2, cvm, cpm, qs
707 #ifdef MULTI_GASES
708  real :: rkx, rdx, rzx, c_air
709 #endif
710  real ri_ref, ri, pt1, pt2, ratio, tv, cv, tmp, q_liq, q_sol
711  real tv1, tv2, g2, h0, mc, fra, rk, rz, rdt, tvd, tv_surf
712  real dh, dq, qsw, dqsdt, tcp3
713  integer i, j, k, kk, n, m, iq, km1, im, kbot
714  real, parameter:: ustar2 = 1.e-4
715  real:: cv_air, xvir
716  integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
717 
718  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
719  rk = cp_air/rdgas + 1.
720  cv = cp_air - rdgas
721 
722  g2 = 0.5*grav
723 
724  rdt = 1./ dt
725  im = ie-is+1
726 
727  if ( present(k_bot) ) then
728  if ( k_bot < 3 ) return
729  kbot = k_bot
730  else
731  kbot = km
732  endif
733 
734  sphum = get_tracer_index(model_atmos, 'sphum')
735  if ( nwat == 0 ) then
736  xvir = 0.
737  rz = 0.
738  else
739  xvir = zvir
740  rz = rvgas - rdgas ! rz = zvir * rdgas
741  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
742  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
743  rainwat = get_tracer_index(model_atmos, 'rainwat')
744  snowwat = get_tracer_index(model_atmos, 'snowwat')
745  graupel = get_tracer_index(model_atmos, 'graupel')
746  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
747  endif
748 
749 !------------------------------------------------------------------------
750 ! The nonhydrostatic pressure changes if there is heating (under constant
751 ! volume and mass is locally conserved).
752 !------------------------------------------------------------------------
753  m = 3
754  fra = dt/real(tau)
755 
756 !$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
757 !$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
758 !$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra,cld_amt, &
759 !$OMP u_dt,rdt,v_dt,xvir,nwat) &
760 !$OMP private(kk,lcp2,icp2,tcp3,dh,dq,den,qs,qsw,dqsdt,qcon,q0, &
761 !$OMP t0,u0,v0,w0,h0,pm,gzh,tvm,tmp,cpm,cvm, q_liq,q_sol,&
762 #ifdef MULTI_GASES
763 !$OMP rkx,rdx,rzx,c_air &
764 #endif
765 !$OMP tv,gz,hd,te,ratio,pt1,pt2,tv1,tv2,ri_ref, ri,mc,km1)
766  do 1000 j=js,je
767 
768  do iq=1, nq
769  do k=1,kbot
770  do i=is,ie
771  q0(i,k,iq) = qa(i,j,k,iq)
772  enddo
773  enddo
774  enddo
775 
776  do k=1,kbot
777  do i=is,ie
778  t0(i,k) = ta(i,j,k)
779 #ifdef MULTI_GASES
780  tvm(i,k) = t0(i,k)*virq(q0(i,k,:))
781 #else
782  tvm(i,k) = t0(i,k)*(1.+xvir*q0(i,k,sphum))
783 #endif
784  u0(i,k) = ua(i,j,k)
785  v0(i,k) = va(i,j,k)
786  pm(i,k) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
787  enddo
788  enddo
789 
790  do i=is,ie
791  gzh(i) = 0.
792  enddo
793 
794  if( hydrostatic ) then
795  do k=kbot, 1,-1
796  do i=is,ie
797  tv = rdgas*tvm(i,k)
798  den(i,k) = pm(i,k)/tv
799  gz(i,k) = gzh(i) + tv*(1.-pe(i,k,j)/pm(i,k))
800  hd(i,k) = cp_air*tvm(i,k)+gz(i,k)+0.5*(u0(i,k)**2+v0(i,k)**2)
801  gzh(i) = gzh(i) + tv*(peln(i,k+1,j)-peln(i,k,j))
802  enddo
803  enddo
804  else
805  do k=kbot, 1, -1
806  if ( nwat == 0 ) then
807  do i=is,ie
808 #ifdef MULTI_GASES
809  cpm(i) = cp_air*vicpqd(q0(i,k,:))
810  cvm(i) = cv_air*vicvqd(q0(i,k,:))
811 #else
812  cpm(i) = cp_air
813  cvm(i) = cv_air
814 #endif
815  enddo
816  elseif ( nwat==1 ) then
817  do i=is,ie
818 #ifdef MULTI_GASES
819  cpm(i) = (1.-q0(i,k,sphum))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor
820  cvm(i) = (1.-q0(i,k,sphum))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap
821 #else
822  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
823  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
824 #endif
825  enddo
826  elseif ( nwat==2 ) then ! GFS
827  do i=is,ie
828  q_sol = q0(i,k,liq_wat)
829 #ifdef MULTI_GASES
830  c_air = cp_air*vicpqd(q0(i,k,:))
831  cpm(i) = c_air + (cp_vapor-c_air)*q0(i,k,sphum)/(1.0-q0(i,k,nwat))
832  c_air = cv_air*vicvqd(q0(i,k,:))
833  cvm(i) = c_air + (cv_vap-c_air)*q0(i,k,sphum)/(1.0-q0(i,k,nwat))
834 #else
835  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
836  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
837 #endif
838  enddo
839  elseif ( nwat==3 ) then
840  do i=is,ie
841  q_liq = q0(i,k,liq_wat)
842  q_sol = q0(i,k,ice_wat)
843 #ifdef MULTI_GASES
844  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
845  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
846 #else
847  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
848  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
849 #endif
850  enddo
851  elseif ( nwat==4 ) then
852  do i=is,ie
853  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
854 #ifdef MULTI_GASES
855  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
856  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq
857 #else
858  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
859  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq
860 #endif
861  enddo
862  elseif ( nwat==5 ) then
863  do i=is,ie
864  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
865  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat)
866 #ifdef MULTI_GASES
867  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
868  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
869 #else
870  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
871  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
872 #endif
873  enddo
874  else
875  do i=is,ie
876  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
877  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
878 #ifdef MULTI_GASES
879  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
880  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
881 #else
882  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
883  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
884 #endif
885  enddo
886  endif
887 
888  do i=is,ie
889  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
890  w0(i,k) = w(i,j,k)
891  gz(i,k) = gzh(i) - g2*delz(i,j,k)
892  tmp = gz(i,k) + 0.5*(u0(i,k)**2+v0(i,k)**2+w0(i,k)**2)
893  hd(i,k) = cpm(i)*t0(i,k) + tmp
894  te(i,k) = cvm(i)*t0(i,k) + tmp
895  gzh(i) = gzh(i) - grav*delz(i,j,k)
896  enddo
897  enddo
898  endif
899 
900  do n=1,m
901 
902  if ( m==3 ) then
903  if ( n==1) ratio = 0.25
904  if ( n==2) ratio = 0.5
905  if ( n==3) ratio = 0.999
906  else
907  ratio = real(n)/real(m)
908  endif
909 
910  do i=is,ie
911  gzh(i) = 0.
912  enddo
913 
914 ! Compute total condensate
915  if ( nwat<2 ) then
916  do k=1,kbot
917  do i=is,ie
918  qcon(i,k) = 0.
919  enddo
920  enddo
921  elseif ( nwat==2 ) then ! GFS_2015
922  do k=1,kbot
923  do i=is,ie
924  qcon(i,k) = q0(i,k,liq_wat)
925  enddo
926  enddo
927  elseif ( nwat==3 ) then
928  do k=1,kbot
929  do i=is,ie
930  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,ice_wat)
931  enddo
932  enddo
933  elseif ( nwat==4 ) then
934  do k=1,kbot
935  do i=is,ie
936  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat)
937  enddo
938  enddo
939  elseif ( nwat==5 ) then
940  do k=1,kbot
941  do i=is,ie
942  qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)
943  enddo
944  else
945  do k=1,kbot
946  do i=is,ie
947  qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)
948  enddo
949  enddo
950  endif
951 
952  do k=kbot, 2, -1
953  km1 = k-1
954 #ifdef TEST_MQ
955 #ifdef MULTI_GASES
956  if(nwat>0) call qsmith(im, km, im, 1, t0(is,km1), pm(is,km1), q0(is,km1,sphum), qs)
957 #else
958  if(nwat>0) call qsmith(im, 1, 1, t0(is,km1), pm(is,km1), q0(is,km1,sphum), qs)
959 #endif
960 #endif
961  do i=is,ie
962 ! Richardson number = g*delz * del_theta/theta / (del_u**2 + del_v**2)
963 ! Use exact form for "density temperature"
964 #ifdef MULTI_GASES
965  tv1 = t0(i,km1)*virq_qpz(q0(i,km1,:),qcon(i,km1))
966  tv2 = t0(i,k )*virq_qpz(q0(i,k ,:),qcon(i,k ))
967 #else
968  tv1 = t0(i,km1)*(1.+xvir*q0(i,km1,sphum)-qcon(i,km1))
969  tv2 = t0(i,k )*(1.+xvir*q0(i,k ,sphum)-qcon(i,k))
970 #endif
971  pt1 = tv1 / pkz(i,j,km1)
972  pt2 = tv2 / pkz(i,j,k )
973  ri = (gz(i,km1)-gz(i,k))*(pt1-pt2)/( 0.5*(pt1+pt2)* &
974  ((u0(i,km1)-u0(i,k))**2+(v0(i,km1)-v0(i,k))**2+ustar2) )
975 
976 ! Adjustment for K-H instability:
977 ! Compute equivalent mass flux: mc
978 ! Add moist 2-dz instability consideration:
979  ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(500.e2,pm(i,k))/250.e2 )
980 #ifdef TEST_MQ
981  if ( nwat > 0 ) then
982  h0 = hd(i,k)-hd(i,km1) + hlv0*(q0(i,k,sphum)-qs(i))
983  if ( h0 > 0. ) ri_ref = 5.*ri_ref
984  endif
985 #endif
986  if ( ri < ri_ref ) then
987  mc = ratio*delp(i,j,km1)*delp(i,j,k)/(delp(i,j,km1)+delp(i,j,k))*(1.-max(0.0,ri/ri_ref))**2
988  do iq=1,nq
989  h0 = mc*(q0(i,k,iq)-q0(i,km1,iq))
990  q0(i,km1,iq) = q0(i,km1,iq) + h0/delp(i,j,km1)
991  q0(i,k ,iq) = q0(i,k ,iq) - h0/delp(i,j,k )
992  enddo
993 ! Recompute qcon
994  if ( nwat<2 ) then
995  qcon(i,km1) = 0.
996  elseif ( nwat==2 ) then ! GFS_2015
997  qcon(i,km1) = q0(i,km1,liq_wat)
998  elseif ( nwat==3 ) then ! AM3/AM4
999  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat)
1000  elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice
1001  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat)
1002  elseif ( nwat==5 ) then
1003  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
1004  q0(i,km1,snowwat) + q0(i,km1,rainwat)
1005  else
1006  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
1007  q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
1008  endif
1009 ! u:
1010  h0 = mc*(u0(i,k)-u0(i,k-1))
1011  u0(i,k-1) = u0(i,k-1) + h0/delp(i,j,k-1)
1012  u0(i,k ) = u0(i,k ) - h0/delp(i,j,k )
1013 ! v:
1014  h0 = mc*(v0(i,k)-v0(i,k-1))
1015  v0(i,k-1) = v0(i,k-1) + h0/delp(i,j,k-1)
1016  v0(i,k ) = v0(i,k ) - h0/delp(i,j,k )
1017 
1018  if ( hydrostatic ) then
1019 ! Static energy
1020  h0 = mc*(hd(i,k)-hd(i,k-1))
1021  hd(i,k-1) = hd(i,k-1) + h0/delp(i,j,k-1)
1022  hd(i,k ) = hd(i,k ) - h0/delp(i,j,k )
1023  else
1024 ! Total energy
1025  h0 = mc*(hd(i,k)-hd(i,k-1))
1026  te(i,k-1) = te(i,k-1) + h0/delp(i,j,k-1)
1027  te(i,k ) = te(i,k ) - h0/delp(i,j,k )
1028 ! w:
1029  h0 = mc*(w0(i,k)-w0(i,k-1))
1030  w0(i,k-1) = w0(i,k-1) + h0/delp(i,j,k-1)
1031  w0(i,k ) = w0(i,k ) - h0/delp(i,j,k )
1032  endif
1033  endif
1034  enddo
1035 
1036 !--------------
1037 ! Retrive Temp:
1038 !--------------
1039  if ( hydrostatic ) then
1040  kk = k
1041  do i=is,ie
1042 #ifdef MULTI_GASES
1043  rkx = cp_air/rdgas * (vicpqd(q0(i,kk,:))/virqd(q0(i,kk,:))) + 1
1044  rdx = rdgas*virqd(q0(i,kk,:))
1045  rzx = rvgas - rdx
1046  rdx = rdx + rzx*q0(i,kk,sphum)/(1.0-sum(q0(i,kk,sphum+1:sphum+nwat-1)))
1047  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
1048  / ( rkx - pe(i,kk,j)/pm(i,kk) )
1049  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
1050  t0(i,kk) = t0(i,kk) / ( rdx + rzx*q0(i,kk,sphum) )
1051 #else
1052  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
1053  / ( rk - pe(i,kk,j)/pm(i,kk) )
1054  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
1055  t0(i,kk) = t0(i,kk) / ( rdgas + rz*q0(i,kk,sphum) )
1056 #endif
1057  enddo
1058  kk = k-1
1059  do i=is,ie
1060 #ifdef MULTI_GASES
1061  rkx = cp_air/rdgas * (vicpqd(q0(i,kk,:))/virqd(q0(i,kk,:))) + 1
1062  rdx = rdgas*virqd(i)
1063  rzx = rvgas - rdx
1064  rdx = rdx + rzx*q0(i,kk,sphum)/(1.0-sum(q0(i,kk,sphum+1:sphum+nwat-1)))
1065  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
1066  / ((rkx-pe(i,kk,j)/pm(i,kk))*rdx)
1067 #else
1068  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
1069  / ((rk-pe(i,kk,j)/pm(i,kk))*(rdgas+rz*q0(i,kk,sphum)))
1070 #endif
1071  enddo
1072  else
1073 ! Non-hydrostatic under constant volume heating/cooling
1074  do kk=k-1,k
1075  if ( nwat == 0 ) then
1076  do i=is,ie
1077 #ifdef MULTI_GASES
1078  cpm(i) = cp_air*vicpqd(q0(i,kk,:))
1079  cvm(i) = cv_air*vicvqd(q0(i,kk,:))
1080 #else
1081  cpm(i) = cp_air
1082  cvm(i) = cv_air
1083 #endif
1084  enddo
1085  elseif ( nwat == 1 ) then
1086  do i=is,ie
1087 #ifdef MULTI_GASES
1088  cpm(i) = (1.-q0(i,kk,sphum))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor
1089  cvm(i) = (1.-q0(i,kk,sphum))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap
1090 #else
1091  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
1092  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
1093 #endif
1094  enddo
1095  elseif ( nwat == 2 ) then
1096  do i=is,ie
1097 #ifdef MULTI_GASES
1098  c_air = cp_air*vicpqd(q0(i,kk,:))
1099  cpm(i) = c_air+(cp_vapor-c_air)*q0(i,kk,sphum)/(1.0-q0(i,kk,sphum+nwat-1))
1100  c_air = cv_air*vicvqd(q0(i,kk,:))
1101  cvm(i) = c_air+(cv_vap -c_air)*q0(i,kk,sphum)/(1.0-q0(i,kk,sphum+nwat-1))
1102 #else
1103  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
1104  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
1105 #endif
1106  enddo
1107  elseif ( nwat == 3 ) then
1108  do i=is,ie
1109  q_liq = q0(i,kk,liq_wat)
1110  q_sol = q0(i,kk,ice_wat)
1111 #ifdef MULTI_GASES
1112  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1113  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1114 #else
1115  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1116  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1117 #endif
1118  enddo
1119  elseif ( nwat == 4 ) then
1120  do i=is,ie
1121  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
1122 #ifdef MULTI_GASES
1123  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
1124  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
1125 #else
1126  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
1127  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
1128 #endif
1129  enddo
1130  elseif ( nwat == 5 ) then
1131  do i=is,ie
1132  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
1133  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat)
1134 #ifdef MULTI_GASES
1135  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1136  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1137 #else
1138  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1139  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1140 #endif
1141  enddo
1142  else
1143  do i=is,ie
1144  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
1145  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
1146 #ifdef MULTI_GASES
1147  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1148  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1149 #else
1150  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1151  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1152 #endif
1153  enddo
1154  endif
1155 
1156 
1157  do i=is,ie
1158  tv = gz(i,kk) + 0.5*(u0(i,kk)**2+v0(i,kk)**2+w0(i,kk)**2)
1159  t0(i,kk) = (te(i,kk)- tv) / cvm(i)
1160  hd(i,kk) = cpm(i)*t0(i,kk) + tv
1161  enddo
1162  enddo
1163  endif
1164  enddo ! k-loop
1165  enddo ! n-loop
1166 
1167 !--------------------
1168  if ( fra < 1. ) then
1169  do k=1, kbot
1170  do i=is,ie
1171  t0(i,k) = ta(i,j,k) + (t0(i,k) - ta(i,j,k))*fra
1172  u0(i,k) = ua(i,j,k) + (u0(i,k) - ua(i,j,k))*fra
1173  v0(i,k) = va(i,j,k) + (v0(i,k) - va(i,j,k))*fra
1174  enddo
1175  enddo
1176 
1177  if ( .not. hydrostatic ) then
1178  do k=1,kbot
1179  do i=is,ie
1180  w0(i,k) = w(i,j,k) + (w0(i,k) - w(i,j,k))*fra
1181  enddo
1182  enddo
1183  endif
1184 
1185  do iq=1,nq
1186  do k=1,kbot
1187  do i=is,ie
1188  q0(i,k,iq) = qa(i,j,k,iq) + (q0(i,k,iq) - qa(i,j,k,iq))*fra
1189  enddo
1190  enddo
1191  enddo
1192  endif
1193 
1194 !----------------------
1195 ! Saturation adjustment
1196 !----------------------
1197 #ifndef GFS_PHYS
1198  if ( nwat > 5 ) then
1199  do k=1, kbot
1200  if ( hydrostatic ) then
1201  do i=is, ie
1202 ! Compute pressure hydrostatically
1203 #ifdef MULTI_GASES
1204  den(i,k) = pm(i,k)/(rdgas*t0(i,k)*virq(q0(i,k,:)))
1205 #else
1206  den(i,k) = pm(i,k)/(rdgas*t0(i,k)*(1.+xvir*q0(i,k,sphum)))
1207 #endif
1208  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
1209  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
1210 #ifdef MULTI_GASES
1211  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1212 #else
1213  cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1214 #endif
1215  lcp2(i) = hlv / cpm(i)
1216  icp2(i) = hlf / cpm(i)
1217  enddo
1218  else
1219  do i=is, ie
1220  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
1221  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
1222  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
1223 #ifdef MULTI_GASES
1224  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1225 #else
1226  cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
1227 #endif
1228  lcp2(i) = (lv0+dc_vap*t0(i,k)) / cvm(i)
1229  icp2(i) = (li0+dc_ice*t0(i,k)) / cvm(i)
1230  enddo
1231  endif
1232 
1233 ! Prevent super saturation over water:
1234  do i=is, ie
1235  qsw = wqs2(t0(i,k), den(i,k), dqsdt)
1236  dq = q0(i,k,sphum) - qsw
1237  if ( dq > 0. ) then ! remove super-saturation
1238  tcp3 = lcp2(i) + icp2(i)*min(1., dim(tice,t0(i,k))/40.)
1239  tmp = dq/(1.+tcp3*dqsdt)
1240  t0(i,k) = t0(i,k) + tmp*lcp2(i)
1241  q0(i,k, sphum) = q0(i,k, sphum) - tmp
1242  q0(i,k,liq_wat) = q0(i,k,liq_wat) + tmp
1243 ! Grid box mean is saturated; 50% or higher cloud cover
1244  if (cld_amt .gt. 0) then
1245  qa(i,j,k,cld_amt) = max(0.5, min(1., qa(i,j,k,cld_amt)+25.*tmp/qsw))
1246  end if
1247  endif
1248 ! Freezing
1249  tmp = tice-40. - t0(i,k)
1250  if( tmp>0.0 .and. q0(i,k,liq_wat)>0. ) then
1251  dh = min( q0(i,k,liq_wat), q0(i,k,liq_wat)*tmp*0.125, tmp/icp2(i) )
1252  q0(i,k,liq_wat) = q0(i,k,liq_wat) - dh
1253  q0(i,k,ice_wat) = q0(i,k,ice_wat) + dh
1254  t0(i,k) = t0(i,k) + dh*icp2(i)
1255  endif
1256  enddo
1257  enddo
1258  endif
1259 #endif
1260 
1261  do k=1,kbot
1262  do i=is,ie
1263  u_dt(i,j,k) = rdt*(u0(i,k) - ua(i,j,k))
1264  v_dt(i,j,k) = rdt*(v0(i,k) - va(i,j,k))
1265  ta(i,j,k) = t0(i,k) ! *** temperature updated ***
1266 #ifdef GFS_PHYS
1267  ua(i,j,k) = u0(i,k)
1268  va(i,j,k) = v0(i,k)
1269 #endif
1270  enddo
1271  do iq=1,nq
1272  if (iq .ne. cld_amt ) then
1273  do i=is,ie
1274  qa(i,j,k,iq) = q0(i,k,iq)
1275  enddo
1276  endif
1277  enddo
1278  enddo
1279 
1280  if ( .not. hydrostatic ) then
1281  do k=1,kbot
1282  do i=is,ie
1283  w(i,j,k) = w0(i,k) ! w updated
1284  enddo
1285  enddo
1286  endif
1287 
1288 1000 continue
1289 
1290 
1291  end subroutine fv_subgrid_z
1292 #endif
1293 
1294 
1295  subroutine qsmith_init
1296  integer, parameter:: length=2621
1297  integer i
1298 
1299  if( .not. allocated(table) ) then
1300 ! Generate es table (dT = 0.1 deg. C)
1301 
1302  allocate ( table(length) )
1303  allocate ( des(length) )
1304 
1305  call qs_table(length, table)
1306 
1307  do i=1,length-1
1308  des(i) = table(i+1) - table(i)
1309  enddo
1310  des(length) = des(length-1)
1311  endif
1312 
1313  end subroutine qsmith_init
1314 
1315 
1316 #ifdef MULTI_GASES
1317  subroutine qsmith(im, km, imx, kmx, t, p, q, qs, dqdt)
1318 #else
1319  subroutine qsmith(im, km, k1, t, p, q, qs, dqdt)
1320 #endif
1321 ! input T in deg K; p (Pa)
1322  real, intent(in),dimension(im,km):: t, p
1323 #ifdef MULTI_GASES
1324  real, intent(in),dimension(im,km,*):: q
1325  integer, intent(in):: im, km, imx, kmx
1326 #else
1327  integer, intent(in):: im, km, k1
1328  real, intent(in),dimension(im,km):: q
1329 #endif
1330  real, intent(out),dimension(im,km):: qs
1331  real, intent(out), optional:: dqdt(im,km)
1332 ! Local:
1333  real es(im,km)
1334  real ap1, eps10
1335  real tmin
1336  integer i, k, it, n
1337 
1338  tmin = tice-160.
1339  eps10 = 10.*esl
1340 
1341  if( .not. allocated(table) ) call qsmith_init
1342 
1343 #ifdef MULTI_GASES
1344  do k=1,kmx
1345  do i=1,imx
1346 #else
1347  do k=k1,km
1348  do i=1,im
1349 #endif
1350  ap1 = 10.*dim(t(i,k), tmin) + 1.
1351  ap1 = min(2621., ap1)
1352  it = ap1
1353  es(i,k) = table(it) + (ap1-it)*des(it)
1354 #ifdef MULTI_GASES
1355  qs(i,k) = esl*es(i,k)*virq(q(i,k,1:num_gas))/p(i,k)
1356 #else
1357  qs(i,k) = esl*es(i,k)*(1.+zvir*q(i,k))/p(i,k)
1358 #endif
1359  enddo
1360  enddo
1361 
1362  if ( present(dqdt) ) then
1363 #ifdef MULTI_GASES
1364  do k=1,kmx
1365  do i=1,imx
1366 #else
1367  do k=k1,km
1368  do i=1,im
1369 #endif
1370  ap1 = 10.*dim(t(i,k), tmin) + 1.
1371  ap1 = min(2621., ap1) - 0.5
1372  it = ap1
1373 #ifdef MULTI_GASES
1374  dqdt(i,k) = eps10*(des(it)+(ap1-it)*(des(it+1)-des(it)))*virq(q(i,k,1:num_gas))/p(i,k)
1375 #else
1376  dqdt(i,k) = eps10*(des(it)+(ap1-it)*(des(it+1)-des(it)))*(1.+zvir*q(i,k))/p(i,k)
1377 #endif
1378  enddo
1379  enddo
1380  endif
1381 
1382  end subroutine qsmith
1383 
1384 
1385  subroutine qs_table(n,table)
1386  integer, intent(in):: n
1387  real table (n)
1388  real:: dt=0.1
1389  real esbasw, tbasw, esbasi, tbasi, Tmin, tem, aa, b, c, d, e, esh20
1390  real wice, wh2o
1391  integer i
1392 ! Constants
1393  esbasw = 1013246.0
1394  tbasw = 373.16
1395  tbasi = 273.16
1396  tmin = tbasi - 160.
1397 ! Compute es over water
1398 ! see smithsonian meteorological tables page 350.
1399  do i=1,n
1400  tem = tmin+dt*real(i-1)
1401  aa = -7.90298*(tbasw/tem-1)
1402  b = 5.02808*alog10(tbasw/tem)
1403  c = -1.3816e-07*(10**((1-tem/tbasw)*11.344)-1)
1404  d = 8.1328e-03*(10**((tbasw/tem-1)*(-3.49149))-1)
1405  e = alog10(esbasw)
1406  table(i) = 0.1*10**(aa+b+c+d+e)
1407  enddo
1408 
1409  end subroutine qs_table
1410 
1411  subroutine qs_table_m(n,table)
1412 ! Mixed (blended) table
1413  integer, intent(in):: n
1414  real table (n)
1415  real esupc(200)
1416  real:: dt=0.1
1417  real esbasw, tbasw, esbasi, tbasi, Tmin, tem, aa, b, c, d, e, esh20
1418  real wice, wh2o
1419  integer i
1420 
1421 ! Constants
1422  esbasw = 1013246.0
1423  tbasw = 373.16
1424  esbasi = 6107.1
1425  tbasi = 273.16
1426 ! ****************************************************
1427 ! Compute es over ice between -160c and 0 c.
1428  tmin = tbasi - 160.
1429 ! see smithsonian meteorological tables page 350.
1430  do i=1,1600
1431  tem = tmin+dt*real(i-1)
1432  aa = -9.09718 *(tbasi/tem-1.0)
1433  b = -3.56654 *alog10(tbasi/tem)
1434  c = 0.876793*(1.0-tem/tbasi)
1435  e = alog10(esbasi)
1436  table(i)=10**(aa+b+c+e)
1437  enddo
1438 ! *****************************************************
1439 ! Compute es over water between -20c and 102c.
1440 ! see smithsonian meteorological tables page 350.
1441  do i=1,1221
1442  tem = 253.16+dt*real(i-1)
1443  aa = -7.90298*(tbasw/tem-1)
1444  b = 5.02808*alog10(tbasw/tem)
1445  c = -1.3816e-07*(10**((1-tem/tbasw)*11.344)-1)
1446  d = 8.1328e-03*(10**((tbasw/tem-1)*(-3.49149))-1)
1447  e = alog10(esbasw)
1448  esh20 = 10**(aa+b+c+d+e)
1449  if (i <= 200) then
1450  esupc(i) = esh20
1451  else
1452  table(i+1400) = esh20
1453  endif
1454  enddo
1455 !********************************************************************
1456 ! Derive blended es over ice and supercooled water between -20c and 0c
1457  do i=1,200
1458  tem = 253.16+dt*real(i-1)
1459  wice = 0.05*(273.16-tem)
1460  wh2o = 0.05*(tem-253.16)
1461  table(i+1400) = wice*table(i+1400)+wh2o*esupc(i)
1462  enddo
1463 
1464  do i=1,n
1465  table(i) = table(i)*0.1
1466  enddo
1467 
1468  end subroutine qs_table_m
1469 
1470  subroutine neg_adj3(is, ie, js, je, ng, kbot, hydrostatic, peln, delz, pt, dp, &
1471 #ifdef MULTI_GASES
1472  qvi, &
1473 #else
1474  qv, &
1475 #endif
1476  ql, qr, qi, qs, qg, qa, check_negative)
1477 
1478 ! This is designed for 6-class micro-physics schemes
1479  integer, intent(in):: is, ie, js, je, ng, kbot
1480  logical, intent(in):: hydrostatic
1481  real, intent(in):: dp(is-ng:ie+ng,js-ng:je+ng,kbot)
1482  real, intent(in):: delz(is-ng:,js-ng:,1:)
1483  real, intent(in):: peln(is:ie,kbot+1,js:je)
1484  logical, intent(in), OPTIONAL :: check_negative
1485 #ifdef MULTI_GASES
1486  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot,*):: qvi
1487 #else
1488  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qv
1489 #endif
1490  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: &
1491  pt, ql, qr, qi, qs, qg
1492  real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa
1493 ! Local:
1494  logical:: sat_adj = .false.
1495  real, parameter :: t48 = tice - 48.
1496 #ifdef MULTI_GASES
1497  real, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qv
1498 #endif
1499  real, dimension(is:ie,kbot):: dpk, q2
1500  real, dimension(is:ie,js:je):: pt2, qv2, ql2, qi2, qs2, qr2, qg2, dp2, p2, icpk, lcpk
1501  real:: cv_air
1502  real:: dq, qsum, dq1, q_liq, q_sol, cpm, sink, qsw, dwsdt
1503  integer i, j, k
1504 
1505  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
1506 
1507 #ifdef MULTI_GASES
1508  qv(:,:,:) = qvi(:,:,:,1)
1509 #endif
1510 
1511  if ( present(check_negative) ) then
1512  if ( check_negative ) then
1513  call prt_negative('Temperature', pt, is, ie, js, je, ng, kbot, 165.)
1514  call prt_negative('sphum', qv, is, ie, js, je, ng, kbot, -1.e-8)
1515  call prt_negative('liq_wat', ql, is, ie, js, je, ng, kbot, -1.e-7)
1516  call prt_negative('rainwat', qr, is, ie, js, je, ng, kbot, -1.e-7)
1517  call prt_negative('ice_wat', qi, is, ie, js, je, ng, kbot, -1.e-7)
1518  call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7)
1519  call prt_negative('graupel', qg, is, ie, js, je, ng, kbot, -1.e-7)
1520 #ifdef MULTI_GASES
1521  qvi(:,:,:,1) = qv(:,:,:)
1522 #endif
1523  endif
1524  endif
1525 
1526  if ( hydrostatic ) then
1527  d0_vap = cp_vapor - c_liq
1528  else
1529  d0_vap = cv_vap - c_liq
1530  endif
1531  lv00 = hlv0 - d0_vap*t_ice
1532 
1533 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,ql,qi,qs,qr,qg,dp,pt, &
1534 #ifdef MULTI_GASES
1535 !$OMP qvi, num_gas, &
1536 #endif
1537 !$OMP lv00, d0_vap,hydrostatic,peln,delz,cv_air,sat_adj) &
1538 !$OMP private(dq,dq1,qsum,dp2,p2,pt2,qv2,ql2,qi2,qs2,qg2,qr2, &
1539 !$OMP lcpk,icpk,qsw,dwsdt,sink,q_liq,q_sol,cpm)
1540  do k=1, kbot
1541  do j=js, je
1542  do i=is, ie
1543  qv2(i,j) = qv(i,j,k)
1544  ql2(i,j) = ql(i,j,k)
1545  qi2(i,j) = qi(i,j,k)
1546  qs2(i,j) = qs(i,j,k)
1547  qr2(i,j) = qr(i,j,k)
1548  qg2(i,j) = qg(i,j,k)
1549  dp2(i,j) = dp(i,j,k)
1550  pt2(i,j) = pt(i,j,k)
1551  enddo
1552  enddo
1553 
1554  if ( hydrostatic ) then
1555  do j=js, je
1556  do i=is, ie
1557  p2(i,j) = dp2(i,j)/(peln(i,k+1,j)-peln(i,k,j))
1558  q_liq = max(0., ql2(i,j) + qr2(i,j))
1559  q_sol = max(0., qi2(i,j) + qs2(i,j))
1560 #ifdef MULTI_GASES
1561  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cp_air*vicpqd_qpz(qvi(i,j,k,1:num_gas),qv2(i,j)+q_liq+q_sol) + &
1562  qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1563 #else
1564  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cp_air + qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1565 #endif
1566  lcpk(i,j) = hlv / cpm
1567  icpk(i,j) = hlf / cpm
1568  enddo
1569  enddo
1570  else
1571  do j=js, je
1572  do i=is, ie
1573  q_liq = max(0., ql2(i,j) + qr2(i,j))
1574  q_sol = max(0., qi2(i,j) + qs2(i,j))
1575 #ifdef MULTI_GASES
1576  p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*virq(qvi(i,j,k,1:num_gas))
1577  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cv_air*vicvqd_qpz(qvi(i,j,k,1:num_gas),qv2(i,j)+q_liq+q_sol) + &
1578  qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice
1579 #else
1580  p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*(1.+zvir*qv2(i,j))
1581  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cv_air + qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice
1582 #endif
1583  lcpk(i,j) = (lv00+d0_vap*pt2(i,j)) / cpm
1584  icpk(i,j) = (li0+dc_ice*pt2(i,j)) / cpm
1585  enddo
1586  enddo
1587  endif
1588 
1589 ! Fix the negatives:
1590 !-----------
1591 ! Ice-phase:
1592 !-----------
1593  do j=js, je
1594  do i=is, ie
1595  qsum = qi2(i,j) + qs2(i,j)
1596  if ( qsum > 0. ) then
1597  if ( qi2(i,j) < 0. ) then
1598  qi2(i,j) = 0.
1599  qs2(i,j) = qsum
1600  elseif ( qs2(i,j) < 0. ) then
1601  qs2(i,j) = 0.
1602  qi2(i,j) = qsum
1603  endif
1604  else
1605 ! borrow froom graupel
1606  qi2(i,j) = 0.
1607  qs2(i,j) = 0.
1608  qg2(i,j) = qg2(i,j) + qsum
1609  endif
1610 
1611 ! At this stage qi and qs should be positive definite
1612 ! If graupel < 0 then borrow from qs then qi
1613  if ( qg2(i,j) < 0. ) then
1614  dq = min( qs2(i,j), -qg2(i,j) )
1615  qs2(i,j) = qs2(i,j) - dq
1616  qg2(i,j) = qg2(i,j) + dq
1617  if ( qg2(i,j) < 0. ) then
1618 ! if qg is still negative
1619  dq = min( qi2(i,j), -qg2(i,j) )
1620  qi2(i,j) = qi2(i,j) - dq
1621  qg2(i,j) = qg2(i,j) + dq
1622  endif
1623  endif
1624 
1625 ! If qg is still negative then borrow from rain water: phase change
1626  if ( qg2(i,j)<0. .and. qr2(i,j)>0. ) then
1627  dq = min( qr2(i,j), -qg2(i,j) )
1628  qg2(i,j) = qg2(i,j) + dq
1629  qr2(i,j) = qr2(i,j) - dq
1630  pt2(i,j) = pt2(i,j) + dq*icpk(i,j) ! conserve total energy
1631  endif
1632 ! If qg is still negative then borrow from cloud water: phase change
1633  if ( qg2(i,j)<0. .and. ql2(i,j)>0. ) then
1634  dq = min( ql2(i,j), -qg2(i,j) )
1635  qg2(i,j) = qg2(i,j) + dq
1636  ql2(i,j) = ql2(i,j) - dq
1637  pt2(i,j) = pt2(i,j) + dq*icpk(i,j)
1638  endif
1639 ! Last resort; borrow from water vapor
1640  if ( qg2(i,j)<0. .and. qv2(i,j)>0. ) then
1641  dq = min( 0.999*qv2(i,j), -qg2(i,j) )
1642  qg2(i,j) = qg2(i,j) + dq
1643  qv2(i,j) = qv2(i,j) - dq
1644  pt2(i,j) = pt2(i,j) + dq*(icpk(i,j)+lcpk(i,j))
1645  endif
1646 
1647 !--------------
1648 ! Liquid phase:
1649 !--------------
1650  qsum = ql2(i,j) + qr2(i,j)
1651  if ( qsum > 0. ) then
1652  if ( qr2(i,j) < 0. ) then
1653  qr2(i,j) = 0.
1654  ql2(i,j) = qsum
1655  elseif ( ql2(i,j) < 0. ) then
1656  ql2(i,j) = 0.
1657  qr2(i,j) = qsum
1658  endif
1659  else
1660  ql2(i,j) = 0.
1661  qr2(i,j) = qsum ! rain water is still negative
1662 ! fill negative rain with qg first
1663  dq = min( max(0.0, qg2(i,j)), -qr2(i,j) )
1664  qr2(i,j) = qr2(i,j) + dq
1665  qg2(i,j) = qg2(i,j) - dq
1666  pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
1667  if ( qr(i,j,k) < 0. ) then
1668 ! fill negative rain with available qi & qs (cooling)
1669  dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) )
1670  qr2(i,j) = qr2(i,j) + dq
1671  dq1 = min( dq, qs2(i,j) )
1672  qs2(i,j) = qs2(i,j) - dq1
1673  qi2(i,j) = qi2(i,j) + dq1 - dq
1674  pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
1675  endif
1676 ! fix negative rain water with available vapor
1677  if ( qr2(i,j)<0. .and. qv2(i,j)>0. ) then
1678  dq = min( 0.999*qv2(i,j), -qr2(i,j) )
1679  qv2(i,j) = qv2(i,j) - dq
1680  qr2(i,j) = qr2(i,j) + dq
1681  pt2(i,j) = pt2(i,j) + dq*lcpk(i,j)
1682  endif
1683  endif
1684  enddo
1685  enddo
1686 
1687 !******************************************
1688 ! Fast moist physics: Saturation adjustment
1689 !******************************************
1690 #ifndef GFS_PHYS
1691  if ( sat_adj ) then
1692 
1693  do j=js, je
1694  do i=is, ie
1695 ! Melting of cloud ice into cloud water ********
1696  if ( qi2(i,j)>1.e-8 .and. pt2(i,j) > tice ) then
1697  sink = min( qi2(i,j), (pt2(i,j)-tice)/icpk(i,j) )
1698  ql2(i,j) = ql2(i,j) + sink
1699  qi2(i,j) = qi2(i,j) - sink
1700  pt2(i,j) = pt2(i,j) - sink*icpk(i,j)
1701  endif
1702 
1703 ! vapor <---> liquid water --------------------------------
1704  qsw = wqsat2_moist(pt2(i,j), qv2(i,j), p2(i,j), dwsdt)
1705  sink = min( ql2(i,j), (qsw-qv2(i,j))/(1.+lcpk(i,j)*dwsdt) )
1706  qv2(i,j) = qv2(i,j) + sink
1707  ql2(i,j) = ql2(i,j) - sink
1708  pt2(i,j) = pt2(i,j) - sink*lcpk(i,j)
1709 !-----------------------------------------------------------
1710 
1711 ! freezing of cloud water ********
1712  if( ql2(i,j)>1.e-8 .and. pt2(i,j) < t48 ) then
1713 ! Enforce complete freezing below t_00 (-48 C)
1714  sink = min( ql2(i,j), (t48-pt2(i,j))/icpk(i,j) )
1715  ql2(i,j) = ql2(i,j) - sink
1716  qi2(i,j) = qi2(i,j) + sink
1717  pt2(i,j) = pt2(i,j) + sink*icpk(i,j)
1718  endif ! significant ql existed
1719  enddo
1720  enddo
1721  endif
1722 #endif
1723 
1724 !----------------------------------------------------------------
1725 ! Update fields:
1726  do j=js, je
1727  do i=is, ie
1728  qv(i,j,k) = qv2(i,j)
1729  ql(i,j,k) = ql2(i,j)
1730  qi(i,j,k) = qi2(i,j)
1731  qs(i,j,k) = qs2(i,j)
1732  qr(i,j,k) = qr2(i,j)
1733  qg(i,j,k) = qg2(i,j)
1734  pt(i,j,k) = pt2(i,j)
1735  enddo
1736  enddo
1737 #ifdef MULTI_GASES
1738  qvi(:,:,k,1) = qv(:,:,k)
1739 #endif
1740 
1741  enddo
1742 
1743 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,dp,qg,qr) &
1744 !$OMP private(dpk, q2)
1745  do j=js, je
1746 ! Graupel:
1747  do k=1,kbot
1748  do i=is,ie
1749  dpk(i,k) = dp(i,j,k)
1750  q2(i,k) = qg(i,j,k)
1751  enddo
1752  enddo
1753  call fillq(ie-is+1, kbot, q2, dpk)
1754  do k=1,kbot
1755  do i=is,ie
1756  qg(i,j,k) = q2(i,k)
1757  enddo
1758  enddo
1759 ! Rain water:
1760  do k=1,kbot
1761  do i=is,ie
1762  q2(i,k) = qr(i,j,k)
1763  enddo
1764  enddo
1765  call fillq(ie-is+1, kbot, q2, dpk)
1766  do k=1,kbot
1767  do i=is,ie
1768  qr(i,j,k) = q2(i,k)
1769  enddo
1770  enddo
1771  enddo
1772 
1773 !-----------------------------------
1774 ! Fix water vapor
1775 !-----------------------------------
1776 ! Top layer: borrow from below
1777  k = 1
1778 !$OMP parallel do default(none) shared(is,ie,js,je,k,qv,dp)
1779  do j=js, je
1780  do i=is, ie
1781  if( qv(i,j,k) < 0. ) then
1782  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1783  qv(i,j,k ) = 0.
1784  endif
1785  enddo
1786  enddo
1787 
1788 ! this OpenMP do-loop cannot be parallelized with recursion on k/k-1
1789 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) &
1790 !$OMP private(dq)
1791  do j=js, je
1792  do k=2,kbot-1
1793  do i=is, ie
1794  if( qv(i,j,k) < 0. .and. qv(i,j,k-1) > 0. ) then
1795  dq = min(-qv(i,j,k)*dp(i,j,k), qv(i,j,k-1)*dp(i,j,k-1))
1796  qv(i,j,k-1) = qv(i,j,k-1) - dq/dp(i,j,k-1)
1797  qv(i,j,k ) = qv(i,j,k ) + dq/dp(i,j,k )
1798  endif
1799  if( qv(i,j,k) < 0. ) then
1800  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1801  qv(i,j,k ) = 0.
1802  endif
1803  enddo
1804  enddo
1805  enddo
1806 
1807 ! Bottom layer; Borrow from above
1808 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) private(dq)
1809  do j=js, je
1810  do i=is, ie
1811  if( qv(i,j,kbot) < 0. ) then
1812  do k=kbot-1,1,-1
1813  if ( qv(i,j,kbot)>=0. ) goto 123
1814  if ( qv(i,j,k) > 0. ) then
1815  dq = min(-qv(i,j,kbot)*dp(i,j,kbot), qv(i,j,k)*dp(i,j,k))
1816  qv(i,j,k ) = qv(i,j,k ) - dq/dp(i,j,k)
1817  qv(i,j,kbot) = qv(i,j,kbot) + dq/dp(i,j,kbot)
1818  endif
1819  enddo ! k-loop
1820 123 continue
1821  endif
1822  enddo ! i-loop
1823  enddo ! j-loop
1824 #ifdef MULTI_GASES
1825  qvi(:,:,:,1) = qv(:,:,:)
1826 #endif
1827 
1828 
1829  if (present(qa)) then
1830 !-----------------------------------
1831 ! Fix negative cloud fraction
1832 !-----------------------------------
1833 ! this OpenMP do-loop cannot be parallelized by the recursion on k/k+1
1834 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qa,dp)
1835  do j=js, je
1836  do k=1,kbot-1
1837  do i=is, ie
1838  if( qa(i,j,k) < 0. ) then
1839  qa(i,j,k+1) = qa(i,j,k+1) + qa(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1840  qa(i,j,k ) = 0.
1841  endif
1842  enddo
1843  enddo
1844  enddo
1845 
1846 ! Bottom layer; Borrow from above
1847 !$OMP parallel do default(none) shared(is,ie,js,je,qa,kbot,dp) &
1848 !$OMP private(dq)
1849  do j=js, je
1850  do i=is, ie
1851  if( qa(i,j,kbot) < 0. .and. qa(i,j,kbot-1)>0.) then
1852  dq = min(-qa(i,j,kbot)*dp(i,j,kbot), qa(i,j,kbot-1)*dp(i,j,kbot-1))
1853  qa(i,j,kbot-1) = qa(i,j,kbot-1) - dq/dp(i,j,kbot-1)
1854  qa(i,j,kbot ) = qa(i,j,kbot ) + dq/dp(i,j,kbot )
1855  endif
1856 ! if qa is still < 0
1857  qa(i,j,kbot) = max(0., qa(i,j,kbot))
1858  enddo
1859  enddo
1860 
1861  endif
1862 
1863  end subroutine neg_adj3
1864 
1865  subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, &
1866  peln, delz, pt, dp, qv, ql, qr, qi, qs, qa, check_negative)
1867 
1868 ! This is designed for 6-class micro-physics schemes
1869  integer, intent(in):: is, ie, js, je, ng, kbot
1870  logical, intent(in):: hydrostatic
1871  real, intent(in):: dp(is-ng:ie+ng,js-ng:je+ng,kbot)
1872  real, intent(in):: delz(is-ng:,js-ng:,1:)
1873  real, intent(in):: peln(is:ie,kbot+1,js:je)
1874  logical, intent(in), OPTIONAL :: check_negative
1875  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: &
1876  pt, qv, ql, qr, qi, qs
1877  real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa
1878 ! Local:
1879  logical:: sat_adj = .false.
1880  real, parameter :: t48 = tice - 48.
1881  real, dimension(is:ie,kbot):: dpk, q2
1882  real, dimension(is:ie,js:je):: pt2, qv2, ql2, qi2, qs2, qr2, dp2, p2, icpk, lcpk
1883  real:: cv_air
1884  real:: dq, qsum, dq1, q_liq, q_sol, oneocpm, sink, qsw, dwsdt, tx1
1885  integer i, j, k
1886 
1887  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
1888 
1889  if ( present(check_negative) ) then
1890  if ( check_negative ) then
1891  call prt_negative('Temperature', pt, is, ie, js, je, ng, kbot, 165.)
1892  call prt_negative('sphum', qv, is, ie, js, je, ng, kbot, -1.e-8)
1893  call prt_negative('liq_wat', ql, is, ie, js, je, ng, kbot, -1.e-7)
1894  call prt_negative('rainwat', qr, is, ie, js, je, ng, kbot, -1.e-7)
1895  call prt_negative('ice_wat', qi, is, ie, js, je, ng, kbot, -1.e-7)
1896  call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7)
1897  endif
1898  endif
1899 
1900  if ( hydrostatic ) then
1901  d0_vap = cp_vapor - c_liq
1902  else
1903  d0_vap = cv_vap - c_liq
1904  endif
1905  lv00 = hlv0 - d0_vap*t_ice
1906 
1907 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,ql,qi,qs,qr,dp,pt, &
1908 !$OMP lv00, d0_vap,hydrostatic,peln,delz,cv_air,sat_adj) &
1909 !$OMP private(dq,dq1,qsum,dp2,p2,pt2,qv2,ql2,qi2,qs2,qr2, &
1910 !$OMP lcpk,icpk,qsw,dwsdt,sink,q_liq,q_sol,oneocpm)
1911  do k=1, kbot
1912  do j=js, je
1913  do i=is, ie
1914  qv2(i,j) = qv(i,j,k)
1915  ql2(i,j) = ql(i,j,k)
1916  qi2(i,j) = qi(i,j,k)
1917  qs2(i,j) = qs(i,j,k)
1918  qr2(i,j) = qr(i,j,k)
1919  dp2(i,j) = dp(i,j,k)
1920  pt2(i,j) = pt(i,j,k)
1921  enddo
1922  enddo
1923 
1924  if ( hydrostatic ) then
1925  do j=js, je
1926  do i=is, ie
1927  p2(i,j) = dp2(i,j)/(peln(i,k+1,j)-peln(i,k,j))
1928  q_liq = max(0., ql2(i,j) + qr2(i,j))
1929  q_sol = max(0., qi2(i,j) + qs2(i,j))
1930  oneocpm = 1.0 / ((1.-(qv2(i,j)+q_liq+q_sol))*cp_air + qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice)
1931  lcpk(i,j) = hlv * oneocpm
1932  icpk(i,j) = hlf * oneocpm
1933  enddo
1934  enddo
1935  else
1936  do j=js, je
1937  do i=is, ie
1938  p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*(1.+zvir*qv2(i,j))
1939  q_liq = max(0., ql2(i,j) + qr2(i,j))
1940  q_sol = max(0., qi2(i,j) + qs2(i,j))
1941  oneocpm = 1.0 / ((1.-(qv2(i,j)+q_liq+q_sol))*cv_air + qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice)
1942  lcpk(i,j) = (lv00+d0_vap*pt2(i,j)) * oneocpm
1943  icpk(i,j) = (li0+dc_ice*pt2(i,j)) * oneocpm
1944  enddo
1945  enddo
1946  endif
1947 
1948 ! Fix the negatives:
1949 !-----------
1950 ! Ice-phase:
1951 !-----------
1952  do j=js, je
1953  do i=is, ie
1954  qsum = qi2(i,j) + qs2(i,j)
1955  if ( qsum > 0. ) then
1956  if ( qi2(i,j) < 0. ) then
1957  qi2(i,j) = 0.
1958  qs2(i,j) = qsum
1959  elseif ( qs2(i,j) < 0. ) then
1960  qs2(i,j) = 0.
1961  qi2(i,j) = qsum
1962  endif
1963  else
1964  qi2(i,j) = 0.
1965  qs2(i,j) = qsum
1966 
1967 ! If qsum is negative then borrow from rain water: phase change
1968  if ( qs2(i,j) < 0. .and. qr2(i,j) > 0. ) then
1969  dq = min( qr2(i,j), -qs2(i,j) )
1970  qs2(i,j) = qs2(i,j) + dq
1971  qr2(i,j) = qr2(i,j) - dq
1972  pt2(i,j) = pt2(i,j) + dq*icpk(i,j) ! conserve total energy
1973  endif
1974 ! If qs2 is still negative then borrow from cloud water: phase change
1975  if ( qs2(i,j) < 0. .and. ql2(i,j) > 0. ) then
1976  dq = min( ql2(i,j), -qs2(i,j) )
1977  qs2(i,j) = qs2(i,j) + dq
1978  ql2(i,j) = ql2(i,j) - dq
1979  pt2(i,j) = pt2(i,j) + dq*icpk(i,j)
1980  endif
1981 ! Last resort; borrow from water vapor
1982  if ( qs2(i,j) < 0. .and. qv2(i,j) > 0. ) then
1983  dq = min( 0.999*qv2(i,j), -qs2(i,j) )
1984  qs2(i,j) = qs2(i,j) + dq
1985  qv2(i,j) = qv2(i,j) - dq
1986  pt2(i,j) = pt2(i,j) + dq*(icpk(i,j)+lcpk(i,j))
1987  endif
1988  endif
1989 
1990 !--------------
1991 ! Liquid phase:
1992 !--------------
1993  qsum = ql2(i,j) + qr2(i,j)
1994  if ( qsum > 0. ) then
1995  if ( qr2(i,j) < 0. ) then
1996  qr2(i,j) = 0.
1997  ql2(i,j) = qsum
1998  elseif ( ql2(i,j) < 0. ) then
1999  ql2(i,j) = 0.
2000  qr2(i,j) = qsum
2001  endif
2002  else
2003  ql2(i,j) = 0.
2004  qr2(i,j) = qsum ! rain water is still negative
2005  if ( qr(i,j,k) < 0. ) then
2006 ! fill negative rain with available qi & qs (cooling)
2007  dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) )
2008  qr2(i,j) = qr2(i,j) + dq
2009  dq1 = min( dq, qs2(i,j) )
2010  qs2(i,j) = qs2(i,j) - dq1
2011  qi2(i,j) = qi2(i,j) + dq1 - dq
2012  pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
2013  endif
2014 ! fix negative rain water with available vapor
2015  if ( qr2(i,j) < 0. .and. qv2(i,j) > 0. ) then
2016  dq = min( 0.999*qv2(i,j), -qr2(i,j) )
2017  qv2(i,j) = qv2(i,j) - dq
2018  qr2(i,j) = qr2(i,j) + dq
2019  pt2(i,j) = pt2(i,j) + dq*lcpk(i,j)
2020  endif
2021  endif
2022  enddo
2023  enddo
2024 
2025 !******************************************
2026 ! Fast moist physics: Saturation adjustment
2027 !******************************************
2028 #ifndef GFS_PHYS
2029  if ( sat_adj ) then
2030 
2031  do j=js, je
2032  do i=is, ie
2033 ! Melting of cloud ice into cloud water ********
2034  if ( qi2(i,j)>1.e-8 .and. pt2(i,j) > tice ) then
2035  sink = min( qi2(i,j), (pt2(i,j)-tice)/icpk(i,j) )
2036  ql2(i,j) = ql2(i,j) + sink
2037  qi2(i,j) = qi2(i,j) - sink
2038  pt2(i,j) = pt2(i,j) - sink*icpk(i,j)
2039  endif
2040 
2041 ! vapor <---> liquid water --------------------------------
2042  qsw = wqsat2_moist(pt2(i,j), qv2(i,j), p2(i,j), dwsdt)
2043  sink = min( ql2(i,j), (qsw-qv2(i,j))/(1.+lcpk(i,j)*dwsdt) )
2044  qv2(i,j) = qv2(i,j) + sink
2045  ql2(i,j) = ql2(i,j) - sink
2046  pt2(i,j) = pt2(i,j) - sink*lcpk(i,j)
2047 !-----------------------------------------------------------
2048 
2049 ! freezing of cloud water ********
2050  if( ql2(i,j)>1.e-8 .and. pt2(i,j) < t48 ) then
2051 ! Enforce complete freezing below t_00 (-48 C)
2052  sink = min( ql2(i,j), (t48-pt2(i,j))/icpk(i,j) )
2053  ql2(i,j) = ql2(i,j) - sink
2054  qi2(i,j) = qi2(i,j) + sink
2055  pt2(i,j) = pt2(i,j) + sink*icpk(i,j)
2056  endif ! significant ql existed
2057  enddo
2058  enddo
2059  endif
2060 #endif
2061 
2062 !----------------------------------------------------------------
2063 ! Update fields:
2064  do j=js, je
2065  do i=is, ie
2066  qv(i,j,k) = qv2(i,j)
2067  ql(i,j,k) = ql2(i,j)
2068  qi(i,j,k) = qi2(i,j)
2069  qs(i,j,k) = qs2(i,j)
2070  qr(i,j,k) = qr2(i,j)
2071  pt(i,j,k) = pt2(i,j)
2072  enddo
2073  enddo
2074 
2075  enddo
2076 
2077 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,dp,qr) &
2078 !$OMP private(dpk, q2)
2079  do j=js, je
2080 ! Rain water:
2081  do k=1,kbot
2082  do i=is,ie
2083  dpk(i,k) = dp(i,j,k)
2084  q2(i,k) = qr(i,j,k)
2085  enddo
2086  enddo
2087  call fillq(ie-is+1, kbot, q2, dpk)
2088  do k=1,kbot
2089  do i=is,ie
2090  qr(i,j,k) = q2(i,k)
2091  enddo
2092  enddo
2093  enddo
2094 
2095 !-----------------------------------
2096 ! Fix water vapor
2097 !-----------------------------------
2098 ! Top layer: borrow from below
2099  k = 1
2100 !$OMP parallel do default(none) shared(is,ie,js,je,k,qv,dp)
2101  do j=js, je
2102  do i=is, ie
2103  if( qv(i,j,k) < 0. ) then
2104  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
2105  qv(i,j,k ) = 0.
2106  endif
2107  enddo
2108  enddo
2109 
2110 ! this OpenMP do-loop cannot be parallelized with recursion on k/k-1
2111 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) &
2112 !$OMP private(dq)
2113  do j=js, je
2114  do k=2,kbot-1
2115  do i=is, ie
2116  if( qv(i,j,k) < 0. .and. qv(i,j,k-1) > 0. ) then
2117  dq = min(-qv(i,j,k)*dp(i,j,k), qv(i,j,k-1)*dp(i,j,k-1))
2118  qv(i,j,k-1) = qv(i,j,k-1) - dq/dp(i,j,k-1)
2119  qv(i,j,k ) = qv(i,j,k ) + dq/dp(i,j,k )
2120  endif
2121  if( qv(i,j,k) < 0. ) then
2122  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
2123  qv(i,j,k ) = 0.
2124  endif
2125  enddo
2126  enddo
2127  enddo
2128 
2129 ! Bottom layer; Borrow from above
2130 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) private(dq,tx1)
2131  do j=js, je
2132  do i=is, ie
2133  if( qv(i,j,kbot) < 0. ) then
2134  tx1 = 1.0 / dp(i,j,kbot)
2135  do k=kbot-1,1,-1
2136  if ( qv(i,j,kbot)>= 0. ) goto 123
2137  if ( qv(i,j,k) > 0. ) then
2138  dq = min(-qv(i,j,kbot)*dp(i,j,kbot), qv(i,j,k)*dp(i,j,k))
2139  qv(i,j,k ) = qv(i,j,k ) - dq / dp(i,j,k)
2140  qv(i,j,kbot) = qv(i,j,kbot) + dq * tx1
2141  endif
2142  enddo ! k-loop
2143 123 continue
2144  endif
2145  enddo ! i-loop
2146  enddo ! j-loop
2147 
2148 
2149  if (present(qa)) then
2150 !-----------------------------------
2151 ! Fix negative cloud fraction
2152 !-----------------------------------
2153 ! this OpenMP do-loop cannot be parallelized by the recursion on k/k+1
2154 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qa,dp)
2155  do j=js, je
2156  do k=1,kbot-1
2157  do i=is, ie
2158  if( qa(i,j,k) < 0. ) then
2159  qa(i,j,k+1) = qa(i,j,k+1) + qa(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
2160  qa(i,j,k ) = 0.
2161  endif
2162  enddo
2163  enddo
2164  enddo
2165 
2166 ! Bottom layer; Borrow from above
2167 !$OMP parallel do default(none) shared(is,ie,js,je,qa,kbot,dp) &
2168 !$OMP private(dq)
2169  do j=js, je
2170  do i=is, ie
2171  if( qa(i,j,kbot) < 0. .and. qa(i,j,kbot-1)>0.) then
2172  dq = min(-qa(i,j,kbot)*dp(i,j,kbot), qa(i,j,kbot-1)*dp(i,j,kbot-1))
2173  qa(i,j,kbot-1) = qa(i,j,kbot-1) - dq/dp(i,j,kbot-1)
2174  qa(i,j,kbot ) = qa(i,j,kbot ) + dq/dp(i,j,kbot )
2175  endif
2176 ! if qa is still < 0
2177  qa(i,j,kbot) = max(0., qa(i,j,kbot))
2178  enddo
2179  enddo
2180 
2181  endif
2182 
2183  end subroutine neg_adj2
2184 
2185  subroutine fillq(im, km, q, dp)
2186 ! Aggresive 1D filling algorithm for qr and qg
2187  integer, intent(in):: im, km
2188  real, intent(inout), dimension(im,km):: q, dp
2189  integer:: i, k
2190  real:: sum1, sum2, dq
2191 
2192  do 500 i=1,im
2193  sum1 = 0.
2194  do k=1,km
2195  if ( q(i,k)>0. ) then
2196  sum1 = sum1 + q(i,k)*dp(i,k)
2197  endif
2198  enddo
2199  if ( sum1<1.e-12 ) goto 500
2200  sum2 = 0.
2201  do k=km,1,-1
2202  if ( q(i,k)<0.0 .and. sum1>0. ) then
2203  dq = min( sum1, -q(i,k)*dp(i,k) )
2204  sum1 = sum1 - dq
2205  sum2 = sum2 + dq
2206  q(i,k) = q(i,k) + dq/dp(i,k)
2207  endif
2208  enddo
2209  do k=km,1,-1
2210  if ( q(i,k)>0.0 .and. sum2>0. ) then
2211  dq = min( sum2, q(i,k)*dp(i,k) )
2212  sum2 = sum2 - dq
2213  q(i,k) = q(i,k) - dq/dp(i,k)
2214  endif
2215  enddo
2216 500 continue
2217 
2218  end subroutine fillq
2219 
2220  subroutine prt_negative(qname, q, is, ie, js, je, n_g, km, threshold)
2221  character(len=*), intent(in):: qname
2222  integer, intent(in):: is, ie, js, je
2223  integer, intent(in):: n_g, km
2224  real, intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
2225  real, intent(in):: threshold
2226  real qmin
2227  integer i,j,k
2228 
2229  qmin = q(is,js,1)
2230  do k=1,km
2231  do j=js,je
2232  do i=is,ie
2233  qmin = min(qmin, q(i,j,k))
2234  enddo
2235  enddo
2236  enddo
2237  call mp_reduce_min(qmin)
2238  if(is_master() .and. qmin<threshold) write(6,*) qname, ' min (negative) = ', qmin
2239 
2240  end subroutine prt_negative
2241 
2242 
2243 end module fv_sg_mod
fv_sg_mod::t3_max
real, parameter t3_max
Definition: fv_sg.F90:97
fv_sg_mod::dc_ice
real, parameter dc_ice
= 2112.
Definition: fv_sg.F90:85
fv_sg_mod::qsmith_init
subroutine qsmith_init
Definition: fv_sg.F90:1296
multi_gases_mod::vicpqd
pure real function, public vicpqd(q)
Definition: multi_gases.F90:239
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
fv_sg_mod::neg_adj3
subroutine, public neg_adj3(is, ie, js, je, ng, kbot, hydrostatic, peln, delz, pt, dp, qv, ql, qr, qi, qs, qg, qa, check_negative)
Definition: fv_sg.F90:1477
fv_sg_mod::li0
real, parameter li0
= -2.431928e5
Definition: fv_sg.F90:99
multi_gases_mod::num_gas
integer, public num_gas
Definition: multi_gases.F90:44
fv_sg_mod::c_ice
real, parameter c_ice
-15 C
Definition: fv_sg.F90:77
multi_gases_mod
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
fv_sg_mod::cv_vap
real, parameter cv_vap
1384.5
Definition: fv_sg.F90:80
multi_gases_mod::vicvqd_qpz
pure real function, public vicvqd_qpz(q, qpz)
Definition: multi_gases.F90:300
multi_gases_mod::virq_qpz
pure real function, public virq_qpz(q, qpz)
Definition: multi_gases.F90:197
fv_sg_mod::lv00
real lv00
Definition: fv_sg.F90:103
fv_sg_mod
The module &#39;fv_sg&#39; performs FV sub-grid mixing.
Definition: fv_sg.F90:54
multi_gases_mod::virqd
pure real function, public virqd(q)
Definition: multi_gases.F90:219
multi_gases_mod::virq
pure real function, public virq(q)
Definition: multi_gases.F90:136
fv_sg_mod::zvir
real, parameter zvir
= 0.607789855
Definition: fv_sg.F90:101
fv_sg_mod::d0_vap
real d0_vap
Definition: fv_sg.F90:103
fv_sg_mod::tice
real, parameter tice
Definition: fv_sg.F90:75
fv_sg_mod::ri_min
real, parameter ri_min
Definition: fv_sg.F90:93
fv_sg_mod::fillq
subroutine fillq(im, km, q, dp)
Definition: fv_sg.F90:2186
fv_sg_mod::qsmith
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
Definition: fv_sg.F90:1320
fv_sg_mod::c_con
real, parameter c_con
Definition: fv_sg.F90:81
fv_sg_mod::t2_min
real, parameter t2_min
Definition: fv_sg.F90:95
fv_sg_mod::hlv0
real, parameter hlv0
Definition: fv_sg.F90:87
fv_sg_mod::c_liq
real, parameter c_liq
GFS.
Definition: fv_sg.F90:78
fv_sg_mod::t1_min
real, parameter t1_min
Definition: fv_sg.F90:94
fv_sg_mod::lv0
real, parameter lv0
= 3.147782e6
Definition: fv_sg.F90:98
fv_sg_mod::t2_max
real, parameter t2_max
Definition: fv_sg.F90:96
fv_sg_mod::hlf0
real, parameter hlf0
Definition: fv_sg.F90:88
fv_sg_mod::prt_negative
subroutine prt_negative(qname, q, is, ie, js, je, n_g, km, threshold)
Definition: fv_sg.F90:2221
fv_sg_mod::t_ice
real, parameter t_ice
Definition: fv_sg.F90:91
fv_sg_mod::fv_subgrid_z
subroutine, public fv_subgrid_z(isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, hydrostatic, w, delz, u_dt, v_dt, t_dt, q_dt, k_bot)
Definition: fv_sg.F90:680
multi_gases_mod::vicpqd_qpz
pure real function, public vicpqd_qpz(q, qpz)
Definition: multi_gases.F90:259
fv_sg_mod::qs_table
subroutine qs_table(n, table)
Definition: fv_sg.F90:1386
fv_sg_mod::table
real, dimension(:), allocatable table
Definition: fv_sg.F90:102
fv_sg_mod::des
real, dimension(:), allocatable des
Definition: fv_sg.F90:102
fv_sg_mod::ri_max
real, parameter ri_max
Definition: fv_sg.F90:92
fv_sg_mod::neg_adj2
subroutine, public neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, peln, delz, pt, dp, qv, ql, qr, qi, qs, qa, check_negative)
Definition: fv_sg.F90:1867
fv_sg_mod::qs_table_m
subroutine qs_table_m(n, table)
Definition: fv_sg.F90:1412
fv_sg_mod::esl
real, parameter esl
Definition: fv_sg.F90:74
multi_gases_mod::vicvqd
pure real function, public vicvqd(q)
Definition: multi_gases.F90:280
fv_sg_mod::dc_vap
real, parameter dc_vap
= -2368.
Definition: fv_sg.F90:84