CCPP SciDoc v7.0.x  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
ecmwf_ngw.F90
1
3
4module ecmwf_ngw
5
6
7contains
8!------------------------------------------------------------------------------------------
9! April 2025 Adding ECMWF Non-stationary gravity wave scheme option by Bo Yang
10!------------------------------------------------------------------------------------------
11! different orientation at vertical
12! 1 is the highest level for ECMWF, 1 is the lowest level for GFS
13! Vertical levels are reversed after entering the subroutine ecmwf_ngw_emc and reversed back before exiting
14! This non-stationary GWD module was obtained by Fanglin Yang from ECMWF, with permission for operational use at NCEP. We would
15! like to thank Andy Brown, Michael Sleigh, Peter Bechtold, and Nils Wedi at ECMWF for their support in porting this code to the
16! UFS.
17!! Original Fortran Code by J. SCINOCCIA
18! Rewritten in IFS format by A. ORR E.C.M.W.F. August 2008
19! PURPOSE
20! -------
21
22! THIS ROUTINE COMPUTES NON-OROGRAPHIC GRAVITY WAVE DRAG
23! AFTER SCINOCCA (2003) AND Mc LANDRESS AND SCINOCCIA (JAS 2005)
24! HYDROSTATIC NON-ROTATIONAL SIMPLIFIED VERSION OF THE
25! WARNER AND MCINTYRE (1996) NON-OROGRAPHIC GRAVITY WAVE PARAMETERIZATION
26! CONSTANTS HAVE BEEN OPTIMIZED FOLLOWING M. ERN ET AL. (ATMOS. CHEM. PHYS. 2006)
27
28! REFERENCE: Orr, A., P. Bechtold, J. Scinoccia, M. Ern, M. Janiskova, 2010:
29! Improved middle atmosphere climate and analysis in the ECMWF forecasting system
30! through a non-orographic gravity wave parametrization. J. Climate., 23, 5905-5926.
31
32! LAUNCH SPECTRUM - GENERALIZED DESAUBIES
33! INCLUDES A CRITICAL-LEVEL CORRECTION THAT PREVENTS THE
34! MOMEMTUM DEPOSITION IN EACH AZIMUTH FROM DRIVING THE FLOW TO SPEEDS FASTER
35! THAN THE PHASE SPEED OF THE WAVES, I.E. WHEN WAVES BREAK THEY DRAG THE MEAN
36! FLOW TOWARDS THEIR PHASE SPEED - NOT PAST IT.
37!!! different orientation for vertical
38!!! 1 is the highest level for ECMWF, 1 is the lowest level for GFS
39!---------------------------------------------------
40! subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, &
41! tau_ngw, tm , um, vm, qm, prsl, prsi, zmet, zmeti, prslk, &
42! xlatd, sinlat, coslat, &
43! pdudt, pdvdt, pdtdt, dked, zngw)
44
45! (C) Copyright 1989- ECMWF.
46!
47! This software is licensed under the terms of the Apache Licence Version 2.0
48! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
49!
50! In applying this licence, ECMWF does not waive the privileges and immunities
51! granted to it by virtue of its status as an intergovernmental organisation
52! nor does it submit to any jurisdiction.
53
54 subroutine ecmwf_ngw_emc(mpi_id, master, KLON, KLEV, kdt, PTSTEP, DX, &
55 tau_ngw, PTM11, PUM11, PVM11, qm1, PAPM11, PAPHM11, PGEO11, zmeti1, prslk1, &
56 xlatd, sinlat, coslat, &
57 PTENU, PTENV, pdtdt, dked, zngw)
58
59
60!
61 use machine, only : kind_phys
62
63
64 use cires_ugwpv1_module,only : psrc => knob_ugwp_palaunch
65
66 use cires_ugwpv1_module,only : maxdudt, maxdtdt, max_eps, dked_min, dked_max
67!! maxdudt=250.e-5; maxdtdt=15.e-2; dked_min=0.01; dked_max=250.0
68!! max_eps=max_kdis*4.e-4=450*4.e-4
69
70 use ugwp_common , only : rgrav, grav, cpd, rd, rv, rcpdl, grav2cpd, &
71 omega2, rcpd, rcpd2, pi, pi2, fv, &
72 rad_to_deg, deg_to_rad, &
73 rdi, gor, grcp, gocp, &
74 bnv2min, bnv2max, dw2min, velmin, gr2, &
75 hpscale, rhp, rh4, grav2, rgrav2, mkzmin, mkz2min
76
77!!! grav=con_g; rgrav=1/grav; cpd=con_cp; rd=con_rd; rv=con_rv; rcpdl=cpd*rgrav=cpd/g
78!!!grav2cpd= grav*grcp=grav*grav*rcpd= g**2/cpd omega=con_omega; omega2=2*omega1
79!! rcpd2=0.5*rcpd=1/(2*cpd)
80!! pi=con_pi; pi2=2*pi
81!! fv=con_fvirt; rad_to_deg=180.0/pi; deg_to_rad=pi2/180.0
82!! rdi= 1.0/rd; gor=grav/rd; grcp= grav*rcpd=g/Cpd; gocp=grcp=grav*rcpd= g/cpd
83!! bnv2min=(pi2/1800.)*(pi2/1800.); bnv2max=(pi2/30.)*(pi2/30.)
84!! dw2min=1.0, velmin=sqrt(dw2min)
85!! gr2=grav*gor=grav*grav/rd=g**2/rd
86!! hpscale=7000; rhp = 1./hpscale = 1/7000; rh4=rhp2*rhp2=(0.5*rhp)**2=1/4*1/7000=1/28000
87!! rgrav2=rgrav*rgrav=1/(g**2); grav2=grav+grav=2*g; mkzmin=pi2/80.0e-3=2*pi/80.0e3
88!! mkz2min= mkzmin*mkzmin=(2*pi/80.0E3)**2
89
90
91 use ugwp_wmsdis_init, only : v_kxw, rv_kxw, v_kxw2, tamp_mpa, tau_min, ucrit, &
92 gw_eff, &
93! zms, &
94 zci4, zci3, zci2, &
95 rimin, sc2, sc2u, ric
96! v_kxw=kxw=pi2/lhmet=pi2/200e3; rv_kxw=200e3/pi2
97! v_kxw2=v_kxw*v_kxw=(pi2/200e3)**2; tamp_mpa=knob_ugwp_tauamp amplitude for GEOS-5/MERRA-2
98! tau_min=min of GW MF 0.25mPa
99! ucrit=cdmin=2e-2/mkzmax=2e-2/((2pie)/500) = 10/(2pie)
100! gw_eff=effac=1.0
101! zci4=(zms*zci(inc))**4; zci2=(zms*zci(inc))**2; zci3(inc)=(zms*zci(inc))**3
102! rimin=-10.0; sc2=lturb*lturb=30m*30m; sc2u=ulturb*ulturb=150*150
103
104
105 implicit none
106
107!in
108!work
109
110
111! integer, intent(in) :: KLAUNCH ! index for launch level
112 integer, intent(in) :: KLEV ! vertical level
113 integer, intent(in) :: KLON ! horiz tiles
114 integer, intent(in) :: mpi_id, master, kdt
115
116 real(kind=kind_phys) ,intent(in) :: ptstep ! model time step
117 real(kind=kind_phys) ,intent(in) :: dx(klon) ! model grid size
118
119 real(kind=kind_phys) ,intent(in) :: tau_ngw(klon)
120
121 real(kind=kind_phys) ,intent(in) :: pvm11(klon,klev) ! meridional wind
122 real(kind=kind_phys) ,intent(in) :: pum11(klon,klev) ! zonal wind
123 real(kind=kind_phys) ,intent(in) :: qm1(klon,klev) ! spec. humidity
124 real(kind=kind_phys) ,intent(in) :: ptm11(klon,klev) ! kinetic temperature
125 real(kind=kind_phys) ,intent(in) :: papm11(klon,klev) ! mid-layer pressure
126 real(kind=kind_phys) ,intent(in) :: paphm11(klon,klev+1) ! interface pressure
127 real(kind=kind_phys) ,intent(in) :: pgeo11(klon,klev) ! full model level geopotential in meters
128
129 real(kind=kind_phys) :: pvm1(klon,klev) ! meridional wind
130 real(kind=kind_phys) :: pum1(klon,klev) ! zonal wind
131 real(kind=kind_phys) :: qm(klon,klev) ! spec. humidity
132 real(kind=kind_phys) :: ptm1(klon,klev) ! kinetic temperature
133 real(kind=kind_phys) :: papm1(klon,klev) ! mid-layer pressure
134 real(kind=kind_phys) :: paphm1(klon,klev+1) ! interface pressure
135 real(kind=kind_phys) :: pgeo1(klon,klev) ! full model level geopotential in meters
136
137! real(kind=kind_phys) :: PGAW(KLON) !normalised gaussian quadrature weight/nb of longitude pts
138 ! local sub-area == 4*RPI*RA**2 * PGAW
139! real(kind=kind_phys) ,intent(in) :: PPRECIP(KLON) ! total surface precipitation
140
141
142
143
144 real(kind=kind_phys) ,intent(in) :: prslk1(klon,klev) ! mid-layer exner function
145 real(kind=kind_phys) :: prslk(klon,klev) ! mid-layer exner function
146! real(kind=kind_phys) ,intent(in) :: zmet(KLON,KLEV) ! meters phil =philg/grav ! use PGEO1 instead
147
148 real(kind=kind_phys) ,intent(in) :: zmeti1(klon,klev+1) ! interface geopi/meters
149 real(kind=kind_phys) :: zmeti(klon,klev+1) ! interface geopi/meters
150 real(kind=kind_phys) ,intent(in) :: xlatd(klon) ! xlat_d in degrees
151 real(kind=kind_phys) ,intent(in) :: sinlat(klon)
152 real(kind=kind_phys) ,intent(in) :: coslat(klon)
153!
154! out-gw effects
155!
156 real(kind=kind_phys) ,intent(out) :: ptenu(klon,klev) ! zonal momentum tendency
157 real(kind=kind_phys) ,intent(out) :: ptenv(klon,klev) ! meridional momentum tendency
158 real(kind=kind_phys) ,intent(out) :: pdtdt(klon,klev) ! gw-heating (u*ax+v*ay)/cp and cooling
159
160
161 real(kind=kind_phys) :: pfluxu(klon,klev+1) ! zonal momentum tendency
162 real(kind=kind_phys) :: pfluxv(klon,klev+1) ! meridional momentum tendency
163
164 real(kind=kind_phys) ,intent(out) :: dked(klon,klev) ! gw-eddy diffusion
165
166 real(kind=kind_phys) ,intent(out) :: zngw(klon) ! launch height
167!
168!work
169 INTEGER, PARAMETER :: IAZIDIM=4 !number of azimuths
170 INTEGER, PARAMETER :: INCDIM=20 !number of discretized c spectral elements in launch spectrum
171
172 REAL(kind=kind_phys), PARAMETER :: ra=6370.e3 !half-model level zonal velocity
173 REAL(kind=kind_phys), PARAMETER :: gptwo=2.0 ! 2p in equation
174
175 REAL(kind=kind_phys) :: zuhm1(klon,klev) !half-model level zonal velocity
176 REAL(kind=kind_phys) :: zvhm1(klon,klev) !half-model level meridional velocity
177
178 REAL(kind=kind_phys) :: zbvfhm1(klon,klev) !half-model level Brunt-Vaisalla frequency
179 REAL(kind=kind_phys) :: zrhohm1(klon,klev) !half-model level density
180 REAL(kind=kind_phys) :: zx(incdim) !coordinate transformation
181 REAL(kind=kind_phys) :: zci(incdim) !phase speed element
182 REAL(kind=kind_phys) :: zdci(incdim)
183 REAL(kind=kind_phys) :: zui(klon,klev,iazidim) !intrinsic velocity
184 REAL(kind=kind_phys) :: zul(klon,iazidim) !velocity in azimuthal direction at launch level
185 REAL(kind=kind_phys) :: zbvfl(klon) !buoyancy at launch level
186 REAL(kind=kind_phys) :: zcosang(iazidim) !cos of azimuth angle
187 REAL(kind=kind_phys) :: zsinang(iazidim) !sin of azimuth angle
188 REAL(kind=kind_phys) :: zfct(klon,klev)
189 REAL(kind=kind_phys) :: zfnorm(klon) !normalisation factor (A)
190 REAL(kind=kind_phys) :: zci_min(klon,iazidim)
191 REAL(kind=kind_phys) :: zthm1(klon,klev) !temperature on half-model levels
192 REAL(kind=kind_phys) :: zflux(klon,incdim,iazidim) !momentum flux at each vertical level and azimuth
193 REAL(kind=kind_phys) :: zpu(klon,klev,iazidim) !momentum flux
194 REAL(kind=kind_phys) :: zdfl(klon,klev,iazidim)
195 REAL(kind=kind_phys) :: zact(klon,incdim,iazidim) !if =1 then critical level encountered
196 REAL(kind=kind_phys) :: zacc(klon,incdim,iazidim)
197 REAL(kind=kind_phys) :: zcrt(klon,klev,iazidim)
198
199 INTEGER :: ILAUNCH !model level from which GW spectrum is launched
200 INTEGER :: INC, JK, JL, IAZI
201
202
203 REAL(KIND=kind_phys) :: zradtodeg, zgelatdeg
204 REAL(KIND=kind_phys) :: zcimin, zcimax
205 REAL(KIND=kind_phys) :: zgam, zpexp, zxmax, zxmin, zxran
206 REAL(KIND=kind_phys) :: zdx
207
208 REAL(KIND=kind_phys) :: zx1, zx2, zdxa, zdxb, zdxs
209 REAL(KIND=kind_phys) :: zang, zaz_fct, znorm, zang1, ztx
210 REAL(KIND=kind_phys) :: zu, zcin, zcpeak
211 REAL(KIND=kind_phys) :: zcin4, zbvfl4, zcin2, zbvfl2, zcin3, zbvfl3, zcinc
212 REAL(KIND=kind_phys) :: zatmp, zfluxs, zdep, zfluxsq, zulm, zdft, ze1, ze2
213 REAL(KIND=kind_phys) :: zms_l,zms, z0p5, z0p0, z50s
214 REAL(KIND=kind_phys) :: zgauss(klon), zfluxlaun(klon), zcngl(klon)
215 REAL(KIND=kind_phys) :: zcons1,zcons2,zdelp,zrgpts
216
217!!! try to assign values for the following
218
219 REAL(KIND=kind_phys) :: gssec
220 REAL(KIND=kind_phys) :: zgaussb,zfluxglob
221
222! REAL(KIND=kind_phys) :: PGAW(KLON) don't need this value, set ZDX directly
223
224
225! REAL(KIND=kind_phys) :: PGELAT(KLON) !!! use xlatd from gfs instead
226 REAL(KIND=kind_phys) :: gcoeff, ggaussa !!! GCOEFF link to precip
227! REAL(KIND=kind_phys) :: GGAUSSB !!! GGAUSSB->ZGAUSS
228 !!!! GGAUSSA
229 REAL(KIND=kind_phys) :: gcstar
230
231
232 LOGICAL :: LGACALC, LGSATL, LOZPR
233
234 integer :: NGAUSS, NSLOPE
235
236
237 nslope=1
238 lgacalc=.false.
239 lgsatl=.false.
240 lozpr=.true.
241! NGAUSS=4
242 ngauss=1
243! GGAUSSA=20._kind_phys
244! GGAUSSA=10._kind_phys
245 ggaussa=5._kind_phys
246! GGAUSSB=1.0_kind_phys
247! ZGAUSSB=0.25_kind_phys
248! ZGAUSSB=0.3_kind_phys
249! ZGAUSSB=0.35_kind_phys
250! ZGAUSSB=0.38_kind_phys
251! ZGAUSSB=-0.25_kind_phys
252 zgaussb=0.5_kind_phys
253! ZGAUSSB=0.3_kind_phys
254 gcstar=1.0_kind_phys
255
256
257
258! GSSEC=(pi2/1800.)*(pi2/1800.)
259 gssec=1.e-24
260 gcoeff=1.0_kind_phys !!do not know the value, but never use it when NGAUSS is not equal 1
261
262
263
264!! LOGIC :: LGINDL !!LGINDL=.true. using standard atm values to calculate!! comment out
265!! REAL(KIND=kind_phys) :: STPHI(KM),STTEM(KM), STPREH(KM)
266!! REAL(KIND=kind_phys) :: ZTHSTD,ZRHOSTD,ZBVFSTD
267
268! Set parameters which are a function of launch height
269!!! need to set the parameter ILAUCH, ZFLUXGLOB,ZGAUSSB,ZMS_L directly
270!ILAUNCH=NLAUNCHL(KLAUNCH)
271!ZFLUXGLOB=GFLUXLAUNL(KLAUNCH)
272!ZGAUSSB=GGAUSSB(KLAUNCH)
273!ZMS_L=GMSTAR_L(KLAUNCH)
274
275!*INPUT PARAMETERS
276!* ----------------
277 zradtodeg=57.29577951_kind_phys
278 zms_l=2.e3_kind_phys
279! ZFLUXGLOB=3.75e-3_kind_phys
280! ZFLUXGLOB=3.7e-3_kind_phys
281! ZFLUXGLOB=3.55e-3_kind_phys
282! ZFLUXGLOB=3.65e-3_kind_phys !standard value
283! ZFLUXGLOB=3.62e-3_kind_phys !standard value
284 zfluxglob=3.60e-3_kind_phys
285
286! ZFLUXGLOB=3.2e-3_kind_phys
287! ZFLUXGLOB=3.0e-3_kind_phys
288! ZFLUXGLOB=5.0e-3_kind_phys
289! ZFLUXGLOB=4.0e-3_kind_phys
290! ZFLUXGLOB=0.0_kind_phys
291
292 zms=2*pi/zms_l
293
294!* INITIALIZE FIELDS TO ZERO
295!* -------------------------
296
297 ptenu(:,:)=0.0_kind_phys
298 ptenv(:,:)=0.0_kind_phys
299 pdtdt(:,:)=0.0_kind_phys
300 dked(:,:)=0.0_kind_phys
301
302 DO jk=1,klev+1
303 DO jl=1,klon
304 pfluxu(jl,jk)=0.0_kind_phys
305 pfluxv(jl,jk)=0.0_kind_phys
306 ENDDO
307 ENDDO
308
309
310 DO iazi=1,iazidim
311 DO jk=1,klev
312 DO jl=1,klon
313 zpu(jl,jk,iazi)=0.0_kind_phys
314 zcrt(jl,jk,iazi)=0.0_kind_phys
315 zdfl(jl,jk,iazi)=0.0_kind_phys
316 ENDDO
317 ENDDO
318 ENDDO
319
320!!!!!!!!!!!!!!!!!!!!!!!!
321! redefine ilaunch
322
323 DO jk=1, klev
324 if (papm11(klon,jk) .LT. psrc) exit
325 ENDDO
326 ilaunch = max(jk-1,3)
327
328 DO jl=1,klon
329 zngw(jl) = pgeo11(jl, ilaunch)
330 ENDDO
331
332 ilaunch = klev + 1 - ilaunch
333
334!!!!!!!!!!!!!!!!!!!!!!!!!!!
335
336
337!* reverse vertical coordinate to ECMWF
338 DO jl=1,klon
339 ptm1(jl,:)=transfer(ptm11(jl,klev:1:-1),ptm11(jl,:))
340 pum1(jl,:)=transfer(pum11(jl,klev:1:-1),pum11(jl,:))
341 pvm1(jl,:)=transfer(pvm11(jl,klev:1:-1),pvm11(jl,:))
342 qm(jl,:)=transfer(qm1(jl,klev:1:-1),qm1(jl,:))
343 papm1(jl,:)=transfer(papm11(jl,klev:1:-1),papm11(jl,:))
344 pgeo1(jl,:)=transfer(pgeo11(jl,klev:1:-1),pgeo11(jl,:))
345 prslk(jl,:)=transfer(prslk1(jl,klev:1:-1),prslk1(jl,:))
346
347 dked(jl,:)=transfer(dked(jl,klev:1:-1),dked(jl,:))
348 ptenu(jl,:)=transfer(ptenu(jl,klev:1:-1),ptenu(jl,:))
349 ptenv(jl,:)=transfer(ptenv(jl,klev:1:-1),ptenv(jl,:))
350 pdtdt(jl,:)=transfer(pdtdt(jl,klev:1:-1),pdtdt(jl,:))
351
352
353 paphm1(jl,:)=transfer(paphm11(jl,klev+1:1:-1),paphm11(jl,:))
354 zmeti(jl,:)=transfer(zmeti1(jl,klev+1:1:-1),zmeti1(jl,:))
355 ENDDO
356!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
357
358
359
360
361!* INITIALIZE PARAMETERS FOR COORDINATE TRANSFORM
362!* ----------------------------------------------
363
364! ZCIMIN,ZCIMAX - min,max intrinsic launch-level phase speed (c-U_o) (m/s)
365! ZGAM - half=width of coordinate stretch
366
367 zcimin=0.50_kind_phys
368 zcimax=100.0_kind_phys
369 zgam=0.25_kind_phys
370
371 zpexp=gptwo/2.0_kind_phys
372
373! set initial min ci in each column and azimuth (used for critical levels)
374
375 DO iazi=1,iazidim
376 DO jl=1,klon
377 zci_min(jl,iazi)=zcimin
378 ENDDO
379 ENDDO
380
381
382!* DEFINE HALF MODEL LEVEL WINDS AND TEMPERATURE
383!* -----------------------------------
384
385 DO jk=2,klev
386 DO jl=1,klon
387 zthm1(jl,jk) =0.5_kind_phys*(ptm1(jl,jk-1)+ptm1(jl,jk))
388 zuhm1(jl,jk) =0.5_kind_phys*(pum1(jl,jk-1)+pum1(jl,jk))
389 zvhm1(jl,jk) =0.5_kind_phys*(pvm1(jl,jk-1)+pvm1(jl,jk))
390 ENDDO
391 ENDDO
392
393 jk=1
394 DO jl=1,klon
395 zthm1(jl,jk)=ptm1(jl,jk)
396 zuhm1(jl,jk)=pum1(jl,jk)
397 zvhm1(jl,jk)=pvm1(jl,jk)
398 ENDDO
399
400
401!* DEFINE STATIC STABILITY AND AIR DENSITY ON HALF MODEL LEVELS
402!* ------------------------------------------------------------
403
404! ZCONS1=1.0_kind_phys/RD
405! ZCONS2=RG**2/RCPD
406
407 zcons1=1.0_kind_phys/rd
408 zcons2=grav2cpd
409
410
411 DO jk=klev,2,-1
412 DO jl=1,klon
413! ZDELP=PAPM1(JL,JK)-PAPM1(JL,JK-1)
414 zdelp=(pgeo1(jl,jk)-pgeo1(jl,jk-1))*grav !! times grav since import from zmet
415 zrhohm1(jl,jk)=paphm1(jl,jk)*zcons1/zthm1(jl,jk)
416 zbvfhm1(jl,jk)=zcons2/zthm1(jl,jk)*&
417 & (1.0_kind_phys+cpd*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdelp)
418! & (1.0_kind_phys-RCPD*ZRHOHM1(JL,JK)*(PTM1(JL,JK)-PTM1(JL,JK-1))/ZDELP)
419 zbvfhm1(jl,jk)=max(zbvfhm1(jl,jk),gssec)
420 zbvfhm1(jl,jk)=sqrt(zbvfhm1(jl,jk))
421 ENDDO
422 ENDDO
423
424!* SET UP AZIMUTH DIRECTIONS AND SOME TRIG FACTORS
425!* -----------------------------------------------
426
427 zang=2*pi/iazidim
428 zaz_fct=1.0_kind_phys
429
430
431! get normalization factor to ensure that the same amount of momentum
432! flux is directed (n,s,e,w) no mater how many azimuths are selected.
433! note, however, the code below assumes a symmetric distribution of
434! of azimuthal directions (ie 4,8,16,32,...)
435
436 znorm=0.0_kind_phys
437 DO iazi=1,iazidim
438 zang1=(iazi-1)*zang
439 zcosang(iazi)=cos(zang1)
440 zsinang(iazi)=sin(zang1)
441 znorm=znorm+abs(zcosang(iazi))
442 ENDDO
443 zaz_fct=2._kind_phys*zaz_fct/znorm
444
445
446!* DEFINE COORDINATE TRANSFORM
447!* -----------------------------------------------
448
449! note that this is expresed in terms of the intrinsic phase speed
450! at launch ci=c-u_o so that the transformation is identical at every
451! launch site.
452! See Eq. 28-30 of Scinocca 2003.
453
454 zxmax=1.0_kind_phys/zcimin
455 zxmin=1.0_kind_phys/zcimax
456
457 zxran=zxmax-zxmin
458 zdx=zxran/real(incdim-1)
459 IF(lgacalc) zgam=(zxmax-zxmin)/log(zxmax/zxmin)
460!! LGACALC=.false. ZGAM =0.25 ZX1=0.0007
461!! LGACALC=.true. ZGAM =0.37559 ZX1=0.01
462
463 zx1=zxran/(exp(zxran/zgam)-1.0_kind_phys)
464 zx2=zxmin-zx1
465
466
467
468
469 DO inc=1,incdim
470 ztx=real(inc-1)*zdx+zxmin
471 zx(inc)=zx1*exp((ztx-zxmin)/zgam)+zx2 !Eq. 29 of Scinocca 2003
472 zci(inc)=1.0_kind_phys/zx(inc) !Eq. 28 of Scinocca 2003
473 zdci(inc)=zci(inc)**2*(zx1/zgam)*exp((ztx-zxmin)/zgam)*zdx !Eq. 30 of Scinocca 2003
474 ENDDO
475
476
477!* DEFINE INTRINSIC VELOCITY (RELATIVE TO LAUNCH LEVEL VELOCITY) U(Z)-U(Zo), AND COEFFICINETS
478!* ------------------------------------------------------------------------------------------
479
480 DO iazi=1,iazidim
481 DO jl=1,klon
482 zul(jl,iazi)=zcosang(iazi)*zuhm1(jl,ilaunch)&
483 & +zsinang(iazi)*zvhm1(jl,ilaunch)
484 ENDDO
485 ENDDO
486 DO jl=1,klon
487 zbvfl(jl)=zbvfhm1(jl,ilaunch)
488 ENDDO
489
490 DO jk=2,ilaunch
491 DO iazi=1,iazidim
492 DO jl=1,klon
493 zu=zcosang(iazi)*zuhm1(jl,jk)+zsinang(iazi)*zvhm1(jl,jk)
494 zui(jl,jk,iazi)=zu-zul(jl,iazi)
495 ENDDO
496 ENDDO
497 ENDDO
498
499!* DEFINE RHO(Zo)/N(Zo)
500!* -------------------
501 DO jk=2,ilaunch
502 DO jl=1,klon
503 zfct(jl,jk)=zrhohm1(jl,jk)/zbvfhm1(jl,jk)
504 ENDDO
505 ENDDO
506
507! Optionally set ZFCT at launch level using standard atmos values, to ensure saturation is
508! independent of location
509! IF (LGINDL) THEN
510! ZCONS1=1.0_kind_phys/rd
511! ZCONS2=grav2cpd
512! ZDELP=STPHI(ILAUNCH)-STPHI(ILAUNCH-1) !!probably need to time grav depending on input
513! THSTD=0.5_kind_phys*(STTEM(ILAUNCH-1)+STTEM(ILAUNCH))
514! ZRHOSTD=STPREH(ILAUNCH-1)*ZCONS1/ZTHSTD
515! ZBVFSTD=ZCONS2/ZTHSTD*(1.0_kind_phys+rcpd*
516! & (STTEM(ILAUNCH)-STTEM(ILAUNCH-1))/ZDELP)
517!
518! ZBVFSTD=MAX(ZBVFSTD,GSSEC)
519! ZBVFSTD=SQRT(ZBVFSTD)
520! DO JL=1,KLON
521! ZFCT(JL,ILAUNCH)=ZRHOSTD/ZBVFSTD
522! ENDDO
523! ENDIF
524
525!* SET LAUNCH MOMENTUM FLUX SPECTRAL DENSITY
526!* -----------------------------------------
527
528! Eq. (25) of Scinocca 2003 (not including the 'A' component), and with U-Uo=0
529! do this for only one azimuth since it is identical to all azimuths, and it will be renormalized
530! Initial spectrum fully saturated if LGSATL
531
532 IF(nslope==1) THEN
533! s=1 case
534 DO inc=1,incdim
535 zcin=zci(inc)
536 zcin4=(zms*zcin)**4
537 DO jl=1,klon
538 zbvfl4=zbvfl(jl)**4
539 IF(lgsatl) THEN
540 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl4*min(zcin/zcin4,&
541 & zcin/zbvfl4)
542 ELSE
543 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl4*zcin/(zbvfl4+zcin4)
544 ENDIF
545 zact(jl,inc,1)=1.0_kind_phys
546 ENDDO
547 ENDDO
548
549 ELSEIF(nslope==2) THEN
550! s=2 case
551 DO inc=1,incdim
552 zcin=zci(inc)
553 zcin4=(zms*zcin)**4
554 DO jl=1,klon
555 zbvfl4=zbvfl(jl)**4
556 zcpeak=zbvfl(jl)/zms
557 IF(lgsatl) THEN
558 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl4*min&
559 & (zcpeak/zcin4,zcin/zbvfl4)
560 ELSE
561 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl4*zcin*zcpeak/(zbvfl4&
562 & *zcpeak+zcin4*zcin)
563 ENDIF
564 zact(jl,inc,1)=1.0_kind_phys
565 ENDDO
566 ENDDO
567
568 ELSEIF(nslope==-1) THEN
569! s=-1 case
570 DO inc=1,incdim
571 zcin=zci(inc)
572 zcin2=(zms*zcin)**2
573 DO jl=1,klon
574 zbvfl2=zbvfl(jl)**2
575 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl2*zcin/(zbvfl2+zcin2)
576 zact(jl,inc,1)=1.0_kind_phys
577 ENDDO
578 ENDDO
579
580 ELSEIF(nslope==0) THEN
581! s=0 case
582 DO inc=1,incdim
583 zcin=zci(inc)
584 zcin3=(zms*zcin)**3
585 DO jl=1,klon
586 zbvfl3=zbvfl(jl)**3
587 zflux(jl,inc,1)=zfct(jl,ilaunch)*zbvfl3*zcin/(zbvfl3+zcin3)
588 zact(jl,inc,1)=1.0_kind_phys
589 zacc(jl,inc,1)=1.0_kind_phys
590 ENDDO
591 ENDDO
592
593 ENDIF
594
595!* NORMALIZE LAUNCH MOMENTUM FLUX
596!* ------------------------------
597
598! (rho x F^H = rho_o x F_p^total)
599
600! integrate (ZFLUX x dX)
601 DO inc=1,incdim
602 zcinc=zdci(inc)
603 DO jl=1,klon
604 zpu(jl,ilaunch,1)=zpu(jl,ilaunch,1)+zflux(jl,inc,1)*zcinc
605 ENDDO
606 ENDDO
607
608
609!* NORMALIZE GFLUXLAUN TO INCLUDE SENSITIVITY TO PRECIPITATION
610!* -----------------------------------------------------------
611
612! Also other options to alter tropical values
613
614! A=ZFNORM in Scinocca 2003. A is independent of height.
615 zdxa=1.0_kind_phys/29.e3_kind_phys
616 zdxb=1.0_kind_phys/3.5e3_kind_phys
617 DO jl=1,klon
618!! ZDX=MAX(1.E2_kind_phys,2*RA*SQRT(pi*PGAW(JL))) !grid resolution (m)
619!! ZDX=50.E3 ! c192 for c192 usage
620!!! ZDX=25.E3 !C384
621!!! ZDX=13.E3 !C768
622
623
624 !Scaling factor for launch flux depending on grid resolution
625 ! smooth reduction below 30 km
626 zdxs=1.0_kind_phys-min(1.0_kind_phys,atan((max(1.0_kind_phys&
627 & /dx(jl),zdxa)-zdxa)/(zdxb-zdxa)))
628 zfluxlaun(jl)=zfluxglob*zdxs
629 zfnorm(jl)=zfluxlaun(jl)/zpu(jl,ilaunch,1)
630 ENDDO
631
632! If LOZPR=TRUE then vary EPLAUNCH over tropics
633 IF (lozpr) THEN
634 IF (ngauss==1) THEN
635 z50s=-50.0_kind_phys
636
637 DO jl=1,klon
638
639! ZFLUXLAUN(JL)=ZFLUXLAUN(JL)*(1.0_kind_phys+MIN&
640! & (0.5_kind_phys,GCOEFF*PPRECIP(JL))) !precip
641! ZFNORM(JL)=ZFLUXLAUN(JL)/ZPU(JL,ILAUNCH,1)
642
643 zgelatdeg=xlatd(jl)-z50s
644 zgauss(jl)=zgaussb*exp((-zgelatdeg*zgelatdeg)&
645! & /(2*GGAUSSA*GGAUSSA))
646 & /(2*20.*20.))
647
648! ZGELATDEG=xlatd(JL)
649! ZGAUSS(JL)=-0.1_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)&
650! & /(2*GGAUSSA*GGAUSSA))+ZGAUSS(JL)
651
652
653! ZGELATDEG=xlatd(JL)
654! ZGAUSS(JL)=-0.05_kind_phys*EXP((-ZGELATDEG*ZGELATDEG)&
655! & /(2*GGAUSSA*GGAUSSA))+ZGAUSS(JL)
656
657 zgelatdeg=xlatd(jl)
658 zgauss(jl)=0.08_kind_phys*exp((-zgelatdeg*zgelatdeg)&
659 & /(2*ggaussa*ggaussa))+zgauss(jl)
660
661 zgelatdeg=xlatd(jl)-50.0_kind_phys
662 zgauss(jl)= 0.1_kind_phys*exp((-zgelatdeg*zgelatdeg)&
663 & /(2*10.*10.))+zgauss(jl)
664
665
666 zfluxlaun(jl)=(1.0_kind_phys+zgauss(jl))*zfluxlaun(jl)
667
668
669
670 zfnorm(jl)=zfluxlaun(jl)/zpu(jl,ilaunch,1)
671
672
673 ENDDO
674 ELSEIF (ngauss==2) THEN
675
676 DO jl=1,klon
677! ZGELATDEG=PGELAT(JL)*ZRADTODEG
678 zgelatdeg=xlatd(jl)
679 zgauss(jl)=zgaussb*exp((-zgelatdeg*zgelatdeg)&
680 & /(2*ggaussa*ggaussa))
681 zfluxlaun(jl)=(1.0_kind_phys+zgauss(jl))*zfluxlaun(jl)
682 zfnorm(jl)=zfluxlaun(jl)/zpu(jl,ilaunch,1)
683 ENDDO
684
685
686
687 ELSEIF (ngauss==4) THEN
688
689! Set latitudinal dependence to optimize stratospheric winds for 36r1
690 z50s=-50.0_kind_phys
691 DO jl=1,klon
692! ZGELATDEG=PGELAT(JL)*ZRADTODEG-Z50S
693 zgelatdeg=xlatd(jl)-z50s
694 zgauss(jl)=zgaussb*exp((-zgelatdeg*zgelatdeg)/(2*ggaussa*ggaussa))
695 zfluxlaun(jl)=(1.0_kind_phys+zgauss(jl))*zfluxlaun(jl)
696 zfnorm(jl)=zfluxlaun(jl)/zpu(jl,ilaunch,1)
697 ENDDO
698
699 ENDIF
700 ENDIF
701
702 DO iazi=1,iazidim
703 DO jl=1,klon
704 zpu(jl,ilaunch,iazi)=zfluxlaun(jl)
705 ENDDO
706 ENDDO
707
708!* ADJUST CONSTANT ZFCT
709!* --------------------
710 DO jk=2,ilaunch
711 DO jl=1,klon
712 zfct(jl,jk)=zfnorm(jl)*zfct(jl,jk)
713 ENDDO
714 ENDDO
715
716!* RENORMALIZE EACH SPECTRAL ELEMENT IN FIRST AZIMUTH
717!* --------------------------------------------------
718 DO inc=1,incdim
719 DO jl=1,klon
720 zflux(jl,inc,1)=zfnorm(jl)*zflux(jl,inc,1)
721 ENDDO
722 ENDDO
723
724
725
726!* COPY ZFLUX INTO ALL OTHER AZIMUTHS
727!* --------------------------------
728
729! ZACT=1 then no critical level
730! ZACT=0 then critical level
731
732 DO iazi=2,iazidim
733 DO inc=1,incdim
734 DO jl=1,klon
735 zflux(jl,inc,iazi)=zflux(jl,inc,1)
736 zact(jl,inc,iazi)=1.0_kind_phys
737 zacc(jl,inc,iazi)=1.0_kind_phys
738 ENDDO
739 ENDDO
740 ENDDO
741
742! -----------------------------------------------------------------------------
743
744!* BEGIN MAIN LOOP OVER LEVELS
745!* ---------------------------
746
747!* begin IAZIDIM do-loop
748!* --------------------
749
750 DO iazi=1,iazidim
751
752!* begin JK do-loop
753!* ----------------
754
755 DO jk=ilaunch-1,2,-1
756
757!* first do critical levels
758!* ------------------------
759
760 DO jl=1,klon
761 zci_min(jl,iazi)=max(zci_min(jl,iazi),zui(jl,jk,iazi))
762 ENDDO
763
764!* set ZACT to zero if critical level encountered
765!* ----------------------------------------------
766
767 z0p5=0.5_kind_phys
768 DO inc=1,incdim
769 zcin=zci(inc)
770 DO jl=1,klon
771 zatmp=z0p5+sign(z0p5,zcin-zci_min(jl,iazi))
772 zacc(jl,inc,iazi)=zact(jl,inc,iazi)-zatmp
773 zact(jl,inc,iazi)=zatmp
774 ENDDO
775 ENDDO
776
777!* integrate to get critical-level contribution to mom deposition on this level, i.e. ZACC=1
778!* ----------------------------------------------------------------------------------------
779
780 DO inc=1,incdim
781 zcinc=zdci(inc)
782 DO jl=1,klon
783 zdfl(jl,jk,iazi)=zdfl(jl,jk,iazi)+&
784 & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zcinc
785 ENDDO
786 ENDDO
787
788!* get weighted average of phase speed in layer
789
790 DO jl=1,klon
791 IF(zdfl(jl,jk,iazi)>0.0_kind_phys) THEN
792 zatmp=zcrt(jl,jk,iazi)
793 DO inc=1,incdim
794 zatmp=zatmp+zci(inc)*&
795 & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zdci(inc)
796 ENDDO
797 zcrt(jl,jk,iazi)=zatmp/zdfl(jl,jk,iazi)
798 ELSE
799 zcrt(jl,jk,iazi)=zcrt(jl,jk+1,iazi)
800 ENDIF
801 ENDDO
802
803!* do saturation (Eq. (26) and (27) of Scinocca 2003)
804!* -------------------------------------------------
805
806 IF(gptwo==3.0_kind_phys) THEN
807 DO inc=1,incdim
808 zcin=zci(inc)
809 zcinc=1.0_kind_phys/zcin
810 DO jl=1,klon
811 ze1=zcin-zui(jl,jk,iazi)
812 ze2=gcstar*zfct(jl,jk)*ze1
813 zfluxsq=ze2*ze2*ze1*zcinc
814 ! ZFLUXSQ=ZE2*ZE2*ZE1/ZCIN
815 zdep=zact(jl,inc,iazi)*(zflux(jl,inc,iazi)**2-zfluxsq)
816 IF(zdep>0.0_kind_phys) THEN
817 zflux(jl,inc,iazi)=sqrt(zfluxsq)
818 ENDIF
819 ENDDO
820 ENDDO
821 ELSEIF(gptwo==2.0_kind_phys) THEN
822 DO inc=1,incdim
823 zcin=zci(inc)
824 zcinc=1.0_kind_phys/zcin
825 DO jl=1,klon
826 zfluxs=gcstar*zfct(jl,jk)*&
827 & (zcin-zui(jl,jk,iazi))**2*zcinc
828 ! ZFLUXS=GCSTAR*ZFCT(JL,JK)*(ZCIN-ZUI(JL,JK,IAZI))**2/ZCIN
829 zdep=zact(jl,inc,iazi)*(zflux(jl,inc,iazi)-zfluxs)
830 IF(zdep>0.0_kind_phys) THEN
831 zflux(jl,inc,iazi)=zfluxs
832 ENDIF
833 ENDDO
834 ENDDO
835 ENDIF
836
837!* integrate spectrum
838!* ------------------
839
840 DO inc=1,incdim
841 zcinc=zdci(inc)
842 DO jl=1,klon
843 zpu(jl,jk,iazi)=zpu(jl,jk,iazi)+&
844 & zact(jl,inc,iazi)*zflux(jl,inc,iazi)*zcinc
845 ENDDO
846 ENDDO
847
848!* end JK do-loop
849!* --------------
850
851 ENDDO
852!* end IAZIDIM do-loop
853!* ---------------
854
855 ENDDO
856
857! -----------------------------------------------------------------------------
858
859!* MAKE CORRECTION FOR CRITICAL-LEVEL MOMENTUM DEPOSITION
860!* ------------------------------------------------------
861
862 z0p0=0._kind_phys
863! ZRGPTS=1.0_kind_phys/(RG*PTSTEP)
864 zrgpts=1.0_kind_phys/(grav*ptstep)
865 DO iazi=1,iazidim
866 DO jl=1,klon
867 zcngl(jl)=0.0_kind_phys
868 ENDDO
869 DO jk=2,ilaunch
870 DO jl=1,klon
871 zulm=zcosang(iazi)*pum1(jl,jk)+zsinang(iazi)*&
872 & pvm1(jl,jk)-zul(jl,iazi)
873 zdfl(jl,jk-1,iazi)=zdfl(jl,jk-1,iazi)+zcngl(jl)
874 zdft=min(zdfl(jl,jk-1,iazi),2.0_kind_phys*(papm1(jl,jk-1)-&
875 & papm1(jl,jk))*(zcrt(jl,jk-1,iazi)-zulm)*zrgpts)
876
877 zdft=max(zdft,z0p0)
878 zcngl(jl)=(zdfl(jl,jk-1,iazi)-zdft)
879 zpu(jl,jk,iazi)=zpu(jl,jk,iazi)-zcngl(jl)
880 ENDDO
881 ENDDO
882 ENDDO
883
884
885!* SUM CONTRIBUTION FOR TOTAL ZONAL AND MERIDIONAL FLUX
886!* ---------------------------------------------------
887
888 DO iazi=1,iazidim
889 DO jk=ilaunch,2,-1
890 DO jl=1,klon
891 pfluxu(jl,jk)=pfluxu(jl,jk)+zpu(jl,jk,iazi)*zaz_fct*zcosang(iazi)
892 pfluxv(jl,jk)=pfluxv(jl,jk)+zpu(jl,jk,iazi)*zaz_fct*zsinang(iazi)
893 ENDDO
894 ENDDO
895 ENDDO
896
897
898!* UPDATE U AND V TENDENCIES
899!* ----------------------------
900
901! ZCONS1=1.0_kind_phys/RCPD
902 zcons1=rcpd
903 DO jk=1,ilaunch
904 DO jl=1, klon
905! ZDELP= RG/(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
906 zdelp= grav/(paphm1(jl,jk+1)-paphm1(jl,jk))
907 ze1=(pfluxu(jl,jk+1)-pfluxu(jl,jk))*zdelp
908 ze2=(pfluxv(jl,jk+1)-pfluxv(jl,jk))*zdelp
909
910 if (abs(ze1) >= maxdudt ) then
911 ze1 = sign(maxdudt, ze1)
912 endif
913 if (abs(ze2) >= maxdudt ) then
914 ze2 = sign(maxdudt, ze2)
915 endif
916
917
918 ptenu(jl,jk)=ze1
919 ptenv(jl,jk)=ze2
920! add the tendency of dT/dt
921 ze2=-(pum1(jl,jk)*ptenu(jl,jk)+pvm1(jl,jk)*ptenv(jl,jk))/cpd
922 if (abs(ze2) >= max_eps) pdtdt(jl,jk) = sign(max_eps, ze2)
923
924! end of the tendency of dT/dt
925 ENDDO
926 ENDDO
927
928
929!* reverse vertical coordinate back to GFS for outbound variables only
930 DO jl=1,klon
931 dked(jl,:)=transfer(dked(jl,klev:1:-1),dked(jl,:))
932 ptenu(jl,:)=transfer(ptenu(jl,klev:1:-1),ptenu(jl,:))
933 ptenv(jl,:)=transfer(ptenv(jl,klev:1:-1),ptenv(jl,:))
934 pdtdt(jl,:)=transfer(pdtdt(jl,klev:1:-1),pdtdt(jl,:))
935 ENDDO
936
937
938! DO JL=1,KLON
939! dked(JL,:)=0.0
940! PTENU(JL,:)=0.0
941! PTENV(JL,:)=0.0
942! pdtdt(JL,:)=0.0
943! ENDDO
944
945!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
946
947
948
949 return
950 end subroutine ecmwf_ngw_emc
951
952
953end module ecmwf_ngw
954
955
956
957
ecmwf_ngw
Definition ecmwf_ngw.F90:4