cires_ugwp_initialize.F90

 
!  initialization of ugwp_common_v0
!  init gw-solvers (1,2) .. no UFS-funds for (3,4) tests
!  init gw-source specifications
!  init gw-background dissipation
!===============================    
 
    module ugwp_common_v0
!
     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_v0
!
!
!===================================================
!
!Part-1 init =>   wave dissipation + RFriction
!
!===================================================
     subroutine init_global_gwdis_v0(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_v0
 
     
! ========================================================================
! Part 2 - sources
!      wave  sources
! ========================================================================
!
!    ugwpv0_oro_init
!
!=========================================================================
     module ugwpv0_oro_init
 
     use ugwp_common_v0, 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_v0(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_v0
!
 
    end module ugwpv0_oro_init
!=============================== end of GW  sources
!
!  init specific  gw-solvers (1,2,3,4)
!
 
!===============================
!  Part -3  init  wave solvers
!===============================
 
  module ugwpv0_lsatdis_init
     implicit none
 
      integer  :: nwav, nazd
      integer  :: nst
      real     :: eff
      integer, parameter  :: incdim = 4, iazdim = 4
!
     contains
 
     subroutine initsolv_lsatdis_v0(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_v0
!
  end module ugwpv0_lsatdis_init
!
!
  module ugwpv0_wmsdis_init
 
    use ugwp_common_v0, 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  :: 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_v0(me, master,  nwaves, nazdir, nstoch, effac, do_physb, 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_v0
!
  end module ugwpv0_wmsdis_init