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