!===============================
! cu-cires ugwp-scheme
! initialization of selected
! init gw-solvers (1,2,3,4)
! init gw-source specifications
! init gw-background dissipation
!==============================
!
! Part-0 specifications of common constants, limiters and "criiical" values
! module oro_state
! integer, parameter :: kind_phys=8
! integer, parameter :: nvaroro=14
! real (kind=kind_phys), allocatable :: oro_stat(:, :)
! contains
! subroutine fill_oro_stat(nx, oc, oa4, clx4, theta, gamm, sigma, elvmax, hprime)
! real (kind=kind_phys),dimension(nx) :: oc, theta, gamm, sigma, elvmax, hprime
! real(kind=kind_phys),dimension(nx,4) :: oa4, clx4
! integer :: i
! do i=1, nx
! oro_stat(i,1) = hprime(i)
! oro_stat(i,2) = oc(i)
! oro_stat(i,3:6) = oa4(i,1:4)
! oro_stat(i,7:10) = clx4(i,1:4)
! oro_stat(i,11) = theta(i)
! oro_stat(i,12) = gamm(i)
! oro_stat(i,13) = sigma(i)
! oro_stat(i,14) = elvmax(i)
! enddo
! end subroutine fill_oro_stat
! end module oro_state
module ugwp_common
!
use machine, only: kind_phys
use physcons, only : pi => con_pi, grav => con_g, rd => con_rd, &
rv => con_rv, cpd => con_cp, fv => con_fvirt,&
arad => con_rerth
implicit none
real(kind=kind_phys), parameter :: grcp = grav/cpd, rgrav = 1.0d0/grav, &
rdi = 1.0d0/rd, &
gor = grav/rd, gr2 = grav*gor, gocp = grav/cpd, &
rcpd = 1./cpd, rcpd2 = 0.5*rcpd, &
pi2 = pi + pi, omega1 = pi2/86400.0, &
omega2 = omega1+omega1, &
rad_to_deg=180.0/pi, deg_to_rad=pi/180.0, &
dw2min=1.0, bnv2min=1.e-6, velmin=sqrt(dw2min)
end module ugwp_common
!
!
!===================================================
!
!Part-1 init => wave dissipation + RFriction
!
!===================================================
subroutine init_global_gwdis(levs, zkm, pmb, kvg, ktg, krad, kion)
implicit none
integer :: levs
real, intent(in) :: zkm(levs), pmb(levs)
real, intent(out), dimension(levs+1) :: kvg, ktg, krad, kion
!
!locals + data
!
integer :: k
real, parameter :: vusurf = 2.e-5
real, parameter :: musurf = vusurf/1.95
real, parameter :: hpmol = 8.5
!
real, parameter :: kzmin = 0.1
real, parameter :: kturbo = 100.
real, parameter :: zturbo = 130.
real, parameter :: zturw = 30.
real, parameter :: inv_pra = 3. !kt/kv =inv_pr
!
real, parameter :: alpha = 1./86400./15.
!
real, parameter :: kdrag = 1./86400./10.
real, parameter :: zdrag = 100.
real, parameter :: zgrow = 50.
!
real :: vumol, mumol, keddy, ion_drag
!
do k=1, levs
vumol = vusurf*exp(-zkm(k)/hpmol)
mumol = musurf*exp(-zkm(k)/hpmol)
keddy = kturbo*exp(-((zkm(k)-zturbo) /zturw)**2)
kvg(k) = vumol + keddy
ktg(k) = mumol + keddy*inv_pra
krad(k) = alpha
!
ion_drag = kdrag
!
kion(k) = ion_drag
enddo
k= levs+1
kion(k) = kion(k-1)
krad(k) = krad(k-1)
kvg(k) = kvg(k-1)
ktg(k) = ktg(k-1)
!
end subroutine init_global_gwdis
!
!
subroutine rf_damp_init(levs, pa_rf, tau_rf, dtp, pmb, rfdis, rfdist, levs_rf)
implicit none
integer :: levs
real :: pa_rf, tau_rf
real :: dtp
real :: pmb(levs)
real :: rfdis(levs), rfdist(levs)
integer :: levs_rf
real :: krf, krfz
integer :: k
!
rfdis(1:levs) = 1.0
rfdist(1:levs) = 0.0
levs_rf = levs
if (tau_rf <= 0.0 .or. pa_rf == 0.0) return
krf = 1.0/(tau_rf*86400.0)
do k=levs, 1, -1
if(pmb(k) < pa_rf ) then ! applied only on constant pressure surfaces fixed pmb in "Pa"
krfz = krf*log(pa_rf/pmb(k))
rfdis(k) = 1.0/(1.+krfz*dtp)
rfdist(k) = (rfdis(k) -1.0)/dtp ! du/dtp
levs_rf = k
endif
enddo
end subroutine rf_damp_init
! ========================================================================
! Part 2 - sources
! wave sources
! ========================================================================
!
! ugwp_oro_init
!
!=========================================================================
module ugwp_oro_init
use ugwp_common, only : bnv2min, grav, grcp, fv, grav, cpd, grcp, pi
implicit none
!
! constants and "crirtical" values to run oro-mtb_gw physics
!
! choice of oro-scheme: strver = 'vay_2018' , 'gfs_2018', 'kdn_2005', 'smc_2000'
!
character(len=8) :: strver = 'gfs_2018'
character(len=8) :: strbase = 'gfs_2018'
real, parameter :: rimin=-10., ric=0.25
!
real, parameter :: efmin=0.5, efmax=10.0
real, parameter :: hpmax=2400.0, hpmin=25.0
real, parameter :: sigma_std=1./100., gamm_std=1.0
real, parameter :: frmax=10., frc =1.0, frmin =0.01
!
real, parameter :: ce=0.8, ceofrc=ce/frc, cg=0.5
real, parameter :: gmax=1.0, veleps=1.0, factop=0.5
!
real, parameter :: rlolev=50000.0
!
real, parameter :: hncrit=9000. ! max value in meters for elvmax
! hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt
real, parameter :: sigfac=4.0 ! mb3a expt test for elvmax factor
real, parameter :: hminmt=50. ! min mtn height (*j*)
real, parameter :: minwnd=1.0 ! min wind component (*j*)
real, parameter :: dpmin=5000.0 ! minimum thickness of the reference layer in pa
real, parameter :: kxoro=6.28e-3/200. !
real, parameter :: coro = 0.0
integer, parameter :: nridge=2
real :: cdmb ! scale factors for mtb
real :: cleff ! scale factors for orogw
integer :: nworo ! number of waves
integer :: nazoro ! number of azimuths
integer :: nstoro ! flag for stochastic launch above SG-peak
integer, parameter :: mdir = 8
real, parameter :: fdir=.5*mdir/pi
integer nwdir(mdir)
data nwdir/6,7,5,8,2,3,1,4/
save nwdir
real, parameter :: odmin = 0.1, odmax = 10.0
!------------------------------------------------------------------------------
! small-scale orography parameters for TOFD of Beljaars et al., 2004, QJRMS
!------------------------------------------------------------------------------
integer, parameter :: n_tofd = 2 ! depth of SSO for TOFD compared with Zpbl
real, parameter :: const_tofd = 0.0759 ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
real, parameter :: ze_tofd = 1500.0 ! BJ's z-decay in meters
real, parameter :: a12_tofd = 0.0002662*0.005363 ! BJ's k-spect const for sigf2 * a1*a2*exp(-[z/zdec]**1.5]
real, parameter :: ztop_tofd = 10.*ze_tofd ! no TOFD > this height too higher 15 km
!------------------------------------------------------------------------------
!
real, parameter :: fcrit_sm = 0.7, fcrit_sm2 = fcrit_sm * fcrit_sm
real, parameter :: fcrit_gfs = 0.7
real, parameter :: fcrit_mtb = 0.7
real, parameter :: lzmax = 18.e3 ! 18 km
real, parameter :: mkzmin = 6.28/lzmax
real, parameter :: mkz2min = mkzmin*mkzmin
real, parameter :: zbr_pi = (3.0/2.0)*pi
real, parameter :: zbr_ifs = 0.5*pi
contains
!
subroutine init_oro_gws(nwaves, nazdir, nstoch, effac, &
lonr, kxw, cdmbgwd )
!
!
integer :: nwaves, nazdir, nstoch
integer :: lonr
real :: cdmbgwd(2) ! scaling factors for MTb (1) & (2) for cleff = cleff * cdmbgwd(2)
! high res-n "larger" MTB and "less-active" cleff in GFS-2018
real :: cdmbX
real :: kxw
real :: effac ! it is analog of cdmbgwd(2) for GWs, off for now
!-----------------------------! GFS-setup for cdmb & cleff
! cdmb = 4.0 * (192.0/IMX)
! cleff = 0.5E-5 / SQRT(IMX/192.0) = 0.5E-5*SQRT(192./IMX)
!
real, parameter :: lonr_refmb = 4.0 * 192.0
real, parameter :: lonr_refgw = 192.0
! copy to "ugwp_oro_init" => nwaves, nazdir, nstoch
nworo = nwaves
nazoro = nazdir
nstoro = nstoch
cdmbX = lonr_refmb/float(lonr)
cdmb = cdmbX
if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
cleff = 0.5e-5 * sqrt(lonr_refgw/float(lonr)) !* effac
!!! cleff = kxw * sqrt(lonr_refgw/float(lonr)) !* effac
if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
!
!....................................................................
! higher res => smaller h' ..&.. higher kx
! flux_gwd ~ 'u'^2*kx/kz ~kxu/n ~1/dx *u/n tau ~ h'*h'*kx*kx = const (h'-less kx-grow)
!....................................................................
!
! print *, ' init_oro_gws 2-1cdmb', cdmbgwd(2), cdmbgwd(1)
end subroutine init_oro_gws
!
end module ugwp_oro_init
! =========================================================================
!
! ugwp_conv_init
!
!=========================================================================
module ugwp_conv_init
implicit none
real :: eff_con ! scale factors for conv GWs
integer :: nwcon ! number of waves
integer :: nazcon ! number of azimuths
integer :: nstcon ! flag for stochastic choice of launch level above Conv-cloud
real :: con_dlength
real :: con_cldf
real, parameter :: cmin = 5 !2.5
real, parameter :: cmax = 95. !82.5
real, parameter :: cmid = 22.5
real, parameter :: cwid = cmid
real, parameter :: bns = 2.e-2, bns2 = bns*bns, bns4=bns2*bns2
real, parameter :: mstar = 6.28e-3/2. ! 2km
real :: dc
real, allocatable :: ch_conv(:), spf_conv(:)
real, allocatable :: xaz_conv(:), yaz_conv(:)
contains
!
subroutine init_conv_gws(nwaves, nazdir, nstoch, effac, &
lonr, kxw, cgwf)
use ugwp_common, only : pi2, arad
implicit none
integer :: nwaves, nazdir, nstoch
integer :: lonr
real :: cgwf(2)
real :: kxw, effac
real :: work1 = 0.5
real :: chk, tn4, snorm
integer :: k
nwcon = nwaves
nazcon = nazdir
nstcon = nstoch
eff_con = effac
con_dlength = pi2*arad/float(lonr)
con_cldf = cgwf(1) * work1 + cgwf(2) *(1.-work1)
!
! allocate & define spectra in "selected direction": "dc" "ch(nwaves)"
!
if (.not. allocated(ch_conv)) allocate (ch_conv(nwaves))
if (.not. allocated(spf_conv)) allocate (spf_conv(nwaves))
if (.not. allocated(xaz_conv)) allocate (xaz_conv(nazdir))
if (.not. allocated(yaz_conv)) allocate (yaz_conv(nazdir))
dc = (cmax-cmin)/float(nwaves-1)
!
! we may use different spectral "shapes"
! for example FVS-93 "Desabeius"
! E(s=1, t=3,m, w, k) ~ m^s/(m*^4 + m^4) ~ m^-3 saturated tail
!
do k = 1,nwaves
chk = cmin + (k-1)*dc
tn4 = (mstar*chk)**4
ch_conv(k) = chk
spf_conv(k) = bns4*chk/(bns4+tn4)
enddo
snorm = sum(spf_conv)
spf_conv = spf_conv/snorm*1.5
call init_nazdir(nazdir, xaz_conv, yaz_conv)
end subroutine init_conv_gws
end module ugwp_conv_init
!=========================================================================
!
! ugwp_fjet_init
!
!=========================================================================
module ugwp_fjet_init
implicit none
real :: eff_fj ! scale factors for conv GWs
integer :: nwfj ! number of waves
integer :: nazfj ! number of azimuths
integer :: nstfj ! flag for stochastic choice of launch level above Conv-cloud
!
real, parameter :: fjet_trig=0. ! if ( abs(frgf) > fjet_trig ) launch GW-packet
real, parameter :: cmin = 2.5
real, parameter :: cmax = 67.5
real :: dc
real, allocatable :: ch_fjet(:) , spf_fjet(:)
real, allocatable :: xaz_fjet(:), yaz_fjet(:)
contains
subroutine init_fjet_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
use ugwp_common, only : pi2, arad
implicit none
integer :: nwaves, nazdir, nstoch
integer :: lonr
real :: kxw, effac , chk
integer :: k
nwfj = nwaves
nazfj = nazdir
nstfj = nstoch
eff_fj = effac
if (.not. allocated(ch_fjet)) allocate (ch_fjet(nwaves))
if (.not. allocated(spf_fjet)) allocate (spf_fjet(nwaves))
if (.not. allocated(xaz_fjet)) allocate (xaz_fjet(nazdir))
if (.not. allocated(yaz_fjet)) allocate (yaz_fjet(nazdir))
dc = (cmax-cmin)/float(nwaves-1)
do k = 1,nwaves
chk = cmin + (k-1)*dc
ch_fjet(k) = chk
spf_fjet(k) = 1.0
enddo
call init_nazdir(nazdir, xaz_fjet, yaz_fjet)
end subroutine init_fjet_gws
end module ugwp_fjet_init
!
!=========================================================================
!
!
module ugwp_okw_init
!=========================================================================
implicit none
real :: eff_okw ! scale factors for conv GWs
integer :: nwokw ! number of waves
integer :: nazokw ! number of azimuths
integer :: nstokw ! flag for stochastic choice of launch level above Conv-cloud
!
real, parameter :: okw_trig=0. ! if ( abs(okwp) > okw_trig ) launch GW-packet
real, parameter :: cmin = 2.5
real, parameter :: cmax = 67.5
real :: dc
real, allocatable :: ch_okwp(:), spf_okwp(:)
real, allocatable :: xaz_okwp(:), yaz_okwp(:)
contains
!
subroutine init_okw_gws(nwaves, nazdir, nstoch, effac, lonr, kxw)
use ugwp_common, only : pi2, arad
implicit none
integer :: nwaves, nazdir, nstoch
integer :: lonr
real :: kxw, effac , chk
integer :: k
nwokw = nwaves
nazokw = nazdir
nstokw = nstoch
eff_okw = effac
if (.not. allocated(ch_okwp)) allocate (ch_okwp(nwaves))
if (.not. allocated(spf_okwp)) allocate (spf_okwp(nwaves))
if (.not. allocated(xaz_okwp)) allocate (xaz_okwp(nazdir))
if (.not. allocated(yaz_okwp)) allocate (yaz_okwp(nazdir))
dc = (cmax-cmin)/float(nwaves-1)
do k = 1,nwaves
chk = cmin + (k-1)*dc
ch_okwp(k) = chk
spf_okwp(k) = 1.
enddo
call init_nazdir(nazdir, xaz_okwp, yaz_okwp)
end subroutine init_okw_gws
end module ugwp_okw_init
!=============================== end of GW sources
!
! init specific gw-solvers (1,2,3,4)
!
!===============================
! Part -3 init wave solvers
!===============================
module ugwp_lsatdis_init
implicit none
integer :: nwav, nazd
integer :: nst
real :: eff
integer, parameter :: incdim = 4, iazdim = 4
!
contains
subroutine initsolv_lsatdis(me, master, nwaves, nazdir, nstoch, effac, do_physb, kxw)
implicit none
!
integer :: me, master
integer :: nwaves, nazdir
integer :: nstoch
real :: effac
logical :: do_physb
real :: kxw
!
!locals: define azimuths and Ch(nwaves) - domain when physics-based soureces
! are not actibve
!
integer :: inc, jk, jl, iazi, i, j, k
if( nwaves == 0 .or. nstoch == 1 ) then
! redefine from the default
nwav = incdim
nazd = iazdim
nst = 0
eff = 1.0
else
! from input_nml multi-wave spectra
nwav = nwaves
nazd = nazdir
nst = nstoch
eff = effac
endif
!
end subroutine initsolv_lsatdis
!
end module ugwp_lsatdis_init
!
!
module ugwp_wmsdis_init
use ugwp_common, only : pi, pi2
implicit none
real, parameter :: maxdudt = 250.e-5
real, parameter :: hpscale= 7000., rhp2 = 0.5/hpscale
real, parameter :: omega2 = 2.*6.28/86400
real, parameter :: gptwo=2.0
real, parameter :: dked_min =0.01
real, parameter :: gssec = (6.28/30.)**2 ! max-value for bn2
real, parameter :: bv2min = (6.28/60./120.)**2 ! min-value for bn2 7.6(-7) 2 hrs
real, parameter :: minvel = 0.5
!
! make parameter list that will be passed to SOLVER
!
real, parameter :: v_kxw = 6.28e-3/200.
real, parameter :: v_kxw2 = v_kxw*v_kxw
real, parameter :: tamp_mpa = 30.e-3
real, parameter :: zfluxglob= 3.75e-3
real , parameter :: nslope=1 ! the GW sprctral slope at small-m
! integer, parameter :: klaunch=55 ! 32 - ~ 1km ;55 - 5.5 km ; 52 4.7km ; 60-7km index for selecting launch level
! integer, parameter :: ilaunch=klaunch
integer , parameter :: iazidim=4 ! number of azimuths
integer , parameter :: incdim=25 ! number of discrete cx - spectral elements in launch spectrum
real , parameter :: ucrit2=0.5
real , parameter :: zcimin = ucrit2
real , parameter :: zcimax = 125.0
real , parameter :: zgam = 0.25
real , parameter :: zms_l = 2000.0, zms = pi2 / zms_l, zmsi = 1.0 / zms
integer :: ilaunch
real :: gw_eff
!===========================================================================
integer :: nwav, nazd, nst
real :: eff
real :: zaz_fct
real, allocatable :: zci(:), zci4(:), zci3(:),zci2(:), zdci(:)
real, allocatable :: zcosang(:), zsinang(:)
contains
!============================================================================
subroutine initsolv_wmsdis(me, master, nwaves, nazdir, nstoch, effac, do_physb, kxw)
! call initsolv_wmsdis(me, master, knob_ugwp_wvspec(2), knob_ugwp_azdir(2), &
! knob_ugwp_stoch(2), knob_ugwp_effac(2), do_physb_gwsrcs, kxw)
!
implicit none
!
!input -control for solvers:
! nwaves, nazdir, nstoch, effac, do_physb, kxw
!
!
integer :: me, master, nwaves, nazdir, nstoch
real :: effac, kxw
logical :: do_physb
!
!locals
!
integer :: inc, jk, jl, iazi
!
real :: zang, zang1, znorm
real :: zx1, zx2, ztx, zdx, zxran, zxmin, zxmax, zx, zpexp
if( nwaves == 0) then
!
! redefine from the deafault
!
nwav = incdim
nazd = iazidim
nst = 0
eff = 1.0
gw_eff = eff
else
!
! from input.nml
!
nwav = nwaves
nazd = nazdir
nst = nstoch
gw_eff = effac
endif
allocate ( zci(nwav), zci4(nwav), zci3(nwav),zci2(nwav), zdci(nwav) )
allocate ( zcosang(nazd), zsinang(nazd) )
if (me == master) then
print *, 'ugwp_v0: init_gw_wmsdis_control '
! print *, 'ugwp_v0: WMSDIS launch layer ', klaunch
print *, 'ugwp_v0: WMSDIS launch layer ', ilaunch
print *, 'ugwp_v0: WMSDID tot_mflux in mpa', tamp_mpa*1000.
endif
zpexp = gptwo * 0.5 ! gptwo=2 , zpexp = 1.
!
! set up azimuth directions and some trig factors
!
!
zang = pi2 / float(nazd)
! get normalization factor to ensure that the same amount of momentum
! flux is directed (n,s,e,w) no mater how many azimuths are selected.
!
znorm = 0.0
do iazi=1, nazd
zang1 = (iazi-1)*zang
zcosang(iazi) = cos(zang1)
zsinang(iazi) = sin(zang1)
znorm = znorm + abs(zcosang(iazi))
enddo
! zaz_fct = 1.0
zaz_fct = 2.0 / znorm ! correction factor for azimuthal sums
! define coordinate transform for "Ch" ....x = 1/c stretching transform
! -----------------------------------------------
! note that this is expresed in terms of the intrinsic phase speed
! at launch ci=c-u_o so that the transformation is identical
! see eq. 28-30 of scinocca 2003. x = 1/c stretching transform
!
zxmax = 1.0 / zcimin
zxmin = 1.0 / zcimax
zxran = zxmax - zxmin
zdx = zxran / real(nwav-1) ! dkz
!
zx1 = zxran/(exp(zxran/zgam)-1.0 ) ! zgam =1./4.
zx2 = zxmin - zx1
!
! computations for zci =1/zx
! if(lgacalc) zgam=(zxmax-zxmin)/log(zxmax/zxmin)
! zx1=zxran/(exp(zxran/zgam)-1.0_jprb)
! zx2=zxmin-zx1
! zms = pi2 / zms_l
do inc=1, nwav
ztx = real(inc-1)*zdx+zxmin
zx = zx1*exp((ztx-zxmin)/zgam)+zx2 !eq. 29 of scinocca 2003
zci(inc) = 1.0 /zx !eq. 28 of scinocca 2003
zdci(inc) = zci(inc)**2*(zx1/zgam)*exp((ztx-zxmin)/zgam)*zdx !eq. 30 of scinocca 2003
zci4(inc) = (zms*zci(inc))**4
zci2(inc) = (zms*zci(inc))**2
zci3(inc) = (zms*zci(inc))**3
enddo
!
!
! all done and print-out
!
!
if (me == master) then
print *
print *, 'ugwp_v0: zcimin=' , zcimin
print *, 'ugwp_v0: zcimax=' , zcimax
print *, 'ugwp_v0: cd_crit=', zgam ! m/s precision for crit-level
print *, 'ugwp_v0: launch_level', ilaunch
print *, ' ugwp_v0 zms_l=', zms_l
print *, ' ugwp_vgw nslope=', nslope
print *
endif
end subroutine initsolv_wmsdis
!
! make a list of all-initilized parameters needed for "gw_solver_wmsdis"
!
end module ugwp_wmsdis_init
!=========================================================================
!
! work TODO for 2-extra WAM-solvers:
! DSPDIS (Hines)+ADODIS (Alexander-Dunkerton-Ortland)
!
!=========================================================================
subroutine init_dspdis
implicit none
end subroutine init_dspdis
subroutine init_adodis
implicit none
end subroutine init_adodis