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