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