CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cires_ugwpv1_initialize.F90
1
3
4!===============================
5! cu-cires ugwp-scheme
6! initialization of selected
7! init gw-solvers (1,2,3,4)
8! init gw-source specifications
9! init gw-background dissipation
10!==============================
11!
12! Part-0 specifications of common constants, limiters and "criiical" values
13!
14!
15
17!
18 use machine, only : kind_phys
19
20 implicit none
21
22 real(kind=kind_phys) :: pi, pi2, pih, rad_to_deg, deg_to_rad
23 real(kind=kind_phys) :: arad, p0s
24 real(kind=kind_phys) :: grav, grav2, rgrav, rgrav2
25 real(kind=kind_phys) :: cpd, rd, rv, fv
26 real(kind=kind_phys) :: rdi, rcpd, rcpd2
27
28 real(kind=kind_phys) :: gor, gr2, grcp, gocp, rcpdl, grav2cpd
29 real(kind=kind_phys) :: bnv2min, bnv2max
30 real(kind=kind_phys) :: dw2min, velmin, minvel
31 real(kind=kind_phys) :: omega1, omega2, omega3
32 real(kind=kind_phys) :: hpscale, rhp, rhp2, rh4, rhp4, khp, hpskm
33 real(kind=kind_phys) :: mkzmin, mkz2min, mkzmax, mkz2max, cdmin
34 real(kind=kind_phys) :: rcpdt
35
36! real(kind=kind_phys), parameter :: grav2 = grav + grav
37! real(kind=kind_phys), parameter :: rgrav = 1.0/grav, rgrav2= rgrav*rgrav
38! real(kind=kind_phys), parameter :: rdi = 1.0 / rd, rcpd = 1./cpd, rcpd2 = 0.5/cpd
39! real(kind=kind_phys), parameter :: gor = grav/rd, rcpdt = 1./(cp*dtp)
40
41! real(kind=kind_phys), parameter :: gr2 = grav*gor
42! real(kind=kind_phys), parameter :: grcp = grav*rcpd, gocp = grcp
43! real(kind=kind_phys), parameter :: rcpdl = cpd*rgrav ! 1/[g/cp] == cp/g
44! real(kind=kind_phys), parameter :: grav2cpd = grav*grcp ! g*(g/cp)= g^2/cp
45! real(kind=kind_phys), parameter :: pi2 = 2.*pi, pih = .5*pi
46! real(kind=kind_phys), parameter :: rad_to_deg=180.0/pi, deg_to_rad=pi/180.0
47!
48! real(kind=kind_phys), parameter :: bnv2min = (pi2/1800.)*(pi2/1800.)
49! real(kind=kind_phys), parameter :: bnv2max = (pi2/30.)*(pi2/30.)
50! real(kind=kind_phys), parameter :: dw2min=1.0, velmin=sqrt(dw2min), minvel = 0.5
51! real(kind=kind_phys), parameter :: omega1 = pi2/86400., omega2 = 2.*omega1, omega3 = 3.*omega1
52!
53! real(kind=kind_phys), parameter :: hpscale= 7000., rhp=1./hpscale, rhp2=.5*rhp, rh4 = 0.25*rhp
54! real(kind=kind_phys), parameter :: mkzmin = pi2/80.0e3, mkz2min = mkzmin*mkzmin
55! real(kind=kind_phys), parameter :: mkzmax = pi2/500., mkz2max = mkzmax*mkzmax
56! real(kind=kind_phys), parameter :: cdmin = 2.e-2/mkzmax
57! real(kind=kind_phys), parameter :: pi = 4.*atan(1.0),
58! real(kind=kind_phys), parameter :: grav =9.81, cpd = 1004.
59! real(kind=kind_phys), parameter :: rd = 287.0 , rv =461.5
60! real(kind=kind_phys), parameter :: fv = rv/rd - 1.0
61! real(kind=kind_phys), parameter :: arad = 6370.e3
62
63 end module ugwp_common
64
65 subroutine init_nazdir(naz, xaz, yaz)
66
67 use machine, only : kind_phys
68 use ugwp_common, only : pi2
69
70 implicit none
71
72 integer :: naz
73 real(kind=kind_phys), dimension(naz) :: xaz, yaz
74 integer :: idir
75 real(kind=kind_phys) :: phic, drad
76
77 drad = pi2/float(naz)
78 if (naz.ne.4) then
79 do idir =1, naz
80 phic = drad*(float(idir)-1.0)
81 xaz(idir) = cos(phic)
82 yaz(idir) = sin(phic)
83 enddo
84 else
85! if (naz.eq.4) then
86 xaz(1) = 1.0 !E
87 yaz(1) = 0.0
88 xaz(2) = 0.0
89 yaz(2) = 1.0 !N
90 xaz(3) =-1.0 !W
91 yaz(3) = 0.0
92 xaz(4) = 0.0
93 yaz(4) =-1.0 !S
94 endif
95 end subroutine init_nazdir
96!
97!
98!===================================================
99!
100!Part-1 init => wave dissipation + RFriction
101!
102!===================================================
103 subroutine init_global_gwdis(levs, zkm, pmb, kvg, ktg, krad, kion, me, master)
104!
105 use machine , only : kind_phys
106 use ugwp_common, only : pih, pi
107
108 implicit none
109 integer , intent(in) :: me, master
110 integer , intent(in) :: levs
111 real(kind=kind_phys), intent(in) :: zkm(levs), pmb(levs) ! in km-Pa
112 real(kind=kind_phys), intent(out), dimension(levs+1) :: kvg, ktg, krad, kion
113!
114!locals + data
115!
116 integer :: k
117 real(kind=kind_phys), parameter :: vusurf = 2.e-5
118 real(kind=kind_phys), parameter :: musurf = vusurf/1.95
119 real(kind=kind_phys), parameter :: hpmol = 7.0
120!
121 real(kind=kind_phys), parameter :: kzmin = 0.1
122 real(kind=kind_phys), parameter :: kturbo = 100.
123 real(kind=kind_phys), parameter :: zturbo = 130.
124 real(kind=kind_phys), parameter :: zturw = 30.
125 real(kind=kind_phys), parameter :: inv_pra = 3. !kt/kv =inv_pr
126!
127 real(kind=kind_phys), parameter :: alpha = 1./86400./15. ! height variable see Zhu-1993 from 60-days => 6 days
128 real(kind=kind_phys) :: pa_alp = 750. ! super-RF parameters from FV3-dycore GFSv17/16 sett
129 real(kind=kind_phys) :: tau_alp = 10. ! days (750 Pa /10days)
130!
131 real(kind=kind_phys), parameter :: kdrag = 1./86400./30. !parametrization for WAM ion drag as e-density function
132 real(kind=kind_phys), parameter :: zdrag = 100.
133 real(kind=kind_phys), parameter :: zgrow = 50.
134!
135 real(kind=kind_phys) :: vumol, mumol, keddy, ion_drag
136 real(kind=kind_phys) :: rf_fv3, rtau_fv3, ptop, pih_dlog
137!
138 real(kind=kind_phys) :: ae1 ,ae2
139!
140
141 ptop = pmb(levs)
142 rtau_fv3 = 1./86400./tau_alp
143 pih_dlog = pih/log(pa_alp/ptop)
144
145 do k=1, levs
146 ae1 = zkm(k)/hpmol
147 vumol = vusurf*exp(ae1)
148 mumol = musurf*exp(ae1)
149 ae2 = -((zkm(k)-zturbo) /zturw)**2
150 keddy = kturbo*exp(ae2)
151
152 kvg(k) = vumol + keddy
153 ktg(k) = mumol + keddy*inv_pra
154
155 krad(k) = alpha
156!
157 ion_drag = kdrag
158!
159 kion(k) = ion_drag!
160! add Rayleigh_Super of FV3 for pmb < pa_alp
161!
162 if (pmb(k) .le. pa_alp) then
163 rf_fv3=rtau_fv3*sin(pih_dlog*log(pa_alp/pmb(k)))**2
164 krad(k) = krad(k) + rf_fv3
165 kion(k) = kion(k) + rf_fv3
166
167 endif
168
169! write(6,132) zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2, kion(k)
170 enddo
171
172 k= levs+1
173 kion(k) = kion(k-1)
174 krad(k) = krad(k-1)
175 kvg(k) = kvg(k-1)
176 ktg(k) = ktg(k-1)
177
178! if (me == master) then
179! write(6, * ) ' zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2, kion(k) ... init_global_gwdis'
180! do k=1, levs, 1
181! write(6,132) zkm(k), kvg(k), kvg(k)*(6.28/5000.)**2, kion(k), pmb(k)
182! enddo
183! endif
184!
185! 132 format( 2x, F8.3,' dis-scales:', 4(2x, E10.3))
186
187 end subroutine init_global_gwdis
188!
189! ========================================================================
190! Part 2 - sources
191! wave sources
192! ========================================================================
193!
194! ugwp_oro_init
195!
196!=========================================================================
198 use machine , only : kind_phys
199 use ugwp_common, only : bnv2min, grav, grcp, fv, grav, cpd, grcp, pi
200 use ugwp_common, only : mkzmin, mkz2min
201 implicit none
202!
203! constants and "crirtical" values to run oro-mtb_gw physics
204!
205!
206 real(kind=kind_phys), parameter :: hncrit=9000. ! max value in meters for elvmax
207 real(kind=kind_phys), parameter :: hminmt=50. ! min mtn height (*j*)
208 real(kind=kind_phys), parameter :: sigfac=4.0 ! mb3a expt test for elvmax factor
209 real(kind=kind_phys), parameter :: hpmax=2500.0
210 real(kind=kind_phys), parameter :: hpmin=25.0
211!
212!
213 real(kind=kind_phys), parameter :: minwnd=1.0 ! min wind component (*j*)
214 real(kind=kind_phys), parameter :: dpmin=5000.0 ! minimum thickness of the reference layer in pa
215
216
217 character(len=8) :: strver = 'gfs_2018'
218 character(len=8) :: strbase = 'gfs_2018'
219
220 real(kind=kind_phys), parameter :: rimin=-10., ric=0.25
221
222 real(kind=kind_phys), parameter :: frmax=10., frc =1.0, frmin =0.01
223 real(kind=kind_phys), parameter :: ce=0.8, ceofrc=ce/frc, cg=0.5
224 real(kind=kind_phys), parameter :: gmax=1.0, veleps=1.0, factop=0.5!
225 real(kind=kind_phys), parameter :: efmin=0.5, efmax=10.0
226
227 real(kind=kind_phys), parameter :: rlolev=50000.0
228 integer, parameter :: mdir = 8
229 real(kind=kind_phys), parameter :: fdir=mdir/(8.*atan(1.0))
230 real(kind=kind_phys), parameter :: zpgeo=2.*atan(1.0)
231
232 integer nwdir(mdir)
233 data nwdir/6,7,5,8,2,3,1,4/
234 save nwdir
235
236 real(kind=kind_phys), parameter :: odmin = 0.1, odmax = 10.0
237 real(kind=kind_phys), parameter :: fcrit_sm = 0.7, fcrit_sm2 = fcrit_sm * fcrit_sm
238 real(kind=kind_phys), parameter :: fcrit_gfs = 0.7, fcrit_v1 = 0.7
239 real(kind=kind_phys), parameter :: fcrit_mtb = 0.7
240
241 real(kind=kind_phys), parameter :: zbr_pi = zpgeo
242 real(kind=kind_phys), parameter :: zbr_ifs = zpgeo
243
244!
245
246 real(kind=kind_phys), parameter :: kxoro=6.28e-3/200. !
247 real(kind=kind_phys), parameter :: coro = 0.0
248 integer,parameter :: nridge=2
249 real(kind=kind_phys), parameter :: sigma_std=1./100., gamm_std=1.0
250
251 real(kind=kind_phys) :: cdmb ! scale factors for mtb
252 real(kind=kind_phys) :: cleff ! scale factors for orogw
253
254 integer :: nworo ! number of waves
255 integer :: nazoro ! number of azimuths
256 integer :: nstoro ! flag for stochastic launch above SG-peak
257
258
259!------------------------------------------------------------------------------
260! small-scale orography parameters for TOFD of Beljaars et al., 2004, QJRMS
261! SA-option can be controlled by Integral limits of fluxes
262! in B2004: klow = 0.003 1/m ~ 2km and kinf ~ 6.28/10/(Z1)~< 1 km => meters
263! these limits can change strength of TOFD... choice of k0tr ~1/10 km (10km ~dx of C768)
264! kmax = kdis_pbl
265!------------------------------------------------------------------------------
266 real(kind=kind_phys), parameter :: kmax = 6.28/(10.*25.) ! max k-tofd
267 real(kind=kind_phys), parameter :: k1tr = 6.28/(2100) ! max k-transition from -1.9/slope to -2.8/slope
268 real(kind=kind_phys), parameter :: kflt = 6.28/(18.e3) !
269 real(kind=kind_phys), parameter :: k0tr = 6.28/(10.e3) ! min k-tofd
270 real(kind=kind_phys), parameter :: nk1tr = 2.8
271 real(kind=kind_phys), parameter :: nk0tr = 1.9
272 real(kind=kind_phys), parameter :: a1_tofd = kflt ** nk1tr *1.e3
273 real(kind=kind_phys), parameter :: a2_tofd = k1tr ** (nk0tr-nk1tr)
274 real(kind=kind_phys), parameter :: fix_tofd = 2.* 0.005 * 12 *0.6 !value= 0.072
275!
276! B2004 scheme is based on the empirical vertical profile of the tofd divergence:
277! Ax_tofd(Z)=exp(-[Z/ze_tofd]^3/2) / Z^1.2.....
278! TOFD-flux/TMS-flux must dissipate due to PBL-diffusion with spectral damping
279! Here we can enhance TOFD-impact by selecting k0tr and kmax limits
280! as functions of resolution and PBL-dissipation
281!
282 integer, parameter :: n_tofd = 2 ! depth of SSO for TOFD compared with Zpbl
283 real(kind=kind_phys), parameter :: const_tofd = 0.0759 ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
284 real(kind=kind_phys), parameter :: ze_tofd = 1500.0 ! BJ's z-decay in meters, 1.5 km
285 real(kind=kind_phys), parameter :: a12_tofd = 0.0002662*0.005363 ! BJ's k-spect const for sigf2 * a1*a2*exp(-[z/zdec]**1.5]
286 real(kind=kind_phys), parameter :: ztop_tofd = 3.*ze_tofd ! no TOFD > this height 4.5 km
287!------------------------------------------------------------------------------
288!
289
290 contains
291!
292 subroutine init_oro_gws(nwaves, nazdir, nstoch, effac, &
293 lonr, kxw )
294!
295!
296 integer :: nwaves, nazdir, nstoch
297 integer :: lonr
298 real(kind=kind_phys) :: cdmbx
299 real(kind=kind_phys) :: kxw
300 real(kind=kind_phys) :: effac ! it is analog of cdmbgwd(2) for GWs, off for now
301!-----------------------------! GFS-setup for cdmb & cleff
302! cdmb = 4.0 * (192.0/IMX)
303! cleff = 0.5E-5 / SQRT(IMX/192.0) = 0.5E-5*SQRT(192./IMX)
304!
305 real(kind=kind_phys), parameter :: lonr_refmb = 4.0 * 192.0
306 real(kind=kind_phys), parameter :: lonr_refgw = 192.0
307 real(kind=kind_phys), parameter :: cleff_ref = 0.5e-5 ! 1256 km = 10 * 125 km ???
308
309! copy to "ugwp_oro_init" => nwaves, nazdir, nstoch
310
311 nworo = nwaves
312 nazoro = nazdir
313 nstoro = nstoch
314
315 cdmbx = lonr_refmb/float(lonr)
316
317 cdmb = cdmbx
318 cleff = cleff_ref * sqrt(lonr_refgw/float(lonr)) !* effac
319!
320 end subroutine init_oro_gws
321!
322
323 end module ugwp_oro_init
324! =========================================================================
325!
326! ugwp_conv_init
327!
328!=========================================================================
330
331 use machine , only : kind_phys
332
333
334 implicit none
335 real(kind=kind_phys) :: eff_con ! scale factors for conv GWs
336 integer :: nwcon ! number of waves
337 integer :: nazcon ! number of azimuths
338 integer :: nstcon ! flag for stochastic choice of launch level above Conv-cloud
339 real(kind=kind_phys) :: con_dlength
340 real(kind=kind_phys) :: con_cldf
341
342 real(kind=kind_phys), parameter :: cmin = 5 !2.5
343 real(kind=kind_phys), parameter :: cmax = 95. !82.5
344 real(kind=kind_phys), parameter :: cmid = 22.5
345 real(kind=kind_phys), parameter :: cwid = cmid
346 real(kind=kind_phys), parameter :: bns = 2.e-2, bns2 = bns*bns, bns4=bns2*bns2
347 real(kind=kind_phys), parameter :: mstar = 6.28e-3/2. ! 2km
348 real(kind=kind_phys) :: dc
349
350 real(kind=kind_phys), allocatable :: ch_conv(:), spf_conv(:)
351 real(kind=kind_phys), allocatable :: xaz_conv(:), yaz_conv(:)
352 contains
353!
354 subroutine init_conv_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
355!
356 use ugwp_common, only : pi2, arad
357
358 implicit none
359
360
361 integer :: nwaves, nazdir, nstoch
362 integer :: lonr
363!
364! ccpp
365!
366
367 real(kind=kind_phys) :: kxw, effac
368 real(kind=kind_phys) :: work1 = 0.5
369 real(kind=kind_phys) :: chk, tn4, snorm
370 integer :: k
371
372 nwcon = nwaves
373 nazcon = nazdir
374 nstcon = nstoch
375 eff_con = effac
376
377 con_dlength = pi2*arad/float(lonr)
378!
379! allocate & define spectra in "selected direction": "dc" "ch(nwaves)"
380!
381 if (.not. allocated(ch_conv)) allocate (ch_conv(nwaves))
382 if (.not. allocated(spf_conv)) allocate (spf_conv(nwaves))
383 if (.not. allocated(xaz_conv)) allocate (xaz_conv(nazdir))
384 if (.not. allocated(yaz_conv)) allocate (yaz_conv(nazdir))
385
386 dc = (cmax-cmin)/float(nwaves-1)
387!
388! we may use different spectral "shapes"
389! for example FVS-93 "Desabeius"
390! E(s=1, t=3,m, w, k) ~ m^s/(m*^4 + m^4) ~ m^-3 saturated tail
391!
392 do k = 1,nwaves
393 chk = cmin + (k-1)*dc
394 tn4 = (mstar*chk)**4
395 ch_conv(k) = chk
396 spf_conv(k) = bns4*chk/(bns4+tn4)
397 enddo
398
399 snorm = sum(spf_conv)
400 spf_conv = spf_conv/snorm*1.5
401
402 call init_nazdir(nazdir, xaz_conv, yaz_conv)
403 end subroutine init_conv_gws
404
405
406 end module ugwp_conv_init
407!=========================================================================
408!
409! ugwp_fjet_init
410!
411!=========================================================================
412
414 use machine , only : kind_phys
415
416
417
418 implicit none
419 real(kind=kind_phys) :: eff_fj ! scale factors for conv GWs
420 integer :: nwfj ! number of waves
421 integer :: nazfj ! number of azimuths
422 integer :: nstfj ! flag for stochastic choice of launch level above Conv-cloud
423!
424 real(kind=kind_phys), parameter :: fjet_trig=0. ! if ( abs(frgf) > fjet_trig ) launch GW-packet
425
426
427 real(kind=kind_phys), parameter :: cmin = 2.5
428 real(kind=kind_phys), parameter :: cmax = 67.5
429 real(kind=kind_phys) :: dc
430 real(kind=kind_phys), allocatable :: ch_fjet(:) , spf_fjet(:)
431 real(kind=kind_phys), allocatable :: xaz_fjet(:), yaz_fjet(:)
432 contains
433
434 subroutine init_fjet_gws(nwaves, nazdir, nstoch, effac,lonr, kxw)
435
436 use ugwp_common, only : pi2, arad
437
438 implicit none
439
440 integer :: nwaves, nazdir, nstoch
441 integer :: lonr
442 real(kind=kind_phys) :: kxw, effac , chk
443
444 integer :: k
445
446 nwfj = nwaves
447 nazfj = nazdir
448 nstfj = nstoch
449 eff_fj = effac
450
451 if (.not. allocated(ch_fjet)) allocate (ch_fjet(nwaves))
452 if (.not. allocated(spf_fjet)) allocate (spf_fjet(nwaves))
453 if (.not. allocated(xaz_fjet)) allocate (xaz_fjet(nazdir))
454 if (.not. allocated(yaz_fjet)) allocate (yaz_fjet(nazdir))
455
456 dc = (cmax-cmin)/float(nwaves-1)
457 do k = 1,nwaves
458 chk = cmin + (k-1)*dc
459 ch_fjet(k) = chk
460 spf_fjet(k) = 1.0
461 enddo
462 call init_nazdir(nazdir, xaz_fjet, yaz_fjet)
463
464 end subroutine init_fjet_gws
465
466 end module ugwp_fjet_init
467!
468!=========================================================================
469!
470!
472!=========================================================================
473 use machine , only : kind_phys
474
475 implicit none
476
477 real(kind=kind_phys) :: eff_okw ! scale factors for conv GWs
478 integer :: nwokw ! number of waves
479 integer :: nazokw ! number of azimuths
480 integer :: nstokw ! flag for stochastic choice of launch level above Conv-cloud
481!
482 real(kind=kind_phys), parameter :: okw_trig=0. ! if ( abs(okwp) > okw_trig ) launch GW-packet
483
484 real(kind=kind_phys), parameter :: cmin = 2.5
485 real(kind=kind_phys), parameter :: cmax = 67.5
486 real(kind=kind_phys) :: dc
487 real(kind=kind_phys), allocatable :: ch_okwp(:), spf_okwp(:)
488 real(kind=kind_phys), allocatable :: xaz_okwp(:), yaz_okwp(:)
489
490 contains
491!
492
493
494 subroutine init_okw_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
495 use ugwp_common, only : pi2, arad
496
497 implicit none
498
499 integer :: nwaves, nazdir, nstoch
500 integer :: lonr
501 real(kind=kind_phys) :: kxw, effac , chk
502
503 integer :: k
504
505 nwokw = nwaves
506 nazokw = nazdir
507 nstokw = nstoch
508 eff_okw = effac
509
510 if (.not. allocated(ch_okwp)) allocate (ch_okwp(nwaves))
511 if (.not. allocated(spf_okwp)) allocate (spf_okwp(nwaves))
512 if (.not. allocated(xaz_okwp)) allocate (xaz_okwp(nazdir))
513 if (.not. allocated(yaz_okwp)) allocate (yaz_okwp(nazdir))
514 dc = (cmax-cmin)/float(nwaves-1)
515 do k = 1,nwaves
516 chk = cmin + (k-1)*dc
517 ch_okwp(k) = chk
518 spf_okwp(k) = 1.
519 enddo
520
521 call init_nazdir(nazdir, xaz_okwp, yaz_okwp)
522!
523 end subroutine init_okw_gws
524
525 end module ugwp_okw_init
526
527!=============================== end of GW sources
528!
529! init specific gw-solvers (1,2,3,4)
530!
531!===============================
532! Part -3 init wave solvers
533!===============================
534
536 use machine , only : kind_phys
537 implicit none
538
539 integer :: nwav, nazd
540 integer :: nst
541 real(kind=kind_phys) :: eff
542 integer, parameter :: incdim = 4, iazdim = 4
543!
544 contains
545
546 subroutine initsolv_lsatdis(me, master, nwaves, nazdir, nstoch, effac, do_physb, kxw)
547
548 implicit none
549!
550 integer :: me, master
551 integer :: nwaves, nazdir
552 integer :: nstoch
553 real(kind=kind_phys) :: effac
554 logical :: do_physb
555 real(kind=kind_phys) :: kxw
556!
557!locals: define azimuths and Ch(nwaves) - domain when physics-based soureces
558! are not actibve
559!
560 integer :: inc, jk, jl, iazi, i, j, k
561
562 if( nwaves == 0 .or. nstoch == 1 ) then
563! redefine from the default
564 nwav = incdim
565 nazd = iazdim
566 nst = 0
567 eff = 1.0
568 else
569! from input_nml multi-wave spectra
570 nwav = nwaves
571 nazd = nazdir
572 nst = nstoch
573 eff = effac
574 endif
575!
576 end subroutine initsolv_lsatdis
577!
578 end module ugwp_lsatdis_init
579!
580!
582
583 use machine , only : kind_phys
584 use ugwp_common, only : arad, pi, pi2, hpscale, rhp, rhp2, rh4, omega2
585 use ugwp_common, only : bnv2max, bnv2min, minvel
586 use ugwp_common, only : mkzmin, mkz2min, mkzmax, mkz2max, ucrit => cdmin
587
588 implicit none
589
590 real(kind=kind_phys), parameter :: maxdudt = 250.e-5, maxdtdt=15.e-2
591 real(kind=kind_phys), parameter :: dked_min =0.01, dked_max=250.0
592
593 real(kind=kind_phys), parameter :: gptwo=2.0
594
595 real(kind=kind_phys) , parameter :: bnfix = 6.28/300., bnfix2= bnfix * bnfix
596 real(kind=kind_phys) , parameter :: bnfix4 = bnfix2 * bnfix2
597 real(kind=kind_phys) , parameter :: bnfix3 = bnfix2 * bnfix
598!
599! make parameter list that will be passed to SOLVER
600!
601 integer , parameter :: iazidim=4 ! number of azimuths
602 integer , parameter :: incdim=25 ! number of discrete cx - spectral elements in launch spectrum
603
604 real(kind=kind_phys) , parameter :: zcimin = 2.5
605 real(kind=kind_phys) , parameter :: zcimax = 125.0
606 real(kind=kind_phys) , parameter :: zgam = 0.25
607!
608! Verical spectra
609!
610 real(kind=kind_phys) , parameter :: pind_wd = 5./3.
611 real(kind=kind_phys) , parameter :: sind_kz = 1.
612 real(kind=kind_phys) , parameter :: tind_kz = 3.
613 real(kind=kind_phys) , parameter :: stind_kz = sind_kz + tind_kz
614!
615! copies from kmob_ugwp namelist
616!
617 real(kind=kind_phys) :: nslope ! the GW sprctral slope at small-m
618 real(kind=kind_phys) :: lzstar
619 real(kind=kind_phys) :: lzmin
620 real(kind=kind_phys) :: lzmax
621 real(kind=kind_phys) :: lhmet
622 real(kind=kind_phys) :: tamp_mpa !amplitude for GEOS-5/MERRA-2
623 real(kind=kind_phys) :: tau_min ! min of GW MF 0.25 mPa
624 integer :: ilaunch
625 real(kind=kind_phys) :: gw_eff
626
627 real(kind=kind_phys) :: v_kxw, rv_kxw, v_kxw2
628
629
630
631!===========================================================================
632 integer :: nwav, nazd, nst
633 real(kind=kind_phys) :: eff
634
635 real(kind=kind_phys) :: zaz_fct, zms
636 real(kind=kind_phys), allocatable :: zci(:), zci4(:), zci3(:),zci2(:), zdci(:)
637 real(kind=kind_phys), allocatable :: zcosang(:), zsinang(:)
638 real(kind=kind_phys), allocatable :: lzmet(:), czmet(:), mkzmet(:), dczmet(:), dmkz(:)
639
640!
641! GW-eddy constants for wave-mode dissipation by background and stability of
642! "final" flow after application of GW-effects
643!
644 real(kind=kind_phys), parameter :: ipr_pt = 0.5
645 real(kind=kind_phys), parameter :: lturb = 30., sc2 = lturb*lturb ! stable on 80-km TL lmix ~ 500 met.
646 real(kind=kind_phys), parameter :: ulturb=150., sc2u = ulturb* ulturb ! unstable
647 real(kind=kind_phys), parameter :: ric =0.25
648 real(kind=kind_phys), parameter :: rimin = -10., prmin = 0.25
649 real(kind=kind_phys), parameter :: prmax = 4.0
650!
651 contains
652!============================================================================
653 subroutine initsolv_wmsdis(me, master, nwaves, nazdir, nstoch, effac, do_physb, kxw, version)
654
655! call initsolv_wmsdis(me, master, knob_ugwp_wvspec(2), knob_ugwp_azdir(2), &
656! knob_ugwp_stoch(2), knob_ugwp_effac(2), do_physb_gwsrcs, kxw,version)
657!
658 implicit none
659!
660!input-control for solvers:
661! nwaves, nazdir, nstoch, effac, do_physb, kxw
662!
663!
664 integer, intent(in) :: me, master, nwaves, nazdir, nstoch
665 integer, intent(in) :: version
666
667 real(kind=kind_phys), intent(in) :: effac, kxw
668 logical, intent(in) :: do_physb
669
670!
671!locals
672!
673 real(kind=kind_phys) :: dlzmet
674 real(kind=kind_phys) :: cstar,rcstar, nslope3, fnorm, zcin
675
676 integer :: inc, jk, jl, iazi
677!
678 real(kind=kind_phys) :: zang, zang1, znorm
679 real(kind=kind_phys) :: zx1, zx2, ztx, zdx, zxran, zxmin, zxmax, zx, zpexp
680 real(kind=kind_phys) :: fpc, fpc_dc
681 real(kind=kind_phys) :: ae1,ae2
682 if( nwaves == 0) then
683!
684! redefine from the deafault
685!
686 nwav = incdim
687 nazd = iazidim
688 nst = 0
689 eff = 1.0
690 gw_eff = eff
691 else
692!
693! from input.nml
694!
695 nwav = nwaves
696 nazd = nazdir
697 nst = nstoch
698 gw_eff = effac
699 endif
700
701
702 v_kxw = kxw ; v_kxw2 = v_kxw*v_kxw
703 rv_kxw = 1./v_kxw
704
705 allocate ( zci(nwav), zci4(nwav), zci3(nwav),zci2(nwav), zdci(nwav) )
706 allocate ( zcosang(nazd), zsinang(nazd) )
707 allocate (lzmet(nwav), czmet(nwav), mkzmet(nwav), dczmet(nwav), dmkz(nwav) )
708
709! if (me == master) then
710! print *, 'ugwp_v1/v0: init_gw_wmsdis_control '
711!
712! print *, 'ugwp_v1/v0: WMS_DIS launch layer ', ilaunch
713! print *, 'ugwp_v1/v0: WMS_DIS tot_mflux in mpa', tamp_mpa*1000.
714! print *, 'ugwp_v1/v0: WMS_DIS lhmet in km ' , lhmet*1.e-3
715! endif
716
717 zpexp = gptwo * 0.5 ! gptwo=2 , zpexp = 1.
718
719!
720! set up azimuth directions and some trig factors
721!
722!
723 zang = pi2 / float(nazd)
724
725! get normalization factor to ensure that the same amount of momentum
726! flux is directed (n,s,e,w) no mater how many azimuths are selected.
727!
728 znorm = 0.0
729 do iazi=1, nazd
730 zang1 = (iazi-1)*zang
731 zcosang(iazi) = cos(zang1)
732 zsinang(iazi) = sin(zang1)
733 znorm = znorm + abs(zcosang(iazi))
734 enddo
735! zaz_fct = 1.0
736 zaz_fct = 2.0 / znorm ! correction factor for azimuthal sums
737
738! define coordinate transform for "Ch" ....x = 1/c stretching transform
739! -----------------------------------------------
740!
741! x=1/Cphase transform
742! Scinocca 2003. x = 1/c stretching transform
743!
744 zxmax = 1.0 / zcimin
745 zxmin = 1.0 / zcimax
746 zxran = zxmax - zxmin
747 zdx = zxran / real(nwav-1) ! dkz
748!
749 ae1=zxran/zgam
750 zx1 = zxran/(exp(ae1)-1.0 ) ! zgam =1./4.
751 zx2 = zxmin - zx1
752
753!
754! computations for zci =1/zx, stretching "accuracy" is not "accurate" spectra transform
755! it represents additional "empirical" redistribution of "spectral" mode in C-space
756!
757 zms = pi2 / lzstar
758
759 do inc=1, nwav
760 ztx = real(inc-1)*zdx+zxmin
761 ae1 = (ztx-zxmin)/zgam
762 zx = zx1*exp(ae1)+zx2 !eq.(29-30),Scinocca-2003
763 zci(inc) = 1.0 /zx !
764 zdci(inc) = zci(inc)**2*(zx1/zgam)*exp(ae1)*zdx !
765 zci4(inc) = (zms*zci(inc))**4
766 zci2(inc) = (zms*zci(inc))**2
767 zci3(inc) = (zms*zci(inc))**3
768 enddo
769!
770!
771! alternatuve lzmax-lzmin without x=1/c transform
772!
773!
774 if (version == 1) then
775
776 dlzmet = (lzmax-lzmin)/ real(nwav-1)
777 do inc=1, nwav
778 lzmet(inc) = lzmin + (inc-1)*dlzmet
779 mkzmet(inc) = pi2/lzmet(inc)
780 zci(inc) =lzmet(inc)/(pi2/bnfix)
781 zci4(inc) = (zms*zci(inc))**4
782 zci2(inc) = (zms*zci(inc))**2
783 zci3(inc) = (zms*zci(inc))**3
784
785 enddo
786
787 zdx = (zci(nwav)-zci(1))/ real(nwav-1)
788 do inc=1, nwav
789 zdci(inc) = zdx
790 enddo
791
792 cstar = bnfix/zms
793 rcstar = 1./cstar
794 ENDIF ! if (version == 1) then
795
796 RETURN
797!=================== Diag prints after return ====================
798 if (me == master) then
799 print *
800 print *, 'ugwp_v0: zcimin=' , zcimin
801 print *, 'ugwp_v0: zcimax=' , zcimax
802 print *, 'ugwp_v0: zgam= ', zgam
803 print *
804
805! print *, ' ugwp_v1 nslope=', nslope
806 print *
807 print *, 'ugwp_v1: zcimin/zci=' , maxval(zci)
808 print *, 'ugwp_v1: zcimax/zci=' , minval(zci)
809 print *, 'ugwp_v1: cd_crit=', ucrit
810 print *, 'ugwp_v1: launch_level', ilaunch
811 print *, ' ugwp_v1 lzstar=', lzstar
812 print *, ' ugwp_v1 nslope=', nslope
813
814 print *
815 nslope3=nslope+3.0
816 do inc=1, nwav
817 zcin =zci(inc)*rcstar
818 fpc = rcstar*(zcin*zcin)/(1.+ zcin**nslope3)
819 fpc_dc = fpc * zdci(inc)
820 write(6,111) inc, zci(inc), zdci(inc),ucrit, fpc, fpc_dc, 6.28e-3/bnfix*zci(inc)
821 enddo
822 endif
823
824
825
826 111 format( 'wms-zci', i4, 7 (3x, f8.3))
827
828 end subroutine initsolv_wmsdis
829!
830!
831 end module ugwp_wmsdis_init