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