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