FV3DYCORE  Version 2.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 SHiELD
58 !!! This routine will be kept here
59 !!! for the time being to not disrupt idealized tests
60 #ifdef USE_VAR_ETA
61  subroutine set_eta(km, ks, ptop, ak, bk, npz_type)
62 ! This is the easy to use version of the set_eta
63  integer, intent(in):: km ! vertical dimension
64  integer, intent(out):: ks ! number of pure p layers
65  real:: a60(61),b60(61)
66 ! Thfollowing L63 setting is the same as NCEP GFS's L64 except the top
67 ! 3 layers
68  data a60/300.0000, 430.00000, 558.00000, &
69  700.00000, 863.05803, 1051.07995, &
70  1265.75194, 1510.71101, 1790.05098, &
71  2108.36604, 2470.78817, 2883.03811, &
72  3351.46002, 3883.05187, 4485.49315, &
73  5167.14603, 5937.04991, 6804.87379, &
74  7780.84698, 8875.64338, 10100.20534, &
75  11264.35673, 12190.64366, 12905.42546, &
76  13430.87867, 13785.88765, 13986.77987, &
77  14047.96335, 13982.46770, 13802.40331, &
78  13519.33841, 13144.59486, 12689.45608, &
79  12165.28766, 11583.57006, 10955.84778, &
80  10293.60402, 9608.08306, 8910.07678, &
81  8209.70131, 7516.18560, 6837.69250, &
82  6181.19473, 5552.39653, 4955.72632, &
83  4394.37629, 3870.38682, 3384.76586, &
84  2937.63489, 2528.37666, 2155.78385, &
85  1818.20722, 1513.68173, 1240.03585, &
86  994.99144, 776.23591, 581.48797, &
87  408.53400, 255.26520, 119.70243, 0. /
88 
89  data b60/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.00000, 0.00000, 0.00000, &
95  0.00000, 0.00000, 0.00000, &
96  0.00201, 0.00792, 0.01755, &
97  0.03079, 0.04751, 0.06761, &
98  0.09097, 0.11746, 0.14690, &
99  0.17911, 0.21382, 0.25076, &
100  0.28960, 0.32994, 0.37140, &
101  0.41353, 0.45589, 0.49806, &
102  0.53961, 0.58015, 0.61935, &
103  0.65692, 0.69261, 0.72625, &
104  0.75773, 0.78698, 0.81398, &
105  0.83876, 0.86138, 0.88192, &
106  0.90050, 0.91722, 0.93223, &
107  0.94565, 0.95762, 0.96827, &
108  0.97771, 0.98608, 0.99347, 1./
109  real, intent(out):: ak(km+1)
110  real, intent(out):: bk(km+1)
111  real, intent(out):: ptop ! model top (Pa)
112  character(24), intent(IN) :: npz_type
113  real pint, stretch_fac
114  integer k
115  real :: s_rate = -1.0 ! dummy value to not use var_les
116 
117  pint = 100.e2
118 
119 !- Notes ---------------------------------
120 ! low-top: ptop = 100. ! ~45 km
121 ! mid-top: ptop = 10. ! ~60 km
122 ! hi -top: ptop = 1. ! ~80 km
123 !-----------------------------------------
124  select case (km)
125 
126 ! Optimal number = 8 * N -1 (for 8 openMP threads)
127 ! 16 * M -1 (for 16 openMP threads)
128 
129 #ifdef HIWPP
130 #ifdef SUPER_K
131  case (20)
132  ptop = 56.e2
133  pint = ptop
134  stretch_fac = 1.03
135  case (24)
136  ptop = 56.e2
137  pint = ptop
138  stretch_fac = 1.03
139  case (30)
140  ptop = 56.e2
141  pint = ptop
142  stretch_fac = 1.03
143  case (40)
144  ptop = 56.e2
145  pint = ptop
146  stretch_fac = 1.03
147  case (50)
148  ptop = 56.e2
149  pint = ptop
150  stretch_fac = 1.03
151  case (60)
152  ptop = 56.e2
153  pint = ptop
154  stretch_fac = 1.03
155  case (80)
156  ptop = 56.e2
157  pint = ptop
158  stretch_fac = 1.03
159 #else
160  case (30) ! For Baroclinic Instability Test
161  ptop = 2.26e2
162  pint = 250.e2
163  stretch_fac = 1.03
164  case (40)
165  ptop = 50.e2 ! For super cell test
166  pint = 300.e2
167  stretch_fac = 1.03
168  case (50) ! Mountain waves?
169  ptop = 30.e2
170  stretch_fac = 1.025
171  case (60) ! For Baroclinic Instability Test
172 #ifdef GFSL60
173  ks = 20
174  ptop = a60(1)
175  pint = a60(ks+1)
176  do k=1,km+1
177  ak(k) = a60(k)
178  bk(k) = b60(k)
179  enddo
180 #else
181  ptop = 3.e2
182 ! pint = 250.E2
183  pint = 300.e2 ! revised for Moist test
184  stretch_fac = 1.03
185 #endif
186 #endif
187  case (64)
188 !!! ptop = 3.e2
189  ptop = 2.0e2
190  pint = 300.e2
191  stretch_fac = 1.03
192 #else
193 ! *Very-low top: for idealized super-cell simulation:
194  case (50)
195  ptop = 50.e2
196  pint = 250.e2
197  stretch_fac = 1.03
198  case (60)
199  ptop = 40.e2
200  pint = 250.e2
201  stretch_fac = 1.03
202  case (90) ! super-duper cell
203  ptop = 40.e2
204  stretch_fac = 1.025
205 #endif
206 ! Low-top:
207  case (31) ! N = 4, M=2
208  ptop = 100.
209  stretch_fac = 1.035
210  case (32) ! N = 4, M=2
211  ptop = 100.
212  stretch_fac = 1.035
213  case (39) ! N = 5
214  ptop = 100.
215  stretch_fac = 1.035
216  case (41)
217  ptop = 100.
218  stretch_fac = 1.035
219  case (47) ! N = 6, M=3
220  ptop = 100.
221  stretch_fac = 1.035
222  case (51)
223  ptop = 100.
224  stretch_fac = 1.03
225  case (52) ! very low top
226  ptop = 30.e2 ! for special DPM RCE experiments
227  stretch_fac = 1.03
228 ! Mid-top:
229  case (55) ! N = 7
230  ptop = 10.
231  stretch_fac = 1.035
232 ! Hi-top:
233  case (63) ! N = 8, M=4
234  ptop = 1.
235  ! c360 or c384
236  stretch_fac = 1.035
237  case (71) ! N = 9
238  ptop = 1.
239  stretch_fac = 1.03
240  case (79) ! N = 10, M=5
241  ptop = 1.
242  stretch_fac = 1.03
243  case (127) ! N = 10, M=5
244  ptop = 1.
245  stretch_fac = 1.03
246  case (151)
247  ptop = 75.e2
248  pint = 500.e2
249  s_rate = 1.01
250  case default
251  ptop = 1.
252  stretch_fac = 1.03
253  end select
254 
255 #ifdef MOUNTAIN_WAVES
256  call mount_waves(km, ak, bk, ptop, ks, pint)
257 #else
258  if (s_rate > 0.) then
259  call var_les(km, ak, bk, ptop, ks, pint, s_rate)
260  else
261  if ( km > 79 ) then
262  call var_hi2(km, ak, bk, ptop, ks, pint, stretch_fac)
263  elseif (km==5 .or. km==10 ) then
264 ! Equivalent Shallow Water: for NGGPS, variable-resolution testing
265  ptop = 500.e2
266  ks = 0
267  do k=1,km+1
268  bk(k) = real(k-1) / real (km)
269  ak(k) = ptop*(1.-bk(k))
270  enddo
271  else
272 #ifndef GFSL60
273  call var_hi(km, ak, bk, ptop, ks, pint, stretch_fac)
274 #endif
275  endif
276 #endif
277  endif
278 
279  ptop = ak(1)
280  pint = ak(ks+1)
281 
282  end subroutine set_eta
283 
284 
285 #else
286  !This is the version of set_eta used in SHiELD and AM4
287  subroutine set_eta(km, ks, ptop, ak, bk, npz_type)
288 
289 !Level definitions are now in this header file
290 #include <tools/fv_eta.h>
291 
292  integer, intent(in):: km ! vertical dimension
293  integer, intent(out):: ks ! number of pure p layers
294  real, intent(out):: ak(km+1)
295  real, intent(out):: bk(km+1)
296  real, intent(out):: ptop ! model top (Pa)
297  character(24), intent(IN) :: npz_type
298 
299  real:: p0=1000.e2
300  real:: pc=200.e2
301 
302  real pt, lnpe, dlnp
303  real press(km+1), pt1(km)
304  integer k
305  integer :: var_fn = 0
306 
307  real :: pint = 100.e2
308  real :: stretch_fac = 1.03
309  integer :: auto_routine = 0
310 
311 
312  ptop = 1.
313 
314  ! Definition: press(i,j,k) = ak(k) + bk(k) * ps(i,j)
315 
316  if (trim(npz_type) == 'superC' .or. trim(npz_type) == 'superK') then
317 
318  auto_routine = 1
319  select case (km)
320  case (20)
321  ptop = 56.e2
322  pint = ptop
323  stretch_fac = 1.03
324  case (24)
325  ptop = 56.e2
326  pint = ptop
327  stretch_fac = 1.03
328  case (30)
329  ptop = 56.e2
330  pint = ptop
331  stretch_fac = 1.03
332  case (40)
333  ptop = 56.e2
334  pint = ptop
335  stretch_fac = 1.03
336  case (50)
337  ptop = 56.e2
338  pint = ptop
339  stretch_fac = 1.03
340  case (60)
341  ptop = 56.e2
342  pint = ptop
343  stretch_fac = 1.03
344  case (80)
345  ptop = 56.e2
346  pint = ptop
347  stretch_fac = 1.03
348  case (90) ! super-duper cell
349  ptop = 40.e2
350  stretch_fac = 1.025
351  auto_routine = 2
352  end select
353 
354  else
355 
356  select case (km)
357 
358  case (5,10) ! does this work????
359 
360  ! Equivalent Shallow Water: for modon test
361  ptop = 500.e2
362  ks = 0
363  do k=1,km+1
364  bk(k) = real(k-1) / real (km)
365  ak(k) = ptop*(1.-bk(k))
366  enddo
367 
368  case (24)
369 
370  ks = 5
371  do k=1,km+1
372  ak(k) = a24(k)
373  bk(k) = b24(k)
374  enddo
375 
376  case (26)
377 
378  ks = 7
379  do k=1,km+1
380  ak(k) = a26(k)
381  bk(k) = b26(k)
382  enddo
383 
384  case (30) ! For Baroclinic Instability Test
385  ptop = 2.26e2
386  pint = 250.e2
387  stretch_fac = 1.03
388  auto_routine = 1
389 
390  case (31) ! N = 4, M=2
391  if (trim(npz_type) == 'lowtop') then
392  ptop = 300.
393  pint = 100.e2
394  stretch_fac = 1.035
395  auto_routine = 5
396  else
397  ptop = 100.
398  stretch_fac = 1.035
399  auto_routine = 1
400  endif
401 
402  case (32)
403 
404  if (trim(npz_type) == 'old32') then
405  ks = 13 ! high-res trop_32 setup
406  do k=1,km+1
407  ak(k) = a32old(k)
408  bk(k) = b32old(k)
409  enddo
410  elseif (trim(npz_type) == 'lowtop') then
411  ptop = 100.
412  stretch_fac = 1.035
413  auto_routine = 1
414  else
415  ks = 7
416  do k=1,km+1
417  ak(k) = a32(k)
418  bk(k) = b32(k)
419  enddo
420  endif
421  !miz
422  case (33)
423  ks = 7
424  do k=1,km+1
425  ak(k) = a33(k)
426  bk(k) = b33(k)
427  enddo
428  !miz
429 
430  case (39) ! N = 5
431  ptop = 100.
432  stretch_fac = 1.035
433  auto_routine = 1
434 
435  case (40)
436  ptop = 50.e2 ! For super cell test
437  pint = 300.e2
438  stretch_fac = 1.03
439  auto_routine = 1
440 
441  case (41)
442  ptop = 100.
443  pint = 100.e2
444  stretch_fac = 1.035
445  auto_routine = 1
446 
447  case (47)
448 
449  if (trim(npz_type) == 'lowtop') then
450  ptop = 100.
451  stretch_fac = 1.035
452  auto_routine = 1
453  else
454  ! ks = 27 ! high-res trop-strat
455  ks = 20 ! Oct 23, 2012
456  do k=1,km+1
457  ak(k) = a47(k)
458  bk(k) = b47(k)
459  enddo
460  endif
461 
462  case (48)
463  ks = 28
464  do k=1,km+1
465  ak(k) = a48(k)
466  bk(k) = b48(k)
467  enddo
468 
469  case (50)
470  ! *Very-low top: for idealized super-cell simulation:
471  ptop = 50.e2
472  pint = 250.e2
473  stretch_fac = 1.03
474  auto_routine = 1
475 
476  case (51)
477  if (trim(npz_type) == 'lowtop') then
478  ptop = 100.
479  stretch_fac = 1.03
480  auto_routine = 1
481  elseif (trim(npz_type) == 'meso') then
482  ptop = 20.e2
483  pint = 100.e2
484  stretch_fac = 1.05
485  auto_routine = 1
486  elseif (trim(npz_type) == 'meso2') then
487  ptop = 1.e2
488  pint = 100.e2
489  stretch_fac = 1.05
490  auto_routine = 6
491  else
492  ptop = 100.
493  pint = 100.e2
494  stretch_fac = 1.035
495  auto_routine = 1
496  endif
497 
498  case (52)
499 
500  if (trim(npz_type) == 'rce') then
501  ptop = 30.e2 ! for special DPM RCE experiments
502  stretch_fac = 1.03
503  auto_routine = 1
504  else
505  ks = 35 ! pint = 223
506  do k=1,km+1
507  ak(k) = a52(k)
508  bk(k) = b52(k)
509  enddo
510  endif
511 
512  case (54)
513  ks = 11 ! pint = 109.4
514  do k=1,km+1
515  ak(k) = a54(k)
516  bk(k) = b54(k)
517  enddo
518 
519  ! Mid-top:
520  case (55) ! N = 7
521  ptop = 10.
522  pint = 100.e2
523  stretch_fac = 1.035
524  auto_routine = 1
525 
526  case (56)
527  ks = 26
528  do k=1,km+1
529  ak(k) = a56(k)
530  bk(k) = b56(k)
531  enddo
532 
533  case (60)
534 
535  if (trim(npz_type) == 'gfs') then
536  ks = 20
537  do k=1,km+1
538  ak(k) = a60gfs(k)
539  bk(k) = b60gfs(k)
540  enddo
541  else if (trim(npz_type) == 'BCwave') then
542  ptop = 3.e2
543  ! pint = 250.E2
544  pint = 300.e2 ! revised for Moist test
545  stretch_fac = 1.03
546  auto_routine = 1
547  else if (trim(npz_type) == 'meso') then
548 
549  ptop = 40.e2
550  pint = 250.e2
551  stretch_fac = 1.03
552  auto_routine = 1
553 
554  else
555  ks = 19
556  do k=1,km+1
557  ak(k) = a60(k)
558  bk(k) = b60(k)
559  enddo
560  endif
561 
562  case (63)
563  if (trim(npz_type) == 'meso') then
564  ks = 11
565  do k=1,km+1
566  ak(k) = a63meso(k)
567  bk(k) = b63meso(k)
568  enddo
569  elseif (trim(npz_type) == 'hitop') then
570  ptop = 1. ! high top
571  pint = 100.e2
572  stretch_fac = 1.035
573  auto_routine = 1
574  else!if (trim(npz_type) == 'gfs') then
575  !Used for SHiELD
576  ! GFS L64 equivalent setting
577  ks = 23
578  do k=1,km+1
579  ak(k) = a63(k)
580  bk(k) = b63(k)
581  enddo
582  endif
583 
584  case (64)
585 
586  if (trim(npz_type) == 'gfs') then
587  ks = 23
588  do k=1,km+1
589  ak(k) = a64gfs(k)
590  bk(k) = b64gfs(k)
591  enddo
592 
593  else
594 
595  ks = 46
596  do k=1,km+1
597  ak(k) = a64(k)
598  bk(k) = b64(k)
599  enddo
600 
601  endif
602  !-->cjg
603  case (68)
604  ks = 27
605  do k=1,km+1
606  ak(k) = a68(k)
607  bk(k) = b68(k)
608  enddo
609 
610  case (71) ! N = 9
611  ptop = 1.
612  stretch_fac = 1.03
613  auto_routine = 1
614  case (75) ! HS-SGO test configuration
615  pint = 100.e2
616  ptop = 10.e2
617  stretch_fac = 1.035
618  auto_routine = 6
619  case (79) ! N = 10, M=5
620  if (trim(npz_type) == 'gcrm') then
621  pint = 100.e2
622  ptop = 3.e2
623  stretch_fac = 1.035
624  auto_routine = 6
625  else
626  ptop = 1.
627  stretch_fac = 1.03
628  auto_routine = 1
629  endif
630  case (90) ! super-duper cell
631  ptop = 40.e2
632  stretch_fac = 1.025
633  auto_routine = 2
634 
635  ! NGGPS_GFS
636  case (91)
637  pint = 100.e2
638  ptop = 40.
639  stretch_fac = 1.029
640  auto_routine = 6
641 
642  case (95)
643  ! Mid-top settings:
644  pint = 100.e2
645  ptop = 20.
646  stretch_fac = 1.029
647  auto_routine = 6
648 
649  case (96)
650  ks = 27
651  do k=1,km+1
652  ak(k) = a96(k)
653  bk(k) = b96(k)
654  enddo
655 
656 
657  case (100)
658  ks = 38
659  do k=1,km+1
660  ak(k) = a100(k)
661  bk(k) = b100(k)
662  enddo
663 
664  case (104)
665  ks = 73
666  do k=1,km+1
667  ak(k) = a104(k)
668  bk(k) = b104(k)
669  enddo
670 
671  ! IFS-like L125
672  case (125)
673  ks = 33
674  ptop = a125(1)
675  pint = a125(ks+1)
676  do k=1,km+1
677  ak(k) = a125(k)
678  bk(k) = b125(k)
679  enddo
680 
681  case (127) ! N = 10, M=5
682  if (trim(npz_type) == 'hitop') then
683  ptop = 1.
684  stretch_fac = 1.03
685  auto_routine = 2
686  else
687  ptop = 1.
688  pint = 75.e2
689  stretch_fac = 1.028
690  auto_routine = 6
691  endif
692  case (151)
693  !LES applications
694  ptop = 75.e2
695  pint = 500.e2
696  stretch_fac = 1.01
697  auto_routine = 3
698 
699  case default
700 
701  if(trim(npz_type) == 'hitop') then
702  ptop = 1.
703  pint = 100.e2
704  elseif(trim(npz_type) == 'midtop') then
705  ptop = 10.
706  pint = 100.e2
707  elseif(trim(npz_type) == 'lowtop') then
708  ptop = 1.e2
709  pint = 100.e2
710  endif
711 
712  if (trim(npz_type) == 'gfs') then
713  auto_routine = 6
714  elseif(trim(npz_type) == 'les') then
715  auto_routine = 3
716  elseif(trim(npz_type) == 'mountain_wave') then
717  auto_routine = 4
718  elseif (km > 79) then
719  auto_routine = 2
720  else
721  auto_routine = 1
722  endif
723 
724  end select
725 
726  endif ! superC/superK
727 
728  select case (auto_routine)
729 
730  case (1)
731  call var_hi(km, ak, bk, ptop, ks, pint, stretch_fac)
732  case (2)
733  call var_hi2(km, ak, bk, ptop, ks, pint, stretch_fac)
734  case (3)
735  call var_les(km, ak, bk, ptop, ks, pint, stretch_fac)
736  case (4)
737  call mount_waves(km, ak, bk, ptop, ks, pint)
738  case (5)
739  call var_dz(km, ak, bk, ptop, ks, pint, stretch_fac)
740  case (6)
741  call var_gfs(km, ak, bk, ptop, ks, pint, stretch_fac)
742  end select
743 
744  ptop = ak(1)
745  pint = ak(ks+1)
746 
747  if (is_master()) then
748  write(*, '(A4, A13, A13, A11)') 'klev', 'ak', 'bk', 'p_ref'
749  do k=1,km+1
750  write(*,'(I4, F13.5, F13.5, F11.2)') k, ak(k), bk(k), 1000.e2*bk(k) + ak(k)
751  enddo
752  endif
753 
754 
755  end subroutine set_eta
756 #endif
757 
761  subroutine set_external_eta(ak, bk, ptop, ks)
762  implicit none
763  real, intent(in) :: ak(:)
764  real, intent(in) :: bk(:)
765  real, intent(out) :: ptop
766  integer, intent(out) :: ks
767  !--- local variables
768  integer :: k
769  real :: eps = 1.d-7
770 
771  ptop = ak(1)
772  ks = 1
773  do k = 1, size(bk(:))
774  if (bk(k).lt.eps) ks = k
775  enddo
776  !--- change ks to layers from levels
777  ks = ks - 1
778  if (is_master()) write(6,*) ' ptop & ks ', ptop, ks
779 
780  end subroutine set_external_eta
781 
782 
783  subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate)
784  implicit none
785  integer, intent(in):: km
786  real, intent(in):: ptop
787  real, intent(in):: s_rate
788  real, intent(out):: ak(km+1), bk(km+1)
789  real, intent(inout):: pint
790  integer, intent(out):: ks
791 ! Local
792  real, parameter:: p00 = 1.e5
793  real, dimension(km+1):: ze, pe1, peln, eta
794  real, dimension(km):: dz, s_fac, dlnp, pm, dp, dk
795  real ztop, t0, dz0, sum1, tmp1
796  real ep, es, alpha, beta, gama
797  real, parameter:: akap = 2./7.
798 !---- Tunable parameters:
799  real:: k_inc = 10
800  real:: s0 = 0.8
801 !-----------------------
802  real:: s_inc
803  integer k
804 
805  pe1(1) = ptop
806  peln(1) = log(pe1(1))
807  pe1(km+1) = p00
808  peln(km+1) = log(pe1(km+1))
809 
810  t0 = 273.
811  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
812 
813  s_inc = (1.-s0) / real(k_inc)
814  s_fac(km) = s0
815 
816  do k=km-1, km-k_inc, -1
817  s_fac(k) = s_fac(k+1) + s_inc
818  enddo
819 
820  s_fac(km-k_inc-1) = 0.5*(s_fac(km-k_inc) + s_rate)
821 
822  do k=km-k_inc-2, 5, -1
823  s_fac(k) = s_rate * s_fac(k+1)
824  enddo
825 
826  s_fac(4) = 0.5*(1.1+s_rate)*s_fac(5)
827  s_fac(3) = 1.1 *s_fac(4)
828  s_fac(2) = 1.1 *s_fac(3)
829  s_fac(1) = 1.1 *s_fac(2)
830 
831  sum1 = 0.
832  do k=1,km
833  sum1 = sum1 + s_fac(k)
834  enddo
835 
836  dz0 = ztop / sum1
837 
838  do k=1,km
839  dz(k) = s_fac(k) * dz0
840  enddo
841 
842  ze(km+1) = 0.
843  do k=km,1,-1
844  ze(k) = ze(k+1) + dz(k)
845  enddo
846 
847 ! Re-scale dz with the stretched ztop
848  do k=1,km
849  dz(k) = dz(k) * (ztop/ze(1))
850  enddo
851 
852  do k=km,1,-1
853  ze(k) = ze(k+1) + dz(k)
854  enddo
855 ! ze(1) = ztop
856 
857  if ( is_master() ) then
858  write(*,*) 'var_les: computed model top (m)=', ztop, ' bottom/top dz=', dz(km), dz(1)
859 ! do k=1,km
860 ! write(*,*) k, s_fac(k)
861 ! enddo
862  endif
863 
864  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
865 
866 ! Given z --> p
867  do k=1,km
868  dz(k) = ze(k) - ze(k+1)
869  dlnp(k) = grav*dz(k) / (rdgas*t0)
870  !write(*,*) k, dz(k)
871  enddo
872  do k=2,km
873  peln(k) = peln(k-1) + dlnp(k-1)
874  pe1(k) = exp(peln(k))
875  enddo
876 
877 
878 ! Pe(k) = ak(k) + bk(k) * PS
879 ! Locate pint and KS
880  ks = 0
881  do k=2,km
882  if ( pint < pe1(k)) then
883  ks = k-1
884  exit
885  endif
886  enddo
887  if ( is_master() ) then
888  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
889  endif
890  pint = pe1(ks+1)
891 
892  do k=1,km+1
893  eta(k) = pe1(k) / pe1(km+1)
894  enddo
895 
896  ep = eta(ks+1)
897  es = eta(km)
898 ! es = 1.
899  alpha = (ep**2-2.*ep*es) / (es-ep)**2
900  beta = 2.*ep*es**2 / (es-ep)**2
901  gama = -(ep*es)**2 / (es-ep)**2
902 
903 ! Pure pressure:
904  do k=1,ks+1
905  ak(k) = eta(k)*1.e5
906  bk(k) = 0.
907  enddo
908 
909  do k=ks+2, km
910  ak(k) = alpha*eta(k) + beta + gama/eta(k)
911  ak(k) = ak(k)*1.e5
912  enddo
913  ak(km+1) = 0.
914 
915  do k=ks+2, km
916  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
917  enddo
918  bk(km+1) = 1.
919 
920  if ( is_master() ) then
921  ! write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
922  ! do k=1,km
923  ! pm(k) = 0.5*(pe1(k)+pe1(k+1))/100.
924  ! write(*,*) k, pm(k), dz(k)
925  ! enddo
926  tmp1 = ak(ks+1)
927  do k=ks+1,km
928  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
929  enddo
930  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
931  write(*,800) (pm(k), k=km,1,-1)
932  endif
933 
934  do k=1,km
935  dp(k) = (pe1(k+1) - pe1(k))/100.
936  dk(k) = pe1(k+1)**akap - pe1(k)**akap
937  enddo
938 
939 800 format(1x,5(1x,f9.4))
940 
941  end subroutine var_les
942 
943 
944 
945  subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate)
946  integer, intent(in):: km
947  real, intent(in):: ptop
948  real, intent(in):: s_rate
949  real, intent(out):: ak(km+1), bk(km+1)
950  real, intent(inout):: pint
951  integer, intent(out):: ks
952 ! Local
953  real, parameter:: p00 = 1.e5
954  real, dimension(km+1):: ze, pe1, peln, eta
955  real, dimension(km):: dz, s_fac, dlnp
956  real ztop, t0, dz0, sum1, tmp1
957  real ep, es, alpha, beta, gama
958 !---- Tunable parameters:
959  integer:: k_inc = 25
960  real:: s0 = 0.13
961 !-----------------------
962  real:: s_inc
963  integer k
964 
965  pe1(1) = ptop
966  peln(1) = log(pe1(1))
967  pe1(km+1) = p00
968  peln(km+1) = log(pe1(km+1))
969 
970  t0 = 270.
971  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
972 
973  s_inc = (1.-s0) / real(k_inc)
974  s_fac(km) = s0
975 
976  do k=km-1, km-k_inc, -1
977  s_fac(k) = s_fac(k+1) + s_inc
978  enddo
979 
980  do k=km-k_inc-1, 9, -1
981  s_fac(k) = s_rate * s_fac(k+1)
982  enddo
983  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
984  s_fac(7) = 1.10*s_fac(8)
985  s_fac(6) = 1.15*s_fac(7)
986  s_fac(5) = 1.20*s_fac(6)
987  s_fac(4) = 1.26*s_fac(5)
988  s_fac(3) = 1.33*s_fac(4)
989  s_fac(2) = 1.41*s_fac(3)
990  s_fac(1) = 1.60*s_fac(2)
991 
992  sum1 = 0.
993  do k=1,km
994  sum1 = sum1 + s_fac(k)
995  enddo
996 
997  dz0 = ztop / sum1
998 
999  do k=1,km
1000  dz(k) = s_fac(k) * dz0
1001  enddo
1002 
1003  ze(km+1) = 0.
1004  do k=km,1,-1
1005  ze(k) = ze(k+1) + dz(k)
1006  enddo
1007 
1008 ! Re-scale dz with the stretched ztop
1009  do k=1,km
1010  dz(k) = dz(k) * (ztop/ze(1))
1011  enddo
1012 
1013  do k=km,1,-1
1014  ze(k) = ze(k+1) + dz(k)
1015  enddo
1016 ! ze(1) = ztop
1017 
1018  if ( is_master() ) then
1019  write(*,*) 'var_gfs: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1020 ! do k=1,km
1021 ! write(*,*) k, s_fac(k)
1022 ! enddo
1023  endif
1024 
1025 ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 3)
1026 
1027 ! Given z --> p
1028  do k=1,km
1029  dz(k) = ze(k) - ze(k+1)
1030  dlnp(k) = grav*dz(k) / (rdgas*t0)
1031  enddo
1032  do k=2,km
1033  peln(k) = peln(k-1) + dlnp(k-1)
1034  pe1(k) = exp(peln(k))
1035  enddo
1036 
1037 ! Pe(k) = ak(k) + bk(k) * PS
1038 ! Locate pint and KS
1039  ks = 0
1040  do k=2,km
1041  if ( pint < pe1(k)) then
1042  ks = k-1
1043  exit
1044  endif
1045  enddo
1046  if ( is_master() ) then
1047  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1048  write(*,*) 'ptop =', ptop
1049  endif
1050  pint = pe1(ks+1)
1051 
1052 #ifdef NO_UKMO_HB
1053  do k=1,ks+1
1054  ak(k) = pe1(k)
1055  bk(k) = 0.
1056  enddo
1057 
1058  do k=ks+2,km+1
1059  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1060  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1061  enddo
1062  bk(km+1) = 1.
1063  ak(km+1) = 0.
1064 #else
1065 ! Problematic for non-hydrostatic
1066  do k=1,km+1
1067  eta(k) = pe1(k) / pe1(km+1)
1068  enddo
1069 
1070  ep = eta(ks+1)
1071  es = eta(km)
1072 ! es = 1.
1073  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1074  beta = 2.*ep*es**2 / (es-ep)**2
1075  gama = -(ep*es)**2 / (es-ep)**2
1076 
1077 ! Pure pressure:
1078  do k=1,ks+1
1079  ak(k) = eta(k)*1.e5
1080  bk(k) = 0.
1081  enddo
1082 
1083  do k=ks+2, km
1084  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1085  ak(k) = ak(k)*1.e5
1086  enddo
1087  ak(km+1) = 0.
1088 
1089  do k=ks+2, km
1090  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1091  enddo
1092  bk(km+1) = 1.
1093 #endif
1094 
1095  if ( is_master() ) then
1096  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1097  do k=1,km
1098  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1099  enddo
1100  tmp1 = ak(ks+1)
1101  do k=ks+1,km
1102  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1103  enddo
1104  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1105  endif
1106 
1107  end subroutine var_gfs
1108 
1109  subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate)
1110  integer, intent(in):: km
1111  real, intent(in):: ptop
1112  real, intent(in):: s_rate
1113  real, intent(out):: ak(km+1), bk(km+1)
1114  real, intent(inout):: pint
1115  integer, intent(out):: ks
1116 ! Local
1117  real, parameter:: p00 = 1.e5
1118  real, dimension(km+1):: ze, pe1, peln, eta
1119  real, dimension(km):: dz, s_fac, dlnp
1120  real ztop, t0, dz0, sum1, tmp1
1121  real ep, es, alpha, beta, gama
1122 !---- Tunable parameters:
1123  integer:: k_inc = 15
1124  real:: s0 = 0.10
1125 !-----------------------
1126  real:: s_inc
1127  integer k
1128 
1129  pe1(1) = ptop
1130  peln(1) = log(pe1(1))
1131  pe1(km+1) = p00
1132  peln(km+1) = log(pe1(km+1))
1133 
1134  t0 = 270.
1135  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1136 
1137  s_inc = (1.-s0) / real(k_inc)
1138  s_fac(km) = s0
1139 
1140  do k=km-1, km-k_inc, -1
1141  s_fac(k) = s_fac(k+1) + s_inc
1142  enddo
1143 
1144  s_fac(km-k_inc-1) = 0.5*(s_fac(km-k_inc) + s_rate)
1145 
1146 #ifdef HIWPP
1147  do k=km-k_inc-2, 4, -1
1148  s_fac(k) = s_rate * s_fac(k+1)
1149  enddo
1150  s_fac(3) = 0.5*(1.15+s_rate)*s_fac(4)
1151  s_fac(2) = 1.15 *s_fac(3)
1152  s_fac(1) = 1.3 *s_fac(2)
1153 #else
1154  do k=km-k_inc-2, 9, -1
1155  s_fac(k) = s_rate * s_fac(k+1)
1156  enddo
1157 
1158  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
1159  s_fac(7) = 1.1 *s_fac(8)
1160  s_fac(6) = 1.15*s_fac(7)
1161  s_fac(5) = 1.2 *s_fac(6)
1162  s_fac(4) = 1.3 *s_fac(5)
1163  s_fac(3) = 1.4 *s_fac(4)
1164  s_fac(2) = 1.45 *s_fac(3)
1165  s_fac(1) = 1.5 *s_fac(2)
1166 #endif
1167 
1168  sum1 = 0.
1169  do k=1,km
1170  sum1 = sum1 + s_fac(k)
1171  enddo
1172 
1173  dz0 = ztop / sum1
1174 
1175  do k=1,km
1176  dz(k) = s_fac(k) * dz0
1177  enddo
1178 
1179  ze(km+1) = 0.
1180  do k=km,1,-1
1181  ze(k) = ze(k+1) + dz(k)
1182  enddo
1183 
1184 ! Re-scale dz with the stretched ztop
1185  do k=1,km
1186  dz(k) = dz(k) * (ztop/ze(1))
1187  enddo
1188 
1189  do k=km,1,-1
1190  ze(k) = ze(k+1) + dz(k)
1191  enddo
1192 ! ze(1) = ztop
1193 
1194  if ( is_master() ) then
1195  write(*,*) 'var_hi: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1196 ! do k=1,km
1197 ! write(*,*) k, s_fac(k)
1198 ! enddo
1199  endif
1200 
1201  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
1202 
1203 ! Given z --> p
1204  do k=1,km
1205  dz(k) = ze(k) - ze(k+1)
1206  dlnp(k) = grav*dz(k) / (rdgas*t0)
1207  enddo
1208  do k=2,km
1209  peln(k) = peln(k-1) + dlnp(k-1)
1210  pe1(k) = exp(peln(k))
1211  enddo
1212 
1213 ! Pe(k) = ak(k) + bk(k) * PS
1214 ! Locate pint and KS
1215  ks = 0
1216  do k=2,km
1217  if ( pint < pe1(k)) then
1218  ks = k-1
1219  exit
1220  endif
1221  enddo
1222  if ( is_master() ) then
1223  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1224  write(*,*) 'ptop =', ptop
1225  endif
1226  pint = pe1(ks+1)
1227 
1228 #ifdef NO_UKMO_HB
1229  do k=1,ks+1
1230  ak(k) = pe1(k)
1231  bk(k) = 0.
1232  enddo
1233 
1234  do k=ks+2,km+1
1235  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1236  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1237  enddo
1238  bk(km+1) = 1.
1239  ak(km+1) = 0.
1240 #else
1241 ! Problematic for non-hydrostatic
1242  do k=1,km+1
1243  eta(k) = pe1(k) / pe1(km+1)
1244  enddo
1245 
1246  ep = eta(ks+1)
1247  es = eta(km)
1248 ! es = 1.
1249  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1250  beta = 2.*ep*es**2 / (es-ep)**2
1251  gama = -(ep*es)**2 / (es-ep)**2
1252 
1253 ! Pure pressure:
1254  do k=1,ks+1
1255  ak(k) = eta(k)*1.e5
1256  bk(k) = 0.
1257  enddo
1258 
1259  do k=ks+2, km
1260  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1261  ak(k) = ak(k)*1.e5
1262  enddo
1263  ak(km+1) = 0.
1264 
1265  do k=ks+2, km
1266  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1267  enddo
1268  bk(km+1) = 1.
1269 #endif
1270 
1271  if ( is_master() ) then
1272  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1273  do k=1,km
1274  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1275  enddo
1276  tmp1 = ak(ks+1)
1277  do k=ks+1,km
1278  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1279  enddo
1280  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1281  endif
1282 
1283 
1284  end subroutine var_hi
1285  subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate)
1286  integer, intent(in):: km
1287  real, intent(in):: ptop
1288  real, intent(in):: s_rate
1289  real, intent(out):: ak(km+1), bk(km+1)
1290  real, intent(inout):: pint
1291  integer, intent(out):: ks
1292 ! Local
1293  real, parameter:: p00 = 1.e5
1294  real, dimension(km+1):: ze, pe1, peln, eta
1295  real, dimension(km):: dz, s_fac, dlnp
1296  real ztop, t0, dz0, sum1, tmp1
1297  real ep, es, alpha, beta, gama
1298  integer k
1299 
1300  pe1(1) = ptop
1301  peln(1) = log(pe1(1))
1302  pe1(km+1) = p00
1303  peln(km+1) = log(pe1(km+1))
1304 
1305  t0 = 270.
1306  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1307 
1308  s_fac(km ) = 0.15
1309  s_fac(km-1) = 0.20
1310  s_fac(km-2) = 0.30
1311  s_fac(km-3) = 0.40
1312  s_fac(km-4) = 0.50
1313  s_fac(km-5) = 0.60
1314  s_fac(km-6) = 0.70
1315  s_fac(km-7) = 0.80
1316  s_fac(km-8) = 0.90
1317  s_fac(km-9) = 0.95
1318  s_fac(km-10) = 0.5*(s_fac(km-9) + s_rate)
1319 
1320  do k=km-11, 8, -1
1321  s_fac(k) = s_rate * s_fac(k+1)
1322  enddo
1323 
1324  s_fac(7) = 0.5*(1.1+s_rate)*s_fac(9)
1325  s_fac(6) = 1.05*s_fac(7)
1326  s_fac(5) = 1.1*s_fac(6)
1327  s_fac(4) = 1.15*s_fac(5)
1328  s_fac(3) = 1.2*s_fac(4)
1329  s_fac(2) = 1.3*s_fac(3)
1330  s_fac(1) = 1.4*s_fac(2)
1331 
1332  sum1 = 0.
1333  do k=1,km
1334  sum1 = sum1 + s_fac(k)
1335  enddo
1336 
1337  dz0 = ztop / sum1
1338 
1339  do k=1,km
1340  dz(k) = s_fac(k) * dz0
1341  enddo
1342 
1343  ze(km+1) = 0.
1344  do k=km,1,-1
1345  ze(k) = ze(k+1) + dz(k)
1346  enddo
1347 
1348 ! Re-scale dz with the stretched ztop
1349  do k=1,km
1350  dz(k) = dz(k) * (ztop/ze(1))
1351  enddo
1352 
1353  do k=km,1,-1
1354  ze(k) = ze(k+1) + dz(k)
1355  enddo
1356 ! ze(1) = ztop
1357 
1358  if ( is_master() ) write(*,*) 'var_hi2: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1359  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
1360 
1361 ! Given z --> p
1362  do k=1,km
1363  dz(k) = ze(k) - ze(k+1)
1364  dlnp(k) = grav*dz(k) / (rdgas*t0)
1365  enddo
1366  do k=2,km
1367  peln(k) = peln(k-1) + dlnp(k-1)
1368  pe1(k) = exp(peln(k))
1369  enddo
1370 
1371 ! Pe(k) = ak(k) + bk(k) * PS
1372 ! Locate pint and KS
1373  ks = 0
1374  do k=2,km
1375  if ( pint < pe1(k)) then
1376  ks = k-1
1377  exit
1378  endif
1379  enddo
1380  if ( is_master() ) then
1381  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1382  endif
1383  pint = pe1(ks+1)
1384 
1385 #ifdef NO_UKMO_HB
1386  do k=1,ks+1
1387  ak(k) = pe1(k)
1388  bk(k) = 0.
1389  enddo
1390 
1391  do k=ks+2,km+1
1392  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1393  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1394  enddo
1395  bk(km+1) = 1.
1396  ak(km+1) = 0.
1397 #else
1398 ! Problematic for non-hydrostatic
1399  do k=1,km+1
1400  eta(k) = pe1(k) / pe1(km+1)
1401  enddo
1402 
1403  ep = eta(ks+1)
1404  es = eta(km)
1405 ! es = 1.
1406  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1407  beta = 2.*ep*es**2 / (es-ep)**2
1408  gama = -(ep*es)**2 / (es-ep)**2
1409 
1410 ! Pure pressure:
1411  do k=1,ks+1
1412  ak(k) = eta(k)*1.e5
1413  bk(k) = 0.
1414  enddo
1415 
1416  do k=ks+2, km
1417  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1418  ak(k) = ak(k)*1.e5
1419  enddo
1420  ak(km+1) = 0.
1421 
1422  do k=ks+2, km
1423  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1424  enddo
1425  bk(km+1) = 1.
1426 #endif
1427 
1428  if ( is_master() ) then
1429  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1430  do k=1,km
1431  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1432  enddo
1433  tmp1 = ak(ks+1)
1434  do k=ks+1,km
1435  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1436  enddo
1437  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1438  endif
1439 
1440 
1441  end subroutine var_hi2
1442 
1443 
1444  subroutine var_dz(km, ak, bk, ptop, ks, pint, s_rate)
1445  integer, intent(in):: km
1446  real, intent(in):: ptop
1447  real, intent(in):: s_rate
1448  real, intent(out):: ak(km+1), bk(km+1)
1449  real, intent(inout):: pint
1450  integer, intent(out):: ks
1451 ! Local
1452  real, parameter:: p00 = 1.e5
1453  real, dimension(km+1):: ze, pe1, peln, eta
1454  real, dimension(km):: dz, s_fac, dlnp
1455  real ztop, t0, dz0, sum1, tmp1
1456  real ep, es, alpha, beta, gama
1457  integer k
1458 
1459  pe1(1) = ptop
1460  peln(1) = log(pe1(1))
1461  pe1(km+1) = p00
1462  peln(km+1) = log(pe1(km+1))
1463 
1464  t0 = 270.
1465  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1466 
1467  s_fac(km ) = 0.10
1468  s_fac(km-1) = 0.20
1469  s_fac(km-2) = 0.30
1470  s_fac(km-3) = 0.40
1471  s_fac(km-4) = 0.50
1472  s_fac(km-5) = 0.60
1473  s_fac(km-6) = 0.70
1474  s_fac(km-7) = 0.80
1475  s_fac(km-8) = 0.90
1476  s_fac(km-9) = 0.95
1477  s_fac(km-10) = 0.5*(s_fac(km-9) + s_rate)
1478 
1479  do k=km-11, 9, -1
1480  s_fac(k) = min(10.0, s_rate * s_fac(k+1) )
1481  enddo
1482 
1483  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
1484  s_fac(7) = 1.1 *s_fac(8)
1485  s_fac(6) = 1.15*s_fac(7)
1486  s_fac(5) = 1.2 *s_fac(6)
1487  s_fac(4) = 1.3 *s_fac(5)
1488  s_fac(3) = 1.4 *s_fac(4)
1489  s_fac(2) = 1.5 *s_fac(3)
1490  s_fac(1) = 1.6 *s_fac(2)
1491 
1492  sum1 = 0.
1493  do k=1,km
1494  sum1 = sum1 + s_fac(k)
1495  enddo
1496 
1497  dz0 = ztop / sum1
1498 
1499  do k=1,km
1500  dz(k) = s_fac(k) * dz0
1501  enddo
1502 
1503  ze(km+1) = 0.
1504  do k=km,1,-1
1505  ze(k) = ze(k+1) + dz(k)
1506  enddo
1507 
1508 ! Re-scale dz with the stretched ztop
1509  do k=1,km
1510  dz(k) = dz(k) * (ztop/ze(1))
1511  enddo
1512 
1513  do k=km,1,-1
1514  ze(k) = ze(k+1) + dz(k)
1515  enddo
1516 ! ze(1) = ztop
1517 
1518  if ( is_master() ) write(*,*) 'var_dz: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1519  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
1520 
1521 ! Given z --> p
1522  do k=1,km
1523  dz(k) = ze(k) - ze(k+1)
1524  dlnp(k) = grav*dz(k) / (rdgas*t0)
1525  enddo
1526  do k=2,km
1527  peln(k) = peln(k-1) + dlnp(k-1)
1528  pe1(k) = exp(peln(k))
1529  enddo
1530 
1531 ! Pe(k) = ak(k) + bk(k) * PS
1532 ! Locate pint and KS
1533  ks = 0
1534  do k=2,km
1535  if ( pint < pe1(k)) then
1536  ks = k-1
1537  exit
1538  endif
1539  enddo
1540  if ( is_master() ) then
1541  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1542  write(*,*) 'ptop =', ptop
1543  endif
1544  pint = pe1(ks+1)
1545 
1546 #ifdef NO_UKMO_HB
1547  do k=1,ks+1
1548  ak(k) = pe1(k)
1549  bk(k) = 0.
1550  enddo
1551 
1552  do k=ks+2,km+1
1553  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1554  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1555  enddo
1556  bk(km+1) = 1.
1557  ak(km+1) = 0.
1558 #else
1559 ! Problematic for non-hydrostatic
1560  do k=1,km+1
1561  eta(k) = pe1(k) / pe1(km+1)
1562  enddo
1563 
1564  ep = eta(ks+1)
1565  es = eta(km)
1566 ! es = 1.
1567  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1568  beta = 2.*ep*es**2 / (es-ep)**2
1569  gama = -(ep*es)**2 / (es-ep)**2
1570 
1571 ! Pure pressure:
1572  do k=1,ks+1
1573  ak(k) = eta(k)*1.e5
1574  bk(k) = 0.
1575  enddo
1576 
1577  do k=ks+2, km
1578  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1579  ak(k) = ak(k)*1.e5
1580  enddo
1581  ak(km+1) = 0.
1582 
1583  do k=ks+2, km
1584  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1585  enddo
1586  bk(km+1) = 1.
1587 #endif
1588 
1589  if ( is_master() ) then
1590  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1591  do k=1,km
1592  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1593  enddo
1594  tmp1 = ak(ks+1)
1595  do k=ks+1,km
1596  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1597  enddo
1598  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1599  endif
1600 
1601 
1602  end subroutine var_dz
1603 
1604 
1605  subroutine var55_dz(km, ak, bk, ptop, ks, pint, s_rate)
1606  integer, intent(in):: km
1607  real, intent(in):: ptop
1608  real, intent(in):: s_rate
1609  real, intent(out):: ak(km+1), bk(km+1)
1610  real, intent(inout):: pint
1611  integer, intent(out):: ks
1612 ! Local
1613  real, parameter:: p00 = 1.e5
1614  real, dimension(km+1):: ze, pe1, peln, eta
1615  real, dimension(km):: dz, s_fac, dlnp
1616  real ztop, t0, dz0, sum1, tmp1
1617  real ep, es, alpha, beta, gama
1618  integer k
1619 
1620  pe1(1) = ptop
1621  peln(1) = log(pe1(1))
1622  pe1(km+1) = p00
1623  peln(km+1) = log(pe1(km+1))
1624 
1625  t0 = 270.
1626  ztop = rdgas/grav*t0*(peln(km+1) - peln(1))
1627 
1628  s_fac(km ) = 0.10
1629  s_fac(km-1) = 0.15
1630  s_fac(km-2) = 0.20
1631  s_fac(km-3) = 0.25
1632  s_fac(km-4) = 0.30
1633  s_fac(km-5) = 0.35
1634  s_fac(km-6) = 0.40
1635  s_fac(km-7) = 0.45
1636  s_fac(km-8) = 0.50
1637  s_fac(km-9) = 0.55
1638  s_fac(km-10) = 0.60
1639  s_fac(km-11) = 0.70
1640  s_fac(km-12) = 0.85
1641  s_fac(km-13) = 1.00
1642 
1643  do k=km-14, 9, -1
1644  s_fac(k) = min(10.0, s_rate * s_fac(k+1) )
1645  enddo
1646 
1647  s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9)
1648  s_fac(7) = 1.1 *s_fac(8)
1649  s_fac(6) = 1.15*s_fac(7)
1650  s_fac(5) = 1.2 *s_fac(6)
1651  s_fac(4) = 1.3 *s_fac(5)
1652  s_fac(3) = 1.4 *s_fac(4)
1653  s_fac(2) = 1.5 *s_fac(3)
1654  s_fac(1) = 1.6 *s_fac(2)
1655 
1656  sum1 = 0.
1657  do k=1,km
1658  sum1 = sum1 + s_fac(k)
1659  enddo
1660 
1661  dz0 = ztop / sum1
1662 
1663  do k=1,km
1664  dz(k) = s_fac(k) * dz0
1665  enddo
1666 
1667  ze(km+1) = 0.
1668  do k=km,1,-1
1669  ze(k) = ze(k+1) + dz(k)
1670  enddo
1671 
1672 ! Re-scale dz with the stretched ztop
1673  do k=1,km
1674  dz(k) = dz(k) * (ztop/ze(1))
1675  enddo
1676 
1677  do k=km,1,-1
1678  ze(k) = ze(k+1) + dz(k)
1679  enddo
1680 ! ze(1) = ztop
1681 
1682  if ( is_master() ) write(*,*) 'var55_dz: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1)
1683  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
1684 
1685 ! Given z --> p
1686  do k=1,km
1687  dz(k) = ze(k) - ze(k+1)
1688  dlnp(k) = grav*dz(k) / (rdgas*t0)
1689  enddo
1690  do k=2,km
1691  peln(k) = peln(k-1) + dlnp(k-1)
1692  pe1(k) = exp(peln(k))
1693  enddo
1694 
1695 ! Pe(k) = ak(k) + bk(k) * PS
1696 ! Locate pint and KS
1697  ks = 0
1698  do k=2,km
1699  if ( pint < pe1(k)) then
1700  ks = k-1
1701  exit
1702  endif
1703  enddo
1704  if ( is_master() ) then
1705  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
1706  endif
1707  pint = pe1(ks+1)
1708 
1709 #ifdef NO_UKMO_HB
1710  do k=1,ks+1
1711  ak(k) = pe1(k)
1712  bk(k) = 0.
1713  enddo
1714 
1715  do k=ks+2,km+1
1716  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
1717  ak(k) = pe1(k) - bk(k) * pe1(km+1)
1718  enddo
1719  bk(km+1) = 1.
1720  ak(km+1) = 0.
1721 #else
1722 ! Problematic for non-hydrostatic
1723  do k=1,km+1
1724  eta(k) = pe1(k) / pe1(km+1)
1725  enddo
1726 
1727  ep = eta(ks+1)
1728  es = eta(km)
1729 ! es = 1.
1730  alpha = (ep**2-2.*ep*es) / (es-ep)**2
1731  beta = 2.*ep*es**2 / (es-ep)**2
1732  gama = -(ep*es)**2 / (es-ep)**2
1733 
1734 ! Pure pressure:
1735  do k=1,ks+1
1736  ak(k) = eta(k)*1.e5
1737  bk(k) = 0.
1738  enddo
1739 
1740  do k=ks+2, km
1741  ak(k) = alpha*eta(k) + beta + gama/eta(k)
1742  ak(k) = ak(k)*1.e5
1743  enddo
1744  ak(km+1) = 0.
1745 
1746  do k=ks+2, km
1747  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
1748  enddo
1749  bk(km+1) = 1.
1750 #endif
1751 
1752  if ( is_master() ) then
1753  write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100.
1754  do k=1,km
1755  write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k)
1756  enddo
1757  tmp1 = ak(ks+1)
1758  do k=ks+1,km
1759  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
1760  enddo
1761  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
1762  endif
1763 
1764 
1765  end subroutine var55_dz
1766 
1767  subroutine hybrid_z_dz(km, dz, ztop, s_rate)
1768  integer, intent(in):: km
1769  real, intent(in):: s_rate
1770  real, intent(in):: ztop
1771  real, intent(out):: dz(km)
1772 ! Local
1773  real, parameter:: p00 = 1.e5
1774  real, dimension(km+1):: ze, pe1, peln, eta
1775  real, dimension(km):: s_fac, dlnp
1776  real t0, dz0, sum1, tmp1
1777  real ep, es, alpha, beta, gama
1778  integer k
1779 
1780  s_fac(km ) = 0.12
1781  s_fac(km-1) = 0.20
1782  s_fac(km-2) = 0.30
1783  s_fac(km-3) = 0.40
1784  s_fac(km-4) = 0.50
1785  s_fac(km-5) = 0.60
1786  s_fac(km-6) = 0.70
1787  s_fac(km-7) = 0.80
1788  s_fac(km-8) = 0.90
1789  s_fac(km-9) = 1.
1790 
1791  do k=km-10, 9, -1
1792  s_fac(k) = min(4.0, s_rate * s_fac(k+1) )
1793  enddo
1794 
1795  s_fac(8) = 1.05*s_fac(9)
1796  s_fac(7) = 1.1 *s_fac(8)
1797  s_fac(6) = 1.15*s_fac(7)
1798  s_fac(5) = 1.2 *s_fac(6)
1799  s_fac(4) = 1.3 *s_fac(5)
1800  s_fac(3) = 1.4 *s_fac(4)
1801  s_fac(2) = 1.5 *s_fac(3)
1802  s_fac(1) = 1.6 *s_fac(2)
1803 
1804  sum1 = 0.
1805  do k=1,km
1806  sum1 = sum1 + s_fac(k)
1807  enddo
1808 
1809  dz0 = ztop / sum1
1810 
1811  do k=1,km
1812  dz(k) = s_fac(k) * dz0
1813  enddo
1814 
1815  ze(km+1) = 0.
1816  do k=km,1,-1
1817  ze(k) = ze(k+1) + dz(k)
1818  enddo
1819 
1820  ze(1) = ztop
1821 
1822  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
1823 
1824  do k=1,km
1825  dz(k) = ze(k) - ze(k+1)
1826  enddo
1827 
1828  end subroutine hybrid_z_dz
1829 
1832  subroutine get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
1833  integer, intent(in) :: npz
1834  real, intent(in) :: p_s
1835  real, intent(in) :: ak(npz+1)
1836  real, intent(in) :: bk(npz+1)
1837  real, intent(in), optional :: pscale
1838  real, intent(out) :: pf(npz)
1839  real, intent(out) :: ph(npz+1)
1840  integer k
1841 
1842  ph(1) = ak(1)
1843  do k=2,npz+1
1844  ph(k) = ak(k) + bk(k)*p_s
1845  enddo
1846 
1847  if ( present(pscale) ) then
1848  do k=1,npz+1
1849  ph(k) = pscale*ph(k)
1850  enddo
1851  endif
1852 
1853  if( ak(1) > 1.e-8 ) then
1854  pf(1) = (ph(2) - ph(1)) / log(ph(2)/ph(1))
1855  else
1856  pf(1) = (ph(2) - ph(1)) * kappa/(kappa+1.)
1857  endif
1858 
1859  do k=2,npz
1860  pf(k) = (ph(k+1) - ph(k)) / log(ph(k+1)/ph(k))
1861  enddo
1862 
1863  end subroutine get_eta_level
1864 
1865 
1866 
1867  subroutine compute_dz(km, ztop, dz)
1868 
1869  integer, intent(in):: km
1870  real, intent(in):: ztop ! try 50.E3
1871  real, intent(out):: dz(km)
1872 !------------------------------
1873  real ze(km+1), dzt(km)
1874  integer k
1875 
1876 
1877 ! ztop = 30.E3
1878  dz(1) = ztop / real(km)
1879  dz(km) = 0.5*dz(1)
1880 
1881  do k=2,km-1
1882  dz(k) = dz(1)
1883  enddo
1884 
1885 ! Top:
1886  dz(1) = 2.*dz(2)
1887 
1888  ze(km+1) = 0.
1889  do k=km,1,-1
1890  ze(k) = ze(k+1) + dz(k)
1891  enddo
1892 
1893  if ( is_master() ) then
1894  write(*,*) 'Hybrid_z: dz, zm'
1895  do k=1,km
1896  dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000.
1897  write(*,*) k, dz(k), dzt(k)
1898  enddo
1899  endif
1900 
1901  end subroutine compute_dz
1902 
1903  subroutine compute_dz_var(km, ztop, dz)
1904 
1905  integer, intent(in):: km
1906  real, intent(in):: ztop ! try 50.E3
1907  real, intent(out):: dz(km)
1908 !------------------------------
1909  real, parameter:: s_rate = 1.0
1910  real ze(km+1)
1911  real s_fac(km)
1912  real sum1, dz0
1913  integer k
1914 
1915  s_fac(km ) = 0.125
1916  s_fac(km-1) = 0.20
1917  s_fac(km-2) = 0.30
1918  s_fac(km-3) = 0.40
1919  s_fac(km-4) = 0.50
1920  s_fac(km-5) = 0.60
1921  s_fac(km-6) = 0.70
1922  s_fac(km-7) = 0.80
1923  s_fac(km-8) = 0.90
1924  s_fac(km-9) = 1.
1925 
1926  do k=km-10, 9, -1
1927  s_fac(k) = s_rate * s_fac(k+1)
1928  enddo
1929 
1930  s_fac(8) = 1.05*s_fac(9)
1931  s_fac(7) = 1.1 *s_fac(8)
1932  s_fac(6) = 1.15*s_fac(7)
1933  s_fac(5) = 1.2 *s_fac(6)
1934  s_fac(4) = 1.3 *s_fac(5)
1935  s_fac(3) = 1.4 *s_fac(4)
1936  s_fac(2) = 1.5 *s_fac(3)
1937  s_fac(1) = 1.6 *s_fac(2)
1938 
1939  sum1 = 0.
1940  do k=1,km
1941  sum1 = sum1 + s_fac(k)
1942  enddo
1943 
1944  dz0 = ztop / sum1
1945 
1946  do k=1,km
1947  dz(k) = s_fac(k) * dz0
1948  enddo
1949 
1950  ze(1) = ztop
1951  ze(km+1) = 0.
1952  do k=km,2,-1
1953  ze(k) = ze(k+1) + dz(k)
1954  enddo
1955 
1956 ! Re-scale dz with the stretched ztop
1957  do k=1,km
1958  dz(k) = dz(k) * (ztop/ze(1))
1959  enddo
1960 
1961  do k=km,2,-1
1962  ze(k) = ze(k+1) + dz(k)
1963  enddo
1964 
1965  call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2)
1966 
1967  do k=1,km
1968  dz(k) = ze(k) - ze(k+1)
1969  enddo
1970 
1971  end subroutine compute_dz_var
1972 
1973  subroutine compute_dz_l32(km, ztop, dz)
1974 
1975  integer, intent(in):: km
1976  real, intent(out):: dz(km)
1977  real, intent(out):: ztop ! try 50.E3
1978 !------------------------------
1979  real dzt(km)
1980  real ze(km+1)
1981  real dz0, dz1, dz2
1982  real z0, z1, z2
1983  integer k, k0, k1, k2, n
1984 
1985 !-------------------
1986  k2 = 8
1987  z2 = 30.e3
1988 !-------------------
1989  k1 = 21
1990  z1 = 10.0e3
1991 !-------------------
1992  k0 = 2
1993  z0 = 0.
1994  dz0 = 75. ! meters
1995 !-------------------
1996 ! Treat the surface layer as a special layer
1997  ze(1) = z0
1998  dz(1) = dz0
1999 
2000  ze(2) = dz(1)
2001  dz0 = 1.5*dz0
2002  dz(2) = dz0
2003 
2004  ze(3) = ze(2) + dz(2)
2005 
2006  dz1 = 2.*(z1-ze(3) - k1*dz0) / (k1*(k1-1))
2007 
2008  do k=k0+1,k0+k1
2009  dz(k) = dz0 + (k-k0)*dz1
2010  ze(k+1) = ze(k) + dz(k)
2011  enddo
2012 
2013  dz0 = dz(k1+k0)
2014  dz2 = 2.*(z2-ze(k0+k1+1)-k2*dz0) / (k2*(k2-1))
2015 
2016  do k=k0+k1+1,k0+k1+k2
2017  dz(k) = dz0 + (k-k0-k1)*dz2
2018  ze(k+1) = ze(k) + dz(k)
2019  enddo
2020 
2021  dz(km) = 2.*dz(km-1)
2022  ztop = ze(km) + dz(km)
2023  ze(km+1) = ze(km) + dz(km)
2024 
2025  call zflip (dz, 1, km)
2026 
2027  ze(km+1) = 0.
2028  do k=km,1,-1
2029  ze(k) = ze(k+1) + dz(k)
2030  enddo
2031 
2032 ! if ( is_master() ) then
2033 ! write(*,*) 'Hybrid_z: dz, zm'
2034 ! do k=1,km
2035 ! dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000.
2036 ! write(*,*) k, dz(k), dzt(k)
2037 ! enddo
2038 ! endif
2039 
2040  end subroutine compute_dz_l32
2041 
2042  subroutine compute_dz_l101(km, ztop, dz)
2043 
2044  integer, intent(in):: km ! km==101
2045  real, intent(out):: dz(km)
2046  real, intent(out):: ztop ! try 30.E3
2047 !------------------------------
2048  real ze(km+1)
2049  real dz0, dz1
2050  real:: stretch_f = 1.16
2051  integer k, k0, k1
2052 
2053  k1 = 2
2054  k0 = 25
2055  dz0 = 40. ! meters
2056 
2057  ze(km+1) = 0.
2058 
2059  do k=km, k0, -1
2060  dz(k) = dz0
2061  ze(k) = ze(k+1) + dz(k)
2062  enddo
2063 
2064  do k=k0+1, k1, -1
2065  dz(k) = stretch_f * dz(k+1)
2066  ze(k) = ze(k+1) + dz(k)
2067  enddo
2068 
2069  dz(1) = 4.0*dz(2)
2070  ze(1) = ze(2) + dz(1)
2071  ztop = ze(1)
2072 
2073  if ( is_master() ) then
2074  write(*,*) 'Hybrid_z: dz, ze'
2075  do k=1,km
2076  write(*,*) k, 0.001*dz(k), 0.001*ze(k)
2077  enddo
2078 ! ztop (km) = 20.2859154
2079  write(*,*) 'ztop (km) =', ztop * 0.001
2080  endif
2081 
2082  end subroutine compute_dz_l101
2083 
2084  subroutine set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
2085 
2086  integer, intent(in):: is, ie, js, je, ng, km
2087  real, intent(in):: rgrav, ztop
2088  real, intent(in):: dz(km)
2089  real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng)
2090  real, intent(inout):: ze(is:ie,js:je,km+1)
2091  real, optional, intent(out):: dz3(is-ng:ie+ng,js-ng:je+ng,km)
2092 ! local
2093  logical:: filter_xy = .false.
2094  real, allocatable:: delz(:,:,:)
2095  integer ntimes
2096  real zint
2097  real:: z1(is:ie,js:je)
2098  real:: z(km+1)
2099  real sig, z_rat
2100  integer ks(is:ie,js:je)
2101  integer i, j, k, ks_min, kint
2102 
2103  z(km+1) = 0.
2104  do k=km,1,-1
2105  z(k) = z(k+1) + dz(k)
2106  enddo
2107 
2108  do j=js,je
2109  do i=is,ie
2110  ze(i,j, 1) = ztop
2111  ze(i,j,km+1) = hs(i,j) * rgrav
2112  enddo
2113  enddo
2114 
2115  do k=2,km
2116  do j=js,je
2117  do i=is,ie
2118  ze(i,j,k) = z(k)
2119  enddo
2120  enddo
2121  enddo
2122 
2123 ! Set interface:
2124 #ifndef USE_VAR_ZINT
2125  zint = 12.0e3
2126  ntimes = 2
2127  kint = 2
2128  do k=2,km
2129  if ( z(k)<=zint ) then
2130  kint = k
2131  exit
2132  endif
2133  enddo
2134 
2135  if ( is_master() ) write(*,*) 'Z_coord interface set at k=',kint, ' ZE=', z(kint)
2136 
2137  do j=js,je
2138  do i=is,ie
2139  z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1))
2140  do k=km,kint+1,-1
2141  ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat
2142  enddo
2143 !--------------------------------------
2144 ! Apply vertical smoother locally to dz
2145 !--------------------------------------
2146  call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2147  enddo
2148  enddo
2149 #else
2150 ! ZINT is a function of local terrain
2151  ntimes = 4
2152  do j=js,je
2153  do i=is,ie
2154  z1(i,j) = dim(ze(i,j,km+1), 2500.) + 7500.
2155  enddo
2156  enddo
2157 
2158  ks_min = km
2159  do j=js,je
2160  do i=is,ie
2161  do k=km,2,-1
2162  if ( z(k)>=z1(i,j) ) then
2163  ks(i,j) = k
2164  ks_min = min(ks_min, k)
2165  go to 555
2166  endif
2167  enddo
2168 555 continue
2169  enddo
2170  enddo
2171 
2172  do j=js,je
2173  do i=is,ie
2174  kint = ks(i,j) + 1
2175  z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1))
2176  do k=km,kint+1,-1
2177  ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat
2178  enddo
2179 !--------------------------------------
2180 ! Apply vertical smoother locally to dz
2181 !--------------------------------------
2182  call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2183  enddo
2184  enddo
2185 #endif
2186 
2187 #ifdef DEV_ETA
2188  if ( filter_xy ) then
2189  allocate (delz(isd:ied, jsd:jed, km) )
2190  ntimes = 2
2191  do k=1,km
2192  do j=js,je
2193  do i=is,ie
2194  delz(i,j,k) = ze(i,j,k+1) - ze(i,j,k)
2195  enddo
2196  enddo
2197  enddo
2198  call del2_cubed(delz, 0.2*da_min, npx, npy, km, ntimes)
2199  do k=km,2,-1
2200  do j=js,je
2201  do i=is,ie
2202  ze(i,j,k) = ze(i,j,k+1) - delz(i,j,k)
2203  enddo
2204  enddo
2205  enddo
2206  deallocate ( delz )
2207  endif
2208 #endif
2209  if ( present(dz3) ) then
2210  do k=1,km
2211  do j=js,je
2212  do i=is,ie
2213  dz3(i,j,k) = ze(i,j,k+1) - ze(i,j,k)
2214  enddo
2215  enddo
2216  enddo
2217  endif
2218 
2219  end subroutine set_hybrid_z
2220 
2221 
2222  subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
2223  integer, intent(in):: is, ie, js, je, km
2224  integer, intent(in):: ntimes, i, j
2225  real, intent(inout):: ze(is:ie,js:je,km+1)
2226 ! local:
2227  real, parameter:: df = 0.25
2228  real dz(km)
2229  real flux(km+1)
2230  integer k, n, k1, k2
2231 
2232  k2 = km-1
2233  do k=1,km
2234  dz(k) = ze(i,j,k+1) - ze(i,j,k)
2235  enddo
2236 
2237  do n=1,ntimes
2238  k1 = 2 + (ntimes-n)
2239 
2240  flux(k1 ) = 0.
2241  flux(k2+1) = 0.
2242  do k=k1+1,k2
2243  flux(k) = df*(dz(k) - dz(k-1))
2244  enddo
2245 
2246  do k=k1,k2
2247  dz(k) = dz(k) - flux(k) + flux(k+1)
2248  enddo
2249  enddo
2250 
2251  do k=km,1,-1
2252  ze(i,j,k) = ze(i,j,k+1) - dz(k)
2253  enddo
2254 
2255  end subroutine sm1_edge
2256 
2257 
2258 
2259  subroutine gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
2260  integer, intent(in):: km
2261  real, intent(in):: p0, ztop
2262  real, intent(inout):: ptop
2263  real, intent(inout):: ak(km+1), bk(km+1)
2264  real, intent(out):: pt1(km)
2265 ! Local
2266  logical:: isothermal
2267  real, dimension(km+1):: ze, pe1, pk1
2268  real, dimension(km):: dz1
2269  real t0, n2, s0
2270  integer k
2271 
2272 ! Set up vertical coordinare with constant del-z spacing:
2273  isothermal = .false.
2274  t0 = 300.
2275 
2276  if ( isothermal ) then
2277  n2 = grav**2/(cp_air*t0)
2278  else
2279  n2 = 0.0001
2280  endif
2281 
2282  s0 = grav*grav / (cp_air*n2)
2283 
2284  ze(km+1) = 0.
2285  do k=km,1,-1
2286  dz1(k) = ztop / real(km)
2287  ze(k) = ze(k+1) + dz1(k)
2288  enddo
2289 
2290 ! Given z --> p
2291  do k=1,km+1
2292  pe1(k) = p0*( (1.-s0/t0) + s0/t0*exp(-n2*ze(k)/grav) )**(1./kappa)
2293  enddo
2294 
2295  ptop = pe1(1)
2296 ! if ( is_master() ) write(*,*) 'GW_1D: computed model top (pa)=', ptop
2297 
2298 ! Set up "sigma" coordinate
2299  ak(1) = pe1(1)
2300  bk(1) = 0.
2301  do k=2,km
2302  bk(k) = (pe1(k) - pe1(1)) / (pe1(km+1)-pe1(1)) ! bk == sigma
2303  ak(k) = pe1(1)*(1.-bk(k))
2304  enddo
2305  ak(km+1) = 0.
2306  bk(km+1) = 1.
2307 
2308  do k=1,km+1
2309  pk1(k) = pe1(k) ** kappa
2310  enddo
2311 
2312 ! Compute volume mean potential temperature with hydrostatic eqn:
2313  do k=1,km
2314  pt1(k) = grav*dz1(k) / ( cp_air*(pk1(k+1)-pk1(k)) )
2315  enddo
2316 
2317  end subroutine gw_1d
2318 
2319  subroutine mount_waves(km, ak, bk, ptop, ks, pint)
2320  integer, intent(in):: km
2321  real, intent(out):: ak(km+1), bk(km+1)
2322  real, intent(out):: ptop, pint
2323  integer, intent(out):: ks
2324 ! Local
2325  real, parameter:: p00 = 1.e5
2326  real, dimension(km+1):: ze, pe1, peln, eta
2327  real, dimension(km):: dz, dlnp
2328  real ztop, t0, dz0, sum1, tmp1
2329  real ep, es, alpha, beta, gama, s_fac
2330  integer k, k500
2331 
2332  pint = 300.e2
2333 ! s_fac = 1.05
2334 ! dz0 = 500.
2335  if ( km <= 60 ) then
2336  s_fac = 1.0
2337  dz0 = 500.
2338  else
2339  s_fac = 1.
2340  dz0 = 250.
2341  endif
2342 
2343 ! Basic parameters for HIWPP mountain waves
2344  t0 = 300.
2345 ! ztop = 20.0e3; 500-m resolution in halft of the vertical domain
2346 ! ztop = real(km-1)*500.
2347 !-----------------------
2348 ! Compute temp ptop based on isothermal atm
2349 ! ptop = p00*exp(-grav*ztop/(rdgas*t0))
2350 
2351 ! Lowest half has constant resolution
2352  ze(km+1) = 0.
2353  do k=km, km-19, -1
2354  ze(k) = ze(k+1) + dz0
2355  enddo
2356 
2357 ! Stretching from 10-km and up:
2358  do k=km-20, 3, -1
2359  dz0 = s_fac * dz0
2360  ze(k) = ze(k+1) + dz0
2361  enddo
2362  ze(2) = ze(3) + sqrt(2.)*dz0
2363  ze(1) = ze(2) + 2.0*dz0
2364 
2365 ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
2366 
2367 ! Given z --> p
2368  do k=1,km
2369  dz(k) = ze(k) - ze(k+1)
2370  dlnp(k) = grav*dz(k) / (rdgas*t0)
2371  enddo
2372 
2373  pe1(km+1) = p00
2374  peln(km+1) = log(p00)
2375  do k=km,1,-1
2376  peln(k) = peln(k+1) - dlnp(k)
2377  pe1(k) = exp(peln(k))
2378  enddo
2379 
2380 ! Comnpute new ptop
2381  ptop = pe1(1)
2382 
2383 ! Pe(k) = ak(k) + bk(k) * PS
2384 ! Locate pint and KS
2385  ks = 0
2386  do k=2,km
2387  if ( pint < pe1(k)) then
2388  ks = k-1
2389  exit
2390  endif
2391  enddo
2392 
2393  if ( is_master() ) then
2394  write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1)
2395  write(*,*) 'Modified ptop =', ptop, ' ztop=', ze(1)/1000.
2396  do k=1,km
2397  write(*,*) k, 'ze =', ze(k)/1000.
2398  enddo
2399  endif
2400  pint = pe1(ks+1)
2401 
2402 #ifdef NO_UKMO_HB
2403  do k=1,ks+1
2404  ak(k) = pe1(k)
2405  bk(k) = 0.
2406  enddo
2407 
2408  do k=ks+2,km+1
2409  bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma
2410  ak(k) = pe1(k) - bk(k) * pe1(km+1)
2411  enddo
2412  bk(km+1) = 1.
2413  ak(km+1) = 0.
2414 #else
2415 ! Problematic for non-hydrostatic
2416  do k=1,km+1
2417  eta(k) = pe1(k) / pe1(km+1)
2418  enddo
2419  ep = eta(ks+1)
2420  es = eta(km)
2421 ! es = 1.
2422  alpha = (ep**2-2.*ep*es) / (es-ep)**2
2423  beta = 2.*ep*es**2 / (es-ep)**2
2424  gama = -(ep*es)**2 / (es-ep)**2
2425 
2426 ! Pure pressure:
2427  do k=1,ks+1
2428  ak(k) = eta(k)*1.e5
2429  bk(k) = 0.
2430  enddo
2431 
2432  do k=ks+2, km
2433  ak(k) = alpha*eta(k) + beta + gama/eta(k)
2434  ak(k) = ak(k)*1.e5
2435  enddo
2436  ak(km+1) = 0.
2437 
2438  do k=ks+2, km
2439  bk(k) = (pe1(k) - ak(k))/pe1(km+1)
2440  enddo
2441  bk(km+1) = 1.
2442 #endif
2443 
2444  if ( is_master() ) then
2445  tmp1 = ak(ks+1)
2446  do k=ks+1,km
2447  tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.e-5, (bk(k+1)-bk(k))) )
2448  enddo
2449  write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
2450  endif
2451 
2452  end subroutine mount_waves
2453 
2454 
2455  subroutine zflip(q, im, km)
2456  integer, intent(in):: im, km
2457  real, intent(inout):: q(im,km)
2458 !---
2459  integer i, k
2460  real qtmp
2461 
2462  do i = 1, im
2463  do k = 1, (km+1)/2
2464  qtmp = q(i,k)
2465  q(i,k) = q(i,km+1-k)
2466  q(i,km+1-k) = qtmp
2467  end do
2468  end do
2469 
2470  end subroutine zflip
2471 
2472 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:2223
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:1768
fv_eta_mod::var_hi2
subroutine var_hi2(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:1286
fv_eta_mod::var_les
subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:784
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:1833
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:2085
fv_eta_mod::gw_1d
subroutine, public gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
Definition: fv_eta.F90:2260
fv_eta_mod::var55_dz
subroutine var55_dz(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:1606
fv_eta_mod::compute_dz
subroutine, public compute_dz(km, ztop, dz)
Definition: fv_eta.F90:1868
fv_eta_mod::compute_dz_l32
subroutine, public compute_dz_l32(km, ztop, dz)
Definition: fv_eta.F90:1974
fv_eta_mod::zflip
subroutine zflip(q, im, km)
Definition: fv_eta.F90:2456
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:762
fv_eta_mod::set_eta
subroutine, public set_eta(km, ks, ptop, ak, bk, npz_type)
Definition: fv_eta.F90:288
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:1904
fv_eta_mod::var_gfs
subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:946
fv_eta_mod::var_dz
subroutine var_dz(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:1445
fv_eta_mod::mount_waves
subroutine mount_waves(km, ak, bk, ptop, ks, pint)
Definition: fv_eta.F90:2320
fv_eta_mod::var_hi
subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate)
Definition: fv_eta.F90:1110
fv_eta_mod::compute_dz_l101
subroutine, public compute_dz_l101(km, ztop, dz)
Definition: fv_eta.F90:2043