FV3 Bundle
fv_eta_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
20 module fv_eta_nlm_mod
21  use constants_mod, only: kappa, grav, cp_air, rdgas
22  use fv_mp_nlm_mod, only: is_master
23  use mpp_mod, only: fatal, mpp_error
24  implicit none
25  private
26  public set_eta, get_eta_level, compute_dz_var, compute_dz_l32, compute_dz_l101, set_hybrid_z, compute_dz, gw_1d, sm1_edge, hybrid_z_dz
27 
28 ! Developer: Shian-Jiann Lin, NOAA/GFDL
29 
30  contains
31 
32 #ifdef USE_VAR_ETA
33  subroutine set_eta(km, ks, ptop, ak, bk)
34 ! This is the easy to use version of the set_eta
35  integer, intent(in):: km ! vertical dimension
36  integer, intent(out):: ks ! number of pure p layers
37  real:: a60(61),b60(61)
38 ! Thfollowing L63 setting is the same as NCEP GFS's L64 except the top
39 ! 3 layers
40  data a60/300.0000, 430.00000, 558.00000, &
41  700.00000, 863.05803, 1051.07995, &
42  1265.75194, 1510.71101, 1790.05098, &
43  2108.36604, 2470.78817, 2883.03811, &
44  3351.46002, 3883.05187, 4485.49315, &
45  5167.14603, 5937.04991, 6804.87379, &
46  7780.84698, 8875.64338, 10100.20534, &
47  11264.35673, 12190.64366, 12905.42546, &
48  13430.87867, 13785.88765, 13986.77987, &
49  14047.96335, 13982.46770, 13802.40331, &
50  13519.33841, 13144.59486, 12689.45608, &
51  12165.28766, 11583.57006, 10955.84778, &
52  10293.60402, 9608.08306, 8910.07678, &
53  8209.70131, 7516.18560, 6837.69250, &
54  6181.19473, 5552.39653, 4955.72632, &
55  4394.37629, 3870.38682, 3384.76586, &
56  2937.63489, 2528.37666, 2155.78385, &
57  1818.20722, 1513.68173, 1240.03585, &
58  994.99144, 776.23591, 581.48797, &
59  408.53400, 255.26520, 119.70243, 0. /
60 
61  data b60/0.00000, 0.00000, 0.00000, &
62  0.00000, 0.00000, 0.00000, &
63  0.00000, 0.00000, 0.00000, &
64  0.00000, 0.00000, 0.00000, &
65  0.00000, 0.00000, 0.00000, &
66  0.00000, 0.00000, 0.00000, &
67  0.00000, 0.00000, 0.00000, &
68  0.00201, 0.00792, 0.01755, &
69  0.03079, 0.04751, 0.06761, &
70  0.09097, 0.11746, 0.14690, &
71  0.17911, 0.21382, 0.25076, &
72  0.28960, 0.32994, 0.37140, &
73  0.41353, 0.45589, 0.49806, &
74  0.53961, 0.58015, 0.61935, &
75  0.65692, 0.69261, 0.72625, &
76  0.75773, 0.78698, 0.81398, &
77  0.83876, 0.86138, 0.88192, &
78  0.90050, 0.91722, 0.93223, &
79  0.94565, 0.95762, 0.96827, &
80  0.97771, 0.98608, 0.99347, 1./
81  real, intent(out):: ak(km+1)
82  real, intent(out):: bk(km+1)
83  real, intent(out):: ptop ! model top (Pa)
84  real pint, stretch_fac
85  integer k
86  real :: s_rate = -1.0 ! dummy value to not use var_les
87 
88  pint = 100.e2
89 
90 !- Notes ---------------------------------
91 ! low-top: ptop = 100. ! ~45 km
92 ! mid-top: ptop = 10. ! ~60 km
93 ! hi -top: ptop = 1. ! ~80 km
94 !-----------------------------------------
95  select case (km)
96 
97 ! Optimal number = 8 * N -1 (for 8 openMP threads)
98 ! 16 * M -1 (for 16 openMP threads)
99 
100 #ifdef HIWPP
101 #ifdef SUPER_K
102  case (20)
103  ptop = 56.e2
104  pint = ptop
105  stretch_fac = 1.03
106  case (24)
107  ptop = 56.e2
108  pint = ptop
109  stretch_fac = 1.03
110  case (30)
111  ptop = 56.e2
112  pint = ptop
113  stretch_fac = 1.03
114  case (40)
115  ptop = 56.e2
116  pint = ptop
117  stretch_fac = 1.03
118  case (50)
119  ptop = 56.e2
120  pint = ptop
121  stretch_fac = 1.03
122  case (60)
123  ptop = 56.e2
124  pint = ptop
125  stretch_fac = 1.03
126  case (80)
127  ptop = 56.e2
128  pint = ptop
129  stretch_fac = 1.03
130 #else
131  case (30) ! For Baroclinic Instability Test
132  ptop = 2.26e2
133  pint = 250.e2
134  stretch_fac = 1.03
135  case (40)
136  ptop = 50.e2 ! For super cell test
137  pint = 300.e2
138  stretch_fac = 1.03
139  case (50) ! Mountain waves?
140  ptop = 30.e2
141  stretch_fac = 1.025
142  case (60) ! For Baroclinic Instability Test
143 #ifdef GFSL60
144  ks = 20
145  ptop = a60(1)
146  pint = a60(ks+1)
147  do k=1,km+1
148  ak(k) = a60(k)
149  bk(k) = b60(k)
150  enddo
151 #else
152  ptop = 3.e2
153 ! pint = 250.E2
154  pint = 300.e2 ! revised for Moist test
155  stretch_fac = 1.03
156 #endif
157 #endif
158  case (64)
159 !!! ptop = 3.e2
160  ptop = 2.0e2
161  pint = 300.e2
162  stretch_fac = 1.03
163 #else
164 ! *Very-low top: for idealized super-cell simulation:
165  case (50)
166  ptop = 50.e2
167  pint = 250.e2
168  stretch_fac = 1.03
169  case (60)
170  ptop = 40.e2
171  pint = 250.e2
172  stretch_fac = 1.03
173  case (90) ! super-duper cell
174  ptop = 40.e2
175  stretch_fac = 1.025
176 #endif
177 ! Low-top:
178  case (31) ! N = 4, M=2
179  ptop = 100.
180  stretch_fac = 1.035
181  case (32) ! N = 4, M=2
182  ptop = 100.
183  stretch_fac = 1.035
184  case (39) ! N = 5
185  ptop = 100.
186  stretch_fac = 1.035
187  case (41)
188  ptop = 100.
189  stretch_fac = 1.035
190  case (47) ! N = 6, M=3
191  ptop = 100.
192  stretch_fac = 1.035
193  case (51)
194  ptop = 100.
195  stretch_fac = 1.03
196  case (52) ! very low top
197  ptop = 30.e2 ! for special DPM RCE experiments
198  stretch_fac = 1.03
199 ! Mid-top:
200  case (55) ! N = 7
201  ptop = 10.
202  stretch_fac = 1.035
203 ! Hi-top:
204  case (63) ! N = 8, M=4
205  ptop = 1.
206  ! c360 or c384
207  stretch_fac = 1.035
208  case (71) ! N = 9
209  ptop = 1.
210  stretch_fac = 1.03
211  case (79) ! N = 10, M=5
212  ptop = 1.
213  stretch_fac = 1.03
214  case (127) ! N = 10, M=5
215  ptop = 1.
216  stretch_fac = 1.03
217  case (151)
218  ptop = 75.e2
219  pint = 500.e2
220  s_rate = 1.01
221  case default
222  ptop = 1.
223  stretch_fac = 1.03
224  end select
225 
226 #ifdef MOUNTAIN_WAVES
227  call mount_waves(km, ak, bk, ptop, ks, pint)
228 #else
229  if (s_rate > 0.) then
230  call var_les(km, ak, bk, ptop, ks, pint, s_rate)
231  else
232  if ( km > 79 ) then
233  call var_hi2(km, ak, bk, ptop, ks, pint, stretch_fac)
234  elseif (km==5 .or. km==10 ) then
235 ! Equivalent Shallow Water: for NGGPS, variable-resolution testing
236  ptop = 500.e2
237  ks = 0
238  do k=1,km+1
239  bk(k) = real(k-1) / real (km)
240  ak(k) = ptop*(1.-bk(k))
241  enddo
242  else
243 #ifndef GFSL60
244  call var_hi(km, ak, bk, ptop, ks, pint, stretch_fac)
245 #endif
246  endif
247 #endif
248  endif
249 
250  ptop = ak(1)
251  pint = ak(ks+1)
252 
253  end subroutine set_eta
254 
255  subroutine mount_waves(km, ak, bk, ptop, ks, pint)
256  integer, intent(in):: km
257  real, intent(out):: ak(km+1), bk(km+1)
258  real, intent(out):: ptop, pint
259  integer, intent(out):: ks
260 ! Local
261  real, parameter:: p00 = 1.e5
262  real, dimension(km+1):: ze, pe1, peln, eta
263  real, dimension(km):: dz, dlnp
264  real ztop, t0, dz0, sum1, tmp1
265  real ep, es, alpha, beta, gama, s_fac
266  integer k, k500
267 
268  pint = 300.e2
269 ! s_fac = 1.05
270 ! dz0 = 500.
271  if ( km <= 60 ) then
272  s_fac = 1.0
273  dz0 = 500.
274  else
275  s_fac = 1.
276  dz0 = 250.
277  endif
278 
279 ! Basic parameters for HIWPP mountain waves
280  t0 = 300.
281 ! ztop = 20.0e3; 500-m resolution in halft of the vertical domain
282 ! ztop = real(km-1)*500.
283 !-----------------------
284 ! Compute temp ptop based on isothermal atm
285 ! ptop = p00*exp(-grav*ztop/(rdgas*t0))
286 
287 ! Lowest half has constant resolution
288  ze(km+1) = 0.
289  do k=km, km-19, -1
290  ze(k) = ze(k+1) + dz0
291  enddo
292 
293 ! Stretching from 10-km and up:
294  do k=km-20, 3, -1
295  dz0 = s_fac * dz0
296  ze(k) = ze(k+1) + dz0
297  enddo
298  ze(2) = ze(3) + sqrt(2.)*dz0
299  ze(1) = ze(2) + 2.0*dz0
300 
301 ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
302 
303 ! Given z --> p
304  do k=1,km
305  dz(k) = ze(k) - ze(k+1)
306  dlnp(k) = grav*dz(k) / (rdgas*t0)
307  enddo
308 
309  pe1(km+1) = p00
310  peln(km+1) = log(p00)
311  do k=km,1,-1
312  peln(k) = peln(k+1) - dlnp(k)
313  pe1(k) = exp(peln(k))
314  enddo
315 
316 ! Comnpute new ptop
317  ptop = pe1(1)
318 
319 ! Pe(k) = ak(k) + bk(k) * PS
320 ! Locate pint and KS
321  ks = 0
322  do k=2,km
323  if ( pint < pe1(k)) then
324  ks = k-1
325  exit
326  endif
327  enddo
328 
329  if ( is_master() ) then
330  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
331  write(*,*) 'Modified ptop =', ptop, ' ztop=', ze(1)/1000.
332  do k=1,km
333  write(*,*) k, 'ze =', ze(k)/1000.
334  enddo
335  endif
336  pint = pe1(ks+1)
337 
338 #ifdef NO_UKMO_HB
339  do k=1,ks+1
340  ak(k) = pe1(k)
341  bk(k) = 0.
342  enddo
343 
344  do k=ks+2,km+1
345  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
346  ak(k) = pe1(k) - bk(k) * pe1(km+1)
347  enddo
348  bk(km+1) = 1.
349  ak(km+1) = 0.
350 #else
351 ! Problematic for non-hydrostatic
352  do k=1,km+1
353  eta(k) = pe1(k) / pe1(km+1)
354  enddo
355  ep = eta(ks+1)
356  es = eta(km)
357 ! es = 1.
358  alpha = (ep**2-2.*ep*es) / (es-ep)**2
359  beta = 2.*ep*es**2 / (es-ep)**2
360  gama = -(ep*es)**2 / (es-ep)**2
361 
362 ! Pure pressure:
363  do k=1,ks+1
364  ak(k) = eta(k)*1.e5
365  bk(k) = 0.
366  enddo
367 
368  do k=ks+2, km
369  ak(k) = alpha*eta(k) + beta + gama/eta(k)
370  ak(k) = ak(k)*1.e5
371  enddo
372  ak(km+1) = 0.
373 
374  do k=ks+2, km
375  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
376  enddo
377  bk(km+1) = 1.
378 #endif
379 
380  if ( is_master() ) then
381  tmp1 = ak(ks+1)
382  do k=ks+1,km
383  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
384  enddo
385  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
386  endif
387 
388  end subroutine mount_waves
389 
390 #else
391  subroutine set_eta(km, ks, ptop, ak, bk)
392 
393  integer, intent(in):: km ! vertical dimension
394  integer, intent(out):: ks ! number of pure p layers
395  real, intent(out):: ak(km+1)
396  real, intent(out):: bk(km+1)
397  real, intent(out):: ptop ! model top (Pa)
398 ! local
399  real a24(25),b24(25) ! GFDL AM2L24
400  real a26(27),b26(27) ! Jablonowski & Williamson 26-level
401  real a32(33),b32(33)
402  real a32w(33),b32w(33)
403  real a47(48),b47(48)
404  real a48(49),b48(49)
405  real a52(53),b52(53)
406  real a54(55),b54(55)
407  real a56(57),b56(57)
408  real a60(61),b60(61)
409  real a63(64),b63(64)
410  real a64(65),b64(65)
411  real a68(69),b68(69) ! cjg: grid with enhanced PBL resolution
412  real a96(97),b96(97) ! cjg: grid with enhanced PBL resolution
413  real a100(101),b100(101)
414  real a104(105),b104(105)
415  real a125(126),b125(126)
416 
417  real:: p0=1000.e2
418  real:: pc=200.e2
419 
420  real pt, pint, lnpe, dlnp
421  real press(km+1), pt1(km)
422  integer k
423 
424 ! Definition: press(i,j,k) = ak(k) + bk(k) * ps(i,j)
425 
426 !-----------------------------------------------
427 ! GFDL AM2-L24: modified by SJL at the model top
428 !-----------------------------------------------
429 ! data a24 / 100.0000, 1050.0000, 3474.7942, 7505.5556, 12787.2428, &
430  data a24 / 100.0000, 903.4465, 3474.7942, 7505.5556, 12787.2428, &
431  19111.3683, 21854.9274, 22884.1866, 22776.3058, 21716.1604, &
432  20073.2963, 18110.5123, 16004.7832, 13877.6253, 11812.5452, &
433  9865.8840, 8073.9726, 6458.0834, 5027.9899, 3784.6085, &
434  2722.0086, 1828.9752, 1090.2396, 487.4595, 0.0000 /
435 
436  data b24 / 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, &
437  0.0000000, 0.0435679, 0.1102275, 0.1922249, 0.2817656, &
438  0.3694997, 0.4532348, 0.5316253, 0.6038733, 0.6695556, &
439  0.7285176, 0.7808017, 0.8265992, 0.8662148, 0.9000406, &
440  0.9285364, 0.9522140, 0.9716252, 0.9873523, 1.0000000 /
441 
442 ! Jablonowski & Williamson 26-level setup
443  data a26 / 219.4067, 489.5209, 988.2418, 1805.2010, 2983.7240, 4462.3340, &
444  6160.5870, 7851.2430, 7731.2710, 7590.1310, 7424.0860, &
445  7228.7440, 6998.9330, 6728.5740, 6410.5090, 6036.3220, &
446  5596.1110, 5078.2250, 4468.9600, 3752.1910, 2908.9490, &
447  2084.739, 1334.443, 708.499, 252.1360, 0.0, 0.0 /
448 
449  data b26 / 0.0, 0.0, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000,&
450  0.0000000, 0.01505309, 0.03276228, 0.05359622, 0.07810627, &
451  0.1069411, 0.1408637, 0.1807720, 0.2277220, 0.2829562, &
452  0.3479364, 0.4243822, 0.5143168, 0.6201202, 0.7235355, &
453  0.8176768, 0.8962153, 0.9534761, 0.9851122, 1.0000000 /
454 
455 
456 ! High-resolution troposphere setup
457 #ifdef OLD_32
458 ! Revised Apr 14, 2004: PINT = 245.027 mb
459  data a32/100.00000, 400.00000, 818.60211, &
460  1378.88653, 2091.79519, 2983.64084, &
461  4121.78960, 5579.22148, 7419.79300, &
462  9704.82578, 12496.33710, 15855.26306, &
463  19839.62499, 24502.73262, 28177.10152, &
464  29525.28447, 29016.34358, 27131.32792, &
465  24406.11225, 21326.04907, 18221.18357, &
466  15275.14642, 12581.67796, 10181.42843, &
467  8081.89816, 6270.86956, 4725.35001, &
468  3417.39199, 2317.75459, 1398.09473, &
469  632.49506, 0.00000, 0.00000 /
470 
471  data b32/0.00000, 0.00000, 0.00000, &
472  0.00000, 0.00000, 0.00000, &
473  0.00000, 0.00000, 0.00000, &
474  0.00000, 0.00000, 0.00000, &
475  0.00000, 0.00000, 0.01711, &
476  0.06479, 0.13730, 0.22693, &
477  0.32416, 0.42058, 0.51105, &
478  0.59325, 0.66628, 0.73011, &
479  0.78516, 0.83217, 0.87197, &
480  0.90546, 0.93349, 0.95685, &
481  0.97624, 0.99223, 1.00000 /
482 #else
483 ! SJL June 26, 2012
484 ! pint= 55.7922
485  data a32/100.00000, 400.00000, 818.60211, &
486  1378.88653, 2091.79519, 2983.64084, &
487  4121.78960, 5579.22148, 6907.19063, &
488  7735.78639, 8197.66476, 8377.95525, &
489  8331.69594, 8094.72213, 7690.85756, &
490  7139.01788, 6464.80251, 5712.35727, &
491  4940.05347, 4198.60465, 3516.63294, &
492  2905.19863, 2366.73733, 1899.19455, &
493  1497.78137, 1156.25252, 867.79199, &
494  625.59324, 423.21322, 254.76613, &
495  115.06646, 0.00000, 0.00000 /
496 
497  data b32/0.00000, 0.00000, 0.00000, &
498  0.00000, 0.00000, 0.00000, &
499  0.00000, 0.00000, 0.00513, &
500  0.01969, 0.04299, 0.07477, &
501  0.11508, 0.16408, 0.22198, &
502  0.28865, 0.36281, 0.44112, &
503  0.51882, 0.59185, 0.65810, &
504  0.71694, 0.76843, 0.81293, &
505  0.85100, 0.88331, 0.91055, &
506  0.93338, 0.95244, 0.96828, &
507  0.98142, 0.99223, 1.00000 /
508 #endif
509 
510 !---------------------
511 ! Wilson's 32L settings:
512 !---------------------
513 ! Top changed to 0.01 mb
514  data a32w/ 1.00, 26.6378, 84.5529, 228.8592, &
515  539.9597, 1131.7087, 2141.8082, 3712.0454, &
516  5963.5317, 8974.1873, 12764.5388, 17294.5911, &
517  20857.7007, 22221.8651, 22892.7202, 22891.1641, &
518  22286.0724, 21176.0846, 19673.0671, 17889.0989, &
519  15927.5060, 13877.6239, 11812.5474, 9865.8830, &
520  8073.9717, 6458.0824, 5027.9893, 3784.6104, &
521  2722.0093, 1828.9741, 1090.2397, 487.4575, &
522  0.0000 /
523 
524  data b32w/ 0.0000, 0.0000, 0.0000, 0.0000, &
525  0.0000, 0.0000, 0.0000, 0.0000, &
526  0.0000, 0.0000, 0.0000, 0.0000, &
527  0.0159, 0.0586, 0.1117, 0.1734, &
528  0.2415, 0.3137, 0.3878, 0.4619, &
529  0.5344, 0.6039, 0.6696, 0.7285, &
530  0.7808, 0.8266, 0.8662, 0.9000, &
531  0.9285, 0.9522, 0.9716, 0.9874, &
532  1.0000 /
533 
534 
535 #ifdef OLD_L47
536 ! QBO setting with ptop = 0.1 mb and p_full=0.17 mb; pint ~ 100 mb
537  data a47/ 10.00000, 24.45365, 48.76776, &
538  85.39458, 133.41983, 191.01402, &
539  257.94919, 336.63306, 431.52741, &
540  548.18995, 692.78825, 872.16512, &
541  1094.18467, 1368.11917, 1704.99489, &
542  2117.91945, 2622.42986, 3236.88281, &
543  3982.89623, 4885.84733, 5975.43260, &
544  7286.29500, 8858.72424, 10739.43477, &
545  12982.41110, 15649.68745, 18811.37629, &
546  22542.71275, 25724.93857, 27314.36781, &
547  27498.59474, 26501.79312, 24605.92991, &
548  22130.51655, 19381.30274, 16601.56419, &
549  13952.53231, 11522.93244, 9350.82303, &
550  7443.47723, 5790.77434, 4373.32696, &
551  3167.47008, 2148.51663, 1293.15510, &
552  581.62429, 0.00000, 0.00000 /
553 
554  data b47/ 0.0000, 0.0000, 0.0000, &
555  0.00000, 0.00000, 0.00000, &
556  0.00000, 0.00000, 0.00000, &
557  0.00000, 0.00000, 0.00000, &
558  0.00000, 0.00000, 0.00000, &
559  0.00000, 0.00000, 0.00000, &
560  0.00000, 0.00000, 0.00000, &
561  0.00000, 0.00000, 0.00000, &
562  0.00000, 0.00000, 0.00000, &
563  0.00000, 0.01188, 0.04650, &
564  0.10170, 0.17401, 0.25832, &
565  0.34850, 0.43872, 0.52448, &
566  0.60307, 0.67328, 0.73492, &
567  0.78834, 0.83418, 0.87320, &
568  0.90622, 0.93399, 0.95723, &
569  0.97650, 0.99223, 1.00000 /
570 #else
571 ! Oct 23, 2012
572 ! QBO setting with ptop = 0.1 mb, pint ~ 60 mb
573  data a47/ 10.00000, 24.45365, 48.76776, &
574  85.39458, 133.41983, 191.01402, &
575  257.94919, 336.63306, 431.52741, &
576  548.18995, 692.78825, 872.16512, &
577  1094.18467, 1368.11917, 1704.99489, &
578  2117.91945, 2622.42986, 3236.88281, &
579  3982.89623, 4885.84733, 5975.43260, &
580  7019.26669, 7796.15848, 8346.60209, &
581  8700.31838, 8878.27554, 8894.27179, &
582  8756.46404, 8469.60171, 8038.92687, &
583  7475.89006, 6803.68067, 6058.68992, &
584  5285.28859, 4526.01565, 3813.00206, &
585  3164.95553, 2589.26318, 2085.96929, &
586  1651.11596, 1278.81205, 962.38875, &
587  695.07046, 470.40784, 282.61654, &
588  126.92745, 0.00000, 0.00000 /
589  data b47/ 0.0000, 0.0000, 0.0000, &
590  0.00000, 0.00000, 0.00000, &
591  0.00000, 0.00000, 0.00000, &
592  0.00000, 0.00000, 0.00000, &
593  0.00000, 0.00000, 0.00000, &
594  0.00000, 0.00000, 0.00000, &
595  0.00000, 0.00000, 0.00000, &
596  0.00267, 0.01063, 0.02393, &
597  0.04282, 0.06771, 0.09917, &
598  0.13786, 0.18444, 0.23925, &
599  0.30193, 0.37100, 0.44379, &
600  0.51695, 0.58727, 0.65236, &
601  0.71094, 0.76262, 0.80757, &
602  0.84626, 0.87930, 0.90731, &
603  0.93094, 0.95077, 0.96733, &
604  0.98105, 0.99223, 1.00000 /
605 #endif
606 
607  data a48/ &
608  1.00000, 2.69722, 5.17136, &
609  8.89455, 14.24790, 22.07157, &
610  33.61283, 50.48096, 74.79993, &
611  109.40055, 158.00460, 225.44108, &
612  317.89560, 443.19350, 611.11558, &
613  833.74392, 1125.83405, 1505.20759, &
614  1993.15829, 2614.86254, 3399.78420, &
615  4382.06240, 5600.87014, 7100.73115, &
616  8931.78242, 11149.97021, 13817.16841, &
617  17001.20930, 20775.81856, 23967.33875, &
618  25527.64563, 25671.22552, 24609.29622, &
619  22640.51220, 20147.13482, 17477.63530, &
620  14859.86462, 12414.92533, 10201.44191, &
621  8241.50255, 6534.43202, 5066.17865, &
622  3815.60705, 2758.60264, 1870.64631, &
623  1128.33931, 510.47983, 0.00000, &
624  0.00000 /
625 
626  data b48/ &
627  0.00000, 0.00000, 0.00000, &
628  0.00000, 0.00000, 0.00000, &
629  0.00000, 0.00000, 0.00000, &
630  0.00000, 0.00000, 0.00000, &
631  0.00000, 0.00000, 0.00000, &
632  0.00000, 0.00000, 0.00000, &
633  0.00000, 0.00000, 0.00000, &
634  0.00000, 0.00000, 0.00000, &
635  0.00000, 0.00000, 0.00000, &
636  0.00000, 0.00000, 0.01253, &
637  0.04887, 0.10724, 0.18455, &
638  0.27461, 0.36914, 0.46103, &
639  0.54623, 0.62305, 0.69099, &
640  0.75016, 0.80110, 0.84453, &
641  0.88127, 0.91217, 0.93803, &
642  0.95958, 0.97747, 0.99223, &
643  1.00000 /
644 
645 ! High PBL resolution with top at 1 mb
646 ! SJL modified May 7, 2013 to ptop ~ 100 mb
647  data a54/100.00000, 254.83931, 729.54278, &
648  1602.41121, 2797.50667, 4100.18977, &
649  5334.87140, 6455.24153, 7511.80944, &
650  8580.26355, 9714.44293, 10938.62253, &
651  12080.36051, 12987.13921, 13692.75084, &
652  14224.92180, 14606.55444, 14856.69953, &
653  14991.32121, 15023.90075, 14965.91493, &
654  14827.21612, 14616.33505, 14340.72252, &
655  14006.94280, 13620.82849, 13187.60470, &
656  12711.98873, 12198.27003, 11650.37451, &
657  11071.91608, 10466.23819, 9836.44706, &
658  9185.43852, 8515.96231, 7831.01080, &
659  7135.14301, 6436.71659, 5749.00215, &
660  5087.67188, 4465.67510, 3889.86419, &
661  3361.63433, 2879.51065, 2441.02496, &
662  2043.41345, 1683.80513, 1359.31122, &
663  1067.09135, 804.40101, 568.62625, &
664  357.32525, 168.33263, 0.00000, &
665  0.00000 /
666 
667  data b54/0.00000, 0.00000, 0.00000, &
668  0.00000, 0.00000, 0.00000, &
669  0.00000, 0.00000, 0.00000, &
670  0.00000, 0.00000, 0.00000, &
671  0.00180, 0.00694, 0.01510, &
672  0.02601, 0.03942, 0.05515, &
673  0.07302, 0.09288, 0.11459, &
674  0.13803, 0.16307, 0.18960, &
675  0.21753, 0.24675, 0.27716, &
676  0.30866, 0.34115, 0.37456, &
677  0.40879, 0.44375, 0.47935, &
678  0.51551, 0.55215, 0.58916, &
679  0.62636, 0.66334, 0.69946, &
680  0.73395, 0.76622, 0.79594, &
681  0.82309, 0.84780, 0.87020, &
682  0.89047, 0.90876, 0.92524, &
683  0.94006, 0.95336, 0.96529, &
684  0.97596, 0.98551, 0.99400, &
685  1.00000 /
686 
687 
688 ! The 56-L setup
689  data a56/ 10.00000, 24.97818, 58.01160, &
690  115.21466, 199.29210, 309.39897, &
691  445.31785, 610.54747, 812.28518, &
692  1059.80882, 1363.07092, 1732.09335, &
693  2176.91502, 2707.68972, 3334.70962, &
694  4068.31964, 4918.76594, 5896.01890, &
695  7009.59166, 8268.36324, 9680.41211, &
696  11252.86491, 12991.76409, 14901.95764, &
697  16987.01313, 19249.15733, 21689.24182, &
698  23845.11055, 25330.63353, 26243.52467, &
699  26663.84998, 26657.94696, 26281.61371, &
700  25583.05256, 24606.03265, 23393.39510, &
701  21990.28845, 20445.82122, 18811.93894, &
702  17139.59660, 15473.90350, 13850.50167, &
703  12294.49060, 10821.62655, 9440.57746, &
704  8155.11214, 6965.72496, 5870.70511, &
705  4866.83822, 3949.90019, 3115.03562, &
706  2357.07879, 1670.87329, 1051.65120, &
707  495.51399, 0.00000, 0.00000 /
708 
709  data b56 /0.00000, 0.00000, 0.00000, &
710  0.00000, 0.00000, 0.00000, &
711  0.00000, 0.00000, 0.00000, &
712  0.00000, 0.00000, 0.00000, &
713  0.00000, 0.00000, 0.00000, &
714  0.00000, 0.00000, 0.00000, &
715  0.00000, 0.00000, 0.00000, &
716  0.00000, 0.00000, 0.00000, &
717  0.00000, 0.00000, 0.00000, &
718  0.00462, 0.01769, 0.03821, &
719  0.06534, 0.09834, 0.13659, &
720  0.17947, 0.22637, 0.27660, &
721  0.32929, 0.38343, 0.43791, &
722  0.49162, 0.54361, 0.59319, &
723  0.63989, 0.68348, 0.72391, &
724  0.76121, 0.79545, 0.82679, &
725  0.85537, 0.88135, 0.90493, &
726  0.92626, 0.94552, 0.96286, &
727  0.97840, 0.99223, 1.00000 /
728 
729  data a60/ 1.7861000000e-01, 1.0805100000e+00, 3.9647100000e+00, &
730  9.7516000000e+00, 1.9816580000e+01, 3.6695950000e+01, &
731  6.2550570000e+01, 9.9199620000e+01, 1.4792505000e+02, &
732  2.0947487000e+02, 2.8422571000e+02, 3.7241721000e+02, &
733  4.7437835000e+02, 5.9070236000e+02, 7.2236063000e+02, &
734  8.7076746000e+02, 1.0378138800e+03, 1.2258877300e+03, &
735  1.4378924600e+03, 1.6772726600e+03, 1.9480506400e+03, &
736  2.2548762700e+03, 2.6030909400e+03, 2.9988059200e+03, &
737  3.4489952300e+03, 3.9616028900e+03, 4.5456641600e+03, &
738  5.2114401700e+03, 5.9705644000e+03, 6.8361981800e+03, &
739  7.8231906000e+03, 8.9482351000e+03, 1.0230010660e+04, &
740  1.1689289750e+04, 1.3348986860e+04, 1.5234111060e+04, &
741  1.7371573230e+04, 1.9789784580e+04, 2.2005564550e+04, &
742  2.3550115120e+04, 2.4468583320e+04, 2.4800548800e+04, &
743  2.4582445070e+04, 2.3849999620e+04, 2.2640519740e+04, &
744  2.0994737150e+04, 1.8957848730e+04, 1.6579413230e+04, &
745  1.4080071030e+04, 1.1753630920e+04, 9.6516996300e+03, &
746  7.7938009300e+03, 6.1769062800e+03, 4.7874276000e+03, &
747  3.6050497500e+03, 2.6059860700e+03, 1.7668328200e+03, &
748  1.0656131200e+03, 4.8226201000e+02, 0.0000000000e+00, &
749  0.0000000000e+00 /
750 
751 
752  data b60/ 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
753  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
754  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
755  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
756  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
757  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
758  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
759  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
760  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
761  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
762  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
763  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
764  0.0000000000e+00, 0.0000000000e+00, 5.0600000000e-03, &
765  2.0080000000e-02, 4.4900000000e-02, 7.9360000000e-02, &
766  1.2326000000e-01, 1.7634000000e-01, 2.3820000000e-01, &
767  3.0827000000e-01, 3.8581000000e-01, 4.6989000000e-01, &
768  5.5393000000e-01, 6.2958000000e-01, 6.9642000000e-01, &
769  7.5458000000e-01, 8.0463000000e-01, 8.4728000000e-01, &
770  8.8335000000e-01, 9.1368000000e-01, 9.3905000000e-01, &
771  9.6020000000e-01, 9.7775000000e-01, 9.9223000000e-01, &
772  1.0000000000e+00 /
773 
774 ! This is activated by USE_GFSL63
775 ! Thfollowing L63 setting is the same as NCEP GFS's L64 except the top
776 ! 3 layers
777  data a63/64.247, 137.790, 221.958, &
778  318.266, 428.434, 554.424, &
779  698.457, 863.05803, 1051.07995, &
780  1265.75194, 1510.71101, 1790.05098, &
781  2108.36604, 2470.78817, 2883.03811, &
782  3351.46002, 3883.05187, 4485.49315, &
783  5167.14603, 5937.04991, 6804.87379, &
784  7780.84698, 8875.64338, 10100.20534, &
785  11264.35673, 12190.64366, 12905.42546, &
786  13430.87867, 13785.88765, 13986.77987, &
787  14047.96335, 13982.46770, 13802.40331, &
788  13519.33841, 13144.59486, 12689.45608, &
789  12165.28766, 11583.57006, 10955.84778, &
790  10293.60402, 9608.08306, 8910.07678, &
791  8209.70131, 7516.18560, 6837.69250, &
792  6181.19473, 5552.39653, 4955.72632, &
793  4394.37629, 3870.38682, 3384.76586, &
794  2937.63489, 2528.37666, 2155.78385, &
795  1818.20722, 1513.68173, 1240.03585, &
796  994.99144, 776.23591, 581.48797, &
797  408.53400, 255.26520, 119.70243, 0. /
798 
799  data b63/0.00000, 0.00000, 0.00000, &
800  0.00000, 0.00000, 0.00000, &
801  0.00000, 0.00000, 0.00000, &
802  0.00000, 0.00000, 0.00000, &
803  0.00000, 0.00000, 0.00000, &
804  0.00000, 0.00000, 0.00000, &
805  0.00000, 0.00000, 0.00000, &
806  0.00000, 0.00000, 0.00000, &
807  0.00201, 0.00792, 0.01755, &
808  0.03079, 0.04751, 0.06761, &
809  0.09097, 0.11746, 0.14690, &
810  0.17911, 0.21382, 0.25076, &
811  0.28960, 0.32994, 0.37140, &
812  0.41353, 0.45589, 0.49806, &
813  0.53961, 0.58015, 0.61935, &
814  0.65692, 0.69261, 0.72625, &
815  0.75773, 0.78698, 0.81398, &
816  0.83876, 0.86138, 0.88192, &
817  0.90050, 0.91722, 0.93223, &
818  0.94565, 0.95762, 0.96827, &
819  0.97771, 0.98608, 0.99347, 1./
820 #ifdef GFSL64
821  data a64/20.00000, 68.00000, 137.79000, &
822  221.95800, 318.26600, 428.43400, &
823  554.42400, 698.45700, 863.05803, &
824  1051.07995, 1265.75194, 1510.71101, &
825  1790.05098, 2108.36604, 2470.78817, &
826  2883.03811, 3351.46002, 3883.05187, &
827  4485.49315, 5167.14603, 5937.04991, &
828  6804.87379, 7780.84698, 8875.64338, &
829  9921.40745, 10760.99844, 11417.88354, &
830  11911.61193, 12258.61668, 12472.89642, &
831  12566.58298, 12550.43517, 12434.26075, &
832  12227.27484, 11938.39468, 11576.46910, &
833  11150.43640, 10669.41063, 10142.69482, &
834  9579.72458, 8989.94947, 8382.67090, &
835  7766.85063, 7150.91171, 6542.55077, &
836  5948.57894, 5374.81094, 4825.99383, &
837  4305.79754, 3816.84622, 3360.78848, &
838  2938.39801, 2549.69756, 2194.08449, &
839  1870.45732, 1577.34218, 1313.00028, &
840  1075.52114, 862.90778, 673.13815, &
841  504.22118, 354.22752, 221.32110, &
842  103.78014, 0./
843  data b64/0.00000, 0.00000, 0.00000, &
844  0.00000, 0.00000, 0.00000, &
845  0.00000, 0.00000, 0.00000, &
846  0.00000, 0.00000, 0.00000, &
847  0.00000, 0.00000, 0.00000, &
848  0.00000, 0.00000, 0.00000, &
849  0.00000, 0.00000, 0.00000, &
850  0.00000, 0.00000, 0.00000, &
851  0.00179, 0.00705, 0.01564, &
852  0.02749, 0.04251, 0.06064, &
853  0.08182, 0.10595, 0.13294, &
854  0.16266, 0.19492, 0.22950, &
855  0.26615, 0.30455, 0.34435, &
856  0.38516, 0.42656, 0.46815, &
857  0.50949, 0.55020, 0.58989, &
858  0.62825, 0.66498, 0.69987, &
859  0.73275, 0.76351, 0.79208, &
860  0.81845, 0.84264, 0.86472, &
861  0.88478, 0.90290, 0.91923, &
862  0.93388, 0.94697, 0.95865, &
863  0.96904, 0.97826, 0.98642, &
864  0.99363, 1./
865 #else
866  data a64/1.00000, 3.90000, 8.70000, &
867  15.42000, 24.00000, 34.50000, &
868  47.00000, 61.50000, 78.60000, &
869  99.13500, 124.12789, 154.63770, &
870  191.69700, 236.49300, 290.38000, &
871  354.91000, 431.82303, 523.09300, &
872  630.92800, 757.79000, 906.45000, &
873  1079.85000, 1281.00000, 1515.00000, &
874  1788.00000, 2105.00000, 2470.00000, &
875  2889.00000, 3362.00000, 3890.00000, &
876  4475.00000, 5120.00000, 5830.00000, &
877  6608.00000, 7461.00000, 8395.00000, &
878  9424.46289, 10574.46880, 11864.80270, &
879  13312.58890, 14937.03710, 16759.70700, &
880  18804.78710, 21099.41210, 23674.03710, &
881  26562.82810, 29804.11720, 32627.31640, &
882  34245.89840, 34722.28910, 34155.19920, &
883  32636.50390, 30241.08200, 27101.44920, &
884  23362.20700, 19317.05270, 15446.17090, &
885  12197.45210, 9496.39941, 7205.66992, &
886  5144.64307, 3240.79346, 1518.62134, &
887  0.00000, 0.00000 /
888 
889  data b64/0.00000, 0.00000, 0.00000, &
890  0.00000, 0.00000, 0.00000, &
891  0.00000, 0.00000, 0.00000, &
892  0.00000, 0.00000, 0.00000, &
893  0.00000, 0.00000, 0.00000, &
894  0.00000, 0.00000, 0.00000, &
895  0.00000, 0.00000, 0.00000, &
896  0.00000, 0.00000, 0.00000, &
897  0.00000, 0.00000, 0.00000, &
898  0.00000, 0.00000, 0.00000, &
899  0.00000, 0.00000, 0.00000, &
900  0.00000, 0.00000, 0.00000, &
901  0.00000, 0.00000, 0.00000, &
902  0.00000, 0.00000, 0.00000, &
903  0.00000, 0.00000, 0.00000, &
904  0.00000, 0.00000, 0.00813, &
905  0.03224, 0.07128, 0.12445, &
906  0.19063, 0.26929, 0.35799, &
907  0.45438, 0.55263, 0.64304, &
908  0.71703, 0.77754, 0.82827, &
909  0.87352, 0.91502, 0.95235, &
910  0.98511, 1.00000 /
911 #endif
912 !-->cjg
913  data a68/1.00000, 2.68881, 5.15524, &
914  8.86683, 14.20349, 22.00278, &
915  33.50807, 50.32362, 74.56680, &
916  109.05958, 157.51214, 224.73844, &
917  316.90481, 441.81219, 609.21090, &
918  831.14537, 1122.32514, 1500.51628, &
919  1986.94617, 2606.71274, 3389.18802, &
920  4368.40473, 5583.41379, 7078.60015, &
921  8903.94455, 11115.21886, 13774.60566, &
922  16936.82070, 20340.47045, 23193.71492, &
923  24870.36141, 25444.59363, 25252.57081, &
924  24544.26211, 23474.29096, 22230.65331, &
925  20918.50731, 19589.96280, 18296.26682, &
926  17038.02866, 15866.85655, 14763.18943, &
927  13736.83624, 12794.11850, 11930.72442, &
928  11137.17217, 10404.78946, 9720.03954, &
929  9075.54055, 8466.72650, 7887.12346, &
930  7333.90490, 6805.43028, 6297.33773, &
931  5805.78227, 5327.94995, 4859.88765, &
932  4398.63854, 3942.81761, 3491.08449, &
933  3043.04531, 2598.71608, 2157.94527, &
934  1720.87444, 1287.52805, 858.02944, &
935  432.71276, 8.10905, 0.00000 /
936 
937  data b68/0.00000, 0.00000, 0.00000, &
938  0.00000, 0.00000, 0.00000, &
939  0.00000, 0.00000, 0.00000, &
940  0.00000, 0.00000, 0.00000, &
941  0.00000, 0.00000, 0.00000, &
942  0.00000, 0.00000, 0.00000, &
943  0.00000, 0.00000, 0.00000, &
944  0.00000, 0.00000, 0.00000, &
945  0.00000, 0.00000, 0.00000, &
946  0.00000, 0.00283, 0.01590, &
947  0.04412, 0.08487, 0.13284, &
948  0.18470, 0.23828, 0.29120, &
949  0.34211, 0.39029, 0.43518, &
950  0.47677, 0.51536, 0.55091, &
951  0.58331, 0.61263, 0.63917, &
952  0.66333, 0.68552, 0.70617, &
953  0.72555, 0.74383, 0.76117, &
954  0.77765, 0.79335, 0.80838, &
955  0.82287, 0.83693, 0.85069, &
956  0.86423, 0.87760, 0.89082, &
957  0.90392, 0.91689, 0.92973, &
958  0.94244, 0.95502, 0.96747, &
959  0.97979, 0.99200, 1.00000 /
960 
961  data a96/1.00000, 2.35408, 4.51347, &
962  7.76300, 12.43530, 19.26365, &
963  29.33665, 44.05883, 65.28397, &
964  95.48274, 137.90344, 196.76073, &
965  277.45330, 386.81095, 533.37018, &
966  727.67600, 982.60677, 1313.71685, &
967  1739.59104, 2282.20281, 2967.26766, &
968  3824.58158, 4888.33404, 6197.38450, &
969  7795.49158, 9731.48414, 11969.71024, &
970  14502.88894, 17304.52434, 20134.76139, &
971  22536.63814, 24252.54459, 25230.65591, &
972  25585.72044, 25539.91412, 25178.87141, &
973  24644.84493, 23978.98781, 23245.49366, &
974  22492.11600, 21709.93990, 20949.64473, &
975  20225.94258, 19513.31158, 18829.32485, &
976  18192.62250, 17589.39396, 17003.45386, &
977  16439.01774, 15903.91204, 15396.39758, &
978  14908.02140, 14430.65897, 13967.88643, &
979  13524.16667, 13098.30227, 12687.56457, &
980  12287.08757, 11894.41553, 11511.54106, &
981  11139.22483, 10776.01912, 10419.75711, &
982  10067.11881, 9716.63489, 9369.61967, &
983  9026.69066, 8687.29884, 8350.04978, &
984  8013.20925, 7677.12187, 7343.12994, &
985  7011.62844, 6681.98102, 6353.09764, &
986  6025.10535, 5699.10089, 5375.54503, &
987  5053.63074, 4732.62740, 4413.38037, &
988  4096.62775, 3781.79777, 3468.45371, &
989  3157.19882, 2848.25306, 2541.19150, &
990  2236.21942, 1933.50628, 1632.83741, &
991  1334.35954, 1038.16655, 744.22318, &
992  452.71094, 194.91899, 0.00000, &
993  0.00000 /
994 
995  data b96/0.00000, 0.00000, 0.00000, &
996  0.00000, 0.00000, 0.00000, &
997  0.00000, 0.00000, 0.00000, &
998  0.00000, 0.00000, 0.00000, &
999  0.00000, 0.00000, 0.00000, &
1000  0.00000, 0.00000, 0.00000, &
1001  0.00000, 0.00000, 0.00000, &
1002  0.00000, 0.00000, 0.00000, &
1003  0.00000, 0.00000, 0.00000, &
1004  0.00000, 0.00000, 0.00193, &
1005  0.00974, 0.02538, 0.04876, &
1006  0.07817, 0.11081, 0.14514, &
1007  0.18007, 0.21486, 0.24866, &
1008  0.28088, 0.31158, 0.34030, &
1009  0.36701, 0.39210, 0.41554, &
1010  0.43733, 0.45774, 0.47707, &
1011  0.49540, 0.51275, 0.52922, &
1012  0.54495, 0.56007, 0.57459, &
1013  0.58850, 0.60186, 0.61471, &
1014  0.62715, 0.63922, 0.65095, &
1015  0.66235, 0.67348, 0.68438, &
1016  0.69510, 0.70570, 0.71616, &
1017  0.72651, 0.73675, 0.74691, &
1018  0.75700, 0.76704, 0.77701, &
1019  0.78690, 0.79672, 0.80649, &
1020  0.81620, 0.82585, 0.83542, &
1021  0.84492, 0.85437, 0.86375, &
1022  0.87305, 0.88229, 0.89146, &
1023  0.90056, 0.90958, 0.91854, &
1024  0.92742, 0.93623, 0.94497, &
1025  0.95364, 0.96223, 0.97074, &
1026  0.97918, 0.98723, 0.99460, &
1027  1.00000 /
1028 !<--cjg
1029 !
1030 ! Ultra high troposphere resolution
1031  data a100/100.00000, 300.00000, 800.00000, &
1032  1762.35235, 3106.43596, 4225.71874, &
1033  4946.40525, 5388.77387, 5708.35540, &
1034  5993.33124, 6277.73673, 6571.49996, &
1035  6877.05339, 7195.14327, 7526.24920, &
1036  7870.82981, 8229.35361, 8602.30193, &
1037  8990.16936, 9393.46399, 9812.70768, &
1038  10248.43625, 10701.19980, 11171.56286, &
1039  11660.10476, 12167.41975, 12694.11735, &
1040  13240.82253, 13808.17600, 14396.83442, &
1041  15007.47066, 15640.77407, 16297.45067, &
1042  16978.22343, 17683.83253, 18415.03554, &
1043  19172.60771, 19957.34218, 20770.05022, &
1044  21559.14829, 22274.03147, 22916.87519, &
1045  23489.70456, 23994.40187, 24432.71365, &
1046  24806.25734, 25116.52754, 25364.90190, &
1047  25552.64670, 25680.92203, 25750.78675, &
1048  25763.20311, 25719.04113, 25619.08274, &
1049  25464.02630, 25254.49482, 24991.06137, &
1050  24674.32737, 24305.11235, 23884.79781, &
1051  23415.77059, 22901.76510, 22347.84738, &
1052  21759.93950, 21144.07284, 20505.73136, &
1053  19849.54271, 19179.31412, 18498.23400, &
1054  17809.06809, 17114.28232, 16416.10343, &
1055  15716.54833, 15017.44246, 14320.43478, &
1056  13627.01116, 12938.50682, 12256.11762, &
1057  11580.91062, 10913.83385, 10255.72526, &
1058  9607.32122, 8969.26427, 8342.11044, &
1059  7726.33606, 7122.34405, 6530.46991, &
1060  5950.98721, 5384.11279, 4830.01153, &
1061  4288.80090, 3760.55514, 3245.30920, &
1062  2743.06250, 2253.78294, 1777.41285, &
1063  1313.88054, 863.12371, 425.13088, &
1064  0.00000, 0.00000 /
1065 
1066 
1067  data b100/0.00000, 0.00000, 0.00000, &
1068  0.00000, 0.00000, 0.00000, &
1069  0.00000, 0.00000, 0.00000, &
1070  0.00000, 0.00000, 0.00000, &
1071  0.00000, 0.00000, 0.00000, &
1072  0.00000, 0.00000, 0.00000, &
1073  0.00000, 0.00000, 0.00000, &
1074  0.00000, 0.00000, 0.00000, &
1075  0.00000, 0.00000, 0.00000, &
1076  0.00000, 0.00000, 0.00000, &
1077  0.00000, 0.00000, 0.00000, &
1078  0.00000, 0.00000, 0.00000, &
1079  0.00000, 0.00000, 0.00000, &
1080  0.00052, 0.00209, 0.00468, &
1081  0.00828, 0.01288, 0.01849, &
1082  0.02508, 0.03266, 0.04121, &
1083  0.05075, 0.06126, 0.07275, &
1084  0.08521, 0.09866, 0.11308, &
1085  0.12850, 0.14490, 0.16230, &
1086  0.18070, 0.20009, 0.22042, &
1087  0.24164, 0.26362, 0.28622, &
1088  0.30926, 0.33258, 0.35605, &
1089  0.37958, 0.40308, 0.42651, &
1090  0.44981, 0.47296, 0.49591, &
1091  0.51862, 0.54109, 0.56327, &
1092  0.58514, 0.60668, 0.62789, &
1093  0.64872, 0.66919, 0.68927, &
1094  0.70895, 0.72822, 0.74709, &
1095  0.76554, 0.78357, 0.80117, &
1096  0.81835, 0.83511, 0.85145, &
1097  0.86736, 0.88286, 0.89794, &
1098  0.91261, 0.92687, 0.94073, &
1099  0.95419, 0.96726, 0.97994, &
1100  0.99223, 1.00000 /
1101 
1102  data a104/ &
1103  1.8827062944e-01, 7.7977549145e-01, 2.1950593583e+00, &
1104  4.9874566624e+00, 9.8041418997e+00, 1.7019717163e+01, &
1105  2.7216579591e+01, 4.0518628401e+01, 5.6749646818e+01, &
1106  7.5513868331e+01, 9.6315093333e+01, 1.1866706195e+02, &
1107  1.4216835396e+02, 1.6653733709e+02, 1.9161605772e+02, &
1108  2.1735580129e+02, 2.4379516604e+02, 2.7103771847e+02, &
1109  2.9923284173e+02, 3.2856100952e+02, 3.5922338766e+02, &
1110  3.9143507908e+02, 4.2542117983e+02, 4.6141487902e+02, &
1111  4.9965698106e+02, 5.4039638379e+02, 5.8389118154e+02, &
1112  6.3041016829e+02, 6.8023459505e+02, 7.3366009144e+02, &
1113  7.9099869949e+02, 8.5258099392e+02, 9.1875827946e+02, &
1114  9.8990486716e+02, 1.0664204381e+03, 1.1487325074e+03, &
1115  1.2372990044e+03, 1.3326109855e+03, 1.4351954993e+03, &
1116  1.5456186222e+03, 1.6644886848e+03, 1.7924597105e+03, &
1117  1.9302350870e+03, 2.0785714934e+03, 2.2382831070e+03, &
1118  2.4102461133e+03, 2.5954035462e+03, 2.7947704856e+03, &
1119  3.0094396408e+03, 3.2405873512e+03, 3.4894800360e+03, &
1120  3.7574811281e+03, 4.0460585279e+03, 4.3567926151e+03, &
1121  4.6913848588e+03, 5.0516670674e+03, 5.4396113207e+03, &
1122  5.8573406270e+03, 6.3071403487e+03, 6.7914704368e+03, &
1123  7.3129785102e+03, 7.8745138115e+03, 8.4791420557e+03, &
1124  9.1301611750e+03, 9.8311179338e+03, 1.0585825354e+04, &
1125  1.1398380836e+04, 1.2273184781e+04, 1.3214959424e+04, &
1126  1.4228767429e+04, 1.5320029596e+04, 1.6494540743e+04, &
1127  1.7758482452e+04, 1.9118430825e+04, 2.0422798801e+04, &
1128  2.1520147587e+04, 2.2416813461e+04, 2.3118184510e+04, &
1129  2.3628790785e+04, 2.3952411814e+04, 2.4092209011e+04, &
1130  2.4050892106e+04, 2.3830930156e+04, 2.3434818358e+04, &
1131  2.2865410898e+04, 2.2126326004e+04, 2.1222420323e+04, &
1132  2.0160313690e+04, 1.8948920926e+04, 1.7599915822e+04, &
1133  1.6128019809e+04, 1.4550987232e+04, 1.2889169132e+04, &
1134  1.1164595563e+04, 9.4227665517e+03, 7.7259097899e+03, &
1135  6.1538244381e+03, 4.7808126007e+03, 3.5967415552e+03, &
1136  2.5886394104e+03, 1.7415964865e+03, 1.0393721271e+03, &
1137  4.6478852032e+02, 7.0308342481e-13, 0.0000000000e+00 /
1138 
1139 
1140  data b104/ &
1141  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1142  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1143  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1144  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1145  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1146  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1147  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1148  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1149  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1150  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1151  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1152  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1153  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1154  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1155  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1156  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1157  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1158  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1159  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1160  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1161  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1162  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1163  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1164  0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, &
1165  0.0000000000e+00, 0.0000000000e+00, 1.5648447298e-03, &
1166  6.2617046389e-03, 1.4104157933e-02, 2.5118187415e-02, &
1167  3.9340510972e-02, 5.6816335609e-02, 7.7596328431e-02, &
1168  1.0173255472e-01, 1.2927309709e-01, 1.6025505622e-01, &
1169  1.9469566981e-01, 2.3258141217e-01, 2.7385520518e-01, &
1170  3.1840233814e-01, 3.6603639170e-01, 4.1648734767e-01, &
1171  4.6939496013e-01, 5.2431098738e-01, 5.8071350676e-01, &
1172  6.3803478105e-01, 6.9495048840e-01, 7.4963750338e-01, &
1173  7.9975208897e-01, 8.4315257576e-01, 8.8034012292e-01, &
1174  9.1184389721e-01, 9.3821231526e-01, 9.6000677644e-01, &
1175  9.7779792223e-01, 9.9216315122e-01, 1.0000000000e+00 /
1176 
1177 ! IFS-like L125(top 12 levels removed from IFSL137)
1178  data a125/ 64., &
1179  86.895882, 107.415741, 131.425507, 159.279404, 191.338562, &
1180  227.968948, 269.539581, 316.420746, 368.982361, 427.592499, 492.616028, &
1181  564.413452, 643.339905, 729.744141, 823.967834, 926.344910, 1037.201172, &
1182  1156.853638, 1285.610352, 1423.770142, 1571.622925, 1729.448975, 1897.519287, &
1183  2076.095947, 2265.431641, 2465.770508, 2677.348145, 2900.391357, 3135.119385, &
1184  3381.743652, 3640.468262, 3911.490479, 4194.930664, 4490.817383, 4799.149414, &
1185  5119.895020, 5452.990723, 5798.344727, 6156.074219, 6526.946777, 6911.870605, &
1186  7311.869141, 7727.412109, 8159.354004, 8608.525391, 9076.400391, 9562.682617, &
1187  10065.978516, 10584.631836, 11116.662109, 11660.067383, 12211.547852, 12766.873047, &
1188  13324.668945, 13881.331055, 14432.139648, 14975.615234, 15508.256836, 16026.115234, &
1189  16527.322266, 17008.789063, 17467.613281, 17901.621094, 18308.433594, 18685.718750, &
1190  19031.289063, 19343.511719, 19620.042969, 19859.390625, 20059.931641, 20219.664063, &
1191  20337.863281, 20412.308594, 20442.078125, 20425.718750, 20361.816406, 20249.511719, &
1192  20087.085938, 19874.025391, 19608.572266, 19290.226563, 18917.460938, 18489.707031, &
1193  18006.925781, 17471.839844, 16888.687500, 16262.046875, 15596.695313, 14898.453125, &
1194  14173.324219, 13427.769531, 12668.257813, 11901.339844, 11133.304688, 10370.175781, &
1195  9617.515625, 8880.453125, 8163.375000, 7470.343750, 6804.421875, 6168.531250, &
1196  5564.382813, 4993.796875, 4457.375000, 3955.960938, 3489.234375, 3057.265625, &
1197  2659.140625, 2294.242188, 1961.500000, 1659.476563, 1387.546875, 1143.250000, &
1198  926.507813, 734.992188, 568.062500, 424.414063, 302.476563, 202.484375, &
1199  122.101563, 62.781250, 22.835938, 3.757813, 0.000000, 0.000000 /
1200 
1201  data b125/ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1202  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1203  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1204  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1205  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1206  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1207  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1208  0.000000, 0.000007, 0.000024, 0.000059, 0.000112, 0.000199, &
1209  0.000340, 0.000562, 0.000890, 0.001353, 0.001992, 0.002857, &
1210  0.003971, 0.005378, 0.007133, 0.009261, 0.011806, 0.014816, &
1211  0.018318, 0.022355, 0.026964, 0.032176, 0.038026, 0.044548, &
1212  0.051773, 0.059728, 0.068448, 0.077958, 0.088286, 0.099462, &
1213  0.111505, 0.124448, 0.138313, 0.153125, 0.168910, 0.185689, &
1214  0.203491, 0.222333, 0.242244, 0.263242, 0.285354, 0.308598, &
1215  0.332939, 0.358254, 0.384363, 0.411125, 0.438391, 0.466003, &
1216  0.493800, 0.521619, 0.549301, 0.576692, 0.603648, 0.630036, &
1217  0.655736, 0.680643, 0.704669, 0.727739, 0.749797, 0.770798, &
1218  0.790717, 0.809536, 0.827256, 0.843881, 0.859432, 0.873929, &
1219  0.887408, 0.899900, 0.911448, 0.922096, 0.931881, 0.940860, &
1220  0.949064, 0.956550, 0.963352, 0.969513, 0.975078, 0.980072, &
1221  0.984542, 0.988500, 0.991984, 0.995003, 0.997630, 1.000000 /
1222 
1223  select case (km)
1224 
1225  case (24)
1226 
1227  ks = 5
1228  do k=1,km+1
1229  ak(k) = a24(k)
1230  bk(k) = b24(k)
1231  enddo
1232 
1233  case (26)
1234 
1235  ks = 7
1236  do k=1,km+1
1237  ak(k) = a26(k)
1238  bk(k) = b26(k)
1239  enddo
1240 
1241  case (32)
1242 #ifdef OLD_32
1243  ks = 13 ! high-res trop_32 setup
1244 #else
1245  ks = 7
1246 #endif
1247  do k=1,km+1
1248  ak(k) = a32(k)
1249  bk(k) = b32(k)
1250  enddo
1251 
1252  case (47)
1253 ! ks = 27 ! high-res trop-strat
1254  ks = 20 ! Oct 23, 2012
1255  do k=1,km+1
1256  ak(k) = a47(k)
1257  bk(k) = b47(k)
1258  enddo
1259 
1260  case (48)
1261  ks = 28
1262  do k=1,km+1
1263  ak(k) = a48(k)
1264  bk(k) = b48(k)
1265  enddo
1266 
1267  case (52)
1268  ks = 35 ! pint = 223
1269  do k=1,km+1
1270  ak(k) = a52(k)
1271  bk(k) = b52(k)
1272  enddo
1273 
1274  case (54)
1275  ks = 11 ! pint = 109.4
1276  do k=1,km+1
1277  ak(k) = a54(k)
1278  bk(k) = b54(k)
1279  enddo
1280 
1281  case (56)
1282  ks = 26
1283  do k=1,km+1
1284  ak(k) = a56(k)
1285  bk(k) = b56(k)
1286  enddo
1287 
1288  case (60)
1289  ks = 37
1290  do k=1,km+1
1291  ak(k) = a60(k)
1292  bk(k) = b60(k)
1293  enddo
1294 
1295 
1296  case (64)
1297 #ifdef GFSL64
1298  ks = 23
1299 #else
1300  ks = 46
1301 #endif
1302  do k=1,km+1
1303  ak(k) = a64(k)
1304  bk(k) = b64(k)
1305  enddo
1306 !-->cjg
1307  case (68)
1308  ks = 27
1309  do k=1,km+1
1310  ak(k) = a68(k)
1311  bk(k) = b68(k)
1312  enddo
1313 
1314  case (96)
1315  ks = 27
1316  do k=1,km+1
1317  ak(k) = a96(k)
1318  bk(k) = b96(k)
1319  enddo
1320 !<--cjg
1321 
1322  case (100)
1323  ks = 38
1324  do k=1,km+1
1325  ak(k) = a100(k)
1326  bk(k) = b100(k)
1327  enddo
1328 
1329  case (104)
1330  ks = 73
1331  do k=1,km+1
1332  ak(k) = a104(k)
1333  bk(k) = b104(k)
1334  enddo
1335 
1336 #ifndef TEST_GWAVES
1337  case (10)
1338 !--------------------------------------------------
1339 ! Pure sigma-coordinate with uniform spacing in "z"
1340 !--------------------------------------------------
1341 !
1342  pt = 2000. ! model top pressure (pascal)
1343 ! pt = 100. ! 1 mb
1344  press(1) = pt
1345  press(km+1) = p0
1346  dlnp = (log(p0) - log(pt)) / real(km)
1347 
1348  lnpe = log(press(km+1))
1349  do k=km,2,-1
1350  lnpe = lnpe - dlnp
1351  press(k) = exp(lnpe)
1352  enddo
1353 
1354 ! Search KS
1355  ks = 0
1356  do k=1,km
1357  if(press(k) >= pc) then
1358  ks = k-1
1359  goto 123
1360  endif
1361  enddo
1362 123 continue
1363 
1364  if(ks /= 0) then
1365  do k=1,ks
1366  ak(k) = press(k)
1367  bk(k) = 0.
1368  enddo
1369  endif
1370 
1371  pint = press(ks+1)
1372  do k=ks+1,km
1373  ak(k) = pint*(press(km)-press(k))/(press(km)-pint)
1374  bk(k) = (press(k) - ak(k)) / press(km+1)
1375  enddo
1376  ak(km+1) = 0.
1377  bk(km+1) = 1.
1378 
1379 ! do k=2,km
1380 ! bk(k) = real(k-1) / real(km)
1381 ! ak(k) = pt * ( 1. - bk(k) )
1382 ! enddo
1383 #endif
1384 
1385 ! The following 4 selections are better for non-hydrostatic
1386 ! Low top:
1387  case (31)
1388  ptop = 300.
1389  pint = 100.e2
1390  call var_dz(km, ak, bk, ptop, ks, pint, 1.035)
1391 #ifndef TEST_GWAVES
1392  case (41)
1393  ptop = 100.
1394  pint = 100.e2
1395  call var_hi(km, ak, bk, ptop, ks, pint, 1.035)
1396 #endif
1397  case (51)
1398  ptop = 100.
1399  pint = 100.e2
1400  call var_hi(km, ak, bk, ptop, ks, pint, 1.035)
1401 ! Mid-top:
1402  case (55)
1403  ptop = 10.
1404  pint = 100.e2
1405 ! call var_dz(km, ak, bk, ptop, ks, pint, 1.035)
1406  call var_hi(km, ak, bk, ptop, ks, pint, 1.035)
1407 #ifdef USE_GFSL63
1408 ! GFS L64 equivalent setting
1409  case (63)
1410  ks = 23
1411  ptop = a63(1)
1412  pint = a63(ks+1)
1413  do k=1,km+1
1414  ak(k) = a63(k)
1415  bk(k) = b63(k)
1416  enddo
1417 #else
1418  case (63)
1419  ptop = 1. ! high top
1420  pint = 100.e2
1421  call var_hi(km, ak, bk, ptop, ks, pint, 1.035)
1422 #endif
1423 ! NGGPS_GFS
1424  case (91)
1425  pint = 100.e2
1426  ptop = 40.
1427  call var_gfs(km, ak, bk, ptop, ks, pint, 1.029)
1428  case (95)
1429 ! Mid-top settings:
1430  pint = 100.e2
1431  ptop = 20.
1432  call var_gfs(km, ak, bk, ptop, ks, pint, 1.028)
1433  case (127)
1434  ptop = 1.
1435  pint = 75.e2
1436  call var_gfs(km, ak, bk, ptop, ks, pint, 1.028)
1437 ! IFS-like L125
1438  case (125)
1439  ks = 33
1440  ptop = a125(1)
1441  pint = a125(ks+1)
1442  do k=1,km+1
1443  ak(k) = a125(k)
1444  bk(k) = b125(k)
1445  enddo
1446  case default
1447 
1448 #ifdef TEST_GWAVES
1449 !--------------------------------------------------
1450 ! Pure sigma-coordinate with uniform spacing in "z"
1451 !--------------------------------------------------
1452  call gw_1d(km, 1000.e2, ak, bk, ptop, 10.e3, pt1)
1453  ks = 0
1454  pint = ak(1)
1455 #else
1456 
1457 !----------------------------------------------------------------
1458 ! Sigma-coordinate with uniform spacing in sigma and ptop = 1 mb
1459 !----------------------------------------------------------------
1460  pt = 100.
1461 ! One pressure layer
1462  ks = 1
1463 ! pint = pt + 0.5*1.E5/real(km) ! SJL: 20120327
1464  pint = pt + 1.e5/real(km)
1465 
1466  ak(1) = pt
1467  bk(1) = 0.
1468  ak(2) = pint
1469  bk(2) = 0.
1470 
1471  do k=3,km+1
1472  bk(k) = real(k-2) / real(km-1)
1473  ak(k) = pint - bk(k)*pint
1474  enddo
1475  ak(km+1) = 0.
1476  bk(km+1) = 1.
1477 #endif
1478  end select
1479  ptop = ak(1)
1480  pint = ak(ks+1)
1481 
1482  end subroutine set_eta
1483 #endif
1484 
1485  subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate)
1486  implicit none
1487  integer, intent(in):: km
1488  real, intent(in):: ptop
1489  real, intent(in):: s_rate ! between [1. 1.1]
1490  real, intent(out):: ak(km+1), bk(km+1)
1491  real, intent(inout):: pint
1492  integer, intent(out):: ks
1493 ! Local
1494  real, parameter:: p00 = 1.e5
1495  real, dimension(km+1):: ze, pe1, peln, eta
1496  real, dimension(km):: dz, s_fac, dlnp, pm, dp, dk
1497  real ztop, t0, dz0, sum1, tmp1
1498  real ep, es, alpha, beta, gama
1499  real, parameter:: akap = 2./7.
1500 !---- Tunable parameters:
1501  integer:: k_inc = 10 ! # of layers from bottom up to near const dz region
1502  real:: s0 = 0.8 ! lowest layer stretch factor
1503 !-----------------------
1504  real:: s_inc
1505  integer k
1506 
1507  pe1(1) = ptop
1508  peln(1) = log(pe1(1))
1509  pe1(km+1) = p00
1510  peln(km+1) = log(pe1(km+1))
1511 
1512  t0 = 273.
1513  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1514 
1515  s_inc = (1.-s0) / real(k_inc)
1516  s_fac(km) = s0
1517 
1518  do k=km-1, km-k_inc, -1
1519  s_fac(k) = s_fac(k+1) + s_inc
1520  enddo
1521 
1522  s_fac(km-k_inc-1) = 0.5*(s_fac(km-k_inc) + s_rate)
1523 
1524  do k=km-k_inc-2, 5, -1
1525  s_fac(k) = s_rate * s_fac(k+1)
1526  enddo
1527 
1528  s_fac(4) = 0.5*(1.1+s_rate)*s_fac(5)
1529  s_fac(3) = 1.1 *s_fac(4)
1530  s_fac(2) = 1.1 *s_fac(3)
1531  s_fac(1) = 1.1 *s_fac(2)
1532 
1533  sum1 = 0.
1534  do k=1,km
1535  sum1 = sum1 + s_fac(k)
1536  enddo
1537 
1538  dz0 = ztop / sum1
1539 
1540  do k=1,km
1541  dz(k) = s_fac(k) * dz0
1542  enddo
1543 
1544  ze(km+1) = 0.
1545  do k=km,1,-1
1546  ze(k) = ze(k+1) + dz(k)
1547  enddo
1548 
1549 ! Re-scale dz with the stretched ztop
1550  do k=1,km
1551  dz(k) = dz(k) * (ztop/ze(1))
1552  enddo
1553 
1554  do k=km,1,-1
1555  ze(k) = ze(k+1) + dz(k)
1556  enddo
1557 ! ze(1) = ztop
1558 
1559  if ( is_master() ) then
1560  write(*,*) 'var_les: computed model top (m)=', ztop, ' bottom/top dz=', dz(km), dz(1)
1561 ! do k=1,km
1562 ! write(*,*) k, s_fac(k)
1563 ! enddo
1564  endif
1565 
1566  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
1567 
1568 ! Given z --> p
1569  do k=1,km
1570  dz(k) = ze(k) - ze(k+1)
1571  dlnp(k) = grav*dz(k) / (rdgas*t0)
1572  !write(*,*) k, dz(k)
1573  enddo
1574  do k=2,km
1575  peln(k) = peln(k-1) + dlnp(k-1)
1576  pe1(k) = exp(peln(k))
1577  enddo
1578 
1579 
1580 ! Pe(k) = ak(k) + bk(k) * PS
1581 ! Locate pint and KS
1582  ks = 0
1583  do k=2,km
1584  if ( pint < pe1(k)) then
1585  ks = k-1
1586  exit
1587  endif
1588  enddo
1589  if ( is_master() ) then
1590  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1591  endif
1592  pint = pe1(ks+1)
1593 
1594  do k=1,km+1
1595  eta(k) = pe1(k) / pe1(km+1)
1596  enddo
1597 
1598  ep = eta(ks+1)
1599  es = eta(km)
1600 ! es = 1.
1601  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1602  beta = 2.*ep*es**2 / (es-ep)**2
1603  gama = -(ep*es)**2 / (es-ep)**2
1604 
1605 ! Pure pressure:
1606  do k=1,ks+1
1607  ak(k) = eta(k)*1.e5
1608  bk(k) = 0.
1609  enddo
1610 
1611  do k=ks+2, km
1612  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1613  ak(k) = ak(k)*1.e5
1614  enddo
1615  ak(km+1) = 0.
1616 
1617  do k=ks+2, km
1618  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1619  enddo
1620  bk(km+1) = 1.
1621 
1622  if ( is_master() ) then
1623  ! write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1624  ! do k=1,km
1625  ! pm(k) = 0.5*(pe1(k)+pe1(k+1))/100.
1626  ! write(*,*) k, pm(k), dz(k)
1627  ! enddo
1628  tmp1 = ak(ks+1)
1629  do k=ks+1,km
1630  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1631  enddo
1632  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1633  write(*,800) (pm(k), k=km,1,-1)
1634  endif
1635 
1636  do k=1,km
1637  dp(k) = (pe1(k+1) - pe1(k))/100.
1638  dk(k) = pe1(k+1)**akap - pe1(k)**akap
1639  enddo
1640 
1641 800 format(1x,5(1x,f9.4))
1642 
1643  end subroutine var_les
1644 
1645 
1646 
1647  subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate)
1648  integer, intent(in):: km
1649  real, intent(in):: ptop
1650  real, intent(in):: s_rate ! between [1. 1.1]
1651  real, intent(out):: ak(km+1), bk(km+1)
1652  real, intent(inout):: pint
1653  integer, intent(out):: ks
1654 ! Local
1655  real, parameter:: p00 = 1.e5
1656  real, dimension(km+1):: ze, pe1, peln, eta
1657  real, dimension(km):: dz, s_fac, dlnp
1658  real ztop, t0, dz0, sum1, tmp1
1659  real ep, es, alpha, beta, gama
1660 !---- Tunable parameters:
1661  integer:: k_inc = 25 ! # of layers from bottom up to near const dz region
1662  real:: s0 = 0.13 ! lowest layer stretch factor
1663 !-----------------------
1664  real:: s_inc
1665  integer k
1666 
1667  pe1(1) = ptop
1668  peln(1) = log(pe1(1))
1669  pe1(km+1) = p00
1670  peln(km+1) = log(pe1(km+1))
1671 
1672  t0 = 270.
1673  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1674 
1675  s_inc = (1.-s0) / real(k_inc)
1676  s_fac(km) = s0
1677 
1678  do k=km-1, km-k_inc, -1
1679  s_fac(k) = s_fac(k+1) + s_inc
1680  enddo
1681 
1682  do k=km-k_inc-1, 9, -1
1683  s_fac(k) = s_rate * s_fac(k+1)
1684  enddo
1685  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
1686  s_fac(7) = 1.10*s_fac(8)
1687  s_fac(6) = 1.15*s_fac(7)
1688  s_fac(5) = 1.20*s_fac(6)
1689  s_fac(4) = 1.26*s_fac(5)
1690  s_fac(3) = 1.33*s_fac(4)
1691  s_fac(2) = 1.41*s_fac(3)
1692  s_fac(1) = 1.60*s_fac(2)
1693 
1694  sum1 = 0.
1695  do k=1,km
1696  sum1 = sum1 + s_fac(k)
1697  enddo
1698 
1699  dz0 = ztop / sum1
1700 
1701  do k=1,km
1702  dz(k) = s_fac(k) * dz0
1703  enddo
1704 
1705  ze(km+1) = 0.
1706  do k=km,1,-1
1707  ze(k) = ze(k+1) + dz(k)
1708  enddo
1709 
1710 ! Re-scale dz with the stretched ztop
1711  do k=1,km
1712  dz(k) = dz(k) * (ztop/ze(1))
1713  enddo
1714 
1715  do k=km,1,-1
1716  ze(k) = ze(k+1) + dz(k)
1717  enddo
1718 ! ze(1) = ztop
1719 
1720  if ( is_master() ) then
1721  write(*,*) 'var_gfs: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1722 ! do k=1,km
1723 ! write(*,*) k, s_fac(k)
1724 ! enddo
1725  endif
1726 
1727 ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 3)
1728 
1729 ! Given z --> p
1730  do k=1,km
1731  dz(k) = ze(k) - ze(k+1)
1732  dlnp(k) = grav*dz(k) / (rdgas*t0)
1733  enddo
1734  do k=2,km
1735  peln(k) = peln(k-1) + dlnp(k-1)
1736  pe1(k) = exp(peln(k))
1737  enddo
1738 
1739 ! Pe(k) = ak(k) + bk(k) * PS
1740 ! Locate pint and KS
1741  ks = 0
1742  do k=2,km
1743  if ( pint < pe1(k)) then
1744  ks = k-1
1745  exit
1746  endif
1747  enddo
1748  if ( is_master() ) then
1749  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1750  write(*,*) 'ptop =', ptop
1751  endif
1752  pint = pe1(ks+1)
1753 
1754 #ifdef NO_UKMO_HB
1755  do k=1,ks+1
1756  ak(k) = pe1(k)
1757  bk(k) = 0.
1758  enddo
1759 
1760  do k=ks+2,km+1
1761  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1762  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1763  enddo
1764  bk(km+1) = 1.
1765  ak(km+1) = 0.
1766 #else
1767 ! Problematic for non-hydrostatic
1768  do k=1,km+1
1769  eta(k) = pe1(k) / pe1(km+1)
1770  enddo
1771 
1772  ep = eta(ks+1)
1773  es = eta(km)
1774 ! es = 1.
1775  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1776  beta = 2.*ep*es**2 / (es-ep)**2
1777  gama = -(ep*es)**2 / (es-ep)**2
1778 
1779 ! Pure pressure:
1780  do k=1,ks+1
1781  ak(k) = eta(k)*1.e5
1782  bk(k) = 0.
1783  enddo
1784 
1785  do k=ks+2, km
1786  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1787  ak(k) = ak(k)*1.e5
1788  enddo
1789  ak(km+1) = 0.
1790 
1791  do k=ks+2, km
1792  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1793  enddo
1794  bk(km+1) = 1.
1795 #endif
1796 
1797  if ( is_master() ) then
1798  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1799  do k=1,km
1800  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1801  enddo
1802  tmp1 = ak(ks+1)
1803  do k=ks+1,km
1804  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1805  enddo
1806  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1807  endif
1808 
1809  end subroutine var_gfs
1810 
1811  subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate)
1812  integer, intent(in):: km
1813  real, intent(in):: ptop
1814  real, intent(in):: s_rate ! between [1. 1.1]
1815  real, intent(out):: ak(km+1), bk(km+1)
1816  real, intent(inout):: pint
1817  integer, intent(out):: ks
1818 ! Local
1819  real, parameter:: p00 = 1.e5
1820  real, dimension(km+1):: ze, pe1, peln, eta
1821  real, dimension(km):: dz, s_fac, dlnp
1822  real ztop, t0, dz0, sum1, tmp1
1823  real ep, es, alpha, beta, gama
1824 !---- Tunable parameters:
1825  integer:: k_inc = 15 ! # of layers from bottom up to near const dz region
1826  real:: s0 = 0.10 ! lowest layer stretch factor
1827 !-----------------------
1828  real:: s_inc
1829  integer k
1830 
1831  pe1(1) = ptop
1832  peln(1) = log(pe1(1))
1833  pe1(km+1) = p00
1834  peln(km+1) = log(pe1(km+1))
1835 
1836  t0 = 270.
1837  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1838 
1839  s_inc = (1.-s0) / real(k_inc)
1840  s_fac(km) = s0
1841 
1842  do k=km-1, km-k_inc, -1
1843  s_fac(k) = s_fac(k+1) + s_inc
1844  enddo
1845 
1846  s_fac(km-k_inc-1) = 0.5*(s_fac(km-k_inc) + s_rate)
1847 
1848 #ifdef HIWPP
1849  do k=km-k_inc-2, 4, -1
1850  s_fac(k) = s_rate * s_fac(k+1)
1851  enddo
1852  s_fac(3) = 0.5*(1.15+s_rate)*s_fac(4)
1853  s_fac(2) = 1.15 *s_fac(3)
1854  s_fac(1) = 1.3 *s_fac(2)
1855 #else
1856  do k=km-k_inc-2, 9, -1
1857  s_fac(k) = s_rate * s_fac(k+1)
1858  enddo
1859 
1860  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
1861  s_fac(7) = 1.1 *s_fac(8)
1862  s_fac(6) = 1.15*s_fac(7)
1863  s_fac(5) = 1.2 *s_fac(6)
1864  s_fac(4) = 1.3 *s_fac(5)
1865  s_fac(3) = 1.4 *s_fac(4)
1866  s_fac(2) = 1.45 *s_fac(3)
1867  s_fac(1) = 1.5 *s_fac(2)
1868 #endif
1869 
1870  sum1 = 0.
1871  do k=1,km
1872  sum1 = sum1 + s_fac(k)
1873  enddo
1874 
1875  dz0 = ztop / sum1
1876 
1877  do k=1,km
1878  dz(k) = s_fac(k) * dz0
1879  enddo
1880 
1881  ze(km+1) = 0.
1882  do k=km,1,-1
1883  ze(k) = ze(k+1) + dz(k)
1884  enddo
1885 
1886 ! Re-scale dz with the stretched ztop
1887  do k=1,km
1888  dz(k) = dz(k) * (ztop/ze(1))
1889  enddo
1890 
1891  do k=km,1,-1
1892  ze(k) = ze(k+1) + dz(k)
1893  enddo
1894 ! ze(1) = ztop
1895 
1896  if ( is_master() ) then
1897  write(*,*) 'var_hi: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1898 ! do k=1,km
1899 ! write(*,*) k, s_fac(k)
1900 ! enddo
1901  endif
1902 
1903  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
1904 
1905 ! Given z --> p
1906  do k=1,km
1907  dz(k) = ze(k) - ze(k+1)
1908  dlnp(k) = grav*dz(k) / (rdgas*t0)
1909  enddo
1910  do k=2,km
1911  peln(k) = peln(k-1) + dlnp(k-1)
1912  pe1(k) = exp(peln(k))
1913  enddo
1914 
1915 ! Pe(k) = ak(k) + bk(k) * PS
1916 ! Locate pint and KS
1917  ks = 0
1918  do k=2,km
1919  if ( pint < pe1(k)) then
1920  ks = k-1
1921  exit
1922  endif
1923  enddo
1924  if ( is_master() ) then
1925  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1926  write(*,*) 'ptop =', ptop
1927  endif
1928  pint = pe1(ks+1)
1929 
1930 #ifdef NO_UKMO_HB
1931  do k=1,ks+1
1932  ak(k) = pe1(k)
1933  bk(k) = 0.
1934  enddo
1935 
1936  do k=ks+2,km+1
1937  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1938  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1939  enddo
1940  bk(km+1) = 1.
1941  ak(km+1) = 0.
1942 #else
1943 ! Problematic for non-hydrostatic
1944  do k=1,km+1
1945  eta(k) = pe1(k) / pe1(km+1)
1946  enddo
1947 
1948  ep = eta(ks+1)
1949  es = eta(km)
1950 ! es = 1.
1951  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1952  beta = 2.*ep*es**2 / (es-ep)**2
1953  gama = -(ep*es)**2 / (es-ep)**2
1954 
1955 ! Pure pressure:
1956  do k=1,ks+1
1957  ak(k) = eta(k)*1.e5
1958  bk(k) = 0.
1959  enddo
1960 
1961  do k=ks+2, km
1962  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1963  ak(k) = ak(k)*1.e5
1964  enddo
1965  ak(km+1) = 0.
1966 
1967  do k=ks+2, km
1968  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1969  enddo
1970  bk(km+1) = 1.
1971 #endif
1972 
1973  if ( is_master() ) then
1974  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1975  do k=1,km
1976  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1977  enddo
1978  tmp1 = ak(ks+1)
1979  do k=ks+1,km
1980  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1981  enddo
1982  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1983  endif
1984 
1985 
1986  end subroutine var_hi
1987  subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate)
1988  integer, intent(in):: km
1989  real, intent(in):: ptop
1990  real, intent(in):: s_rate ! between [1. 1.1]
1991  real, intent(out):: ak(km+1), bk(km+1)
1992  real, intent(inout):: pint
1993  integer, intent(out):: ks
1994 ! Local
1995  real, parameter:: p00 = 1.e5
1996  real, dimension(km+1):: ze, pe1, peln, eta
1997  real, dimension(km):: dz, s_fac, dlnp
1998  real ztop, t0, dz0, sum1, tmp1
1999  real ep, es, alpha, beta, gama
2000  integer k
2001 
2002  pe1(1) = ptop
2003  peln(1) = log(pe1(1))
2004  pe1(km+1) = p00
2005  peln(km+1) = log(pe1(km+1))
2006 
2007  t0 = 270.
2008  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
2009 
2010  s_fac(km ) = 0.15
2011  s_fac(km-1) = 0.20
2012  s_fac(km-2) = 0.30
2013  s_fac(km-3) = 0.40
2014  s_fac(km-4) = 0.50
2015  s_fac(km-5) = 0.60
2016  s_fac(km-6) = 0.70
2017  s_fac(km-7) = 0.80
2018  s_fac(km-8) = 0.90
2019  s_fac(km-9) = 0.95
2020  s_fac(km-10) = 0.5*(s_fac(km-9) + s_rate)
2021 
2022  do k=km-11, 8, -1
2023  s_fac(k) = s_rate * s_fac(k+1)
2024  enddo
2025 
2026  s_fac(7) = 0.5*(1.1+s_rate)*s_fac(9)
2027  s_fac(6) = 1.05*s_fac(7)
2028  s_fac(5) = 1.1*s_fac(6)
2029  s_fac(4) = 1.15*s_fac(5)
2030  s_fac(3) = 1.2*s_fac(4)
2031  s_fac(2) = 1.3*s_fac(3)
2032  s_fac(1) = 1.4*s_fac(2)
2033 
2034  sum1 = 0.
2035  do k=1,km
2036  sum1 = sum1 + s_fac(k)
2037  enddo
2038 
2039  dz0 = ztop / sum1
2040 
2041  do k=1,km
2042  dz(k) = s_fac(k) * dz0
2043  enddo
2044 
2045  ze(km+1) = 0.
2046  do k=km,1,-1
2047  ze(k) = ze(k+1) + dz(k)
2048  enddo
2049 
2050 ! Re-scale dz with the stretched ztop
2051  do k=1,km
2052  dz(k) = dz(k) * (ztop/ze(1))
2053  enddo
2054 
2055  do k=km,1,-1
2056  ze(k) = ze(k+1) + dz(k)
2057  enddo
2058 ! ze(1) = ztop
2059 
2060  if ( is_master() ) write(*,*) 'var_hi2: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
2061  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
2062 
2063 ! Given z --> p
2064  do k=1,km
2065  dz(k) = ze(k) - ze(k+1)
2066  dlnp(k) = grav*dz(k) / (rdgas*t0)
2067  enddo
2068  do k=2,km
2069  peln(k) = peln(k-1) + dlnp(k-1)
2070  pe1(k) = exp(peln(k))
2071  enddo
2072 
2073 ! Pe(k) = ak(k) + bk(k) * PS
2074 ! Locate pint and KS
2075  ks = 0
2076  do k=2,km
2077  if ( pint < pe1(k)) then
2078  ks = k-1
2079  exit
2080  endif
2081  enddo
2082  if ( is_master() ) then
2083  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
2084  endif
2085  pint = pe1(ks+1)
2086 
2087 #ifdef NO_UKMO_HB
2088  do k=1,ks+1
2089  ak(k) = pe1(k)
2090  bk(k) = 0.
2091  enddo
2092 
2093  do k=ks+2,km+1
2094  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
2095  ak(k) = pe1(k) - bk(k) * pe1(km+1)
2096  enddo
2097  bk(km+1) = 1.
2098  ak(km+1) = 0.
2099 #else
2100 ! Problematic for non-hydrostatic
2101  do k=1,km+1
2102  eta(k) = pe1(k) / pe1(km+1)
2103  enddo
2104 
2105  ep = eta(ks+1)
2106  es = eta(km)
2107 ! es = 1.
2108  alpha = (ep**2-2.*ep*es) / (es-ep)**2
2109  beta = 2.*ep*es**2 / (es-ep)**2
2110  gama = -(ep*es)**2 / (es-ep)**2
2111 
2112 ! Pure pressure:
2113  do k=1,ks+1
2114  ak(k) = eta(k)*1.e5
2115  bk(k) = 0.
2116  enddo
2117 
2118  do k=ks+2, km
2119  ak(k) = alpha*eta(k) + beta + gama/eta(k)
2120  ak(k) = ak(k)*1.e5
2121  enddo
2122  ak(km+1) = 0.
2123 
2124  do k=ks+2, km
2125  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
2126  enddo
2127  bk(km+1) = 1.
2128 #endif
2129 
2130  if ( is_master() ) then
2131  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
2132  do k=1,km
2133  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
2134  enddo
2135  tmp1 = ak(ks+1)
2136  do k=ks+1,km
2137  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
2138  enddo
2139  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
2140  endif
2141 
2142 
2143  end subroutine var_hi2
2144 
2145 
2146  subroutine var_dz(km, ak, bk, ptop, ks, pint, s_rate)
2147  integer, intent(in):: km
2148  real, intent(in):: ptop
2149  real, intent(in):: s_rate ! between [1. 1.1]
2150  real, intent(out):: ak(km+1), bk(km+1)
2151  real, intent(inout):: pint
2152  integer, intent(out):: ks
2153 ! Local
2154  real, parameter:: p00 = 1.e5
2155  real, dimension(km+1):: ze, pe1, peln, eta
2156  real, dimension(km):: dz, s_fac, dlnp
2157  real ztop, t0, dz0, sum1, tmp1
2158  real ep, es, alpha, beta, gama
2159  integer k
2160 
2161  pe1(1) = ptop
2162  peln(1) = log(pe1(1))
2163  pe1(km+1) = p00
2164  peln(km+1) = log(pe1(km+1))
2165 
2166  t0 = 270.
2167  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
2168 
2169  s_fac(km ) = 0.10
2170  s_fac(km-1) = 0.20
2171  s_fac(km-2) = 0.30
2172  s_fac(km-3) = 0.40
2173  s_fac(km-4) = 0.50
2174  s_fac(km-5) = 0.60
2175  s_fac(km-6) = 0.70
2176  s_fac(km-7) = 0.80
2177  s_fac(km-8) = 0.90
2178  s_fac(km-9) = 0.95
2179  s_fac(km-10) = 0.5*(s_fac(km-9) + s_rate)
2180 
2181  do k=km-11, 9, -1
2182  s_fac(k) = min(10.0, s_rate * s_fac(k+1) )
2183  enddo
2184 
2185  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
2186  s_fac(7) = 1.1 *s_fac(8)
2187  s_fac(6) = 1.15*s_fac(7)
2188  s_fac(5) = 1.2 *s_fac(6)
2189  s_fac(4) = 1.3 *s_fac(5)
2190  s_fac(3) = 1.4 *s_fac(4)
2191  s_fac(2) = 1.5 *s_fac(3)
2192  s_fac(1) = 1.6 *s_fac(2)
2193 
2194  sum1 = 0.
2195  do k=1,km
2196  sum1 = sum1 + s_fac(k)
2197  enddo
2198 
2199  dz0 = ztop / sum1
2200 
2201  do k=1,km
2202  dz(k) = s_fac(k) * dz0
2203  enddo
2204 
2205  ze(km+1) = 0.
2206  do k=km,1,-1
2207  ze(k) = ze(k+1) + dz(k)
2208  enddo
2209 
2210 ! Re-scale dz with the stretched ztop
2211  do k=1,km
2212  dz(k) = dz(k) * (ztop/ze(1))
2213  enddo
2214 
2215  do k=km,1,-1
2216  ze(k) = ze(k+1) + dz(k)
2217  enddo
2218 ! ze(1) = ztop
2219 
2220  if ( is_master() ) write(*,*) 'var_dz: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
2221  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
2222 
2223 ! Given z --> p
2224  do k=1,km
2225  dz(k) = ze(k) - ze(k+1)
2226  dlnp(k) = grav*dz(k) / (rdgas*t0)
2227  enddo
2228  do k=2,km
2229  peln(k) = peln(k-1) + dlnp(k-1)
2230  pe1(k) = exp(peln(k))
2231  enddo
2232 
2233 ! Pe(k) = ak(k) + bk(k) * PS
2234 ! Locate pint and KS
2235  ks = 0
2236  do k=2,km
2237  if ( pint < pe1(k)) then
2238  ks = k-1
2239  exit
2240  endif
2241  enddo
2242  if ( is_master() ) then
2243  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
2244  write(*,*) 'ptop =', ptop
2245  endif
2246  pint = pe1(ks+1)
2247 
2248 #ifdef NO_UKMO_HB
2249  do k=1,ks+1
2250  ak(k) = pe1(k)
2251  bk(k) = 0.
2252  enddo
2253 
2254  do k=ks+2,km+1
2255  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
2256  ak(k) = pe1(k) - bk(k) * pe1(km+1)
2257  enddo
2258  bk(km+1) = 1.
2259  ak(km+1) = 0.
2260 #else
2261 ! Problematic for non-hydrostatic
2262  do k=1,km+1
2263  eta(k) = pe1(k) / pe1(km+1)
2264  enddo
2265 
2266  ep = eta(ks+1)
2267  es = eta(km)
2268 ! es = 1.
2269  alpha = (ep**2-2.*ep*es) / (es-ep)**2
2270  beta = 2.*ep*es**2 / (es-ep)**2
2271  gama = -(ep*es)**2 / (es-ep)**2
2272 
2273 ! Pure pressure:
2274  do k=1,ks+1
2275  ak(k) = eta(k)*1.e5
2276  bk(k) = 0.
2277  enddo
2278 
2279  do k=ks+2, km
2280  ak(k) = alpha*eta(k) + beta + gama/eta(k)
2281  ak(k) = ak(k)*1.e5
2282  enddo
2283  ak(km+1) = 0.
2284 
2285  do k=ks+2, km
2286  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
2287  enddo
2288  bk(km+1) = 1.
2289 #endif
2290 
2291  if ( is_master() ) then
2292  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
2293  do k=1,km
2294  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
2295  enddo
2296  tmp1 = ak(ks+1)
2297  do k=ks+1,km
2298  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
2299  enddo
2300  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
2301  endif
2302 
2303 
2304  end subroutine var_dz
2305 
2306 
2307  subroutine var55_dz(km, ak, bk, ptop, ks, pint, s_rate)
2308  integer, intent(in):: km
2309  real, intent(in):: ptop
2310  real, intent(in):: s_rate ! between [1. 1.1]
2311  real, intent(out):: ak(km+1), bk(km+1)
2312  real, intent(inout):: pint
2313  integer, intent(out):: ks
2314 ! Local
2315  real, parameter:: p00 = 1.e5
2316  real, dimension(km+1):: ze, pe1, peln, eta
2317  real, dimension(km):: dz, s_fac, dlnp
2318  real ztop, t0, dz0, sum1, tmp1
2319  real ep, es, alpha, beta, gama
2320  integer k
2321 
2322  pe1(1) = ptop
2323  peln(1) = log(pe1(1))
2324  pe1(km+1) = p00
2325  peln(km+1) = log(pe1(km+1))
2326 
2327  t0 = 270.
2328  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
2329 
2330  s_fac(km ) = 0.10
2331  s_fac(km-1) = 0.15
2332  s_fac(km-2) = 0.20
2333  s_fac(km-3) = 0.25
2334  s_fac(km-4) = 0.30
2335  s_fac(km-5) = 0.35
2336  s_fac(km-6) = 0.40
2337  s_fac(km-7) = 0.45
2338  s_fac(km-8) = 0.50
2339  s_fac(km-9) = 0.55
2340  s_fac(km-10) = 0.60
2341  s_fac(km-11) = 0.70
2342  s_fac(km-12) = 0.85
2343  s_fac(km-13) = 1.00
2344 
2345  do k=km-14, 9, -1
2346  s_fac(k) = min(10.0, s_rate * s_fac(k+1) )
2347  enddo
2348 
2349  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
2350  s_fac(7) = 1.1 *s_fac(8)
2351  s_fac(6) = 1.15*s_fac(7)
2352  s_fac(5) = 1.2 *s_fac(6)
2353  s_fac(4) = 1.3 *s_fac(5)
2354  s_fac(3) = 1.4 *s_fac(4)
2355  s_fac(2) = 1.5 *s_fac(3)
2356  s_fac(1) = 1.6 *s_fac(2)
2357 
2358  sum1 = 0.
2359  do k=1,km
2360  sum1 = sum1 + s_fac(k)
2361  enddo
2362 
2363  dz0 = ztop / sum1
2364 
2365  do k=1,km
2366  dz(k) = s_fac(k) * dz0
2367  enddo
2368 
2369  ze(km+1) = 0.
2370  do k=km,1,-1
2371  ze(k) = ze(k+1) + dz(k)
2372  enddo
2373 
2374 ! Re-scale dz with the stretched ztop
2375  do k=1,km
2376  dz(k) = dz(k) * (ztop/ze(1))
2377  enddo
2378 
2379  do k=km,1,-1
2380  ze(k) = ze(k+1) + dz(k)
2381  enddo
2382 ! ze(1) = ztop
2383 
2384  if ( is_master() ) write(*,*) 'var55_dz: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
2385  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
2386 
2387 ! Given z --> p
2388  do k=1,km
2389  dz(k) = ze(k) - ze(k+1)
2390  dlnp(k) = grav*dz(k) / (rdgas*t0)
2391  enddo
2392  do k=2,km
2393  peln(k) = peln(k-1) + dlnp(k-1)
2394  pe1(k) = exp(peln(k))
2395  enddo
2396 
2397 ! Pe(k) = ak(k) + bk(k) * PS
2398 ! Locate pint and KS
2399  ks = 0
2400  do k=2,km
2401  if ( pint < pe1(k)) then
2402  ks = k-1
2403  exit
2404  endif
2405  enddo
2406  if ( is_master() ) then
2407  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
2408  endif
2409  pint = pe1(ks+1)
2410 
2411 #ifdef NO_UKMO_HB
2412  do k=1,ks+1
2413  ak(k) = pe1(k)
2414  bk(k) = 0.
2415  enddo
2416 
2417  do k=ks+2,km+1
2418  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
2419  ak(k) = pe1(k) - bk(k) * pe1(km+1)
2420  enddo
2421  bk(km+1) = 1.
2422  ak(km+1) = 0.
2423 #else
2424 ! Problematic for non-hydrostatic
2425  do k=1,km+1
2426  eta(k) = pe1(k) / pe1(km+1)
2427  enddo
2428 
2429  ep = eta(ks+1)
2430  es = eta(km)
2431 ! es = 1.
2432  alpha = (ep**2-2.*ep*es) / (es-ep)**2
2433  beta = 2.*ep*es**2 / (es-ep)**2
2434  gama = -(ep*es)**2 / (es-ep)**2
2435 
2436 ! Pure pressure:
2437  do k=1,ks+1
2438  ak(k) = eta(k)*1.e5
2439  bk(k) = 0.
2440  enddo
2441 
2442  do k=ks+2, km
2443  ak(k) = alpha*eta(k) + beta + gama/eta(k)
2444  ak(k) = ak(k)*1.e5
2445  enddo
2446  ak(km+1) = 0.
2447 
2448  do k=ks+2, km
2449  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
2450  enddo
2451  bk(km+1) = 1.
2452 #endif
2453 
2454  if ( is_master() ) then
2455  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
2456  do k=1,km
2457  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
2458  enddo
2459  tmp1 = ak(ks+1)
2460  do k=ks+1,km
2461  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
2462  enddo
2463  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
2464  endif
2465 
2466 
2467  end subroutine var55_dz
2468 
2469  subroutine hybrid_z_dz(km, dz, ztop, s_rate)
2470  integer, intent(in):: km
2471  real, intent(in):: s_rate ! between [1. 1.1]
2472  real, intent(in):: ztop
2473  real, intent(out):: dz(km)
2474 ! Local
2475  real, parameter:: p00 = 1.e5
2476  real, dimension(km+1):: ze, pe1, peln, eta
2477  real, dimension(km):: s_fac, dlnp
2478  real t0, dz0, sum1, tmp1
2479  real ep, es, alpha, beta, gama
2480  integer k
2481 
2482  s_fac(km ) = 0.12
2483  s_fac(km-1) = 0.20
2484  s_fac(km-2) = 0.30
2485  s_fac(km-3) = 0.40
2486  s_fac(km-4) = 0.50
2487  s_fac(km-5) = 0.60
2488  s_fac(km-6) = 0.70
2489  s_fac(km-7) = 0.80
2490  s_fac(km-8) = 0.90
2491  s_fac(km-9) = 1.
2492 
2493  do k=km-10, 9, -1
2494  s_fac(k) = min(4.0, s_rate * s_fac(k+1) )
2495  enddo
2496 
2497  s_fac(8) = 1.05*s_fac(9)
2498  s_fac(7) = 1.1 *s_fac(8)
2499  s_fac(6) = 1.15*s_fac(7)
2500  s_fac(5) = 1.2 *s_fac(6)
2501  s_fac(4) = 1.3 *s_fac(5)
2502  s_fac(3) = 1.4 *s_fac(4)
2503  s_fac(2) = 1.5 *s_fac(3)
2504  s_fac(1) = 1.6 *s_fac(2)
2505 
2506  sum1 = 0.
2507  do k=1,km
2508  sum1 = sum1 + s_fac(k)
2509  enddo
2510 
2511  dz0 = ztop / sum1
2512 
2513  do k=1,km
2514  dz(k) = s_fac(k) * dz0
2515  enddo
2516 
2517  ze(km+1) = 0.
2518  do k=km,1,-1
2519  ze(k) = ze(k+1) + dz(k)
2520  enddo
2521 
2522  ze(1) = ztop
2523 
2524  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
2525 
2526  do k=1,km
2527  dz(k) = ze(k) - ze(k+1)
2528  enddo
2529 
2530  end subroutine hybrid_z_dz
2531 
2532 
2533 
2534  subroutine get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
2535  integer, intent(in) :: npz
2536  real, intent(in) :: p_s ! unit: pascal
2537  real, intent(in) :: ak(npz+1)
2538  real, intent(in) :: bk(npz+1)
2539  real, intent(in), optional :: pscale
2540  real, intent(out) :: pf(npz)
2541  real, intent(out) :: ph(npz+1)
2542  integer k
2543 
2544  ph(1) = ak(1)
2545  do k=2,npz+1
2546  ph(k) = ak(k) + bk(k)*p_s
2547  enddo
2548 
2549  if ( present(pscale) ) then
2550  do k=1,npz+1
2551  ph(k) = pscale*ph(k)
2552  enddo
2553  endif
2554 
2555  if( ak(1) > 1.e-8 ) then
2556  pf(1) = (ph(2) - ph(1)) / log(ph(2)/ph(1))
2557  else
2558  pf(1) = (ph(2) - ph(1)) * kappa/(kappa+1.)
2559  endif
2560 
2561  do k=2,npz
2562  pf(k) = (ph(k+1) - ph(k)) / log(ph(k+1)/ph(k))
2563  enddo
2564 
2565  end subroutine get_eta_level
2566 
2567 
2568 
2569  subroutine compute_dz(km, ztop, dz)
2570 
2571  integer, intent(in):: km
2572  real, intent(in):: ztop ! try 50.E3
2573  real, intent(out):: dz(km)
2574 !------------------------------
2575  real ze(km+1), dzt(km)
2576  integer k
2577 
2578 
2579 ! ztop = 30.E3
2580  dz(1) = ztop / real(km)
2581  dz(km) = 0.5*dz(1)
2582 
2583  do k=2,km-1
2584  dz(k) = dz(1)
2585  enddo
2586 
2587 ! Top:
2588  dz(1) = 2.*dz(2)
2589 
2590  ze(km+1) = 0.
2591  do k=km,1,-1
2592  ze(k) = ze(k+1) + dz(k)
2593  enddo
2594 
2595  if ( is_master() ) then
2596  write(*,*) 'Hybrid_z: dz, zm'
2597  do k=1,km
2598  dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000.
2599  write(*,*) k, dz(k), dzt(k)
2600  enddo
2601  endif
2602 
2603  end subroutine compute_dz
2604 
2605  subroutine compute_dz_var(km, ztop, dz)
2606 
2607  integer, intent(in):: km
2608  real, intent(in):: ztop ! try 50.E3
2609  real, intent(out):: dz(km)
2610 !------------------------------
2611  real, parameter:: s_rate = 1.0
2612  real ze(km+1)
2613  real s_fac(km)
2614  real sum1, dz0
2615  integer k
2616 
2617  s_fac(km ) = 0.125
2618  s_fac(km-1) = 0.20
2619  s_fac(km-2) = 0.30
2620  s_fac(km-3) = 0.40
2621  s_fac(km-4) = 0.50
2622  s_fac(km-5) = 0.60
2623  s_fac(km-6) = 0.70
2624  s_fac(km-7) = 0.80
2625  s_fac(km-8) = 0.90
2626  s_fac(km-9) = 1.
2627 
2628  do k=km-10, 9, -1
2629  s_fac(k) = s_rate * s_fac(k+1)
2630  enddo
2631 
2632  s_fac(8) = 1.05*s_fac(9)
2633  s_fac(7) = 1.1 *s_fac(8)
2634  s_fac(6) = 1.15*s_fac(7)
2635  s_fac(5) = 1.2 *s_fac(6)
2636  s_fac(4) = 1.3 *s_fac(5)
2637  s_fac(3) = 1.4 *s_fac(4)
2638  s_fac(2) = 1.5 *s_fac(3)
2639  s_fac(1) = 1.6 *s_fac(2)
2640 
2641  sum1 = 0.
2642  do k=1,km
2643  sum1 = sum1 + s_fac(k)
2644  enddo
2645 
2646  dz0 = ztop / sum1
2647 
2648  do k=1,km
2649  dz(k) = s_fac(k) * dz0
2650  enddo
2651 
2652  ze(1) = ztop
2653  ze(km+1) = 0.
2654  do k=km,2,-1
2655  ze(k) = ze(k+1) + dz(k)
2656  enddo
2657 
2658 ! Re-scale dz with the stretched ztop
2659  do k=1,km
2660  dz(k) = dz(k) * (ztop/ze(1))
2661  enddo
2662 
2663  do k=km,2,-1
2664  ze(k) = ze(k+1) + dz(k)
2665  enddo
2666 
2667  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
2668 
2669  do k=1,km
2670  dz(k) = ze(k) - ze(k+1)
2671  enddo
2672 
2673  end subroutine compute_dz_var
2674 
2675  subroutine compute_dz_l32(km, ztop, dz)
2676 
2677  integer, intent(in):: km
2678  real, intent(out):: dz(km)
2679  real, intent(out):: ztop ! try 50.E3
2680 !------------------------------
2681  real dzt(km)
2682  real ze(km+1)
2683  real dz0, dz1, dz2
2684  real z0, z1, z2
2685  integer k, k0, k1, k2, n
2686 
2687 !-------------------
2688  k2 = 8
2689  z2 = 30.e3
2690 !-------------------
2691  k1 = 21
2692  z1 = 10.0e3
2693 !-------------------
2694  k0 = 2
2695  z0 = 0.
2696  dz0 = 75. ! meters
2697 !-------------------
2698 ! Treat the surface layer as a special layer
2699  ze(1) = z0
2700  dz(1) = dz0
2701 
2702  ze(2) = dz(1)
2703  dz0 = 1.5*dz0
2704  dz(2) = dz0
2705 
2706  ze(3) = ze(2) + dz(2)
2707 
2708  dz1 = 2.*(z1-ze(3) - k1*dz0) / (k1*(k1-1))
2709 
2710  do k=k0+1,k0+k1
2711  dz(k) = dz0 + (k-k0)*dz1
2712  ze(k+1) = ze(k) + dz(k)
2713  enddo
2714 
2715  dz0 = dz(k1+k0)
2716  dz2 = 2.*(z2-ze(k0+k1+1)-k2*dz0) / (k2*(k2-1))
2717 
2718  do k=k0+k1+1,k0+k1+k2
2719  dz(k) = dz0 + (k-k0-k1)*dz2
2720  ze(k+1) = ze(k) + dz(k)
2721  enddo
2722 
2723  dz(km) = 2.*dz(km-1)
2724  ztop = ze(km) + dz(km)
2725  ze(km+1) = ze(km) + dz(km)
2726 
2727  call zflip (dz, 1, km)
2728 
2729  ze(km+1) = 0.
2730  do k=km,1,-1
2731  ze(k) = ze(k+1) + dz(k)
2732  enddo
2733 
2734 ! if ( is_master() ) then
2735 ! write(*,*) 'Hybrid_z: dz, zm'
2736 ! do k=1,km
2737 ! dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000.
2738 ! write(*,*) k, dz(k), dzt(k)
2739 ! enddo
2740 ! endif
2741 
2742  end subroutine compute_dz_l32
2743 
2744  subroutine compute_dz_l101(km, ztop, dz)
2745 
2746  integer, intent(in):: km ! km==101
2747  real, intent(out):: dz(km)
2748  real, intent(out):: ztop ! try 30.E3
2749 !------------------------------
2750  real ze(km+1)
2751  real dz0, dz1
2752  real:: stretch_f = 1.16
2753  integer k, k0, k1
2754 
2755  k1 = 2
2756  k0 = 25
2757  dz0 = 40. ! meters
2758 
2759  ze(km+1) = 0.
2760 
2761  do k=km, k0, -1
2762  dz(k) = dz0
2763  ze(k) = ze(k+1) + dz(k)
2764  enddo
2765 
2766  do k=k0+1, k1, -1
2767  dz(k) = stretch_f * dz(k+1)
2768  ze(k) = ze(k+1) + dz(k)
2769  enddo
2770 
2771  dz(1) = 4.0*dz(2)
2772  ze(1) = ze(2) + dz(1)
2773  ztop = ze(1)
2774 
2775  if ( is_master() ) then
2776  write(*,*) 'Hybrid_z: dz, ze'
2777  do k=1,km
2778  write(*,*) k, 0.001*dz(k), 0.001*ze(k)
2779  enddo
2780 ! ztop (km) = 20.2859154
2781  write(*,*) 'ztop (km) =', ztop * 0.001
2782  endif
2783 
2784  end subroutine compute_dz_l101
2785 
2786  subroutine set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
2787 
2788  integer, intent(in):: is, ie, js, je, ng, km
2789  real, intent(in):: rgrav, ztop
2790  real, intent(in):: dz(km) ! Reference vertical resolution for zs=0
2791  real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng)
2792  real, intent(inout):: ze(is:ie,js:je,km+1)
2793  real, optional, intent(out):: dz3(is-ng:ie+ng,js-ng:je+ng,km)
2794 ! local
2795  logical:: filter_xy = .false.
2796  real, allocatable:: delz(:,:,:)
2797  integer ntimes
2798  real zint
2799  real:: z1(is:ie,js:je)
2800  real:: z(km+1)
2801  real sig, z_rat
2802  integer ks(is:ie,js:je)
2803  integer i, j, k, ks_min, kint
2804 
2805  z(km+1) = 0.
2806  do k=km,1,-1
2807  z(k) = z(k+1) + dz(k)
2808  enddo
2809 
2810  do j=js,je
2811  do i=is,ie
2812  ze(i,j, 1) = ztop
2813  ze(i,j,km+1) = hs(i,j) * rgrav
2814  enddo
2815  enddo
2816 
2817  do k=2,km
2818  do j=js,je
2819  do i=is,ie
2820  ze(i,j,k) = z(k)
2821  enddo
2822  enddo
2823  enddo
2824 
2825 ! Set interface:
2826 #ifndef USE_VAR_ZINT
2827  zint = 12.0e3
2828  ntimes = 2
2829  kint = 2
2830  do k=2,km
2831  if ( z(k)<=zint ) then
2832  kint = k
2833  exit
2834  endif
2835  enddo
2836 
2837  if ( is_master() ) write(*,*) 'Z_coord interface set at k=',kint, ' ZE=', z(kint)
2838 
2839  do j=js,je
2840  do i=is,ie
2841  z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1))
2842  do k=km,kint+1,-1
2843  ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat
2844  enddo
2845 !--------------------------------------
2846 ! Apply vertical smoother locally to dz
2847 !--------------------------------------
2848  call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2849  enddo
2850  enddo
2851 #else
2852 ! ZINT is a function of local terrain
2853  ntimes = 4
2854  do j=js,je
2855  do i=is,ie
2856  z1(i,j) = dim(ze(i,j,km+1), 2500.) + 7500.
2857  enddo
2858  enddo
2859 
2860  ks_min = km
2861  do j=js,je
2862  do i=is,ie
2863  do k=km,2,-1
2864  if ( z(k)>=z1(i,j) ) then
2865  ks(i,j) = k
2866  ks_min = min(ks_min, k)
2867  go to 555
2868  endif
2869  enddo
2870 555 continue
2871  enddo
2872  enddo
2873 
2874  do j=js,je
2875  do i=is,ie
2876  kint = ks(i,j) + 1
2877  z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1))
2878  do k=km,kint+1,-1
2879  ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat
2880  enddo
2881 !--------------------------------------
2882 ! Apply vertical smoother locally to dz
2883 !--------------------------------------
2884  call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2885  enddo
2886  enddo
2887 #endif
2888 
2889 #ifdef DEV_ETA
2890  if ( filter_xy ) then
2891  allocate (delz(isd:ied, jsd:jed, km) )
2892  ntimes = 2
2893  do k=1,km
2894  do j=js,je
2895  do i=is,ie
2896  delz(i,j,k) = ze(i,j,k+1) - ze(i,j,k)
2897  enddo
2898  enddo
2899  enddo
2900  call del2_cubed(delz, 0.2*da_min, npx, npy, km, ntimes)
2901  do k=km,2,-1
2902  do j=js,je
2903  do i=is,ie
2904  ze(i,j,k) = ze(i,j,k+1) - delz(i,j,k)
2905  enddo
2906  enddo
2907  enddo
2908  deallocate ( delz )
2909  endif
2910 #endif
2911  if ( present(dz3) ) then
2912  do k=1,km
2913  do j=js,je
2914  do i=is,ie
2915  dz3(i,j,k) = ze(i,j,k+1) - ze(i,j,k)
2916  enddo
2917  enddo
2918  enddo
2919  endif
2920 
2921  end subroutine set_hybrid_z
2922 
2923 
2924  subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2925  integer, intent(in):: is, ie, js, je, km
2926  integer, intent(in):: ntimes, i, j
2927  real, intent(inout):: ze(is:ie,js:je,km+1)
2928 ! local:
2929  real, parameter:: df = 0.25
2930  real dz(km)
2931  real flux(km+1)
2932  integer k, n, k1, k2
2933 
2934  k2 = km-1
2935  do k=1,km
2936  dz(k) = ze(i,j,k+1) - ze(i,j,k)
2937  enddo
2938 
2939  do n=1,ntimes
2940  k1 = 2 + (ntimes-n)
2941 
2942  flux(k1 ) = 0.
2943  flux(k2+1) = 0.
2944  do k=k1+1,k2
2945  flux(k) = df*(dz(k) - dz(k-1))
2946  enddo
2947 
2948  do k=k1,k2
2949  dz(k) = dz(k) - flux(k) + flux(k+1)
2950  enddo
2951  enddo
2952 
2953  do k=km,1,-1
2954  ze(i,j,k) = ze(i,j,k+1) - dz(k)
2955  enddo
2956 
2957  end subroutine sm1_edge
2958 
2959 
2960 
2961  subroutine gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
2962  integer, intent(in):: km
2963  real, intent(in):: p0, ztop
2964  real, intent(inout):: ptop
2965  real, intent(inout):: ak(km+1), bk(km+1)
2966  real, intent(out):: pt1(km)
2967 ! Local
2968  logical:: isothermal
2969  real, dimension(km+1):: ze, pe1, pk1
2970  real, dimension(km):: dz1
2971  real t0, n2, s0
2972  integer k
2973 
2974 ! Set up vertical coordinare with constant del-z spacing:
2975  isothermal = .false.
2976  t0 = 300.
2977 
2978  if ( isothermal ) then
2979  n2 = grav**2/(cp_air*t0)
2980  else
2981  n2 = 0.0001
2982  endif
2983 
2984  s0 = grav*grav / (cp_air*n2)
2985 
2986  ze(km+1) = 0.
2987  do k=km,1,-1
2988  dz1(k) = ztop / real(km)
2989  ze(k) = ze(k+1) + dz1(k)
2990  enddo
2991 
2992 ! Given z --> p
2993  do k=1,km+1
2994  pe1(k) = p0*( (1.-s0/t0) + s0/t0*exp(-n2*ze(k)/grav) )**(1./kappa)
2995  enddo
2996 
2997  ptop = pe1(1)
2998 ! if ( is_master() ) write(*,*) 'GW_1D: computed model top (pa)=', ptop
2999 
3000 ! Set up "sigma" coordinate
3001  ak(1) = pe1(1)
3002  bk(1) = 0.
3003  do k=2,km
3004  bk(k) = (pe1(k) - pe1(1)) / (pe1(km+1)-pe1(1)) ! bk == sigma
3005  ak(k) = pe1(1)*(1.-bk(k))
3006  enddo
3007  ak(km+1) = 0.
3008  bk(km+1) = 1.
3009 
3010  do k=1,km+1
3011  pk1(k) = pe1(k) ** kappa
3012  enddo
3013 
3014 ! Compute volume mean potential temperature with hydrostatic eqn:
3015  do k=1,km
3016  pt1(k) = grav*dz1(k) / ( cp_air*(pk1(k+1)-pk1(k)) )
3017  enddo
3018 
3019  end subroutine gw_1d
3020 
3021 
3022 
3023  subroutine zflip(q, im, km)
3024  integer, intent(in):: im, km
3025  real, intent(inout):: q(im,km)
3026 !---
3027  integer i, k
3028  real qtmp
3029 
3030  do i = 1, im
3031  do k = 1, (km+1)/2
3032  qtmp = q(i,k)
3033  q(i,k) = q(i,km+1-k)
3034  q(i,km+1-k) = qtmp
3035  end do
3036  end do
3037 
3038  end subroutine zflip
3039 
3040 end module fv_eta_nlm_mod
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
constants_mod
Definition: constants.F90:24
fv_eta_nlm_mod::sm1_edge
subroutine, public sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
Definition: fv_eta_nlm.F90:2925
fv_eta_nlm_mod::var_les
subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:1486
fv_eta_nlm_mod::set_eta
subroutine, public set_eta(km, ks, ptop, ak, bk)
Definition: fv_eta_nlm.F90:392
fv_eta_nlm_mod::var_gfs
subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:1648
fv_eta_nlm_mod::var_dz
subroutine var_dz(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:2147
constants_mod::rdgas
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
Definition: constants.F90:77
fv_eta_nlm_mod::compute_dz
subroutine, public compute_dz(km, ztop, dz)
Definition: fv_eta_nlm.F90:2570
fv_eta_nlm_mod::gw_1d
subroutine, public gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
Definition: fv_eta_nlm.F90:2962
fv_eta_nlm_mod::get_eta_level
subroutine, public get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
Definition: fv_eta_nlm.F90:2535
mpp_mod
Definition: mpp.F90:39
fv_eta_nlm_mod::hybrid_z_dz
subroutine, public hybrid_z_dz(km, dz, ztop, s_rate)
Definition: fv_eta_nlm.F90:2470
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_eta_nlm_mod::var55_dz
subroutine var55_dz(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:2308
fv_eta_nlm_mod::var_hi
subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:1812
fv_eta_nlm_mod
Definition: fv_eta_nlm.F90:20
fv_eta_nlm_mod::compute_dz_var
subroutine, public compute_dz_var(km, ztop, dz)
Definition: fv_eta_nlm.F90:2606
constants_mod::grav
real, parameter, public grav
Acceleration due to gravity [m/s^2].
Definition: constants.F90:76
mpp_mod::mpp_error
Definition: mpp.F90:402
fv_eta_nlm_mod::var_hi2
subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta_nlm.F90:1988
max
#define max(a, b)
Definition: mosaic_util.h:33
fv_eta_nlm_mod::set_hybrid_z
subroutine, public set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
Definition: fv_eta_nlm.F90:2787
fv_eta_nlm_mod::compute_dz_l101
subroutine, public compute_dz_l101(km, ztop, dz)
Definition: fv_eta_nlm.F90:2745
fv_eta_nlm_mod::compute_dz_l32
subroutine, public compute_dz_l32(km, ztop, dz)
Definition: fv_eta_nlm.F90:2676
min
#define min(a, b)
Definition: mosaic_util.h:32
constants_mod::kappa
real, parameter, public kappa
RDGAS / CP_AIR [dimensionless].
Definition: constants.F90:82
fv_eta_nlm_mod::zflip
subroutine zflip(q, im, km)
Definition: fv_eta_nlm.F90:3024