FV3DYCORE  Version 2.0.0
nh_utils.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 
25 
26 module nh_utils_mod
27 
28 ! Modules Included:
29 ! <table>
30 ! <tr>
31 ! <th>Module Name</th>
32 ! <th>Functions Included</th>
33 ! </tr>
34 ! <tr>
35 ! <td>constants_mod</td>
36 ! <td>rdgas, cp_air, grav</td>
37 ! </tr>
38 ! <tr>
39 ! <td>fv_arrays_mod</td>
40 ! <td>fv_grid_bounds_type, fv_grid_type</td>
41 ! </tr>
42 ! <tr>
43 ! <tr>
44 ! <td>sw_core_mod</td>
45 ! <td>fill_4corners, del6_vt_flux</td>
46 ! </tr>
47 ! <tr>
48 ! <td>tp_core_mod</td>
49 ! <td>fv_tp_2d</td>
50 ! </tr>
51 ! </table>
52 
53  use constants_mod, only: rdgas, cp_air, grav
54  use tp_core_mod, only: fv_tp_2d
55  use sw_core_mod, only: fill_4corners, del6_vt_flux
56  use fv_arrays_mod, only: fv_grid_bounds_type, fv_grid_type,fv_nest_bc_type_3d
57 #ifdef MULTI_GASES
58  use multi_gases_mod, only: vicpqd, vicvqd
59 #endif
60 
61  implicit none
62  private
63 
64  public update_dz_c, update_dz_d, nh_bc
65  public sim_solver, sim1_solver, sim3_solver
66  public sim3p0_solver, rim_2d
67  public riem_solver_c
68 
69  real, parameter:: dz_min = 6.
70  real, parameter:: r3 = 1./3.
71 
72 CONTAINS
73 
74  subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, &
75  npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type)
76 ! !INPUT PARAMETERS:
77  type(fv_grid_bounds_type), intent(IN) :: bd
78  integer, intent(in):: is, ie, js, je, ng, km, npx, npy, grid_type
79  logical, intent(IN):: sw_corner, se_corner, ne_corner, nw_corner
80  real, intent(in):: dt
81  real, intent(in):: dp0(km)
82  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: ut, vt
83  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng):: area
84  real, intent(inout):: gz(is-ng:ie+ng,js-ng:je+ng,km+1)
85  real, intent(in ):: zs(is-ng:ie+ng, js-ng:je+ng)
86  real, intent( out):: ws(is-ng:ie+ng, js-ng:je+ng)
87 ! Local Work array:
88  real:: gz2(is-ng:ie+ng,js-ng:je+ng)
89  real, dimension(is-1:ie+2,js-1:je+1):: xfx, fx
90  real, dimension(is-1:ie+1,js-1:je+2):: yfx, fy
91  real, parameter:: r14 = 1./14.
92  integer i, j, k
93  integer:: is1, ie1, js1, je1
94  integer:: ie2, je2
95  real:: rdt, top_ratio, bot_ratio, int_ratio
96 !--------------------------------------------------------------------
97 
98  rdt = 1. / dt
99 
100  top_ratio = dp0(1 ) / (dp0( 1)+dp0(2 ))
101  bot_ratio = dp0(km) / (dp0(km-1)+dp0(km))
102 
103  is1 = is - 1
104  js1 = js - 1
105 
106  ie1 = ie + 1
107  je1 = je + 1
108 
109  ie2 = ie + 2
110  je2 = je + 2
111 
112 !$OMP parallel do default(none) shared(js1,je1,is1,ie2,km,je2,ie1,ut,top_ratio,vt, &
113 !$OMP bot_ratio,dp0,js,je,ng,is,ie,gz,grid_type, &
114 !$OMP bd,npx,npy,sw_corner,se_corner,ne_corner, &
115 !$OMP nw_corner,area) &
116 !$OMP private(gz2, xfx, yfx, fx, fy, int_ratio)
117  do 6000 k=1,km+1
118 
119  if ( k==1 ) then
120  do j=js1, je1
121  do i=is1, ie2
122  xfx(i,j) = ut(i,j,1) + (ut(i,j,1)-ut(i,j,2))*top_ratio
123  enddo
124  enddo
125  do j=js1, je2
126  do i=is1, ie1
127  yfx(i,j) = vt(i,j,1) + (vt(i,j,1)-vt(i,j,2))*top_ratio
128  enddo
129  enddo
130  elseif ( k==km+1 ) then
131 ! Bottom extrapolation
132  do j=js1, je1
133  do i=is1, ie2
134  xfx(i,j) = ut(i,j,km) + (ut(i,j,km)-ut(i,j,km-1))*bot_ratio
135 ! xfx(i,j) = r14*(3.*ut(i,j,km-2)-13.*ut(i,j,km-1)+24.*ut(i,j,km))
136 ! if ( xfx(i,j)*ut(i,j,km)<0. ) xfx(i,j) = 0.
137  enddo
138  enddo
139  do j=js1, je2
140  do i=is1, ie1
141  yfx(i,j) = vt(i,j,km) + (vt(i,j,km)-vt(i,j,km-1))*bot_ratio
142 ! yfx(i,j) = r14*(3.*vt(i,j,km-2)-13.*vt(i,j,km-1)+24.*vt(i,j,km))
143 ! if ( yfx(i,j)*vt(i,j,km)<0. ) yfx(i,j) = 0.
144  enddo
145  enddo
146  else
147  int_ratio = 1./(dp0(k-1)+dp0(k))
148  do j=js1, je1
149  do i=is1, ie2
150  xfx(i,j) = (dp0(k)*ut(i,j,k-1)+dp0(k-1)*ut(i,j,k))*int_ratio
151  enddo
152  enddo
153  do j=js1, je2
154  do i=is1, ie1
155  yfx(i,j) = (dp0(k)*vt(i,j,k-1)+dp0(k-1)*vt(i,j,k))*int_ratio
156  enddo
157  enddo
158  endif
159 
160  do j=js-ng, je+ng
161  do i=is-ng, ie+ng
162  gz2(i,j) = gz(i,j,k)
163  enddo
164  enddo
165 
166  if (grid_type < 3) call fill_4corners(gz2, 1, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
167  do j=js1, je1
168  do i=is1, ie2
169  if( xfx(i,j) > 0. ) then
170  fx(i,j) = gz2(i-1,j)
171  else
172  fx(i,j) = gz2(i ,j)
173  endif
174  fx(i,j) = xfx(i,j)*fx(i,j)
175  enddo
176  enddo
177 
178  if (grid_type < 3) call fill_4corners(gz2, 2, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
179  do j=js1,je2
180  do i=is1,ie1
181  if( yfx(i,j) > 0. ) then
182  fy(i,j) = gz2(i,j-1)
183  else
184  fy(i,j) = gz2(i,j)
185  endif
186  fy(i,j) = yfx(i,j)*fy(i,j)
187  enddo
188  enddo
189 
190  do j=js1, je1
191  do i=is1,ie1
192  gz(i,j,k) = (gz2(i,j)*area(i,j) + fx(i,j)- fx(i+1,j)+ fy(i,j)- fy(i,j+1)) &
193  / ( area(i,j) + xfx(i,j)-xfx(i+1,j)+yfx(i,j)-yfx(i,j+1))
194  enddo
195  enddo
196 6000 continue
197 
198 ! Enforce monotonicity of height to prevent blowup
199 !$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km)
200  do j=js1, je1
201  do k=2, km+1
202  do i=is1, ie1
203  gz(i,j,k) = min( gz(i,j,k), gz(i,j,k-1) - dz_min )
204  enddo
205  enddo
206  do i=is1, ie1
207  ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
208  enddo
209  enddo
210 
211  end subroutine update_dz_c
212 
213 
214  subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea, &
215  dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac)
216 
217  type(fv_grid_bounds_type), intent(IN) :: bd
218  integer, intent(in):: is, ie, js, je, ng, km, npx, npy
219  integer, intent(in):: hord
220  real, intent(in) :: rdt
221  real, intent(in) :: dp0(km)
222  real, intent(in) :: area(is-ng:ie+ng,js-ng:je+ng)
223  real, intent(in) :: rarea(is-ng:ie+ng,js-ng:je+ng)
224  real, intent(inout):: damp(km+1)
225  integer, intent(inout):: ndif(km+1)
226  real, intent(in ) :: zs(is-ng:ie+ng,js-ng:je+ng)
227  real, intent(inout) :: zh(is-ng:ie+ng,js-ng:je+ng,km+1)
228  real, intent(inout), dimension(is:ie+1,js-ng:je+ng,km):: crx, xfx
229  real, intent(inout), dimension(is-ng:ie+ng,js:je+1,km):: cry, yfx
230  real, intent(out) :: ws(is:ie,js:je)
231  type(fv_grid_type), intent(IN), target :: gridstruct
232  real, intent(in) :: lim_fac
233 !-----------------------------------------------------
234 ! Local array:
235  real, dimension(is: ie+1, js-ng:je+ng,km+1):: crx_adv, xfx_adv
236  real, dimension(is-ng:ie+ng,js: je+1,km+1 ):: cry_adv, yfx_adv
237  real, dimension(is:ie+1,js:je ):: fx
238  real, dimension(is:ie ,js:je+1):: fy
239  real, dimension(is-ng:ie+ng+1,js-ng:je+ng ):: fx2
240  real, dimension(is-ng:ie+ng ,js-ng:je+ng+1):: fy2
241  real, dimension(is-ng:ie+ng ,js-ng:je+ng ):: wk2, z2
242  real:: ra_x(is:ie,js-ng:je+ng)
243  real:: ra_y(is-ng:ie+ng,js:je)
244 !--------------------------------------------------------------------
245  integer i, j, k, isd, ied, jsd, jed
246  logical:: uniform_grid
247 
248  uniform_grid = .false.
249 
250  damp(km+1) = damp(km)
251  ndif(km+1) = ndif(km)
252 
253  isd = is - ng; ied = ie + ng
254  jsd = js - ng; jed = je + ng
255 
256 !$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, &
257 !$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv)
258  do j=jsd,jed
259  call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, &
260  dp0, uniform_grid, 0)
261  if ( j<=je+1 .and. j>=js ) &
262  call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, &
263  dp0, uniform_grid, 0)
264  enddo
265 
266 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,area,xfx_adv,yfx_adv, &
267 !$OMP damp,zh,crx_adv,cry_adv,npx,npy,hord,gridstruct,bd, &
268 !$OMP ndif,rarea,lim_fac) &
269 !$OMP private(z2, fx2, fy2, ra_x, ra_y, fx, fy,wk2)
270  do k=1,km+1
271 
272  do j=jsd,jed
273  do i=is,ie
274  ra_x(i,j) = area(i,j) + xfx_adv(i,j,k) - xfx_adv(i+1,j,k)
275  enddo
276  enddo
277  do j=js,je
278  do i=isd,ied
279  ra_y(i,j) = area(i,j) + yfx_adv(i,j,k) - yfx_adv(i,j+1,k)
280  enddo
281  enddo
282 
283  if ( damp(k)>1.e-5 ) then
284  do j=jsd,jed
285  do i=isd,ied
286  z2(i,j) = zh(i,j,k)
287  enddo
288  enddo
289  call fv_tp_2d(z2, crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, &
290  fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac)
291  call del6_vt_flux(ndif(k), npx, npy, damp(k), z2, wk2, fx2, fy2, gridstruct, bd)
292  do j=js,je
293  do i=is,ie
294  zh(i,j,k) = (z2(i,j)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) &
295  / (ra_x(i,j)+ra_y(i,j)-area(i,j)) + (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*rarea(i,j)
296  enddo
297  enddo
298  else
299  call fv_tp_2d(zh(isd,jsd,k), crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, &
300  fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac)
301  do j=js,je
302  do i=is,ie
303  zh(i,j,k) = (zh(i,j,k)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) &
304  / (ra_x(i,j) + ra_y(i,j) - area(i,j))
305 ! zh(i,j,k) = rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) &
306 ! + zh(i,j,k)*(3.-rarea(i,j)*(ra_x(i,j) + ra_y(i,j)))
307  enddo
308  enddo
309  endif
310 
311  enddo
312 
313 !$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt)
314  do j=js, je
315  do k=2, km+1
316  do i=is, ie
317 ! Enforce monotonicity of height to prevent blowup
318  zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min )
319  enddo
320  enddo
321  do i=is,ie
322  ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
323  enddo
324  enddo
325 
326  end subroutine update_dz_d
327 
328 
329  subroutine riem_solver_c(ms, dt, is, ie, js, je, km, ng, &
330  akap, cappa, cp, &
331 #ifdef MULTI_GASES
332  kapad, &
333 #endif
334  ptop, hs, w3, pt, q_con, &
335  delp, gz, pef, ws, p_fac, a_imp, scale_m)
336 
337  integer, intent(in):: is, ie, js, je, ng, km
338  integer, intent(in):: ms
339  real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m
340  real, intent(in):: ws(is-ng:ie+ng,js-ng:je+ng)
341  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delp
342  real, intent(in), dimension(is-ng:,js-ng:,1:):: q_con, cappa
343 #ifdef MULTI_GASES
344  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: kapad
345 #endif
346  real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng)
347  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: w3
348 ! OUTPUT PARAMETERS
349  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: gz
350  real, intent( out), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: pef
351 ! Local:
352  real, dimension(is-1:ie+1,km ):: dm, dz2, w2, pm2, gm2, cp2
353  real, dimension(is-1:ie+1,km+1):: pem, pe2, peg
354 #ifdef MULTI_GASES
355  real, dimension(is-1:ie+1,km ):: kapad2
356 #endif
357  real gama, rgrav
358  integer i, j, k
359  integer is1, ie1
360 
361  gama = 1./(1.-akap)
362  rgrav = 1./grav
363 
364  is1 = is - 1
365  ie1 = ie + 1
366 
367 !$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, &
368 #ifdef MULTI_GASES
369 !$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad) &
370 !$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg, kapad2)
371 #else
372 !$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) &
373 !$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg)
374 #endif
375  do 2000 j=js-1, je+1
376 
377  do k=1,km
378  do i=is1, ie1
379  dm(i,k) = delp(i,j,k)
380  enddo
381  enddo
382 
383  do i=is1, ie1
384  pef(i,j,1) = ptop ! full pressure at top
385  pem(i,1) = ptop
386 #ifdef USE_COND
387  peg(i,1) = ptop
388 #endif
389  enddo
390 
391  do k=2,km+1
392  do i=is1, ie1
393  pem(i,k) = pem(i,k-1) + dm(i,k-1)
394 #ifdef USE_COND
395 ! Excluding contribution from condensates:
396  peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1))
397 #endif
398  enddo
399  enddo
400 
401  do k=1,km
402  do i=is1, ie1
403  dz2(i,k) = gz(i,j,k+1) - gz(i,j,k)
404 #ifdef USE_COND
405  pm2(i,k) = (peg(i,k+1)-peg(i,k))/log(peg(i,k+1)/peg(i,k))
406 
407 #ifdef MOIST_CAPPA
408  cp2(i,k) = cappa(i,j,k)
409  gm2(i,k) = 1. / (1.-cp2(i,k))
410 #endif
411 
412 #else
413  pm2(i,k) = dm(i,k)/log(pem(i,k+1)/pem(i,k))
414 #endif
415 #ifdef MULTI_GASES
416  kapad2(i,k) = kapad(i,j,k)
417 #endif
418  dm(i,k) = dm(i,k) * rgrav
419  w2(i,k) = w3(i,j,k)
420  enddo
421  enddo
422 
423 
424  if ( a_imp < -0.01 ) then
425  call sim3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, &
426 #ifdef MULTI_GASES
427  kapad2, &
428 #endif
429  pe2, dm, &
430  pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, scale_m)
431  elseif ( a_imp <= 0.5 ) then
432  call rim_2d(ms, dt, is1, ie1, km, rdgas, gama, gm2, &
433 #ifdef MULTI_GASES
434  kapad2, &
435 #endif
436  pe2, &
437  dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.)
438  else
439  call sim1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, &
440 #ifdef MULTI_GASES
441  kapad2, &
442 #endif
443  pe2, &
444  dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac)
445  endif
446 
447 
448  do k=2,km+1
449  do i=is1, ie1
450  pef(i,j,k) = pe2(i,k) + pem(i,k) ! add hydrostatic full-component
451  enddo
452  enddo
453 
454 ! Compute Height * grav (for p-gradient computation)
455  do i=is1, ie1
456  gz(i,j,km+1) = hs(i,j)
457  enddo
458 
459  do k=km,1,-1
460  do i=is1, ie1
461  gz(i,j,k) = gz(i,j,k+1) - dz2(i,k)*grav
462  enddo
463  enddo
464 
465 2000 continue
466 
467  end subroutine riem_solver_c
468 
469 
472  subroutine riem_solver3test(ms, dt, is, ie, js, je, km, ng, &
473  isd, ied, jsd, jed, akap, cappa, cp, &
474 #ifdef MULTI_GASES
475  kapad, &
476 #endif
477  ptop, zs, q_con, w, delz, pt, &
478  delp, zh, pe, ppe, pk3, pk, peln, &
479  ws, scale_m, p_fac, a_imp, &
480  use_logp, last_call, fp_out)
481 !--------------------------------------------
482 ! !OUTPUT PARAMETERS
483 ! Ouput: gz: grav*height at edges
484 ! pe: full hydrostatic pressure
485 ! ppe: non-hydrostatic pressure perturbation
486 !--------------------------------------------
487  integer, intent(in):: ms, is, ie, js, je, km, ng
488  integer, intent(in):: isd, ied, jsd, jed
489  real, intent(in):: dt ! the BIG horizontal Lagrangian time step
490  real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m
491  real, intent(in):: zs(isd:ied,jsd:jed)
492  logical, intent(in):: last_call, use_logp, fp_out
493  real, intent(in):: ws(is:ie,js:je)
494  real, intent(in), dimension(isd:,jsd:,1:):: q_con, cappa
495  real, intent(in), dimension(isd:ied,jsd:jed,km):: delp, pt
496 #ifdef MULTI_GASES
497  real, intent(in), dimension(isd:ied,jsd:jed,km):: kapad
498 #endif
499  real, intent(inout), dimension(isd:ied,jsd:jed,km+1):: zh
500  real, intent(inout), dimension(isd:ied,jsd:jed,km):: w
501  real, intent(inout):: pe(is-1:ie+1,km+1,js-1:je+1)
502  real, intent(out):: peln(is:ie,km+1,js:je) ! ln(pe)
503  real, intent(out), dimension(isd:ied,jsd:jed,km+1):: ppe
504  real, intent(out):: delz(is:ie,js:je,km)
505  real, intent(out):: pk(is:ie,js:je,km+1)
506  real, intent(out):: pk3(isd:ied,jsd:jed,km+1)
507 ! Local:
508  real, dimension(is:ie,km):: dm, dz2, pm2, w2, gm2, cp2
509  real, dimension(is:ie,km+1)::pem, pe2, peln2, peg, pelng
510 #ifdef MULTI_GASES
511  real, dimension(is:ie,km):: kapad2
512 #endif
513  real gama, rgrav, ptk, peln1
514  integer i, j, k
515 
516  gama = 1./(1.-akap)
517  rgrav = 1./grav
518  peln1 = log(ptop)
519  ptk = exp(akap*peln1)
520 
521 !$OMP parallel do default(none) shared(is,ie,js,je,km,delp,ptop,peln1,pk3,ptk,akap,rgrav,zh,pt, &
522 !$OMP w,a_imp,dt,gama,ws,p_fac,scale_m,ms,delz,last_call, &
523 #ifdef MULTI_GASES
524 !$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad ) &
525 !$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2,kapad2)
526 #else
527 !$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con ) &
528 !$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2)
529 #endif
530  do 2000 j=js, je
531 
532  do k=1,km
533  do i=is, ie
534  dm(i,k) = delp(i,j,k)
535 #ifdef MOIST_CAPPA
536  cp2(i,k) = cappa(i,j,k)
537 #endif
538 #ifdef MULTI_GASES
539  kapad2(i,k) = kapad(i,j,k)
540 #endif
541  enddo
542  enddo
543 
544  do i=is,ie
545  pem(i,1) = ptop
546  peln2(i,1) = peln1
547  pk3(i,j,1) = ptk
548 #ifdef USE_COND
549  peg(i,1) = ptop
550  pelng(i,1) = peln1
551 #endif
552  enddo
553  do k=2,km+1
554  do i=is,ie
555  pem(i,k) = pem(i,k-1) + dm(i,k-1)
556  peln2(i,k) = log(pem(i,k))
557 #ifdef USE_COND
558 ! Excluding contribution from condensates:
559 ! peln used during remap; pk3 used only for p_grad
560  peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1))
561  pelng(i,k) = log(peg(i,k))
562 #endif
563  pk3(i,j,k) = exp(akap*peln2(i,k))
564  enddo
565  enddo
566 
567  do k=1,km
568  do i=is, ie
569 #ifdef USE_COND
570  pm2(i,k) = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
571 
572 #ifdef MOIST_CAPPA
573  gm2(i,k) = 1. / (1.-cp2(i,k))
574 #endif
575 
576 #else
577  pm2(i,k) = dm(i,k)/(peln2(i,k+1)-peln2(i,k))
578 #endif
579  dm(i,k) = dm(i,k) * rgrav
580  dz2(i,k) = zh(i,j,k+1) - zh(i,j,k)
581  w2(i,k) = w(i,j,k)
582  enddo
583  enddo
584 
585 
586  if ( a_imp < -0.999 ) then
587  call sim3p0_solver(dt, is, ie, km, rdgas, gama, akap, &
588 #ifdef MULTI_GASES
589  kapad2, &
590 #endif
591  pe2, dm, &
592  pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, scale_m )
593  elseif ( a_imp < -0.5 ) then
594  call sim3_solver(dt, is, ie, km, rdgas, gama, akap, &
595 #ifdef MULTI_GASES
596  kapad2, &
597 #endif
598  pe2, dm, &
599  pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), abs(a_imp), p_fac, scale_m)
600  elseif ( a_imp <= 0.5 ) then
601  call rim_2d(ms, dt, is, ie, km, rdgas, gama, gm2, &
602 #ifdef MULTI_GASES
603  kapad2, &
604 #endif
605  pe2, &
606  dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.)
607  elseif ( a_imp > 0.999 ) then
608  call sim1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
609 #ifdef MULTI_GASES
610  kapad2, &
611 #endif
612  pe2, dm, &
613  pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac)
614  else
615  call sim_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
616 #ifdef MULTI_GASES
617  kapad2, &
618 #endif
619  pe2, dm, &
620  pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), &
621  a_imp, p_fac, scale_m)
622  endif
623 
624 
625  do k=1, km
626  do i=is, ie
627  w(i,j,k) = w2(i,k)
628  delz(i,j,k) = dz2(i,k)
629  enddo
630  enddo
631 
632  if ( last_call ) then
633  do k=1,km+1
634  do i=is,ie
635  peln(i,k,j) = peln2(i,k)
636  pk(i,j,k) = pk3(i,j,k)
637  pe(i,k,j) = pem(i,k)
638  enddo
639  enddo
640  endif
641 
642  if( fp_out ) then
643  do k=1,km+1
644  do i=is, ie
645  ppe(i,j,k) = pe2(i,k) + pem(i,k)
646  enddo
647  enddo
648  else
649  do k=1,km+1
650  do i=is, ie
651  ppe(i,j,k) = pe2(i,k)
652  enddo
653  enddo
654  endif
655 
656  if ( use_logp ) then
657  do k=2,km+1
658  do i=is, ie
659  pk3(i,j,k) = peln2(i,k)
660  enddo
661  enddo
662  endif
663 
664  do i=is, ie
665  zh(i,j,km+1) = zs(i,j)
666  enddo
667  do k=km,1,-1
668  do i=is, ie
669  zh(i,j,k) = zh(i,j,k+1) - dz2(i,k)
670  enddo
671  enddo
672 
673 2000 continue
674 
675  end subroutine riem_solver3test
676 
677  subroutine imp_diff_w(j, is, ie, js, je, ng, km, cd, delz, ws, w, w3)
678  integer, intent(in) :: j, is, ie, js, je, km, ng
679  real, intent(in) :: cd
680  real, intent(in) :: delz(is:ie, km)
681  real, intent(in) :: w(is:ie, km)
682  real, intent(in) :: ws(is:ie)
683  real, intent(out) :: w3(is-ng:ie+ng,js-ng:je+ng,km)
684 ! Local:
685  real, dimension(is:ie,km):: c, gam, dz, wt
686  real:: bet(is:ie)
687  real:: a
688  integer:: i, k
689 
690  do k=2,km
691  do i=is,ie
692  dz(i,k) = 0.5*(delz(i,k-1)+delz(i,k))
693  enddo
694  enddo
695  do k=1,km-1
696  do i=is,ie
697  c(i,k) = -cd/(dz(i,k+1)*delz(i,k))
698  enddo
699  enddo
700 
701 ! model top:
702  do i=is,ie
703  bet(i) = 1. - c(i,1) ! bet(i) = b
704  wt(i,1) = w(i,1) / bet(i)
705  enddo
706 
707 ! Interior:
708  do k=2,km-1
709  do i=is,ie
710  gam(i,k) = c(i,k-1)/bet(i)
711  a = cd/(dz(i,k)*delz(i,k))
712  bet(i) = (1.+a-c(i,k)) + a*gam(i,k)
713  wt(i,k) = (w(i,k) + a*wt(i,k-1)) / bet(i)
714  enddo
715  enddo
716 
717 ! Bottom:
718  do i=is,ie
719  gam(i,km) = c(i,km-1) / bet(i)
720  a = cd/(dz(i,km)*delz(i,km))
721  wt(i,km) = (w(i,km) + 2.*ws(i)*cd/delz(i,km)**2 &
722  + a*wt(i,km-1))/(1. + a + (cd+cd)/delz(i,km)**2 + a*gam(i,km))
723  enddo
724 
725  do k=km-1,1,-1
726  do i=is,ie
727  wt(i,k) = wt(i,k) - gam(i,k+1)*wt(i,k+1)
728  enddo
729  enddo
730 
731  do k=1,km
732  do i=is,ie
733  w3(i,j,k) = wt(i,k)
734  enddo
735  enddo
736 
737  end subroutine imp_diff_w
738 
739 
740  subroutine rim_2d(ms, bdt, is, ie, km, rgas, gama, gm2, &
741 #ifdef MULTI_GASES
742  kapad2, &
743 #endif
744  pe2, dm2, pm2, w2, dz2, pt2, ws, c_core )
745 
746  integer, intent(in):: ms, is, ie, km
747  real, intent(in):: bdt, gama, rgas
748  real, intent(in), dimension(is:ie,km):: dm2, pm2, gm2
749  logical, intent(in):: c_core
750  real, intent(in ) :: pt2(is:ie,km)
751  real, intent(in ) :: ws(is:ie)
752 ! IN/OUT:
753  real, intent(inout):: dz2(is:ie,km)
754  real, intent(inout):: w2(is:ie,km)
755  real, intent(out ):: pe2(is:ie,km+1)
756 #ifdef MULTI_GASES
757  real, intent(inout), dimension(is:ie,km):: kapad2
758 #endif
759 ! Local:
760  real:: ws2(is:ie)
761  real, dimension(km+1):: m_bot, m_top, r_bot, r_top, pe1, pbar, wbar
762  real, dimension(km):: r_hi, r_lo, dz, wm, dm, dts
763  real, dimension(km):: pf1, wc, cm , pp, pt1
764  real:: dt, rdt, grg, z_frac, ptmp1, rden, pf, time_left
765  real:: m_surf
766 #ifdef MULTI_GASES
767  real gamax
768 #endif
769  integer:: i, k, n, ke, kt1, ktop
770  integer:: ks0, ks1
771 
772  grg = gama * rgas
773  rdt = 1. / bdt
774  dt = bdt / real(ms)
775 
776  pbar(:) = 0.
777  wbar(:) = 0.
778 
779  do i=is,ie
780  ws2(i) = 2.*ws(i)
781  enddo
782 
783  do 6000 i=is,ie
784 
785  do k=1,km
786  dz(k) = dz2(i,k)
787  dm(k) = dm2(i,k)
788  wm(k) = w2(i,k)*dm(k)
789  pt1(k) = pt2(i,k)
790  enddo
791 
792  pe1(:) = 0.
793  wbar(km+1) = ws(i)
794 
795  ks0 = 1
796  if ( ms > 1 .and. ms < 8 ) then
797 ! Continuity of (pbar, wbar) is maintained
798  do k=1, km
799  rden = -rgas*dm(k)/dz(k)
800 #ifdef MOIST_CAPPA
801  pf1(k) = exp( gm2(i,k)*log(rden*pt1(k)) )
802 ! dts(k) = -dz(k)/sqrt(gm2(i,k)*rgas*pf1(k)/rden)
803  dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden)
804 #else
805 #ifdef MULTI_GASES
806  gamax = 1./(1.-kapad2(i,k))
807  pf1(k) = exp( gamax*log(rden*pt1(k)) )
808 #else
809  pf1(k) = exp( gama*log(rden*pt1(k)) )
810 #endif
811  dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden)
812 #endif
813  if ( bdt > dts(k) ) then
814  ks0 = k-1
815  goto 222
816  endif
817  enddo
818  ks0 = km
819 222 if ( ks0==1 ) goto 244
820 
821  do k=1, ks0
822  cm(k) = dm(k) / dts(k)
823  wc(k) = wm(k) / dts(k)
824  pp(k) = pf1(k) - pm2(i,k)
825  enddo
826 
827  wbar(1) = (wc(1)+pp(1)) / cm(1)
828  do k=2, ks0
829  wbar(k) = (wc(k-1)+wc(k) + pp(k)-pp(k-1)) / (cm(k-1)+cm(k))
830  pbar(k) = bdt*(cm(k-1)*wbar(k) - wc(k-1) + pp(k-1))
831  pe1(k) = pbar(k)
832  enddo
833 
834  if ( ks0 == km ) then
835  pbar(km+1) = bdt*(cm(km)*wbar(km+1) - wc(km) + pp(km))
836  if ( c_core ) then
837  do k=1,km
838  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
839  enddo
840  else
841  do k=1,km
842  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
843  w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k)
844  enddo
845  endif
846  pe2(i,1) = 0.
847  do k=2,km+1
848  pe2(i,k) = pbar(k)*rdt
849  enddo
850  goto 6000 ! next i
851  else
852  if ( c_core ) then
853  do k=1, ks0-1
854  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
855  enddo
856  else
857  do k=1, ks0-1
858  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
859  w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k)
860  enddo
861  endif
862  pbar(ks0) = pbar(ks0) / real(ms)
863  endif
864  endif
865 244 ks1 = ks0
866 
867  do 5000 n=1, ms
868 
869  do k=ks1, km
870  rden = -rgas*dm(k)/dz(k)
871 #ifdef MOIST_CAPPA
872  pf = exp( gm2(i,k)*log(rden*pt1(k)) )
873 ! dts(k) = -dz(k) / sqrt( gm2(i,k)*rgas*pf/rden )
874  dts(k) = -dz(k) / sqrt( grg*pf/rden )
875 #else
876 #ifdef MULTI_GASES
877  gamax = 1./(1.-kapad2(i,k))
878  pf = exp( gamax*log(rden*pt1(k)) )
879 #else
880  pf = exp( gama*log(rden*pt1(k)) )
881 #endif
882  dts(k) = -dz(k) / sqrt( grg*pf/rden )
883 #endif
884  ptmp1 = dts(k)*(pf - pm2(i,k))
885  r_lo(k) = wm(k) + ptmp1
886  r_hi(k) = wm(k) - ptmp1
887  enddo
888 
889  ktop = ks1
890  do k=ks1, km
891  if( dt > dts(k) ) then
892  ktop = k-1
893  goto 333
894  endif
895  enddo
896  ktop = km
897 333 continue
898 
899  if ( ktop >= ks1 ) then
900  do k=ks1, ktop
901  z_frac = dt/dts(k)
902  r_bot(k ) = z_frac*r_lo(k)
903  r_top(k+1) = z_frac*r_hi(k)
904  m_bot(k ) = z_frac* dm(k)
905  m_top(k+1) = m_bot(k)
906  enddo
907  if ( ktop == km ) goto 666
908  endif
909 
910  do k=ktop+2, km+1
911  m_top(k) = 0.
912  r_top(k) = 0.
913  enddo
914 
915  kt1 = max(1, ktop)
916  do 444 ke=km+1, ktop+2, -1
917  time_left = dt
918  do k=ke-1, kt1, -1
919  if ( time_left > dts(k) ) then
920  time_left = time_left - dts(k)
921  m_top(ke) = m_top(ke) + dm(k)
922  r_top(ke) = r_top(ke) + r_hi(k)
923  else
924  z_frac = time_left/dts(k)
925  m_top(ke) = m_top(ke) + z_frac*dm(k)
926  r_top(ke) = r_top(ke) + z_frac*r_hi(k)
927  go to 444 ! next level
928  endif
929  enddo
930 444 continue
931 
932  do k=ktop+1, km
933  m_bot(k) = 0.
934  r_bot(k) = 0.
935  enddo
936 
937  do 4000 ke=ktop+1, km
938  time_left = dt
939  do k=ke, km
940  if ( time_left > dts(k) ) then
941  time_left = time_left - dts(k)
942  m_bot(ke) = m_bot(ke) + dm(k)
943  r_bot(ke) = r_bot(ke) + r_lo(k)
944  else
945  z_frac = time_left/dts(k)
946  m_bot(ke) = m_bot(ke) + z_frac* dm(k)
947  r_bot(ke) = r_bot(ke) + z_frac*r_lo(k)
948  go to 4000 ! next interface
949  endif
950  enddo
951  m_surf = m_bot(ke)
952  do k=km, kt1, -1
953  if ( time_left > dts(k) ) then
954  time_left = time_left - dts(k)
955  m_bot(ke) = m_bot(ke) + dm(k)
956  r_bot(ke) = r_bot(ke) - r_hi(k)
957  else
958  z_frac = time_left/dts(k)
959  m_bot(ke) = m_bot(ke) + z_frac* dm(k)
960  r_bot(ke) = r_bot(ke) - z_frac*r_hi(k) + (m_bot(ke)-m_surf)*ws2(i)
961  go to 4000 ! next interface
962  endif
963  enddo
964 4000 continue
965 
966 666 if ( ks1==1 ) wbar(1) = r_bot(1) / m_bot(1)
967  do k=ks1+1, km
968  wbar(k) = (r_bot(k)+r_top(k)) / (m_top(k)+m_bot(k))
969  enddo
970 ! pbar here is actually dt*pbar
971  do k=ks1+1, km+1
972  pbar(k) = m_top(k)*wbar(k) - r_top(k)
973  pe1(k) = pe1(k) + pbar(k)
974  enddo
975 
976  if ( n==ms ) then
977  if ( c_core ) then
978  do k=ks1, km
979  dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k))
980  enddo
981  else
982  do k=ks1, km
983  dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k))
984  w2(i,k) = ( wm(k) + pbar(k+1) - pbar(k) ) / dm(k)
985  enddo
986  endif
987  else
988  do k=ks1, km
989  dz(k) = dz(k) + dt*(wbar(k+1)-wbar(k))
990  wm(k) = wm(k) + pbar(k+1) - pbar(k)
991  enddo
992  endif
993 
994 5000 continue
995  pe2(i,1) = 0.
996  do k=2,km+1
997  pe2(i,k) = pe1(k)*rdt
998  enddo
999 
1000 6000 continue ! end i-loop
1001 
1002  end subroutine rim_2d
1003 
1004  subroutine sim3_solver(dt, is, ie, km, rgas, gama, kappa, &
1005 #ifdef MULTI_GASES
1006  kapad2, &
1007 #endif
1008  pe2, dm, &
1009  pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
1010  integer, intent(in):: is, ie, km
1011  real, intent(in):: dt, rgas, gama, kappa, alpha, p_fac, scale_m
1012  real, intent(in ), dimension(is:ie,km):: dm, pt2
1013  real, intent(in ):: ws(is:ie)
1014  real, intent(in ), dimension(is:ie,km+1):: pem
1015  real, intent(out):: pe2(is:ie,km+1)
1016  real, intent(inout), dimension(is:ie,km):: dz2, w2
1017 #ifdef MULTI_GASES
1018  real, intent(inout), dimension(is:ie,km):: kapad2
1019 #endif
1020 ! Local
1021  real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam
1022  real, dimension(is:ie,km+1):: pp
1023  real, dimension(is:ie):: p1, wk1, bet
1024  real beta, t2, t1g, rdt, ra, capa1, r2g, r6g
1025 #ifdef MULTI_GASES
1026  real gamax, capa1x, t1gx
1027 #endif
1028  integer i, k
1029 
1030  beta = 1. - alpha
1031  ra = 1. / alpha
1032  t2 = beta / alpha
1033  t1g = gama * 2.*(alpha*dt)**2
1034  rdt = 1. / dt
1035  capa1 = kappa - 1.
1036  r2g = grav / 2.
1037  r6g = grav / 6.
1038 
1039 
1040  do k=1,km
1041  do i=is, ie
1042  w1(i,k) = w2(i,k)
1043 ! Full pressure at center
1044 #ifdef MULTI_GASES
1045  gamax = 1. / (1.-kapad2(i,k))
1046  aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1047 #else
1048  aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1049 #endif
1050  enddo
1051  enddo
1052 
1053  do k=1,km-1
1054  do i=is, ie
1055  g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction
1056  bb(i,k) = 2.*(1.+g_rat(i,k))
1057  dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1))
1058  enddo
1059  enddo
1060 
1061 ! pe2 is full p at edges
1062  do i=is, ie
1063 ! Top:
1064  bet(i) = bb(i,1)
1065  pe2(i,1) = pem(i,1)
1066  pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i)
1067 ! Bottom:
1068  bb(i,km) = 2.
1069  dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km)
1070  enddo
1071 
1072  do k=2,km
1073  do i=is, ie
1074  gam(i,k) = g_rat(i,k-1) / bet(i)
1075  bet(i) = bb(i,k) - gam(i,k)
1076  pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i)
1077  enddo
1078  enddo
1079 
1080  do k=km, 2, -1
1081  do i=is, ie
1082  pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1)
1083  enddo
1084  enddo
1085 ! done reconstruction of full:
1086 
1087 ! pp is pert. p at edges
1088  do k=1, km+1
1089  do i=is, ie
1090  pp(i,k) = pe2(i,k) - pem(i,k)
1091  enddo
1092  enddo
1093 
1094  do k=2, km
1095  do i=is, ie
1096 #ifdef MULTI_GASES
1097  gamax = 1./(1.-kapad2(i,k))
1098  t1gx = gamax*2.*(alpha*dt)**2
1099  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1100 #else
1101  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1102 #endif
1103  wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k))
1104  aa(i,k) = aa(i,k) - scale_m*dm(i,1)
1105  enddo
1106  enddo
1107  do i=is, ie
1108  bet(i) = dm(i,1) - aa(i,2)
1109  w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2) + wk(i,2)) / bet(i)
1110  enddo
1111  do k=2,km-1
1112  do i=is, ie
1113  gam(i,k) = aa(i,k) / bet(i)
1114  bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1115  w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) &
1116  - aa(i,k)*w2(i,k-1)) / bet(i)
1117  enddo
1118  enddo
1119  do i=is, ie
1120 #ifdef MULTI_GASES
1121  gamax = 1./(1.-kapad2(i,km))
1122  t1gx = gamax*2.*(alpha*dt)**2
1123  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1124 #else
1125  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1126 #endif
1127  gam(i,km) = aa(i,km) / bet(i)
1128  bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1129  w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk(i,km) + &
1130  wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i)
1131  enddo
1132  do k=km-1, 1, -1
1133  do i=is, ie
1134  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1135  enddo
1136  enddo
1137 
1138 ! pe2 is updated perturbation p at edges
1139  do i=is, ie
1140  pe2(i,1) = 0.
1141  enddo
1142  do k=1,km
1143  do i=is, ie
1144  pe2(i,k+1) = pe2(i,k) + ( dm(i,k)*(w2(i,k)-w1(i,k))*rdt &
1145  - beta*(pp(i,k+1)-pp(i,k)) )*ra
1146  enddo
1147  enddo
1148 
1149 ! Full non-hydro pressure at edges:
1150  do i=is, ie
1151  pe2(i,1) = pem(i,1)
1152  enddo
1153  do k=2,km+1
1154  do i=is, ie
1155  pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k))
1156  enddo
1157  enddo
1158 
1159  do i=is, ie
1160 ! Recover cell-averaged pressure
1161  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km)
1162 #ifdef MULTI_GASES
1163  capa1x = kapad2(i,km) - 1.
1164  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) )
1165 #else
1166  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) )
1167 #endif
1168  enddo
1169 
1170  do k=km-1, 1, -1
1171  do i=is, ie
1172  p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i)
1173 #ifdef MULTI_GASES
1174  capa1x = kapad2(i,k) - 1.
1175  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) )
1176 #else
1177  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) )
1178 #endif
1179  enddo
1180  enddo
1181 
1182  do k=1,km+1
1183  do i=is, ie
1184  pe2(i,k) = pe2(i,k) - pem(i,k)
1185  pe2(i,k) = pe2(i,k) + beta*(pp(i,k) - pe2(i,k))
1186  enddo
1187  enddo
1188 
1189  end subroutine sim3_solver
1190 
1191  subroutine sim3p0_solver(dt, is, ie, km, rgas, gama, kappa, &
1192 #ifdef MULTI_GASES
1193  kapad2, &
1194 #endif
1195  pe2, dm, &
1196  pem, w2, dz2, pt2, ws, p_fac, scale_m)
1197 ! Sa SIM3, but for beta==0
1198  integer, intent(in):: is, ie, km
1199  real, intent(in):: dt, rgas, gama, kappa, p_fac, scale_m
1200  real, intent(in ), dimension(is:ie,km):: dm, pt2
1201  real, intent(in ):: ws(is:ie)
1202  real, intent(in ):: pem(is:ie,km+1)
1203  real, intent(out):: pe2(is:ie,km+1)
1204  real, intent(inout), dimension(is:ie,km):: dz2, w2
1205 #ifdef MULTI_GASES
1206  real, intent(inout), dimension(is:ie,km):: kapad2
1207 #endif
1208 ! Local
1209  real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam
1210  real, dimension(is:ie,km+1):: pp
1211  real, dimension(is:ie):: p1, wk1, bet
1212  real t1g, rdt, capa1, r2g, r6g
1213 #ifdef MULTI_GASES
1214  real gamax, capa1x, t1gx
1215 #endif
1216  integer i, k
1217 
1218  t1g = 2.*gama*dt**2
1219  rdt = 1. / dt
1220  capa1 = kappa - 1.
1221  r2g = grav / 2.
1222  r6g = grav / 6.
1223 
1224  do k=1,km
1225  do i=is, ie
1226  w1(i,k) = w2(i,k)
1227 ! Full pressure at center
1228 #ifdef MULTI_GASES
1229  gamax = 1. / ( 1. - kapad2(i,k) )
1230  aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1231 #else
1232  aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1233 #endif
1234  enddo
1235  enddo
1236 
1237  do k=1,km-1
1238  do i=is, ie
1239  g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction
1240  bb(i,k) = 2.*(1.+g_rat(i,k))
1241  dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1))
1242  enddo
1243  enddo
1244 
1245 ! pe2 is full p at edges
1246  do i=is, ie
1247 ! Top:
1248  bet(i) = bb(i,1)
1249  pe2(i,1) = pem(i,1)
1250  pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i)
1251 ! Bottom:
1252  bb(i,km) = 2.
1253  dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km)
1254  enddo
1255 
1256  do k=2,km
1257  do i=is, ie
1258  gam(i,k) = g_rat(i,k-1) / bet(i)
1259  bet(i) = bb(i,k) - gam(i,k)
1260  pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i)
1261  enddo
1262  enddo
1263 
1264  do k=km, 2, -1
1265  do i=is, ie
1266  pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1)
1267  enddo
1268  enddo
1269 ! done reconstruction of full:
1270 
1271 ! pp is pert. p at edges
1272  do k=1, km+1
1273  do i=is, ie
1274  pp(i,k) = pe2(i,k) - pem(i,k)
1275  enddo
1276  enddo
1277 
1278  do k=2, km
1279  do i=is, ie
1280 #ifdef MULTI_GASES
1281  gamax = 1. / (1.-kapad2(i,k))
1282  t1gx = 2.*gamax*dt**2
1283  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1)
1284 #else
1285  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1)
1286 #endif
1287  enddo
1288  enddo
1289  do i=is, ie
1290  bet(i) = dm(i,1) - aa(i,2)
1291  w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2)) / bet(i)
1292  enddo
1293  do k=2,km-1
1294  do i=is, ie
1295  gam(i,k) = aa(i,k) / bet(i)
1296  bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1297  w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1))/bet(i)
1298  enddo
1299  enddo
1300  do i=is, ie
1301 #ifdef MULTI_GASES
1302  gamax = 1. / (1.-kapad2(i,km))
1303  t1gx = 2.*gamax*dt**2
1304  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1305 #else
1306  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1307 #endif
1308  gam(i,km) = aa(i,km) / bet(i)
1309  bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1310  w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk1(i)*ws(i) - &
1311  aa(i,km)*w2(i,km-1)) / bet(i)
1312  enddo
1313  do k=km-1, 1, -1
1314  do i=is, ie
1315  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1316  enddo
1317  enddo
1318 
1319 ! pe2 is updated perturbation p at edges
1320  do i=is, ie
1321  pe2(i,1) = 0.
1322  enddo
1323  do k=1,km
1324  do i=is, ie
1325  pe2(i,k+1) = pe2(i,k) + dm(i,k)*(w2(i,k)-w1(i,k))*rdt
1326  enddo
1327  enddo
1328 
1329 ! Full non-hydro pressure at edges:
1330  do i=is, ie
1331  pe2(i,1) = pem(i,1)
1332  enddo
1333  do k=2,km+1
1334  do i=is, ie
1335  pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k))
1336  enddo
1337  enddo
1338 
1339  do i=is, ie
1340 ! Recover cell-averaged pressure
1341  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km)
1342 #ifdef MULTI_GASES
1343  capa1x = kapad2(i,km) - 1.
1344  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) )
1345 #else
1346  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) )
1347 #endif
1348  enddo
1349 
1350  do k=km-1, 1, -1
1351  do i=is, ie
1352  p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3-g_rat(i,k)*p1(i)
1353 #ifdef MULTI_GASES
1354  capa1x = kapad2(i,k) - 1.
1355  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) )
1356 #else
1357  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) )
1358 #endif
1359  enddo
1360  enddo
1361 
1362  do k=1,km+1
1363  do i=is, ie
1364  pe2(i,k) = pe2(i,k) - pem(i,k)
1365  enddo
1366  enddo
1367 
1368  end subroutine sim3p0_solver
1369 
1370 
1371  subroutine sim1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
1372 #ifdef MULTI_GASES
1373  kapad2, &
1374 #endif
1375  pe, dm2, &
1376  pm2, pem, w2, dz2, pt2, ws, p_fac)
1377  integer, intent(in):: is, ie, km
1378  real, intent(in):: dt, rgas, gama, kappa, p_fac
1379  real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
1380  real, intent(in ):: ws(is:ie)
1381  real, intent(in ), dimension(is:ie,km+1):: pem
1382  real, intent(out):: pe(is:ie,km+1)
1383  real, intent(inout), dimension(is:ie,km):: dz2, w2
1384 #ifdef MULTI_GASES
1385  real, intent(inout), dimension(is:ie,km):: kapad2
1386 #endif
1387 ! Local
1388  real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam
1389  real, dimension(is:ie,km+1):: pp
1390  real, dimension(is:ie):: p1, bet
1391  real t1g, rdt, capa1
1392 #ifdef MULTI_GASES
1393  real gamax, capa1x, t1gx
1394 #endif
1395  integer i, k
1396 
1397 #ifdef MOIST_CAPPA
1398  t1g = 2.*dt*dt
1399 #else
1400  t1g = gama * 2.*dt*dt
1401 #endif
1402  rdt = 1. / dt
1403  capa1 = kappa - 1.
1404 
1405  do k=1,km
1406  do i=is, ie
1407 #ifdef MOIST_CAPPA
1408  pe(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1409 #else
1410 #ifdef MULTI_GASES
1411  gamax = 1. / ( 1. - kapad2(i,k) )
1412  pe(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1413 #else
1414  pe(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1415 #endif
1416 #endif
1417  w1(i,k) = w2(i,k)
1418  enddo
1419  enddo
1420 
1421  do k=1,km-1
1422  do i=is, ie
1423  g_rat(i,k) = dm2(i,k)/dm2(i,k+1)
1424  bb(i,k) = 2.*(1.+g_rat(i,k))
1425  dd(i,k) = 3.*(pe(i,k) + g_rat(i,k)*pe(i,k+1))
1426  enddo
1427  enddo
1428 
1429  do i=is, ie
1430  bet(i) = bb(i,1)
1431  pp(i,1) = 0.
1432  pp(i,2) = dd(i,1) / bet(i)
1433  bb(i,km) = 2.
1434  dd(i,km) = 3.*pe(i,km)
1435  enddo
1436 
1437  do k=2,km
1438  do i=is, ie
1439  gam(i,k) = g_rat(i,k-1) / bet(i)
1440  bet(i) = bb(i,k) - gam(i,k)
1441  pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i)
1442  enddo
1443  enddo
1444 
1445  do k=km, 2, -1
1446  do i=is, ie
1447  pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1)
1448  enddo
1449  enddo
1450 
1451 ! Start the w-solver
1452  do k=2, km
1453  do i=is, ie
1454 #ifdef MOIST_CAPPA
1455  aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
1456 #else
1457 #ifdef MULTI_GASES
1458  gamax = 1./(1.-kapad2(i,k))
1459  t1gx = gamax * 2.*dt*dt
1460  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
1461 #else
1462  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
1463 #endif
1464 #endif
1465  enddo
1466  enddo
1467  do i=is, ie
1468  bet(i) = dm2(i,1) - aa(i,2)
1469  w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2)) / bet(i)
1470  enddo
1471  do k=2,km-1
1472  do i=is, ie
1473  gam(i,k) = aa(i,k) / bet(i)
1474  bet(i) = dm2(i,k) - (aa(i,k) + aa(i,k+1) + aa(i,k)*gam(i,k))
1475  w2(i,k) = (dm2(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k))-aa(i,k)*w2(i,k-1)) / bet(i)
1476  enddo
1477  enddo
1478  do i=is, ie
1479 #ifdef MOIST_CAPPA
1480  p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1481 #else
1482 #ifdef MULTI_GASES
1483  gamax = 1./(1.-kapad2(i,km))
1484  t1gx = gamax * 2.*dt*dt
1485  p1(i) = t1gx/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1486 #else
1487  p1(i) = t1g/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1488 #endif
1489 #endif
1490  gam(i,km) = aa(i,km) / bet(i)
1491  bet(i) = dm2(i,km) - (aa(i,km)+p1(i) + aa(i,km)*gam(i,km))
1492  w2(i,km) = (dm2(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-p1(i)*ws(i)-aa(i,km)*w2(i,km-1))/bet(i)
1493  enddo
1494  do k=km-1, 1, -1
1495  do i=is, ie
1496  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1497  enddo
1498  enddo
1499 
1500  do i=is, ie
1501  pe(i,1) = 0.
1502  enddo
1503  do k=1,km
1504  do i=is, ie
1505  pe(i,k+1) = pe(i,k) + dm2(i,k)*(w2(i,k)-w1(i,k))*rdt
1506  enddo
1507  enddo
1508 
1509  do i=is, ie
1510  p1(i) = ( pe(i,km) + 2.*pe(i,km+1) )*r3
1511 #ifdef MOIST_CAPPA
1512  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1513 #else
1514 #ifdef MULTI_GASES
1515  capa1x = kapad2(i,km)-1.
1516  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1517 #else
1518  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1519 #endif
1520 #endif
1521  enddo
1522 
1523  do k=km-1, 1, -1
1524  do i=is, ie
1525  p1(i) = (pe(i,k) + bb(i,k)*pe(i,k+1) + g_rat(i,k)*pe(i,k+2))*r3 - g_rat(i,k)*p1(i)
1526 #ifdef MOIST_CAPPA
1527  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1528 
1529 #else
1530 #ifdef MULTI_GASES
1531  capa1x = kapad2(i,k)-1.
1532  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1533 #else
1534  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1535 #endif
1536 #endif
1537  enddo
1538  enddo
1539 
1540  end subroutine sim1_solver
1541 
1542  subroutine sim_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
1543 #ifdef MULTI_GASES
1544  kapad2, &
1545 #endif
1546  pe2, dm2, &
1547  pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
1548  integer, intent(in):: is, ie, km
1549  real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m
1550  real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
1551  real, intent(in ):: ws(is:ie)
1552  real, intent(in ), dimension(is:ie,km+1):: pem
1553  real, intent(out):: pe2(is:ie,km+1)
1554  real, intent(inout), dimension(is:ie,km):: dz2, w2
1555 #ifdef MULTI_GASES
1556  real, intent(inout), dimension(is:ie,km):: kapad2
1557 #endif
1558 ! Local
1559  real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam
1560  real, dimension(is:ie,km+1):: pp
1561  real, dimension(is:ie):: p1, wk1, bet
1562  real beta, t2, t1g, rdt, ra, capa1
1563 #ifdef MULTI_GASES
1564  real gamax, capa1x, t1gx
1565 #endif
1566  integer i, k
1567 
1568  beta = 1. - alpha
1569  ra = 1. / alpha
1570  t2 = beta / alpha
1571 #ifdef MOIST_CAPPA
1572  t1g = 2.*(alpha*dt)**2
1573 #else
1574  t1g = 2.*gama*(alpha*dt)**2
1575 #endif
1576  rdt = 1. / dt
1577  capa1 = kappa - 1.
1578 
1579  do k=1,km
1580  do i=is, ie
1581  w1(i,k) = w2(i,k)
1582 ! P_g perturbation
1583 #ifdef MOIST_CAPPA
1584  pe2(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1585 #else
1586 #ifdef MULTI_GASES
1587  gamax = 1./(1.-kapad2(i,k))
1588  pe2(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1589 #else
1590  pe2(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1591 #endif
1592 #endif
1593  enddo
1594  enddo
1595 
1596  do k=1,km-1
1597  do i=is, ie
1598  g_rat(i,k) = dm2(i,k)/dm2(i,k+1)
1599  bb(i,k) = 2.*(1.+g_rat(i,k))
1600  dd(i,k) = 3.*(pe2(i,k) + g_rat(i,k)*pe2(i,k+1))
1601  enddo
1602  enddo
1603 
1604  do i=is, ie
1605  bet(i) = bb(i,1)
1606  pp(i,1) = 0.
1607  pp(i,2) = dd(i,1) / bet(i)
1608  bb(i,km) = 2.
1609  dd(i,km) = 3.*pe2(i,km)
1610  enddo
1611 
1612  do k=2,km
1613  do i=is, ie
1614  gam(i,k) = g_rat(i,k-1) / bet(i)
1615  bet(i) = bb(i,k) - gam(i,k)
1616  pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i)
1617  enddo
1618  enddo
1619 
1620  do k=km, 2, -1
1621  do i=is, ie
1622  pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1)
1623  enddo
1624  enddo
1625 
1626  do k=1, km+1
1627  do i=is, ie
1628 ! pe2 is Full p
1629  pe2(i,k) = pem(i,k) + pp(i,k)
1630  enddo
1631  enddo
1632 
1633  do k=2, km
1634  do i=is, ie
1635 #ifdef MOIST_CAPPA
1636  aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1637 #else
1638 #ifdef MULTI_GASES
1639  gamax = 1./(1.-kapad2(i,k))
1640  t1gx = 2.*gamax*(alpha*dt)**2
1641  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1642 #else
1643  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1644 #endif
1645 #endif
1646  wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k))
1647  aa(i,k) = aa(i,k) - scale_m*dm2(i,1)
1648  enddo
1649  enddo
1650 ! Top:
1651  do i=is, ie
1652  bet(i) = dm2(i,1) - aa(i,2)
1653  w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2) + wk(i,2)) / bet(i)
1654  enddo
1655 ! Interior:
1656  do k=2,km-1
1657  do i=is, ie
1658  gam(i,k) = aa(i,k) / bet(i)
1659  bet(i) = dm2(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1660  w2(i,k) = (dm2(i,k)*w1(i,k) + dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) &
1661  - aa(i,k)*w2(i,k-1)) / bet(i)
1662  enddo
1663  enddo
1664 ! Bottom: k=km
1665  do i=is, ie
1666 #ifdef MOIST_CAPPA
1667  wk1(i) = t1g*gm2(i,km)/dz2(i,km)*pe2(i,km+1)
1668 #else
1669 #ifdef MULTI_GASES
1670  gamax = 1./(1.-kapad2(i,km))
1671  t1gx = 2.*gamax*(alpha*dt)**2
1672  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1673 #else
1674  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1675 #endif
1676 #endif
1677  gam(i,km) = aa(i,km) / bet(i)
1678  bet(i) = dm2(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1679  w2(i,km) = (dm2(i,km)*w1(i,km) + dt*(pp(i,km+1)-pp(i,km)) - wk(i,km) + &
1680  wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i)
1681  enddo
1682  do k=km-1, 1, -1
1683  do i=is, ie
1684  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1685  enddo
1686  enddo
1687 
1688  do i=is, ie
1689  pe2(i,1) = 0.
1690  enddo
1691  do k=1,km
1692  do i=is, ie
1693  pe2(i,k+1) = pe2(i,k) + ( dm2(i,k)*(w2(i,k)-w1(i,k))*rdt &
1694  - beta*(pp(i,k+1)-pp(i,k)) ) * ra
1695  enddo
1696  enddo
1697 
1698  do i=is, ie
1699  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3
1700 #ifdef MOIST_CAPPA
1701  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp((cp2(i,km)-1.)*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1702 #else
1703 #ifdef MULTI_GASES
1704  capa1x = kapad2(i,km)-1.
1705  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1706 #else
1707  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1708 #endif
1709 #endif
1710  enddo
1711 
1712  do k=km-1, 1, -1
1713  do i=is, ie
1714  p1(i) = (pe2(i,k)+bb(i,k)*pe2(i,k+1)+g_rat(i,k)*pe2(i,k+2))*r3 - g_rat(i,k)*p1(i)
1715 ! delz = -dm*R*T_m / p_gas
1716 #ifdef MOIST_CAPPA
1717  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp((cp2(i,k)-1.)*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1718 #else
1719 #ifdef MULTI_GASES
1720  capa1x = kapad2(i,k)-1.
1721  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1722 #else
1723  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1724 #endif
1725 #endif
1726  enddo
1727  enddo
1728 
1729  do k=1, km+1
1730  do i=is, ie
1731  pe2(i,k) = pe2(i,k) + beta*(pp(i,k)-pe2(i,k))
1732  enddo
1733  enddo
1734 
1735  end subroutine sim_solver
1736 
1737 
1738  subroutine edge_scalar(q1, qe, i1, i2, km, id)
1739 ! Optimized for wind profile reconstruction:
1740  integer, intent(in):: i1, i2, km
1741  integer, intent(in):: id ! 0: pp 1: wind
1742  real, intent(in ), dimension(i1:i2,km):: q1
1743  real, intent(out), dimension(i1:i2,km+1):: qe
1744 !-----------------------------------------------------------------------
1745  real, parameter:: r2o3 = 2./3.
1746  real, parameter:: r4o3 = 4./3.
1747  real gak(km)
1748  real bet
1749  integer i, k
1750 
1751 !------------------------------------------------
1752 ! Optimized coding for uniform grid: SJL Apr 2007
1753 !------------------------------------------------
1754 
1755  if ( id==1 ) then
1756  do i=i1,i2
1757  qe(i,1) = r4o3*q1(i,1) + r2o3*q1(i,2)
1758  enddo
1759  else
1760  do i=i1,i2
1761  qe(i,1) = 1.e30
1762  enddo
1763  endif
1764 
1765  gak(1) = 7./3.
1766  do k=2,km
1767  gak(k) = 1. / (4. - gak(k-1))
1768  do i=i1,i2
1769  qe(i,k) = (3.*(q1(i,k-1) + q1(i,k)) - qe(i,k-1)) * gak(k)
1770  enddo
1771  enddo
1772 
1773  bet = 1. / (1.5 - 3.5*gak(km))
1774  do i=i1,i2
1775  qe(i,km+1) = (4.*q1(i,km) + q1(i,km-1) - 3.5*qe(i,km)) * bet
1776  enddo
1777 
1778  do k=km,1,-1
1779  do i=i1,i2
1780  qe(i,k) = qe(i,k) - gak(k)*qe(i,k+1)
1781  enddo
1782  enddo
1783 
1784  end subroutine edge_scalar
1785 
1786 
1787 
1788  subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter)
1789 ! Optimized for wind profile reconstruction:
1790  integer, intent(in):: i1, i2, j1, j2
1791  integer, intent(in):: j, km
1792  integer, intent(in):: limiter
1793  logical, intent(in):: uniform_grid
1794  real, intent(in):: dp0(km)
1795  real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2
1796  real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e
1797 !-----------------------------------------------------------------------
1798  real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values
1799  real gak(km)
1800  real bet, r2o3, r4o3
1801  real g0, gk, xt1, xt2, a_bot
1802  integer i, k
1803 
1804  if ( uniform_grid ) then
1805 !------------------------------------------------
1806 ! Optimized coding for uniform grid: SJL Apr 2007
1807 !------------------------------------------------
1808  r2o3 = 2./3.
1809  r4o3 = 4./3.
1810  do i=i1,i2
1811  qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2)
1812  qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2)
1813  enddo
1814 
1815  gak(1) = 7./3.
1816  do k=2,km
1817  gak(k) = 1. / (4. - gak(k-1))
1818  do i=i1,i2
1819  qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k)
1820  qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k)
1821  enddo
1822  enddo
1823 
1824  bet = 1. / (1.5 - 3.5*gak(km))
1825  do i=i1,i2
1826  qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet
1827  qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet
1828  enddo
1829 
1830  do k=km,1,-1
1831  do i=i1,i2
1832  qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1)
1833  qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1)
1834  enddo
1835  enddo
1836  else
1837 ! Assuming grid varying in vertical only
1838  g0 = dp0(2) / dp0(1)
1839  xt1 = 2.*g0*(g0+1. )
1840  bet = g0*(g0+0.5)
1841  do i=i1,i2
1842  qe1(i,1) = ( xt1*q1(i,j,1) + q1(i,j,2) ) / bet
1843  qe2(i,1) = ( xt1*q2(i,j,1) + q2(i,j,2) ) / bet
1844  gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet
1845  enddo
1846 
1847  do k=2,km
1848  gk = dp0(k-1) / dp0(k)
1849  do i=i1,i2
1850  bet = 2. + 2.*gk - gam(i,k-1)
1851  qe1(i,k) = ( 3.*(q1(i,j,k-1)+gk*q1(i,j,k)) - qe1(i,k-1) ) / bet
1852  qe2(i,k) = ( 3.*(q2(i,j,k-1)+gk*q2(i,j,k)) - qe2(i,k-1) ) / bet
1853  gam(i,k) = gk / bet
1854  enddo
1855  enddo
1856 
1857  a_bot = 1. + gk*(gk+1.5)
1858  xt1 = 2.*gk*(gk+1.)
1859  do i=i1,i2
1860  xt2 = gk*(gk+0.5) - a_bot*gam(i,km)
1861  qe1(i,km+1) = ( xt1*q1(i,j,km) + q1(i,j,km-1) - a_bot*qe1(i,km) ) / xt2
1862  qe2(i,km+1) = ( xt1*q2(i,j,km) + q2(i,j,km-1) - a_bot*qe2(i,km) ) / xt2
1863  enddo
1864 
1865  do k=km,1,-1
1866  do i=i1,i2
1867  qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1)
1868  qe2(i,k) = qe2(i,k) - gam(i,k)*qe2(i,k+1)
1869  enddo
1870  enddo
1871  endif
1872 
1873 !------------------
1874 ! Apply constraints
1875 !------------------
1876  if ( limiter/=0 ) then ! limit the top & bottom winds
1877  do i=i1,i2
1878 ! Top
1879  if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0.
1880  if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0.
1881 ! Surface:
1882  if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0.
1883  if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0.
1884  enddo
1885  endif
1886 
1887  do k=1,km+1
1888  do i=i1,i2
1889  q1e(i,j,k) = qe1(i,k)
1890  q2e(i,j,k) = qe2(i,k)
1891  enddo
1892  enddo
1893 
1894  end subroutine edge_profile
1895 
1896 !TODO LMH 25may18: do not need delz defined on full compute domain; pass appropriate BCs instead
1897  subroutine nh_bc(ptop, grav, kappa, cp, delp, delzBC, pt, phis, &
1898 #ifdef MULTI_GASES
1899  q , &
1900 #endif
1901 #ifdef USE_COND
1902  q_con, &
1903 #ifdef MOIST_CAPPA
1904  cappa, &
1905 #endif
1906 #endif
1907  pkc, gz, pk3, &
1908  BC_step, BC_split, &
1909  npx, npy, npz, bounded_domain, pkc_pertn, computepk3, fullhalo, bd)
1910 
1911  !INPUT: delp, delz (BC), pt
1912  !OUTPUT: gz, pkc, pk3 (optional)
1913  integer, intent(IN) :: npx, npy, npz
1914  logical, intent(IN) :: pkc_pertn, computepk3, fullhalo, bounded_domain
1915  real, intent(IN) :: ptop, kappa, cp, grav, BC_step, BC_split
1916  type(fv_grid_bounds_type), intent(IN) :: bd
1917  real, intent(IN) :: phis(bd%isd:bd%ied,bd%jsd:bd%jed)
1918  real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: pt, delp
1919  type(fv_nest_bc_type_3d), intent(IN) :: delzBC
1920 #ifdef MULTI_GASES
1921  real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,*):: q
1922 #endif
1923 #ifdef USE_COND
1924  real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: q_con
1925 #ifdef MOIST_CAPPA
1926  real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: cappa
1927 #endif
1928 #endif
1929  real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1):: gz, pkc, pk3
1930 
1931  integer :: i,j,k
1932  real :: gama !'gamma'
1933  real :: ptk, rgrav, rkap, peln1, rdg
1934 
1935  integer :: istart, iend
1936 
1937  integer :: is, ie, js, je
1938  integer :: isd, ied, jsd, jed
1939 
1940  is = bd%is
1941  ie = bd%ie
1942  js = bd%js
1943  je = bd%je
1944  isd = bd%isd
1945  ied = bd%ied
1946  jsd = bd%jsd
1947  jed = bd%jed
1948 
1949  if (.not. bounded_domain) return
1950 
1951  if (is == 1) then
1952 
1953  call nh_bc_k(ptop, grav, kappa, cp, delp, delzbc%west_t0, delzbc%west_t1, pt, phis, &
1954 #ifdef MULTI_GASES
1955  q , &
1956 #endif
1957 #ifdef USE_COND
1958  q_con, &
1959 #ifdef MOIST_CAPPA
1960  cappa, &
1961 #endif
1962 #endif
1963  pkc, gz, pk3, &
1964  bc_step, bc_split, &
1965  pkc_pertn, computepk3, isd, ied, isd, 0, isd, 0, jsd, jed, jsd, jed, npz)
1966 
1967  endif
1968 
1969  if (ie == npx-1) then
1970 
1971  call nh_bc_k(ptop, grav, kappa, cp, delp, delzbc%east_t0, delzbc%east_t1, pt, phis, &
1972 #ifdef MULTI_GASES
1973  q , &
1974 #endif
1975 #ifdef USE_COND
1976  q_con, &
1977 #ifdef MOIST_CAPPA
1978  cappa, &
1979 #endif
1980 #endif
1981  pkc, gz, pk3, &
1982  bc_step, bc_split, &
1983  pkc_pertn, computepk3, isd, ied, npx, ied, npx, ied, jsd, jed, jsd, jed, npz)
1984 
1985  endif
1986 
1987  if (is == 1) then
1988  istart = is
1989  else
1990  istart = isd
1991  end if
1992  if (ie == npx-1) then
1993  iend = ie
1994  else
1995  iend = ied
1996  end if
1997 
1998  if (js == 1) then
1999 
2000  call nh_bc_k(ptop, grav, kappa, cp, delp, delzbc%south_t0, delzbc%south_t1, pt, phis, &
2001 #ifdef MULTI_GASES
2002  q , &
2003 #endif
2004 #ifdef USE_COND
2005  q_con, &
2006 #ifdef MOIST_CAPPA
2007  cappa, &
2008 #endif
2009 #endif
2010  pkc, gz, pk3, &
2011  bc_step, bc_split, &
2012  pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, jsd, 0, npz)
2013 
2014  end if
2015 
2016  if (je == npy-1) then
2017 
2018  call nh_bc_k(ptop, grav, kappa, cp, delp, delzbc%north_t0, delzbc%north_t1, pt, phis, &
2019 #ifdef MULTI_GASES
2020  q , &
2021 #endif
2022 #ifdef USE_COND
2023  q_con, &
2024 #ifdef MOIST_CAPPA
2025  cappa, &
2026 #endif
2027 #endif
2028  pkc, gz, pk3, &
2029  bc_step, bc_split, &
2030  pkc_pertn, computepk3, isd, ied, isd, ied, istart, iend, jsd, jed, npy, jed, npz)
2031  endif
2032 
2033 end subroutine nh_bc
2034 
2035 subroutine nh_bc_k(ptop, grav, kappa, cp, delp, delzBC_t0, delzBC_t1, pt, phis, &
2036 #ifdef MULTI_GASES
2037  q , &
2038 #endif
2039 #ifdef USE_COND
2040  q_con, &
2041 #ifdef MOIST_CAPPA
2042  cappa, &
2043 #endif
2044 #endif
2045  pkc, gz, pk3, &
2046  BC_step, BC_split, &
2047  pkc_pertn, computepk3, isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz)
2048 
2049  integer, intent(IN) :: isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz
2050  real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: delzBC_t0, delzBC_t1
2051  real, intent(IN) :: BC_step, BC_split
2052 
2053  logical, intent(IN) :: pkc_pertn, computepk3
2054  real, intent(IN) :: ptop, kappa, cp, grav
2055  real, intent(IN) :: phis(isd:ied,jsd:jed)
2056  real, intent(IN), dimension(isd:ied,jsd:jed,npz):: pt, delp
2057 #ifdef MULTI_GASES
2058  real, intent(IN), dimension(isd:ied,jsd:jed,npz,*):: q
2059 #endif
2060 #ifdef USE_COND
2061  real, intent(IN), dimension(isd:ied,jsd:jed,npz):: q_con
2062 #ifdef MOIST_CAPPA
2063  real, intent(INOUT), dimension(isd:ied,jsd:jed,npz):: cappa
2064 #endif
2065 #endif
2066  real, intent(INOUT), dimension(isd:ied,jsd:jed,npz+1):: gz, pkc, pk3
2067 
2068  integer :: i,j,k
2069  real :: gama !'gamma'
2070  real :: ptk, rgrav, rkap, peln1, rdg, denom
2071 
2072  real, dimension(istart:iend, npz+1, jstart:jend ) :: pe, peln
2073 #ifdef USE_COND
2074  real, dimension(istart:iend, npz+1 ) :: peg, pelng
2075 #endif
2076  real, dimension(istart:iend, npz) :: gam, bb, dd, pkz
2077  real, dimension(istart:iend, npz-1) :: g_rat
2078  real, dimension(istart:iend) :: bet
2079  real :: pm, delz_int
2080 
2081 
2082  real :: pealn, pebln, rpkz
2083 #ifdef MULTI_GASES
2084  real gamax
2085 #endif
2086  rgrav = 1./grav
2087  gama = 1./(1.-kappa)
2088  ptk = ptop ** kappa
2089  rkap = 1./kappa
2090  peln1 = log(ptop)
2091  rdg = - rdgas * rgrav
2092  denom = 1./bc_split
2093 
2094  do j=jstart,jend
2095 
2096  !GZ
2097  do i=istart,iend
2098  gz(i,j,npz+1) = phis(i,j)
2099  enddo
2100  do k=npz,1,-1
2101  do i=istart,iend
2102  delz_int = (delzbc_t0(i,j,k)*(bc_split-bc_step) + bc_step*delzbc_t1(i,j,k))*denom
2103  gz(i,j,k) = gz(i,j,k+1) - delz_int*grav
2104  enddo
2105  enddo
2106 
2107  !Hydrostatic interface pressure
2108  do i=istart,iend
2109  pe(i,1,j) = ptop
2110  peln(i,1,j) = peln1
2111 #ifdef USE_COND
2112  peg(i,1) = ptop
2113  pelng(i,1) = peln1
2114 #endif
2115  enddo
2116  do k=2,npz+1
2117  do i=istart,iend
2118  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2119  peln(i,k,j) = log(pe(i,k,j))
2120 #ifdef USE_COND
2121  peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1))
2122  pelng(i,k) = log(peg(i,k))
2123 #endif
2124  enddo
2125  enddo
2126 
2127  !Perturbation nonhydro layer-mean pressure (NOT to the kappa)
2128  do k=1,npz
2129  do i=istart,iend
2130  delz_int = (delzbc_t0(i,j,k)*(bc_split-bc_step) + bc_step*delzbc_t1(i,j,k))*denom
2131 
2132  !Full p
2133 #ifdef MOIST_CAPPA
2134  pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz_int*pt(i,j,k)))
2135 #else
2136 #ifdef MULTI_GASES
2137  gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:)))
2138  pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k)))
2139 #else
2140  pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz_int*rdgas*pt(i,j,k)))
2141 #endif
2142 #endif
2143  !hydro
2144 #ifdef USE_COND
2145  pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
2146 #else
2147  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2148 #endif
2149  !Remove hydro cell-mean pressure
2150  pkz(i,k) = pkz(i,k) - pm
2151  enddo
2152  enddo
2153 
2154  !pressure solver
2155  do k=1,npz-1
2156  do i=istart,iend
2157  g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1)
2158  bb(i,k) = 2.*(1. + g_rat(i,k))
2159  dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1))
2160  enddo
2161  enddo
2162 
2163  do i=istart,iend
2164  bet(i) = bb(i,1)
2165  pkc(i,j,1) = 0.
2166  pkc(i,j,2) = dd(i,1)/bet(i)
2167  bb(i,npz) = 2.
2168  dd(i,npz) = 3.*pkz(i,npz)
2169  enddo
2170  do k=2,npz
2171  do i=istart,iend
2172  gam(i,k) = g_rat(i,k-1)/bet(i)
2173  bet(i) = bb(i,k) - gam(i,k)
2174  pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i)
2175  enddo
2176  enddo
2177  do k=npz,2,-1
2178  do i=istart,iend
2179  pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1)
2180 #ifdef NHNEST_DEBUG
2181  if (abs(pkc(i,j,k)) > 1.e5) then
2182  print*, mpp_pe(), i,j,k, 'PKC: ', pkc(i,j,k)
2183  endif
2184 #endif
2185  enddo
2186  enddo
2187 
2188 
2189  enddo
2190 
2191  do j=jstart,jend
2192 
2193  if (.not. pkc_pertn) then
2194  do k=npz+1,1,-1
2195  do i=istart,iend
2196  pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j)
2197  enddo
2198  enddo
2199  endif
2200 
2201  !pk3 if necessary; doesn't require condenstate loading calculation
2202  if (computepk3) then
2203  do i=istart,iend
2204  pk3(i,j,1) = ptk
2205  enddo
2206  do k=2,npz+1
2207  do i=istart,iend
2208  pk3(i,j,k) = exp(kappa*log(pe(i,k,j)))
2209  enddo
2210  enddo
2211  endif
2212 
2213  enddo
2214 
2215 end subroutine nh_bc_k
2216 
2217 
2218 end module nh_utils_mod
nh_utils_mod::sim_solver
subroutine, public sim_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe2, dm2, pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
Definition: nh_utils.F90:1548
multi_gases_mod::vicpqd
pure real function, public vicpqd(q)
Definition: multi_gases.F90:281
nh_utils_mod::sim3_solver
subroutine, public sim3_solver(dt, is, ie, km, rgas, gama, kappa, pe2, dm, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
Definition: nh_utils.F90:1010
fv_arrays_mod::fv_nest_bc_type_3d
Definition: fv_arrays.F90:1043
tp_core_mod::fv_tp_2d
subroutine, public fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, lim_fac, mfx, mfy, mass, nord, damp_c)
The subroutine &#39;fv_tp_2d&#39; contains the FV advection scheme .
Definition: tp_core.F90:109
fv_arrays_mod::fv_grid_type
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:123
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1187
sw_core_mod
The module &#39;sw_core&#39; advances the forward step of the Lagrangian dynamics as described by ...
Definition: sw_core.F90:26
multi_gases_mod
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
nh_utils_mod::r3
real, parameter r3
Definition: nh_utils.F90:70
nh_utils_mod::riem_solver3test
subroutine riem_solver3test(ms, dt, is, ie, js, je, km, ng, isd, ied, jsd, jed, akap, cappa, cp, ptop, zs, q_con, w, delz, pt, delp, zh, pe, ppe, pk3, pk, peln, ws, scale_m, p_fac, a_imp, use_logp, last_call, fp_out)
GFDL - This routine will not give absoulte reproducibility when compiled with -fast-transcendentals. GFDL - It is now inside of nh_core.F90 and being compiled without -fast-transcendentals.
Definition: nh_utils.F90:481
nh_utils_mod::edge_profile
subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter)
Definition: nh_utils.F90:1789
sw_core_mod::del6_vt_flux
subroutine, public del6_vt_flux(nord, npx, npy, damp, q, d2, fx2, fy2, gridstruct, bd)
The subroutine &#39;del6_vt_flux&#39; applies 2nd, 4th, or 6th-order damping to fluxes ("vorticity damping") ...
Definition: sw_core.F90:1581
nh_utils_mod::nh_bc
subroutine, public nh_bc(ptop, grav, kappa, cp, delp, delzBC, pt, phis, pkc, gz, pk3, BC_step, BC_split, npx, npy, npz, bounded_domain, pkc_pertn, computepk3, fullhalo, bd)
Definition: nh_utils.F90:1910
nh_utils_mod
The module &#39;nh_utils&#39; peforms non-hydrostatic computations.
Definition: nh_utils.F90:26
tp_core_mod
The module &#39;tp_core&#39; is a collection of routines to support FV transport.
Definition: tp_core.F90:24
nh_utils_mod::update_dz_c
subroutine, public update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type)
Definition: nh_utils.F90:76
nh_utils_mod::nh_bc_k
subroutine nh_bc_k(ptop, grav, kappa, cp, delp, delzBC_t0, delzBC_t1, pt, phis, pkc, gz, pk3, BC_step, BC_split, pkc_pertn, computepk3, isd, ied, isd_BC, ied_BC, istart, iend, jsd, jed, jstart, jend, npz)
Definition: nh_utils.F90:2048
sw_core_mod::fill_4corners
subroutine, public fill_4corners(q, dir, bd, npx, npy, sw_corner, se_corner, ne_corner, nw_corner)
The subroutine &#39;fill_4corners&#39; fills the 4 corners of the scalar fields only as needed by c_core...
Definition: sw_core.F90:3313
nh_utils_mod::imp_diff_w
subroutine imp_diff_w(j, is, ie, js, je, ng, km, cd, delz, ws, w, w3)
Definition: nh_utils.F90:678
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
nh_utils_mod::sim1_solver
subroutine, public sim1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2, pm2, pem, w2, dz2, pt2, ws, p_fac)
Definition: nh_utils.F90:1377
nh_utils_mod::sim3p0_solver
subroutine, public sim3p0_solver(dt, is, ie, km, rgas, gama, kappa, pe2, dm, pem, w2, dz2, pt2, ws, p_fac, scale_m)
Definition: nh_utils.F90:1197
nh_utils_mod::rim_2d
subroutine, public rim_2d(ms, bdt, is, ie, km, rgas, gama, gm2, pe2, dm2, pm2, w2, dz2, pt2, ws, c_core)
Definition: nh_utils.F90:745
nh_utils_mod::edge_scalar
subroutine edge_scalar(q1, qe, i1, i2, km, id)
Definition: nh_utils.F90:1739
multi_gases_mod::vicvqd
pure real function, public vicvqd(q)
Definition: multi_gases.F90:322
nh_utils_mod::update_dz_d
subroutine, public update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea, dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac)
Definition: nh_utils.F90:216
nh_utils_mod::dz_min
real, parameter dz_min
Definition: nh_utils.F90:69
nh_utils_mod::riem_solver_c
subroutine, public riem_solver_c(ms, dt, is, ie, js, je, km, ng, akap, cappa, cp, ptop, hs, w3, pt, q_con, delp, gz, pef, ws, p_fac, a_imp, scale_m)
Definition: nh_utils.F90:336