FV3DYCORE  Version1.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
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, nest_halo_nh
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 = 2.
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 i=is1, ie1
202  ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
203  enddo
204  do k=km, 1, -1
205  do i=is1, ie1
206  gz(i,j,k) = max( gz(i,j,k), gz(i,j,k+1) + dz_min )
207  enddo
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, delz, ws, rdt, gridstruct, bd, lim_fac, regional)
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( out) ::delz(is-ng:ie+ng,js-ng:je+ng,km)
229  real, intent(inout), dimension(is:ie+1,js-ng:je+ng,km):: crx, xfx
230  real, intent(inout), dimension(is-ng:ie+ng,js:je+1,km):: cry, yfx
231  real, intent(out) :: ws(is:ie,js:je)
232  type(fv_grid_type), intent(IN), target :: gridstruct
233  real, intent(in) :: lim_fac
234  logical,intent(in) :: regional
235 !-----------------------------------------------------
236 ! Local array:
237  real, dimension(is: ie+1, js-ng:je+ng,km+1):: crx_adv, xfx_adv
238  real, dimension(is-ng:ie+ng,js: je+1,km+1 ):: cry_adv, yfx_adv
239  real, dimension(is:ie+1,js:je ):: fx
240  real, dimension(is:ie ,js:je+1):: fy
241  real, dimension(is-ng:ie+ng+1,js-ng:je+ng ):: fx2
242  real, dimension(is-ng:ie+ng ,js-ng:je+ng+1):: fy2
243  real, dimension(is-ng:ie+ng ,js-ng:je+ng ):: wk2, z2
244  real:: ra_x(is:ie,js-ng:je+ng)
245  real:: ra_y(is-ng:ie+ng,js:je)
246 !--------------------------------------------------------------------
247  integer i, j, k, isd, ied, jsd, jed
248  logical:: uniform_grid
249 
250  uniform_grid = .false.
251 
252  damp(km+1) = damp(km)
253  ndif(km+1) = ndif(km)
254 
255  isd = is - ng; ied = ie + ng
256  jsd = js - ng; jed = je + ng
257 
258 !$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, &
259 !$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv)
260  do j=jsd,jed
261  call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, &
262  dp0, uniform_grid, 0)
263  if ( j<=je+1 .and. j>=js ) &
264  call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, &
265  dp0, uniform_grid, 0)
266  enddo
267 
268 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,area,xfx_adv,yfx_adv, &
269 !$OMP damp,zh,crx_adv,cry_adv,npx,npy,hord,gridstruct,bd, &
270 !$OMP ndif,rarea,lim_fac,regional) &
271 !$OMP private(z2, fx2, fy2, ra_x, ra_y, fx, fy,wk2)
272  do k=1,km+1
273 
274  do j=jsd,jed
275  do i=is,ie
276  ra_x(i,j) = area(i,j) + xfx_adv(i,j,k) - xfx_adv(i+1,j,k)
277  enddo
278  enddo
279  do j=js,je
280  do i=isd,ied
281  ra_y(i,j) = area(i,j) + yfx_adv(i,j,k) - yfx_adv(i,j+1,k)
282  enddo
283  enddo
284 
285  if ( damp(k)>1.e-5 ) then
286  do j=jsd,jed
287  do i=isd,ied
288  z2(i,j) = zh(i,j,k)
289  enddo
290  enddo
291  call fv_tp_2d(z2, crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, &
292  fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac, regional)
293  call del6_vt_flux(ndif(k), npx, npy, damp(k), z2, wk2, fx2, fy2, gridstruct, bd)
294  do j=js,je
295  do i=is,ie
296  zh(i,j,k) = (z2(i,j)*area(i,j)+fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) &
297  / (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)
298  enddo
299  enddo
300  else
301  call fv_tp_2d(zh(isd,jsd,k), crx_adv(is,jsd,k), cry_adv(isd,js,k), npx, npy, hord, &
302  fx, fy, xfx_adv(is,jsd,k), yfx_adv(isd,js,k), gridstruct, bd, ra_x, ra_y, lim_fac, regional)
303  do j=js,je
304  do i=is,ie
305  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)) &
306  / (ra_x(i,j) + ra_y(i,j) - area(i,j))
307 ! zh(i,j,k) = rarea(i,j)*(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1)) &
308 ! + zh(i,j,k)*(3.-rarea(i,j)*(ra_x(i,j) + ra_y(i,j)))
309  enddo
310  enddo
311  endif
312 
313  enddo
314 
315 !$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt)
316  do j=js, je
317  do i=is,ie
318  ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
319  enddo
320  do k=km, 1, -1
321  do i=is, ie
322 ! Enforce monotonicity of height to prevent blowup
323  zh(i,j,k) = max( zh(i,j,k), zh(i,j,k+1) + dz_min )
324  enddo
325  enddo
326  enddo
327 
328  end subroutine update_dz_d
329 
330 
331  subroutine riem_solver_c(ms, dt, is, ie, js, je, km, ng, &
332  akap, cappa, cp, &
333 #ifdef MULTI_GASES
334  kapad, &
335 #endif
336  ptop, hs, w3, pt, q_con, &
337  delp, gz, pef, ws, p_fac, a_imp, scale_m)
338 
339  integer, intent(in):: is, ie, js, je, ng, km
340  integer, intent(in):: ms
341  real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m
342  real, intent(in):: ws(is-ng:ie+ng,js-ng:je+ng)
343  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delp
344  real, intent(in), dimension(is-ng:,js-ng:,1:):: q_con, cappa
345 #ifdef MULTI_GASES
346  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: kapad
347 #endif
348  real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng)
349  real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: w3
350 ! OUTPUT PARAMETERS
351  real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: gz
352  real, intent( out), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: pef
353 ! Local:
354  real, dimension(is-1:ie+1,km ):: dm, dz2, w2, pm2, gm2, cp2
355  real, dimension(is-1:ie+1,km+1):: pem, pe2, peg
356 #ifdef MULTI_GASES
357  real, dimension(is-1:ie+1,km ):: kapad2
358 #endif
359  real gama, rgrav
360  integer i, j, k
361  integer is1, ie1
362 
363  gama = 1./(1.-akap)
364  rgrav = 1./grav
365 
366  is1 = is - 1
367  ie1 = ie + 1
368 
369 !$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, &
370 #ifdef MULTI_GASES
371 !$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad) &
372 !$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg, kapad2)
373 #else
374 !$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) &
375 !$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg)
376 #endif
377  do 2000 j=js-1, je+1
378 
379  do k=1,km
380  do i=is1, ie1
381  dm(i,k) = delp(i,j,k)
382  enddo
383  enddo
384 
385  do i=is1, ie1
386  pef(i,j,1) = ptop ! full pressure at top
387  pem(i,1) = ptop
388 #ifdef USE_COND
389  peg(i,1) = ptop
390 #endif
391  enddo
392 
393  do k=2,km+1
394  do i=is1, ie1
395  pem(i,k) = pem(i,k-1) + dm(i,k-1)
396 #ifdef USE_COND
397 ! Excluding contribution from condensates:
398  peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1))
399 #endif
400  enddo
401  enddo
402 
403  do k=1,km
404  do i=is1, ie1
405  dz2(i,k) = gz(i,j,k+1) - gz(i,j,k)
406 #ifdef USE_COND
407  pm2(i,k) = (peg(i,k+1)-peg(i,k))/log(peg(i,k+1)/peg(i,k))
408 
409 #ifdef MOIST_CAPPA
410  cp2(i,k) = cappa(i,j,k)
411  gm2(i,k) = 1. / (1.-cp2(i,k))
412 #endif
413 
414 #else
415  pm2(i,k) = dm(i,k)/log(pem(i,k+1)/pem(i,k))
416 #endif
417 #ifdef MULTI_GASES
418  kapad2(i,k) = kapad(i,j,k)
419 #endif
420  dm(i,k) = dm(i,k) * rgrav
421  w2(i,k) = w3(i,j,k)
422  enddo
423  enddo
424 
425 
426  if ( a_imp < -0.01 ) then
427  call sim3p0_solver(dt, is1, ie1, km, rdgas, gama, akap, &
428 #ifdef MULTI_GASES
429  kapad2, &
430 #endif
431  pe2, dm, &
432  pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, scale_m)
433  elseif ( a_imp <= 0.5 ) then
434  call rim_2d(ms, dt, is1, ie1, km, rdgas, gama, gm2, &
435 #ifdef MULTI_GASES
436  kapad2, &
437 #endif
438  pe2, &
439  dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.)
440  else
441  call sim1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, &
442 #ifdef MULTI_GASES
443  kapad2, &
444 #endif
445  pe2, &
446  dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac)
447  endif
448 
449 
450  do k=2,km+1
451  do i=is1, ie1
452  pef(i,j,k) = pe2(i,k) + pem(i,k) ! add hydrostatic full-component
453  enddo
454  enddo
455 
456 ! Compute Height * grav (for p-gradient computation)
457  do i=is1, ie1
458  gz(i,j,km+1) = hs(i,j)
459  enddo
460 
461  do k=km,1,-1
462  do i=is1, ie1
463  gz(i,j,k) = gz(i,j,k+1) - dz2(i,k)*grav
464  enddo
465  enddo
466 
467 2000 continue
468 
469  end subroutine riem_solver_c
470 
471 
474  subroutine riem_solver3test(ms, dt, is, ie, js, je, km, ng, &
475  isd, ied, jsd, jed, akap, cappa, cp, &
476 #ifdef MULTI_GASES
477  kapad, &
478 #endif
479  ptop, zs, q_con, w, delz, pt, &
480  delp, zh, pe, ppe, pk3, pk, peln, &
481  ws, scale_m, p_fac, a_imp, &
482  use_logp, last_call, fp_out)
483 !--------------------------------------------
484 ! !OUTPUT PARAMETERS
485 ! Ouput: gz: grav*height at edges
486 ! pe: full hydrostatic pressure
487 ! ppe: non-hydrostatic pressure perturbation
488 !--------------------------------------------
489  integer, intent(in):: ms, is, ie, js, je, km, ng
490  integer, intent(in):: isd, ied, jsd, jed
491  real, intent(in):: dt ! the BIG horizontal Lagrangian time step
492  real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m
493  real, intent(in):: zs(isd:ied,jsd:jed)
494  logical, intent(in):: last_call, use_logp, fp_out
495  real, intent(in):: ws(is:ie,js:je)
496  real, intent(in), dimension(isd:,jsd:,1:):: q_con, cappa
497  real, intent(in), dimension(isd:ied,jsd:jed,km):: delp, pt
498 #ifdef MULTI_GASES
499  real, intent(in), dimension(isd:ied,jsd:jed,km):: kapad
500 #endif
501  real, intent(inout), dimension(isd:ied,jsd:jed,km+1):: zh
502  real, intent(inout), dimension(isd:ied,jsd:jed,km):: w
503  real, intent(inout):: pe(is-1:ie+1,km+1,js-1:je+1)
504  real, intent(out):: peln(is:ie,km+1,js:je) ! ln(pe)
505  real, intent(out), dimension(isd:ied,jsd:jed,km+1):: ppe
506  real, intent(out):: delz(is-ng:ie+ng,js-ng:je+ng,km)
507  real, intent(out):: pk(is:ie,js:je,km+1)
508  real, intent(out):: pk3(isd:ied,jsd:jed,km+1)
509 ! Local:
510  real, dimension(is:ie,km):: dm, dz2, pm2, w2, gm2, cp2
511  real, dimension(is:ie,km+1)::pem, pe2, peln2, peg, pelng
512 #ifdef MULTI_GASES
513  real, dimension(is:ie,km):: kapad2
514 #endif
515  real gama, rgrav, ptk, peln1
516  integer i, j, k
517 
518  gama = 1./(1.-akap)
519  rgrav = 1./grav
520  peln1 = log(ptop)
521  ptk = exp(akap*peln1)
522 
523 !$OMP parallel do default(none) shared(is,ie,js,je,km,delp,ptop,peln1,pk3,ptk,akap,rgrav,zh,pt, &
524 !$OMP w,a_imp,dt,gama,ws,p_fac,scale_m,ms,delz,last_call, &
525 #ifdef MULTI_GASES
526 !$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad ) &
527 !$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2,kapad2)
528 #else
529 !$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con ) &
530 !$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2)
531 #endif
532  do 2000 j=js, je
533 
534  do k=1,km
535  do i=is, ie
536  dm(i,k) = delp(i,j,k)
537 #ifdef MOIST_CAPPA
538  cp2(i,k) = cappa(i,j,k)
539 #endif
540 #ifdef MULTI_GASES
541  kapad2(i,k) = kapad(i,j,k)
542 #endif
543  enddo
544  enddo
545 
546  do i=is,ie
547  pem(i,1) = ptop
548  peln2(i,1) = peln1
549  pk3(i,j,1) = ptk
550 #ifdef USE_COND
551  peg(i,1) = ptop
552  pelng(i,1) = peln1
553 #endif
554  enddo
555  do k=2,km+1
556  do i=is,ie
557  pem(i,k) = pem(i,k-1) + dm(i,k-1)
558  peln2(i,k) = log(pem(i,k))
559 #ifdef USE_COND
560 ! Excluding contribution from condensates:
561 ! peln used during remap; pk3 used only for p_grad
562  peg(i,k) = peg(i,k-1) + dm(i,k-1)*(1.-q_con(i,j,k-1))
563  pelng(i,k) = log(peg(i,k))
564 #endif
565  pk3(i,j,k) = exp(akap*peln2(i,k))
566  enddo
567  enddo
568 
569  do k=1,km
570  do i=is, ie
571 #ifdef USE_COND
572  pm2(i,k) = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
573 
574 #ifdef MOIST_CAPPA
575  gm2(i,k) = 1. / (1.-cp2(i,k))
576 #endif
577 
578 #else
579  pm2(i,k) = dm(i,k)/(peln2(i,k+1)-peln2(i,k))
580 #endif
581  dm(i,k) = dm(i,k) * rgrav
582  dz2(i,k) = zh(i,j,k+1) - zh(i,j,k)
583  w2(i,k) = w(i,j,k)
584  enddo
585  enddo
586 
587 
588  if ( a_imp < -0.999 ) then
589  call sim3p0_solver(dt, is, ie, km, rdgas, gama, akap, &
590 #ifdef MULTI_GASES
591  kapad2, &
592 #endif
593  pe2, dm, &
594  pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, scale_m )
595  elseif ( a_imp < -0.5 ) then
596  call sim3_solver(dt, is, ie, km, rdgas, gama, akap, &
597 #ifdef MULTI_GASES
598  kapad2, &
599 #endif
600  pe2, dm, &
601  pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), abs(a_imp), p_fac, scale_m)
602  elseif ( a_imp <= 0.5 ) then
603  call rim_2d(ms, dt, is, ie, km, rdgas, gama, gm2, &
604 #ifdef MULTI_GASES
605  kapad2, &
606 #endif
607  pe2, &
608  dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.)
609  elseif ( a_imp > 0.999 ) then
610  call sim1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
611 #ifdef MULTI_GASES
612  kapad2, &
613 #endif
614  pe2, dm, &
615  pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac)
616  else
617  call sim_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
618 #ifdef MULTI_GASES
619  kapad2, &
620 #endif
621  pe2, dm, &
622  pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), &
623  a_imp, p_fac, scale_m)
624  endif
625 
626 
627  do k=1, km
628  do i=is, ie
629  w(i,j,k) = w2(i,k)
630  delz(i,j,k) = dz2(i,k)
631  enddo
632  enddo
633 
634  if ( last_call ) then
635  do k=1,km+1
636  do i=is,ie
637  peln(i,k,j) = peln2(i,k)
638  pk(i,j,k) = pk3(i,j,k)
639  pe(i,k,j) = pem(i,k)
640  enddo
641  enddo
642  endif
643 
644  if( fp_out ) then
645  do k=1,km+1
646  do i=is, ie
647  ppe(i,j,k) = pe2(i,k) + pem(i,k)
648  enddo
649  enddo
650  else
651  do k=1,km+1
652  do i=is, ie
653  ppe(i,j,k) = pe2(i,k)
654  enddo
655  enddo
656  endif
657 
658  if ( use_logp ) then
659  do k=2,km+1
660  do i=is, ie
661  pk3(i,j,k) = peln2(i,k)
662  enddo
663  enddo
664  endif
665 
666  do i=is, ie
667  zh(i,j,km+1) = zs(i,j)
668  enddo
669  do k=km,1,-1
670  do i=is, ie
671  zh(i,j,k) = zh(i,j,k+1) - dz2(i,k)
672  enddo
673  enddo
674 
675 2000 continue
676 
677  end subroutine riem_solver3test
678 
679  subroutine imp_diff_w(j, is, ie, js, je, ng, km, cd, delz, ws, w, w3)
680  integer, intent(in) :: j, is, ie, js, je, km, ng
681  real, intent(in) :: cd
682  real, intent(in) :: delz(is-ng:ie+ng, km)
683  real, intent(in) :: w(is:ie, km)
684  real, intent(in) :: ws(is:ie)
685  real, intent(out) :: w3(is-ng:ie+ng,js-ng:je+ng,km)
686 ! Local:
687  real, dimension(is:ie,km):: c, gam, dz, wt
688  real:: bet(is:ie)
689  real:: a
690  integer:: i, k
691 
692  do k=2,km
693  do i=is,ie
694  dz(i,k) = 0.5*(delz(i,k-1)+delz(i,k))
695  enddo
696  enddo
697  do k=1,km-1
698  do i=is,ie
699  c(i,k) = -cd/(dz(i,k+1)*delz(i,k))
700  enddo
701  enddo
702 
703 ! model top:
704  do i=is,ie
705  bet(i) = 1. - c(i,1) ! bet(i) = b
706  wt(i,1) = w(i,1) / bet(i)
707  enddo
708 
709 ! Interior:
710  do k=2,km-1
711  do i=is,ie
712  gam(i,k) = c(i,k-1)/bet(i)
713  a = cd/(dz(i,k)*delz(i,k))
714  bet(i) = (1.+a-c(i,k)) + a*gam(i,k)
715  wt(i,k) = (w(i,k) + a*wt(i,k-1)) / bet(i)
716  enddo
717  enddo
718 
719 ! Bottom:
720  do i=is,ie
721  gam(i,km) = c(i,km-1) / bet(i)
722  a = cd/(dz(i,km)*delz(i,km))
723  wt(i,km) = (w(i,km) + 2.*ws(i)*cd/delz(i,km)**2 &
724  + a*wt(i,km-1))/(1. + a + (cd+cd)/delz(i,km)**2 + a*gam(i,km))
725  enddo
726 
727  do k=km-1,1,-1
728  do i=is,ie
729  wt(i,k) = wt(i,k) - gam(i,k+1)*wt(i,k+1)
730  enddo
731  enddo
732 
733  do k=1,km
734  do i=is,ie
735  w3(i,j,k) = wt(i,k)
736  enddo
737  enddo
738 
739  end subroutine imp_diff_w
740 
741 
742  subroutine rim_2d(ms, bdt, is, ie, km, rgas, gama, gm2, &
743 #ifdef MULTI_GASES
744  kapad2, &
745 #endif
746  pe2, dm2, pm2, w2, dz2, pt2, ws, c_core )
747 
748  integer, intent(in):: ms, is, ie, km
749  real, intent(in):: bdt, gama, rgas
750  real, intent(in), dimension(is:ie,km):: dm2, pm2, gm2
751  logical, intent(in):: c_core
752  real, intent(in ) :: pt2(is:ie,km)
753  real, intent(in ) :: ws(is:ie)
754 ! IN/OUT:
755  real, intent(inout):: dz2(is:ie,km)
756  real, intent(inout):: w2(is:ie,km)
757  real, intent(out ):: pe2(is:ie,km+1)
758 #ifdef MULTI_GASES
759  real, intent(inout), dimension(is:ie,km):: kapad2
760 #endif
761 ! Local:
762  real:: ws2(is:ie)
763  real, dimension(km+1):: m_bot, m_top, r_bot, r_top, pe1, pbar, wbar
764  real, dimension(km):: r_hi, r_lo, dz, wm, dm, dts
765  real, dimension(km):: pf1, wc, cm , pp, pt1
766  real:: dt, rdt, grg, z_frac, ptmp1, rden, pf, time_left
767  real:: m_surf
768 #ifdef MULTI_GASES
769  real gamax
770 #endif
771  integer:: i, k, n, ke, kt1, ktop
772  integer:: ks0, ks1
773 
774  grg = gama * rgas
775  rdt = 1. / bdt
776  dt = bdt / real(ms)
777 
778  pbar(:) = 0.
779  wbar(:) = 0.
780 
781  do i=is,ie
782  ws2(i) = 2.*ws(i)
783  enddo
784 
785  do 6000 i=is,ie
786 
787  do k=1,km
788  dz(k) = dz2(i,k)
789  dm(k) = dm2(i,k)
790  wm(k) = w2(i,k)*dm(k)
791  pt1(k) = pt2(i,k)
792  enddo
793 
794  pe1(:) = 0.
795  wbar(km+1) = ws(i)
796 
797  ks0 = 1
798  if ( ms > 1 .and. ms < 8 ) then
799 ! Continuity of (pbar, wbar) is maintained
800  do k=1, km
801  rden = -rgas*dm(k)/dz(k)
802 #ifdef MOIST_CAPPA
803  pf1(k) = exp( gm2(i,k)*log(rden*pt1(k)) )
804 ! dts(k) = -dz(k)/sqrt(gm2(i,k)*rgas*pf1(k)/rden)
805  dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden)
806 #else
807 #ifdef MULTI_GASES
808  gamax = 1./(1.-kapad2(i,k))
809  pf1(k) = exp( gamax*log(rden*pt1(k)) )
810 #else
811  pf1(k) = exp( gama*log(rden*pt1(k)) )
812 #endif
813  dts(k) = -dz(k)/sqrt(grg*pf1(k)/rden)
814 #endif
815  if ( bdt > dts(k) ) then
816  ks0 = k-1
817  goto 222
818  endif
819  enddo
820  ks0 = km
821 222 if ( ks0==1 ) goto 244
822 
823  do k=1, ks0
824  cm(k) = dm(k) / dts(k)
825  wc(k) = wm(k) / dts(k)
826  pp(k) = pf1(k) - pm2(i,k)
827  enddo
828 
829  wbar(1) = (wc(1)+pp(1)) / cm(1)
830  do k=2, ks0
831  wbar(k) = (wc(k-1)+wc(k) + pp(k)-pp(k-1)) / (cm(k-1)+cm(k))
832  pbar(k) = bdt*(cm(k-1)*wbar(k) - wc(k-1) + pp(k-1))
833  pe1(k) = pbar(k)
834  enddo
835 
836  if ( ks0 == km ) then
837  pbar(km+1) = bdt*(cm(km)*wbar(km+1) - wc(km) + pp(km))
838  if ( c_core ) then
839  do k=1,km
840  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
841  enddo
842  else
843  do k=1,km
844  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
845  w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k)
846  enddo
847  endif
848  pe2(i,1) = 0.
849  do k=2,km+1
850  pe2(i,k) = pbar(k)*rdt
851  enddo
852  goto 6000 ! next i
853  else
854  if ( c_core ) then
855  do k=1, ks0-1
856  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
857  enddo
858  else
859  do k=1, ks0-1
860  dz2(i,k) = dz(k) + bdt*(wbar(k+1) - wbar(k))
861  w2(i,k) = (wm(k)+pbar(k+1)-pbar(k))/dm(k)
862  enddo
863  endif
864  pbar(ks0) = pbar(ks0) / real(ms)
865  endif
866  endif
867 244 ks1 = ks0
868 
869  do 5000 n=1, ms
870 
871  do k=ks1, km
872  rden = -rgas*dm(k)/dz(k)
873 #ifdef MOIST_CAPPA
874  pf = exp( gm2(i,k)*log(rden*pt1(k)) )
875 ! dts(k) = -dz(k) / sqrt( gm2(i,k)*rgas*pf/rden )
876  dts(k) = -dz(k) / sqrt( grg*pf/rden )
877 #else
878 #ifdef MULTI_GASES
879  gamax = 1./(1.-kapad2(i,k))
880  pf = exp( gamax*log(rden*pt1(k)) )
881 #else
882  pf = exp( gama*log(rden*pt1(k)) )
883 #endif
884  dts(k) = -dz(k) / sqrt( grg*pf/rden )
885 #endif
886  ptmp1 = dts(k)*(pf - pm2(i,k))
887  r_lo(k) = wm(k) + ptmp1
888  r_hi(k) = wm(k) - ptmp1
889  enddo
890 
891  ktop = ks1
892  do k=ks1, km
893  if( dt > dts(k) ) then
894  ktop = k-1
895  goto 333
896  endif
897  enddo
898  ktop = km
899 333 continue
900 
901  if ( ktop >= ks1 ) then
902  do k=ks1, ktop
903  z_frac = dt/dts(k)
904  r_bot(k ) = z_frac*r_lo(k)
905  r_top(k+1) = z_frac*r_hi(k)
906  m_bot(k ) = z_frac* dm(k)
907  m_top(k+1) = m_bot(k)
908  enddo
909  if ( ktop == km ) goto 666
910  endif
911 
912  do k=ktop+2, km+1
913  m_top(k) = 0.
914  r_top(k) = 0.
915  enddo
916 
917  kt1 = max(1, ktop)
918  do 444 ke=km+1, ktop+2, -1
919  time_left = dt
920  do k=ke-1, kt1, -1
921  if ( time_left > dts(k) ) then
922  time_left = time_left - dts(k)
923  m_top(ke) = m_top(ke) + dm(k)
924  r_top(ke) = r_top(ke) + r_hi(k)
925  else
926  z_frac = time_left/dts(k)
927  m_top(ke) = m_top(ke) + z_frac*dm(k)
928  r_top(ke) = r_top(ke) + z_frac*r_hi(k)
929  go to 444 ! next level
930  endif
931  enddo
932 444 continue
933 
934  do k=ktop+1, km
935  m_bot(k) = 0.
936  r_bot(k) = 0.
937  enddo
938 
939  do 4000 ke=ktop+1, km
940  time_left = dt
941  do k=ke, km
942  if ( time_left > dts(k) ) then
943  time_left = time_left - dts(k)
944  m_bot(ke) = m_bot(ke) + dm(k)
945  r_bot(ke) = r_bot(ke) + r_lo(k)
946  else
947  z_frac = time_left/dts(k)
948  m_bot(ke) = m_bot(ke) + z_frac* dm(k)
949  r_bot(ke) = r_bot(ke) + z_frac*r_lo(k)
950  go to 4000 ! next interface
951  endif
952  enddo
953  m_surf = m_bot(ke)
954  do k=km, kt1, -1
955  if ( time_left > dts(k) ) then
956  time_left = time_left - dts(k)
957  m_bot(ke) = m_bot(ke) + dm(k)
958  r_bot(ke) = r_bot(ke) - r_hi(k)
959  else
960  z_frac = time_left/dts(k)
961  m_bot(ke) = m_bot(ke) + z_frac* dm(k)
962  r_bot(ke) = r_bot(ke) - z_frac*r_hi(k) + (m_bot(ke)-m_surf)*ws2(i)
963  go to 4000 ! next interface
964  endif
965  enddo
966 4000 continue
967 
968 666 if ( ks1==1 ) wbar(1) = r_bot(1) / m_bot(1)
969  do k=ks1+1, km
970  wbar(k) = (r_bot(k)+r_top(k)) / (m_top(k)+m_bot(k))
971  enddo
972 ! pbar here is actually dt*pbar
973  do k=ks1+1, km+1
974  pbar(k) = m_top(k)*wbar(k) - r_top(k)
975  pe1(k) = pe1(k) + pbar(k)
976  enddo
977 
978  if ( n==ms ) then
979  if ( c_core ) then
980  do k=ks1, km
981  dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k))
982  enddo
983  else
984  do k=ks1, km
985  dz2(i,k) = dz(k) + dt*(wbar(k+1)-wbar(k))
986  w2(i,k) = ( wm(k) + pbar(k+1) - pbar(k) ) / dm(k)
987  enddo
988  endif
989  else
990  do k=ks1, km
991  dz(k) = dz(k) + dt*(wbar(k+1)-wbar(k))
992  wm(k) = wm(k) + pbar(k+1) - pbar(k)
993  enddo
994  endif
995 
996 5000 continue
997  pe2(i,1) = 0.
998  do k=2,km+1
999  pe2(i,k) = pe1(k)*rdt
1000  enddo
1001 
1002 6000 continue ! end i-loop
1003 
1004  end subroutine rim_2d
1005 
1006  subroutine sim3_solver(dt, is, ie, km, rgas, gama, kappa, &
1007 #ifdef MULTI_GASES
1008  kapad2, &
1009 #endif
1010  pe2, dm, &
1011  pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
1012  integer, intent(in):: is, ie, km
1013  real, intent(in):: dt, rgas, gama, kappa, alpha, p_fac, scale_m
1014  real, intent(in ), dimension(is:ie,km):: dm, pt2
1015  real, intent(in ):: ws(is:ie)
1016  real, intent(in ), dimension(is:ie,km+1):: pem
1017  real, intent(out):: pe2(is:ie,km+1)
1018  real, intent(inout), dimension(is:ie,km):: dz2, w2
1019 #ifdef MULTI_GASES
1020  real, intent(inout), dimension(is:ie,km):: kapad2
1021 #endif
1022 ! Local
1023  real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam
1024  real, dimension(is:ie,km+1):: pp
1025  real, dimension(is:ie):: p1, wk1, bet
1026  real beta, t2, t1g, rdt, ra, capa1, r2g, r6g
1027 #ifdef MULTI_GASES
1028  real gamax, capa1x, t1gx
1029 #endif
1030  integer i, k
1031 
1032  beta = 1. - alpha
1033  ra = 1. / alpha
1034  t2 = beta / alpha
1035  t1g = gama * 2.*(alpha*dt)**2
1036  rdt = 1. / dt
1037  capa1 = kappa - 1.
1038  r2g = grav / 2.
1039  r6g = grav / 6.
1040 
1041 
1042  do k=1,km
1043  do i=is, ie
1044  w1(i,k) = w2(i,k)
1045 ! Full pressure at center
1046 #ifdef MULTI_GASES
1047  gamax = 1. / (1.-kapad2(i,k))
1048  aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1049 #else
1050  aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1051 #endif
1052  enddo
1053  enddo
1054 
1055  do k=1,km-1
1056  do i=is, ie
1057  g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction
1058  bb(i,k) = 2.*(1.+g_rat(i,k))
1059  dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1))
1060  enddo
1061  enddo
1062 
1063 ! pe2 is full p at edges
1064  do i=is, ie
1065 ! Top:
1066  bet(i) = bb(i,1)
1067  pe2(i,1) = pem(i,1)
1068  pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i)
1069 ! Bottom:
1070  bb(i,km) = 2.
1071  dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km)
1072  enddo
1073 
1074  do k=2,km
1075  do i=is, ie
1076  gam(i,k) = g_rat(i,k-1) / bet(i)
1077  bet(i) = bb(i,k) - gam(i,k)
1078  pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i)
1079  enddo
1080  enddo
1081 
1082  do k=km, 2, -1
1083  do i=is, ie
1084  pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1)
1085  enddo
1086  enddo
1087 ! done reconstruction of full:
1088 
1089 ! pp is pert. p at edges
1090  do k=1, km+1
1091  do i=is, ie
1092  pp(i,k) = pe2(i,k) - pem(i,k)
1093  enddo
1094  enddo
1095 
1096  do k=2, km
1097  do i=is, ie
1098 #ifdef MULTI_GASES
1099  gamax = 1./(1.-kapad2(i,k))
1100  t1gx = gamax*2.*(alpha*dt)**2
1101  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1102 #else
1103  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1104 #endif
1105  wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k))
1106  aa(i,k) = aa(i,k) - scale_m*dm(i,1)
1107  enddo
1108  enddo
1109  do i=is, ie
1110  bet(i) = dm(i,1) - aa(i,2)
1111  w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2) + wk(i,2)) / bet(i)
1112  enddo
1113  do k=2,km-1
1114  do i=is, ie
1115  gam(i,k) = aa(i,k) / bet(i)
1116  bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1117  w2(i,k) = (dm(i,k)*w1(i,k)+dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) &
1118  - aa(i,k)*w2(i,k-1)) / bet(i)
1119  enddo
1120  enddo
1121  do i=is, ie
1122 #ifdef MULTI_GASES
1123  gamax = 1./(1.-kapad2(i,km))
1124  t1gx = gamax*2.*(alpha*dt)**2
1125  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1126 #else
1127  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1128 #endif
1129  gam(i,km) = aa(i,km) / bet(i)
1130  bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1131  w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk(i,km) + &
1132  wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i)
1133  enddo
1134  do k=km-1, 1, -1
1135  do i=is, ie
1136  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1137  enddo
1138  enddo
1139 
1140 ! pe2 is updated perturbation p at edges
1141  do i=is, ie
1142  pe2(i,1) = 0.
1143  enddo
1144  do k=1,km
1145  do i=is, ie
1146  pe2(i,k+1) = pe2(i,k) + ( dm(i,k)*(w2(i,k)-w1(i,k))*rdt &
1147  - beta*(pp(i,k+1)-pp(i,k)) )*ra
1148  enddo
1149  enddo
1150 
1151 ! Full non-hydro pressure at edges:
1152  do i=is, ie
1153  pe2(i,1) = pem(i,1)
1154  enddo
1155  do k=2,km+1
1156  do i=is, ie
1157  pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k))
1158  enddo
1159  enddo
1160 
1161  do i=is, ie
1162 ! Recover cell-averaged pressure
1163  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km)
1164 #ifdef MULTI_GASES
1165  capa1x = kapad2(i,km) - 1.
1166  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) )
1167 #else
1168  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) )
1169 #endif
1170  enddo
1171 
1172  do k=km-1, 1, -1
1173  do i=is, ie
1174  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)
1175 #ifdef MULTI_GASES
1176  capa1x = kapad2(i,k) - 1.
1177  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) )
1178 #else
1179  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) )
1180 #endif
1181  enddo
1182  enddo
1183 
1184  do k=1,km+1
1185  do i=is, ie
1186  pe2(i,k) = pe2(i,k) - pem(i,k)
1187  pe2(i,k) = pe2(i,k) + beta*(pp(i,k) - pe2(i,k))
1188  enddo
1189  enddo
1190 
1191  end subroutine sim3_solver
1192 
1193  subroutine sim3p0_solver(dt, is, ie, km, rgas, gama, kappa, &
1194 #ifdef MULTI_GASES
1195  kapad2, &
1196 #endif
1197  pe2, dm, &
1198  pem, w2, dz2, pt2, ws, p_fac, scale_m)
1199 ! Sa SIM3, but for beta==0
1200  integer, intent(in):: is, ie, km
1201  real, intent(in):: dt, rgas, gama, kappa, p_fac, scale_m
1202  real, intent(in ), dimension(is:ie,km):: dm, pt2
1203  real, intent(in ):: ws(is:ie)
1204  real, intent(in ):: pem(is:ie,km+1)
1205  real, intent(out):: pe2(is:ie,km+1)
1206  real, intent(inout), dimension(is:ie,km):: dz2, w2
1207 #ifdef MULTI_GASES
1208  real, intent(inout), dimension(is:ie,km):: kapad2
1209 #endif
1210 ! Local
1211  real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam
1212  real, dimension(is:ie,km+1):: pp
1213  real, dimension(is:ie):: p1, wk1, bet
1214  real t1g, rdt, capa1, r2g, r6g
1215 #ifdef MULTI_GASES
1216  real gamax, capa1x, t1gx
1217 #endif
1218  integer i, k
1219 
1220  t1g = 2.*gama*dt**2
1221  rdt = 1. / dt
1222  capa1 = kappa - 1.
1223  r2g = grav / 2.
1224  r6g = grav / 6.
1225 
1226  do k=1,km
1227  do i=is, ie
1228  w1(i,k) = w2(i,k)
1229 ! Full pressure at center
1230 #ifdef MULTI_GASES
1231  gamax = 1. / ( 1. - kapad2(i,k) )
1232  aa(i,k) = exp(gamax*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1233 #else
1234  aa(i,k) = exp(gama*log(-dm(i,k)/dz2(i,k)*rgas*pt2(i,k)))
1235 #endif
1236  enddo
1237  enddo
1238 
1239  do k=1,km-1
1240  do i=is, ie
1241  g_rat(i,k) = dm(i,k)/dm(i,k+1) ! for profile reconstruction
1242  bb(i,k) = 2.*(1.+g_rat(i,k))
1243  dd(i,k) = 3.*(aa(i,k) + g_rat(i,k)*aa(i,k+1))
1244  enddo
1245  enddo
1246 
1247 ! pe2 is full p at edges
1248  do i=is, ie
1249 ! Top:
1250  bet(i) = bb(i,1)
1251  pe2(i,1) = pem(i,1)
1252  pe2(i,2) = (dd(i,1)-pem(i,1)) / bet(i)
1253 ! Bottom:
1254  bb(i,km) = 2.
1255  dd(i,km) = 3.*aa(i,km) + r2g*dm(i,km)
1256  enddo
1257 
1258  do k=2,km
1259  do i=is, ie
1260  gam(i,k) = g_rat(i,k-1) / bet(i)
1261  bet(i) = bb(i,k) - gam(i,k)
1262  pe2(i,k+1) = (dd(i,k) - pe2(i,k) ) / bet(i)
1263  enddo
1264  enddo
1265 
1266  do k=km, 2, -1
1267  do i=is, ie
1268  pe2(i,k) = pe2(i,k) - gam(i,k)*pe2(i,k+1)
1269  enddo
1270  enddo
1271 ! done reconstruction of full:
1272 
1273 ! pp is pert. p at edges
1274  do k=1, km+1
1275  do i=is, ie
1276  pp(i,k) = pe2(i,k) - pem(i,k)
1277  enddo
1278  enddo
1279 
1280  do k=2, km
1281  do i=is, ie
1282 #ifdef MULTI_GASES
1283  gamax = 1. / (1.-kapad2(i,k))
1284  t1gx = 2.*gamax*dt**2
1285  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1)
1286 #else
1287  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k) - scale_m*dm(i,1)
1288 #endif
1289  enddo
1290  enddo
1291  do i=is, ie
1292  bet(i) = dm(i,1) - aa(i,2)
1293  w2(i,1) = (dm(i,1)*w1(i,1)+dt*pp(i,2)) / bet(i)
1294  enddo
1295  do k=2,km-1
1296  do i=is, ie
1297  gam(i,k) = aa(i,k) / bet(i)
1298  bet(i) = dm(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1299  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)
1300  enddo
1301  enddo
1302  do i=is, ie
1303 #ifdef MULTI_GASES
1304  gamax = 1. / (1.-kapad2(i,km))
1305  t1gx = 2.*gamax*dt**2
1306  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1307 #else
1308  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1309 #endif
1310  gam(i,km) = aa(i,km) / bet(i)
1311  bet(i) = dm(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1312  w2(i,km) = (dm(i,km)*w1(i,km)+dt*(pp(i,km+1)-pp(i,km))-wk1(i)*ws(i) - &
1313  aa(i,km)*w2(i,km-1)) / bet(i)
1314  enddo
1315  do k=km-1, 1, -1
1316  do i=is, ie
1317  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1318  enddo
1319  enddo
1320 
1321 ! pe2 is updated perturbation p at edges
1322  do i=is, ie
1323  pe2(i,1) = 0.
1324  enddo
1325  do k=1,km
1326  do i=is, ie
1327  pe2(i,k+1) = pe2(i,k) + dm(i,k)*(w2(i,k)-w1(i,k))*rdt
1328  enddo
1329  enddo
1330 
1331 ! Full non-hydro pressure at edges:
1332  do i=is, ie
1333  pe2(i,1) = pem(i,1)
1334  enddo
1335  do k=2,km+1
1336  do i=is, ie
1337  pe2(i,k) = max(p_fac*pem(i,k), pe2(i,k)+pem(i,k))
1338  enddo
1339  enddo
1340 
1341  do i=is, ie
1342 ! Recover cell-averaged pressure
1343  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3 - r6g*dm(i,km)
1344 #ifdef MULTI_GASES
1345  capa1x = kapad2(i,km) - 1.
1346  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1x*log(p1(i)) )
1347 #else
1348  dz2(i,km) = -dm(i,km)*rgas*pt2(i,km)*exp( capa1*log(p1(i)) )
1349 #endif
1350  enddo
1351 
1352  do k=km-1, 1, -1
1353  do i=is, ie
1354  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)
1355 #ifdef MULTI_GASES
1356  capa1x = kapad2(i,k) - 1.
1357  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1x*log(p1(i)) )
1358 #else
1359  dz2(i,k) = -dm(i,k)*rgas*pt2(i,k)*exp( capa1*log(p1(i)) )
1360 #endif
1361  enddo
1362  enddo
1363 
1364  do k=1,km+1
1365  do i=is, ie
1366  pe2(i,k) = pe2(i,k) - pem(i,k)
1367  enddo
1368  enddo
1369 
1370  end subroutine sim3p0_solver
1371 
1372 
1373  subroutine sim1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
1374 #ifdef MULTI_GASES
1375  kapad2, &
1376 #endif
1377  pe, dm2, &
1378  pm2, pem, w2, dz2, pt2, ws, p_fac)
1379  integer, intent(in):: is, ie, km
1380  real, intent(in):: dt, rgas, gama, kappa, p_fac
1381  real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
1382  real, intent(in ):: ws(is:ie)
1383  real, intent(in ), dimension(is:ie,km+1):: pem
1384  real, intent(out):: pe(is:ie,km+1)
1385  real, intent(inout), dimension(is:ie,km):: dz2, w2
1386 #ifdef MULTI_GASES
1387  real, intent(inout), dimension(is:ie,km):: kapad2
1388 #endif
1389 ! Local
1390  real, dimension(is:ie,km ):: aa, bb, dd, w1, g_rat, gam
1391  real, dimension(is:ie,km+1):: pp
1392  real, dimension(is:ie):: p1, bet
1393  real t1g, rdt, capa1
1394 #ifdef MULTI_GASES
1395  real gamax, capa1x, t1gx
1396 #endif
1397  integer i, k
1398 
1399 #ifdef MOIST_CAPPA
1400  t1g = 2.*dt*dt
1401 #else
1402  t1g = gama * 2.*dt*dt
1403 #endif
1404  rdt = 1. / dt
1405  capa1 = kappa - 1.
1406 
1407  do k=1,km
1408  do i=is, ie
1409 #ifdef MOIST_CAPPA
1410  pe(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1411 #else
1412 #ifdef MULTI_GASES
1413  gamax = 1. / ( 1. - kapad2(i,k) )
1414  pe(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1415 #else
1416  pe(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1417 #endif
1418 #endif
1419  w1(i,k) = w2(i,k)
1420  enddo
1421  enddo
1422 
1423  do k=1,km-1
1424  do i=is, ie
1425  g_rat(i,k) = dm2(i,k)/dm2(i,k+1)
1426  bb(i,k) = 2.*(1.+g_rat(i,k))
1427  dd(i,k) = 3.*(pe(i,k) + g_rat(i,k)*pe(i,k+1))
1428  enddo
1429  enddo
1430 
1431  do i=is, ie
1432  bet(i) = bb(i,1)
1433  pp(i,1) = 0.
1434  pp(i,2) = dd(i,1) / bet(i)
1435  bb(i,km) = 2.
1436  dd(i,km) = 3.*pe(i,km)
1437  enddo
1438 
1439  do k=2,km
1440  do i=is, ie
1441  gam(i,k) = g_rat(i,k-1) / bet(i)
1442  bet(i) = bb(i,k) - gam(i,k)
1443  pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i)
1444  enddo
1445  enddo
1446 
1447  do k=km, 2, -1
1448  do i=is, ie
1449  pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1)
1450  enddo
1451  enddo
1452 
1453 ! Start the w-solver
1454  do k=2, km
1455  do i=is, ie
1456 #ifdef MOIST_CAPPA
1457  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))
1458 #else
1459 #ifdef MULTI_GASES
1460  gamax = 1./(1.-kapad2(i,k))
1461  t1gx = gamax * 2.*dt*dt
1462  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
1463 #else
1464  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
1465 #endif
1466 #endif
1467  enddo
1468  enddo
1469  do i=is, ie
1470  bet(i) = dm2(i,1) - aa(i,2)
1471  w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2)) / bet(i)
1472  enddo
1473  do k=2,km-1
1474  do i=is, ie
1475  gam(i,k) = aa(i,k) / bet(i)
1476  bet(i) = dm2(i,k) - (aa(i,k) + aa(i,k+1) + aa(i,k)*gam(i,k))
1477  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)
1478  enddo
1479  enddo
1480  do i=is, ie
1481 #ifdef MOIST_CAPPA
1482  p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1483 #else
1484 #ifdef MULTI_GASES
1485  gamax = 1./(1.-kapad2(i,km))
1486  t1gx = gamax * 2.*dt*dt
1487  p1(i) = t1gx/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1488 #else
1489  p1(i) = t1g/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
1490 #endif
1491 #endif
1492  gam(i,km) = aa(i,km) / bet(i)
1493  bet(i) = dm2(i,km) - (aa(i,km)+p1(i) + aa(i,km)*gam(i,km))
1494  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)
1495  enddo
1496  do k=km-1, 1, -1
1497  do i=is, ie
1498  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1499  enddo
1500  enddo
1501 
1502  do i=is, ie
1503  pe(i,1) = 0.
1504  enddo
1505  do k=1,km
1506  do i=is, ie
1507  pe(i,k+1) = pe(i,k) + dm2(i,k)*(w2(i,k)-w1(i,k))*rdt
1508  enddo
1509  enddo
1510 
1511  do i=is, ie
1512  p1(i) = ( pe(i,km) + 2.*pe(i,km+1) )*r3
1513 #ifdef MOIST_CAPPA
1514  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))))
1515 #else
1516 #ifdef MULTI_GASES
1517  capa1x = kapad2(i,km)-1.
1518  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1519 #else
1520  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1521 #endif
1522 #endif
1523  enddo
1524 
1525  do k=km-1, 1, -1
1526  do i=is, ie
1527  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)
1528 #ifdef MOIST_CAPPA
1529  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))))
1530 
1531 #else
1532 #ifdef MULTI_GASES
1533  capa1x = kapad2(i,k)-1.
1534  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1535 #else
1536  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1537 #endif
1538 #endif
1539  enddo
1540  enddo
1541 
1542  end subroutine sim1_solver
1543 
1544  subroutine sim_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
1545 #ifdef MULTI_GASES
1546  kapad2, &
1547 #endif
1548  pe2, dm2, &
1549  pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
1550  integer, intent(in):: is, ie, km
1551  real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m
1552  real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
1553  real, intent(in ):: ws(is:ie)
1554  real, intent(in ), dimension(is:ie,km+1):: pem
1555  real, intent(out):: pe2(is:ie,km+1)
1556  real, intent(inout), dimension(is:ie,km):: dz2, w2
1557 #ifdef MULTI_GASES
1558  real, intent(inout), dimension(is:ie,km):: kapad2
1559 #endif
1560 ! Local
1561  real, dimension(is:ie,km ):: aa, bb, dd, w1, wk, g_rat, gam
1562  real, dimension(is:ie,km+1):: pp
1563  real, dimension(is:ie):: p1, wk1, bet
1564  real beta, t2, t1g, rdt, ra, capa1
1565 #ifdef MULTI_GASES
1566  real gamax, capa1x, t1gx
1567 #endif
1568  integer i, k
1569 
1570  beta = 1. - alpha
1571  ra = 1. / alpha
1572  t2 = beta / alpha
1573 #ifdef MOIST_CAPPA
1574  t1g = 2.*(alpha*dt)**2
1575 #else
1576  t1g = 2.*gama*(alpha*dt)**2
1577 #endif
1578  rdt = 1. / dt
1579  capa1 = kappa - 1.
1580 
1581  do k=1,km
1582  do i=is, ie
1583  w1(i,k) = w2(i,k)
1584 ! P_g perturbation
1585 #ifdef MOIST_CAPPA
1586  pe2(i,k) = exp(gm2(i,k)*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1587 #else
1588 #ifdef MULTI_GASES
1589  gamax = 1./(1.-kapad2(i,k))
1590  pe2(i,k) = exp(gamax*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1591 #else
1592  pe2(i,k) = exp(gama*log(-dm2(i,k)/dz2(i,k)*rgas*pt2(i,k))) - pm2(i,k)
1593 #endif
1594 #endif
1595  enddo
1596  enddo
1597 
1598  do k=1,km-1
1599  do i=is, ie
1600  g_rat(i,k) = dm2(i,k)/dm2(i,k+1)
1601  bb(i,k) = 2.*(1.+g_rat(i,k))
1602  dd(i,k) = 3.*(pe2(i,k) + g_rat(i,k)*pe2(i,k+1))
1603  enddo
1604  enddo
1605 
1606  do i=is, ie
1607  bet(i) = bb(i,1)
1608  pp(i,1) = 0.
1609  pp(i,2) = dd(i,1) / bet(i)
1610  bb(i,km) = 2.
1611  dd(i,km) = 3.*pe2(i,km)
1612  enddo
1613 
1614  do k=2,km
1615  do i=is, ie
1616  gam(i,k) = g_rat(i,k-1) / bet(i)
1617  bet(i) = bb(i,k) - gam(i,k)
1618  pp(i,k+1) = (dd(i,k) - pp(i,k) ) / bet(i)
1619  enddo
1620  enddo
1621 
1622  do k=km, 2, -1
1623  do i=is, ie
1624  pp(i,k) = pp(i,k) - gam(i,k)*pp(i,k+1)
1625  enddo
1626  enddo
1627 
1628  do k=1, km+1
1629  do i=is, ie
1630 ! pe2 is Full p
1631  pe2(i,k) = pem(i,k) + pp(i,k)
1632  enddo
1633  enddo
1634 
1635  do k=2, km
1636  do i=is, ie
1637 #ifdef MOIST_CAPPA
1638  aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1639 #else
1640 #ifdef MULTI_GASES
1641  gamax = 1./(1.-kapad2(i,k))
1642  t1gx = 2.*gamax*(alpha*dt)**2
1643  aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1644 #else
1645  aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k))*pe2(i,k)
1646 #endif
1647 #endif
1648  wk(i,k) = t2*aa(i,k)*(w1(i,k-1)-w1(i,k))
1649  aa(i,k) = aa(i,k) - scale_m*dm2(i,1)
1650  enddo
1651  enddo
1652 ! Top:
1653  do i=is, ie
1654  bet(i) = dm2(i,1) - aa(i,2)
1655  w2(i,1) = (dm2(i,1)*w1(i,1) + dt*pp(i,2) + wk(i,2)) / bet(i)
1656  enddo
1657 ! Interior:
1658  do k=2,km-1
1659  do i=is, ie
1660  gam(i,k) = aa(i,k) / bet(i)
1661  bet(i) = dm2(i,k) - (aa(i,k)+aa(i,k+1) + aa(i,k)*gam(i,k))
1662  w2(i,k) = (dm2(i,k)*w1(i,k) + dt*(pp(i,k+1)-pp(i,k)) + wk(i,k+1)-wk(i,k) &
1663  - aa(i,k)*w2(i,k-1)) / bet(i)
1664  enddo
1665  enddo
1666 ! Bottom: k=km
1667  do i=is, ie
1668 #ifdef MOIST_CAPPA
1669  wk1(i) = t1g*gm2(i,km)/dz2(i,km)*pe2(i,km+1)
1670 #else
1671 #ifdef MULTI_GASES
1672  gamax = 1./(1.-kapad2(i,km))
1673  t1gx = 2.*gamax*(alpha*dt)**2
1674  wk1(i) = t1gx/dz2(i,km)*pe2(i,km+1)
1675 #else
1676  wk1(i) = t1g/dz2(i,km)*pe2(i,km+1)
1677 #endif
1678 #endif
1679  gam(i,km) = aa(i,km) / bet(i)
1680  bet(i) = dm2(i,km) - (aa(i,km)+wk1(i) + aa(i,km)*gam(i,km))
1681  w2(i,km) = (dm2(i,km)*w1(i,km) + dt*(pp(i,km+1)-pp(i,km)) - wk(i,km) + &
1682  wk1(i)*(t2*w1(i,km)-ra*ws(i)) - aa(i,km)*w2(i,km-1)) / bet(i)
1683  enddo
1684  do k=km-1, 1, -1
1685  do i=is, ie
1686  w2(i,k) = w2(i,k) - gam(i,k+1)*w2(i,k+1)
1687  enddo
1688  enddo
1689 
1690  do i=is, ie
1691  pe2(i,1) = 0.
1692  enddo
1693  do k=1,km
1694  do i=is, ie
1695  pe2(i,k+1) = pe2(i,k) + ( dm2(i,k)*(w2(i,k)-w1(i,k))*rdt &
1696  - beta*(pp(i,k+1)-pp(i,k)) ) * ra
1697  enddo
1698  enddo
1699 
1700  do i=is, ie
1701  p1(i) = (pe2(i,km)+ 2.*pe2(i,km+1))*r3
1702 #ifdef MOIST_CAPPA
1703  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))))
1704 #else
1705 #ifdef MULTI_GASES
1706  capa1x = kapad2(i,km)-1.
1707  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1x*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1708 #else
1709  dz2(i,km) = -dm2(i,km)*rgas*pt2(i,km)*exp(capa1*log(max(p_fac*pm2(i,km),p1(i)+pm2(i,km))))
1710 #endif
1711 #endif
1712  enddo
1713 
1714  do k=km-1, 1, -1
1715  do i=is, ie
1716  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)
1717 ! delz = -dm*R*T_m / p_gas
1718 #ifdef MOIST_CAPPA
1719  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))))
1720 #else
1721 #ifdef MULTI_GASES
1722  capa1x = kapad2(i,k)-1.
1723  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1x*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1724 #else
1725  dz2(i,k) = -dm2(i,k)*rgas*pt2(i,k)*exp(capa1*log(max(p_fac*pm2(i,k),p1(i)+pm2(i,k))))
1726 #endif
1727 #endif
1728  enddo
1729  enddo
1730 
1731  do k=1, km+1
1732  do i=is, ie
1733  pe2(i,k) = pe2(i,k) + beta*(pp(i,k)-pe2(i,k))
1734  enddo
1735  enddo
1736 
1737  end subroutine sim_solver
1738 
1739 
1740  subroutine edge_scalar(q1, qe, i1, i2, km, id)
1741 ! Optimized for wind profile reconstruction:
1742  integer, intent(in):: i1, i2, km
1743  integer, intent(in):: id ! 0: pp 1: wind
1744  real, intent(in ), dimension(i1:i2,km):: q1
1745  real, intent(out), dimension(i1:i2,km+1):: qe
1746 !-----------------------------------------------------------------------
1747  real, parameter:: r2o3 = 2./3.
1748  real, parameter:: r4o3 = 4./3.
1749  real gak(km)
1750  real bet
1751  integer i, k
1752 
1753 !------------------------------------------------
1754 ! Optimized coding for uniform grid: SJL Apr 2007
1755 !------------------------------------------------
1756 
1757  if ( id==1 ) then
1758  do i=i1,i2
1759  qe(i,1) = r4o3*q1(i,1) + r2o3*q1(i,2)
1760  enddo
1761  else
1762  do i=i1,i2
1763  qe(i,1) = 1.e30
1764  enddo
1765  endif
1766 
1767  gak(1) = 7./3.
1768  do k=2,km
1769  gak(k) = 1. / (4. - gak(k-1))
1770  do i=i1,i2
1771  qe(i,k) = (3.*(q1(i,k-1) + q1(i,k)) - qe(i,k-1)) * gak(k)
1772  enddo
1773  enddo
1774 
1775  bet = 1. / (1.5 - 3.5*gak(km))
1776  do i=i1,i2
1777  qe(i,km+1) = (4.*q1(i,km) + q1(i,km-1) - 3.5*qe(i,km)) * bet
1778  enddo
1779 
1780  do k=km,1,-1
1781  do i=i1,i2
1782  qe(i,k) = qe(i,k) - gak(k)*qe(i,k+1)
1783  enddo
1784  enddo
1785 
1786  end subroutine edge_scalar
1787 
1788 
1789 
1790  subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter)
1791 ! Optimized for wind profile reconstruction:
1792  integer, intent(in):: i1, i2, j1, j2
1793  integer, intent(in):: j, km
1794  integer, intent(in):: limiter
1795  logical, intent(in):: uniform_grid
1796  real, intent(in):: dp0(km)
1797  real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2
1798  real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e
1799 !-----------------------------------------------------------------------
1800  real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values
1801  real gak(km)
1802  real bet, r2o3, r4o3
1803  real g0, gk, xt1, xt2, a_bot
1804  integer i, k
1805 
1806  if ( uniform_grid ) then
1807 !------------------------------------------------
1808 ! Optimized coding for uniform grid: SJL Apr 2007
1809 !------------------------------------------------
1810  r2o3 = 2./3.
1811  r4o3 = 4./3.
1812  do i=i1,i2
1813  qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2)
1814  qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2)
1815  enddo
1816 
1817  gak(1) = 7./3.
1818  do k=2,km
1819  gak(k) = 1. / (4. - gak(k-1))
1820  do i=i1,i2
1821  qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k)
1822  qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k)
1823  enddo
1824  enddo
1825 
1826  bet = 1. / (1.5 - 3.5*gak(km))
1827  do i=i1,i2
1828  qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet
1829  qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet
1830  enddo
1831 
1832  do k=km,1,-1
1833  do i=i1,i2
1834  qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1)
1835  qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1)
1836  enddo
1837  enddo
1838  else
1839 ! Assuming grid varying in vertical only
1840  g0 = dp0(2) / dp0(1)
1841  xt1 = 2.*g0*(g0+1. )
1842  bet = g0*(g0+0.5)
1843  do i=i1,i2
1844  qe1(i,1) = ( xt1*q1(i,j,1) + q1(i,j,2) ) / bet
1845  qe2(i,1) = ( xt1*q2(i,j,1) + q2(i,j,2) ) / bet
1846  gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet
1847  enddo
1848 
1849  do k=2,km
1850  gk = dp0(k-1) / dp0(k)
1851  do i=i1,i2
1852  bet = 2. + 2.*gk - gam(i,k-1)
1853  qe1(i,k) = ( 3.*(q1(i,j,k-1)+gk*q1(i,j,k)) - qe1(i,k-1) ) / bet
1854  qe2(i,k) = ( 3.*(q2(i,j,k-1)+gk*q2(i,j,k)) - qe2(i,k-1) ) / bet
1855  gam(i,k) = gk / bet
1856  enddo
1857  enddo
1858 
1859  a_bot = 1. + gk*(gk+1.5)
1860  xt1 = 2.*gk*(gk+1.)
1861  do i=i1,i2
1862  xt2 = gk*(gk+0.5) - a_bot*gam(i,km)
1863  qe1(i,km+1) = ( xt1*q1(i,j,km) + q1(i,j,km-1) - a_bot*qe1(i,km) ) / xt2
1864  qe2(i,km+1) = ( xt1*q2(i,j,km) + q2(i,j,km-1) - a_bot*qe2(i,km) ) / xt2
1865  enddo
1866 
1867  do k=km,1,-1
1868  do i=i1,i2
1869  qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1)
1870  qe2(i,k) = qe2(i,k) - gam(i,k)*qe2(i,k+1)
1871  enddo
1872  enddo
1873  endif
1874 
1875 !------------------
1876 ! Apply constraints
1877 !------------------
1878  if ( limiter/=0 ) then ! limit the top & bottom winds
1879  do i=i1,i2
1880 ! Top
1881  if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0.
1882  if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0.
1883 ! Surface:
1884  if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0.
1885  if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0.
1886  enddo
1887  endif
1888 
1889  do k=1,km+1
1890  do i=i1,i2
1891  q1e(i,j,k) = qe1(i,k)
1892  q2e(i,j,k) = qe2(i,k)
1893  enddo
1894  enddo
1895 
1896  end subroutine edge_profile
1897 
1898  subroutine nest_halo_nh(ptop, grav, kappa, cp, delp, delz, pt, phis, &
1899 #ifdef MULTI_GASES
1900  q , &
1901 #endif
1902 #ifdef USE_COND
1903  q_con, &
1904 #ifdef MOIST_CAPPA
1905  cappa, &
1906 #endif
1907 #endif
1908  pkc, gz, pk3, &
1909  npx, npy, npz, nested, pkc_pertn, computepk3, fullhalo, bd, regional)
1910 
1911  !INPUT: delp, delz, pt
1912  !OUTPUT: gz, pkc, pk3 (optional)
1913  integer, intent(IN) :: npx, npy, npz
1914  logical, intent(IN) :: pkc_pertn, computepk3, fullhalo, nested, regional
1915  real, intent(IN) :: ptop, kappa, cp, grav
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, delz
1919 #ifdef MULTI_GASES
1920  real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,*):: q
1921 #endif
1922 #ifdef USE_COND
1923  real, intent(IN), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: q_con
1924 #ifdef MOIST_CAPPA
1925  real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: cappa
1926 #endif
1927 #endif
1928  real, intent(INOUT), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1):: gz, pkc, pk3
1929 
1930  integer :: i,j,k
1931  real :: gama !'gamma'
1932  real :: ptk, rgrav, rkap, peln1, rdg
1933 
1934  real, dimension(bd%isd:bd%ied, npz+1, bd%jsd:bd%jed ) :: pe, peln
1935 #ifdef USE_COND
1936  real, dimension(bd%isd:bd%ied, npz+1 ) :: peg, pelng
1937 #endif
1938  real, dimension(bd%isd:bd%ied, npz) :: gam, bb, dd, pkz
1939  real, dimension(bd%isd:bd%ied, npz-1) :: g_rat
1940  real, dimension(bd%isd:bd%ied) :: bet
1941  real :: pm
1942 #ifdef MULTI_GASES
1943  real gamax
1944 #endif
1945 
1946  integer :: ifirst, ilast, jfirst, jlast
1947 
1948  integer :: is, ie, js, je
1949  integer :: isd, ied, jsd, jed
1950 
1951  is = bd%is
1952  ie = bd%ie
1953  js = bd%js
1954  je = bd%je
1955  isd = bd%isd
1956  ied = bd%ied
1957  jsd = bd%jsd
1958  jed = bd%jed
1959 
1960  if (.not. (nested .or. regional)) return
1961  ifirst = isd
1962  jfirst = jsd
1963  ilast = ied
1964  jlast = jed
1965 
1966  !Remember we want to compute these in the HALO. Note also this routine
1967  !requires an appropriate
1968 
1969  rgrav = 1./grav
1970  gama = 1./(1.-kappa)
1971  ptk = ptop ** kappa
1972  rkap = 1./kappa
1973  peln1 = log(ptop)
1974  rdg = - rdgas * rgrav
1975 
1976  !NOTE: Compiler does NOT like this sort of nested-grid BC code. Is it trying to do some ugly optimization?
1977 
1978  if (is == 1) then
1979 
1980  do j=jfirst,jlast
1981 
1982  !GZ
1983  do i=ifirst,0
1984  gz(i,j,npz+1) = phis(i,j)
1985  enddo
1986  do k=npz,1,-1
1987  do i=ifirst,0
1988  gz(i,j,k) = gz(i,j,k+1) - delz(i,j,k)*grav
1989  enddo
1990  enddo
1991 
1992  !Hydrostatic interface pressure
1993  do i=ifirst,0
1994  pe(i,1,j) = ptop
1995  peln(i,1,j) = peln1
1996 #ifdef USE_COND
1997  peg(i,1) = ptop
1998  pelng(i,1) = peln1
1999 #endif
2000  enddo
2001  do k=2,npz+1
2002  do i=ifirst,0
2003  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2004  peln(i,k,j) = log(pe(i,k,j))
2005 #ifdef USE_COND
2006  peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1))
2007  pelng(i,k) = log(peg(i,k))
2008 #endif
2009  enddo
2010  enddo
2011 
2012  !Perturbation nonhydro layer-mean pressure (NOT to the kappa)
2013  do k=1,npz
2014  do i=ifirst,0
2015  !Full p
2016 #ifdef MOIST_CAPPA
2017  pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
2018 #else
2019 #ifdef MULTI_GASES
2020  gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:)))
2021  pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2022 #else
2023  pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2024 #endif
2025 #endif
2026 ! !hydro
2027 #ifdef USE_COND
2028  pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
2029 #else
2030  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2031 #endif
2032  !Remove hydro cell-mean pressure
2033  pkz(i,k) = pkz(i,k) - pm
2034  enddo
2035  enddo
2036 
2037  !pressure solver
2038  do k=1,npz-1
2039  do i=ifirst,0
2040  g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1)
2041  bb(i,k) = 2.*(1. + g_rat(i,k))
2042  dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1))
2043  enddo
2044  enddo
2045 
2046  do i=ifirst,0
2047  bet(i) = bb(i,1)
2048  pkc(i,j,1) = 0.
2049  pkc(i,j,2) = dd(i,1)/bet(i)
2050  bb(i,npz) = 2.
2051  dd(i,npz) = 3.*pkz(i,npz)
2052  enddo
2053  do k=2,npz
2054  do i=ifirst,0
2055  gam(i,k) = g_rat(i,k-1)/bet(i)
2056  bet(i) = bb(i,k) - gam(i,k)
2057  pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i)
2058  enddo
2059  enddo
2060  do k=npz,2,-1
2061  do i=ifirst,0
2062  pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1)
2063 #ifdef NHNEST_DEBUG
2064  if (abs(pkc(i,j,k)) > 1.e5) then
2065  print*, mpp_pe(), i,j,k, 'PKC: ', pkc(i,j,k)
2066  endif
2067 #endif
2068  enddo
2069  enddo
2070 
2071  enddo
2072 
2073  do j=jfirst,jlast
2074 
2075  if (.not. pkc_pertn) then
2076  do k=npz+1,1,-1
2077  do i=ifirst,0
2078  pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j)
2079  enddo
2080  enddo
2081  endif
2082 
2083  !pk3 if necessary; doesn't require condenstate loading calculation
2084  if (computepk3) then
2085  do i=ifirst,0
2086  pk3(i,j,1) = ptk
2087  enddo
2088  do k=2,npz+1
2089  do i=ifirst,0
2090  pk3(i,j,k) = exp(kappa*log(pe(i,k,j)))
2091  enddo
2092  enddo
2093  endif
2094 
2095  enddo
2096 
2097  endif
2098 
2099  if (ie == npx-1) then
2100 
2101  do j=jfirst,jlast
2102 
2103  !GZ
2104  do i=npx,ilast
2105  gz(i,j,npz+1) = phis(i,j)
2106  enddo
2107  do k=npz,1,-1
2108  do i=npx,ilast
2109  gz(i,j,k) = gz(i,j,k+1) - delz(i,j,k)*grav
2110  enddo
2111  enddo
2112 
2113  !Hydrostatic interface pressure
2114  do i=npx,ilast
2115  pe(i,1,j) = ptop
2116  peln(i,1,j) = peln1
2117 #ifdef USE_COND
2118  peg(i,1) = ptop
2119  pelng(i,1) = peln1
2120 #endif
2121  enddo
2122  do k=2,npz+1
2123  do i=npx,ilast
2124  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2125  peln(i,k,j) = log(pe(i,k,j))
2126 #ifdef USE_COND
2127  peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1))
2128  pelng(i,k) = log(peg(i,k))
2129 #endif
2130  enddo
2131  enddo
2132 
2133  !Perturbation nonhydro layer-mean pressure (NOT to the kappa)
2134  do k=1,npz
2135  do i=npx,ilast
2136  !Full p
2137 #ifdef MOIST_CAPPA
2138  pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
2139 #else
2140 #ifdef MULTI_GASES
2141  gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:)))
2142  pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2143 #else
2144  pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2145 #endif
2146 #endif
2147  !hydro
2148 #ifdef USE_COND
2149  pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
2150 #else
2151  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2152 #endif
2153  !Remove hydro cell-mean pressure
2154  pkz(i,k) = pkz(i,k) - pm
2155  enddo
2156  enddo
2157 
2158  !pressure solver
2159  do k=1,npz-1
2160  do i=npx,ilast
2161  g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1)
2162  bb(i,k) = 2.*(1. + g_rat(i,k))
2163  dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1))
2164  enddo
2165  enddo
2166 
2167  do i=npx,ilast
2168  bet(i) = bb(i,1)
2169  pkc(i,j,1) = 0.
2170  pkc(i,j,2) = dd(i,1)/bet(i)
2171  bb(i,npz) = 2.
2172  dd(i,npz) = 3.*pkz(i,npz)
2173  enddo
2174  do k=2,npz
2175  do i=npx,ilast
2176  gam(i,k) = g_rat(i,k-1)/bet(i)
2177  bet(i) = bb(i,k) - gam(i,k)
2178  pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i)
2179  enddo
2180  enddo
2181  do k=npz,2,-1
2182  do i=npx,ilast
2183  pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1)
2184  enddo
2185  enddo
2186 
2187 
2188  enddo
2189 
2190  do j=jfirst,jlast
2191 
2192  if (.not. pkc_pertn) then
2193  do k=npz+1,1,-1
2194  do i=npx,ilast
2195  pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j)
2196  enddo
2197  enddo
2198  endif
2199 
2200  !pk3 if necessary
2201  if (computepk3) then
2202  do i=npx,ilast
2203  pk3(i,j,1) = ptk
2204  enddo
2205  do k=2,npz+1
2206  do i=npx,ilast
2207  pk3(i,j,k) = exp(kappa*log(pe(i,k,j)))
2208  enddo
2209  enddo
2210  endif
2211 
2212  enddo
2213 
2214  endif
2215 
2216  if (js == 1) then
2217 
2218  do j=jfirst,0
2219 
2220  !GZ
2221  do i=ifirst,ilast
2222  gz(i,j,npz+1) = phis(i,j)
2223  enddo
2224  do k=npz,1,-1
2225  do i=ifirst,ilast
2226  gz(i,j,k) = gz(i,j,k+1) - delz(i,j,k)*grav
2227  enddo
2228  enddo
2229 
2230  !Hydrostatic interface pressure
2231  do i=ifirst,ilast
2232  pe(i,1,j) = ptop
2233  peln(i,1,j) = peln1
2234 #ifdef USE_COND
2235  peg(i,1) = ptop
2236  pelng(i,1) = peln1
2237 #endif
2238  enddo
2239  do k=2,npz+1
2240  do i=ifirst,ilast
2241  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2242  peln(i,k,j) = log(pe(i,k,j))
2243 #ifdef USE_COND
2244  peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1))
2245  pelng(i,k) = log(peg(i,k))
2246 #endif
2247  enddo
2248  enddo
2249 
2250  !Perturbation nonhydro layer-mean pressure (NOT to the kappa)
2251  do k=1,npz
2252  do i=ifirst,ilast
2253  !Full p
2254 #ifdef MOIST_CAPPA
2255  pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
2256 #else
2257 #ifdef MULTI_GASES
2258  gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:)))
2259  pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2260 #else
2261  pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2262 #endif
2263 #endif
2264  !hydro
2265 #ifdef USE_COND
2266  pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
2267 #else
2268  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2269 #endif
2270  !hydro
2271  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2272  !Remove hydro cell-mean pressure
2273  pkz(i,k) = pkz(i,k) - pm
2274  enddo
2275  enddo
2276 
2277  !pressure solver
2278  do k=1,npz-1
2279  do i=ifirst,ilast
2280  g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1)
2281  bb(i,k) = 2.*(1. + g_rat(i,k))
2282  dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1))
2283  enddo
2284  enddo
2285 
2286  do i=ifirst,ilast
2287  bet(i) = bb(i,1)
2288  pkc(i,j,1) = 0.
2289  pkc(i,j,2) = dd(i,1)/bet(i)
2290  bb(i,npz) = 2.
2291  dd(i,npz) = 3.*pkz(i,npz)
2292  enddo
2293  do k=2,npz
2294  do i=ifirst,ilast
2295  gam(i,k) = g_rat(i,k-1)/bet(i)
2296  bet(i) = bb(i,k) - gam(i,k)
2297  pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i)
2298  enddo
2299  enddo
2300  do k=npz,2,-1
2301  do i=ifirst,ilast
2302  pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1)
2303 #ifdef NHNEST_DEBUG
2304  if (abs(pkc(i,j,k)) > 1.e5) then
2305  print*, mpp_pe(), i,j,k, 'PKC: ', pkc(i,j,k)
2306  endif
2307 #endif
2308  enddo
2309  enddo
2310 
2311  enddo
2312 
2313  do j=jfirst,0
2314 
2315  if (.not. pkc_pertn) then
2316  do k=npz+1,1,-1
2317  do i=ifirst,ilast
2318  pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j)
2319  enddo
2320  enddo
2321  endif
2322 
2323  !pk3 if necessary
2324  if (computepk3) then
2325  do i=ifirst,ilast
2326  pk3(i,j,1) = ptk
2327  enddo
2328  do k=2,npz+1
2329  do i=ifirst,ilast
2330  pk3(i,j,k) = exp(kappa*log(pe(i,k,j)))
2331  enddo
2332  enddo
2333  endif
2334 
2335  enddo
2336 
2337  endif
2338 
2339  if (je == npy-1) then
2340 
2341  do j=npy,jlast
2342 
2343  !GZ
2344  do i=ifirst,ilast
2345  gz(i,j,npz+1) = phis(i,j)
2346  enddo
2347  do k=npz,1,-1
2348  do i=ifirst,ilast
2349  gz(i,j,k) = gz(i,j,k+1) - delz(i,j,k)*grav
2350  enddo
2351  enddo
2352 
2353  !Hydrostatic interface pressure
2354  do i=ifirst,ilast
2355  pe(i,1,j) = ptop
2356  peln(i,1,j) = peln1
2357 #ifdef USE_COND
2358  peg(i,1) = ptop
2359  pelng(i,1) = peln1
2360 #endif
2361  enddo
2362  do k=2,npz+1
2363  do i=ifirst,ilast
2364  pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
2365  peln(i,k,j) = log(pe(i,k,j))
2366 #ifdef USE_COND
2367  peg(i,k) = peg(i,k-1) + delp(i,j,k-1)*(1.-q_con(i,j,k-1))
2368  pelng(i,k) = log(peg(i,k))
2369 #endif
2370  enddo
2371  enddo
2372 
2373  !Perturbation nonhydro layer-mean pressure (NOT to the kappa)
2374  do k=1,npz
2375  do i=ifirst,ilast
2376  !Full p
2377 #ifdef MOIST_CAPPA
2378  pkz(i,k) = exp(1./(1.-cappa(i,j,k))*log(rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
2379 #else
2380 #ifdef MULTI_GASES
2381  gamax = gama * (vicpqd(q(i,j,k,:))/vicvqd(q(i,j,k,:)))
2382  pkz(i,k) = exp(gamax*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2383 #else
2384  pkz(i,k) = exp(gama*log(-delp(i,j,k)*rgrav/delz(i,j,k)*rdgas*pt(i,j,k)))
2385 #endif
2386 #endif
2387  !hydro
2388 #ifdef USE_COND
2389  pm = (peg(i,k+1)-peg(i,k))/(pelng(i,k+1)-pelng(i,k))
2390 #else
2391  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2392 #endif
2393  !hydro
2394  pm = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
2395  !Remove hydro cell-mean pressure
2396  pkz(i,k) = pkz(i,k) - pm
2397  enddo
2398  enddo
2399 
2400  !Reversible interpolation on layer NH pressure perturbation
2401  ! to recover lastge NH pressure perturbation
2402  do k=1,npz-1
2403  do i=ifirst,ilast
2404  g_rat(i,k) = delp(i,j,k)/delp(i,j,k+1)
2405  bb(i,k) = 2.*(1. + g_rat(i,k))
2406  dd(i,k) = 3.*(pkz(i,k) + g_rat(i,k)*pkz(i,k+1))
2407  enddo
2408  enddo
2409 
2410  do i=ifirst,ilast
2411  bet(i) = bb(i,1)
2412  pkc(i,j,1) = 0.
2413  pkc(i,j,2) = dd(i,1)/bet(i)
2414  bb(i,npz) = 2.
2415  dd(i,npz) = 3.*pkz(i,npz)
2416  enddo
2417  do k=2,npz
2418  do i=ifirst,ilast
2419  gam(i,k) = g_rat(i,k-1)/bet(i)
2420  bet(i) = bb(i,k) - gam(i,k)
2421  pkc(i,j,k+1) = (dd(i,k) - pkc(i,j,k))/bet(i)
2422  enddo
2423  enddo
2424  do k=npz,2,-1
2425  do i=ifirst,ilast
2426  pkc(i,j,k) = pkc(i,j,k) - gam(i,k)*pkc(i,j,k+1)
2427  enddo
2428  enddo
2429 
2430 
2431  enddo
2432 
2433  do j=npy,jlast
2434 
2435  if (.not. pkc_pertn) then
2436  do k=npz+1,1,-1
2437  do i=ifirst,ilast
2438  pkc(i,j,k) = pkc(i,j,k) + pe(i,k,j)
2439  enddo
2440  enddo
2441  endif
2442 
2443  !pk3 if necessary
2444  if (computepk3) then
2445  do i=ifirst,ilast
2446  pk3(i,j,1) = ptk
2447  enddo
2448  do k=2,npz+1
2449  do i=ifirst,ilast
2450  pk3(i,j,k) = exp(kappa*log(pe(i,k,j)))
2451  enddo
2452  enddo
2453  endif
2454 
2455  enddo
2456 
2457  endif
2458 
2459 end subroutine nest_halo_nh
2460 
2461 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:1550
multi_gases_mod::vicpqd
pure real function, public vicpqd(q)
Definition: multi_gases.F90:239
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:1012
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:120
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1137
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:483
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:1791
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:1586
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
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:3320
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:680
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:1379
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:1199
nh_utils_mod::nest_halo_nh
subroutine, public nest_halo_nh(ptop, grav, kappa, cp, delp, delz, pt, phis, pkc, gz, pk3, npx, npy, npz, nested, pkc_pertn, computepk3, fullhalo, bd, regional)
Definition: nh_utils.F90:1910
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, regional, mfx, mfy, mass, nord, damp_c)
The subroutine &#39;fv_tp_2d&#39; contains the FV advection scheme .
Definition: tp_core.F90:110
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:747
nh_utils_mod::edge_scalar
subroutine edge_scalar(q1, qe, i1, i2, km, id)
Definition: nh_utils.F90:1741
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, delz, ws, rdt, gridstruct, bd, lim_fac, regional)
Definition: nh_utils.F90:216
multi_gases_mod::vicvqd
pure real function, public vicvqd(q)
Definition: multi_gases.F90:280
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:338