!*********************************************************************** !* GNU Lesser General Public License !* !* This file is part of the FV3 dynamical core. !* !* The FV3 dynamical core is free software: you can redistribute it !* and/or modify it under the terms of the !* GNU Lesser General Public License as published by the !* Free Software Foundation, either version 3 of the License, or !* (at your option) any later version. !* !* The FV3 dynamical core is distributed in the hope that it will be !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. !* See the GNU General Public License for more details. !* !* You should have received a copy of the GNU Lesser General Public !* License along with the FV3 dynamical core. !* If not, see . !*********************************************************************** !>@brief The module 'fv_eta' contains routine to set up the reference !! (Eulerian) pressure coordinate module fv_eta_mod ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
Module NameFunctions Included
constants_modkappa, grav, cp_air, rdgas
fv_mp_modis_master
mpp_modmpp_error, FATAL
use constants_mod, only: kappa, grav, cp_air, rdgas use fv_mp_mod, only: is_master use mpp_mod, only: FATAL, mpp_error implicit none private public set_eta, set_external_eta, get_eta_level, compute_dz_var, & compute_dz_L32, compute_dz_L101, set_hybrid_z, compute_dz, & gw_1d, sm1_edge, hybrid_z_dz contains !!!NOTE: USE_VAR_ETA not used in fvGFS #ifdef USE_VAR_ETA subroutine set_eta(km, ks, ptop, ak, bk) ! This is the easy to use version of the set_eta integer, intent(in):: km ! vertical dimension integer, intent(out):: ks ! number of pure p layers real:: a60(61),b60(61) ! Thfollowing L63 setting is the same as NCEP GFS's L64 except the top ! 3 layers data a60/300.0000, 430.00000, 558.00000, & 700.00000, 863.05803, 1051.07995, & 1265.75194, 1510.71101, 1790.05098, & 2108.36604, 2470.78817, 2883.03811, & 3351.46002, 3883.05187, 4485.49315, & 5167.14603, 5937.04991, 6804.87379, & 7780.84698, 8875.64338, 10100.20534, & 11264.35673, 12190.64366, 12905.42546, & 13430.87867, 13785.88765, 13986.77987, & 14047.96335, 13982.46770, 13802.40331, & 13519.33841, 13144.59486, 12689.45608, & 12165.28766, 11583.57006, 10955.84778, & 10293.60402, 9608.08306, 8910.07678, & 8209.70131, 7516.18560, 6837.69250, & 6181.19473, 5552.39653, 4955.72632, & 4394.37629, 3870.38682, 3384.76586, & 2937.63489, 2528.37666, 2155.78385, & 1818.20722, 1513.68173, 1240.03585, & 994.99144, 776.23591, 581.48797, & 408.53400, 255.26520, 119.70243, 0. / data b60/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00201, 0.00792, 0.01755, & 0.03079, 0.04751, 0.06761, & 0.09097, 0.11746, 0.14690, & 0.17911, 0.21382, 0.25076, & 0.28960, 0.32994, 0.37140, & 0.41353, 0.45589, 0.49806, & 0.53961, 0.58015, 0.61935, & 0.65692, 0.69261, 0.72625, & 0.75773, 0.78698, 0.81398, & 0.83876, 0.86138, 0.88192, & 0.90050, 0.91722, 0.93223, & 0.94565, 0.95762, 0.96827, & 0.97771, 0.98608, 0.99347, 1./ real, intent(out):: ak(km+1) real, intent(out):: bk(km+1) real, intent(out):: ptop ! model top (Pa) real pint, stretch_fac integer k real :: s_rate = -1.0 ! dummy value to not use var_les pint = 100.E2 !- Notes --------------------------------- ! low-top: ptop = 100. ! ~45 km ! mid-top: ptop = 10. ! ~60 km ! hi -top: ptop = 1. ! ~80 km !----------------------------------------- select case (km) ! Optimal number = 8 * N -1 (for 8 openMP threads) ! 16 * M -1 (for 16 openMP threads) #ifdef HIWPP #ifdef SUPER_K case (20) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (24) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (30) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (40) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (50) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (60) ptop = 56.e2 pint = ptop stretch_fac = 1.03 case (80) ptop = 56.e2 pint = ptop stretch_fac = 1.03 #else case (30) ! For Baroclinic Instability Test ptop = 2.26e2 pint = 250.E2 stretch_fac = 1.03 case (40) ptop = 50.e2 ! For super cell test pint = 300.E2 stretch_fac = 1.03 case (50) ! Mountain waves? ptop = 30.e2 stretch_fac = 1.025 case (60) ! For Baroclinic Instability Test #ifdef GFSL60 ks = 20 ptop = a60(1) pint = a60(ks+1) do k=1,km+1 ak(k) = a60(k) bk(k) = b60(k) enddo #else ptop = 3.e2 ! pint = 250.E2 pint = 300.E2 ! revised for Moist test stretch_fac = 1.03 #endif #endif case (64) !!! ptop = 3.e2 ptop = 2.0e2 pint = 300.E2 stretch_fac = 1.03 #else ! *Very-low top: for idealized super-cell simulation: case (50) ptop = 50.e2 pint = 250.E2 stretch_fac = 1.03 case (60) ptop = 40.e2 pint = 250.E2 stretch_fac = 1.03 case (90) ! super-duper cell ptop = 40.e2 stretch_fac = 1.025 #endif ! Low-top: case (31) ! N = 4, M=2 ptop = 100. stretch_fac = 1.035 case (32) ! N = 4, M=2 ptop = 100. stretch_fac = 1.035 case (39) ! N = 5 ptop = 100. stretch_fac = 1.035 case (41) ptop = 100. stretch_fac = 1.035 case (47) ! N = 6, M=3 ptop = 100. stretch_fac = 1.035 case (51) ptop = 100. stretch_fac = 1.03 case (52) ! very low top ptop = 30.e2 ! for special DPM RCE experiments stretch_fac = 1.03 ! Mid-top: case (55) ! N = 7 ptop = 10. stretch_fac = 1.035 ! Hi-top: case (63) ! N = 8, M=4 ptop = 1. ! c360 or c384 stretch_fac = 1.035 case (71) ! N = 9 ptop = 1. stretch_fac = 1.03 case (79) ! N = 10, M=5 ptop = 1. stretch_fac = 1.03 case (127) ! N = 10, M=5 ptop = 1. stretch_fac = 1.03 case (151) ptop = 75.e2 pint = 500.E2 s_rate = 1.01 case default ptop = 1. stretch_fac = 1.03 end select #ifdef MOUNTAIN_WAVES call mount_waves(km, ak, bk, ptop, ks, pint) #else if (s_rate > 0.) then call var_les(km, ak, bk, ptop, ks, pint, s_rate) else if ( km > 79 ) then call var_hi2(km, ak, bk, ptop, ks, pint, stretch_fac) elseif (km==5 .or. km==10 ) then ! Equivalent Shallow Water: for NGGPS, variable-resolution testing ptop = 500.e2 ks = 0 do k=1,km+1 bk(k) = real(k-1) / real (km) ak(k) = ptop*(1.-bk(k)) enddo else #ifndef GFSL60 call var_hi(km, ak, bk, ptop, ks, pint, stretch_fac) #endif endif #endif endif ptop = ak(1) pint = ak(ks+1) end subroutine set_eta subroutine mount_waves(km, ak, bk, ptop, ks, pint) integer, intent(in):: km real, intent(out):: ak(km+1), bk(km+1) real, intent(out):: ptop, pint integer, intent(out):: ks ! Local real, parameter:: p00 = 1.E5 real, dimension(km+1):: ze, pe1, peln, eta real, dimension(km):: dz, dlnp real ztop, t0, dz0, sum1, tmp1 real ep, es, alpha, beta, gama, s_fac integer k, k500 pint = 300.e2 ! s_fac = 1.05 ! dz0 = 500. if ( km <= 60 ) then s_fac = 1.0 dz0 = 500. else s_fac = 1. dz0 = 250. endif ! Basic parameters for HIWPP mountain waves t0 = 300. ! ztop = 20.0e3; 500-m resolution in halft of the vertical domain ! ztop = real(km-1)*500. !----------------------- ! Compute temp ptop based on isothermal atm ! ptop = p00*exp(-grav*ztop/(rdgas*t0)) ! Lowest half has constant resolution ze(km+1) = 0. do k=km, km-19, -1 ze(k) = ze(k+1) + dz0 enddo ! Stretching from 10-km and up: do k=km-20, 3, -1 dz0 = s_fac * dz0 ze(k) = ze(k+1) + dz0 enddo ze(2) = ze(3) + sqrt(2.)*dz0 ze(1) = ze(2) + 2.0*dz0 ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1) ! Given z --> p do k=1,km dz(k) = ze(k) - ze(k+1) dlnp(k) = grav*dz(k) / (rdgas*t0) enddo pe1(km+1) = p00 peln(km+1) = log(p00) do k=km,1,-1 peln(k) = peln(k+1) - dlnp(k) pe1(k) = exp(peln(k)) enddo ! Comnpute new ptop ptop = pe1(1) ! Pe(k) = ak(k) + bk(k) * PS ! Locate pint and KS ks = 0 do k=2,km if ( pint < pe1(k)) then ks = k-1 exit endif enddo if ( is_master() ) then write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1) write(*,*) 'Modified ptop =', ptop, ' ztop=', ze(1)/1000. do k=1,km write(*,*) k, 'ze =', ze(k)/1000. enddo endif pint = pe1(ks+1) #ifdef NO_UKMO_HB do k=1,ks+1 ak(k) = pe1(k) bk(k) = 0. enddo do k=ks+2,km+1 bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma ak(k) = pe1(k) - bk(k) * pe1(km+1) enddo bk(km+1) = 1. ak(km+1) = 0. #else ! Problematic for non-hydrostatic do k=1,km+1 eta(k) = pe1(k) / pe1(km+1) enddo ep = eta(ks+1) es = eta(km) ! es = 1. alpha = (ep**2-2.*ep*es) / (es-ep)**2 beta = 2.*ep*es**2 / (es-ep)**2 gama = -(ep*es)**2 / (es-ep)**2 ! Pure pressure: do k=1,ks+1 ak(k) = eta(k)*1.e5 bk(k) = 0. enddo do k=ks+2, km ak(k) = alpha*eta(k) + beta + gama/eta(k) ak(k) = ak(k)*1.e5 enddo ak(km+1) = 0. do k=ks+2, km bk(k) = (pe1(k) - ak(k))/pe1(km+1) enddo bk(km+1) = 1. #endif if ( is_master() ) then tmp1 = ak(ks+1) do k=ks+1,km tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) ) enddo write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100. endif end subroutine mount_waves #else !>This is the version of set_eta used in fvGFS and AM4 !>@note 01/2018: 'set_eta' is being cleaned up. subroutine set_eta(km, ks, ptop, ak, bk) integer, intent(in):: km !< vertical dimension integer, intent(out):: ks !< number of pure p layers real, intent(out):: ak(km+1) real, intent(out):: bk(km+1) real, intent(out):: ptop !< model top (Pa) ! local real a24(25),b24(25) !< GFDL AM2L24 real a26(27),b26(27) !< Jablonowski & Williamson 26-level real a32(33),b32(33) real a32w(33),b32w(33) real a47(48),b47(48) real a48(49),b48(49) real a52(53),b52(53) real a54(55),b54(55) real a56(57),b56(57) real a60(61),b60(61) real a63(64),b63(64) real a64(65),b64(65) real a68(69),b68(69) !< cjg: grid with enhanced PBL resolution real a96(97),b96(97) !< cjg: grid with enhanced PBL resolution real a100(101),b100(101) real a104(105),b104(105) real a125(126),b125(126) real:: p0=1000.E2 real:: pc=200.E2 real pt, pint, lnpe, dlnp real press(km+1), pt1(km) integer k ! Definition: press(i,j,k) = ak(k) + bk(k) * ps(i,j) !----------------------------------------------- ! GFDL AM2-L24: modified by SJL at the model top !----------------------------------------------- ! data a24 / 100.0000, 1050.0000, 3474.7942, 7505.5556, 12787.2428, & data a24 / 100.0000, 903.4465, 3474.7942, 7505.5556, 12787.2428, & 19111.3683, 21854.9274, 22884.1866, 22776.3058, 21716.1604, & 20073.2963, 18110.5123, 16004.7832, 13877.6253, 11812.5452, & 9865.8840, 8073.9726, 6458.0834, 5027.9899, 3784.6085, & 2722.0086, 1828.9752, 1090.2396, 487.4595, 0.0000 / data b24 / 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, & 0.0000000, 0.0435679, 0.1102275, 0.1922249, 0.2817656, & 0.3694997, 0.4532348, 0.5316253, 0.6038733, 0.6695556, & 0.7285176, 0.7808017, 0.8265992, 0.8662148, 0.9000406, & 0.9285364, 0.9522140, 0.9716252, 0.9873523, 1.0000000 / ! Jablonowski & Williamson 26-level setup data a26 / 219.4067, 489.5209, 988.2418, 1805.2010, 2983.7240, 4462.3340, & 6160.5870, 7851.2430, 7731.2710, 7590.1310, 7424.0860, & 7228.7440, 6998.9330, 6728.5740, 6410.5090, 6036.3220, & 5596.1110, 5078.2250, 4468.9600, 3752.1910, 2908.9490, & 2084.739, 1334.443, 708.499, 252.1360, 0.0, 0.0 / data b26 / 0.0, 0.0, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000,& 0.0000000, 0.01505309, 0.03276228, 0.05359622, 0.07810627, & 0.1069411, 0.1408637, 0.1807720, 0.2277220, 0.2829562, & 0.3479364, 0.4243822, 0.5143168, 0.6201202, 0.7235355, & 0.8176768, 0.8962153, 0.9534761, 0.9851122, 1.0000000 / ! High-resolution troposphere setup #ifdef OLD_32 ! Revised Apr 14, 2004: PINT = 245.027 mb data a32/100.00000, 400.00000, 818.60211, & 1378.88653, 2091.79519, 2983.64084, & 4121.78960, 5579.22148, 7419.79300, & 9704.82578, 12496.33710, 15855.26306, & 19839.62499, 24502.73262, 28177.10152, & 29525.28447, 29016.34358, 27131.32792, & 24406.11225, 21326.04907, 18221.18357, & 15275.14642, 12581.67796, 10181.42843, & 8081.89816, 6270.86956, 4725.35001, & 3417.39199, 2317.75459, 1398.09473, & 632.49506, 0.00000, 0.00000 / data b32/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.01711, & 0.06479, 0.13730, 0.22693, & 0.32416, 0.42058, 0.51105, & 0.59325, 0.66628, 0.73011, & 0.78516, 0.83217, 0.87197, & 0.90546, 0.93349, 0.95685, & 0.97624, 0.99223, 1.00000 / #else ! SJL June 26, 2012 ! pint= 55.7922 data a32/100.00000, 400.00000, 818.60211, & 1378.88653, 2091.79519, 2983.64084, & 4121.78960, 5579.22148, 6907.19063, & 7735.78639, 8197.66476, 8377.95525, & 8331.69594, 8094.72213, 7690.85756, & 7139.01788, 6464.80251, 5712.35727, & 4940.05347, 4198.60465, 3516.63294, & 2905.19863, 2366.73733, 1899.19455, & 1497.78137, 1156.25252, 867.79199, & 625.59324, 423.21322, 254.76613, & 115.06646, 0.00000, 0.00000 / data b32/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00513, & 0.01969, 0.04299, 0.07477, & 0.11508, 0.16408, 0.22198, & 0.28865, 0.36281, 0.44112, & 0.51882, 0.59185, 0.65810, & 0.71694, 0.76843, 0.81293, & 0.85100, 0.88331, 0.91055, & 0.93338, 0.95244, 0.96828, & 0.98142, 0.99223, 1.00000 / #endif !--------------------- ! Wilson's 32L settings: !--------------------- ! Top changed to 0.01 mb data a32w/ 1.00, 26.6378, 84.5529, 228.8592, & 539.9597, 1131.7087, 2141.8082, 3712.0454, & 5963.5317, 8974.1873, 12764.5388, 17294.5911, & 20857.7007, 22221.8651, 22892.7202, 22891.1641, & 22286.0724, 21176.0846, 19673.0671, 17889.0989, & 15927.5060, 13877.6239, 11812.5474, 9865.8830, & 8073.9717, 6458.0824, 5027.9893, 3784.6104, & 2722.0093, 1828.9741, 1090.2397, 487.4575, & 0.0000 / data b32w/ 0.0000, 0.0000, 0.0000, 0.0000, & 0.0000, 0.0000, 0.0000, 0.0000, & 0.0000, 0.0000, 0.0000, 0.0000, & 0.0159, 0.0586, 0.1117, 0.1734, & 0.2415, 0.3137, 0.3878, 0.4619, & 0.5344, 0.6039, 0.6696, 0.7285, & 0.7808, 0.8266, 0.8662, 0.9000, & 0.9285, 0.9522, 0.9716, 0.9874, & 1.0000 / #ifdef OLD_L47 ! QBO setting with ptop = 0.1 mb and p_full=0.17 mb; pint ~ 100 mb data a47/ 10.00000, 24.45365, 48.76776, & 85.39458, 133.41983, 191.01402, & 257.94919, 336.63306, 431.52741, & 548.18995, 692.78825, 872.16512, & 1094.18467, 1368.11917, 1704.99489, & 2117.91945, 2622.42986, 3236.88281, & 3982.89623, 4885.84733, 5975.43260, & 7286.29500, 8858.72424, 10739.43477, & 12982.41110, 15649.68745, 18811.37629, & 22542.71275, 25724.93857, 27314.36781, & 27498.59474, 26501.79312, 24605.92991, & 22130.51655, 19381.30274, 16601.56419, & 13952.53231, 11522.93244, 9350.82303, & 7443.47723, 5790.77434, 4373.32696, & 3167.47008, 2148.51663, 1293.15510, & 581.62429, 0.00000, 0.00000 / data b47/ 0.0000, 0.0000, 0.0000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.01188, 0.04650, & 0.10170, 0.17401, 0.25832, & 0.34850, 0.43872, 0.52448, & 0.60307, 0.67328, 0.73492, & 0.78834, 0.83418, 0.87320, & 0.90622, 0.93399, 0.95723, & 0.97650, 0.99223, 1.00000 / #else ! Oct 23, 2012 ! QBO setting with ptop = 0.1 mb, pint ~ 60 mb data a47/ 10.00000, 24.45365, 48.76776, & 85.39458, 133.41983, 191.01402, & 257.94919, 336.63306, 431.52741, & 548.18995, 692.78825, 872.16512, & 1094.18467, 1368.11917, 1704.99489, & 2117.91945, 2622.42986, 3236.88281, & 3982.89623, 4885.84733, 5975.43260, & 7019.26669, 7796.15848, 8346.60209, & 8700.31838, 8878.27554, 8894.27179, & 8756.46404, 8469.60171, 8038.92687, & 7475.89006, 6803.68067, 6058.68992, & 5285.28859, 4526.01565, 3813.00206, & 3164.95553, 2589.26318, 2085.96929, & 1651.11596, 1278.81205, 962.38875, & 695.07046, 470.40784, 282.61654, & 126.92745, 0.00000, 0.00000 / data b47/ 0.0000, 0.0000, 0.0000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00267, 0.01063, 0.02393, & 0.04282, 0.06771, 0.09917, & 0.13786, 0.18444, 0.23925, & 0.30193, 0.37100, 0.44379, & 0.51695, 0.58727, 0.65236, & 0.71094, 0.76262, 0.80757, & 0.84626, 0.87930, 0.90731, & 0.93094, 0.95077, 0.96733, & 0.98105, 0.99223, 1.00000 / #endif data a48/ & 1.00000, 2.69722, 5.17136, & 8.89455, 14.24790, 22.07157, & 33.61283, 50.48096, 74.79993, & 109.40055, 158.00460, 225.44108, & 317.89560, 443.19350, 611.11558, & 833.74392, 1125.83405, 1505.20759, & 1993.15829, 2614.86254, 3399.78420, & 4382.06240, 5600.87014, 7100.73115, & 8931.78242, 11149.97021, 13817.16841, & 17001.20930, 20775.81856, 23967.33875, & 25527.64563, 25671.22552, 24609.29622, & 22640.51220, 20147.13482, 17477.63530, & 14859.86462, 12414.92533, 10201.44191, & 8241.50255, 6534.43202, 5066.17865, & 3815.60705, 2758.60264, 1870.64631, & 1128.33931, 510.47983, 0.00000, & 0.00000 / data b48/ & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.01253, & 0.04887, 0.10724, 0.18455, & 0.27461, 0.36914, 0.46103, & 0.54623, 0.62305, 0.69099, & 0.75016, 0.80110, 0.84453, & 0.88127, 0.91217, 0.93803, & 0.95958, 0.97747, 0.99223, & 1.00000 / ! High PBL resolution with top at 1 mb ! SJL modified May 7, 2013 to ptop ~ 100 mb data a54/100.00000, 254.83931, 729.54278, & 1602.41121, 2797.50667, 4100.18977, & 5334.87140, 6455.24153, 7511.80944, & 8580.26355, 9714.44293, 10938.62253, & 12080.36051, 12987.13921, 13692.75084, & 14224.92180, 14606.55444, 14856.69953, & 14991.32121, 15023.90075, 14965.91493, & 14827.21612, 14616.33505, 14340.72252, & 14006.94280, 13620.82849, 13187.60470, & 12711.98873, 12198.27003, 11650.37451, & 11071.91608, 10466.23819, 9836.44706, & 9185.43852, 8515.96231, 7831.01080, & 7135.14301, 6436.71659, 5749.00215, & 5087.67188, 4465.67510, 3889.86419, & 3361.63433, 2879.51065, 2441.02496, & 2043.41345, 1683.80513, 1359.31122, & 1067.09135, 804.40101, 568.62625, & 357.32525, 168.33263, 0.00000, & 0.00000 / data b54/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00180, 0.00694, 0.01510, & 0.02601, 0.03942, 0.05515, & 0.07302, 0.09288, 0.11459, & 0.13803, 0.16307, 0.18960, & 0.21753, 0.24675, 0.27716, & 0.30866, 0.34115, 0.37456, & 0.40879, 0.44375, 0.47935, & 0.51551, 0.55215, 0.58916, & 0.62636, 0.66334, 0.69946, & 0.73395, 0.76622, 0.79594, & 0.82309, 0.84780, 0.87020, & 0.89047, 0.90876, 0.92524, & 0.94006, 0.95336, 0.96529, & 0.97596, 0.98551, 0.99400, & 1.00000 / ! The 56-L setup data a56/ 10.00000, 24.97818, 58.01160, & 115.21466, 199.29210, 309.39897, & 445.31785, 610.54747, 812.28518, & 1059.80882, 1363.07092, 1732.09335, & 2176.91502, 2707.68972, 3334.70962, & 4068.31964, 4918.76594, 5896.01890, & 7009.59166, 8268.36324, 9680.41211, & 11252.86491, 12991.76409, 14901.95764, & 16987.01313, 19249.15733, 21689.24182, & 23845.11055, 25330.63353, 26243.52467, & 26663.84998, 26657.94696, 26281.61371, & 25583.05256, 24606.03265, 23393.39510, & 21990.28845, 20445.82122, 18811.93894, & 17139.59660, 15473.90350, 13850.50167, & 12294.49060, 10821.62655, 9440.57746, & 8155.11214, 6965.72496, 5870.70511, & 4866.83822, 3949.90019, 3115.03562, & 2357.07879, 1670.87329, 1051.65120, & 495.51399, 0.00000, 0.00000 / data b56 /0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00462, 0.01769, 0.03821, & 0.06534, 0.09834, 0.13659, & 0.17947, 0.22637, 0.27660, & 0.32929, 0.38343, 0.43791, & 0.49162, 0.54361, 0.59319, & 0.63989, 0.68348, 0.72391, & 0.76121, 0.79545, 0.82679, & 0.85537, 0.88135, 0.90493, & 0.92626, 0.94552, 0.96286, & 0.97840, 0.99223, 1.00000 / data a60/ 1.7861000000e-01, 1.0805100000e+00, 3.9647100000e+00, & 9.7516000000e+00, 1.9816580000e+01, 3.6695950000e+01, & 6.2550570000e+01, 9.9199620000e+01, 1.4792505000e+02, & 2.0947487000e+02, 2.8422571000e+02, 3.7241721000e+02, & 4.7437835000e+02, 5.9070236000e+02, 7.2236063000e+02, & 8.7076746000e+02, 1.0378138800e+03, 1.2258877300e+03, & 1.4378924600e+03, 1.6772726600e+03, 1.9480506400e+03, & 2.2548762700e+03, 2.6030909400e+03, 2.9988059200e+03, & 3.4489952300e+03, 3.9616028900e+03, 4.5456641600e+03, & 5.2114401700e+03, 5.9705644000e+03, 6.8361981800e+03, & 7.8231906000e+03, 8.9482351000e+03, 1.0230010660e+04, & 1.1689289750e+04, 1.3348986860e+04, 1.5234111060e+04, & 1.7371573230e+04, 1.9789784580e+04, 2.2005564550e+04, & 2.3550115120e+04, 2.4468583320e+04, 2.4800548800e+04, & 2.4582445070e+04, 2.3849999620e+04, 2.2640519740e+04, & 2.0994737150e+04, 1.8957848730e+04, 1.6579413230e+04, & 1.4080071030e+04, 1.1753630920e+04, 9.6516996300e+03, & 7.7938009300e+03, 6.1769062800e+03, 4.7874276000e+03, & 3.6050497500e+03, 2.6059860700e+03, 1.7668328200e+03, & 1.0656131200e+03, 4.8226201000e+02, 0.0000000000e+00, & 0.0000000000e+00 / data b60/ 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 5.0600000000e-03, & 2.0080000000e-02, 4.4900000000e-02, 7.9360000000e-02, & 1.2326000000e-01, 1.7634000000e-01, 2.3820000000e-01, & 3.0827000000e-01, 3.8581000000e-01, 4.6989000000e-01, & 5.5393000000e-01, 6.2958000000e-01, 6.9642000000e-01, & 7.5458000000e-01, 8.0463000000e-01, 8.4728000000e-01, & 8.8335000000e-01, 9.1368000000e-01, 9.3905000000e-01, & 9.6020000000e-01, 9.7775000000e-01, 9.9223000000e-01, & 1.0000000000e+00 / ! This is activated by USE_GFSL63 ! Thfollowing L63 setting is the same as NCEP GFS's L64 except the top ! 3 layers data a63/64.247, 137.790, 221.958, & 318.266, 428.434, 554.424, & 698.457, 863.05803, 1051.07995, & 1265.75194, 1510.71101, 1790.05098, & 2108.36604, 2470.78817, 2883.03811, & 3351.46002, 3883.05187, 4485.49315, & 5167.14603, 5937.04991, 6804.87379, & 7780.84698, 8875.64338, 10100.20534, & 11264.35673, 12190.64366, 12905.42546, & 13430.87867, 13785.88765, 13986.77987, & 14047.96335, 13982.46770, 13802.40331, & 13519.33841, 13144.59486, 12689.45608, & 12165.28766, 11583.57006, 10955.84778, & 10293.60402, 9608.08306, 8910.07678, & 8209.70131, 7516.18560, 6837.69250, & 6181.19473, 5552.39653, 4955.72632, & 4394.37629, 3870.38682, 3384.76586, & 2937.63489, 2528.37666, 2155.78385, & 1818.20722, 1513.68173, 1240.03585, & 994.99144, 776.23591, 581.48797, & 408.53400, 255.26520, 119.70243, 0. / data b63/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00201, 0.00792, 0.01755, & 0.03079, 0.04751, 0.06761, & 0.09097, 0.11746, 0.14690, & 0.17911, 0.21382, 0.25076, & 0.28960, 0.32994, 0.37140, & 0.41353, 0.45589, 0.49806, & 0.53961, 0.58015, 0.61935, & 0.65692, 0.69261, 0.72625, & 0.75773, 0.78698, 0.81398, & 0.83876, 0.86138, 0.88192, & 0.90050, 0.91722, 0.93223, & 0.94565, 0.95762, 0.96827, & 0.97771, 0.98608, 0.99347, 1./ #ifdef GFSL64 data a64/20.00000, 68.00000, 137.79000, & 221.95800, 318.26600, 428.43400, & 554.42400, 698.45700, 863.05803, & 1051.07995, 1265.75194, 1510.71101, & 1790.05098, 2108.36604, 2470.78817, & 2883.03811, 3351.46002, 3883.05187, & 4485.49315, 5167.14603, 5937.04991, & 6804.87379, 7780.84698, 8875.64338, & 9921.40745, 10760.99844, 11417.88354, & 11911.61193, 12258.61668, 12472.89642, & 12566.58298, 12550.43517, 12434.26075, & 12227.27484, 11938.39468, 11576.46910, & 11150.43640, 10669.41063, 10142.69482, & 9579.72458, 8989.94947, 8382.67090, & 7766.85063, 7150.91171, 6542.55077, & 5948.57894, 5374.81094, 4825.99383, & 4305.79754, 3816.84622, 3360.78848, & 2938.39801, 2549.69756, 2194.08449, & 1870.45732, 1577.34218, 1313.00028, & 1075.52114, 862.90778, 673.13815, & 504.22118, 354.22752, 221.32110, & 103.78014, 0./ data b64/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00179, 0.00705, 0.01564, & 0.02749, 0.04251, 0.06064, & 0.08182, 0.10595, 0.13294, & 0.16266, 0.19492, 0.22950, & 0.26615, 0.30455, 0.34435, & 0.38516, 0.42656, 0.46815, & 0.50949, 0.55020, 0.58989, & 0.62825, 0.66498, 0.69987, & 0.73275, 0.76351, 0.79208, & 0.81845, 0.84264, 0.86472, & 0.88478, 0.90290, 0.91923, & 0.93388, 0.94697, 0.95865, & 0.96904, 0.97826, 0.98642, & 0.99363, 1./ #else data a64/1.00000, 3.90000, 8.70000, & 15.42000, 24.00000, 34.50000, & 47.00000, 61.50000, 78.60000, & 99.13500, 124.12789, 154.63770, & 191.69700, 236.49300, 290.38000, & 354.91000, 431.82303, 523.09300, & 630.92800, 757.79000, 906.45000, & 1079.85000, 1281.00000, 1515.00000, & 1788.00000, 2105.00000, 2470.00000, & 2889.00000, 3362.00000, 3890.00000, & 4475.00000, 5120.00000, 5830.00000, & 6608.00000, 7461.00000, 8395.00000, & 9424.46289, 10574.46880, 11864.80270, & 13312.58890, 14937.03710, 16759.70700, & 18804.78710, 21099.41210, 23674.03710, & 26562.82810, 29804.11720, 32627.31640, & 34245.89840, 34722.28910, 34155.19920, & 32636.50390, 30241.08200, 27101.44920, & 23362.20700, 19317.05270, 15446.17090, & 12197.45210, 9496.39941, 7205.66992, & 5144.64307, 3240.79346, 1518.62134, & 0.00000, 0.00000 / data b64/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00813, & 0.03224, 0.07128, 0.12445, & 0.19063, 0.26929, 0.35799, & 0.45438, 0.55263, 0.64304, & 0.71703, 0.77754, 0.82827, & 0.87352, 0.91502, 0.95235, & 0.98511, 1.00000 / #endif !-->cjg data a68/1.00000, 2.68881, 5.15524, & 8.86683, 14.20349, 22.00278, & 33.50807, 50.32362, 74.56680, & 109.05958, 157.51214, 224.73844, & 316.90481, 441.81219, 609.21090, & 831.14537, 1122.32514, 1500.51628, & 1986.94617, 2606.71274, 3389.18802, & 4368.40473, 5583.41379, 7078.60015, & 8903.94455, 11115.21886, 13774.60566, & 16936.82070, 20340.47045, 23193.71492, & 24870.36141, 25444.59363, 25252.57081, & 24544.26211, 23474.29096, 22230.65331, & 20918.50731, 19589.96280, 18296.26682, & 17038.02866, 15866.85655, 14763.18943, & 13736.83624, 12794.11850, 11930.72442, & 11137.17217, 10404.78946, 9720.03954, & 9075.54055, 8466.72650, 7887.12346, & 7333.90490, 6805.43028, 6297.33773, & 5805.78227, 5327.94995, 4859.88765, & 4398.63854, 3942.81761, 3491.08449, & 3043.04531, 2598.71608, 2157.94527, & 1720.87444, 1287.52805, 858.02944, & 432.71276, 8.10905, 0.00000 / data b68/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00283, 0.01590, & 0.04412, 0.08487, 0.13284, & 0.18470, 0.23828, 0.29120, & 0.34211, 0.39029, 0.43518, & 0.47677, 0.51536, 0.55091, & 0.58331, 0.61263, 0.63917, & 0.66333, 0.68552, 0.70617, & 0.72555, 0.74383, 0.76117, & 0.77765, 0.79335, 0.80838, & 0.82287, 0.83693, 0.85069, & 0.86423, 0.87760, 0.89082, & 0.90392, 0.91689, 0.92973, & 0.94244, 0.95502, 0.96747, & 0.97979, 0.99200, 1.00000 / data a96/1.00000, 2.35408, 4.51347, & 7.76300, 12.43530, 19.26365, & 29.33665, 44.05883, 65.28397, & 95.48274, 137.90344, 196.76073, & 277.45330, 386.81095, 533.37018, & 727.67600, 982.60677, 1313.71685, & 1739.59104, 2282.20281, 2967.26766, & 3824.58158, 4888.33404, 6197.38450, & 7795.49158, 9731.48414, 11969.71024, & 14502.88894, 17304.52434, 20134.76139, & 22536.63814, 24252.54459, 25230.65591, & 25585.72044, 25539.91412, 25178.87141, & 24644.84493, 23978.98781, 23245.49366, & 22492.11600, 21709.93990, 20949.64473, & 20225.94258, 19513.31158, 18829.32485, & 18192.62250, 17589.39396, 17003.45386, & 16439.01774, 15903.91204, 15396.39758, & 14908.02140, 14430.65897, 13967.88643, & 13524.16667, 13098.30227, 12687.56457, & 12287.08757, 11894.41553, 11511.54106, & 11139.22483, 10776.01912, 10419.75711, & 10067.11881, 9716.63489, 9369.61967, & 9026.69066, 8687.29884, 8350.04978, & 8013.20925, 7677.12187, 7343.12994, & 7011.62844, 6681.98102, 6353.09764, & 6025.10535, 5699.10089, 5375.54503, & 5053.63074, 4732.62740, 4413.38037, & 4096.62775, 3781.79777, 3468.45371, & 3157.19882, 2848.25306, 2541.19150, & 2236.21942, 1933.50628, 1632.83741, & 1334.35954, 1038.16655, 744.22318, & 452.71094, 194.91899, 0.00000, & 0.00000 / data b96/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00193, & 0.00974, 0.02538, 0.04876, & 0.07817, 0.11081, 0.14514, & 0.18007, 0.21486, 0.24866, & 0.28088, 0.31158, 0.34030, & 0.36701, 0.39210, 0.41554, & 0.43733, 0.45774, 0.47707, & 0.49540, 0.51275, 0.52922, & 0.54495, 0.56007, 0.57459, & 0.58850, 0.60186, 0.61471, & 0.62715, 0.63922, 0.65095, & 0.66235, 0.67348, 0.68438, & 0.69510, 0.70570, 0.71616, & 0.72651, 0.73675, 0.74691, & 0.75700, 0.76704, 0.77701, & 0.78690, 0.79672, 0.80649, & 0.81620, 0.82585, 0.83542, & 0.84492, 0.85437, 0.86375, & 0.87305, 0.88229, 0.89146, & 0.90056, 0.90958, 0.91854, & 0.92742, 0.93623, 0.94497, & 0.95364, 0.96223, 0.97074, & 0.97918, 0.98723, 0.99460, & 1.00000 / !<--cjg ! ! Ultra high troposphere resolution data a100/100.00000, 300.00000, 800.00000, & 1762.35235, 3106.43596, 4225.71874, & 4946.40525, 5388.77387, 5708.35540, & 5993.33124, 6277.73673, 6571.49996, & 6877.05339, 7195.14327, 7526.24920, & 7870.82981, 8229.35361, 8602.30193, & 8990.16936, 9393.46399, 9812.70768, & 10248.43625, 10701.19980, 11171.56286, & 11660.10476, 12167.41975, 12694.11735, & 13240.82253, 13808.17600, 14396.83442, & 15007.47066, 15640.77407, 16297.45067, & 16978.22343, 17683.83253, 18415.03554, & 19172.60771, 19957.34218, 20770.05022, & 21559.14829, 22274.03147, 22916.87519, & 23489.70456, 23994.40187, 24432.71365, & 24806.25734, 25116.52754, 25364.90190, & 25552.64670, 25680.92203, 25750.78675, & 25763.20311, 25719.04113, 25619.08274, & 25464.02630, 25254.49482, 24991.06137, & 24674.32737, 24305.11235, 23884.79781, & 23415.77059, 22901.76510, 22347.84738, & 21759.93950, 21144.07284, 20505.73136, & 19849.54271, 19179.31412, 18498.23400, & 17809.06809, 17114.28232, 16416.10343, & 15716.54833, 15017.44246, 14320.43478, & 13627.01116, 12938.50682, 12256.11762, & 11580.91062, 10913.83385, 10255.72526, & 9607.32122, 8969.26427, 8342.11044, & 7726.33606, 7122.34405, 6530.46991, & 5950.98721, 5384.11279, 4830.01153, & 4288.80090, 3760.55514, 3245.30920, & 2743.06250, 2253.78294, 1777.41285, & 1313.88054, 863.12371, 425.13088, & 0.00000, 0.00000 / data b100/0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00000, 0.00000, 0.00000, & 0.00052, 0.00209, 0.00468, & 0.00828, 0.01288, 0.01849, & 0.02508, 0.03266, 0.04121, & 0.05075, 0.06126, 0.07275, & 0.08521, 0.09866, 0.11308, & 0.12850, 0.14490, 0.16230, & 0.18070, 0.20009, 0.22042, & 0.24164, 0.26362, 0.28622, & 0.30926, 0.33258, 0.35605, & 0.37958, 0.40308, 0.42651, & 0.44981, 0.47296, 0.49591, & 0.51862, 0.54109, 0.56327, & 0.58514, 0.60668, 0.62789, & 0.64872, 0.66919, 0.68927, & 0.70895, 0.72822, 0.74709, & 0.76554, 0.78357, 0.80117, & 0.81835, 0.83511, 0.85145, & 0.86736, 0.88286, 0.89794, & 0.91261, 0.92687, 0.94073, & 0.95419, 0.96726, 0.97994, & 0.99223, 1.00000 / data a104/ & 1.8827062944e-01, 7.7977549145e-01, 2.1950593583e+00, & 4.9874566624e+00, 9.8041418997e+00, 1.7019717163e+01, & 2.7216579591e+01, 4.0518628401e+01, 5.6749646818e+01, & 7.5513868331e+01, 9.6315093333e+01, 1.1866706195e+02, & 1.4216835396e+02, 1.6653733709e+02, 1.9161605772e+02, & 2.1735580129e+02, 2.4379516604e+02, 2.7103771847e+02, & 2.9923284173e+02, 3.2856100952e+02, 3.5922338766e+02, & 3.9143507908e+02, 4.2542117983e+02, 4.6141487902e+02, & 4.9965698106e+02, 5.4039638379e+02, 5.8389118154e+02, & 6.3041016829e+02, 6.8023459505e+02, 7.3366009144e+02, & 7.9099869949e+02, 8.5258099392e+02, 9.1875827946e+02, & 9.8990486716e+02, 1.0664204381e+03, 1.1487325074e+03, & 1.2372990044e+03, 1.3326109855e+03, 1.4351954993e+03, & 1.5456186222e+03, 1.6644886848e+03, 1.7924597105e+03, & 1.9302350870e+03, 2.0785714934e+03, 2.2382831070e+03, & 2.4102461133e+03, 2.5954035462e+03, 2.7947704856e+03, & 3.0094396408e+03, 3.2405873512e+03, 3.4894800360e+03, & 3.7574811281e+03, 4.0460585279e+03, 4.3567926151e+03, & 4.6913848588e+03, 5.0516670674e+03, 5.4396113207e+03, & 5.8573406270e+03, 6.3071403487e+03, 6.7914704368e+03, & 7.3129785102e+03, 7.8745138115e+03, 8.4791420557e+03, & 9.1301611750e+03, 9.8311179338e+03, 1.0585825354e+04, & 1.1398380836e+04, 1.2273184781e+04, 1.3214959424e+04, & 1.4228767429e+04, 1.5320029596e+04, 1.6494540743e+04, & 1.7758482452e+04, 1.9118430825e+04, 2.0422798801e+04, & 2.1520147587e+04, 2.2416813461e+04, 2.3118184510e+04, & 2.3628790785e+04, 2.3952411814e+04, 2.4092209011e+04, & 2.4050892106e+04, 2.3830930156e+04, 2.3434818358e+04, & 2.2865410898e+04, 2.2126326004e+04, 2.1222420323e+04, & 2.0160313690e+04, 1.8948920926e+04, 1.7599915822e+04, & 1.6128019809e+04, 1.4550987232e+04, 1.2889169132e+04, & 1.1164595563e+04, 9.4227665517e+03, 7.7259097899e+03, & 6.1538244381e+03, 4.7808126007e+03, 3.5967415552e+03, & 2.5886394104e+03, 1.7415964865e+03, 1.0393721271e+03, & 4.6478852032e+02, 7.0308342481e-13, 0.0000000000e+00 / data b104/ & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 0.0000000000e+00, & 0.0000000000e+00, 0.0000000000e+00, 1.5648447298e-03, & 6.2617046389e-03, 1.4104157933e-02, 2.5118187415e-02, & 3.9340510972e-02, 5.6816335609e-02, 7.7596328431e-02, & 1.0173255472e-01, 1.2927309709e-01, 1.6025505622e-01, & 1.9469566981e-01, 2.3258141217e-01, 2.7385520518e-01, & 3.1840233814e-01, 3.6603639170e-01, 4.1648734767e-01, & 4.6939496013e-01, 5.2431098738e-01, 5.8071350676e-01, & 6.3803478105e-01, 6.9495048840e-01, 7.4963750338e-01, & 7.9975208897e-01, 8.4315257576e-01, 8.8034012292e-01, & 9.1184389721e-01, 9.3821231526e-01, 9.6000677644e-01, & 9.7779792223e-01, 9.9216315122e-01, 1.0000000000e+00 / ! IFS-like L125(top 12 levels removed from IFSL137) data a125/ 64., & 86.895882, 107.415741, 131.425507, 159.279404, 191.338562, & 227.968948, 269.539581, 316.420746, 368.982361, 427.592499, 492.616028, & 564.413452, 643.339905, 729.744141, 823.967834, 926.344910, 1037.201172, & 1156.853638, 1285.610352, 1423.770142, 1571.622925, 1729.448975, 1897.519287, & 2076.095947, 2265.431641, 2465.770508, 2677.348145, 2900.391357, 3135.119385, & 3381.743652, 3640.468262, 3911.490479, 4194.930664, 4490.817383, 4799.149414, & 5119.895020, 5452.990723, 5798.344727, 6156.074219, 6526.946777, 6911.870605, & 7311.869141, 7727.412109, 8159.354004, 8608.525391, 9076.400391, 9562.682617, & 10065.978516, 10584.631836, 11116.662109, 11660.067383, 12211.547852, 12766.873047, & 13324.668945, 13881.331055, 14432.139648, 14975.615234, 15508.256836, 16026.115234, & 16527.322266, 17008.789063, 17467.613281, 17901.621094, 18308.433594, 18685.718750, & 19031.289063, 19343.511719, 19620.042969, 19859.390625, 20059.931641, 20219.664063, & 20337.863281, 20412.308594, 20442.078125, 20425.718750, 20361.816406, 20249.511719, & 20087.085938, 19874.025391, 19608.572266, 19290.226563, 18917.460938, 18489.707031, & 18006.925781, 17471.839844, 16888.687500, 16262.046875, 15596.695313, 14898.453125, & 14173.324219, 13427.769531, 12668.257813, 11901.339844, 11133.304688, 10370.175781, & 9617.515625, 8880.453125, 8163.375000, 7470.343750, 6804.421875, 6168.531250, & 5564.382813, 4993.796875, 4457.375000, 3955.960938, 3489.234375, 3057.265625, & 2659.140625, 2294.242188, 1961.500000, 1659.476563, 1387.546875, 1143.250000, & 926.507813, 734.992188, 568.062500, 424.414063, 302.476563, 202.484375, & 122.101563, 62.781250, 22.835938, 3.757813, 0.000000, 0.000000 / data b125/ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000007, 0.000024, 0.000059, 0.000112, 0.000199, & 0.000340, 0.000562, 0.000890, 0.001353, 0.001992, 0.002857, & 0.003971, 0.005378, 0.007133, 0.009261, 0.011806, 0.014816, & 0.018318, 0.022355, 0.026964, 0.032176, 0.038026, 0.044548, & 0.051773, 0.059728, 0.068448, 0.077958, 0.088286, 0.099462, & 0.111505, 0.124448, 0.138313, 0.153125, 0.168910, 0.185689, & 0.203491, 0.222333, 0.242244, 0.263242, 0.285354, 0.308598, & 0.332939, 0.358254, 0.384363, 0.411125, 0.438391, 0.466003, & 0.493800, 0.521619, 0.549301, 0.576692, 0.603648, 0.630036, & 0.655736, 0.680643, 0.704669, 0.727739, 0.749797, 0.770798, & 0.790717, 0.809536, 0.827256, 0.843881, 0.859432, 0.873929, & 0.887408, 0.899900, 0.911448, 0.922096, 0.931881, 0.940860, & 0.949064, 0.956550, 0.963352, 0.969513, 0.975078, 0.980072, & 0.984542, 0.988500, 0.991984, 0.995003, 0.997630, 1.000000 / select case (km) case (24) ks = 5 do k=1,km+1 ak(k) = a24(k) bk(k) = b24(k) enddo case (26) ks = 7 do k=1,km+1 ak(k) = a26(k) bk(k) = b26(k) enddo case (32) #ifdef OLD_32 ks = 13 ! high-res trop_32 setup #else ks = 7 #endif do k=1,km+1 ak(k) = a32(k) bk(k) = b32(k) enddo case (47) ! ks = 27 ! high-res trop-strat ks = 20 ! Oct 23, 2012 do k=1,km+1 ak(k) = a47(k) bk(k) = b47(k) enddo case (48) ks = 28 do k=1,km+1 ak(k) = a48(k) bk(k) = b48(k) enddo case (52) ks = 35 ! pint = 223 do k=1,km+1 ak(k) = a52(k) bk(k) = b52(k) enddo case (54) ks = 11 ! pint = 109.4 do k=1,km+1 ak(k) = a54(k) bk(k) = b54(k) enddo case (56) ks = 26 do k=1,km+1 ak(k) = a56(k) bk(k) = b56(k) enddo case (60) ks = 37 do k=1,km+1 ak(k) = a60(k) bk(k) = b60(k) enddo case (64) #ifdef GFSL64 ks = 23 #else ks = 46 #endif do k=1,km+1 ak(k) = a64(k) bk(k) = b64(k) enddo !-->cjg case (68) ks = 27 do k=1,km+1 ak(k) = a68(k) bk(k) = b68(k) enddo case (96) ks = 27 do k=1,km+1 ak(k) = a96(k) bk(k) = b96(k) enddo !<--cjg case (100) ks = 38 do k=1,km+1 ak(k) = a100(k) bk(k) = b100(k) enddo case (104) ks = 73 do k=1,km+1 ak(k) = a104(k) bk(k) = b104(k) enddo #ifndef TEST_GWAVES case (10) !-------------------------------------------------- ! Pure sigma-coordinate with uniform spacing in "z" !-------------------------------------------------- ! pt = 2000. ! model top pressure (pascal) ! pt = 100. ! 1 mb press(1) = pt press(km+1) = p0 dlnp = (log(p0) - log(pt)) / real(km) lnpe = log(press(km+1)) do k=km,2,-1 lnpe = lnpe - dlnp press(k) = exp(lnpe) enddo ! Search KS ks = 0 do k=1,km if(press(k) >= pc) then ks = k-1 goto 123 endif enddo 123 continue if(ks /= 0) then do k=1,ks ak(k) = press(k) bk(k) = 0. enddo endif pint = press(ks+1) do k=ks+1,km ak(k) = pint*(press(km)-press(k))/(press(km)-pint) bk(k) = (press(k) - ak(k)) / press(km+1) enddo ak(km+1) = 0. bk(km+1) = 1. ! do k=2,km ! bk(k) = real(k-1) / real(km) ! ak(k) = pt * ( 1. - bk(k) ) ! enddo #endif ! The following 4 selections are better for non-hydrostatic ! Low top: case (31) ptop = 300. pint = 100.E2 call var_dz(km, ak, bk, ptop, ks, pint, 1.035) #ifndef TEST_GWAVES case (41) ptop = 100. pint = 100.E2 call var_hi(km, ak, bk, ptop, ks, pint, 1.035) #endif case (51) ptop = 100. pint = 100.E2 call var_hi(km, ak, bk, ptop, ks, pint, 1.035) ! Mid-top: case (55) ptop = 10. pint = 100.E2 ! call var_dz(km, ak, bk, ptop, ks, pint, 1.035) call var_hi(km, ak, bk, ptop, ks, pint, 1.035) #ifdef USE_GFSL63 ! GFS L64 equivalent setting case (63) ks = 23 ptop = a63(1) pint = a63(ks+1) do k=1,km+1 ak(k) = a63(k) bk(k) = b63(k) enddo #else case (63) ptop = 1. ! high top pint = 100.E2 call var_hi(km, ak, bk, ptop, ks, pint, 1.035) #endif ! NGGPS_GFS case (91) pint = 100.E2 ptop = 40. call var_gfs(km, ak, bk, ptop, ks, pint, 1.029) ! call var_gfs(km, ak, bk, ptop, ks, pint, 1.03) case (95) ! Mid-top settings: pint = 100.E2 ptop = 20. call var_gfs(km, ak, bk, ptop, ks, pint, 1.028) case (127) ptop = 1. pint = 75.E2 call var_gfs(km, ak, bk, ptop, ks, pint, 1.028) ! IFS-like L125 case (125) ks = 33 ptop = a125(1) pint = a125(ks+1) do k=1,km+1 ak(k) = a125(k) bk(k) = b125(k) enddo case default #ifdef TEST_GWAVES !-------------------------------------------------- ! Pure sigma-coordinate with uniform spacing in "z" !-------------------------------------------------- call gw_1d(km, 1000.E2, ak, bk, ptop, 10.E3, pt1) ks = 0 pint = ak(1) #else !---------------------------------------------------------------- ! Sigma-coordinate with uniform spacing in sigma and ptop = 1 mb !---------------------------------------------------------------- pt = 100. ! One pressure layer ks = 1 ! pint = pt + 0.5*1.E5/real(km) ! SJL: 20120327 pint = pt + 1.E5/real(km) ak(1) = pt bk(1) = 0. ak(2) = pint bk(2) = 0. do k=3,km+1 bk(k) = real(k-2) / real(km-1) ak(k) = pint - bk(k)*pint enddo ak(km+1) = 0. bk(km+1) = 1. #endif end select ptop = ak(1) pint = ak(ks+1) end subroutine set_eta #endif !>@brief The subroutine 'set_external_eta' sets 'ptop' (model top) and !! 'ks' (first level of pure pressure coordinates given the coefficients !! 'ak' and 'bk' subroutine set_external_eta(ak, bk, ptop, ks) implicit none real, intent(in) :: ak(:) real, intent(in) :: bk(:) real, intent(out) :: ptop !< model top (Pa) integer, intent(out) :: ks !< number of pure p layers !--- local variables integer :: k real :: eps = 1.d-7 ptop = ak(1) ks = 1 do k = 1, size(bk(:)) if (bk(k).lt.eps) ks = k enddo !--- change ks to layers from levels ks = ks - 1 if (is_master()) write(6,*) ' ptop & ks ', ptop, ks end subroutine set_external_eta subroutine var_les(km, ak, bk, ptop, ks, pint, s_rate) implicit none integer, intent(in):: km real, intent(in):: ptop real, intent(in):: s_rate !< between [1. 1.1] real, intent(out):: ak(km+1), bk(km+1) real, intent(inout):: pint integer, intent(out):: ks ! Local real, parameter:: p00 = 1.E5 real, dimension(km+1):: ze, pe1, peln, eta real, dimension(km):: dz, s_fac, dlnp, pm, dp, dk real ztop, t0, dz0, sum1, tmp1 real ep, es, alpha, beta, gama real, parameter:: akap = 2./7. !---- Tunable parameters: real:: k_inc = 10 !< number of layers from bottom up to near const dz region real:: s0 = 0.8 !< lowest layer stretch factor !----------------------- real:: s_inc integer k pe1(1) = ptop peln(1) = log(pe1(1)) pe1(km+1) = p00 peln(km+1) = log(pe1(km+1)) t0 = 273. ztop = rdgas/grav*t0*(peln(km+1) - peln(1)) s_inc = (1.-s0) / real(k_inc) s_fac(km) = s0 do k=km-1, km-k_inc, -1 s_fac(k) = s_fac(k+1) + s_inc enddo s_fac(km-k_inc-1) = 0.5*(s_fac(km-k_inc) + s_rate) do k=km-k_inc-2, 5, -1 s_fac(k) = s_rate * s_fac(k+1) enddo s_fac(4) = 0.5*(1.1+s_rate)*s_fac(5) s_fac(3) = 1.1 *s_fac(4) s_fac(2) = 1.1 *s_fac(3) s_fac(1) = 1.1 *s_fac(2) sum1 = 0. do k=1,km sum1 = sum1 + s_fac(k) enddo dz0 = ztop / sum1 do k=1,km dz(k) = s_fac(k) * dz0 enddo ze(km+1) = 0. do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo ! Re-scale dz with the stretched ztop do k=1,km dz(k) = dz(k) * (ztop/ze(1)) enddo do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo ! ze(1) = ztop if ( is_master() ) then write(*,*) 'var_les: computed model top (m)=', ztop, ' bottom/top dz=', dz(km), dz(1) ! do k=1,km ! write(*,*) k, s_fac(k) ! enddo endif call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2) ! Given z --> p do k=1,km dz(k) = ze(k) - ze(k+1) dlnp(k) = grav*dz(k) / (rdgas*t0) !write(*,*) k, dz(k) enddo do k=2,km peln(k) = peln(k-1) + dlnp(k-1) pe1(k) = exp(peln(k)) enddo ! Pe(k) = ak(k) + bk(k) * PS ! Locate pint and KS ks = 0 do k=2,km if ( pint < pe1(k)) then ks = k-1 exit endif enddo if ( is_master() ) then write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1) endif pint = pe1(ks+1) do k=1,km+1 eta(k) = pe1(k) / pe1(km+1) enddo ep = eta(ks+1) es = eta(km) ! es = 1. alpha = (ep**2-2.*ep*es) / (es-ep)**2 beta = 2.*ep*es**2 / (es-ep)**2 gama = -(ep*es)**2 / (es-ep)**2 ! Pure pressure: do k=1,ks+1 ak(k) = eta(k)*1.e5 bk(k) = 0. enddo do k=ks+2, km ak(k) = alpha*eta(k) + beta + gama/eta(k) ak(k) = ak(k)*1.e5 enddo ak(km+1) = 0. do k=ks+2, km bk(k) = (pe1(k) - ak(k))/pe1(km+1) enddo bk(km+1) = 1. if ( is_master() ) then ! write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100. ! do k=1,km ! pm(k) = 0.5*(pe1(k)+pe1(k+1))/100. ! write(*,*) k, pm(k), dz(k) ! enddo tmp1 = ak(ks+1) do k=ks+1,km tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) ) enddo write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100. write(*,800) (pm(k), k=km,1,-1) endif do k=1,km dp(k) = (pe1(k+1) - pe1(k))/100. dk(k) = pe1(k+1)**akap - pe1(k)**akap enddo 800 format(1x,5(1x,f9.4)) end subroutine var_les subroutine var_gfs(km, ak, bk, ptop, ks, pint, s_rate) integer, intent(in):: km real, intent(in):: ptop real, intent(in):: s_rate !< between [1. 1.1] real, intent(out):: ak(km+1), bk(km+1) real, intent(inout):: pint integer, intent(out):: ks ! Local real, parameter:: p00 = 1.E5 real, dimension(km+1):: ze, pe1, peln, eta real, dimension(km):: dz, s_fac, dlnp real ztop, t0, dz0, sum1, tmp1 real ep, es, alpha, beta, gama !---- Tunable parameters: integer:: k_inc = 25 !< number of layers from bottom up to near const dz region real:: s0 = 0.13 !< lowest layer stretch factor !----------------------- real:: s_inc integer k pe1(1) = ptop peln(1) = log(pe1(1)) pe1(km+1) = p00 peln(km+1) = log(pe1(km+1)) t0 = 270. ztop = rdgas/grav*t0*(peln(km+1) - peln(1)) s_inc = (1.-s0) / real(k_inc) s_fac(km) = s0 do k=km-1, km-k_inc, -1 s_fac(k) = s_fac(k+1) + s_inc enddo do k=km-k_inc-1, 9, -1 s_fac(k) = s_rate * s_fac(k+1) enddo s_fac(8) = 0.5*(1.1+s_rate)*s_fac(9) s_fac(7) = 1.10*s_fac(8) s_fac(6) = 1.15*s_fac(7) s_fac(5) = 1.20*s_fac(6) s_fac(4) = 1.26*s_fac(5) s_fac(3) = 1.33*s_fac(4) s_fac(2) = 1.41*s_fac(3) s_fac(1) = 1.60*s_fac(2) sum1 = 0. do k=1,km sum1 = sum1 + s_fac(k) enddo dz0 = ztop / sum1 do k=1,km dz(k) = s_fac(k) * dz0 enddo ze(km+1) = 0. do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo ! Re-scale dz with the stretched ztop do k=1,km dz(k) = dz(k) * (ztop/ze(1)) enddo do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo ! ze(1) = ztop if ( is_master() ) then write(*,*) 'var_gfs: computed model top (m)=', ztop*0.001, ' bottom/top dz=', dz(km), dz(1) ! do k=1,km ! write(*,*) k, s_fac(k) ! enddo endif ! call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 3) ! Given z --> p do k=1,km dz(k) = ze(k) - ze(k+1) dlnp(k) = grav*dz(k) / (rdgas*t0) enddo do k=2,km peln(k) = peln(k-1) + dlnp(k-1) pe1(k) = exp(peln(k)) enddo ! Pe(k) = ak(k) + bk(k) * PS ! Locate pint and KS ks = 0 do k=2,km if ( pint < pe1(k)) then ks = k-1 exit endif enddo if ( is_master() ) then write(*,*) 'For (input) PINT=', 0.01*pint, ' KS=', ks, 'pint(computed)=', 0.01*pe1(ks+1) write(*,*) 'ptop =', ptop endif pint = pe1(ks+1) #ifdef NO_UKMO_HB do k=1,ks+1 ak(k) = pe1(k) bk(k) = 0. enddo do k=ks+2,km+1 bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint) ! bk == sigma ak(k) = pe1(k) - bk(k) * pe1(km+1) enddo bk(km+1) = 1. ak(km+1) = 0. #else ! Problematic for non-hydrostatic do k=1,km+1 eta(k) = pe1(k) / pe1(km+1) enddo ep = eta(ks+1) es = eta(km) ! es = 1. alpha = (ep**2-2.*ep*es) / (es-ep)**2 beta = 2.*ep*es**2 / (es-ep)**2 gama = -(ep*es)**2 / (es-ep)**2 ! Pure pressure: do k=1,ks+1 ak(k) = eta(k)*1.e5 bk(k) = 0. enddo do k=ks+2, km ak(k) = alpha*eta(k) + beta + gama/eta(k) ak(k) = ak(k)*1.e5 enddo ak(km+1) = 0. do k=ks+2, km bk(k) = (pe1(k) - ak(k))/pe1(km+1) enddo bk(km+1) = 1. #endif if ( is_master() ) then write(*,*) 'KS=', ks, 'PINT (mb)=', pint/100. do k=1,km write(*,*) k, 0.5*(pe1(k)+pe1(k+1))/100., dz(k) enddo tmp1 = ak(ks+1) do k=ks+1,km tmp1 = max(tmp1, (ak(k)-ak(k+1))/max(1.E-5, (bk(k+1)-bk(k))) ) enddo write(*,*) 'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100. endif end subroutine var_gfs subroutine var_hi(km, ak, bk, ptop, ks, pint, s_rate) integer, intent(in):: km real, intent(in):: ptop real, intent(in):: s_rate !< between [1. 1.1] real, intent(out):: ak(km+1), bk(km+1) real, intent(inout):: pint integer, intent(out):: ks ! Local real, parameter:: p00 = 1.E5 real, dimension(km+1):: ze, pe1, peln, eta real, dimension(km):: dz, s_fac, dlnp real ztop, t0, dz0, sum1, tmp1 real ep, es, alpha, beta, gama !---- Tunable parameters: integer:: k_inc = 15 ! 1.E-8 ) then pf(1) = (ph(2) - ph(1)) / log(ph(2)/ph(1)) else pf(1) = (ph(2) - ph(1)) * kappa/(kappa+1.) endif do k=2,npz pf(k) = (ph(k+1) - ph(k)) / log(ph(k+1)/ph(k)) enddo end subroutine get_eta_level subroutine compute_dz(km, ztop, dz) integer, intent(in):: km real, intent(in):: ztop ! try 50.E3 real, intent(out):: dz(km) !------------------------------ real ze(km+1), dzt(km) integer k ! ztop = 30.E3 dz(1) = ztop / real(km) dz(km) = 0.5*dz(1) do k=2,km-1 dz(k) = dz(1) enddo ! Top: dz(1) = 2.*dz(2) ze(km+1) = 0. do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo if ( is_master() ) then write(*,*) 'Hybrid_z: dz, zm' do k=1,km dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000. write(*,*) k, dz(k), dzt(k) enddo endif end subroutine compute_dz subroutine compute_dz_var(km, ztop, dz) integer, intent(in):: km real, intent(in):: ztop ! try 50.E3 real, intent(out):: dz(km) !------------------------------ real, parameter:: s_rate = 1.0 real ze(km+1) real s_fac(km) real sum1, dz0 integer k s_fac(km ) = 0.125 s_fac(km-1) = 0.20 s_fac(km-2) = 0.30 s_fac(km-3) = 0.40 s_fac(km-4) = 0.50 s_fac(km-5) = 0.60 s_fac(km-6) = 0.70 s_fac(km-7) = 0.80 s_fac(km-8) = 0.90 s_fac(km-9) = 1. do k=km-10, 9, -1 s_fac(k) = s_rate * s_fac(k+1) enddo s_fac(8) = 1.05*s_fac(9) s_fac(7) = 1.1 *s_fac(8) s_fac(6) = 1.15*s_fac(7) s_fac(5) = 1.2 *s_fac(6) s_fac(4) = 1.3 *s_fac(5) s_fac(3) = 1.4 *s_fac(4) s_fac(2) = 1.5 *s_fac(3) s_fac(1) = 1.6 *s_fac(2) sum1 = 0. do k=1,km sum1 = sum1 + s_fac(k) enddo dz0 = ztop / sum1 do k=1,km dz(k) = s_fac(k) * dz0 enddo ze(1) = ztop ze(km+1) = 0. do k=km,2,-1 ze(k) = ze(k+1) + dz(k) enddo ! Re-scale dz with the stretched ztop do k=1,km dz(k) = dz(k) * (ztop/ze(1)) enddo do k=km,2,-1 ze(k) = ze(k+1) + dz(k) enddo call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 2) do k=1,km dz(k) = ze(k) - ze(k+1) enddo end subroutine compute_dz_var subroutine compute_dz_L32(km, ztop, dz) integer, intent(in):: km real, intent(out):: dz(km) real, intent(out):: ztop ! try 50.E3 !------------------------------ real dzt(km) real ze(km+1) real dz0, dz1, dz2 real z0, z1, z2 integer k, k0, k1, k2, n !------------------- k2 = 8 z2 = 30.E3 !------------------- k1 = 21 z1 = 10.0E3 !------------------- k0 = 2 z0 = 0. dz0 = 75. ! meters !------------------- ! Treat the surface layer as a special layer ze(1) = z0 dz(1) = dz0 ze(2) = dz(1) dz0 = 1.5*dz0 dz(2) = dz0 ze(3) = ze(2) + dz(2) dz1 = 2.*(z1-ze(3) - k1*dz0) / (k1*(k1-1)) do k=k0+1,k0+k1 dz(k) = dz0 + (k-k0)*dz1 ze(k+1) = ze(k) + dz(k) enddo dz0 = dz(k1+k0) dz2 = 2.*(z2-ze(k0+k1+1)-k2*dz0) / (k2*(k2-1)) do k=k0+k1+1,k0+k1+k2 dz(k) = dz0 + (k-k0-k1)*dz2 ze(k+1) = ze(k) + dz(k) enddo dz(km) = 2.*dz(km-1) ztop = ze(km) + dz(km) ze(km+1) = ze(km) + dz(km) call zflip (dz, 1, km) ze(km+1) = 0. do k=km,1,-1 ze(k) = ze(k+1) + dz(k) enddo ! if ( is_master() ) then ! write(*,*) 'Hybrid_z: dz, zm' ! do k=1,km ! dzt(k) = 0.5*(ze(k)+ze(k+1)) / 1000. ! write(*,*) k, dz(k), dzt(k) ! enddo ! endif end subroutine compute_dz_L32 subroutine compute_dz_L101(km, ztop, dz) integer, intent(in):: km ! km==101 real, intent(out):: dz(km) real, intent(out):: ztop ! try 30.E3 !------------------------------ real ze(km+1) real dz0, dz1 real:: stretch_f = 1.16 integer k, k0, k1 k1 = 2 k0 = 25 dz0 = 40. ! meters ze(km+1) = 0. do k=km, k0, -1 dz(k) = dz0 ze(k) = ze(k+1) + dz(k) enddo do k=k0+1, k1, -1 dz(k) = stretch_f * dz(k+1) ze(k) = ze(k+1) + dz(k) enddo dz(1) = 4.0*dz(2) ze(1) = ze(2) + dz(1) ztop = ze(1) if ( is_master() ) then write(*,*) 'Hybrid_z: dz, ze' do k=1,km write(*,*) k, 0.001*dz(k), 0.001*ze(k) enddo ! ztop (km) = 20.2859154 write(*,*) 'ztop (km) =', ztop * 0.001 endif end subroutine compute_dz_L101 subroutine set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3) integer, intent(in):: is, ie, js, je, ng, km real, intent(in):: rgrav, ztop real, intent(in):: dz(km) !< Reference vertical resolution for zs=0 real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng) real, intent(inout):: ze(is:ie,js:je,km+1) real, optional, intent(out):: dz3(is-ng:ie+ng,js-ng:je+ng,km) ! local logical:: filter_xy = .false. real, allocatable:: delz(:,:,:) integer ntimes real zint real:: z1(is:ie,js:je) real:: z(km+1) real sig, z_rat integer ks(is:ie,js:je) integer i, j, k, ks_min, kint z(km+1) = 0. do k=km,1,-1 z(k) = z(k+1) + dz(k) enddo do j=js,je do i=is,ie ze(i,j, 1) = ztop ze(i,j,km+1) = hs(i,j) * rgrav enddo enddo do k=2,km do j=js,je do i=is,ie ze(i,j,k) = z(k) enddo enddo enddo ! Set interface: #ifndef USE_VAR_ZINT zint = 12.0E3 ntimes = 2 kint = 2 do k=2,km if ( z(k)<=zint ) then kint = k exit endif enddo if ( is_master() ) write(*,*) 'Z_coord interface set at k=',kint, ' ZE=', z(kint) do j=js,je do i=is,ie z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1)) do k=km,kint+1,-1 ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat enddo !-------------------------------------- ! Apply vertical smoother locally to dz !-------------------------------------- call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes) enddo enddo #else ! ZINT is a function of local terrain ntimes = 4 do j=js,je do i=is,ie z1(i,j) = dim(ze(i,j,km+1), 2500.) + 7500. enddo enddo ks_min = km do j=js,je do i=is,ie do k=km,2,-1 if ( z(k)>=z1(i,j) ) then ks(i,j) = k ks_min = min(ks_min, k) go to 555 endif enddo 555 continue enddo enddo do j=js,je do i=is,ie kint = ks(i,j) + 1 z_rat = (ze(i,j,kint)-ze(i,j,km+1)) / (z(kint)-z(km+1)) do k=km,kint+1,-1 ze(i,j,k) = ze(i,j,k+1) + dz(k)*z_rat enddo !-------------------------------------- ! Apply vertical smoother locally to dz !-------------------------------------- call sm1_edge(is, ie, js, je, km, i, j, ze, ntimes) enddo enddo #endif #ifdef DEV_ETA if ( filter_xy ) then allocate (delz(isd:ied, jsd:jed, km) ) ntimes = 2 do k=1,km do j=js,je do i=is,ie delz(i,j,k) = ze(i,j,k+1) - ze(i,j,k) enddo enddo enddo call del2_cubed(delz, 0.2*da_min, npx, npy, km, ntimes) do k=km,2,-1 do j=js,je do i=is,ie ze(i,j,k) = ze(i,j,k+1) - delz(i,j,k) enddo enddo enddo deallocate ( delz ) endif #endif if ( present(dz3) ) then do k=1,km do j=js,je do i=is,ie dz3(i,j,k) = ze(i,j,k+1) - ze(i,j,k) enddo enddo enddo endif end subroutine set_hybrid_z subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes) integer, intent(in):: is, ie, js, je, km integer, intent(in):: ntimes, i, j real, intent(inout):: ze(is:ie,js:je,km+1) ! local: real, parameter:: df = 0.25 real dz(km) real flux(km+1) integer k, n, k1, k2 k2 = km-1 do k=1,km dz(k) = ze(i,j,k+1) - ze(i,j,k) enddo do n=1,ntimes k1 = 2 + (ntimes-n) flux(k1 ) = 0. flux(k2+1) = 0. do k=k1+1,k2 flux(k) = df*(dz(k) - dz(k-1)) enddo do k=k1,k2 dz(k) = dz(k) - flux(k) + flux(k+1) enddo enddo do k=km,1,-1 ze(i,j,k) = ze(i,j,k+1) - dz(k) enddo end subroutine sm1_edge subroutine gw_1d(km, p0, ak, bk, ptop, ztop, pt1) integer, intent(in):: km real, intent(in):: p0, ztop real, intent(inout):: ptop real, intent(inout):: ak(km+1), bk(km+1) real, intent(out):: pt1(km) ! Local logical:: isothermal real, dimension(km+1):: ze, pe1, pk1 real, dimension(km):: dz1 real t0, n2, s0 integer k ! Set up vertical coordinare with constant del-z spacing: isothermal = .false. t0 = 300. if ( isothermal ) then n2 = grav**2/(cp_air*t0) else n2 = 0.0001 endif s0 = grav*grav / (cp_air*n2) ze(km+1) = 0. do k=km,1,-1 dz1(k) = ztop / real(km) ze(k) = ze(k+1) + dz1(k) enddo ! Given z --> p do k=1,km+1 pe1(k) = p0*( (1.-s0/t0) + s0/t0*exp(-n2*ze(k)/grav) )**(1./kappa) enddo ptop = pe1(1) ! if ( is_master() ) write(*,*) 'GW_1D: computed model top (pa)=', ptop ! Set up "sigma" coordinate ak(1) = pe1(1) bk(1) = 0. do k=2,km bk(k) = (pe1(k) - pe1(1)) / (pe1(km+1)-pe1(1)) ! bk == sigma ak(k) = pe1(1)*(1.-bk(k)) enddo ak(km+1) = 0. bk(km+1) = 1. do k=1,km+1 pk1(k) = pe1(k) ** kappa enddo ! Compute volume mean potential temperature with hydrostatic eqn: do k=1,km pt1(k) = grav*dz1(k) / ( cp_air*(pk1(k+1)-pk1(k)) ) enddo end subroutine gw_1d subroutine zflip(q, im, km) integer, intent(in):: im, km real, intent(inout):: q(im,km) !--- integer i, k real qtmp do i = 1, im do k = 1, (km+1)/2 qtmp = q(i,k) q(i,k) = q(i,km+1-k) q(i,km+1-k) = qtmp end do end do end subroutine zflip end module fv_eta_mod