FV3DYCORE  Version 1.1.0
fv_mapz.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 
26 
27 module fv_mapz_mod
28 
29 ! Modules Included:
30 ! <table>
31 ! <tr>
32 ! <th>Module Name</th>
33 ! <th>Functions Included</th>
34 ! </tr>
35 ! <table>
36 ! <tr>
37 ! <td>constants_mod</td>
38 ! <td>radius, pi=>pi_8, rvgas, rdgas, grav, hlv, hlf, cp_air, cp_vapor</td>
39 ! </tr>
40 ! <td>field_manager_mod</td>
41 ! <td>MODEL_ATMOS</td>
42 ! </tr>
43 ! <tr>
44 ! <td>fv_arrays_mod</td>
45 ! <td>fv_grid_type</td>
46 ! </tr>
47 ! <tr>
48 ! <td>fv_cmp_mod</td>
49 ! <td>qs_init, fv_sat_adj</td>
50 ! </tr>
51 ! <tr>
52 ! <td>fv_fill_mod</td>
53 ! <td>fillz</td>
54 ! </tr>
55 ! <tr>
56 ! <td>fv_grid_utils_mod</td>
57 ! <td>g_sum, ptop_min</td>
58 ! </tr>
59 ! <tr>
60 ! <td>fv_mp_mod</td>
61 ! <td>is_master</td>
62 ! </tr>
63 ! <tr>
64 ! <td>fv_timing_mod</td>
65 ! <td>timing_on, timing_off</td>
66 ! </tr>
67 ! <tr>
68 ! <td>fv_tracer2d_mod</td>
69 ! <td>tracer_2d, tracer_2d_1L, tracer_2d_nested</td>
70 ! </tr>
71 ! <tr>
72 ! <td>mpp_mod/td>
73 ! <td>NOTE, mpp_error, get_unit, mpp_root_pe, mpp_pe</td>
74 ! </tr>
75 ! <tr>
76 ! <td>mpp_domains_mod/td>
77 ! <td> mpp_update_domains, domain2d</td>
78 ! </tr>
79 ! <tr>
80 ! <td>tracer_manager_mod</td>
81 ! <td>get_tracer_index</td>
82 ! </tr>
83 ! </table>
84 
85  use constants_mod, only: radius, pi=>pi_8, rvgas, rdgas, grav, hlv, hlf, cp_air, cp_vapor
86  use tracer_manager_mod,only: get_tracer_index
87  use field_manager_mod, only: model_atmos
88  use fv_grid_utils_mod, only: g_sum, ptop_min
89  use fv_fill_mod, only: fillz
90  use mpp_domains_mod, only: mpp_update_domains, domain2d
91  use mpp_mod, only: note, fatal, mpp_error, get_unit, mpp_root_pe, mpp_pe
92  use fv_arrays_mod, only: fv_grid_type
93  use fv_timing_mod, only: timing_on, timing_off
94  use fv_mp_mod, only: is_master
95 #ifndef CCPP
96  use fv_cmp_mod, only: qs_init, fv_sat_adj
97 #else
98 #ifdef STATIC
99 ! For static builds, the ccpp_physics_{init,run,finalize} calls
100 ! are not pointing to code in the CCPP framework, but to auto-generated
101 ! ccpp_suite_cap and ccpp_group_*_cap modules behind a ccpp_static_api
102  use ccpp_api, only: ccpp_initialized
103  use ccpp_static_api, only: ccpp_physics_run
104  use ccpp_data, only: ccpp_suite
105 #else
106  use ccpp_api, only: ccpp_initialized, ccpp_physics_run
107 #endif
108  use ccpp_data, only: cdata => cdata_tile, ccpp_interstitial
109 #endif
110 #ifdef MULTI_GASES
111  use multi_gases_mod, only: virq, virqd, vicpqd, vicvqd, num_gas
112 #endif
113 
114  implicit none
115  real, parameter:: consv_min= 0.001
116  real, parameter:: t_min= 184.
117  real, parameter:: r3 = 1./3., r23 = 2./3., r12 = 1./12.
118  real, parameter:: cv_vap = 3.*rvgas
119  real, parameter:: cv_air = cp_air - rdgas
120 ! real, parameter:: c_ice = 2106. !< heat capacity of ice at 0.C
121  real, parameter:: c_ice = 1972.
122  real, parameter:: c_liq = 4.1855e+3
123 ! real, parameter:: c_liq = 4218. !< ECMWF-IFS
124  real, parameter:: cp_vap = cp_vapor
125  real, parameter:: tice = 273.16
126 
127  real(kind=4) :: e_flux = 0.
128  private
129 
130  public compute_total_energy, lagrangian_to_eulerian, moist_cv, moist_cp, &
131  rst_remap, mappm, e_flux
132 
133 contains
134 
135 
138  subroutine lagrangian_to_eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
139  mdt, pdt, km, is,ie,js,je, isd,ied,jsd,jed, &
140  nq, nwat, sphum, q_con, u, v, w, delz, pt, q, hs, r_vir, cp, &
141  akap, cappa, kord_mt, kord_wz, kord_tr, kord_tm, peln, te0_2d, &
142  ng, ua, va, omga, te, ws, fill, reproduce_sum, out_dt, dtdt, &
143  ptop, ak, bk, pfull, gridstruct, domain, do_sat_adj, &
144  hydrostatic, hybrid_z, do_omega, adiabatic, do_adiabatic_init)
145  logical, intent(in):: last_step
146  real, intent(in):: mdt
147  real, intent(in):: pdt
148  integer, intent(in):: km
149  integer, intent(in):: nq
150  integer, intent(in):: nwat
151  integer, intent(in):: sphum
152  integer, intent(in):: ng
153  integer, intent(in):: is,ie,isd,ied
154  integer, intent(in):: js,je,jsd,jed
155  integer, intent(in):: kord_mt
156  integer, intent(in):: kord_wz
157  integer, intent(in):: kord_tr(nq)
158  integer, intent(in):: kord_tm
159 
160  real, intent(in):: consv
161  real, intent(in):: r_vir
162  real, intent(in):: cp
163  real, intent(in):: akap
164  real, intent(in):: hs(isd:ied,jsd:jed)
165  real, intent(inout):: te0_2d(is:ie,js:je)
166  real, intent(in):: ws(is:ie,js:je)
167 
168  logical, intent(in):: do_sat_adj
169  logical, intent(in):: fill
170  logical, intent(in):: reproduce_sum
171  logical, intent(in):: do_omega, adiabatic, do_adiabatic_init
172  real, intent(in) :: ptop
173  real, intent(in) :: ak(km+1)
174  real, intent(in) :: bk(km+1)
175  real, intent(in):: pfull(km)
176  type(fv_grid_type), intent(IN), target :: gridstruct
177  type(domain2d), intent(INOUT) :: domain
178 
179 ! INPUT/OUTPUT
180  real, intent(inout):: pk(is:ie,js:je,km+1)
181  real, intent(inout):: q(isd:ied,jsd:jed,km,*)
182  real, intent(inout):: delp(isd:ied,jsd:jed,km)
183  real, intent(inout):: pe(is-1:ie+1,km+1,js-1:je+1)
184  real, intent(inout):: ps(isd:ied,jsd:jed)
185 
186 ! u-wind will be ghosted one latitude to the north upon exit
187  real, intent(inout):: u(isd:ied ,jsd:jed+1,km)
188  real, intent(inout):: v(isd:ied+1,jsd:jed ,km)
189  real, intent(inout):: w(isd: ,jsd: ,1:)
190  real, intent(inout):: pt(isd:ied ,jsd:jed ,km)
192  real, intent(inout), dimension(isd:,jsd:,1:)::delz, q_con, cappa
193  logical, intent(in):: hydrostatic
194  logical, intent(in):: hybrid_z
195  logical, intent(in):: out_dt
196 
197  real, intent(inout):: ua(isd:ied,jsd:jed,km)
198  real, intent(inout):: va(isd:ied,jsd:jed,km)
199  real, intent(inout):: omga(isd:ied,jsd:jed,km)
200  real, intent(inout):: peln(is:ie,km+1,js:je)
201  real, intent(inout):: dtdt(is:ie,js:je,km)
202  real, intent(out):: pkz(is:ie,js:je,km)
203  real, intent(out):: te(isd:ied,jsd:jed,km)
204 
205 ! !DESCRIPTION:
206 !
207 ! !REVISION HISTORY:
208 ! SJL 03.11.04: Initial version for partial remapping
209 !
210 !-----------------------------------------------------------------------
211 #ifdef CCPP
212  real, dimension(is:ie,js:je):: te_2d, zsum0, zsum1
213 #else
214  real, dimension(is:ie,js:je):: te_2d, zsum0, zsum1, dpln
215 #endif
216  real, dimension(is:ie,km) :: q2, dp2
217  real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis
218  real, dimension(is:ie+1,km+1):: pe0, pe3
219  real, dimension(is:ie):: gz, cvm, qv
220  real rcp, rg, rrg, bkh, dtmp, k1k
221 #ifndef CCPP
222  logical:: fast_mp_consv
223 #endif
224  integer:: i,j,k
225  integer:: kdelz
226 #ifdef CCPP
227  integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, iq, n, kp, k_next
228  integer :: ierr
229 #else
230  integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, iq, n, kmp, kp, k_next
231 #endif
232 
233 #ifdef CCPP
234  ccpp_associate: associate( fast_mp_consv => ccpp_interstitial%fast_mp_consv, &
235  kmp => ccpp_interstitial%kmp )
236 #endif
237 
238  k1k = rdgas/cv_air ! akap / (1.-akap) = rg/Cv=0.4
239  rg = rdgas
240  rcp = 1./ cp
241  rrg = -rdgas/grav
242 
243  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
244  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
245  rainwat = get_tracer_index(model_atmos, 'rainwat')
246  snowwat = get_tracer_index(model_atmos, 'snowwat')
247  graupel = get_tracer_index(model_atmos, 'graupel')
248  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
249 
250  if ( do_sat_adj ) then
251  fast_mp_consv = (.not.do_adiabatic_init) .and. consv>consv_min
252 #ifndef CCPP
253  do k=1,km
254  kmp = k
255  if ( pfull(k) > 10.e2 ) exit
256  enddo
257  call qs_init(kmp)
258 #endif
259  endif
260 
261 !$OMP parallel do default(none) shared(is,ie,js,je,km,pe,ptop,kord_tm,hydrostatic, &
262 !$OMP pt,pk,rg,peln,q,nwat,liq_wat,rainwat,ice_wat,snowwat, &
263 !$OMP graupel,q_con,sphum,cappa,r_vir,rcp,k1k,delp, &
264 !$OMP delz,akap,pkz,te,u,v,ps, gridstruct, last_step, &
265 !$OMP ak,bk,nq,isd,ied,jsd,jed,kord_tr,fill, adiabatic, &
266 #ifdef MULTI_GASES
267 !$OMP num_gas, &
268 #endif
269 !$OMP hs,w,ws,kord_wz,do_omega,omga,rrg,kord_mt,ua) &
270 !$OMP private(qv,gz,cvm,kp,k_next,bkh,dp2, &
271 !$OMP pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2)
272  do 1000 j=js,je+1
273 
274  do k=1,km+1
275  do i=is,ie
276  pe1(i,k) = pe(i,k,j)
277  enddo
278  enddo
279 
280  do i=is,ie
281  pe2(i, 1) = ptop
282  pe2(i,km+1) = pe(i,km+1,j)
283  enddo
284 
285  if ( j /= (je+1) ) then
286  if ( kord_tm < 0 ) then
287 ! Note: pt at this stage is Theta_v
288  if ( hydrostatic ) then
289 ! Transform virtual pt to virtual Temp
290  do k=1,km
291  do i=is,ie
292 #ifdef MULTI_GASES
293  pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(akap*(peln(i,k+1,j)-peln(i,k,j)))
294  pkz(i,j,k) = exp(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(pkz(i,j,k)))
295  pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
296 #else
297  pt(i,j,k) = pt(i,j,k)*(pk(i,j,k+1)-pk(i,j,k))/(akap*(peln(i,k+1,j)-peln(i,k,j)))
298 #endif
299  enddo
300  enddo
301  else
302 ! Transform "density pt" to "density temp"
303  do k=1,km
304 #ifdef MOIST_CAPPA
305  call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
306  ice_wat, snowwat, graupel, q, gz, cvm)
307  do i=is,ie
308  q_con(i,j,k) = gz(i)
309 #ifdef MULTI_GASES
310  cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/virq(q(i,j,k,1:num_gas)) )
311 #else
312  cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) )
313 #endif
314  pt(i,j,k) = pt(i,j,k)*exp(cappa(i,j,k)/(1.-cappa(i,j,k))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
315  enddo
316 #else
317  do i=is,ie
318 #ifdef MULTI_GASES
319  pt(i,j,k) = pt(i,j,k)*exp(k1k*(virqd(q(i,j,k,1:num_gas))/vicvqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
320 #else
321  pt(i,j,k) = pt(i,j,k)*exp(k1k*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
322 #endif
323 ! Using dry pressure for the definition of the virtual potential temperature
324 ! pt(i,j,k) = pt(i,j,k)*exp(k1k*log(rrg*(1.-q(i,j,k,sphum))*delp(i,j,k)/delz(i,j,k)* &
325 ! pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum))))
326  enddo
327 #endif
328  enddo
329  endif ! hydro test
330  elseif ( hydrostatic ) then
331  call pkez(km, is, ie, js, je, j, pe, pk, akap, peln, pkz, ptop)
332 ! Compute cp_air*Tm + KE
333  do k=1,km
334  do i=is,ie
335 #ifdef MULTI_GASES
336  pkz(i,j,k) = exp(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(pkz(i,j,k)))
337 #endif
338  te(i,j,k) = 0.25*gridstruct%rsin2(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 + &
339  v(i,j,k)**2+v(i+1,j,k)**2 - &
340  (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j)) &
341  + cp_air*pt(i,j,k)*pkz(i,j,k)
342  enddo
343  enddo
344  endif
345 
346  if ( .not. hydrostatic ) then
347  do k=1,km
348  do i=is,ie
349  delz(i,j,k) = -delz(i,j,k) / delp(i,j,k) ! ="specific volume"/grav
350  enddo
351  enddo
352  endif
353 
354 ! update ps
355  do i=is,ie
356  ps(i,j) = pe1(i,km+1)
357  enddo
358 !
359 ! Hybrid sigma-P coordinate:
360 !
361  do k=2,km
362  do i=is,ie
363  pe2(i,k) = ak(k) + bk(k)*pe(i,km+1,j)
364  enddo
365  enddo
366  do k=1,km
367  do i=is,ie
368  dp2(i,k) = pe2(i,k+1) - pe2(i,k)
369  enddo
370  enddo
371 
372 !------------
373 ! update delp
374 !------------
375  do k=1,km
376  do i=is,ie
377  delp(i,j,k) = dp2(i,k)
378  enddo
379  enddo
380 
381 !------------------
382 ! Compute p**Kappa
383 !------------------
384  do k=1,km+1
385  do i=is,ie
386  pk1(i,k) = pk(i,j,k)
387  enddo
388  enddo
389 
390  do i=is,ie
391  pn2(i, 1) = peln(i, 1,j)
392  pn2(i,km+1) = peln(i,km+1,j)
393  pk2(i, 1) = pk1(i, 1)
394  pk2(i,km+1) = pk1(i,km+1)
395  enddo
396 
397  do k=2,km
398  do i=is,ie
399  pn2(i,k) = log(pe2(i,k))
400  pk2(i,k) = exp(akap*pn2(i,k))
401  enddo
402  enddo
403 
404  if ( kord_tm<0 ) then
405 !----------------------------------
406 ! Map t using logp
407 !----------------------------------
408  call map_scalar(km, peln(is,1,j), pt, gz, &
409  km, pn2, pt, &
410  is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm), t_min)
411  else
412 ! Map pt using pe
413  call map1_ppm (km, pe1, pt, gz, &
414  km, pe2, pt, &
415  is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm))
416  endif
417 
418 !----------------
419 ! Map constituents
420 !----------------
421  if( nq > 5 ) then
422  call mapn_tracer(nq, km, pe1, pe2, q, dp2, kord_tr, j, &
423  is, ie, isd, ied, jsd, jed, 0., fill)
424  elseif ( nq > 0 ) then
425 ! Remap one tracer at a time
426  do iq=1,nq
427  call map1_q2(km, pe1, q(isd,jsd,1,iq), &
428  km, pe2, q2, dp2, &
429  is, ie, 0, kord_tr(iq), j, isd, ied, jsd, jed, 0.)
430  if (fill) call fillz(ie-is+1, km, 1, q2, dp2)
431  do k=1,km
432  do i=is,ie
433  q(i,j,k,iq) = q2(i,k)
434  enddo
435  enddo
436  enddo
437  endif
438 
439  if ( .not. hydrostatic ) then
440 ! Remap vertical wind:
441  call map1_ppm (km, pe1, w, ws(is,j), &
442  km, pe2, w, &
443  is, ie, j, isd, ied, jsd, jed, -2, kord_wz)
444 ! Remap delz for hybrid sigma-p coordinate
445  call map1_ppm (km, pe1, delz, gz, &
446  km, pe2, delz, &
447  is, ie, j, isd, ied, jsd, jed, 1, abs(kord_tm))
448  do k=1,km
449  do i=is,ie
450  delz(i,j,k) = -delz(i,j,k)*dp2(i,k)
451  enddo
452  enddo
453  endif
454 
455 !----------
456 ! Update pk
457 !----------
458  do k=1,km+1
459  do i=is,ie
460  pk(i,j,k) = pk2(i,k)
461  enddo
462  enddo
463 
464 !----------------
465  if ( do_omega ) then
466 ! Start do_omega
467 ! Copy omega field to pe3
468  do i=is,ie
469  pe3(i,1) = 0.
470  enddo
471  do k=2,km+1
472  do i=is,ie
473  pe3(i,k) = omga(i,j,k-1)
474  enddo
475  enddo
476  endif
477 
478  do k=1,km+1
479  do i=is,ie
480  pe0(i,k) = peln(i,k,j)
481  peln(i,k,j) = pn2(i,k)
482  enddo
483  enddo
484 
485 !------------
486 ! Compute pkz
487 !hmhj pk is pe**kappa(=rgas/cp_air), but pkz=plyr**kappa(=r/cp)
488 !------------
489  if ( hydrostatic ) then
490  do k=1,km
491  do i=is,ie
492  pkz(i,j,k) = (pk2(i,k+1)-pk2(i,k))/(akap*(peln(i,k+1,j)-peln(i,k,j)))
493 #ifdef MULTI_GASES
494  pkz(i,j,k) = exp(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(pkz(i,j,k)))
495 #endif
496  enddo
497  enddo
498  else
499 ! Note: pt at this stage is T_v or T_m
500  do k=1,km
501 #ifdef MOIST_CAPPA
502  call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
503  ice_wat, snowwat, graupel, q, gz, cvm)
504  do i=is,ie
505  q_con(i,j,k) = gz(i)
506 #ifdef MULTI_GASES
507  cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/virq(q(i,j,k,1:num_gas)) )
508 #else
509  cappa(i,j,k) = rdgas / ( rdgas + cvm(i)/(1.+r_vir*q(i,j,k,sphum)) )
510 #endif
511  pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
512  enddo
513 #else
514  if ( kord_tm < 0 ) then
515  do i=is,ie
516 #ifdef MULTI_GASES
517  pkz(i,j,k) = exp(akap*virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
518 #else
519  pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
520 #endif
521 ! Using dry pressure for the definition of the virtual potential temperature
522 ! pkz(i,j,k) = exp(akap*log(rrg*(1.-q(i,j,k,sphum))*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum))))
523  enddo
524  else
525  do i=is,ie
526 #ifdef MULTI_GASES
527  pkz(i,j,k) = exp(k1k*virqd(q(i,j,k,1:num_gas))/vicvqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
528 #else
529  pkz(i,j,k) = exp(k1k*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
530 #endif
531 ! Using dry pressure for the definition of the virtual potential temperature
532 ! pkz(i,j,k) = exp(k1k*log(rrg*(1.-q(i,j,k,sphum))*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum))))
533  enddo
534  if ( last_step .and. (.not.adiabatic) ) then
535  do i=is,ie
536  pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
537  enddo
538  endif
539  endif
540 #endif
541  enddo
542  endif
543 
544 ! Interpolate omega/pe3 (defined at pe0) to remapped cell center (dp2)
545  if ( do_omega ) then
546  do k=1,km
547  do i=is,ie
548  dp2(i,k) = 0.5*(peln(i,k,j) + peln(i,k+1,j))
549  enddo
550  enddo
551  do i=is,ie
552  k_next = 1
553  do n=1,km
554  kp = k_next
555  do k=kp,km
556  if( dp2(i,n) <= pe0(i,k+1) .and. dp2(i,n) >= pe0(i,k) ) then
557  omga(i,j,n) = pe3(i,k) + (pe3(i,k+1) - pe3(i,k)) * &
558  (dp2(i,n)-pe0(i,k)) / (pe0(i,k+1)-pe0(i,k) )
559  k_next = k
560  exit
561  endif
562  enddo
563  enddo
564  enddo
565  endif ! end do_omega
566 
567  endif !(j < je+1)
568 
569  do i=is,ie+1
570  pe0(i,1) = pe(i,1,j)
571  enddo
572 !------
573 ! map u
574 !------
575  do k=2,km+1
576  do i=is,ie
577  pe0(i,k) = 0.5*(pe(i,k,j-1)+pe1(i,k))
578  enddo
579  enddo
580 
581  do k=1,km+1
582  bkh = 0.5*bk(k)
583  do i=is,ie
584  pe3(i,k) = ak(k) + bkh*(pe(i,km+1,j-1)+pe1(i,km+1))
585  enddo
586  enddo
587 
588  call map1_ppm( km, pe0(is:ie,:), u, gz, &
589  km, pe3(is:ie,:), u, &
590  is, ie, j, isd, ied, jsd, jed+1, -1, kord_mt)
591 
592  if (j < je+1) then
593 !------
594 ! map v
595 !------
596  do i=is,ie+1
597  pe3(i,1) = ak(1)
598  enddo
599 
600  do k=2,km+1
601  bkh = 0.5*bk(k)
602  do i=is,ie+1
603  pe0(i,k) = 0.5*(pe(i-1,k, j)+pe(i,k, j))
604  pe3(i,k) = ak(k) + bkh*(pe(i-1,km+1,j)+pe(i,km+1,j))
605  enddo
606  enddo
607 
608  call map1_ppm (km, pe0, v, gz, &
609  km, pe3, v, is, ie+1, &
610  j, isd, ied+1, jsd, jed, -1, kord_mt)
611  endif ! (j < je+1)
612 
613  do k=1,km
614  do i=is,ie
615  ua(i,j,k) = pe2(i,k+1)
616  enddo
617  enddo
618 
619 1000 continue
620 
621 #if defined(CCPP) && defined(__GFORTRAN__)
622 !$OMP parallel default(none) shared(is,ie,js,je,km,ptop,u,v,pe,ua,isd,ied,jsd,jed,kord_mt, &
623 !$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
624 !$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
625 !$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
626 !$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
627 !$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
628 !$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, &
629 !$OMP kord_tm,cdata,CCPP_interstitial) &
630 #ifdef STATIC
631 !$OMP shared(ccpp_suite) &
632 #endif
633 #ifdef MULTI_GASES
634 !$OMP shared(num_gas) &
635 #endif
636 !$OMP private(pe0,pe1,pe2,pe3,qv,cvm,gz,phis,kdelz,ierr)
637 #elif defined(CCPP)
638 !$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,isd,ied,jsd,jed,kord_mt, &
639 !$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
640 !$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
641 !$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
642 !$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
643 !$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
644 !$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, &
645 !$OMP fast_mp_consv,kord_tm,cdata, CCPP_interstitial) &
646 #ifdef STATIC
647 !$OMP shared(ccpp_suite) &
648 #endif
649 #ifdef MULTI_GASES
650 !$OMP shared(num_gas) &
651 #endif
652 !$OMP private(pe0,pe1,pe2,pe3,qv,cvm,gz,phis,kdelz,ierr)
653 #else
654 !$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,isd,ied,jsd,jed,kord_mt, &
655 !$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
656 !$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
657 !$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
658 !$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
659 !$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
660 !$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, &
661 !$OMP fast_mp_consv,kord_tm) &
662 #ifdef MULTI_GASES
663 !$OMP shared(num_gas) &
664 #endif
665 !$OMP private(pe0,pe1,pe2,pe3,qv,cvm,gz,phis,kdelz,dpln)
666 #endif
667 
668 !$OMP do
669  do k=2,km
670  do j=js,je
671  do i=is,ie
672  pe(i,k,j) = ua(i,j,k-1)
673  enddo
674  enddo
675  enddo
676 
677 dtmp = 0.
678 if( last_step .and. (.not.do_adiabatic_init) ) then
679 
680  if ( consv > consv_min ) then
681 
682 !$OMP do
683  do j=js,je
684  if ( hydrostatic ) then
685  do i=is,ie
686  gz(i) = hs(i,j)
687  do k=1,km
688  gz(i) = gz(i) + rg*pt(i,j,k)*(peln(i,k+1,j)-peln(i,k,j))
689  enddo
690  enddo
691  do i=is,ie
692  te_2d(i,j) = pe(i,km+1,j)*hs(i,j) - pe(i,1,j)*gz(i)
693  enddo
694 
695  do k=1,km
696  do i=is,ie
697  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cp*pt(i,j,k) + &
698  0.25*gridstruct%rsin2(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 + &
699  v(i,j,k)**2+v(i+1,j,k)**2 - &
700  (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j)))
701  enddo
702  enddo
703  else
704  do i=is,ie
705  te_2d(i,j) = 0.
706  phis(i,km+1) = hs(i,j)
707  enddo
708  do k=km,1,-1
709  do i=is,ie
710  phis(i,k) = phis(i,k+1) - grav*delz(i,j,k)
711  enddo
712  enddo
713 
714  do k=1,km
715 #ifdef MOIST_CAPPA
716  call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
717  ice_wat, snowwat, graupel, q, gz, cvm)
718  do i=is,ie
719 ! KE using 3D winds:
720  q_con(i,j,k) = gz(i)
721 #ifdef MULTI_GASES
722  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cvm(i)*pt(i,j,k)/virq(q(i,j,k,1:num_gas)) + &
723 #else
724  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cvm(i)*pt(i,j,k)/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i))) + &
725 #endif
726  0.5*(phis(i,k)+phis(i,k+1) + w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*( &
727  u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 - &
728  (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j))))
729  enddo
730 #else
731  do i=is,ie
732 #ifdef MULTI_GASES
733  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cv_air*pt(i,j,k)/virq(q(i,j,k,1:num_gas)) + &
734 #else
735  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cv_air*pt(i,j,k)/(1.+r_vir*q(i,j,k,sphum)) + &
736 #endif
737  0.5*(phis(i,k)+phis(i,k+1) + w(i,j,k)**2 + 0.5*gridstruct%rsin2(i,j)*( &
738  u(i,j,k)**2+u(i,j+1,k)**2 + v(i,j,k)**2+v(i+1,j,k)**2 - &
739  (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*gridstruct%cosa_s(i,j))))
740  enddo
741 #endif
742  enddo ! k-loop
743  endif ! end non-hydro
744 
745  do i=is,ie
746  te_2d(i,j) = te0_2d(i,j) - te_2d(i,j)
747  zsum1(i,j) = pkz(i,j,1)*delp(i,j,1)
748  enddo
749  do k=2,km
750  do i=is,ie
751  zsum1(i,j) = zsum1(i,j) + pkz(i,j,k)*delp(i,j,k)
752  enddo
753  enddo
754  if ( hydrostatic ) then
755  do i=is,ie
756  zsum0(i,j) = ptop*(pk(i,j,1)-pk(i,j,km+1)) + zsum1(i,j)
757  enddo
758  endif
759 
760  enddo ! j-loop
761 
762 !$OMP single
763  dtmp = consv*g_sum(domain, te_2d, is, ie, js, je, ng, gridstruct%area_64, 0, reproduce=.true.)
764  e_flux = dtmp / (grav*pdt*4.*pi*radius**2) ! unit: W/m**2
765  ! Note pdt is "phys" time step
766  if ( hydrostatic ) then
767  dtmp = dtmp / (cp* g_sum(domain, zsum0, is, ie, js, je, ng, gridstruct%area_64, 0, reproduce=.true.))
768  else
769  dtmp = dtmp / (cv_air*g_sum(domain, zsum1, is, ie, js, je, ng, gridstruct%area_64, 0, reproduce=.true.))
770  endif
771 !$OMP end single
772 
773  elseif ( consv < -consv_min ) then
774 
775 !$OMP do
776  do j=js,je
777  do i=is,ie
778  zsum1(i,j) = pkz(i,j,1)*delp(i,j,1)
779  enddo
780  do k=2,km
781  do i=is,ie
782  zsum1(i,j) = zsum1(i,j) + pkz(i,j,k)*delp(i,j,k)
783  enddo
784  enddo
785  if ( hydrostatic ) then
786  do i=is,ie
787  zsum0(i,j) = ptop*(pk(i,j,1)-pk(i,j,km+1)) + zsum1(i,j)
788  enddo
789  endif
790  enddo
791 
792  e_flux = consv
793 !$OMP single
794  if ( hydrostatic ) then
795  dtmp = e_flux*(grav*pdt*4.*pi*radius**2) / &
796  (cp*g_sum(domain, zsum0, is, ie, js, je, ng, gridstruct%area_64, 0, reproduce=.true.))
797  else
798  dtmp = e_flux*(grav*pdt*4.*pi*radius**2) / &
799  (cv_air*g_sum(domain, zsum1, is, ie, js, je, ng, gridstruct%area_64, 0, reproduce=.true.))
800  endif
801 
802 !$OMP end single
803  endif ! end consv check
804 endif ! end last_step check
805 
806 ! Note: pt at this stage is T_v
807 ! if ( (.not.do_adiabatic_init) .and. do_sat_adj ) then
808  if ( do_sat_adj ) then
809  call timing_on('sat_adj2')
810 #ifdef CCPP
811  if (ccpp_initialized(cdata)) then
812 #ifdef STATIC
813  call ccpp_physics_run(cdata, suite_name=trim(ccpp_suite), group_name='fast_physics', ierr=ierr)
814 #else
815  call ccpp_physics_run(cdata, group_name='fast_physics', ierr=ierr)
816 #endif
817  if (ierr/=0) call mpp_error(fatal, "Call to ccpp_physics_run for group 'fast_physics' failed")
818  else
819  call mpp_error (fatal, 'Lagrangian_to_Eulerian: can not call CCPP fast physics because cdata not initialized')
820  endif
821 #else
822 !$OMP do
823  do k=kmp,km
824  do j=js,je
825  do i=is,ie
826  dpln(i,j) = peln(i,k+1,j) - peln(i,k,j)
827  enddo
828  enddo
829  if (hydrostatic) then
830  kdelz = 1
831  else
832  kdelz = k
833  end if
834  call fv_sat_adj(abs(mdt), r_vir, is, ie, js, je, ng, hydrostatic, fast_mp_consv, &
835  te(isd,jsd,k), &
836 #ifdef MULTI_GASES
837  km, &
838 #endif
839  q(isd,jsd,k,sphum), q(isd,jsd,k,liq_wat), &
840  q(isd,jsd,k,ice_wat), q(isd,jsd,k,rainwat), &
841  q(isd,jsd,k,snowwat), q(isd,jsd,k,graupel), &
842  hs ,dpln, delz(isd:,jsd:,kdelz), pt(isd,jsd,k), delp(isd,jsd,k), q_con(isd:,jsd:,k), &
843 
844  cappa(isd:,jsd:,k), gridstruct%area_64, dtdt(is,js,k), out_dt, last_step, cld_amt>0, q(isd,jsd,k,cld_amt))
845 
846  if ( .not. hydrostatic ) then
847  do j=js,je
848  do i=is,ie
849 #ifdef MOIST_CAPPA
850  pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
851 #else
852 #ifdef MULTI_GASES
853  pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
854 #else
855  pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
856 #endif
857 #endif
858  enddo
859  enddo
860  endif
861  enddo ! OpenMP k-loop
862 
863 
864  if ( fast_mp_consv ) then
865 !$OMP do
866  do j=js,je
867  do i=is,ie
868  do k=kmp,km
869  te0_2d(i,j) = te0_2d(i,j) + te(i,j,k)
870  enddo
871  enddo
872  enddo
873  endif
874 #endif
875  call timing_off('sat_adj2')
876  endif ! do_sat_adj
877 
878 
879  if ( last_step ) then
880  ! Output temperature if last_step
881 !!! if ( is_master() ) write(*,*) 'dtmp=', dtmp, nwat
882 !$OMP do
883  do k=1,km
884  do j=js,je
885 #ifdef USE_COND
886  if ( nwat==2 ) then
887  do i=is,ie
888  gz(i) = max(0., q(i,j,k,liq_wat))
889  qv(i) = max(0., q(i,j,k,sphum))
890 #ifdef MULTI_GASES
891  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / virq(q(i,j,k,1:num_gas))
892 #else
893  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / ((1.+r_vir*qv(i))*(1.-gz(i)))
894 #endif
895  enddo
896  elseif ( nwat==6 ) then
897  do i=is,ie
898  gz(i) = q(i,j,k,liq_wat)+q(i,j,k,rainwat)+q(i,j,k,ice_wat)+q(i,j,k,snowwat)+q(i,j,k,graupel)
899 #ifdef MULTI_GASES
900  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/ virq(q(i,j,k,1:num_gas))
901 #else
902  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i)))
903 #endif
904  enddo
905  else
906  call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
907  ice_wat, snowwat, graupel, q, gz, cvm)
908  do i=is,ie
909 #ifdef MULTI_GASES
910  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / virq(q(i,j,k,1:num_gas))
911 #else
912  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / ((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i)))
913 #endif
914  enddo
915  endif
916 #else
917  if ( .not. adiabatic ) then
918  do i=is,ie
919 #ifdef MULTI_GASES
920  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / virq(q(i,j,k,1:num_gas))
921 #else
922  pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / (1.+r_vir*q(i,j,k,sphum))
923 #endif
924  enddo
925  endif
926 #endif
927  enddo ! j-loop
928  enddo ! k-loop
929  else ! not last_step
930  if ( kord_tm < 0 ) then
931 !$OMP do
932  do k=1,km
933  do j=js,je
934  do i=is,ie
935  pt(i,j,k) = pt(i,j,k)/pkz(i,j,k)
936  enddo
937  enddo
938  enddo
939  endif
940  endif
941 !$OMP end parallel
942 
943 #ifdef CCPP
944  end associate ccpp_associate
945 #endif
946 
947  end subroutine lagrangian_to_eulerian
948 
949 
952  subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
953  u, v, w, delz, pt, delp, q, qc, pe, peln, hs, &
954  rsin2_l, cosa_s_l, &
955  r_vir, cp, rg, hlv, te_2d, ua, va, teq, &
956  moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
957 !------------------------------------------------------
958 ! Compute vertically integrated total energy per column
959 !------------------------------------------------------
960 ! !INPUT PARAMETERS:
961  integer, intent(in):: km, is, ie, js, je, isd, ied, jsd, jed, id_te
962  integer, intent(in):: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, nwat
963  real, intent(inout), dimension(isd:ied,jsd:jed,km):: ua, va
964  real, intent(in), dimension(isd:ied,jsd:jed,km):: pt, delp
965  real, intent(in), dimension(isd:ied,jsd:jed,km,*):: q
966  real, intent(in), dimension(isd:ied,jsd:jed,km):: qc
967  real, intent(inout):: u(isd:ied, jsd:jed+1,km)
968  real, intent(inout):: v(isd:ied+1,jsd:jed, km)
969  real, intent(in):: w(isd:,jsd:,1:)
970  real, intent(in):: delz(isd:,jsd:,1:)
971  real, intent(in):: hs(isd:ied,jsd:jed)
972  real, intent(in):: pe(is-1:ie+1,km+1,js-1:je+1)
973  real, intent(in):: peln(is:ie,km+1,js:je)
974  real, intent(in):: cp, rg, r_vir, hlv
975  real, intent(in) :: rsin2_l(isd:ied, jsd:jed)
976  real, intent(in) :: cosa_s_l(isd:ied, jsd:jed)
977  logical, intent(in):: moist_phys, hydrostatic
978 !! Output:
979  real, intent(out):: te_2d(is:ie,js:je)
980  real, intent(out):: teq(is:ie,js:je)
982  real, dimension(is:ie,km):: tv
983  real phiz(is:ie,km+1)
984  real cvm(is:ie), qd(is:ie)
985  integer i, j, k
986 
987 !----------------------
988 ! Output lat-lon winds:
989 !----------------------
990 ! call cubed_to_latlon(u, v, ua, va, dx, dy, rdxa, rdya, km, flagstruct%c2l_ord)
991 
992 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,hydrostatic,hs,pt,qc,rg,peln,te_2d, &
993 !$OMP pe,delp,cp,rsin2_l,u,v,cosa_s_l,delz,moist_phys,w, &
994 #ifdef MULTI_GASES
995 !$OMP num_gas, &
996 #endif
997 !$OMP q,nwat,liq_wat,rainwat,ice_wat,snowwat,graupel,sphum) &
998 !$OMP private(phiz, tv, cvm, qd)
999  do j=js,je
1000 
1001  if ( hydrostatic ) then
1002 
1003  do i=is,ie
1004  phiz(i,km+1) = hs(i,j)
1005  enddo
1006  do k=km,1,-1
1007  do i=is,ie
1008  tv(i,k) = pt(i,j,k)*(1.+qc(i,j,k))
1009  phiz(i,k) = phiz(i,k+1) + rg*tv(i,k)*(peln(i,k+1,j)-peln(i,k,j))
1010  enddo
1011  enddo
1012 
1013  do i=is,ie
1014  te_2d(i,j) = pe(i,km+1,j)*phiz(i,km+1) - pe(i,1,j)*phiz(i,1)
1015  enddo
1016 
1017  do k=1,km
1018  do i=is,ie
1019  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*(cp*tv(i,k) + &
1020  0.25*rsin2_l(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 + &
1021  v(i,j,k)**2+v(i+1,j,k)**2 - &
1022  (u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*cosa_s_l(i,j)))
1023  enddo
1024  enddo
1025 
1026  else
1027 !-----------------
1028 ! Non-hydrostatic:
1029 !-----------------
1030  do i=is,ie
1031  phiz(i,km+1) = hs(i,j)
1032  do k=km,1,-1
1033  phiz(i,k) = phiz(i,k+1) - grav*delz(i,j,k)
1034  enddo
1035  enddo
1036  do i=is,ie
1037  te_2d(i,j) = 0.
1038  enddo
1039  if ( moist_phys ) then
1040  do k=1,km
1041 #ifdef MOIST_CAPPA
1042  call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
1043  ice_wat, snowwat, graupel, q, qd, cvm)
1044 #endif
1045  do i=is,ie
1046 #ifdef MOIST_CAPPA
1047  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*( cvm(i)*pt(i,j,k) + &
1048 #else
1049 #ifdef MULTI_GASES
1050  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*( cv_air*vicvqd(q(i,j,k,1:num_gas) )*pt(i,j,k) + &
1051 #else
1052  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*( cv_air*pt(i,j,k) + &
1053 #endif
1054 #endif
1055  0.5*(phiz(i,k)+phiz(i,k+1)+w(i,j,k)**2+0.5*rsin2_l(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 + &
1056  v(i,j,k)**2+v(i+1,j,k)**2-(u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*cosa_s_l(i,j))))
1057  enddo
1058  enddo
1059  else
1060  do k=1,km
1061  do i=is,ie
1062 #ifdef MULTI_GASES
1063  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*( cv_air*vicvqd(q(i,j,k,1:num_gas))*pt(i,j,k) + &
1064 #else
1065  te_2d(i,j) = te_2d(i,j) + delp(i,j,k)*( cv_air*pt(i,j,k) + &
1066 #endif
1067  0.5*(phiz(i,k)+phiz(i,k+1)+w(i,j,k)**2+0.5*rsin2_l(i,j)*(u(i,j,k)**2+u(i,j+1,k)**2 + &
1068  v(i,j,k)**2+v(i+1,j,k)**2-(u(i,j,k)+u(i,j+1,k))*(v(i,j,k)+v(i+1,j,k))*cosa_s_l(i,j))))
1069  enddo
1070  enddo
1071  endif
1072  endif
1073  enddo
1074 
1075 !-------------------------------------
1076 ! Diganostics computation for moist TE
1077 !-------------------------------------
1078  if( id_te>0 ) then
1079 !$OMP parallel do default(none) shared(is,ie,js,je,teq,te_2d,moist_phys,km,hlv,sphum,q,delp)
1080  do j=js,je
1081  do i=is,ie
1082  teq(i,j) = te_2d(i,j)
1083  enddo
1084  if ( moist_phys ) then
1085  do k=1,km
1086  do i=is,ie
1087  teq(i,j) = teq(i,j) + hlv*q(i,j,k,sphum)*delp(i,j,k)
1088  enddo
1089  enddo
1090  endif
1091  enddo
1092  endif
1093 
1094  end subroutine compute_total_energy
1095 
1096  subroutine pkez(km, ifirst, ilast, jfirst, jlast, j, &
1097  pe, pk, akap, peln, pkz, ptop)
1098 
1099 ! INPUT PARAMETERS:
1100  integer, intent(in):: km, j
1101  integer, intent(in):: ifirst, ilast
1102  integer, intent(in):: jfirst, jlast
1103  real, intent(in):: akap
1104  real, intent(in):: pe(ifirst-1:ilast+1,km+1,jfirst-1:jlast+1)
1105  real, intent(in):: pk(ifirst:ilast,jfirst:jlast,km+1)
1106  real, intent(IN):: ptop
1107 ! OUTPUT
1108  real, intent(out):: pkz(ifirst:ilast,jfirst:jlast,km)
1109  real, intent(inout):: peln(ifirst:ilast, km+1, jfirst:jlast)
1110 ! Local
1111  real pk2(ifirst:ilast, km+1)
1112  real pek
1113  real lnp
1114  real ak1
1115  integer i, k
1116 
1117  ak1 = (akap + 1.) / akap
1118 
1119  pek = pk(ifirst,j,1)
1120  do i=ifirst, ilast
1121  pk2(i,1) = pek
1122  enddo
1123 
1124  do k=2,km+1
1125  do i=ifirst, ilast
1126 ! peln(i,k,j) = log(pe(i,k,j))
1127  pk2(i,k) = pk(i,j,k)
1128  enddo
1129  enddo
1130 
1131 !---- GFDL modification
1132  if( ptop < ptop_min ) then
1133  do i=ifirst, ilast
1134  peln(i,1,j) = peln(i,2,j) - ak1
1135  enddo
1136  else
1137  lnp = log( ptop )
1138  do i=ifirst, ilast
1139  peln(i,1,j) = lnp
1140  enddo
1141  endif
1142 !---- GFDL modification
1143 
1144  do k=1,km
1145  do i=ifirst, ilast
1146  pkz(i,j,k) = (pk2(i,k+1) - pk2(i,k) ) / &
1147  (akap*(peln(i,k+1,j) - peln(i,k,j)) )
1148  enddo
1149  enddo
1150 
1151  end subroutine pkez
1152 
1153 
1154 
1155  subroutine remap_z(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord)
1156 
1157 ! INPUT PARAMETERS:
1158  integer, intent(in) :: i1
1159  integer, intent(in) :: i2
1160  integer, intent(in) :: kord
1161  integer, intent(in) :: km
1162  integer, intent(in) :: kn
1163  integer, intent(in) :: iv
1164 
1165  real, intent(in) :: pe1(i1:i2,km+1)
1166  real, intent(in) :: pe2(i1:i2,kn+1)
1167  real, intent(in) :: q1(i1:i2,km)
1168 
1169 ! INPUT/OUTPUT PARAMETERS:
1170  real, intent(inout):: q2(i1:i2,kn)
1171 
1172 ! LOCAL VARIABLES:
1173  real qs(i1:i2)
1174  real dp1( i1:i2,km)
1175  real q4(4,i1:i2,km)
1176  real pl, pr, qsum, delp, esl
1177  integer i, k, l, m, k0
1178 
1179  do k=1,km
1180  do i=i1,i2
1181  dp1(i,k) = pe1(i,k+1) - pe1(i,k) ! negative
1182  q4(1,i,k) = q1(i,k)
1183  enddo
1184  enddo
1185 
1186 ! Compute vertical subgrid distribution
1187  if ( kord >7 ) then
1188  call cs_profile( qs, q4, dp1, km, i1, i2, iv, kord )
1189  else
1190  call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
1191  endif
1192 
1193 ! Mapping
1194  do 1000 i=i1,i2
1195  k0 = 1
1196  do 555 k=1,kn
1197  do 100 l=k0,km
1198 ! locate the top edge: pe2(i,k)
1199  if(pe2(i,k) <= pe1(i,l) .and. pe2(i,k) >= pe1(i,l+1)) then
1200  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1201  if(pe2(i,k+1) >= pe1(i,l+1)) then
1202 ! entire new grid is within the original grid
1203  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1204  q2(i,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
1205  *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
1206  k0 = l
1207  goto 555
1208  else
1209 ! Fractional area...
1210  qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+ &
1211  q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)* &
1212  (r3*(1.+pl*(1.+pl))))
1213  do m=l+1,km
1214 ! locate the bottom edge: pe2(i,k+1)
1215  if(pe2(i,k+1) < pe1(i,m+1) ) then
1216 ! Whole layer..
1217  qsum = qsum + dp1(i,m)*q4(1,i,m)
1218  else
1219  delp = pe2(i,k+1)-pe1(i,m)
1220  esl = delp / dp1(i,m)
1221  qsum = qsum + delp*(q4(2,i,m)+0.5*esl* &
1222  (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.-r23*esl)))
1223  k0 = m
1224  goto 123
1225  endif
1226  enddo
1227  goto 123
1228  endif
1229  endif
1230 100 continue
1231 123 q2(i,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
1232 555 continue
1233 1000 continue
1234 
1235  end subroutine remap_z
1236 
1237  subroutine map_scalar( km, pe1, q1, qs, &
1238  kn, pe2, q2, i1, i2, &
1239  j, ibeg, iend, jbeg, jend, iv, kord, q_min)
1240 ! iv=1
1241  integer, intent(in) :: i1
1242  integer, intent(in) :: i2
1243  integer, intent(in) :: iv
1244  integer, intent(in) :: kord
1245  integer, intent(in) :: j
1246  integer, intent(in) :: ibeg, iend, jbeg, jend
1247  integer, intent(in) :: km
1248  integer, intent(in) :: kn
1249  real, intent(in) :: qs(i1:i2)
1250  real, intent(in) :: pe1(i1:i2,km+1)
1251  real, intent(in) :: pe2(i1:i2,kn+1)
1252  real, intent(in) :: q1(ibeg:iend,jbeg:jend,km)
1253 ! INPUT/OUTPUT PARAMETERS:
1254  real, intent(inout):: q2(ibeg:iend,jbeg:jend,kn)
1255  real, intent(in):: q_min
1256 
1257 ! DESCRIPTION:
1258 ! IV = 0: constituents
1259 ! pe1: pressure at layer edges (from model top to bottom surface)
1260 ! in the original vertical coordinate
1261 ! pe2: pressure at layer edges (from model top to bottom surface)
1262 ! in the new vertical coordinate
1263 ! LOCAL VARIABLES:
1264  real dp1(i1:i2,km)
1265  real q4(4,i1:i2,km)
1266  real pl, pr, qsum, dp, esl
1267  integer i, k, l, m, k0
1268 
1269  do k=1,km
1270  do i=i1,i2
1271  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
1272  q4(1,i,k) = q1(i,j,k)
1273  enddo
1274  enddo
1275 
1276 ! Compute vertical subgrid distribution
1277  if ( kord >7 ) then
1278  call scalar_profile( qs, q4, dp1, km, i1, i2, iv, kord, q_min )
1279  else
1280  call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
1281  endif
1282 
1283  do i=i1,i2
1284  k0 = 1
1285  do 555 k=1,kn
1286  do l=k0,km
1287 ! locate the top edge: pe2(i,k)
1288  if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
1289  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1290  if( pe2(i,k+1) <= pe1(i,l+1) ) then
1291 ! entire new grid is within the original grid
1292  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1293  q2(i,j,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
1294  *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
1295  k0 = l
1296  goto 555
1297  else
1298 ! Fractional area...
1299  qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+ &
1300  q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)* &
1301  (r3*(1.+pl*(1.+pl))))
1302  do m=l+1,km
1303 ! locate the bottom edge: pe2(i,k+1)
1304  if( pe2(i,k+1) > pe1(i,m+1) ) then
1305 ! Whole layer
1306  qsum = qsum + dp1(i,m)*q4(1,i,m)
1307  else
1308  dp = pe2(i,k+1)-pe1(i,m)
1309  esl = dp / dp1(i,m)
1310  qsum = qsum + dp*(q4(2,i,m)+0.5*esl* &
1311  (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.-r23*esl)))
1312  k0 = m
1313  goto 123
1314  endif
1315  enddo
1316  goto 123
1317  endif
1318  endif
1319  enddo
1320 123 q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
1321 555 continue
1322  enddo
1323 
1324  end subroutine map_scalar
1325 
1326 
1327  subroutine map1_ppm( km, pe1, q1, qs, &
1328  kn, pe2, q2, i1, i2, &
1329  j, ibeg, iend, jbeg, jend, iv, kord)
1330  integer, intent(in) :: i1
1331  integer, intent(in) :: i2
1332  integer, intent(in) :: iv
1333  integer, intent(in) :: kord
1334  integer, intent(in) :: j
1335  integer, intent(in) :: ibeg, iend, jbeg, jend
1336  integer, intent(in) :: km
1337  integer, intent(in) :: kn
1338  real, intent(in) :: qs(i1:i2)
1339  real, intent(in) :: pe1(i1:i2,km+1)
1340  real, intent(in) :: pe2(i1:i2,kn+1)
1341  real, intent(in) :: q1(ibeg:iend,jbeg:jend,km)
1342 ! INPUT/OUTPUT PARAMETERS:
1343  real, intent(inout):: q2(ibeg:iend,jbeg:jend,kn)
1344 
1345 ! DESCRIPTION:
1346 ! IV = 0: constituents
1347 ! pe1: pressure at layer edges (from model top to bottom surface)
1348 ! in the original vertical coordinate
1349 ! pe2: pressure at layer edges (from model top to bottom surface)
1350 ! in the new vertical coordinate
1351 
1352 ! LOCAL VARIABLES:
1353  real dp1(i1:i2,km)
1354  real q4(4,i1:i2,km)
1355  real pl, pr, qsum, dp, esl
1356  integer i, k, l, m, k0
1357 
1358  do k=1,km
1359  do i=i1,i2
1360  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
1361  q4(1,i,k) = q1(i,j,k)
1362  enddo
1363  enddo
1364 
1365 ! Compute vertical subgrid distribution
1366  if ( kord >7 ) then
1367  call cs_profile( qs, q4, dp1, km, i1, i2, iv, kord )
1368  else
1369  call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
1370  endif
1371 
1372  do i=i1,i2
1373  k0 = 1
1374  do 555 k=1,kn
1375  do l=k0,km
1376 ! locate the top edge: pe2(i,k)
1377  if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
1378  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1379  if( pe2(i,k+1) <= pe1(i,l+1) ) then
1380 ! entire new grid is within the original grid
1381  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1382  q2(i,j,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
1383  *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
1384  k0 = l
1385  goto 555
1386  else
1387 ! Fractional area...
1388  qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+ &
1389  q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)* &
1390  (r3*(1.+pl*(1.+pl))))
1391  do m=l+1,km
1392 ! locate the bottom edge: pe2(i,k+1)
1393  if( pe2(i,k+1) > pe1(i,m+1) ) then
1394 ! Whole layer
1395  qsum = qsum + dp1(i,m)*q4(1,i,m)
1396  else
1397  dp = pe2(i,k+1)-pe1(i,m)
1398  esl = dp / dp1(i,m)
1399  qsum = qsum + dp*(q4(2,i,m)+0.5*esl* &
1400  (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.-r23*esl)))
1401  k0 = m
1402  goto 123
1403  endif
1404  enddo
1405  goto 123
1406  endif
1407  endif
1408  enddo
1409 123 q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
1410 555 continue
1411  enddo
1412 
1413  end subroutine map1_ppm
1414 
1415 
1416  subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, &
1417  i1, i2, isd, ied, jsd, jed, q_min, fill)
1418 ! INPUT PARAMETERS:
1419  integer, intent(in):: km
1420  integer, intent(in):: j, nq, i1, i2
1421  integer, intent(in):: isd, ied, jsd, jed
1422  integer, intent(in):: kord(nq)
1423  real, intent(in):: pe1(i1:i2,km+1)
1424  real, intent(in):: pe2(i1:i2,km+1)
1425  real, intent(in):: dp2(i1:i2,km)
1426  real, intent(in):: q_min
1427  logical, intent(in):: fill
1428  real, intent(inout):: q1(isd:ied,jsd:jed,km,nq) ! Field input
1429 ! LOCAL VARIABLES:
1430  real:: q4(4,i1:i2,km,nq)
1431  real:: q2(i1:i2,km,nq)
1432  real:: qsum(nq)
1433  real:: dp1(i1:i2,km)
1434  real:: qs(i1:i2)
1435  real:: pl, pr, dp, esl, fac1, fac2
1436  integer:: i, k, l, m, k0, iq
1437 
1438  do k=1,km
1439  do i=i1,i2
1440  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
1441  enddo
1442  enddo
1443 
1444  do iq=1,nq
1445  do k=1,km
1446  do i=i1,i2
1447  q4(1,i,k,iq) = q1(i,j,k,iq)
1448  enddo
1449  enddo
1450  call scalar_profile( qs, q4(1,i1,1,iq), dp1, km, i1, i2, 0, kord(iq), q_min )
1451  enddo
1452 
1453 ! Mapping
1454  do 1000 i=i1,i2
1455  k0 = 1
1456  do 555 k=1,km
1457  do 100 l=k0,km
1458 ! locate the top edge: pe2(i,k)
1459  if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then
1460  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1461  if(pe2(i,k+1) <= pe1(i,l+1)) then
1462 ! entire new grid is within the original grid
1463  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1464  fac1 = pr + pl
1465  fac2 = r3*(pr*fac1 + pl*pl)
1466  fac1 = 0.5*fac1
1467  do iq=1,nq
1468  q2(i,k,iq) = q4(2,i,l,iq) + (q4(4,i,l,iq)+q4(3,i,l,iq)-q4(2,i,l,iq))*fac1 &
1469  - q4(4,i,l,iq)*fac2
1470  enddo
1471  k0 = l
1472  goto 555
1473  else
1474 ! Fractional area...
1475  dp = pe1(i,l+1) - pe2(i,k)
1476  fac1 = 1. + pl
1477  fac2 = r3*(1.+pl*fac1)
1478  fac1 = 0.5*fac1
1479  do iq=1,nq
1480  qsum(iq) = dp*(q4(2,i,l,iq) + (q4(4,i,l,iq)+ &
1481  q4(3,i,l,iq) - q4(2,i,l,iq))*fac1 - q4(4,i,l,iq)*fac2)
1482  enddo
1483  do m=l+1,km
1484 ! locate the bottom edge: pe2(i,k+1)
1485  if(pe2(i,k+1) > pe1(i,m+1) ) then
1486  ! Whole layer..
1487  do iq=1,nq
1488  qsum(iq) = qsum(iq) + dp1(i,m)*q4(1,i,m,iq)
1489  enddo
1490  else
1491  dp = pe2(i,k+1)-pe1(i,m)
1492  esl = dp / dp1(i,m)
1493  fac1 = 0.5*esl
1494  fac2 = 1.-r23*esl
1495  do iq=1,nq
1496  qsum(iq) = qsum(iq) + dp*( q4(2,i,m,iq) + fac1*( &
1497  q4(3,i,m,iq)-q4(2,i,m,iq)+q4(4,i,m,iq)*fac2 ) )
1498  enddo
1499  k0 = m
1500  goto 123
1501  endif
1502  enddo
1503  goto 123
1504  endif
1505  endif
1506 100 continue
1507 123 continue
1508  do iq=1,nq
1509  q2(i,k,iq) = qsum(iq) / dp2(i,k)
1510  enddo
1511 555 continue
1512 1000 continue
1513 
1514  if (fill) call fillz(i2-i1+1, km, nq, q2, dp2)
1515 
1516  do iq=1,nq
1517 ! if (fill) call fillz(i2-i1+1, km, 1, q2(i1,1,iq), dp2)
1518  do k=1,km
1519  do i=i1,i2
1520  q1(i,j,k,iq) = q2(i,k,iq)
1521  enddo
1522  enddo
1523  enddo
1524 
1525  end subroutine mapn_tracer
1526 
1527 
1528  subroutine map1_q2(km, pe1, q1, &
1529  kn, pe2, q2, dp2, &
1530  i1, i2, iv, kord, j, &
1531  ibeg, iend, jbeg, jend, q_min )
1532 
1533 
1534 ! INPUT PARAMETERS:
1535  integer, intent(in) :: j
1536  integer, intent(in) :: i1, i2
1537  integer, intent(in) :: ibeg, iend, jbeg, jend
1538  integer, intent(in) :: iv
1539  integer, intent(in) :: kord
1540  integer, intent(in) :: km
1541  integer, intent(in) :: kn
1542 
1543  real, intent(in) :: pe1(i1:i2,km+1)
1544  real, intent(in) :: pe2(i1:i2,kn+1)
1545  real, intent(in) :: q1(ibeg:iend,jbeg:jend,km)
1546  real, intent(in) :: dp2(i1:i2,kn)
1547  real, intent(in) :: q_min
1548 ! INPUT/OUTPUT PARAMETERS:
1549  real, intent(inout):: q2(i1:i2,kn)
1550 ! LOCAL VARIABLES:
1551  real qs(i1:i2)
1552  real dp1(i1:i2,km)
1553  real q4(4,i1:i2,km)
1554  real pl, pr, qsum, dp, esl
1555 
1556  integer i, k, l, m, k0
1557 
1558  do k=1,km
1559  do i=i1,i2
1560  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
1561  q4(1,i,k) = q1(i,j,k)
1562  enddo
1563  enddo
1564 
1565 ! Compute vertical subgrid distribution
1566  if ( kord >7 ) then
1567  call scalar_profile( qs, q4, dp1, km, i1, i2, iv, kord, q_min )
1568  else
1569  call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
1570  endif
1571 
1572 ! Mapping
1573  do 1000 i=i1,i2
1574  k0 = 1
1575  do 555 k=1,kn
1576  do 100 l=k0,km
1577 ! locate the top edge: pe2(i,k)
1578  if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then
1579  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1580  if(pe2(i,k+1) <= pe1(i,l+1)) then
1581 ! entire new grid is within the original grid
1582  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1583  q2(i,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
1584  *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
1585  k0 = l
1586  goto 555
1587  else
1588 ! Fractional area...
1589  qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+ &
1590  q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)* &
1591  (r3*(1.+pl*(1.+pl))))
1592  do m=l+1,km
1593 ! locate the bottom edge: pe2(i,k+1)
1594  if(pe2(i,k+1) > pe1(i,m+1) ) then
1595  ! Whole layer..
1596  qsum = qsum + dp1(i,m)*q4(1,i,m)
1597  else
1598  dp = pe2(i,k+1)-pe1(i,m)
1599  esl = dp / dp1(i,m)
1600  qsum = qsum + dp*(q4(2,i,m)+0.5*esl* &
1601  (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.-r23*esl)))
1602  k0 = m
1603  goto 123
1604  endif
1605  enddo
1606  goto 123
1607  endif
1608  endif
1609 100 continue
1610 123 q2(i,k) = qsum / dp2(i,k)
1611 555 continue
1612 1000 continue
1613 
1614  end subroutine map1_q2
1615 
1616 
1617 
1618  subroutine remap_2d(km, pe1, q1, &
1619  kn, pe2, q2, &
1620  i1, i2, iv, kord)
1621  integer, intent(in):: i1, i2
1622  integer, intent(in):: iv
1623  integer, intent(in):: kord
1624  integer, intent(in):: km
1625  integer, intent(in):: kn
1626  real, intent(in):: pe1(i1:i2,km+1)
1627  real, intent(in):: pe2(i1:i2,kn+1)
1628  real, intent(in) :: q1(i1:i2,km)
1629  real, intent(out):: q2(i1:i2,kn)
1630 ! LOCAL VARIABLES:
1631  real qs(i1:i2)
1632  real dp1(i1:i2,km)
1633  real q4(4,i1:i2,km)
1634  real pl, pr, qsum, dp, esl
1635  integer i, k, l, m, k0
1636 
1637  do k=1,km
1638  do i=i1,i2
1639  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
1640  q4(1,i,k) = q1(i,k)
1641  enddo
1642  enddo
1643 
1644 ! Compute vertical subgrid distribution
1645  if ( kord >7 ) then
1646  call cs_profile( qs, q4, dp1, km, i1, i2, iv, kord )
1647  else
1648  call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
1649  endif
1650 
1651  do i=i1,i2
1652  k0 = 1
1653  do 555 k=1,kn
1654 #ifdef OLD_TOP_EDGE
1655  if( pe2(i,k+1) <= pe1(i,1) ) then
1656 ! Entire grid above old ptop
1657  q2(i,k) = q4(2,i,1)
1658  elseif( pe2(i,k) < pe1(i,1) .and. pe2(i,k+1)>pe1(i,1) ) then
1659 ! Partially above old ptop:
1660  q2(i,k) = q1(i,1)
1661 #else
1662  if( pe2(i,k) <= pe1(i,1) ) then
1663 ! above old ptop:
1664  q2(i,k) = q1(i,1)
1665 #endif
1666  else
1667  do l=k0,km
1668 ! locate the top edge: pe2(i,k)
1669  if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
1670  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
1671  if(pe2(i,k+1) <= pe1(i,l+1)) then
1672 ! entire new grid is within the original grid
1673  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
1674  q2(i,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
1675  *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
1676  k0 = l
1677  goto 555
1678  else
1679 ! Fractional area...
1680  qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+ &
1681  q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)* &
1682  (r3*(1.+pl*(1.+pl))))
1683  do m=l+1,km
1684 ! locate the bottom edge: pe2(i,k+1)
1685  if(pe2(i,k+1) > pe1(i,m+1) ) then
1686  ! Whole layer..
1687  qsum = qsum + dp1(i,m)*q4(1,i,m)
1688  else
1689  dp = pe2(i,k+1)-pe1(i,m)
1690  esl = dp / dp1(i,m)
1691  qsum = qsum + dp*(q4(2,i,m)+0.5*esl* &
1692  (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.-r23*esl)))
1693  k0 = m
1694  goto 123
1695  endif
1696  enddo
1697  goto 123
1698  endif
1699  endif
1700  enddo
1701 123 q2(i,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
1702  endif
1703 555 continue
1704  enddo
1705 
1706  end subroutine remap_2d
1707 
1708 
1709  subroutine scalar_profile(qs, a4, delp, km, i1, i2, iv, kord, qmin)
1710 ! Optimized vertical profile reconstruction:
1711 ! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL
1712  integer, intent(in):: i1, i2
1713  integer, intent(in):: km
1714  integer, intent(in):: iv
1715  integer, intent(in):: kord
1716  real, intent(in) :: qs(i1:i2)
1717  real, intent(in) :: delp(i1:i2,km)
1718  real, intent(inout):: a4(4,i1:i2,km)
1719  real, intent(in):: qmin
1720 !-----------------------------------------------------------------------
1721  logical, dimension(i1:i2,km):: extm, ext5, ext6
1722  real gam(i1:i2,km)
1723  real q(i1:i2,km+1)
1724  real d4(i1:i2)
1725  real bet, a_bot, grat
1726  real pmp_1, lac_1, pmp_2, lac_2, x0, x1
1727  integer i, k, im
1728 
1729  if ( iv .eq. -2 ) then
1730  do i=i1,i2
1731  gam(i,2) = 0.5
1732  q(i,1) = 1.5*a4(1,i,1)
1733  enddo
1734  do k=2,km-1
1735  do i=i1, i2
1736  grat = delp(i,k-1) / delp(i,k)
1737  bet = 2. + grat + grat - gam(i,k)
1738  q(i,k) = (3.*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
1739  gam(i,k+1) = grat / bet
1740  enddo
1741  enddo
1742  do i=i1,i2
1743  grat = delp(i,km-1) / delp(i,km)
1744  q(i,km) = (3.*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / &
1745  (2. + grat + grat - gam(i,km))
1746  q(i,km+1) = qs(i)
1747  enddo
1748  do k=km-1,1,-1
1749  do i=i1,i2
1750  q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
1751  enddo
1752  enddo
1753  else
1754  do i=i1,i2
1755  grat = delp(i,2) / delp(i,1) ! grid ratio
1756  bet = grat*(grat+0.5)
1757  q(i,1) = ( (grat+grat)*(grat+1.)*a4(1,i,1) + a4(1,i,2) ) / bet
1758  gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
1759  enddo
1760 
1761  do k=2,km
1762  do i=i1,i2
1763  d4(i) = delp(i,k-1) / delp(i,k)
1764  bet = 2. + d4(i) + d4(i) - gam(i,k-1)
1765  q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
1766  gam(i,k) = d4(i) / bet
1767  enddo
1768  enddo
1769 
1770  do i=i1,i2
1771  a_bot = 1. + d4(i)*(d4(i)+1.5)
1772  q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
1773  / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
1774  enddo
1775 
1776  do k=km,1,-1
1777  do i=i1,i2
1778  q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
1779  enddo
1780  enddo
1781  endif
1782 
1783 !----- Perfectly linear scheme --------------------------------
1784  if ( abs(kord) > 16 ) then
1785  do k=1,km
1786  do i=i1,i2
1787  a4(2,i,k) = q(i,k )
1788  a4(3,i,k) = q(i,k+1)
1789  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1790  enddo
1791  enddo
1792  return
1793  endif
1794 !----- Perfectly linear scheme --------------------------------
1795 !------------------
1796 ! Apply constraints
1797 !------------------
1798  im = i2 - i1 + 1
1799 
1800 ! Apply *large-scale* constraints
1801  do i=i1,i2
1802  q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) )
1803  q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) )
1804  enddo
1805 
1806  do k=2,km
1807  do i=i1,i2
1808  gam(i,k) = a4(1,i,k) - a4(1,i,k-1)
1809  enddo
1810  enddo
1811 
1812 ! Interior:
1813  do k=3,km-1
1814  do i=i1,i2
1815  if ( gam(i,k-1)*gam(i,k+1)>0. ) then
1816 ! Apply large-scale constraint to ALL fields if not local max/min
1817  q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) )
1818  q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) )
1819  else
1820  if ( gam(i,k-1) > 0. ) then
1821 ! There exists a local max
1822  q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k)))
1823  else
1824 ! There exists a local min
1825  q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k)))
1826  if ( iv==0 ) q(i,k) = max(0., q(i,k))
1827  endif
1828  endif
1829  enddo
1830  enddo
1831 
1832 ! Bottom:
1833  do i=i1,i2
1834  q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) )
1835  q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) )
1836  enddo
1837 
1838  do k=1,km
1839  do i=i1,i2
1840  a4(2,i,k) = q(i,k )
1841  a4(3,i,k) = q(i,k+1)
1842  enddo
1843  enddo
1844 
1845  do k=1,km
1846  if ( k==1 .or. k==km ) then
1847  do i=i1,i2
1848  extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0.
1849  enddo
1850  else
1851  do i=i1,i2
1852  extm(i,k) = gam(i,k)*gam(i,k+1) < 0.
1853  enddo
1854  endif
1855  if ( abs(kord) > 9 ) then
1856  do i=i1,i2
1857  x0 = 2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))
1858  x1 = abs(a4(2,i,k)-a4(3,i,k))
1859  a4(4,i,k) = 3.*x0
1860  ext5(i,k) = abs(x0) > x1
1861  ext6(i,k) = abs(a4(4,i,k)) > x1
1862  enddo
1863  endif
1864  enddo
1865 
1866 !---------------------------
1867 ! Apply subgrid constraints:
1868 !---------------------------
1869 ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
1870 ! Top 2 and bottom 2 layers always use monotonic mapping
1871 
1872  if ( iv==0 ) then
1873  do i=i1,i2
1874  a4(2,i,1) = max(0., a4(2,i,1))
1875  enddo
1876  elseif ( iv==-1 ) then
1877  do i=i1,i2
1878  if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0.
1879  enddo
1880  elseif ( iv==2 ) then
1881  do i=i1,i2
1882  a4(2,i,1) = a4(1,i,1)
1883  a4(3,i,1) = a4(1,i,1)
1884  a4(4,i,1) = 0.
1885  enddo
1886  endif
1887 
1888  if ( iv/=2 ) then
1889  do i=i1,i2
1890  a4(4,i,1) = 3.*(2.*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1)))
1891  enddo
1892  call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1)
1893  endif
1894 
1895 ! k=2
1896  do i=i1,i2
1897  a4(4,i,2) = 3.*(2.*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2)))
1898  enddo
1899  call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2)
1900 
1901 !-------------------------------------
1902 ! Huynh's 2nd constraint for interior:
1903 !-------------------------------------
1904  do k=3,km-2
1905  if ( abs(kord)<9 ) then
1906  do i=i1,i2
1907 ! Left edges
1908  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
1909  lac_1 = pmp_1 + 1.5*gam(i,k+2)
1910  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1911  max(a4(1,i,k), pmp_1, lac_1) )
1912 ! Right edges
1913  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
1914  lac_2 = pmp_2 - 1.5*gam(i,k-1)
1915  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1916  max(a4(1,i,k), pmp_2, lac_2) )
1917 
1918  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1919  enddo
1920 
1921  elseif ( abs(kord)==9 ) then
1922  do i=i1,i2
1923  if ( extm(i,k) .and. extm(i,k-1) ) then
1924 ! grid-scale 2-delta-z wave detected
1925  a4(2,i,k) = a4(1,i,k)
1926  a4(3,i,k) = a4(1,i,k)
1927  a4(4,i,k) = 0.
1928  else if ( extm(i,k) .and. extm(i,k+1) ) then
1929 ! grid-scale 2-delta-z wave detected
1930  a4(2,i,k) = a4(1,i,k)
1931  a4(3,i,k) = a4(1,i,k)
1932  a4(4,i,k) = 0.
1933  else if ( extm(i,k) .and. a4(1,i,k)<qmin ) then
1934 ! grid-scale 2-delta-z wave detected
1935  a4(2,i,k) = a4(1,i,k)
1936  a4(3,i,k) = a4(1,i,k)
1937  a4(4,i,k) = 0.
1938  else
1939  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1940 ! Check within the smooth region if subgrid profile is non-monotonic
1941  if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
1942  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
1943  lac_1 = pmp_1 + 1.5*gam(i,k+2)
1944  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1945  max(a4(1,i,k), pmp_1, lac_1) )
1946  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
1947  lac_2 = pmp_2 - 1.5*gam(i,k-1)
1948  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1949  max(a4(1,i,k), pmp_2, lac_2) )
1950  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1951  endif
1952  endif
1953  enddo
1954  elseif ( abs(kord)==10 ) then
1955  do i=i1,i2
1956  if( ext5(i,k) ) then
1957  if( ext5(i,k-1) .or. ext5(i,k+1) ) then
1958  a4(2,i,k) = a4(1,i,k)
1959  a4(3,i,k) = a4(1,i,k)
1960  elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
1961  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
1962  lac_1 = pmp_1 + 1.5*gam(i,k+2)
1963  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1964  max(a4(1,i,k), pmp_1, lac_1) )
1965  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
1966  lac_2 = pmp_2 - 1.5*gam(i,k-1)
1967  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1968  max(a4(1,i,k), pmp_2, lac_2) )
1969  endif
1970  elseif( ext6(i,k) ) then
1971  if( ext5(i,k-1) .or. ext5(i,k+1) ) then
1972  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
1973  lac_1 = pmp_1 + 1.5*gam(i,k+2)
1974  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1975  max(a4(1,i,k), pmp_1, lac_1) )
1976  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
1977  lac_2 = pmp_2 - 1.5*gam(i,k-1)
1978  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1979  max(a4(1,i,k), pmp_2, lac_2) )
1980  endif
1981  endif
1982  enddo
1983  do i=i1,i2
1984  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1985  enddo
1986  elseif ( abs(kord)==12 ) then
1987  do i=i1,i2
1988  if( extm(i,k) ) then
1989  a4(2,i,k) = a4(1,i,k)
1990  a4(3,i,k) = a4(1,i,k)
1991  a4(4,i,k) = 0.
1992  else ! not a local extremum
1993  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
1994 ! Check within the smooth region if subgrid profile is non-monotonic
1995  if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
1996  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
1997  lac_1 = pmp_1 + 1.5*gam(i,k+2)
1998  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1999  max(a4(1,i,k), pmp_1, lac_1) )
2000  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2001  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2002  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2003  max(a4(1,i,k), pmp_2, lac_2) )
2004  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
2005  endif
2006  endif
2007  enddo
2008  elseif ( abs(kord)==13 ) then
2009  do i=i1,i2
2010  if( ext6(i,k) ) then
2011  if ( ext6(i,k-1) .and. ext6(i,k+1) ) then
2012 ! grid-scale 2-delta-z wave detected
2013  a4(2,i,k) = a4(1,i,k)
2014  a4(3,i,k) = a4(1,i,k)
2015  endif
2016  endif
2017  enddo
2018  do i=i1,i2
2019  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2020  enddo
2021  elseif ( abs(kord)==14 ) then
2022 
2023  do i=i1,i2
2024  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2025  enddo
2026 
2027  elseif ( abs(kord)==15 ) then ! Revised abs(kord)=9 scheme
2028  do i=i1,i2
2029  if ( ext5(i,k) .and. ext5(i,k-1) ) then
2030  a4(2,i,k) = a4(1,i,k)
2031  a4(3,i,k) = a4(1,i,k)
2032  else if ( ext5(i,k) .and. ext5(i,k+1) ) then
2033  a4(2,i,k) = a4(1,i,k)
2034  a4(3,i,k) = a4(1,i,k)
2035  else if ( ext5(i,k) .and. a4(1,i,k)<qmin ) then
2036  a4(2,i,k) = a4(1,i,k)
2037  a4(3,i,k) = a4(1,i,k)
2038  elseif( ext6(i,k) ) then
2039  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2040  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2041  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2042  max(a4(1,i,k), pmp_1, lac_1) )
2043  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2044  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2045  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2046  max(a4(1,i,k), pmp_2, lac_2) )
2047  endif
2048  enddo
2049  do i=i1,i2
2050  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2051  enddo
2052  elseif ( abs(kord)==16 ) then
2053  do i=i1,i2
2054  if( ext5(i,k) ) then
2055  if ( ext5(i,k-1) .or. ext5(i,k+1) ) then
2056  a4(2,i,k) = a4(1,i,k)
2057  a4(3,i,k) = a4(1,i,k)
2058  elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
2059  ! Left edges
2060  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2061  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2062  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2063  max(a4(1,i,k), pmp_1, lac_1) )
2064  ! Right edges
2065  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2066  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2067  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2068  max(a4(1,i,k), pmp_2, lac_2) )
2069  endif
2070  endif
2071  enddo
2072  do i=i1,i2
2073  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2074  enddo
2075  else ! kord = 11, 13
2076  do i=i1,i2
2077  if ( ext5(i,k) .and. (ext5(i,k-1).or.ext5(i,k+1).or.a4(1,i,k)<qmin) ) then
2078 ! Noisy region:
2079  a4(2,i,k) = a4(1,i,k)
2080  a4(3,i,k) = a4(1,i,k)
2081  a4(4,i,k) = 0.
2082  else
2083  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2084  endif
2085  enddo
2086  endif
2087 
2088 ! Additional constraint to ensure positivity
2089  if ( iv==0 ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 0)
2090 
2091  enddo ! k-loop
2092 
2093 !----------------------------------
2094 ! Bottom layer subgrid constraints:
2095 !----------------------------------
2096  if ( iv==0 ) then
2097  do i=i1,i2
2098  a4(3,i,km) = max(0., a4(3,i,km))
2099  enddo
2100  elseif ( iv .eq. -1 ) then
2101  do i=i1,i2
2102  if ( a4(3,i,km)*a4(1,i,km) <= 0. ) a4(3,i,km) = 0.
2103  enddo
2104  endif
2105 
2106  do k=km-1,km
2107  do i=i1,i2
2108  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2109  enddo
2110  if(k==(km-1)) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 2)
2111  if(k== km ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1)
2112  enddo
2113 
2114  end subroutine scalar_profile
2115 
2116  subroutine cs_profile(qs, a4, delp, km, i1, i2, iv, kord)
2117 ! Optimized vertical profile reconstruction:
2118 ! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL
2119  integer, intent(in):: i1, i2
2120  integer, intent(in):: km
2121  integer, intent(in):: iv
2124  integer, intent(in):: kord
2125  real, intent(in) :: qs(i1:i2)
2126  real, intent(in) :: delp(i1:i2,km)
2127  real, intent(inout):: a4(4,i1:i2,km)
2128 !-----------------------------------------------------------------------
2129  logical, dimension(i1:i2,km):: extm, ext5, ext6
2130  real gam(i1:i2,km)
2131  real q(i1:i2,km+1)
2132  real d4(i1:i2)
2133  real bet, a_bot, grat
2134  real pmp_1, lac_1, pmp_2, lac_2, x0, x1
2135  integer i, k, im
2136 
2137  if ( iv .eq. -2 ) then
2138  do i=i1,i2
2139  gam(i,2) = 0.5
2140  q(i,1) = 1.5*a4(1,i,1)
2141  enddo
2142  do k=2,km-1
2143  do i=i1, i2
2144  grat = delp(i,k-1) / delp(i,k)
2145  bet = 2. + grat + grat - gam(i,k)
2146  q(i,k) = (3.*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
2147  gam(i,k+1) = grat / bet
2148  enddo
2149  enddo
2150  do i=i1,i2
2151  grat = delp(i,km-1) / delp(i,km)
2152  q(i,km) = (3.*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / &
2153  (2. + grat + grat - gam(i,km))
2154  q(i,km+1) = qs(i)
2155  enddo
2156  do k=km-1,1,-1
2157  do i=i1,i2
2158  q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
2159  enddo
2160  enddo
2161  else
2162  do i=i1,i2
2163  grat = delp(i,2) / delp(i,1) ! grid ratio
2164  bet = grat*(grat+0.5)
2165  q(i,1) = ( (grat+grat)*(grat+1.)*a4(1,i,1) + a4(1,i,2) ) / bet
2166  gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
2167  enddo
2168 
2169  do k=2,km
2170  do i=i1,i2
2171  d4(i) = delp(i,k-1) / delp(i,k)
2172  bet = 2. + d4(i) + d4(i) - gam(i,k-1)
2173  q(i,k) = ( 3.*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
2174  gam(i,k) = d4(i) / bet
2175  enddo
2176  enddo
2177 
2178  do i=i1,i2
2179  a_bot = 1. + d4(i)*(d4(i)+1.5)
2180  q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
2181  / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
2182  enddo
2183 
2184  do k=km,1,-1
2185  do i=i1,i2
2186  q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
2187  enddo
2188  enddo
2189  endif
2190 !----- Perfectly linear scheme --------------------------------
2191  if ( abs(kord) > 16 ) then
2192  do k=1,km
2193  do i=i1,i2
2194  a4(2,i,k) = q(i,k )
2195  a4(3,i,k) = q(i,k+1)
2196  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2197  enddo
2198  enddo
2199  return
2200  endif
2201 !----- Perfectly linear scheme --------------------------------
2202 
2203 !------------------
2204 ! Apply constraints
2205 !------------------
2206  im = i2 - i1 + 1
2207 
2208 ! Apply *large-scale* constraints
2209  do i=i1,i2
2210  q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) )
2211  q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) )
2212  enddo
2213 
2214  do k=2,km
2215  do i=i1,i2
2216  gam(i,k) = a4(1,i,k) - a4(1,i,k-1)
2217  enddo
2218  enddo
2219 
2220 ! Interior:
2221  do k=3,km-1
2222  do i=i1,i2
2223  if ( gam(i,k-1)*gam(i,k+1)>0. ) then
2224 ! Apply large-scale constraint to ALL fields if not local max/min
2225  q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) )
2226  q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) )
2227  else
2228  if ( gam(i,k-1) > 0. ) then
2229 ! There exists a local max
2230  q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k)))
2231  else
2232 ! There exists a local min
2233  q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k)))
2234  if ( iv==0 ) q(i,k) = max(0., q(i,k))
2235  endif
2236  endif
2237  enddo
2238  enddo
2239 
2240 ! Bottom:
2241  do i=i1,i2
2242  q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) )
2243  q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) )
2244  enddo
2245 
2246  do k=1,km
2247  do i=i1,i2
2248  a4(2,i,k) = q(i,k )
2249  a4(3,i,k) = q(i,k+1)
2250  enddo
2251  enddo
2252 
2253  do k=1,km
2254  if ( k==1 .or. k==km ) then
2255  do i=i1,i2
2256  extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0.
2257  enddo
2258  else
2259  do i=i1,i2
2260  extm(i,k) = gam(i,k)*gam(i,k+1) < 0.
2261  enddo
2262  endif
2263  if ( abs(kord) > 9 ) then
2264  do i=i1,i2
2265  x0 = 2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))
2266  x1 = abs(a4(2,i,k)-a4(3,i,k))
2267  a4(4,i,k) = 3.*x0
2268  ext5(i,k) = abs(x0) > x1
2269  ext6(i,k) = abs(a4(4,i,k)) > x1
2270  enddo
2271  endif
2272  enddo
2273 
2274 !---------------------------
2275 ! Apply subgrid constraints:
2276 !---------------------------
2277 ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
2278 ! Top 2 and bottom 2 layers always use monotonic mapping
2279 
2280  if ( iv==0 ) then
2281  do i=i1,i2
2282  a4(2,i,1) = max(0., a4(2,i,1))
2283  enddo
2284  elseif ( iv==-1 ) then
2285  do i=i1,i2
2286  if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0.
2287  enddo
2288  elseif ( iv==2 ) then
2289  do i=i1,i2
2290  a4(2,i,1) = a4(1,i,1)
2291  a4(3,i,1) = a4(1,i,1)
2292  a4(4,i,1) = 0.
2293  enddo
2294  endif
2295 
2296  if ( iv/=2 ) then
2297  do i=i1,i2
2298  a4(4,i,1) = 3.*(2.*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1)))
2299  enddo
2300  call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1)
2301  endif
2302 
2303 ! k=2
2304  do i=i1,i2
2305  a4(4,i,2) = 3.*(2.*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2)))
2306  enddo
2307  call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2)
2308 
2309 !-------------------------------------
2310 ! Huynh's 2nd constraint for interior:
2311 !-------------------------------------
2312  do k=3,km-2
2313  if ( abs(kord)<9 ) then
2314  do i=i1,i2
2315 ! Left edges
2316  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2317  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2318  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2319  max(a4(1,i,k), pmp_1, lac_1) )
2320 ! Right edges
2321  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2322  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2323  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2324  max(a4(1,i,k), pmp_2, lac_2) )
2325 
2326  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2327  enddo
2328 
2329  elseif ( abs(kord)==9 ) then
2330  do i=i1,i2
2331  if ( extm(i,k) .and. extm(i,k-1) ) then ! c90_mp122
2332 ! grid-scale 2-delta-z wave detected
2333  a4(2,i,k) = a4(1,i,k)
2334  a4(3,i,k) = a4(1,i,k)
2335  a4(4,i,k) = 0.
2336  else if ( extm(i,k) .and. extm(i,k+1) ) then ! c90_mp122
2337 ! grid-scale 2-delta-z wave detected
2338  a4(2,i,k) = a4(1,i,k)
2339  a4(3,i,k) = a4(1,i,k)
2340  a4(4,i,k) = 0.
2341  else
2342  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
2343 ! Check within the smooth region if subgrid profile is non-monotonic
2344  if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
2345  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2346  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2347  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2348  max(a4(1,i,k), pmp_1, lac_1) )
2349  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2350  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2351  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2352  max(a4(1,i,k), pmp_2, lac_2) )
2353  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
2354  endif
2355  endif
2356  enddo
2357  elseif ( abs(kord)==10 ) then
2358  do i=i1,i2
2359  if( ext5(i,k) ) then
2360  if( ext5(i,k-1) .or. ext5(i,k+1) ) then
2361  a4(2,i,k) = a4(1,i,k)
2362  a4(3,i,k) = a4(1,i,k)
2363  elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
2364  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2365  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2366  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2367  max(a4(1,i,k), pmp_1, lac_1) )
2368  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2369  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2370  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2371  max(a4(1,i,k), pmp_2, lac_2) )
2372  endif
2373  elseif( ext6(i,k) ) then
2374  if( ext5(i,k-1) .or. ext5(i,k+1) ) then
2375  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2376  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2377  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2378  max(a4(1,i,k), pmp_1, lac_1) )
2379  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2380  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2381  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2382  max(a4(1,i,k), pmp_2, lac_2) )
2383  endif
2384  endif
2385  enddo
2386  do i=i1,i2
2387  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2388  enddo
2389  elseif ( abs(kord)==12 ) then
2390  do i=i1,i2
2391  if( extm(i,k) ) then
2392 ! grid-scale 2-delta-z wave detected
2393  a4(2,i,k) = a4(1,i,k)
2394  a4(3,i,k) = a4(1,i,k)
2395  a4(4,i,k) = 0.
2396  else ! not a local extremum
2397  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
2398 ! Check within the smooth region if subgrid profile is non-monotonic
2399  if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
2400  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2401  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2402  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2403  max(a4(1,i,k), pmp_1, lac_1) )
2404  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2405  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2406  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2407  max(a4(1,i,k), pmp_2, lac_2) )
2408  a4(4,i,k) = 6.*a4(1,i,k) - 3.*(a4(2,i,k)+a4(3,i,k))
2409  endif
2410  endif
2411  enddo
2412  elseif ( abs(kord)==13 ) then
2413  do i=i1,i2
2414  if( ext6(i,k) ) then
2415  if ( ext6(i,k-1) .and. ext6(i,k+1) ) then
2416 ! grid-scale 2-delta-z wave detected
2417  a4(2,i,k) = a4(1,i,k)
2418  a4(3,i,k) = a4(1,i,k)
2419  endif
2420  endif
2421  enddo
2422  do i=i1,i2
2423  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2424  enddo
2425  elseif ( abs(kord)==14 ) then
2426 
2427  do i=i1,i2
2428  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2429  enddo
2430 
2431  elseif ( abs(kord)==15 ) then ! revised kord=9 scehem
2432  do i=i1,i2
2433  if ( ext5(i,k) ) then ! c90_mp122
2434  if ( ext5(i,k-1) .or. ext5(i,k+1) ) then ! c90_mp122
2435 ! grid-scale 2-delta-z wave detected
2436  a4(2,i,k) = a4(1,i,k)
2437  a4(3,i,k) = a4(1,i,k)
2438  endif
2439  elseif( ext6(i,k) ) then
2440 ! Check within the smooth region if subgrid profile is non-monotonic
2441  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2442  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2443  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2444  max(a4(1,i,k), pmp_1, lac_1) )
2445  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2446  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2447  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2448  max(a4(1,i,k), pmp_2, lac_2) )
2449  endif
2450  enddo
2451  do i=i1,i2
2452  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2453  enddo
2454  elseif ( abs(kord)==16 ) then
2455  do i=i1,i2
2456  if( ext5(i,k) ) then
2457  if ( ext5(i,k-1) .or. ext5(i,k+1) ) then
2458  a4(2,i,k) = a4(1,i,k)
2459  a4(3,i,k) = a4(1,i,k)
2460  elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
2461  ! Left edges
2462  pmp_1 = a4(1,i,k) - 2.*gam(i,k+1)
2463  lac_1 = pmp_1 + 1.5*gam(i,k+2)
2464  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
2465  max(a4(1,i,k), pmp_1, lac_1) )
2466  ! Right edges
2467  pmp_2 = a4(1,i,k) + 2.*gam(i,k)
2468  lac_2 = pmp_2 - 1.5*gam(i,k-1)
2469  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
2470  max(a4(1,i,k), pmp_2, lac_2) )
2471  endif
2472  endif
2473  enddo
2474  do i=i1,i2
2475  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2476  enddo
2477  else ! kord = 11
2478  do i=i1,i2
2479  if ( ext5(i,k) .and. (ext5(i,k-1) .or. ext5(i,k+1)) ) then
2480 ! Noisy region:
2481  a4(2,i,k) = a4(1,i,k)
2482  a4(3,i,k) = a4(1,i,k)
2483  a4(4,i,k) = 0.
2484  else
2485  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2486  endif
2487  enddo
2488  endif
2489 
2490 ! Additional constraint to ensure positivity
2491  if ( iv==0 ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 0)
2492 
2493  enddo ! k-loop
2494 
2495 !----------------------------------
2496 ! Bottom layer subgrid constraints:
2497 !----------------------------------
2498  if ( iv==0 ) then
2499  do i=i1,i2
2500  a4(3,i,km) = max(0., a4(3,i,km))
2501  enddo
2502  elseif ( iv .eq. -1 ) then
2503  do i=i1,i2
2504  if ( a4(3,i,km)*a4(1,i,km) <= 0. ) a4(3,i,km) = 0.
2505  enddo
2506  endif
2507 
2508  do k=km-1,km
2509  do i=i1,i2
2510  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2511  enddo
2512  if(k==(km-1)) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 2)
2513  if(k== km ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1)
2514  enddo
2515 
2516  end subroutine cs_profile
2517 
2518 
2519  subroutine cs_limiters(im, extm, a4, iv)
2520  integer, intent(in) :: im
2521  integer, intent(in) :: iv
2522  logical, intent(in) :: extm(im)
2523  real , intent(inout) :: a4(4,im)
2524 ! LOCAL VARIABLES:
2525  real da1, da2, a6da
2526  integer i
2527 
2528  if ( iv==0 ) then
2529 ! Positive definite constraint
2530  do i=1,im
2531  if( a4(1,i)<=0.) then
2532  a4(2,i) = a4(1,i)
2533  a4(3,i) = a4(1,i)
2534  a4(4,i) = 0.
2535  else
2536  if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
2537  if( (a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12) < 0. ) then
2538 ! local minimum is negative
2539  if( a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i) ) then
2540  a4(3,i) = a4(1,i)
2541  a4(2,i) = a4(1,i)
2542  a4(4,i) = 0.
2543  elseif( a4(3,i) > a4(2,i) ) then
2544  a4(4,i) = 3.*(a4(2,i)-a4(1,i))
2545  a4(3,i) = a4(2,i) - a4(4,i)
2546  else
2547  a4(4,i) = 3.*(a4(3,i)-a4(1,i))
2548  a4(2,i) = a4(3,i) - a4(4,i)
2549  endif
2550  endif
2551  endif
2552  endif
2553  enddo
2554  elseif ( iv==1 ) then
2555  do i=1,im
2556  if( (a4(1,i)-a4(2,i))*(a4(1,i)-a4(3,i))>=0. ) then
2557  a4(2,i) = a4(1,i)
2558  a4(3,i) = a4(1,i)
2559  a4(4,i) = 0.
2560  else
2561  da1 = a4(3,i) - a4(2,i)
2562  da2 = da1**2
2563  a6da = a4(4,i)*da1
2564  if(a6da < -da2) then
2565  a4(4,i) = 3.*(a4(2,i)-a4(1,i))
2566  a4(3,i) = a4(2,i) - a4(4,i)
2567  elseif(a6da > da2) then
2568  a4(4,i) = 3.*(a4(3,i)-a4(1,i))
2569  a4(2,i) = a4(3,i) - a4(4,i)
2570  endif
2571  endif
2572  enddo
2573  else
2574 ! Standard PPM constraint
2575  do i=1,im
2576  if( extm(i) ) then
2577  a4(2,i) = a4(1,i)
2578  a4(3,i) = a4(1,i)
2579  a4(4,i) = 0.
2580  else
2581  da1 = a4(3,i) - a4(2,i)
2582  da2 = da1**2
2583  a6da = a4(4,i)*da1
2584  if(a6da < -da2) then
2585  a4(4,i) = 3.*(a4(2,i)-a4(1,i))
2586  a4(3,i) = a4(2,i) - a4(4,i)
2587  elseif(a6da > da2) then
2588  a4(4,i) = 3.*(a4(3,i)-a4(1,i))
2589  a4(2,i) = a4(3,i) - a4(4,i)
2590  endif
2591  endif
2592  enddo
2593  endif
2594  end subroutine cs_limiters
2595 
2596 
2597 
2598  subroutine ppm_profile(a4, delp, km, i1, i2, iv, kord)
2599 
2600 ! INPUT PARAMETERS:
2601  integer, intent(in):: iv
2602  integer, intent(in):: i1
2603  integer, intent(in):: i2
2604  integer, intent(in):: km
2605  integer, intent(in):: kord
2606  !
2607  real , intent(in):: delp(i1:i2,km)
2608 
2609 ! !INPUT/OUTPUT PARAMETERS:
2610  real , intent(inout):: a4(4,i1:i2,km)
2611 
2612 ! DESCRIPTION:
2613 !
2614 ! Perform the piecewise parabolic reconstruction
2615 !
2616 ! !REVISION HISTORY:
2617 ! S.-J. Lin revised at GFDL 2007
2618 !-----------------------------------------------------------------------
2619 ! local arrays:
2620  real dc(i1:i2,km)
2621  real h2(i1:i2,km)
2622  real delq(i1:i2,km)
2623  real df2(i1:i2,km)
2624  real d4(i1:i2,km)
2625 
2626 ! local scalars:
2627  integer i, k, km1, lmt, it
2628  real fac
2629  real a1, a2, c1, c2, c3, d1, d2
2630  real qm, dq, lac, qmp, pmp
2631 
2632  km1 = km - 1
2633  it = i2 - i1 + 1
2634 
2635  do k=2,km
2636  do i=i1,i2
2637  delq(i,k-1) = a4(1,i,k) - a4(1,i,k-1)
2638  d4(i,k ) = delp(i,k-1) + delp(i,k)
2639  enddo
2640  enddo
2641 
2642  do k=2,km1
2643  do i=i1,i2
2644  c1 = (delp(i,k-1)+0.5*delp(i,k))/d4(i,k+1)
2645  c2 = (delp(i,k+1)+0.5*delp(i,k))/d4(i,k)
2646  df2(i,k) = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / &
2647  (d4(i,k)+delp(i,k+1))
2648  dc(i,k) = sign( min(abs(df2(i,k)), &
2649  max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))-a4(1,i,k), &
2650  a4(1,i,k)-min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))), df2(i,k) )
2651  enddo
2652  enddo
2653 
2654 !-----------------------------------------------------------
2655 ! 4th order interpolation of the provisional cell edge value
2656 !-----------------------------------------------------------
2657 
2658  do k=3,km1
2659  do i=i1,i2
2660  c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k)
2661  a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1))
2662  a2 = d4(i,k+1) / (d4(i,k) + delp(i,k))
2663  a4(2,i,k) = a4(1,i,k-1) + c1 + 2./(d4(i,k-1)+d4(i,k+1)) * &
2664  ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - &
2665  delp(i,k-1)*a1*dc(i,k ) )
2666  enddo
2667  enddo
2668 
2669 ! if(km>8 .and. kord>4) call steepz(i1, i2, km, a4, df2, dc, delq, delp, d4)
2670 
2671 ! Area preserving cubic with 2nd deriv. = 0 at the boundaries
2672 ! Top
2673  do i=i1,i2
2674  d1 = delp(i,1)
2675  d2 = delp(i,2)
2676  qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2)
2677  dq = 2.*(a4(1,i,2)-a4(1,i,1)) / (d1+d2)
2678  c1 = 4.*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) )
2679  c3 = dq - 0.5*c1*(d2*(5.*d1+d2)-3.*d1*d1)
2680  a4(2,i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1)
2681 ! Top edge:
2682 !-------------------------------------------------------
2683  a4(2,i,1) = d1*(2.*c1*d1**2-c3) + a4(2,i,2)
2684 !-------------------------------------------------------
2685 ! a4(2,i,1) = (12./7.)*a4(1,i,1)-(13./14.)*a4(1,i,2)+(3./14.)*a4(1,i,3)
2686 !-------------------------------------------------------
2687 ! No over- and undershoot condition
2688  a4(2,i,2) = max( a4(2,i,2), min(a4(1,i,1), a4(1,i,2)) )
2689  a4(2,i,2) = min( a4(2,i,2), max(a4(1,i,1), a4(1,i,2)) )
2690  dc(i,1) = 0.5*(a4(2,i,2) - a4(1,i,1))
2691  enddo
2692 
2693 ! Enforce monotonicity within the top layer
2694 
2695  if( iv==0 ) then
2696  do i=i1,i2
2697  a4(2,i,1) = max(0., a4(2,i,1))
2698  a4(2,i,2) = max(0., a4(2,i,2))
2699  enddo
2700  elseif( iv==-1 ) then
2701  do i=i1,i2
2702  if ( a4(2,i,1)*a4(1,i,1) <= 0. ) a4(2,i,1) = 0.
2703  enddo
2704  elseif( abs(iv)==2 ) then
2705  do i=i1,i2
2706  a4(2,i,1) = a4(1,i,1)
2707  a4(3,i,1) = a4(1,i,1)
2708  enddo
2709  endif
2710 
2711 ! Bottom
2712 ! Area preserving cubic with 2nd deriv. = 0 at the surface
2713  do i=i1,i2
2714  d1 = delp(i,km)
2715  d2 = delp(i,km1)
2716  qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2)
2717  dq = 2.*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2)
2718  c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
2719  c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1*d1)
2720  a4(2,i,km) = qm - c1*d1*d2*(d2+3.*d1)
2721 ! Bottom edge:
2722 !-----------------------------------------------------
2723  a4(3,i,km) = d1*(8.*c1*d1**2-c3) + a4(2,i,km)
2724 ! dc(i,km) = 0.5*(a4(3,i,km) - a4(1,i,km))
2725 !-----------------------------------------------------
2726 ! a4(3,i,km) = (12./7.)*a4(1,i,km)-(13./14.)*a4(1,i,km-1)+(3./14.)*a4(1,i,km-2)
2727 ! No over- and under-shoot condition
2728  a4(2,i,km) = max( a4(2,i,km), min(a4(1,i,km), a4(1,i,km1)) )
2729  a4(2,i,km) = min( a4(2,i,km), max(a4(1,i,km), a4(1,i,km1)) )
2730  dc(i,km) = 0.5*(a4(1,i,km) - a4(2,i,km))
2731  enddo
2732 
2733 
2734 ! Enforce constraint on the "slope" at the surface
2735 
2736 #ifdef BOT_MONO
2737  do i=i1,i2
2738  a4(4,i,km) = 0
2739  if( a4(3,i,km) * a4(1,i,km) <= 0. ) a4(3,i,km) = 0.
2740  d1 = a4(1,i,km) - a4(2,i,km)
2741  d2 = a4(3,i,km) - a4(1,i,km)
2742  if ( d1*d2 < 0. ) then
2743  a4(2,i,km) = a4(1,i,km)
2744  a4(3,i,km) = a4(1,i,km)
2745  else
2746  dq = sign(min(abs(d1),abs(d2),0.5*abs(delq(i,km-1))), d1)
2747  a4(2,i,km) = a4(1,i,km) - dq
2748  a4(3,i,km) = a4(1,i,km) + dq
2749  endif
2750  enddo
2751 #else
2752  if( iv==0 ) then
2753  do i=i1,i2
2754  a4(2,i,km) = max(0.,a4(2,i,km))
2755  a4(3,i,km) = max(0.,a4(3,i,km))
2756  enddo
2757  elseif( iv<0 ) then
2758  do i=i1,i2
2759  if( a4(1,i,km)*a4(3,i,km) <= 0. ) a4(3,i,km) = 0.
2760  enddo
2761  endif
2762 #endif
2763 
2764  do k=1,km1
2765  do i=i1,i2
2766  a4(3,i,k) = a4(2,i,k+1)
2767  enddo
2768  enddo
2769 
2770 !-----------------------------------------------------------
2771 ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
2772 !-----------------------------------------------------------
2773 ! Top 2 and bottom 2 layers always use monotonic mapping
2774  do k=1,2
2775  do i=i1,i2
2776  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2777  enddo
2778  call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
2779  enddo
2780 
2781  if(kord >= 7) then
2782 !-----------------------
2783 ! Huynh's 2nd constraint
2784 !-----------------------
2785  do k=2,km1
2786  do i=i1,i2
2787 ! Method#1
2788 ! h2(i,k) = delq(i,k) - delq(i,k-1)
2789 ! Method#2 - better
2790  h2(i,k) = 2.*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1)) &
2791  / ( delp(i,k)+0.5*(delp(i,k-1)+delp(i,k+1)) ) &
2792  * delp(i,k)**2
2793 ! Method#3
2794 !!! h2(i,k) = dc(i,k+1) - dc(i,k-1)
2795  enddo
2796  enddo
2797 
2798  fac = 1.5 ! original quasi-monotone
2799 
2800  do k=3,km-2
2801  do i=i1,i2
2802 ! Right edges
2803 ! qmp = a4(1,i,k) + 2.0*delq(i,k-1)
2804 ! lac = a4(1,i,k) + fac*h2(i,k-1) + 0.5*delq(i,k-1)
2805 !
2806  pmp = 2.*dc(i,k)
2807  qmp = a4(1,i,k) + pmp
2808  lac = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k)
2809  a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), qmp, lac)), &
2810  max(a4(1,i,k), qmp, lac) )
2811 ! Left edges
2812 ! qmp = a4(1,i,k) - 2.0*delq(i,k)
2813 ! lac = a4(1,i,k) + fac*h2(i,k+1) - 0.5*delq(i,k)
2814 !
2815  qmp = a4(1,i,k) - pmp
2816  lac = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k)
2817  a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), qmp, lac)), &
2818  max(a4(1,i,k), qmp, lac))
2819 !-------------
2820 ! Recompute A6
2821 !-------------
2822  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2823  enddo
2824 ! Additional constraint to ensure positivity when kord=7
2825  if (iv == 0 .and. kord >= 6 ) &
2826  call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 2)
2827  enddo
2828 
2829  else
2830 
2831  lmt = kord - 3
2832  lmt = max(0, lmt)
2833  if (iv == 0) lmt = min(2, lmt)
2834 
2835  do k=3,km-2
2836  if( kord /= 4) then
2837  do i=i1,i2
2838  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2839  enddo
2840  endif
2841  if(kord/=6) call ppm_limiters(dc(i1,k), a4(1,i1,k), it, lmt)
2842  enddo
2843  endif
2844 
2845  do k=km1,km
2846  do i=i1,i2
2847  a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
2848  enddo
2849  call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
2850  enddo
2851 
2852  end subroutine ppm_profile
2853 
2854 
2855  subroutine ppm_limiters(dm, a4, itot, lmt)
2856 
2857 ! INPUT PARAMETERS:
2858  real , intent(in):: dm(*)
2859  integer, intent(in) :: itot
2860  integer, intent(in) :: lmt
2863 ! INPUT/OUTPUT PARAMETERS:
2864  real , intent(inout) :: a4(4,*)
2865 ! LOCAL VARIABLES:
2866  real qmp
2867  real da1, da2, a6da
2868  real fmin
2869  integer i
2870 
2871 ! Developer: S.-J. Lin
2872 
2873  if ( lmt == 3 ) return
2874 
2875  if(lmt == 0) then
2876 ! Standard PPM constraint
2877  do i=1,itot
2878  if(dm(i) == 0.) then
2879  a4(2,i) = a4(1,i)
2880  a4(3,i) = a4(1,i)
2881  a4(4,i) = 0.
2882  else
2883  da1 = a4(3,i) - a4(2,i)
2884  da2 = da1**2
2885  a6da = a4(4,i)*da1
2886  if(a6da < -da2) then
2887  a4(4,i) = 3.*(a4(2,i)-a4(1,i))
2888  a4(3,i) = a4(2,i) - a4(4,i)
2889  elseif(a6da > da2) then
2890  a4(4,i) = 3.*(a4(3,i)-a4(1,i))
2891  a4(2,i) = a4(3,i) - a4(4,i)
2892  endif
2893  endif
2894  enddo
2895 
2896  elseif (lmt == 1) then
2897 
2898 ! Improved full monotonicity constraint (Lin 2004)
2899 ! Note: no need to provide first guess of A6 <-- a4(4,i)
2900  do i=1, itot
2901  qmp = 2.*dm(i)
2902  a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp)
2903  a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp)
2904  a4(4,i) = 3.*( 2.*a4(1,i) - (a4(2,i)+a4(3,i)) )
2905  enddo
2906 
2907  elseif (lmt == 2) then
2908 
2909 ! Positive definite constraint
2910  do i=1,itot
2911  if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
2912  fmin = a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12
2913  if( fmin < 0. ) then
2914  if(a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i)) then
2915  a4(3,i) = a4(1,i)
2916  a4(2,i) = a4(1,i)
2917  a4(4,i) = 0.
2918  elseif(a4(3,i) > a4(2,i)) then
2919  a4(4,i) = 3.*(a4(2,i)-a4(1,i))
2920  a4(3,i) = a4(2,i) - a4(4,i)
2921  else
2922  a4(4,i) = 3.*(a4(3,i)-a4(1,i))
2923  a4(2,i) = a4(3,i) - a4(4,i)
2924  endif
2925  endif
2926  endif
2927  enddo
2928 
2929  endif
2930 
2931  end subroutine ppm_limiters
2932 
2933 
2934 
2935  subroutine steepz(i1, i2, km, a4, df2, dm, dq, dp, d4)
2936  integer, intent(in) :: km, i1, i2
2937  real , intent(in) :: dp(i1:i2,km)
2938  real , intent(in) :: dq(i1:i2,km)
2939  real , intent(in) :: d4(i1:i2,km)
2940  real , intent(in) :: df2(i1:i2,km)
2941  real , intent(in) :: dm(i1:i2,km)
2942 ! INPUT/OUTPUT PARAMETERS:
2943  real , intent(inout) :: a4(4,i1:i2,km)
2944 ! LOCAL VARIABLES:
2945  integer i, k
2946  real alfa(i1:i2,km)
2947  real f(i1:i2,km)
2948  real rat(i1:i2,km)
2949  real dg2
2950 
2951 ! Compute ratio of dq/dp
2952  do k=2,km
2953  do i=i1,i2
2954  rat(i,k) = dq(i,k-1) / d4(i,k)
2955  enddo
2956  enddo
2957 
2958 ! Compute F
2959  do k=2,km-1
2960  do i=i1,i2
2961  f(i,k) = (rat(i,k+1) - rat(i,k)) &
2962  / ( dp(i,k-1)+dp(i,k)+dp(i,k+1) )
2963  enddo
2964  enddo
2965 
2966  do k=3,km-2
2967  do i=i1,i2
2968  if(f(i,k+1)*f(i,k-1)<0. .and. df2(i,k)/=0.) then
2969  dg2 = (f(i,k+1)-f(i,k-1))*((dp(i,k+1)-dp(i,k-1))**2 &
2970  + d4(i,k)*d4(i,k+1) )
2971  alfa(i,k) = max(0., min(0.5, -0.1875*dg2/df2(i,k)))
2972  else
2973  alfa(i,k) = 0.
2974  endif
2975  enddo
2976  enddo
2977 
2978  do k=4,km-2
2979  do i=i1,i2
2980  a4(2,i,k) = (1.-alfa(i,k-1)-alfa(i,k)) * a4(2,i,k) + &
2981  alfa(i,k-1)*(a4(1,i,k)-dm(i,k)) + &
2982  alfa(i,k)*(a4(1,i,k-1)+dm(i,k-1))
2983  enddo
2984  enddo
2985 
2986  end subroutine steepz
2987 
2991  subroutine rst_remap(km, kn, is,ie,js,je, isd,ied,jsd,jed, nq, ntp, &
2992  delp_r, u_r, v_r, w_r, delz_r, pt_r, q_r, qdiag_r, &
2993  delp, u, v, w, delz, pt, q, qdiag, &
2994  ak_r, bk_r, ptop, ak, bk, hydrostatic, make_nh, &
2995  domain, square_domain)
2996 !------------------------------------
2997 ! Assuming hybrid sigma-P coordinate:
2998 !------------------------------------
2999 ! INPUT PARAMETERS:
3000  integer, intent(in):: km
3001  integer, intent(in):: kn
3002  integer, intent(in):: nq, ntp
3003  integer, intent(in):: is,ie,isd,ied
3004  integer, intent(in):: js,je,jsd,jed
3005  logical, intent(in):: hydrostatic, make_nh, square_domain
3006  real, intent(IN) :: ptop
3007  real, intent(in) :: ak_r(km+1)
3008  real, intent(in) :: bk_r(km+1)
3009  real, intent(in) :: ak(kn+1)
3010  real, intent(in) :: bk(kn+1)
3011  real, intent(in):: delp_r(is:ie,js:je,km)
3012  real, intent(in):: u_r(is:ie, js:je+1,km)
3013  real, intent(in):: v_r(is:ie+1,js:je ,km)
3014  real, intent(inout):: pt_r(is:ie,js:je,km)
3015  real, intent(in):: w_r(is:ie,js:je,km)
3016  real, intent(in):: q_r(is:ie,js:je,km,1:ntp)
3017  real, intent(in):: qdiag_r(is:ie,js:je,km,ntp+1:nq)
3018  real, intent(inout)::delz_r(is:ie,js:je,km)
3019  type(domain2d), intent(INOUT) :: domain
3020 ! Output:
3021  real, intent(out):: delp(isd:ied,jsd:jed,kn)
3022  real, intent(out):: u(isd:ied ,jsd:jed+1,kn)
3023  real, intent(out):: v(isd:ied+1,jsd:jed ,kn)
3024  real, intent(out):: w(isd: ,jsd: ,1:)
3025  real, intent(out):: pt(isd:ied ,jsd:jed ,kn)
3026  real, intent(out):: q(isd:ied,jsd:jed,kn,1:ntp)
3027  real, intent(out):: qdiag(isd:ied,jsd:jed,kn,ntp+1:nq)
3028  real, intent(out):: delz(isd:,jsd:,1:)
3029 !-----------------------------------------------------------------------
3030  real r_vir, rgrav
3031  real ps(isd:ied,jsd:jed)
3032  real pe1(is:ie,km+1)
3033  real pe2(is:ie,kn+1)
3034  real pv1(is:ie+1,km+1)
3035  real pv2(is:ie+1,kn+1)
3036 
3037  integer i,j,k , iq
3038  integer, parameter:: kord=4
3039 
3040 #ifdef HYDRO_DELZ_REMAP
3041  if (is_master() .and. .not. hydrostatic) then
3042  print*, ''
3043  print*, ' REMAPPING IC: INITIALIZING DELZ WITH HYDROSTATIC STATE '
3044  print*, ''
3045  endif
3046 #endif
3047 
3048 #ifdef HYDRO_DELZ_EXTRAP
3049  if (is_master() .and. .not. hydrostatic) then
3050  print*, ''
3051  print*, ' REMAPPING IC: INITIALIZING DELZ WITH HYDROSTATIC STATE ABOVE INPUT MODEL TOP '
3052  print*, ''
3053  endif
3054 #endif
3055 
3056 #ifdef ZERO_W_EXTRAP
3057  if (is_master() .and. .not. hydrostatic) then
3058  print*, ''
3059  print*, ' REMAPPING IC: INITIALIZING W TO ZERO ABOVE INPUT MODEL TOP '
3060  print*, ''
3061  endif
3062 #endif
3063 
3064  r_vir = rvgas/rdgas - 1.
3065  rgrav = 1./grav
3066 
3067 !$OMP parallel do default(none) shared(is,ie,js,je,ps,ak_r)
3068  do j=js,je
3069  do i=is,ie
3070  ps(i,j) = ak_r(1)
3071  enddo
3072  enddo
3073 
3074 ! this OpenMP do-loop setup cannot work in it's current form....
3075 !$OMP parallel do default(none) shared(is,ie,js,je,km,ps,delp_r)
3076  do j=js,je
3077  do k=1,km
3078  do i=is,ie
3079  ps(i,j) = ps(i,j) + delp_r(i,j,k)
3080  enddo
3081  enddo
3082  enddo
3083 
3084 ! only one cell is needed
3085  if ( square_domain ) then
3086  call mpp_update_domains(ps, domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.true.)
3087  else
3088  call mpp_update_domains(ps, domain, complete=.true.)
3089  endif
3090 
3091 ! Compute virtual Temp
3092 !$OMP parallel do default(none) shared(is,ie,js,je,km,pt_r,r_vir,q_r)
3093  do k=1,km
3094  do j=js,je
3095  do i=is,ie
3096 #ifdef MULTI_GASES
3097  pt_r(i,j,k) = pt_r(i,j,k) * virq(q_r(i,j,k,:))
3098 #else
3099  pt_r(i,j,k) = pt_r(i,j,k) * (1.+r_vir*q_r(i,j,k,1))
3100 #endif
3101  enddo
3102  enddo
3103  enddo
3104 
3105 !$OMP parallel do default(none) shared(is,ie,js,je,km,ak_r,bk_r,ps,kn,ak,bk,u_r,u,delp, &
3106 !$OMP ntp,nq,hydrostatic,make_nh,w_r,w,delz_r,delp_r,delz, &
3107 !$OMP pt_r,pt,v_r,v,q,q_r,qdiag,qdiag_r) &
3108 !$OMP private(pe1, pe2, pv1, pv2)
3109  do 1000 j=js,je+1
3110 !------
3111 ! map u
3112 !------
3113  do k=1,km+1
3114  do i=is,ie
3115  pe1(i,k) = ak_r(k) + 0.5*bk_r(k)*(ps(i,j-1)+ps(i,j))
3116  enddo
3117  enddo
3118 
3119  do k=1,kn+1
3120  do i=is,ie
3121  pe2(i,k) = ak(k) + 0.5*bk(k)*(ps(i,j-1)+ps(i,j))
3122  enddo
3123  enddo
3124 
3125  call remap_2d(km, pe1, u_r(is:ie,j:j,1:km), &
3126  kn, pe2, u(is:ie,j:j,1:kn), &
3127  is, ie, -1, kord)
3128 
3129  if ( j /= (je+1) ) then
3130 
3131 !---------------
3132 ! Hybrid sigma-p
3133 !---------------
3134  do k=1,km+1
3135  do i=is,ie
3136  pe1(i,k) = ak_r(k) + bk_r(k)*ps(i,j)
3137  enddo
3138  enddo
3139 
3140  do k=1,kn+1
3141  do i=is,ie
3142  pe2(i,k) = ak(k) + bk(k)*ps(i,j)
3143  enddo
3144  enddo
3145 
3146 !-------------
3147 ! Compute delp
3148 !-------------
3149  do k=1,kn
3150  do i=is,ie
3151  delp(i,j,k) = pe2(i,k+1) - pe2(i,k)
3152  enddo
3153  enddo
3154 
3155 !----------------
3156 ! Map constituents
3157 !----------------
3158  if( nq /= 0 ) then
3159  do iq=1,ntp
3160  call remap_2d(km, pe1, q_r(is:ie,j:j,1:km,iq:iq), &
3161  kn, pe2, q(is:ie,j:j,1:kn,iq:iq), &
3162  is, ie, 0, kord)
3163  enddo
3164  do iq=ntp+1,nq
3165  call remap_2d(km, pe1, qdiag_r(is:ie,j:j,1:km,iq:iq), &
3166  kn, pe2, qdiag(is:ie,j:j,1:kn,iq:iq), &
3167  is, ie, 0, kord)
3168  enddo
3169  endif
3170 
3171  if ( .not. hydrostatic .and. .not. make_nh) then
3172 ! Remap vertical wind:
3173  call remap_2d(km, pe1, w_r(is:ie,j:j,1:km), &
3174  kn, pe2, w(is:ie,j:j,1:kn), &
3175  is, ie, -1, kord)
3176 
3177 #ifdef ZERO_W_EXTRAP
3178  do k=1,kn
3179  do i=is,ie
3180  if (pe2(i,k) < pe1(i,1)) then
3181  w(i,j,k) = 0.
3182  endif
3183  enddo
3184  enddo
3185 #endif
3186 
3187 #ifndef HYDRO_DELZ_REMAP
3188 ! Remap delz for hybrid sigma-p coordinate
3189  do k=1,km
3190  do i=is,ie
3191  delz_r(i,j,k) = -delz_r(i,j,k)/delp_r(i,j,k) ! ="specific volume"/grav
3192  enddo
3193  enddo
3194  call remap_2d(km, pe1, delz_r(is:ie,j:j,1:km), &
3195  kn, pe2, delz(is:ie,j:j,1:kn), &
3196  is, ie, 1, kord)
3197  do k=1,kn
3198  do i=is,ie
3199  delz(i,j,k) = -delz(i,j,k)*delp(i,j,k)
3200  enddo
3201  enddo
3202 #endif
3203  endif
3204 
3205 ! Geopotential conserving remap of virtual temperature:
3206  do k=1,km+1
3207  do i=is,ie
3208  pe1(i,k) = log(pe1(i,k))
3209  enddo
3210  enddo
3211  do k=1,kn+1
3212  do i=is,ie
3213  pe2(i,k) = log(pe2(i,k))
3214  enddo
3215  enddo
3216 
3217  call remap_2d(km, pe1, pt_r(is:ie,j:j,1:km), &
3218  kn, pe2, pt(is:ie,j:j,1:kn), &
3219  is, ie, 1, kord)
3220 
3221 #ifdef HYDRO_DELZ_REMAP
3222  !initialize delz from the hydrostatic state
3223  do k=1,kn
3224  do i=is,ie
3225  delz(i,j,k) = (rdgas*rgrav)*pt(i,j,k)*(pe2(i,k)-pe2(i,k+1))
3226  enddo
3227  enddo
3228 #endif
3229 #ifdef HYDRO_DELZ_EXTRAP
3230  !initialize delz from the hydrostatic state
3231  do k=1,kn
3232  do i=is,ie
3233  if (pe2(i,k) < pe1(i,1)) then
3234  delz(i,j,k) = (rdgas*rgrav)*pt(i,j,k)*(pe2(i,k)-pe2(i,k+1))
3235  endif
3236  enddo
3237  enddo
3238 #endif
3239 !------
3240 ! map v
3241 !------
3242  do k=1,km+1
3243  do i=is,ie+1
3244  pv1(i,k) = ak_r(k) + 0.5*bk_r(k)*(ps(i-1,j)+ps(i,j))
3245  enddo
3246  enddo
3247  do k=1,kn+1
3248  do i=is,ie+1
3249  pv2(i,k) = ak(k) + 0.5*bk(k)*(ps(i-1,j)+ps(i,j))
3250  enddo
3251  enddo
3252 
3253  call remap_2d(km, pv1, v_r(is:ie+1,j:j,1:km), &
3254  kn, pv2, v(is:ie+1,j:j,1:kn), &
3255  is, ie+1, -1, kord)
3256 
3257  endif !(j < je+1)
3258 1000 continue
3259 
3260 !$OMP parallel do default(none) shared(is,ie,js,je,kn,pt,r_vir,q)
3261  do k=1,kn
3262  do j=js,je
3263  do i=is,ie
3264 #ifdef MULTI_GASES
3265  pt(i,j,k) = pt(i,j,k) / virq(q(i,j,k,:))
3266 #else
3267  pt(i,j,k) = pt(i,j,k) / (1.+r_vir*q(i,j,k,1))
3268 #endif
3269  enddo
3270  enddo
3271  enddo
3272 
3273  end subroutine rst_remap
3274 
3277  subroutine mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
3278 
3279 ! IV = 0: constituents
3280 ! IV = 1: potential temp
3281 ! IV =-1: winds
3282 
3283 ! Mass flux preserving mapping: q1(im,km) -> q2(im,kn)
3284 
3285 ! pe1: pressure at layer edges (from model top to bottom surface)
3286 ! in the original vertical coordinate
3287 ! pe2: pressure at layer edges (from model top to bottom surface)
3288 ! in the new vertical coordinate
3289 
3290  integer, intent(in):: i1, i2, km, kn, kord, iv
3291  real, intent(in ):: pe1(i1:i2,km+1), pe2(i1:i2,kn+1)
3295 ! Mass flux preserving mapping: q1(im,km) -> q2(im,kn)
3296  real, intent(in ):: q1(i1:i2,km)
3297  real, intent(out):: q2(i1:i2,kn)
3298  real, intent(IN) :: ptop
3299 ! local
3300  real qs(i1:i2)
3301  real dp1(i1:i2,km)
3302  real a4(4,i1:i2,km)
3303  integer i, k, l
3304  integer k0, k1
3305  real pl, pr, tt, delp, qsum, dpsum, esl
3306 
3307  do k=1,km
3308  do i=i1,i2
3309  dp1(i,k) = pe1(i,k+1) - pe1(i,k)
3310  a4(1,i,k) = q1(i,k)
3311  enddo
3312  enddo
3313 
3314  if ( kord >7 ) then
3315  call cs_profile( qs, a4, dp1, km, i1, i2, iv, kord )
3316  else
3317  call ppm_profile( a4, dp1, km, i1, i2, iv, kord )
3318  endif
3319 
3320 !------------------------------------
3321 ! Lowest layer: constant distribution
3322 !------------------------------------
3323 #ifdef NGGPS_SUBMITTED
3324  do i=i1,i2
3325  a4(2,i,km) = q1(i,km)
3326  a4(3,i,km) = q1(i,km)
3327  a4(4,i,km) = 0.
3328  enddo
3329 #endif
3330 
3331  do 5555 i=i1,i2
3332  k0 = 1
3333  do 555 k=1,kn
3334 
3335  if(pe2(i,k) .le. pe1(i,1)) then
3336 ! above old ptop
3337  q2(i,k) = q1(i,1)
3338  elseif(pe2(i,k) .ge. pe1(i,km+1)) then
3339 ! Entire grid below old ps
3340 #ifdef NGGPS_SUBMITTED
3341  q2(i,k) = a4(3,i,km) ! this is not good.
3342 #else
3343  q2(i,k) = q1(i,km)
3344 #endif
3345  else
3346 
3347  do 45 l=k0,km
3348 ! locate the top edge at pe2(i,k)
3349  if( pe2(i,k) .ge. pe1(i,l) .and. &
3350  pe2(i,k) .le. pe1(i,l+1) ) then
3351  k0 = l
3352  pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
3353  if(pe2(i,k+1) .le. pe1(i,l+1)) then
3354 
3355 ! entire new grid is within the original grid
3356  pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
3357  tt = r3*(pr*(pr+pl)+pl**2)
3358  q2(i,k) = a4(2,i,l) + 0.5*(a4(4,i,l)+a4(3,i,l) &
3359  - a4(2,i,l))*(pr+pl) - a4(4,i,l)*tt
3360  goto 555
3361  else
3362 ! Fractional area...
3363  delp = pe1(i,l+1) - pe2(i,k)
3364  tt = r3*(1.+pl*(1.+pl))
3365  qsum = delp*(a4(2,i,l)+0.5*(a4(4,i,l)+ &
3366  a4(3,i,l)-a4(2,i,l))*(1.+pl)-a4(4,i,l)*tt)
3367  dpsum = delp
3368  k1 = l + 1
3369  goto 111
3370  endif
3371  endif
3372 45 continue
3373 
3374 111 continue
3375  do 55 l=k1,km
3376  if( pe2(i,k+1) .gt. pe1(i,l+1) ) then
3377 
3378 ! Whole layer..
3379 
3380  qsum = qsum + dp1(i,l)*q1(i,l)
3381  dpsum = dpsum + dp1(i,l)
3382  else
3383  delp = pe2(i,k+1)-pe1(i,l)
3384  esl = delp / dp1(i,l)
3385  qsum = qsum + delp * (a4(2,i,l)+0.5*esl* &
3386  (a4(3,i,l)-a4(2,i,l)+a4(4,i,l)*(1.-r23*esl)) )
3387  dpsum = dpsum + delp
3388  k0 = l
3389  goto 123
3390  endif
3391 55 continue
3392  delp = pe2(i,k+1) - pe1(i,km+1)
3393  if(delp > 0.) then
3394 ! Extended below old ps
3395 #ifdef NGGPS_SUBMITTED
3396  qsum = qsum + delp * a4(3,i,km) ! not good.
3397 #else
3398  qsum = qsum + delp * q1(i,km)
3399 #endif
3400  dpsum = dpsum + delp
3401  endif
3402 123 q2(i,k) = qsum / dpsum
3403  endif
3404 555 continue
3405 5555 continue
3406 
3407  end subroutine mappm
3408 
3409 
3413  subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
3414  ice_wat, snowwat, graupel, q, qd, cvm, t1)
3415  integer, intent(in):: is, ie, isd,ied, jsd,jed, km, nwat, j, k
3416  integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
3417  real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q
3418  real, intent(out), dimension(is:ie):: cvm, qd
3419  real, intent(in), optional:: t1(is:ie)
3420 !
3421  real, parameter:: t_i0 = 15.
3422  real, dimension(is:ie):: qv, ql, qs
3423  integer:: i
3424 
3425  select case (nwat)
3426 
3427  case(2)
3428  if ( present(t1) ) then ! Special case for GFS physics
3429  do i=is,ie
3430  qd(i) = max(0., q(i,j,k,liq_wat))
3431  if ( t1(i) > tice ) then
3432  qs(i) = 0.
3433  elseif ( t1(i) < tice-t_i0 ) then
3434  qs(i) = qd(i)
3435  else
3436  qs(i) = qd(i)*(tice-t1(i))/t_i0
3437  endif
3438  ql(i) = qd(i) - qs(i)
3439  qv(i) = max(0.,q(i,j,k,sphum))
3440 #ifdef MULTI_GASES
3441  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3442 #else
3443  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3444 #endif
3445  enddo
3446  else
3447  do i=is,ie
3448  qv(i) = max(0.,q(i,j,k,sphum))
3449  qs(i) = max(0.,q(i,j,k,liq_wat))
3450  qd(i) = qs(i)
3451 #ifdef MULTI_GASES
3452  cvm(i) = (1.-qv(i))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap
3453 #else
3454  cvm(i) = (1.-qv(i))*cv_air + qv(i)*cv_vap
3455 #endif
3456  enddo
3457  endif
3458  case (3)
3459  do i=is,ie
3460  qv(i) = q(i,j,k,sphum)
3461  ql(i) = q(i,j,k,liq_wat)
3462  qs(i) = q(i,j,k,ice_wat)
3463  qd(i) = ql(i) + qs(i)
3464  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3465  enddo
3466  case(4) ! K_warm_rain with fake ice
3467  do i=is,ie
3468  qv(i) = q(i,j,k,sphum)
3469  qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3470 #ifdef MULTI_GASES
3471  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + qd(i)*c_liq
3472 #else
3473  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + qd(i)*c_liq
3474 #endif
3475  enddo
3476  case(5)
3477  do i=is,ie
3478  qv(i) = q(i,j,k,sphum)
3479  ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3480  qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat)
3481  qd(i) = ql(i) + qs(i)
3482 #ifdef MULTI_GASES
3483  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3484 #else
3485  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3486 #endif
3487  enddo
3488  case(6)
3489  do i=is,ie
3490  qv(i) = q(i,j,k,sphum)
3491  ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3492  qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat) + q(i,j,k,graupel)
3493  qd(i) = ql(i) + qs(i)
3494 #ifdef MULTI_GASES
3495  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3496 #else
3497  cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
3498 #endif
3499  enddo
3500  case default
3501  call mpp_error (note, 'fv_mapz::moist_cv - using default cv_air')
3502  do i=is,ie
3503  qd(i) = 0.
3504 #ifdef MULTI_GASES
3505  cvm(i) = cv_air*vicvqd(q(i,j,k,1:num_gas))
3506 #else
3507  cvm(i) = cv_air
3508 #endif
3509  enddo
3510  end select
3511 
3512  end subroutine moist_cv
3513 
3516  subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
3517  ice_wat, snowwat, graupel, q, qd, cpm, t1)
3518 
3519  integer, intent(in):: is, ie, isd,ied, jsd,jed, km, nwat, j, k
3520  integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
3521  real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q
3522  real, intent(out), dimension(is:ie):: cpm, qd
3523  real, intent(in), optional:: t1(is:ie)
3524 !
3525  real, parameter:: t_i0 = 15.
3526  real, dimension(is:ie):: qv, ql, qs
3527  integer:: i
3528 
3529  select case (nwat)
3530 
3531  case(2)
3532  if ( present(t1) ) then ! Special case for GFS physics
3533  do i=is,ie
3534  qd(i) = max(0., q(i,j,k,liq_wat))
3535  if ( t1(i) > tice ) then
3536  qs(i) = 0.
3537  elseif ( t1(i) < tice-t_i0 ) then
3538  qs(i) = qd(i)
3539  else
3540  qs(i) = qd(i)*(tice-t1(i))/t_i0
3541  endif
3542  ql(i) = qd(i) - qs(i)
3543  qv(i) = max(0.,q(i,j,k,sphum))
3544 #ifdef MULTI_GASES
3545  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air * vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3546 #else
3547  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3548 #endif
3549  enddo
3550  else
3551  do i=is,ie
3552  qv(i) = max(0.,q(i,j,k,sphum))
3553  qs(i) = max(0.,q(i,j,k,liq_wat))
3554  qd(i) = qs(i)
3555 #ifdef MULTI_GASES
3556  cpm(i) = (1.-qv(i))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor
3557 #else
3558  cpm(i) = (1.-qv(i))*cp_air + qv(i)*cp_vapor
3559 #endif
3560  enddo
3561  endif
3562 
3563  case(3)
3564  do i=is,ie
3565  qv(i) = q(i,j,k,sphum)
3566  ql(i) = q(i,j,k,liq_wat)
3567  qs(i) = q(i,j,k,ice_wat)
3568  qd(i) = ql(i) + qs(i)
3569 #ifdef MULTI_GASES
3570  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3571 #else
3572  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3573 #endif
3574  enddo
3575  case(4) ! K_warm_rain scheme with fake ice
3576  do i=is,ie
3577  qv(i) = q(i,j,k,sphum)
3578  qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3579 #ifdef MULTI_GASES
3580  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + qd(i)*c_liq
3581 #else
3582  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + qd(i)*c_liq
3583 #endif
3584  enddo
3585  case(5)
3586  do i=is,ie
3587  qv(i) = q(i,j,k,sphum)
3588  ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3589  qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat)
3590  qd(i) = ql(i) + qs(i)
3591 #ifdef MULTI_GASES
3592  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3593 #else
3594  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3595 #endif
3596  enddo
3597  case(6)
3598  do i=is,ie
3599  qv(i) = q(i,j,k,sphum)
3600  ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
3601  qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat) + q(i,j,k,graupel)
3602  qd(i) = ql(i) + qs(i)
3603 #ifdef MULTI_GASES
3604  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3605 #else
3606  cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
3607 #endif
3608  enddo
3609  case default
3610  call mpp_error (note, 'fv_mapz::moist_cp - using default cp_air')
3611  do i=is,ie
3612  qd(i) = 0.
3613 #ifdef MULTI_GASES
3614  cpm(i) = cp_air*vicpqd(q(i,j,k,:))
3615 #else
3616  cpm(i) = cp_air
3617 #endif
3618  enddo
3619  end select
3620 
3621  end subroutine moist_cp
3622 
3623 end module fv_mapz_mod
fv_mapz_mod::c_liq
real, parameter c_liq
GFS: heat capacity of water at 0C.
Definition: fv_mapz.F90:122
fv_mapz_mod::cp_vap
real, parameter cp_vap
Definition: fv_mapz.F90:124
multi_gases_mod::vicpqd
pure real function, public vicpqd(q)
Definition: multi_gases.F90:239
fv_mapz_mod::mapn_tracer
subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, i1, i2, isd, ied, jsd, jed, q_min, fill)
Definition: fv_mapz.F90:1418
fv_mapz_mod::r12
real, parameter r12
Definition: fv_mapz.F90:117
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
fv_timing_mod::timing_off
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
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_cmp_mod::fv_sat_adj
subroutine, public fv_sat_adj(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te, te0, qv, ql, qi, qr, qs, qg, hs, dpln, delz, pt, dp, q_con, cappa, area, dtdt, out_dt, last_step, do_qa, qa)
The subroutine &#39;fv_sat_adj&#39; performs the fast processes in the GFDL microphysics. ...
Definition: fv_cmp.F90:130
fv_mapz_mod::moist_cv
subroutine, public moist_cv(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cvm, t1)
The subroutine &#39;moist_cv&#39; computes the FV3-consistent moist heat capacity under constant volume...
Definition: fv_mapz.F90:3415
fv_mapz_mod::mappm
subroutine, public mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
The subroutine &#39;mappm&#39; is a general-purpose routine for remapping one set of vertical levels to anoth...
Definition: fv_mapz.F90:3278
multi_gases_mod::num_gas
integer, public num_gas
Definition: multi_gases.F90:44
multi_gases_mod
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
fv_mapz_mod::steepz
subroutine steepz(i1, i2, km, a4, df2, dm, dq, dp, d4)
Definition: fv_mapz.F90:2936
fv_mapz_mod::pkez
subroutine pkez(km, ifirst, ilast, jfirst, jlast, j, pe, pk, akap, peln, pkz, ptop)
Definition: fv_mapz.F90:1098
fv_mapz_mod::cs_profile
subroutine cs_profile(qs, a4, delp, km, i1, i2, iv, kord)
Definition: fv_mapz.F90:2117
fv_grid_utils_mod::g_sum
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function &#39;g_sum&#39; is the fast version of &#39;globalsum&#39;.
Definition: fv_grid_utils.F90:2902
fv_mapz_mod::compute_total_energy
subroutine, public compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, u, v, w, delz, pt, delp, q, qc, pe, peln, hs, rsin2_l, cosa_s_l, r_vir, cp, rg, hlv, te_2d, ua, va, teq, moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
The subroutine &#39;compute_total_energy&#39; performs the FV3-consistent computation of the global total ene...
Definition: fv_mapz.F90:957
fv_mapz_mod::ppm_limiters
subroutine ppm_limiters(dm, a4, itot, lmt)
Definition: fv_mapz.F90:2856
fv_mapz_mod::r23
real, parameter r23
Definition: fv_mapz.F90:117
fv_mapz_mod::c_ice
real, parameter c_ice
heat capacity of ice at -15.C
Definition: fv_mapz.F90:121
multi_gases_mod::virqd
pure real function, public virqd(q)
Definition: multi_gases.F90:219
fv_mapz_mod::remap_2d
subroutine remap_2d(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord)
Definition: fv_mapz.F90:1621
fv_mapz_mod::consv_min
real, parameter consv_min
below which no correction applies
Definition: fv_mapz.F90:115
fv_timing_mod
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
multi_gases_mod::virq
pure real function, public virq(q)
Definition: multi_gases.F90:136
fv_mapz_mod::cv_vap
real, parameter cv_vap
1384.5
Definition: fv_mapz.F90:118
fv_mapz_mod::lagrangian_to_eulerian
subroutine, public lagrangian_to_eulerian(last_step, consv, ps, pe, delp, pkz, pk, mdt, pdt, km, is, ie, js, je, isd, ied, jsd, jed, nq, nwat, sphum, q_con, u, v, w, delz, pt, q, hs, r_vir, cp, akap, cappa, kord_mt, kord_wz, kord_tr, kord_tm, peln, te0_2d, ng, ua, va, omga, te, ws, fill, reproduce_sum, out_dt, dtdt, ptop, ak, bk, pfull, gridstruct, domain, do_sat_adj, hydrostatic, hybrid_z, do_omega, adiabatic, do_adiabatic_init)
The subroutine &#39;Lagrangian_to_Eulerian&#39; remaps deformed Lagrangian layers back to the reference Euler...
Definition: fv_mapz.F90:145
fv_mapz_mod::cs_limiters
subroutine cs_limiters(im, extm, a4, iv)
Definition: fv_mapz.F90:2520
fv_mapz_mod::t_min
real, parameter t_min
below which applies stricter constraint
Definition: fv_mapz.F90:116
fv_mapz_mod::tice
real, parameter tice
Definition: fv_mapz.F90:125
fv_grid_utils_mod::ptop_min
real, parameter, public ptop_min
Definition: fv_grid_utils.F90:109
fv_mapz_mod
The module &#39;fv_mapz&#39; contains the vertical mapping routines .
Definition: fv_mapz.F90:27
fv_fill_mod::fillz
subroutine, public fillz(im, km, nq, q, dp)
The subroutine &#39;fillz&#39; is for mass-conservative filling of nonphysical negative values in the tracers...
Definition: fv_fill.F90:52
fv_mapz_mod::rst_remap
subroutine, public rst_remap(km, kn, is, ie, js, je, isd, ied, jsd, jed, nq, ntp, delp_r, u_r, v_r, w_r, delz_r, pt_r, q_r, qdiag_r, delp, u, v, w, delz, pt, q, qdiag, ak_r, bk_r, ptop, ak, bk, hydrostatic, make_nh, domain, square_domain)
The subroutine &#39;rst_remap&#39; remaps all variables required for a restart.
Definition: fv_mapz.F90:2996
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
fv_mapz_mod::map_scalar
subroutine map_scalar(km, pe1, q1, qs, kn, pe2, q2, i1, i2, j, ibeg, iend, jbeg, jend, iv, kord, q_min)
Definition: fv_mapz.F90:1240
fv_mapz_mod::map1_q2
subroutine map1_q2(km, pe1, q1, kn, pe2, q2, dp2, i1, i2, iv, kord, j, ibeg, iend, jbeg, jend, q_min)
Definition: fv_mapz.F90:1532
fv_mapz_mod::moist_cp
subroutine, public moist_cp(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cpm, t1)
The subroutine &#39;moist_cp&#39; computes the FV3-consistent moist heat capacity under constant pressure...
Definition: fv_mapz.F90:3518
fv_timing_mod::timing_on
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
fv_mapz_mod::remap_z
subroutine remap_z(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord)
Definition: fv_mapz.F90:1156
fv_grid_utils_mod
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
Definition: fv_grid_utils.F90:26
fv_mapz_mod::ppm_profile
subroutine ppm_profile(a4, delp, km, i1, i2, iv, kord)
Definition: fv_mapz.F90:2599
fv_mapz_mod::cv_air
real, parameter cv_air
= rdgas * (7/2-1) = 2.5*rdgas=717.68
Definition: fv_mapz.F90:119
fv_mapz_mod::e_flux
real(kind=4), public e_flux
Definition: fv_mapz.F90:127
fv_cmp_mod::qs_init
subroutine, public qs_init(kmp)
The subroutine &#39;qs_init&#39; initializes lookup tables for the saturation mixing ratio.
Definition: fv_cmp.F90:1029
fv_mapz_mod::scalar_profile
subroutine scalar_profile(qs, a4, delp, km, i1, i2, iv, kord, qmin)
Definition: fv_mapz.F90:1710
fv_cmp_mod
The module &#39;fv_cmp&#39; implements the fast procesesses in the GFDL microphysics >
Definition: fv_cmp.F90:29
fv_mapz_mod::map1_ppm
subroutine map1_ppm(km, pe1, q1, qs, kn, pe2, q2, i1, i2, j, ibeg, iend, jbeg, jend, iv, kord)
Definition: fv_mapz.F90:1330
multi_gases_mod::vicvqd
pure real function, public vicvqd(q)
Definition: multi_gases.F90:280
fv_fill_mod
Definition: fv_fill.F90:21
fv_mapz_mod::r3
real, parameter r3
Definition: fv_mapz.F90:117