!>\file cires_ugwpv1_module.F90 !! module cires_ugwpv1_module ! ! driver is called after pbl & before chem-parameterizations ! it uses ugwp_common (like phys_cons) and some module-param od solvers/sources init-modules !.................................................................................... ! order = dry-adj=>conv=mp-aero=>radiation -sfc/land- chem -> vertdiff-> [rf-gws]=> ion-re !................................................................................... ! ! use machine, only : kind_phys use ugwp_common, only : arad, pi, pi2, hpscale, rhp, rhp2, rh4, rhp4, khp, hpskm use ugwp_wmsdis_init, only : ilaunch, nslope, lhmet, lzmax, lzmin, lzstar use ugwp_wmsdis_init, only : tau_min, tamp_mpa implicit none logical :: module_is_initialized character(len=8) :: strsolver='pss-1986' logical :: do_physb_gwsrcs = .false. ! control for physics-based GW-sources logical :: do_rfdamp = .false. ! control for Rayleigh friction inside ugwp_driver integer, parameter :: idebug_gwrms=0 ! control for diag computaions pw wind-temp GW-rms and MF fluxs logical, parameter :: do_adjoro = .false. real(kind=kind_phys), parameter :: max_kdis = 450. ! 400 m2/s real(kind=kind_phys), parameter :: max_axyz = 450.e-5 ! 400 m/s/day real(kind=kind_phys), parameter :: max_eps = max_kdis*4.e-4 ! max_kdis*BN2 real(kind=kind_phys), parameter :: maxdudt = max_axyz real(kind=kind_phys), parameter :: maxdtdt = max_eps*1.e-3 ! max_kdis*BN2/cp real(kind=kind_phys), parameter :: dked_min = 0.01 real(kind=kind_phys), parameter :: dked_max = max_kdis ! ! ! Pr = Kv/Kt < 1 for upper layers; Pr_mol = 1./1.95 check it ! real(kind=kind_phys), parameter :: Pr_kvkt = 1./1. ! kv/kt = 1./3. real(kind=kind_phys), parameter :: Pr_kdis = Pr_kvkt/(1.+Pr_kvkt) real(kind=kind_phys), parameter :: iPr_ktgw =1./3., iPr_spgw=iPr_ktgw real(kind=kind_phys), parameter :: iPr_turb =1./3., iPr_mol =1.95 real(kind=kind_phys), parameter :: cd_ulim = 1.0 ! critical level precision or Lz ~ 0 ~dz of model real(kind=kind_phys), parameter :: linsat = 1.00 real(kind=kind_phys), parameter :: linsat2 = linsat*linsat real(kind=kind_phys), parameter :: ricrit = 0.25 real(kind=kind_phys), parameter :: frcrit = 0.50 integer :: knob_ugwp_version = 1 integer :: knob_ugwp_solver=1 ! 1, 2, 3, 4 - (linsat, ifs_2010, ad_gfdl, dsp_dis) integer, dimension(4) :: knob_ugwp_source=(/1,0,1,0/) ! [1,0,1,1] - (oro, fronts, conv, imbf-owp] integer, dimension(4) :: knob_ugwp_wvspec=(/1,32,32,32/) ! number of waves for- (oro, fronts, conv, imbf-owp] integer, dimension(4) :: knob_ugwp_azdir=(/2,4,4,4/) ! number of wave azimuths for-(oro, fronts, conv, imbf-owp] integer, dimension(4) :: knob_ugwp_stoch=(/0,0,0,0/) ! 0 - deterministic ; 1 - stochastic real(kind=kind_phys), dimension(4) :: knob_ugwp_effac=(/1.,1.,1.,1./) ! efficiency factors for- (oro, fronts, conv, imbf-owp] integer :: knob_ugwp_doaxyz=1 ! 1 -gwdrag integer :: knob_ugwp_doheat=1 ! 1 -gwheat integer :: knob_ugwp_dokdis=0 ! 1 -gwmixing integer :: knob_ugwp_ndx4lh = 2 ! n-number of "unresolved" "n*dx" for lh_gw integer :: knob_ugwp_nslope = 1 ! spectral"growth" S-slope of GW-energy spectra mkz^S real(kind=kind_phys) :: knob_ugwp_palaunch = 500.e2 ! fixed pressure layer in Pa for "launch" of NGWs real(kind=kind_phys) :: knob_ugwp_lzmax = 12.5e3 ! 12.5 km max-VERT-WL of GW-spectra real(kind=kind_phys) :: knob_ugwp_lzstar = 2.0e3 ! UTLS mstar = 6.28/lzstar 2-2.5 km real(kind=kind_phys) :: knob_ugwp_lzmin = 1.5e3 ! 1.5 km min-VERT-WL of GW-spectra real(kind=kind_phys) :: knob_ugwp_taumin = 0.25e-3 real(kind=kind_phys) :: knob_ugwp_tauamp = 7.75e-3 ! range from 30.e-3 to 3.e-3 ( space-borne values) real(kind=kind_phys) :: knob_ugwp_lhmet = 200.e3 ! 200 km logical :: knob_ugwp_tlimb = .true. character(len=8) :: knob_ugwp_orosolv='pss-1986' real(kind=kind_phys) :: kxw = 6.28/200.e3 ! single horizontal wavenumber of ugwp schemes ! integer :: ugwp_azdir integer :: ugwp_stoch integer :: ugwp_src integer :: ugwp_nws real(kind=kind_phys) :: ugwp_effac ! integer :: launch_level = 55 ! namelist /cires_ugwp_nml/ knob_ugwp_solver, knob_ugwp_source,knob_ugwp_wvspec, knob_ugwp_azdir, & knob_ugwp_stoch, knob_ugwp_effac,knob_ugwp_doaxyz, knob_ugwp_doheat, knob_ugwp_dokdis, & knob_ugwp_ndx4lh, knob_ugwp_version, knob_ugwp_palaunch, knob_ugwp_nslope, knob_ugwp_lzmax, & knob_ugwp_lzmin, knob_ugwp_lzstar, knob_ugwp_lhmet, knob_ugwp_tauamp, knob_ugwp_taumin, & knob_ugwp_tlimb, knob_ugwp_orosolv ! ! allocatable arrays, initilized during "cires_ugwp_init" & ! released during "cires_ugwp_finalize" ! real(kind=kind_phys), allocatable :: kvg(:), ktg(:), krad(:), kion(:) real(kind=kind_phys), allocatable :: zkm(:), pmb(:) real(kind=kind_phys), allocatable :: rfdis(:), rfdist(:) ! ! RF-not active now ! integer :: levs_rf real(kind=kind_phys) :: pa_rf, tau_rf ! ! simple modulation of tau_ngw by the total rain/precip strength ! real(kind=kind_phys), parameter :: rain_max=8.e-5, rain_lat=41.0, rain_lim=1.e-5 real(kind=kind_phys), parameter :: w_merra = 1.0, w_nomerra = 1.-w_merra, w_rain =1. real(kind=kind_phys), parameter :: mtau_rain = 1.e-3, ft_min =0.5, ft_max=2 real(kind=kind_phys), parameter :: tau_ngw_max = 20.e-3 ! 20 mPa real(kind=kind_phys), parameter :: tau_ngw_min = .20e-3 ! .2 mPa ! ! Bushell et al. (2015) tau = tau_rainum (~3.8 km) x sqrt(Precip/base_rainum) ! real(kind=kind_phys), parameter :: tau_rainum = 0.7488e-3 ! 0.74 mPa real(kind=kind_phys), parameter :: base_rainum = 0.1e-5 ! ~0.1 mm/day real(kind=kind_phys), parameter :: pbase_um =1./sqrt(base_rainum) * tau_rainum ! integer, parameter :: metoum_rain = 0 !================================================================= ! switches that can ba activated for NGW physics include/omit ! ! rotational, non-hydrostatic and eddy-dissipative ! F_coriol F_nonhyd F_kds !=================================================== real(kind=kind_phys), parameter :: F_coriol=1.0 ! Coriolis effects real(kind=kind_phys), parameter :: F_nonhyd=1.0 ! Nonhydrostatic waves real(kind=kind_phys), parameter :: F_kds =0.0 ! Eddy mixing due to GW-unstable below contains ! !----------------------------------------------------------------------- ! ! init of cires_ugwp (_init) called from CCPP cap file ! ! --------------------------------------------------------------------------------- ! non-ccpp .... ! ! subroutine cires_ugwp_init_v1 (me, master, nlunit, logunit, jdat_gfs, fn_nml2, & ! lonr, latr, levs, ak, bk, pref, dtp) !----------------------------------------------------------------------------------- subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & con_rerth, fn_nml2, input_nml_file, lonr, latr, levs, ak, bk, & pref, dtp, errmsg, errflg) ! ! input_nml_file ='input.nml'=fn_nml ..... OLD_namelist and cdmvgwd(4) Corrected Bug Oct 4 ! use netcdf use ugwp_oro_init, only : init_oro_gws use ugwp_conv_init, only : init_conv_gws use ugwp_fjet_init, only : init_fjet_gws use ugwp_okw_init, only : init_okw_gws use ugwp_lsatdis_init, only : initsolv_lsatdis use ugwp_wmsdis_init, only : initsolv_wmsdis use ugwp_wmsdis_init, only : ilaunch, nslope, lhmet, lzmax, lzmin, lzstar use ugwp_wmsdis_init, only : tau_min, tamp_mpa implicit none integer, intent (in) :: me integer, intent (in) :: master integer, intent (in) :: nlunit integer, intent (in) :: logunit integer, intent (in) :: lonr integer, intent (in) :: levs integer, intent (in) :: latr integer, intent (in) :: jdat_gfs(8) real(kind=kind_phys), intent (in) :: ak(levs+1), bk(levs+1), pref real(kind=kind_phys), intent (in) :: dtp ! ! consider to retire them ! real(kind=kind_phys), intent (in) :: con_pi, con_rerth character(len=64), intent (in) :: fn_nml2 character(len=*), intent (in) :: input_nml_file(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! character, intent (in) :: input_nml_file ! integer :: ios logical :: exists integer :: ncid, iernc, vid, dimid, status integer :: k integer :: ddd_ugwp, curday_ugwp ! integer :: version ! Initialize CCPP error handling variables errmsg = '' errflg = 0 #ifdef INTERNAL_FILE_NML read (input_nml_file, nml = cires_ugwp_nml) #else if (me == master) print *, trim (fn_nml2), ' GW-namelist file ' inquire (file =trim (fn_nml2) , exist = exists) ! if (.not. exists) then write(errmsg,'(3a)') 'cires_ugwpv1_init: namelist file: ', trim (fn_nml2), ' does not exist' errflg = 1 return else open (unit = nlunit, file = trim(fn_nml2), action = 'read', status = 'old', iostat = ios) endif rewind (nlunit) read (nlunit, nml = cires_ugwp_nml) close (nlunit) #endif strsolver= knob_ugwp_orosolv curday_ugwp = jdat_gfs(1)*10000 + jdat_gfs(2)*100 +jdat_gfs(3) call calendar_ugwp(jdat_gfs(1), jdat_gfs(2), jdat_gfs(3), ddd_ugwp) ! write version number and namelist to log file if (me == master) then write (logunit, *) " ================================================================== " write (logunit, *) "CCPP cires_ugwp_namelist_extended_v1" write (logunit, nml = cires_ugwp_nml) write (logunit, *) " ================================================================== " write (6, *) " ================================================================== " write (6, *) "CCPP cires_ugwp_namelist_extended_v1" write (6, nml = cires_ugwp_nml) write (6, *) " ================================================================== " write (6, *) "calendar_ugwp ddd_ugwp=", ddd_ugwp write (6, *) "calendar_ugwp curday_ugwp=", curday_ugwp write (6, *) " ================================================================== " write (6, *) ddd_ugwp, ' jdat_gfs ddd of year ' endif ! ! effective kxw - resolution-aware ! ! kxw = pi2/knob_ugwp_lhmet ! ! ! init global background dissipation for ugwp -> 4d-variable for fv3wam linked with pbl-vert_diff ! ! ! allocate(fcor(latr), fcor2(latr) ) ! allocate( kvg(levs+1), ktg(levs+1) ) allocate( krad(levs+1), kion(levs+1) ) allocate( zkm(levs), pmb(levs) ) ! ! ak -pa bk-dimensionless from surf => tol_lid_pressure =0 ! do k=1, levs pmb(k) = ak(k) + pref*bk(k) ! Pa -unit Pref = 1.e5, pmb = Pa zkm(k) = -hpskm*alog(pmb(k)/pref) enddo ! ! find ilaunch ! if (knob_ugwp_palaunch > 900.e2) then write(errmsg,'(a,e16.7)') 'cires_ugwpv1_init: unrealistic value of knob_ugwp_palaunch', knob_ugwp_palaunch errflg = 1 return endif do k=levs, 1, -1 if (pmb(k) .gt. knob_ugwp_palaunch ) exit enddo launch_level = max(k-1, 5) ! above 5-layers from the surface if (me == master) then print *, 'cires_ugwpv1 klev_ngw =', launch_level, nint(pmb(launch_level)) endif ! ! Part-1 :init_global_gwdis again "damn"-con_p ! call init_global_gwdis(levs, zkm, pmb, kvg, ktg, krad, kion, me, master) ! ! Part-2 :init_SOURCES_gws ! ! ! call init-solver for "stationary" multi-wave spectra and sub-grid oro ! call init_oro_gws( knob_ugwp_wvspec(1), knob_ugwp_azdir(1), & knob_ugwp_stoch(1), knob_ugwp_effac(1), lonr, kxw ) ! ! call init-sources for "non-sationary" multi-wave spectra ! do_physb_gwsrcs=.true. IF (do_physb_gwsrcs) THEN ! if (me == master) print *, ' do_physb_gwsrcs ', do_physb_gwsrcs, ' in cires_ugwp_init_modv1 ' if (knob_ugwp_wvspec(4) > 0) then ! okw call init_okw_gws(knob_ugwp_wvspec(4), knob_ugwp_azdir(4), & knob_ugwp_stoch(4), knob_ugwp_effac(4), & lonr, kxw ) ! if (me == master) print *, ' init_okw_gws ' endif if (knob_ugwp_wvspec(3) > 0) then ! fronts call init_fjet_gws(knob_ugwp_wvspec(3), knob_ugwp_azdir(3), & knob_ugwp_stoch(3), knob_ugwp_effac(3), & lonr, kxw ) ! if (me == master) print *, ' init_fjet_gws ' endif if (knob_ugwp_wvspec(2) > 0) then ! conv : con_pi, con_rerth, call init_conv_gws(knob_ugwp_wvspec(2), knob_ugwp_azdir(2), & knob_ugwp_stoch(2), knob_ugwp_effac(2), & lonr, kxw ) ! if (me == master) & ! print *, ' init_convective GWs ', knob_ugwp_wvspec(2), knob_ugwp_azdir(2) endif ENDIF !IF (do_physb_gwsrcs) !====================== ! Part-3 :init_SOLVERS ! ===================== ! ! call init-solvers for "broad" non-stationary multi-wave spectra ! if (knob_ugwp_solver==1) then ! kxw = pi2/lhmet call initsolv_lsatdis(me, master, knob_ugwp_wvspec(2), knob_ugwp_azdir(2), & knob_ugwp_stoch(2), knob_ugwp_effac(2), do_physb_gwsrcs, kxw ) endif if (knob_ugwp_solver == 2) then ! ! re-assign from namelists ! nslope = knob_ugwp_nslope ! the GW sprctral slope at small-m lzstar = knob_ugwp_lzstar lzmax = knob_ugwp_lzmax lzmin = knob_ugwp_lzmin lhmet = knob_ugwp_lhmet tamp_mpa =knob_ugwp_tauamp !amplitude for GEOS-5/MERRA-2 tau_min =knob_ugwp_taumin ! min of GW MF 0.25 mPa ilaunch = launch_level kxw = pi2/lhmet 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, knob_ugwp_version) endif !====================== module_is_initialized = .true. if (me == master) print *, ' CIRES_ugwpV1 is initialized ', module_is_initialized end subroutine cires_ugwpv1_init !============================================= subroutine cires_ugwp_advance !----------------------------------------------------------------------- ! FV3-dycore and CCPP-physics has limited options to ! add "horizontal" gradients of winds and temp-re to ! compute GW-triggers: reserved option if it will be funded ...... ! ! the day-to-day variable sources/spectra and diagnostics for stochastic "triggers" ! ! diagnose GW-source functions * FGF + OKWP + SGO/CONV from IAU-fields ! and use for stochastic GWP-sources "memory" ! ! this option is not active due to "weak" flexibility ! in communication between "ccpp/gfsphysics" and FV3-dycore ! extension of State%in is needed to pass horizontal gradients ! winds and temperature to compute "spontatneous" GW triggers !----------------------------------------------------------------------- implicit none ! ! update GW sources and dissipation ! a) physics-based GW triggers eliminated from cires_ugwpv1_triggers.F90 ! b) stochastic-based spectra and amplitudes is not considered ! c) use "memory" on GW-spectra from previous time-step is not considered ! d) update "background" dissipation of GWs as needed (option for FV3WAM) ! end subroutine cires_ugwp_advance ! ! ----------------------------------------------------------------------- ! finalize of cires_ugwp_dealloc ! ----------------------------------------------------------------------- subroutine cires_ugwp_dealloc ! ! deallocate sources/spectra & some diagnostics need to find where "deaalocate them" ! before "end" of the FV3GFS ! implicit none ! ! deallocate arrays employed in: ! cires_ugwp_advance / cires_ugwp_driver / cires_ugwp_init ! if (allocated (kvg)) deallocate (kvg) if (allocated (ktg)) deallocate (ktg) if (allocated (krad)) deallocate (krad) if (allocated (kion)) deallocate (kion) if (allocated (zkm)) deallocate (zkm) if (allocated (pmb)) deallocate (pmb) ! if (allocated (ugwp_taulat)) deallocate(ugwp_taulat) ! if (allocated (tau_limb)) deallocate (tau_limb) ! if (allocated (days_limb)) deallocate(days_limb) end subroutine cires_ugwp_dealloc ! ! subroutine calendar_ugwp(yr, mm, dd, ddd_ugwp) ! ! computes day of year to get tau_limb forcing written with 1-day precision ! implicit none integer, intent(in) :: yr, mm, dd integer :: ddd_ugwp integer :: iw3jdn integer :: jd1, jddd jd1 = iw3jdn(yr,1,1) jddd = iw3jdn(yr,mm,dd) ddd_ugwp = jddd-jd1+1 end subroutine calendar_ugwp subroutine ngwflux_update(me, master, im, levs, kdt, ddd, curdate, & tau_ddd, xlatd, sinlat,coslat, rain, tau_ngw) use machine, only: kind_phys implicit none !input integer, intent(in) :: me, master !, jdat(8) integer, intent(in) :: im, levs, kdt integer, intent(in) :: ddd, curdate ! integer, intent(in), dimension(im) :: j1_tau, j2_tau ! real(kind=kind_phys), intent(in), dimension(im) :: ddy_j2tau, ddy_j1tau real(kind=kind_phys), intent(in), dimension(im) :: xlatd, sinlat,coslat real(kind=kind_phys), intent(in), dimension(im) :: rain, tau_ddd real(kind=kind_phys), intent(inout), dimension(im) :: tau_ngw ! ! locals ! integer :: i, j1, j2, k, it1, it2, iday real(kind=kind_phys) :: tem, tx1, tx2, w1, w2, wlat, rw1, rw2 real(kind=kind_phys) :: tau_rain, flat_rain, tau_3dt ! ! code below inside cires_tauamf_data.F90 ! it1 = 2 ! do iday=1, ntau_d2t ! if (float(ddd) .lt. days_limb(iday) ) then ! it2 = iday ! exit ! endif ! enddo ! it2 = min(it2,ntau_d2t) ! it1 = max(it2-1,1) ! if (it2 > ntau_d2t ) then ! print *, ' it1, it2, ntau_d2t ', it1, it2, ntau_d2t ! stop ! endif ! w2 = (float(ddd)-days_limb(it1))/(days_limb(it2)-days_limb(it1)) ! w1 = 1.0-w2 ! do i=1, im ! j1 = j1_tau(i) ! j2 = j2_tau(i) ! tx1 = tau_limb(j1, it1)*ddy_j1tau(i)+tau_limb(j2, it1)*ddy_j2tau(i) ! tx2 = tau_limb(j1, it2)*ddy_j1tau(i)+tau_limb(j2, it2)*ddy_j2tau(i) ! tau_ddd(i) = tx1*w1 + w2*tx2 ! ! add modulattion by the total "rain"-strength Yudin et al.(2020-FV3GFS) and Bushell et al. (2015-UM/METO) ! do i=1, im tau_3dt = tau_ngw(i) * w_merra + w_nomerra *tau_ddd(i) if (w_rain > 0. .and. rain(i) > 0.) then wlat = abs(xlatd(i)) if (wlat <= rain_lat .and. rain(i) > rain_lim) then flat_rain = wlat/rain_lat rw1 = 0.75 * flat_rain ; rw2 = 1.-rw1 tau_rain = tau_3dt * rw1 + rw2 * mtau_rain*min(rain_max, rain(i))/rain_lim tau_rain = tau_3dt*(1.-w_rain) + w_rain* tau_rain ! ! restict variations from the "tau_ngw" without precip-impact ! ! real, parameter :: ft_min =0.5*tau_g5 < tau_rain < ft_max =2. *tau_g5 ! if (tau_rain < ft_min *tau_3dt) tau_rain = ft_min *tau_3dt if (tau_rain > ft_max *tau_3dt) tau_rain = ft_max *tau_3dt tau_3dt = tau_rain endif if (metoum_rain == 1) then tau_rain = min( sqrt(rain(i))*pbase_um, tau_ngw_max) tau_3dt = max(tau_ngw_min, tau_rain) endif endif tau_ngw(i) = tau_3dt enddo end subroutine ngwflux_update ! end module cires_ugwpv1_module