FV3 Bundle
fv_sg_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
20 module fv_sg_nlm_mod
21 
22 !-----------------------------------------------------------------------
23 ! FV sub-grid mixing
24 !-----------------------------------------------------------------------
25  use constants_mod, only: rdgas, rvgas, cp_air, cp_vapor, hlv, hlf, kappa, grav
26  use tracer_manager_mod, only: get_tracer_index
27  use field_manager_mod, only: model_atmos
28  use lin_cld_microphys_mod, only: wqs2, wqsat2_moist
29  use fv_mp_nlm_mod, only: mp_reduce_min, is_master
30 
31 implicit none
32 private
33 
34 public fv_subgrid_z, qsmith, neg_adj3
35 
36  real, parameter:: esl = 0.621971831
37  real, parameter:: tice = 273.16
38 ! real, parameter:: c_ice = 2106. ! Emanuel table, page 566
39  real, parameter:: c_ice = 1972. ! -15 C
40  real, parameter:: c_liq = 4.1855e+3 ! GFS
41 ! real, parameter:: c_liq = 4218. ! ECMWF-IFS
42  real, parameter:: cv_vap = cp_vapor - rvgas ! 1384.5
43  real, parameter:: c_con = c_ice
44 
45 ! real, parameter:: dc_vap = cp_vapor - c_liq ! = -2368.
46  real, parameter:: dc_vap = cv_vap - c_liq ! = -2368.
47  real, parameter:: dc_ice = c_liq - c_ice ! = 2112.
48 ! Values at 0 Deg C
49  real, parameter:: hlv0 = 2.5e6
50  real, parameter:: hlf0 = 3.3358e5
51 ! real, parameter:: hlv0 = 2.501e6 ! Emanual Appendix-2
52 ! real, parameter:: hlf0 = 3.337e5 ! Emanual
53  real, parameter:: t_ice = 273.16
54  real, parameter:: ri_max = 1.
55  real, parameter:: ri_min = 0.25
56  real, parameter:: t1_min = 160.
57  real, parameter:: t2_min = 165.
58  real, parameter:: t2_max = 315.
59  real, parameter:: t3_max = 325.
60  real, parameter:: lv0 = hlv0 - dc_vap*t_ice ! = 3.147782e6
61  real, parameter:: li0 = hlf0 - dc_ice*t_ice ! = -2.431928e5
62 
63  real, parameter:: zvir = rvgas/rdgas - 1. ! = 0.607789855
64  real, allocatable:: table(:),des(:)
65  real:: lv00, d0_vap
66 
67 contains
68 
69 #if defined(GFS_PHYS) || defined(MAPL_MODE)
70  subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
71  tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, &
72  hydrostatic, w, delz, u_dt, v_dt, t_dt, k_bot )
73 ! Dry convective adjustment-mixing
74 !-------------------------------------------
75  integer, intent(in):: is, ie, js, je, km, nq, nwat
76  integer, intent(in):: isd, ied, jsd, jed
77  integer, intent(in):: tau ! Relaxation time scale
78  real, intent(in):: dt ! model time step
79  real, intent(in):: pe(is-1:ie+1,km+1,js-1:je+1)
80  real, intent(in):: peln(is :ie, km+1,js :je)
81  real, intent(in):: delp(isd:ied,jsd:jed,km) ! Delta p at each model level
82  real, intent(in):: delz(isd:,jsd:,1:) ! Delta z at each model level
83  real, intent(in):: pkz(is:ie,js:je,km)
84  logical, intent(in):: hydrostatic
85  integer, intent(in), optional:: k_bot
86 !
87  real, intent(inout):: ua(isd:ied,jsd:jed,km)
88  real, intent(inout):: va(isd:ied,jsd:jed,km)
89  real, intent(inout):: w(isd:,jsd:,1:)
90  real, intent(inout):: ta(isd:ied,jsd:jed,km) ! Temperature
91  real, intent(inout):: qa(isd:ied,jsd:jed,km,nq) ! Specific humidity & tracers
92  real, intent(inout):: u_dt(isd:ied,jsd:jed,km)
93  real, intent(inout):: v_dt(isd:ied,jsd:jed,km)
94  real, intent(inout):: t_dt(is:ie,js:je,km)
95 !---------------------------Local variables-----------------------------
96  real, dimension(is:ie,km):: u0, v0, w0, t0, hd, te, gz, tvm, pm, den
97  real q0(is:ie,km,nq), qcon(is:ie,km)
98  real, dimension(is:ie):: gzh, lcp2, icp2, cvm, cpm, qs
99  real ri_ref, ri, pt1, pt2, ratio, tv, cv, tmp, q_liq, q_sol
100  real tv1, tv2, g2, h0, mc, fra, rk, rz, rdt, tvd, tv_surf
101  real dh, dq, qsw, dqsdt, tcp3, t_max, t_min
102  integer i, j, k, kk, n, m, iq, km1, im, kbot
103  real, parameter:: ustar2 = 1.e-4
104  real:: cv_air, xvir
105  integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
106 
107  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
108  rk = cp_air/rdgas + 1.
109  cv = cp_air - rdgas
110 
111  g2 = 0.5*grav
112 
113  rdt = 1./ dt
114  im = ie-is+1
115 
116  if ( present(k_bot) ) then
117  if ( k_bot < 3 ) return
118  kbot = k_bot
119  else
120  kbot = km
121  endif
122  if ( pe(is,1,js) < 2. ) then
123  t_min = t1_min
124  else
125  t_min = t2_min
126  endif
127 
128  if ( k_bot < min(km,24) ) then
129  t_max = t2_max
130  else
131  t_max = t3_max
132  endif
133 
134 #ifdef MAPL_MODE
135  sphum = 1
136  if ( nwat == 0 ) then
137  xvir = 0.
138  rz = 0.
139  else
140  xvir = zvir
141  rz = rvgas - rdgas ! rz = zvir * rdgas
142  if ( nwat == 3) then
143  liq_wat = 2
144  ice_wat = 3
145  endif
146  endif
147 #else
148  sphum = get_tracer_index(model_atmos, 'sphum')
149  if ( nwat == 0 ) then
150  xvir = 0.
151  rz = 0.
152  else
153  xvir = zvir
154  rz = rvgas - rdgas ! rz = zvir * rdgas
155  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
156  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
157  rainwat = get_tracer_index(model_atmos, 'rainwat')
158  snowwat = get_tracer_index(model_atmos, 'snowwat')
159  graupel = get_tracer_index(model_atmos, 'graupel')
160  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
161  endif
162 #endif
163 
164 !------------------------------------------------------------------------
165 ! The nonhydrostatic pressure changes if there is heating (under constant
166 ! volume and mass is locally conserved).
167 !------------------------------------------------------------------------
168  m = 3
169  fra = dt/real(tau)
170 
171 !$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
172 !$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
173 !$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra, t_max, t_min, &
174 !$OMP u_dt,rdt,v_dt,xvir,nwat) &
175 !$OMP private(kk,lcp2,icp2,tcp3,dh,dq,den,qs,qsw,dqsdt,qcon,q0, &
176 !$OMP t0,u0,v0,w0,h0,pm,gzh,tvm,tmp,cpm,cvm,q_liq,q_sol, &
177 !$OMP tv,gz,hd,te,ratio,pt1,pt2,tv1,tv2,ri_ref, ri,mc,km1)
178  do 1000 j=js,je
179 
180  do iq=1, nq
181  do k=1,kbot
182  do i=is,ie
183  q0(i,k,iq) = qa(i,j,k,iq)
184  enddo
185  enddo
186  enddo
187 
188  do k=1,kbot
189  do i=is,ie
190  t0(i,k) = ta(i,j,k)
191  tvm(i,k) = t0(i,k)*(1.+xvir*q0(i,k,sphum))
192  u0(i,k) = ua(i,j,k)
193  v0(i,k) = va(i,j,k)
194  pm(i,k) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
195  enddo
196  enddo
197 
198  do i=is,ie
199  gzh(i) = 0.
200  enddo
201 
202  if( hydrostatic ) then
203  do k=kbot, 1,-1
204  do i=is,ie
205  tv = rdgas*tvm(i,k)
206  den(i,k) = pm(i,k)/tv
207  gz(i,k) = gzh(i) + tv*(1.-pe(i,k,j)/pm(i,k))
208  hd(i,k) = cp_air*tvm(i,k)+gz(i,k)+0.5*(u0(i,k)**2+v0(i,k)**2)
209  gzh(i) = gzh(i) + tv*(peln(i,k+1,j)-peln(i,k,j))
210  enddo
211  enddo
212  else
213  do k=kbot, 1, -1
214  if ( nwat == 0 ) then
215  do i=is,ie
216  cpm(i) = cp_air
217  cvm(i) = cv_air
218  enddo
219  elseif ( nwat==1 ) then
220  do i=is,ie
221  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
222  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
223  enddo
224  elseif ( nwat==2 ) then ! GFS
225  do i=is,ie
226  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
227  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
228  enddo
229  elseif ( nwat==3 ) then
230  do i=is,ie
231  q_liq = q0(i,k,liq_wat)
232  q_sol = q0(i,k,ice_wat)
233  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
234  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
235  enddo
236  elseif ( nwat==4 ) then
237  do i=is,ie
238  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
239  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
240  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq
241  enddo
242  else
243  do i=is,ie
244  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
245  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
246  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
247  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
248  enddo
249  endif
250 
251  do i=is,ie
252  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
253  w0(i,k) = w(i,j,k)
254  gz(i,k) = gzh(i) - g2*delz(i,j,k)
255  tmp = gz(i,k) + 0.5*(u0(i,k)**2+v0(i,k)**2+w0(i,k)**2)
256  hd(i,k) = cpm(i)*t0(i,k) + tmp
257  te(i,k) = cvm(i)*t0(i,k) + tmp
258  gzh(i) = gzh(i) - grav*delz(i,j,k)
259  enddo
260  enddo
261  endif
262 
263  do n=1,m
264 
265  if ( m==3 ) then
266  if ( n==1) ratio = 0.25
267  if ( n==2) ratio = 0.5
268  if ( n==3) ratio = 0.999
269  else
270  ratio = real(n)/real(m)
271  endif
272 
273  do i=is,ie
274  gzh(i) = 0.
275  enddo
276 
277 ! Compute total condensate
278  if ( nwat<2 ) then
279  do k=1,kbot
280  do i=is,ie
281  qcon(i,k) = 0.
282  enddo
283  enddo
284  elseif ( nwat==2 ) then ! GFS_2015
285  do k=1,kbot
286  do i=is,ie
287  qcon(i,k) = q0(i,k,liq_wat)
288  enddo
289  enddo
290  elseif ( nwat==3 ) then
291  do k=1,kbot
292  do i=is,ie
293  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,ice_wat)
294  enddo
295  enddo
296  elseif ( nwat==4 ) then
297  do k=1,kbot
298  do i=is,ie
299  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat)
300  enddo
301  enddo
302  else
303  do k=1,kbot
304  do i=is,ie
305  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)
306  enddo
307  enddo
308  endif
309 
310  do k=kbot, 2, -1
311  km1 = k-1
312  do i=is,ie
313 ! Richardson number = g*delz * del_theta/theta / (del_u**2 + del_v**2)
314 ! Use exact form for "density temperature"
315  tv1 = t0(i,km1)*(1.+xvir*q0(i,km1,sphum)-qcon(i,km1))
316  tv2 = t0(i,k )*(1.+xvir*q0(i,k ,sphum)-qcon(i,k))
317  pt1 = tv1 / pkz(i,j,km1)
318  pt2 = tv2 / pkz(i,j,k )
319 !
320  ri = (gz(i,km1)-gz(i,k))*(pt1-pt2)/( 0.5*(pt1+pt2)* &
321  ((u0(i,km1)-u0(i,k))**2+(v0(i,km1)-v0(i,k))**2+ustar2) )
322  if ( tv1>t_max .and. tv1>tv2 ) then
323 ! top layer unphysically warm
324  ri = 0.
325  elseif ( tv2<t_min ) then
326  ri = min(ri, 0.2)
327  endif
328 ! Adjustment for K-H instability:
329 ! Compute equivalent mass flux: mc
330 ! Add moist 2-dz instability consideration:
331 !!! ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(500.e2,pm(i,k))/250.e2 )
332  ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(400.e2,pm(i,k))/200.e2 )
333 ! Enhancing mixing at the model top
334  if ( k==2 ) then
335  ri_ref = 4.*ri_ref
336  elseif ( k==3 ) then
337  ri_ref = 2.*ri_ref
338  elseif ( k==4 ) then
339  ri_ref = 1.5*ri_ref
340  endif
341 
342  if ( ri < ri_ref ) then
343  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
344  do iq=1,nq
345  h0 = mc*(q0(i,k,iq)-q0(i,km1,iq))
346  q0(i,km1,iq) = q0(i,km1,iq) + h0/delp(i,j,km1)
347  q0(i,k ,iq) = q0(i,k ,iq) - h0/delp(i,j,k )
348  enddo
349 ! Recompute qcon
350  if ( nwat<2 ) then
351  qcon(i,km1) = 0.
352  elseif ( nwat==2 ) then ! GFS_2015
353  qcon(i,km1) = q0(i,km1,liq_wat)
354  elseif ( nwat==3 ) then ! AM3/AM4
355  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat)
356  elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice
357  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat)
358  else
359  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
360  q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
361  endif
362 ! u:
363  h0 = mc*(u0(i,k)-u0(i,k-1))
364  u0(i,k-1) = u0(i,k-1) + h0/delp(i,j,k-1)
365  u0(i,k ) = u0(i,k ) - h0/delp(i,j,k )
366 ! v:
367  h0 = mc*(v0(i,k)-v0(i,k-1))
368  v0(i,k-1) = v0(i,k-1) + h0/delp(i,j,k-1)
369  v0(i,k ) = v0(i,k ) - h0/delp(i,j,k )
370 
371  if ( hydrostatic ) then
372 ! Static energy
373  h0 = mc*(hd(i,k)-hd(i,k-1))
374  hd(i,k-1) = hd(i,k-1) + h0/delp(i,j,k-1)
375  hd(i,k ) = hd(i,k ) - h0/delp(i,j,k )
376  else
377 ! Total energy
378  h0 = mc*(hd(i,k)-hd(i,k-1))
379  te(i,k-1) = te(i,k-1) + h0/delp(i,j,k-1)
380  te(i,k ) = te(i,k ) - h0/delp(i,j,k )
381 ! w:
382  h0 = mc*(w0(i,k)-w0(i,k-1))
383  w0(i,k-1) = w0(i,k-1) + h0/delp(i,j,k-1)
384  w0(i,k ) = w0(i,k ) - h0/delp(i,j,k )
385  endif
386  endif
387  enddo
388 
389 !--------------
390 ! Retrive Temp:
391 !--------------
392  if ( hydrostatic ) then
393  kk = k
394  do i=is,ie
395  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
396  / ( rk - pe(i,kk,j)/pm(i,kk) )
397  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
398  t0(i,kk) = t0(i,kk) / ( rdgas + rz*q0(i,kk,sphum) )
399  enddo
400  kk = k-1
401  do i=is,ie
402  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
403  / ((rk-pe(i,kk,j)/pm(i,kk))*(rdgas+rz*q0(i,kk,sphum)))
404  enddo
405  else
406 ! Non-hydrostatic under constant volume heating/cooling
407  do kk=k-1,k
408  if ( nwat == 0 ) then
409  do i=is,ie
410  cpm(i) = cp_air
411  cvm(i) = cv_air
412  enddo
413  elseif ( nwat == 1 ) then
414  do i=is,ie
415  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
416  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
417  enddo
418  elseif ( nwat == 2 ) then
419  do i=is,ie
420  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
421  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
422  enddo
423  elseif ( nwat == 3 ) then
424  do i=is,ie
425  q_liq = q0(i,kk,liq_wat)
426  q_sol = q0(i,kk,ice_wat)
427  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
428  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
429  enddo
430  elseif ( nwat == 4 ) then
431  do i=is,ie
432  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
433  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
434  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
435  enddo
436  else
437  do i=is,ie
438  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
439  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
440  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
441  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
442  enddo
443  endif
444 
445  do i=is,ie
446  tv = gz(i,kk) + 0.5*(u0(i,kk)**2+v0(i,kk)**2+w0(i,kk)**2)
447  t0(i,kk) = (te(i,kk)- tv) / cvm(i)
448  hd(i,kk) = cpm(i)*t0(i,kk) + tv
449  enddo
450  enddo
451  endif
452  enddo ! k-loop
453  enddo ! n-loop
454 
455 !--------------------
456  if ( fra < 1. ) then
457  do k=1, kbot
458  do i=is,ie
459  t0(i,k) = ta(i,j,k) + (t0(i,k) - ta(i,j,k))*fra
460  u0(i,k) = ua(i,j,k) + (u0(i,k) - ua(i,j,k))*fra
461  v0(i,k) = va(i,j,k) + (v0(i,k) - va(i,j,k))*fra
462  enddo
463  enddo
464 
465  if ( .not. hydrostatic ) then
466  do k=1,kbot
467  do i=is,ie
468  w0(i,k) = w(i,j,k) + (w0(i,k) - w(i,j,k))*fra
469  enddo
470  enddo
471  endif
472 
473  do iq=1,nq
474  do k=1,kbot
475  do i=is,ie
476  q0(i,k,iq) = qa(i,j,k,iq) + (q0(i,k,iq) - qa(i,j,k,iq))*fra
477  enddo
478  enddo
479  enddo
480  endif
481 
482  do k=1,kbot
483  do i=is,ie
484  u_dt(i,j,k) = rdt*(u0(i,k) - ua(i,j,k))
485  v_dt(i,j,k) = rdt*(v0(i,k) - va(i,j,k))
486  ta(i,j,k) = t0(i,k) ! *** temperature updated ***
487 #if defined(GFS_PHYS) || defined(MAPL_MODE)
488  ua(i,j,k) = u0(i,k)
489  va(i,j,k) = v0(i,k)
490 #endif
491  enddo
492  do iq=1,nq
493  do i=is,ie
494  qa(i,j,k,iq) = q0(i,k,iq)
495  enddo
496  enddo
497  enddo
498 
499  if ( .not. hydrostatic ) then
500  do k=1,kbot
501  do i=is,ie
502  w(i,j,k) = w0(i,k) ! w updated
503  enddo
504  enddo
505  endif
506 
507 1000 continue
508 
509  end subroutine fv_subgrid_z
510 
511 #else
512  subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
513  tau, nwat, delp, pe, peln, pkz, ta, qa, ua, va, &
514  hydrostatic, w, delz, u_dt, v_dt, t_dt, q_dt, k_bot )
515 ! Dry convective adjustment-mixing
516 !-------------------------------------------
517  integer, intent(in):: is, ie, js, je, km, nq, nwat
518  integer, intent(in):: isd, ied, jsd, jed
519  integer, intent(in):: tau ! Relaxation time scale
520  real, intent(in):: dt ! model time step
521  real, intent(in):: pe(is-1:ie+1,km+1,js-1:je+1)
522  real, intent(in):: peln(is :ie, km+1,js :je)
523  real, intent(in):: delp(isd:ied,jsd:jed,km) ! Delta p at each model level
524  real, intent(in):: delz(isd:,jsd:,1:) ! Delta z at each model level
525  real, intent(in):: pkz(is:ie,js:je,km)
526  logical, intent(in):: hydrostatic
527  integer, intent(in), optional:: k_bot
528 !
529  real, intent(inout):: ua(isd:ied,jsd:jed,km)
530  real, intent(inout):: va(isd:ied,jsd:jed,km)
531  real, intent(inout):: w(isd:,jsd:,1:)
532  real, intent(inout):: ta(isd:ied,jsd:jed,km) ! Temperature
533  real, intent(inout):: qa(isd:ied,jsd:jed,km,nq) ! Specific humidity & tracers
534  real, intent(inout):: u_dt(isd:ied,jsd:jed,km)
535  real, intent(inout):: v_dt(isd:ied,jsd:jed,km)
536  real, intent(inout):: t_dt(is:ie,js:je,km)
537  real, intent(inout):: q_dt(is:ie,js:je,km,nq)
538 !---------------------------Local variables-----------------------------
539  real, dimension(is:ie,km):: u0, v0, w0, t0, hd, te, gz, tvm, pm, den
540  real q0(is:ie,km,nq), qcon(is:ie,km)
541  real, dimension(is:ie):: gzh, lcp2, icp2, cvm, cpm, qs
542  real ri_ref, ri, pt1, pt2, ratio, tv, cv, tmp, q_liq, q_sol
543  real tv1, tv2, g2, h0, mc, fra, rk, rz, rdt, tvd, tv_surf
544  real dh, dq, qsw, dqsdt, tcp3
545  integer i, j, k, kk, n, m, iq, km1, im, kbot
546  real, parameter:: ustar2 = 1.e-4
547  real:: cv_air, xvir
548  integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
549 
550  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
551  rk = cp_air/rdgas + 1.
552  cv = cp_air - rdgas
553 
554  g2 = 0.5*grav
555 
556  rdt = 1./ dt
557  im = ie-is+1
558 
559  if ( present(k_bot) ) then
560  if ( k_bot < 3 ) return
561  kbot = k_bot
562  else
563  kbot = km
564  endif
565 
566  sphum = get_tracer_index(model_atmos, 'sphum')
567  if ( nwat == 0 ) then
568  xvir = 0.
569  rz = 0.
570  else
571  xvir = zvir
572  rz = rvgas - rdgas ! rz = zvir * rdgas
573  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
574  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
575  rainwat = get_tracer_index(model_atmos, 'rainwat')
576  snowwat = get_tracer_index(model_atmos, 'snowwat')
577  graupel = get_tracer_index(model_atmos, 'graupel')
578  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
579  endif
580 
581 !------------------------------------------------------------------------
582 ! The nonhydrostatic pressure changes if there is heating (under constant
583 ! volume and mass is locally conserved).
584 !------------------------------------------------------------------------
585  m = 3
586  fra = dt/real(tau)
587 
588 !$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
589 !$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
590 !$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra,cld_amt, &
591 !$OMP u_dt,rdt,v_dt,xvir,nwat) &
592 !$OMP private(kk,lcp2,icp2,tcp3,dh,dq,den,qs,qsw,dqsdt,qcon,q0, &
593 !$OMP t0,u0,v0,w0,h0,pm,gzh,tvm,tmp,cpm,cvm, q_liq,q_sol,&
594 !$OMP tv,gz,hd,te,ratio,pt1,pt2,tv1,tv2,ri_ref, ri,mc,km1)
595  do 1000 j=js,je
596 
597  do iq=1, nq
598  do k=1,kbot
599  do i=is,ie
600  q0(i,k,iq) = qa(i,j,k,iq)
601  enddo
602  enddo
603  enddo
604 
605  do k=1,kbot
606  do i=is,ie
607  t0(i,k) = ta(i,j,k)
608  tvm(i,k) = t0(i,k)*(1.+xvir*q0(i,k,sphum))
609  u0(i,k) = ua(i,j,k)
610  v0(i,k) = va(i,j,k)
611  pm(i,k) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
612  enddo
613  enddo
614 
615  do i=is,ie
616  gzh(i) = 0.
617  enddo
618 
619  if( hydrostatic ) then
620  do k=kbot, 1,-1
621  do i=is,ie
622  tv = rdgas*tvm(i,k)
623  den(i,k) = pm(i,k)/tv
624  gz(i,k) = gzh(i) + tv*(1.-pe(i,k,j)/pm(i,k))
625  hd(i,k) = cp_air*tvm(i,k)+gz(i,k)+0.5*(u0(i,k)**2+v0(i,k)**2)
626  gzh(i) = gzh(i) + tv*(peln(i,k+1,j)-peln(i,k,j))
627  enddo
628  enddo
629  else
630  do k=kbot, 1, -1
631  if ( nwat == 0 ) then
632  do i=is,ie
633  cpm(i) = cp_air
634  cvm(i) = cv_air
635  enddo
636  elseif ( nwat==1 ) then
637  do i=is,ie
638  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
639  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
640  enddo
641  elseif ( nwat==2 ) then ! GFS
642  do i=is,ie
643  q_sol = q0(i,k,liq_wat)
644  cpm(i) = (1.-q0(i,k,sphum))*cp_air + q0(i,k,sphum)*cp_vapor
645  cvm(i) = (1.-q0(i,k,sphum))*cv_air + q0(i,k,sphum)*cv_vap
646  enddo
647  elseif ( nwat==3 ) then
648  do i=is,ie
649  q_liq = q0(i,k,liq_wat)
650  q_sol = q0(i,k,ice_wat)
651  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
652  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
653  enddo
654  elseif ( nwat==4 ) then
655  do i=is,ie
656  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
657  cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq
658  cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq
659  enddo
660  else
661  do i=is,ie
662  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
663  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
664  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
665  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
666  enddo
667  endif
668 
669  do i=is,ie
670  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
671  w0(i,k) = w(i,j,k)
672  gz(i,k) = gzh(i) - g2*delz(i,j,k)
673  tmp = gz(i,k) + 0.5*(u0(i,k)**2+v0(i,k)**2+w0(i,k)**2)
674  hd(i,k) = cpm(i)*t0(i,k) + tmp
675  te(i,k) = cvm(i)*t0(i,k) + tmp
676  gzh(i) = gzh(i) - grav*delz(i,j,k)
677  enddo
678  enddo
679  endif
680 
681  do n=1,m
682 
683  if ( m==3 ) then
684  if ( n==1) ratio = 0.25
685  if ( n==2) ratio = 0.5
686  if ( n==3) ratio = 0.999
687  else
688  ratio = real(n)/real(m)
689  endif
690 
691  do i=is,ie
692  gzh(i) = 0.
693  enddo
694 
695 ! Compute total condensate
696  if ( nwat<2 ) then
697  do k=1,kbot
698  do i=is,ie
699  qcon(i,k) = 0.
700  enddo
701  enddo
702  elseif ( nwat==2 ) then ! GFS_2015
703  do k=1,kbot
704  do i=is,ie
705  qcon(i,k) = q0(i,k,liq_wat)
706  enddo
707  enddo
708  elseif ( nwat==3 ) then
709  do k=1,kbot
710  do i=is,ie
711  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,ice_wat)
712  enddo
713  enddo
714  elseif ( nwat==4 ) then
715  do k=1,kbot
716  do i=is,ie
717  qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat)
718  enddo
719  enddo
720  else
721  do k=1,kbot
722  do i=is,ie
723  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)
724  enddo
725  enddo
726  endif
727 
728  do k=kbot, 2, -1
729  km1 = k-1
730 #ifdef TEST_MQ
731  if(nwat>0) call qsmith(im, 1, 1, t0(is,km1), pm(is,km1), q0(is,km1,sphum), qs)
732 #endif
733  do i=is,ie
734 ! Richardson number = g*delz * del_theta/theta / (del_u**2 + del_v**2)
735 ! Use exact form for "density temperature"
736  tv1 = t0(i,km1)*(1.+xvir*q0(i,km1,sphum)-qcon(i,km1))
737  tv2 = t0(i,k )*(1.+xvir*q0(i,k ,sphum)-qcon(i,k))
738  pt1 = tv1 / pkz(i,j,km1)
739  pt2 = tv2 / pkz(i,j,k )
740  ri = (gz(i,km1)-gz(i,k))*(pt1-pt2)/( 0.5*(pt1+pt2)* &
741  ((u0(i,km1)-u0(i,k))**2+(v0(i,km1)-v0(i,k))**2+ustar2) )
742 
743 ! Adjustment for K-H instability:
744 ! Compute equivalent mass flux: mc
745 ! Add moist 2-dz instability consideration:
746  ri_ref = min(ri_max, ri_min + (ri_max-ri_min)*dim(500.e2,pm(i,k))/250.e2 )
747 #ifdef TEST_MQ
748  if ( nwat > 0 ) then
749  h0 = hd(i,k)-hd(i,km1) + hlv0*(q0(i,k,sphum)-qs(i))
750  if ( h0 > 0. ) ri_ref = 5.*ri_ref
751  endif
752 #endif
753  if ( ri < ri_ref ) then
754  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
755  do iq=1,nq
756  h0 = mc*(q0(i,k,iq)-q0(i,km1,iq))
757  q0(i,km1,iq) = q0(i,km1,iq) + h0/delp(i,j,km1)
758  q0(i,k ,iq) = q0(i,k ,iq) - h0/delp(i,j,k )
759  enddo
760 ! Recompute qcon
761  if ( nwat<2 ) then
762  qcon(i,km1) = 0.
763  elseif ( nwat==2 ) then ! GFS_2015
764  qcon(i,km1) = q0(i,km1,liq_wat)
765  elseif ( nwat==3 ) then ! AM3/AM4
766  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat)
767  elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice
768  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat)
769  else
770  qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
771  q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
772  endif
773 ! u:
774  h0 = mc*(u0(i,k)-u0(i,k-1))
775  u0(i,k-1) = u0(i,k-1) + h0/delp(i,j,k-1)
776  u0(i,k ) = u0(i,k ) - h0/delp(i,j,k )
777 ! v:
778  h0 = mc*(v0(i,k)-v0(i,k-1))
779  v0(i,k-1) = v0(i,k-1) + h0/delp(i,j,k-1)
780  v0(i,k ) = v0(i,k ) - h0/delp(i,j,k )
781 
782  if ( hydrostatic ) then
783 ! Static energy
784  h0 = mc*(hd(i,k)-hd(i,k-1))
785  hd(i,k-1) = hd(i,k-1) + h0/delp(i,j,k-1)
786  hd(i,k ) = hd(i,k ) - h0/delp(i,j,k )
787  else
788 ! Total energy
789  h0 = mc*(hd(i,k)-hd(i,k-1))
790  te(i,k-1) = te(i,k-1) + h0/delp(i,j,k-1)
791  te(i,k ) = te(i,k ) - h0/delp(i,j,k )
792 ! w:
793  h0 = mc*(w0(i,k)-w0(i,k-1))
794  w0(i,k-1) = w0(i,k-1) + h0/delp(i,j,k-1)
795  w0(i,k ) = w0(i,k ) - h0/delp(i,j,k )
796  endif
797  endif
798  enddo
799 
800 !--------------
801 ! Retrive Temp:
802 !--------------
803  if ( hydrostatic ) then
804  kk = k
805  do i=is,ie
806  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
807  / ( rk - pe(i,kk,j)/pm(i,kk) )
808  gzh(i) = gzh(i) + t0(i,kk)*(peln(i,kk+1,j)-peln(i,kk,j))
809  t0(i,kk) = t0(i,kk) / ( rdgas + rz*q0(i,kk,sphum) )
810  enddo
811  kk = k-1
812  do i=is,ie
813  t0(i,kk) = (hd(i,kk)-gzh(i)-0.5*(u0(i,kk)**2+v0(i,kk)**2)) &
814  / ((rk-pe(i,kk,j)/pm(i,kk))*(rdgas+rz*q0(i,kk,sphum)))
815  enddo
816  else
817 ! Non-hydrostatic under constant volume heating/cooling
818  do kk=k-1,k
819  if ( nwat == 0 ) then
820  do i=is,ie
821  cpm(i) = cp_air
822  cvm(i) = cv_air
823  enddo
824  elseif ( nwat == 1 ) then
825  do i=is,ie
826  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
827  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
828  enddo
829  elseif ( nwat == 2 ) then
830  do i=is,ie
831  cpm(i) = (1.-q0(i,kk,sphum))*cp_air + q0(i,kk,sphum)*cp_vapor
832  cvm(i) = (1.-q0(i,kk,sphum))*cv_air + q0(i,kk,sphum)*cv_vap
833  enddo
834  elseif ( nwat == 3 ) then
835  do i=is,ie
836  q_liq = q0(i,kk,liq_wat)
837  q_sol = q0(i,kk,ice_wat)
838  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
839  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
840  enddo
841  elseif ( nwat == 4 ) then
842  do i=is,ie
843  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
844  cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq
845  cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq
846  enddo
847  else
848  do i=is,ie
849  q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
850  q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
851  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
852  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
853  enddo
854  endif
855 
856 
857  do i=is,ie
858  tv = gz(i,kk) + 0.5*(u0(i,kk)**2+v0(i,kk)**2+w0(i,kk)**2)
859  t0(i,kk) = (te(i,kk)- tv) / cvm(i)
860  hd(i,kk) = cpm(i)*t0(i,kk) + tv
861  enddo
862  enddo
863  endif
864  enddo ! k-loop
865  enddo ! n-loop
866 
867 !--------------------
868  if ( fra < 1. ) then
869  do k=1, kbot
870  do i=is,ie
871  t0(i,k) = ta(i,j,k) + (t0(i,k) - ta(i,j,k))*fra
872  u0(i,k) = ua(i,j,k) + (u0(i,k) - ua(i,j,k))*fra
873  v0(i,k) = va(i,j,k) + (v0(i,k) - va(i,j,k))*fra
874  enddo
875  enddo
876 
877  if ( .not. hydrostatic ) then
878  do k=1,kbot
879  do i=is,ie
880  w0(i,k) = w(i,j,k) + (w0(i,k) - w(i,j,k))*fra
881  enddo
882  enddo
883  endif
884 
885  do iq=1,nq
886  do k=1,kbot
887  do i=is,ie
888  q0(i,k,iq) = qa(i,j,k,iq) + (q0(i,k,iq) - qa(i,j,k,iq))*fra
889  enddo
890  enddo
891  enddo
892  endif
893 
894 !----------------------
895 ! Saturation adjustment
896 !----------------------
897 #ifndef GFS_PHYS
898  if ( nwat > 5 ) then
899  do k=1, kbot
900  if ( hydrostatic ) then
901  do i=is, ie
902 ! Compute pressure hydrostatically
903  den(i,k) = pm(i,k)/(rdgas*t0(i,k)*(1.+xvir*q0(i,k,sphum)))
904  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
905  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
906  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
907  lcp2(i) = hlv / cpm(i)
908  icp2(i) = hlf / cpm(i)
909  enddo
910  else
911  do i=is, ie
912  den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
913  q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
914  q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
915  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
916  lcp2(i) = (lv0+dc_vap*t0(i,k)) / cvm(i)
917  icp2(i) = (li0+dc_ice*t0(i,k)) / cvm(i)
918  enddo
919  endif
920 
921 ! Prevent super saturation over water:
922  do i=is, ie
923  qsw = wqs2(t0(i,k), den(i,k), dqsdt)
924  dq = q0(i,k,sphum) - qsw
925  if ( dq > 0. ) then ! remove super-saturation
926  tcp3 = lcp2(i) + icp2(i)*min(1., dim(tice,t0(i,k))/40.)
927  tmp = dq/(1.+tcp3*dqsdt)
928  t0(i,k) = t0(i,k) + tmp*lcp2(i)
929  q0(i,k, sphum) = q0(i,k, sphum) - tmp
930  q0(i,k,liq_wat) = q0(i,k,liq_wat) + tmp
931 ! Grid box mean is saturated; 50% or higher cloud cover
932  if (cld_amt .gt. 0) then
933  qa(i,j,k,cld_amt) = max(0.5, min(1., qa(i,j,k,cld_amt)+25.*tmp/qsw))
934  end if
935  endif
936 ! Freezing
937  tmp = tice-40. - t0(i,k)
938  if( tmp>0.0 .and. q0(i,k,liq_wat)>0. ) then
939  dh = min( q0(i,k,liq_wat), q0(i,k,liq_wat)*tmp*0.125, tmp/icp2(i) )
940  q0(i,k,liq_wat) = q0(i,k,liq_wat) - dh
941  q0(i,k,ice_wat) = q0(i,k,ice_wat) + dh
942  t0(i,k) = t0(i,k) + dh*icp2(i)
943  endif
944  enddo
945  enddo
946  endif
947 #endif
948 
949  do k=1,kbot
950  do i=is,ie
951  u_dt(i,j,k) = rdt*(u0(i,k) - ua(i,j,k))
952  v_dt(i,j,k) = rdt*(v0(i,k) - va(i,j,k))
953  ta(i,j,k) = t0(i,k) ! *** temperature updated ***
954 #if defined(GFS_PHYS) || defined(MAPL_MODE)
955  ua(i,j,k) = u0(i,k)
956  va(i,j,k) = v0(i,k)
957 #endif
958  enddo
959  do iq=1,nq
960  if (iq .ne. cld_amt ) then
961  do i=is,ie
962  qa(i,j,k,iq) = q0(i,k,iq)
963  enddo
964  endif
965  enddo
966  enddo
967 
968  if ( .not. hydrostatic ) then
969  do k=1,kbot
970  do i=is,ie
971  w(i,j,k) = w0(i,k) ! w updated
972  enddo
973  enddo
974  endif
975 
976 1000 continue
977 
978 
979  end subroutine fv_subgrid_z
980 #endif
981 
982 
983  subroutine qsmith_init
984  integer, parameter:: length=2621
985  integer i
986 
987  if( .not. allocated(table) ) then
988 ! Generate es table (dT = 0.1 deg. C)
989 
990  allocate ( table(length) )
991  allocate ( des(length) )
992 
993  call qs_table(length, table)
994 
995  do i=1,length-1
996  des(i) = table(i+1) - table(i)
997  enddo
998  des(length) = des(length-1)
999  endif
1000 
1001  end subroutine qsmith_init
1002 
1003 
1004  subroutine qsmith(im, km, k1, t, p, q, qs, dqdt)
1005 ! input T in deg K; p (Pa)
1006  integer, intent(in):: im, km, k1
1007  real, intent(in),dimension(im,km):: t, p, q
1008  real, intent(out),dimension(im,km):: qs
1009  real, intent(out), optional:: dqdt(im,km)
1010 ! Local:
1011  real es(im,km)
1012  real ap1, eps10
1013  real tmin
1014  integer i, k, it
1015 
1016  tmin = tice-160.
1017  eps10 = 10.*esl
1018 
1019  if( .not. allocated(table) ) call qsmith_init
1020 
1021  do k=k1,km
1022  do i=1,im
1023  ap1 = 10.*dim(t(i,k), tmin) + 1.
1024  ap1 = min(2621., ap1)
1025  it = ap1
1026  es(i,k) = table(it) + (ap1-it)*des(it)
1027  qs(i,k) = esl*es(i,k)*(1.+zvir*q(i,k))/p(i,k)
1028  enddo
1029  enddo
1030 
1031  if ( present(dqdt) ) then
1032  do k=k1,km
1033  do i=1,im
1034  ap1 = 10.*dim(t(i,k), tmin) + 1.
1035  ap1 = min(2621., ap1) - 0.5
1036  it = ap1
1037  dqdt(i,k) = eps10*(des(it)+(ap1-it)*(des(it+1)-des(it)))*(1.+zvir*q(i,k))/p(i,k)
1038  enddo
1039  enddo
1040  endif
1041 
1042  end subroutine qsmith
1043 
1044 
1045  subroutine qs_table(n,table)
1046  integer, intent(in):: n
1047  real table (n)
1048  real:: dt=0.1
1049  real esbasw, tbasw, esbasi, tbasi, Tmin, tem, aa, b, c, d, e, esh20
1050  real wice, wh2o
1051  integer i
1052 ! Constants
1053  esbasw = 1013246.0
1054  tbasw = 373.16
1055  tbasi = 273.16
1056  tmin = tbasi - 160.
1057 ! Compute es over water
1058 ! see smithsonian meteorological tables page 350.
1059  do i=1,n
1060  tem = tmin+dt*real(i-1)
1061  aa = -7.90298*(tbasw/tem-1)
1062  b = 5.02808*alog10(tbasw/tem)
1063  c = -1.3816e-07*(10**((1-tem/tbasw)*11.344)-1)
1064  d = 8.1328e-03*(10**((tbasw/tem-1)*(-3.49149))-1)
1065  e = alog10(esbasw)
1066  table(i) = 0.1*10**(aa+b+c+d+e)
1067  enddo
1068 
1069  end subroutine qs_table
1070 
1071  subroutine qs_table_m(n,table)
1072 ! Mixed (blended) table
1073  integer, intent(in):: n
1074  real table (n)
1075  real esupc(200)
1076  real:: dt=0.1
1077  real esbasw, tbasw, esbasi, tbasi, Tmin, tem, aa, b, c, d, e, esh20
1078  real wice, wh2o
1079  integer i
1080 
1081 ! Constants
1082  esbasw = 1013246.0
1083  tbasw = 373.16
1084  esbasi = 6107.1
1085  tbasi = 273.16
1086 ! ****************************************************
1087 ! Compute es over ice between -160c and 0 c.
1088  tmin = tbasi - 160.
1089 ! see smithsonian meteorological tables page 350.
1090  do i=1,1600
1091  tem = tmin+dt*real(i-1)
1092  aa = -9.09718 *(tbasi/tem-1.0)
1093  b = -3.56654 *alog10(tbasi/tem)
1094  c = 0.876793*(1.0-tem/tbasi)
1095  e = alog10(esbasi)
1096  table(i)=10**(aa+b+c+e)
1097  enddo
1098 ! *****************************************************
1099 ! Compute es over water between -20c and 102c.
1100 ! see smithsonian meteorological tables page 350.
1101  do i=1,1221
1102  tem = 253.16+dt*real(i-1)
1103  aa = -7.90298*(tbasw/tem-1)
1104  b = 5.02808*alog10(tbasw/tem)
1105  c = -1.3816e-07*(10**((1-tem/tbasw)*11.344)-1)
1106  d = 8.1328e-03*(10**((tbasw/tem-1)*(-3.49149))-1)
1107  e = alog10(esbasw)
1108  esh20 = 10**(aa+b+c+d+e)
1109  if (i <= 200) then
1110  esupc(i) = esh20
1111  else
1112  table(i+1400) = esh20
1113  endif
1114  enddo
1115 !********************************************************************
1116 ! Derive blended es over ice and supercooled water between -20c and 0c
1117  do i=1,200
1118  tem = 253.16+dt*real(i-1)
1119  wice = 0.05*(273.16-tem)
1120  wh2o = 0.05*(tem-253.16)
1121  table(i+1400) = wice*table(i+1400)+wh2o*esupc(i)
1122  enddo
1123 
1124  do i=1,n
1125  table(i) = table(i)*0.1
1126  enddo
1127 
1128  end subroutine qs_table_m
1129 
1130  subroutine neg_adj3(is, ie, js, je, ng, kbot, hydrostatic, &
1131  peln, delz, pt, dp, qv, ql, qr, qi, qs, qg, qa, check_negative)
1132 
1133 ! This is designed for 6-class micro-physics schemes
1134  integer, intent(in):: is, ie, js, je, ng, kbot
1135  logical, intent(in):: hydrostatic
1136  real, intent(in):: dp(is-ng:ie+ng,js-ng:je+ng,kbot) ! total delp-p
1137  real, intent(in):: delz(is-ng:,js-ng:,1:)
1138  real, intent(in):: peln(is:ie,kbot+1,js:je) ! ln(pe)
1139  logical, intent(in), OPTIONAL :: check_negative
1140  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: &
1141  pt, qv, ql, qr, qi, qs, qg
1142  real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa
1143 ! Local:
1144  logical:: sat_adj = .false.
1145  real, parameter :: t48 = tice - 48.
1146  real, dimension(is:ie,kbot):: dpk, q2
1147 real, dimension(is:ie,js:je):: pt2, qv2, ql2, qi2, qs2, qr2, qg2, dp2, p2, icpk, lcpk
1148  real:: cv_air
1149  real:: dq, qsum, dq1, q_liq, q_sol, cpm, sink, qsw, dwsdt
1150  integer i, j, k
1151 
1152  cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
1153 
1154  if ( present(check_negative) ) then
1155  if ( check_negative ) then
1156  call prt_negative('Temperature', pt, is, ie, js, je, ng, kbot, 165.)
1157  call prt_negative('sphum', qv, is, ie, js, je, ng, kbot, -1.e-8)
1158  call prt_negative('liq_wat', ql, is, ie, js, je, ng, kbot, -1.e-7)
1159  call prt_negative('rainwat', qr, is, ie, js, je, ng, kbot, -1.e-7)
1160  call prt_negative('ice_wat', qi, is, ie, js, je, ng, kbot, -1.e-7)
1161  call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7)
1162  call prt_negative('graupel', qg, is, ie, js, je, ng, kbot, -1.e-7)
1163  endif
1164  endif
1165 
1166  if ( hydrostatic ) then
1167  d0_vap = cp_vapor - c_liq
1168  else
1169  d0_vap = cv_vap - c_liq
1170  endif
1171  lv00 = hlv0 - d0_vap*t_ice
1172 
1173 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,ql,qi,qs,qr,qg,dp,pt, &
1174 !$OMP lv00, d0_vap,hydrostatic,peln,delz,cv_air,sat_adj) &
1175 !$OMP private(dq,dq1,qsum,dp2,p2,pt2,qv2,ql2,qi2,qs2,qg2,qr2, &
1176 !$OMP lcpk,icpk,qsw,dwsdt,sink,q_liq,q_sol,cpm)
1177  do k=1, kbot
1178  do j=js, je
1179  do i=is, ie
1180  qv2(i,j) = qv(i,j,k)
1181  ql2(i,j) = ql(i,j,k)
1182  qi2(i,j) = qi(i,j,k)
1183  qs2(i,j) = qs(i,j,k)
1184  qr2(i,j) = qr(i,j,k)
1185  qg2(i,j) = qg(i,j,k)
1186  dp2(i,j) = dp(i,j,k)
1187  pt2(i,j) = pt(i,j,k)
1188  enddo
1189  enddo
1190 
1191  if ( hydrostatic ) then
1192  do j=js, je
1193  do i=is, ie
1194  p2(i,j) = dp2(i,j)/(peln(i,k+1,j)-peln(i,k,j))
1195  q_liq = max(0., ql2(i,j) + qr2(i,j))
1196  q_sol = max(0., qi2(i,j) + qs2(i,j))
1197  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cp_air + qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice
1198  lcpk(i,j) = hlv / cpm
1199  icpk(i,j) = hlf / cpm
1200  enddo
1201  enddo
1202  else
1203  do j=js, je
1204  do i=is, ie
1205  p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*(1.+zvir*qv2(i,j))
1206  q_liq = max(0., ql2(i,j) + qr2(i,j))
1207  q_sol = max(0., qi2(i,j) + qs2(i,j))
1208  cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cv_air + qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice
1209  lcpk(i,j) = (lv00+d0_vap*pt2(i,j)) / cpm
1210  icpk(i,j) = (li0+dc_ice*pt2(i,j)) / cpm
1211  enddo
1212  enddo
1213  endif
1214 
1215 ! Fix the negatives:
1216 !-----------
1217 ! Ice-phase:
1218 !-----------
1219  do j=js, je
1220  do i=is, ie
1221  qsum = qi2(i,j) + qs2(i,j)
1222  if ( qsum > 0. ) then
1223  if ( qi2(i,j) < 0. ) then
1224  qi2(i,j) = 0.
1225  qs2(i,j) = qsum
1226  elseif ( qs2(i,j) < 0. ) then
1227  qs2(i,j) = 0.
1228  qi2(i,j) = qsum
1229  endif
1230  else
1231 ! borrow froom graupel
1232  qi2(i,j) = 0.
1233  qs2(i,j) = 0.
1234  qg2(i,j) = qg2(i,j) + qsum
1235  endif
1236 
1237 ! At this stage qi and qs should be positive definite
1238 ! If graupel < 0 then borrow from qs then qi
1239  if ( qg2(i,j) < 0. ) then
1240  dq = min( qs2(i,j), -qg2(i,j) )
1241  qs2(i,j) = qs2(i,j) - dq
1242  qg2(i,j) = qg2(i,j) + dq
1243  if ( qg2(i,j) < 0. ) then
1244 ! if qg is still negative
1245  dq = min( qi2(i,j), -qg2(i,j) )
1246  qi2(i,j) = qi2(i,j) - dq
1247  qg2(i,j) = qg2(i,j) + dq
1248  endif
1249  endif
1250 
1251 ! If qg is still negative then borrow from rain water: phase change
1252  if ( qg2(i,j)<0. .and. qr2(i,j)>0. ) then
1253  dq = min( qr2(i,j), -qg2(i,j) )
1254  qg2(i,j) = qg2(i,j) + dq
1255  qr2(i,j) = qr2(i,j) - dq
1256  pt2(i,j) = pt2(i,j) + dq*icpk(i,j) ! conserve total energy
1257  endif
1258 ! If qg is still negative then borrow from cloud water: phase change
1259  if ( qg2(i,j)<0. .and. ql2(i,j)>0. ) then
1260  dq = min( ql2(i,j), -qg2(i,j) )
1261  qg2(i,j) = qg2(i,j) + dq
1262  ql2(i,j) = ql2(i,j) - dq
1263  pt2(i,j) = pt2(i,j) + dq*icpk(i,j)
1264  endif
1265 ! Last resort; borrow from water vapor
1266  if ( qg2(i,j)<0. .and. qv2(i,j)>0. ) then
1267  dq = min( 0.999*qv2(i,j), -qg2(i,j) )
1268  qg2(i,j) = qg2(i,j) + dq
1269  qv2(i,j) = qv2(i,j) - dq
1270  pt2(i,j) = pt2(i,j) + dq*(icpk(i,j)+lcpk(i,j))
1271  endif
1272 
1273 !--------------
1274 ! Liquid phase:
1275 !--------------
1276  qsum = ql2(i,j) + qr2(i,j)
1277  if ( qsum > 0. ) then
1278  if ( qr2(i,j) < 0. ) then
1279  qr2(i,j) = 0.
1280  ql2(i,j) = qsum
1281  elseif ( ql2(i,j) < 0. ) then
1282  ql2(i,j) = 0.
1283  qr2(i,j) = qsum
1284  endif
1285  else
1286  ql2(i,j) = 0.
1287  qr2(i,j) = qsum ! rain water is still negative
1288 ! fill negative rain with qg first
1289  dq = min( max(0.0, qg2(i,j)), -qr2(i,j) )
1290  qr2(i,j) = qr2(i,j) + dq
1291  qg2(i,j) = qg2(i,j) - dq
1292  pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
1293  if ( qr(i,j,k) < 0. ) then
1294 ! fill negative rain with available qi & qs (cooling)
1295  dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) )
1296  qr2(i,j) = qr2(i,j) + dq
1297  dq1 = min( dq, qs2(i,j) )
1298  qs2(i,j) = qs2(i,j) - dq1
1299  qi2(i,j) = qi2(i,j) + dq1 - dq
1300  pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
1301  endif
1302 ! fix negative rain water with available vapor
1303  if ( qr2(i,j)<0. .and. qv2(i,j)>0. ) then
1304  dq = min( 0.999*qv2(i,j), -qr2(i,j) )
1305  qv2(i,j) = qv2(i,j) - dq
1306  qr2(i,j) = qr2(i,j) + dq
1307  pt2(i,j) = pt2(i,j) + dq*lcpk(i,j)
1308  endif
1309  endif
1310  enddo
1311  enddo
1312 
1313 !******************************************
1314 ! Fast moist physics: Saturation adjustment
1315 !******************************************
1316 #ifndef GFS_PHYS
1317  if ( sat_adj ) then
1318 
1319  do j=js, je
1320  do i=is, ie
1321 ! Melting of cloud ice into cloud water ********
1322  if ( qi2(i,j)>1.e-8 .and. pt2(i,j) > tice ) then
1323  sink = min( qi2(i,j), (pt2(i,j)-tice)/icpk(i,j) )
1324  ql2(i,j) = ql2(i,j) + sink
1325  qi2(i,j) = qi2(i,j) - sink
1326  pt2(i,j) = pt2(i,j) - sink*icpk(i,j)
1327  endif
1328 
1329 ! vapor <---> liquid water --------------------------------
1330  qsw = wqsat2_moist(pt2(i,j), qv2(i,j), p2(i,j), dwsdt)
1331  sink = min( ql2(i,j), (qsw-qv2(i,j))/(1.+lcpk(i,j)*dwsdt) )
1332  qv2(i,j) = qv2(i,j) + sink
1333  ql2(i,j) = ql2(i,j) - sink
1334  pt2(i,j) = pt2(i,j) - sink*lcpk(i,j)
1335 !-----------------------------------------------------------
1336 
1337 ! freezing of cloud water ********
1338  if( ql2(i,j)>1.e-8 .and. pt2(i,j) < t48 ) then
1339 ! Enforce complete freezing below t_00 (-48 C)
1340  sink = min( ql2(i,j), (t48-pt2(i,j))/icpk(i,j) )
1341  ql2(i,j) = ql2(i,j) - sink
1342  qi2(i,j) = qi2(i,j) + sink
1343  pt2(i,j) = pt2(i,j) + sink*icpk(i,j)
1344  endif ! significant ql existed
1345  enddo
1346  enddo
1347  endif
1348 #endif
1349 
1350 !----------------------------------------------------------------
1351 ! Update fields:
1352  do j=js, je
1353  do i=is, ie
1354  qv(i,j,k) = qv2(i,j)
1355  ql(i,j,k) = ql2(i,j)
1356  qi(i,j,k) = qi2(i,j)
1357  qs(i,j,k) = qs2(i,j)
1358  qr(i,j,k) = qr2(i,j)
1359  qg(i,j,k) = qg2(i,j)
1360  pt(i,j,k) = pt2(i,j)
1361  enddo
1362  enddo
1363 
1364  enddo
1365 
1366 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,dp,qg,qr) &
1367 !$OMP private(dpk, q2)
1368  do j=js, je
1369 ! Graupel:
1370  do k=1,kbot
1371  do i=is,ie
1372  dpk(i,k) = dp(i,j,k)
1373  q2(i,k) = qg(i,j,k)
1374  enddo
1375  enddo
1376  call fillq(ie-is+1, kbot, q2, dpk)
1377  do k=1,kbot
1378  do i=is,ie
1379  qg(i,j,k) = q2(i,k)
1380  enddo
1381  enddo
1382 ! Rain water:
1383  do k=1,kbot
1384  do i=is,ie
1385  q2(i,k) = qr(i,j,k)
1386  enddo
1387  enddo
1388  call fillq(ie-is+1, kbot, q2, dpk)
1389  do k=1,kbot
1390  do i=is,ie
1391  qr(i,j,k) = q2(i,k)
1392  enddo
1393  enddo
1394  enddo
1395 
1396 !-----------------------------------
1397 ! Fix water vapor
1398 !-----------------------------------
1399 ! Top layer: borrow from below
1400  k = 1
1401 !$OMP parallel do default(none) shared(is,ie,js,je,k,qv,dp)
1402  do j=js, je
1403  do i=is, ie
1404  if( qv(i,j,k) < 0. ) then
1405  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1406  qv(i,j,k ) = 0.
1407  endif
1408  enddo
1409  enddo
1410 
1411 ! this OpenMP do-loop cannot be parallelized with recursion on k/k-1
1412 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) &
1413 !$OMP private(dq)
1414  do j=js, je
1415  do k=2,kbot-1
1416  do i=is, ie
1417  if( qv(i,j,k) < 0. .and. qv(i,j,k-1) > 0. ) then
1418  dq = min(-qv(i,j,k)*dp(i,j,k), qv(i,j,k-1)*dp(i,j,k-1))
1419  qv(i,j,k-1) = qv(i,j,k-1) - dq/dp(i,j,k-1)
1420  qv(i,j,k ) = qv(i,j,k ) + dq/dp(i,j,k )
1421  endif
1422  if( qv(i,j,k) < 0. ) then
1423  qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1424  qv(i,j,k ) = 0.
1425  endif
1426  enddo
1427  enddo
1428  enddo
1429 
1430 ! Bottom layer; Borrow from above
1431 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) private(dq)
1432  do j=js, je
1433  do i=is, ie
1434  if( qv(i,j,kbot) < 0. ) then
1435  do k=kbot-1,1,-1
1436  if ( qv(i,j,kbot)>=0. ) goto 123
1437  if ( qv(i,j,k) > 0. ) then
1438  dq = min(-qv(i,j,kbot)*dp(i,j,kbot), qv(i,j,k)*dp(i,j,k))
1439  qv(i,j,k ) = qv(i,j,k ) - dq/dp(i,j,k)
1440  qv(i,j,kbot) = qv(i,j,kbot) + dq/dp(i,j,kbot)
1441  endif
1442  enddo ! k-loop
1443 123 continue
1444  endif
1445  enddo ! i-loop
1446  enddo ! j-loop
1447 
1448 
1449  if (present(qa)) then
1450 !-----------------------------------
1451 ! Fix negative cloud fraction
1452 !-----------------------------------
1453 ! this OpenMP do-loop cannot be parallelized by the recursion on k/k+1
1454 !$OMP parallel do default(none) shared(is,ie,js,je,kbot,qa,dp)
1455  do j=js, je
1456  do k=1,kbot-1
1457  do i=is, ie
1458  if( qa(i,j,k) < 0. ) then
1459  qa(i,j,k+1) = qa(i,j,k+1) + qa(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
1460  qa(i,j,k ) = 0.
1461  endif
1462  enddo
1463  enddo
1464  enddo
1465 
1466 ! Bottom layer; Borrow from above
1467 !$OMP parallel do default(none) shared(is,ie,js,je,qa,kbot,dp) &
1468 !$OMP private(dq)
1469  do j=js, je
1470  do i=is, ie
1471  if( qa(i,j,kbot) < 0. .and. qa(i,j,kbot-1)>0.) then
1472  dq = min(-qa(i,j,kbot)*dp(i,j,kbot), qa(i,j,kbot-1)*dp(i,j,kbot-1))
1473  qa(i,j,kbot-1) = qa(i,j,kbot-1) - dq/dp(i,j,kbot-1)
1474  qa(i,j,kbot ) = qa(i,j,kbot ) + dq/dp(i,j,kbot )
1475  endif
1476 ! if qa is still < 0
1477  qa(i,j,kbot) = max(0., qa(i,j,kbot))
1478  enddo
1479  enddo
1480 
1481  endif
1482 
1483  end subroutine neg_adj3
1484 
1485  subroutine fillq(im, km, q, dp)
1486 ! Aggresive 1D filling algorithm for qr and qg
1487  integer, intent(in):: im, km
1488  real, intent(inout), dimension(im,km):: q, dp
1489  integer:: i, k
1490  real:: sum1, sum2, dq
1491 
1492  do 500 i=1,im
1493  sum1 = 0.
1494  do k=1,km
1495  if ( q(i,k)>0. ) then
1496  sum1 = sum1 + q(i,k)*dp(i,k)
1497  endif
1498  enddo
1499  if ( sum1<1.e-12 ) goto 500
1500  sum2 = 0.
1501  do k=km,1,-1
1502  if ( q(i,k)<0.0 .and. sum1>0. ) then
1503  dq = min( sum1, -q(i,k)*dp(i,k) )
1504  sum1 = sum1 - dq
1505  sum2 = sum2 + dq
1506  q(i,k) = q(i,k) + dq/dp(i,k)
1507  endif
1508  enddo
1509  do k=km,1,-1
1510  if ( q(i,k)>0.0 .and. sum2>0. ) then
1511  dq = min( sum2, q(i,k)*dp(i,k) )
1512  sum2 = sum2 - dq
1513  q(i,k) = q(i,k) - dq/dp(i,k)
1514  endif
1515  enddo
1516 500 continue
1517 
1518  end subroutine fillq
1519 
1520  subroutine prt_negative(qname, q, is, ie, js, je, n_g, km, threshold)
1521  character(len=*), intent(in):: qname
1522  integer, intent(in):: is, ie, js, je
1523  integer, intent(in):: n_g, km
1524  real, intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
1525  real, intent(in):: threshold
1526  real qmin
1527  integer i,j,k
1528 
1529  qmin = q(is,js,1)
1530  do k=1,km
1531  do j=js,je
1532  do i=is,ie
1533  qmin = min(qmin, q(i,j,k))
1534  enddo
1535  enddo
1536  enddo
1537  call mp_reduce_min(qmin)
1538  if(is_master() .and. qmin<threshold) write(6,*) qname, ' min (negative) = ', qmin
1539 
1540  end subroutine prt_negative
1541 
1542 
1543 end module fv_sg_nlm_mod
tracer_manager_mod::get_tracer_index
Definition: tracer_manager.F90:124
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
fv_sg_nlm_mod::t2_min
real, parameter t2_min
Definition: fv_sg_nlm.F90:57
constants_mod
Definition: constants.F90:24
lin_cld_microphys_mod::wqs2
real function, public wqs2(ta, den, dqdt)
Definition: lin_cloud_microphys_nlm.F90:3569
field_manager_mod::model_atmos
integer, parameter, public model_atmos
Definition: field_manager.F90:299
field_manager_mod
Definition: field_manager.F90:20
fv_sg_nlm_mod::d0_vap
real d0_vap
Definition: fv_sg_nlm.F90:65
fv_sg_nlm_mod::t1_min
real, parameter t1_min
Definition: fv_sg_nlm.F90:56
constants_mod::hlv
real, parameter, public hlv
Latent heat of evaporation [J/kg].
Definition: constants.F90:80
fv_sg_nlm_mod::t2_max
real, parameter t2_max
Definition: fv_sg_nlm.F90:58
fv_sg_nlm_mod::ri_min
real, parameter ri_min
Definition: fv_sg_nlm.F90:55
fv_sg_nlm_mod::hlf0
real, parameter hlf0
Definition: fv_sg_nlm.F90:50
fv_sg_nlm_mod
Definition: fv_sg_nlm.F90:20
constants_mod::rdgas
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
Definition: constants.F90:77
constants_mod::cp_vapor
real, parameter, public cp_vapor
Specific heat capacity of water vapor at constant pressure [J/kg/deg].
Definition: constants.F90:89
fv_sg_nlm_mod::prt_negative
subroutine prt_negative(qname, q, is, ie, js, je, n_g, km, threshold)
Definition: fv_sg_nlm.F90:1521
fv_sg_nlm_mod::lv00
real lv00
Definition: fv_sg_nlm.F90:65
fv_sg_nlm_mod::table
real, dimension(:), allocatable table
Definition: fv_sg_nlm.F90:64
fv_sg_nlm_mod::dc_ice
real, parameter dc_ice
Definition: fv_sg_nlm.F90:47
lin_cld_microphys_mod
Definition: lin_cloud_microphys_nlm.F90:7
fv_sg_nlm_mod::tice
real, parameter tice
Definition: fv_sg_nlm.F90:37
fv_sg_nlm_mod::zvir
real, parameter zvir
Definition: fv_sg_nlm.F90:63
fv_sg_nlm_mod::qs_table_m
subroutine qs_table_m(n, table)
Definition: fv_sg_nlm.F90:1072
constants_mod::rvgas
real, parameter, public rvgas
Gas constant for water vapor [J/kg/deg].
Definition: constants.F90:78
fv_sg_nlm_mod::cv_vap
real, parameter cv_vap
Definition: fv_sg_nlm.F90:42
fv_sg_nlm_mod::ri_max
real, parameter ri_max
Definition: fv_sg_nlm.F90:54
constants_mod::cp_air
real, parameter, public cp_air
Specific heat capacity of dry air at constant pressure [J/kg/deg].
Definition: constants.F90:83
fv_sg_nlm_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_nlm.F90:515
fv_sg_nlm_mod::qs_table
subroutine qs_table(n, table)
Definition: fv_sg_nlm.F90:1046
fv_sg_nlm_mod::qsmith
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
Definition: fv_sg_nlm.F90:1005
fv_sg_nlm_mod::c_con
real, parameter c_con
Definition: fv_sg_nlm.F90:43
fv_sg_nlm_mod::hlv0
real, parameter hlv0
Definition: fv_sg_nlm.F90:49
constants_mod::hlf
real, parameter, public hlf
Latent heat of fusion [J/kg].
Definition: constants.F90:81
constants_mod::grav
real, parameter, public grav
Acceleration due to gravity [m/s^2].
Definition: constants.F90:76
tracer_manager_mod
Definition: tracer_manager.F90:20
fv_sg_nlm_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_nlm.F90:1132
fv_sg_nlm_mod::esl
real, parameter esl
Definition: fv_sg_nlm.F90:36
fv_sg_nlm_mod::dc_vap
real, parameter dc_vap
Definition: fv_sg_nlm.F90:46
max
#define max(a, b)
Definition: mosaic_util.h:33
fv_sg_nlm_mod::qsmith_init
subroutine qsmith_init
Definition: fv_sg_nlm.F90:984
fv_sg_nlm_mod::t_ice
real, parameter t_ice
Definition: fv_sg_nlm.F90:53
fv_sg_nlm_mod::c_ice
real, parameter c_ice
Definition: fv_sg_nlm.F90:39
fv_sg_nlm_mod::lv0
real, parameter lv0
Definition: fv_sg_nlm.F90:60
min
#define min(a, b)
Definition: mosaic_util.h:32
lin_cld_microphys_mod::wqsat2_moist
real function, public wqsat2_moist(ta, qv, pa, dqdt)
Definition: lin_cloud_microphys_nlm.F90:3666
fv_sg_nlm_mod::t3_max
real, parameter t3_max
Definition: fv_sg_nlm.F90:59
fv_sg_nlm_mod::fillq
subroutine fillq(im, km, q, dp)
Definition: fv_sg_nlm.F90:1486
fv_sg_nlm_mod::des
real, dimension(:), allocatable des
Definition: fv_sg_nlm.F90:64
fv_sg_nlm_mod::li0
real, parameter li0
Definition: fv_sg_nlm.F90:61
constants_mod::kappa
real, parameter, public kappa
RDGAS / CP_AIR [dimensionless].
Definition: constants.F90:82
fv_sg_nlm_mod::c_liq
real, parameter c_liq
Definition: fv_sg_nlm.F90:40
qg
The namespace for the qg model.
Definition: ioda/src/ioda/random_f.cc:12