MODULE nonoro_gwd_mix_mod IMPLICIT NONE REAL,allocatable,save :: du_eddymix_gwd(:,:) ! Zonal wind tendency due to GWD REAL,allocatable,save :: dv_eddymix_gwd(:,:) ! Meridional wind tendency due to GWD REAL,allocatable,save :: dh_eddymix_gwd(:,:) ! PT tendency due to GWD REAL,allocatable,save :: dq_eddymix_gwd(:,:,:) ! tracers tendency due to GWD REAL,allocatable,save :: de_eddymix_rto(:,:) ! Meridional wind tendency due to GWD REAL,allocatable,save :: df_eddymix_flx(:,:) ! Meridional wind tendency due to GWD !REAL,ALLOCATABLE,SAVE :: east_gwstress(:,:) ! Profile of eastward stress !REAL,ALLOCATABLE,SAVE :: west_gwstress(:,:) ! Profile of westward stress LOGICAL, save :: calljliu_gwimix ! flag for using non-orographic GW-induced mixing !$OMP THREADPRIVATE(du_eddymix_gwd,dv_eddymix_gwd,dh_eddymix_gwd,dq_eddymix_gwd,de_eddymix_rto,df_eddymix_flx,calljliu_gwimix) !,east_gwstress,west_gwstress) CONTAINS SUBROUTINE NONORO_GWD_MIX(ngrid,nlayer,DTIME,nq,cpnew, rnew, pp, & zmax_therm, pt, pu, pv, pq, pht, pdt, pdu, pdv, pdq, pdht, & d_pq, d_t, d_u, d_v) !-------------------------------------------------------------------------------- ! Parametrization of the eddy diffusion coefficient due to a discrete ! number of gravity waves. ! J.LIU ! Version 01, Gaussian distribution of the source ! LMDz model online version ! !--------------------------------------------------------------------------------- use comcstfi_h, only: g, pi, r,rcp USE ioipsl_getin_p_mod, ONLY : getin_p use vertical_layers_mod, only : presnivs use geometry_mod, only: cell_area use write_output_mod, only: write_output #ifdef CPP_XIOS use xios_output_mod, only: send_xios_field #endif implicit none include "callkeys.h" CHARACTER (LEN=20) :: modname='NONORO_GWD_MIX' CHARACTER (LEN=80) :: abort_message ! 0. DECLARATIONS: ! 0.1 INPUTS INTEGER, intent(in):: ngrid ! number of atmospheric columns INTEGER, intent(in):: nlayer ! number of atmospheric layers INTEGER, INTENT(IN) :: nq ! number of tracer species in traceurs.def ! integer, parameter:: nesp =42 ! number of traceurs in chemistry REAL, intent(in):: DTIME ! Time step of the Physics(s) REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m) ! REAL, intent(in):: loss(nesp) ! Chemical reaction loss rate REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere REAL, intent(in):: pp(ngrid,nlayer) ! Pressure at full levels(Pa) REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels(K) REAL, intent(in):: pu(ngrid,nlayer) ! Zonal wind at full levels(m/s) REAL, intent(in):: pv(ngrid,nlayer) ! Meridional winds at full levels(m/s) REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq REAL, INTENT(IN) :: pht(ngrid,nlayer) ! advected field of potential temperature REAL, INTENT(IN) :: pdht(ngrid,nlayer) ! tendancy of potential temperature REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq)! tendancy field pq REAL,INTENT(in) :: pdt(ngrid,nlayer) ! Tendency on temperature (K/s) REAL,INTENT(in) :: pdu(ngrid,nlayer) ! Tendency on zonal wind (m/s/s) REAL,INTENT(in) :: pdv(ngrid,nlayer) ! Tendency on meridional wind (m/s/s) ! 0.2 OUTPUTS ! REAL, intent(out):: zustr(ngrid) ! Zonal surface stress ! REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature (K/s) due to gravity waves (not used set to zero) REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind (m/s/s) due to gravity waves REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves REAL, INTENT(out) :: d_pq(ngrid,nlayer,nq)! tendancy field pq REAL :: d_h(ngrid, nlayer) ! Tendency on PT (T/s/s) due to gravity waves mixing ! 0.3 INTERNAL ARRAYS REAL :: TT(ngrid, nlayer) ! Temperature at full levels REAL :: RHO(ngrid, nlayer) ! Mass density at full levels REAL :: UU(ngrid, nlayer) ! Zonal wind at full levels REAL :: VV(ngrid, nlayer) ! Meridional winds at full levels REAL :: HH(ngrid, nlayer) ! potential temperature at full levels REAL :: BVLOW(ngrid) ! initial N at each grid (not used) INTEGER II, JJ, LL, QQ ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED REAL, parameter:: DELTAT = 24. * 3600. ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves INTEGER JK, JP, JO, JW ! Loop index for NK,NP,NO, and NW REAL, save :: kmax ! Max horizontal wavenumber (lambda min,lambda=2*pi/kmax),kmax=N/u, u(=30~50) zonal wind velocity at launch altitude !$OMP THREADPRIVATE(kmax) REAL, save :: kmin ! Min horizontal wavenumber (lambda max = 314 km,lambda=2*pi/kmin) !$OMP THREADPRIVATE(kmin) REAL kstar ! Min horizontal wavenumber constrained by horizontal resolution REAL :: max_k(ngrid) ! max_k=max(kstar,kmin) REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity (not used) REAL CPHA(ngrid) ! absolute PHASE VELOCITY frequency REAL ZK(NW, ngrid) ! Horizontal wavenumber amplitude REAL ZP(NW, ngrid) ! Horizontal wavenumber angle REAL ZO(NW, ngrid) ! Absolute frequency REAL maxp(NW,ngrid) ! wave saturation index REAL maxs(NW,ngrid) ! wave saturation index integer LLSATURATION(NW,ngrid) ! layer where the gravity waves break integer LLZCRITICAL(NW,ngrid) ! layer where the gravity waves depleting integer ZHSATURATION(NW,ngrid) ! altitude of the layer where the gravity waves break INTEGER ll_zb(ngrid) ! layer where the gravity waves break INTEGER ll_zc(ngrid) ! layer where the gravity waves break integer ll_zb_ii,ll_zc_ii integer ll_zb_max, ll_zb_max_r REAL d_eddy_mix_p(NW,ngrid) ! Diffusion coefficients where ll> = ll_zb REAL d_eddy_mix_m(NW,ngrid) ! Diffusion coefficients where ll< ll_zb REAL d_eddy_mix_s(NW,ngrid) ! Diffusion coefficients where ll< ll_zb REAL lambda_img(NW,ngrid) REAL d_eddy_mix_p_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll> = ll_zb REAL d_eddy_mix_m_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll< ll_zb REAL d_eddy_mix_tot_ll(nlayer,NW,ngrid) ! Diffusion coefficients where ll> = ll_zb REAL d_eddy_mix_tot(ngrid, nlayer+1) REAL d_eddy_mix(NW,ngrid) ! Comprehensive Diffusion coefficients REAL u_eddy_mix_p(NW, ngrid) ! Zonal Diffusion coefficients REAL v_eddy_mix_p(NW, ngrid) ! Meridional Diffusion coefficients REAL h_eddy_mix_p(NW, ngrid) ! potential temperature DC Real u_eddy_mix_tot(ngrid, nlayer+1) ! Total zonal D Real v_eddy_mix_tot(ngrid, nlayer+1) ! Total meridional D Real h_eddy_mix_tot(ngrid, nlayer+1) ! Total PT D REAL U_shear(ngrid,nlayer) Real wwp_vertical_tot(nlayer+1, NW, ngrid) ! Total meridional D Real wwp_vertical_ll(nlayer+1) real eddy_mix_ll(nlayer) real eddy_mix_tot_ll(nlayer) REAL pq_eddy_mix_p(NW,ngrid,nq) REAL pq_eddy_mix_tot(ngrid, nlayer+1,nq) REAL zq(ngrid,nlayer,nq) ! advected field nq REAL, save:: eff !$OMP THREADPRIVATE(eff) REAL, save:: eff1 !$OMP THREADPRIVATE(eff1) REAL, save:: vdl !$OMP THREADPRIVATE(vdl) REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP) REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level REAL wwpsat(nw,ngrid) ! Wave EP-fluxes of saturation REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP) REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP) REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW) REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW) REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0) REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable) !$OMP THREADPRIVATE(epflux_max) INTEGER LAUNCH ! Launching altitude REAL, save :: xlaunch ! Control the launching altitude by pressure !$OMP THREADPRIVATE(xlaunch) REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx. not used) REAL cmax(ngrid,nlayer) ! Max horizontal absolute phase velocity (at the maginitide of zonal wind u at the launch altitude) ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS REAL, save :: sat ! saturation parameter(tunable) !$OMP THREADPRIVATE(sat) REAL, save :: rdiss ! dissipation coefficient (tunable) !$OMP THREADPRIVATE(rdiss) REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic frequency ! 0.3.3 Background flow at 1/2 levels and vertical coordinate REAL H0bis(ngrid, nlayer) ! Characteristic Height of the atmosphere (specific locations) REAL, save:: H0 ! Characteristic Height of the atmosphere (constant) !$OMP THREADPRIVATE(H0) REAL, parameter:: pr = 250 ! Reference pressure [Pa] REAL, parameter:: tr = 220. ! Reference temperature [K] REAL ZH(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) REAL ZHbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H0bis) REAL UH(ngrid, nlayer + 1) ! zonal wind at 1/2 levels REAL VH(ngrid, nlayer + 1) ! meridional wind at 1/2 levels REAL PH(ngrid, nlayer + 1) ! Pressure at 1/2 levels REAL, parameter:: psec = 1.e-20 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure) REAL BV(ngrid, nlayer + 1) ! Brunt Vaisala freq. (N^2) at 1/2 levels REAL, parameter:: bvsec = 1.e-3 ! Security to avoid negative BV REAL HREF(nlayer + 1) ! Reference pressure for launching altitude REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 INTEGER first_breaking_flag(ngrid) INTEGER first_satuatio_flag(ngrid) LOGICAL,SAVE :: firstcall = .true. !$OMP THREADPRIVATE(firstcall) !----------------------------------------------------------------------------------------------------------------------- ! 1. INITIALISATIONS !----------------------------------------------------------------------------------------------------------------------- IF (firstcall) THEN write(*,*) "nonoro_gwd_mix: non-oro GW-induced mixing is active!" epflux_max = 5.E-4 ! Mars' value !! call getin_p("nonoro_gwd_epflux_max", epflux_max) write(*,*) "NONORO_GWD_MIX: epflux_max=", epflux_max sat = 1.5 ! default gravity waves saturation value !! call getin_p("nonoro_gwd_sat", sat) write(*,*) "NONORO_GWD_MIX: sat=", sat ! cmax = 50. ! default gravity waves phase velocity value !! ! call getin_p("nonoro_gwd_cmax", cmax) ! write(*,*) "NONORO_GWD_MIX: cmax=", cmax rdiss=0.07 ! default coefficient of dissipation !! call getin_p("nonoro_gwd_rdiss", rdiss) write(*,*) "NONORO_GWD_MIX: rdiss=", rdiss kmax=1.E-4 ! default Max horizontal wavenumber !! call getin_p("nonoro_gwd_kmax", kmax) write(*,*) "NONORO_GWD_MIX: kmax=", kmax kmin=7.E-6 ! default Max horizontal wavenumber !! call getin_p("nonoro_gwd_kmin", kmin) write(*,*) "NONORO_GWD_MIX: kmin=", kmin xlaunch=0.6 ! default GW luanch altitude controller !! call getin_p("nonoro_gwd_xlaunch", xlaunch) write(*,*) "NONORO_GWD_MIX: xlaunch=", xlaunch eff=0.1 ! Diffusion effective factor !! call getin_p("nonoro_gwimixing_eff", eff) write(*,*) "NONORO_GWD_MIX: eff=", eff eff1=0.1 ! Diffusion effective factor !! call getin_p("nonoro_gwimixing_eff1", eff1) write(*,*) "NONORO_GWD_MIX: eff1=", eff1 vdl=1.5 ! Diffusion effective factor !! call getin_p("nonoro_gwimixing_vdl", vdl) write(*,*) "NONORO_GWD_MIX: vdl=", vdl ! Characteristic vertical scale height H0 = r * tr / g ! Control if (deltat .LT. dtime) THEN call abort_physic("NONORO_GWD_MIX","gwd random: deltat lower than dtime!",1) endif if (nlayer .LT. nw) THEN call abort_physic("NONORO_GWD_MIX","gwd random: nlayer lower than nw!",1) endif firstcall = .false. ENDIF ! Compute current values of temperature and winds tt(:,:)=pt(:,:)+dtime*pdt(:,:) uu(:,:)=pu(:,:)+dtime*pdu(:,:) vv(:,:)=pv(:,:)+dtime*pdv(:,:) zq(:,:,:)=pq(:,:,:)+dtime*pdq(:,:,:) hh(:,:)=pht(:,:)+dtime*pdht(:,:) ! Compute the real mass density by rho=p/R(T)T DO ll=1,nlayer DO ii=1,ngrid rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll)) ENDDO ENDDO ! print*,'epflux_max just after firstcall:', epflux_max !----------------------------------------------------------------------------------------------------------------------- ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS !----------------------------------------------------------------------------------------------------------------------- ! Pressure and inverse of pressure at 1/2 level DO LL = 2, nlayer PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.) end DO PH(:, nlayer + 1) = 0. PH(:, 1) = 2. * PP(:, 1) - PH(:, 2) ! call write_output('nonoro_pp','nonoro_pp', 'm',PP(:,:)) ! call write_output('nonoro_ph','nonoro_ph', 'm',PH(:,:)) ! Launching level for reproductible case !Pour revenir a la version non reproductible en changeant le nombre de !process ! Reprend la formule qui calcule PH en fonction de PP=play DO LL = 2, nlayer HREF(LL) = EXP((LOG(presnivs(LL))+ LOG(presnivs(LL - 1))) / 2.) end DO HREF(nlayer + 1) = 0. HREF(1) = 2. * presnivs(1) - HREF(2) LAUNCH=0 DO LL = 1, nlayer IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL ENDDO ! Log pressure vert. coordinate ZH(:,1) = 0. DO LL = 2, nlayer + 1 !ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC)) H0bis(:, LL-1) = r * TT(:, LL-1) / g ZH(:, LL) = ZH(:, LL-1) & + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1) end DO ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC)) ! call write_output('nonoro_zh','nonoro_zh', 'm',ZH(:,2:nlayer+1)) ! Winds and Brunt Vaisala frequency DO LL = 2, nlayer UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind ! Brunt Vaisala frequency (=g/T*[dT/dz + g/cp] ) BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1))) & *((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+ & G / cpnew(:,LL)) BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small end DO BV(:, 1) = BV(:, 2) UH(:, 1) = 0. VH(:, 1) = 0. BV(:, nlayer + 1) = BV(:, nlayer) UH(:, nlayer + 1) = UU(:, nlayer) VH(:, nlayer + 1) = VV(:, nlayer) cmax(:,launch)=UU(:, launch) DO ii=1,ngrid KSTAR = PI/SQRT(cell_area(II)) MAX_K(II)=MAX(kmin,kstar) ENDDO call write_output('nonoro_bv','Brunt Vaisala frequency in nonoro', 'Hz',BV(:,2:nlayer+1)) !----------------------------------------------------------------------------------------------------------------------- ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY !----------------------------------------------------------------------------------------------------------------------- ! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way DO JW = 1, NW ! Angle DO II = 1, ngrid ! Angle (0 or PI so far) RAN_NUM_1=MOD(TT(II, JW) * 10., 1.) RAN_NUM_2= MOD(TT(II, JW) * 100., 1.) ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.)* PI / 2. ! Horizontal wavenumber amplitude ! From Venus model: TN+GG April/June 2020 (rev2381) ! "Individual waves are not supposed to occupy the entire domain, but only a fraction of it" (Lott+2012) ! ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2 KSTAR = PI/SQRT(cell_area(II)) ZK(JW, II) = MAX_K(II) + (KMAX - MAX_K(II)) *RAN_NUM_2 ! Horizontal phase speed ! this computation allows to get a gaussian distribution centered on 0 (right ?) ! then cmin is not useful, and it favors small values. CPHA(:) = 0. DO JJ = 1, NA RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.) CPHA(ii) = CPHA(ii) + 2.*CMAX(ii,launch)* & (RAN_NUM_3 -0.5)* & SQRT(3.)/SQRT(NA*1.) END DO IF (CPHA(ii).LT.0.) THEN CPHA(ii) = -1.*CPHA(ii) ZP(JW,II) = ZP(JW,II) + PI ENDIF ! otherwise, with the computation below, we get a uniform distribution between cmin and cmax. ! ran_num_3 = mod(tt(ii, jw)**2, 1.) ! cpha = cmin + (cmax - cmin) * ran_num_3 ! Intrinsic frequency ZO(JW, II) = CPHA(II) * ZK(JW, II) ! Intrinsic frequency is imposed ZO(JW, II) = ZO(JW, II) & + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) & + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH) ! Momentum flux at launch level epflux_0(JW, II) = epflux_max & * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.) ENDDO end DO ! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 !----------------------------------------------------------------------------------------------------------------------- ! 4. COMPUTE THE FLUXES !----------------------------------------------------------------------------------------------------------------------- ! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes. !------------------------------------------------------ DO JW = 1, NW ! Evaluate intrinsic frequency at launching altitude: intr_freq_p(JW, :) = ZO(JW, :) & - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) end DO ! VERSION WITHOUT CONVECTIVE SOURCE ! Vertical velocity at launch level, value to ensure the imposed ! mom flux: DO JW = 1, NW ! WW is directly a flux, here, not vertical velocity anymore WWP(JW, :) = epflux_0(JW,:) u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :) v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :) end DO ! 4.2 Initial flux at launching altitude !------------------------------------------------------ u_epflux_tot(:, LAUNCH) = 0 v_epflux_tot(:, LAUNCH) = 0 DO JW = 1, NW u_epflux_tot(:, LAUNCH) = u_epflux_tot(:, LAUNCH) + u_epflux_p(JW, :) v_epflux_tot(:, LAUNCH) = v_epflux_tot(:, LAUNCH) + v_epflux_p(JW, :) end DO ! 4.3 Loop over altitudes, with passage from one level to the next done by: !---------------------------------------------------------------------------- ! i) conserving the EP flux, ! ii) dissipating a little, ! iii) testing critical levels, ! iv) testing the breaking. !---------------------------------------------------------------------------- wwp_vertical_tot(:, :,:) =0. DO LL = LAUNCH, nlayer - 1 ! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL TO THE NEXT) DO JW = 1, NW intr_freq_m(JW, :) = intr_freq_p(JW, :) WWM(JW, :) = WWP(JW, :) ! Intrinsic Frequency intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) & - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1) ! Vertical velocity in flux formulation wwpsat(JW,:) = ABS(intr_freq_p(JW, :))**3 / (2.*BV(:, LL+1)) & * rho(:,launch)*exp(-zh(:, ll + 1) / H0) & * SAT**2 *KMIN**2 / ZK(JW, :)**4 WWP(JW, :) = MIN( & ! No breaking (Eq.6) WWM(JW, :) & ! Dissipation (Eq. 8)(real rho used here rather than pressure ! parametration because the original code has a bug if the density of ! the planet at the launch altitude not approximate 1): ! * EXP(- RDISS*2./rho(:, LL + 1) & * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 & / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 & * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), & ! Critical levels (forced to zero if intrinsic frequency ! changes sign) MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) & ! Saturation (Eq. 12) (rho at launch altitude is imposed by J.Liu. ! Same reason with Eq. 8) * WWPSAT(JW,:)) end DO ! END OF W(KB)ARNING ! first_breaking_flag(:)=0 !mixing start at first breaking ! first_satuatio_flag(:)=0 DO JW=1,NW wwp_vertical_tot(ll, JW, :) = WWP(JW,:)/rho(:,ll) ENDDO end DO ! DO LL = LAUNCH, nlayer - 1 ! print*, "this line is for tunning" DO II=1, ngrid DO JW=1,NW wwp_vertical_ll(:)=wwp_vertical_tot(:, JW, II) ll_zb_max = MAXloc(wwp_vertical_ll(:),1) !if (LLSATURATION(JW,II).ne.-1000) then LLSATURATION(JW,II)=ll_zb_max !endif ! print*, "this line is for tunning" if (ll_zb_max .gt. 1) then !ll_zb_max_r =MAXloc(wwp_vertical_ll(nlayer:1:-1),1) !if (ll_zb_max_r .gt. ll_zb_max) LLSATURATION(JW,II)=ll_zb_max_r DO ll = nlayer,ll_zb_max,-1 if (wwp_vertical_ll(ll)-wwp_vertical_ll(ll-1).gt. 0) then LLSATURATION(JW,II)=ll goto 119 endif ENDDO 119 continue endif ! if (ll_zb_max .gt. launch) then ! DO LL = (nlayer + 1), ll_zb_max,-1 ! if (abs(wwp_vertical_ll(ll)).le. 1.e-9) LLZCRITICAL(JW,II)=LL ! ENDDO ! else ! LLZCRITICAL(JW,II)=1 ! endif !print*, "this line is for tunning" ENDDO ENDDO !print*, "this line is for tunning" DO JW = 1, NW ! Evaluate intrinsic frequency at launching altitude: intr_freq_p(JW, :) = ZO(JW, :) & - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) end DO d_eddy_mix_p_ll(:,:,:)=0. d_eddy_mix_m_ll(:,:,:)=0. d_eddy_mix_p(:,:)=0. d_eddy_mix_m(:,:)=0. d_eddy_mix_s(:,:)=0. lambda_img(:,:) = 0. u_shear(:,:)= 0. DO LL = LAUNCH, nlayer - 1 U_shear(:,ll)=(UU(:, LL+1)-UU(:, LL)) /(ZH(:, LL+1) - ZH(:, LL)) DO II=1,ngrid ! all the eddy diffusion parameters are culculated at here DO JW = 1, NW intr_freq_m(JW, II) = intr_freq_p(JW, II) ! Intrinsic Frequency intr_freq_p(JW, II) = ZO(JW, II) - ZK(JW, II) * COS(ZP(JW,II)) * UH(II, LL + 1) & - ZK(JW, II) * SIN(ZP(JW,II)) * VH(II, LL + 1) ll_zb(II)= LLSATURATION(JW,II) ll_zb_ii = ll_zb(II) ! ll_zc(II)= LLZCRITICAL(JW,II) ! ll_zc_ii = ll_zc(II) ! If (ll_zb_ii.le. launch .or. ll .gt. ll_zc_ii) then If (ll .lt. ll_zb_ii .or. ll_zb_ii.lt. launch) then d_eddy_mix_p(JW,II)=0. d_eddy_mix_m(JW,II)=0. d_eddy_mix_p_ll(ll,JW,ii)=0. d_eddy_mix_m_ll(ll,JW,ii)=0. endif IF (LL.GE.ll_zb_ii .and. ll_zb_ii.ge. launch) THEN lambda_img(JW,II) = 0.5 / H0bis(II,LL) & - 1.5 /((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.) & *(intr_freq_p(JW, II) - intr_freq_m(JW, II)) & /(ZH(II, LL+1) - ZH(II, LL)) & + 1.5/((BV(II,LL)+BV(II, LL+1))/2.) & *(BV(II, LL+1)-BV(II, LL)) /(ZH(II, LL+1) - ZH(II, LL)) !lambda_img(JW,II) = max(0.5 / H0bis(II,LL), abs(lambda_img(JW,II))) if (lambda_img(JW,II).lt. 0) lambda_img(JW,II) = 0. !if (lambda_img(JW,II).gt. 0.5/H0bis(II,LL)) lambda_img(JW,II) = 0.5/H0bis(II,LL) d_eddy_mix_s(JW,II) = eff*MAX(ABS(intr_freq_p(JW, II) + intr_freq_m(JW, II)) / 2., & ZOISEC)**4 / (ZK(JW, II)**3 * ((BV(II,LL)+BV(II, LL+1))/2.)**3)& *lambda_img(JW,II) ! *abs(0.5 / H0bis(II,LL) - 1.5*ZK(JW, II) & ! /abs((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.) & ! *(UU(II, LL+1)-UU(II, LL)) /(ZH(II, LL+1) - ZH(II, LL)) & ! + 1.5/((BV(II,LL)+BV(II, LL+1))/2.) & ! *(BV(II, LL+1)-BV(II, LL)) /(ZH(II, LL+1) - ZH(II, LL))) ! *MAX(0., sign(1., lambda_img(JW,II))) d_eddy_mix_p(JW,II) = Min( d_eddy_mix_s(JW,II) , & -((wwp_vertical_tot(ll+1, JW, II)-wwp_vertical_tot(ll, JW, II)) & /(ZH(II, LL+1) - ZH(II, LL))) & *((intr_freq_p(JW, II)+ intr_freq_m(JW, II)) / 2.)*eff1 & /(((BV(II, LL+1) + BV(II,LL)) / 2.)**2 * ZK(JW, II))) if (d_eddy_mix_p(jw,ii) .lt. 0.) then print*, "this line is for tunning" endif ENDIF end DO !JW = 1, NW end DO !II=1,ngrid ! print*, "this line is for tunning" d_eddy_mix_p_ll(ll,:,:)=d_eddy_mix_p(:,:) ! print*, "this line is for tunning" ! if (ll.ge.55) then ! print*, "this line is for tunning" ! endif end DO ! DO LL = LAUNCH, nlayer - 1 call write_output('zonal_shear','u shear', 's-1',u_shear(:,:)) DO II=1,ngrid DO JW = 1, NW ll_zb(II)= LLSATURATION(JW,II) ll_zb_ii = ll_zb(II) do LL=launch, nlayer-1 IF (LL.eq.ll_zb_ii .and. ll_zb_ii .gt.launch) THEN d_eddy_mix_m(JW,II) = d_eddy_mix_p_ll(ll,JW,ii) ! if (ii.eq. 848) then ! print*, "this line is for tunning" ! endif ENDIF enddo ENDDO ENDDO DO II=1,ngrid DO JW = 1, NW ll_zb(II)= LLSATURATION(JW,II) ll_zb_ii = ll_zb(II) DO LL= launch, (ll_zb_ii - 1) if (ll_zb_ii .gt. launch) then d_eddy_mix_m_ll(ll,JW,II) = d_eddy_mix_m(JW,II) & * exp(vdl*(ZH(II, LL)-ZH(II, ll_zb_ii) ) / H0) else d_eddy_mix_m_ll(ll,JW,II) = 0. endif endDO ENDDO ENDDO ! print*, "this line is for tunning" DO II=1, ngrid DO JW=1,NW eddy_mix_ll(:) = d_eddy_mix_p_ll(:,JW,II)+ d_eddy_mix_m_ll(:,JW,II) ! print*, "this line is for tunning" enddo ENDDO DO JW = 1, NW ! Evaluate intrinsic frequency at launching altitude: intr_freq_p(JW, :) = ZO(JW, :) & - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) end DO u_eddy_mix_tot(:, :) = 0. v_eddy_mix_tot(:, :) = 0. h_eddy_mix_tot(:, :) = 0. u_eddy_mix_p(:, :)=0. v_eddy_mix_p(:, :)=0. h_eddy_mix_p(:, :)=0. d_eddy_mix_tot_ll(:,:,:)=0. pq_eddy_mix_p(:,:,:)=0. pq_eddy_mix_tot(:, :,:)=0. d_eddy_mix(:,:)=0. d_eddy_mix_p_ll(nlayer,:,:)=d_eddy_mix_p_ll(nlayer-1,:,:) d_eddy_mix_tot(:, :) =0. DO LL = LAUNCH, nlayer - 1 d_eddy_mix(:,:) = d_eddy_mix_m_ll(ll,:,:) + d_eddy_mix_p_ll(ll,:,:) DO JW = 1, NW u_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(UU(:, LL + 1) - UU(:, LL)) & /(ZH(:, LL + 1) - ZH(:, LL)) & *SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) v_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(VV(:, LL + 1) - VV(:, LL)) & /(ZH(:, LL + 1) - ZH(:, LL)) & *SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) h_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(HH(:, LL + 1) - HH(:, LL)) & /(ZH(:, LL + 1) - ZH(:, LL)) & *SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) ENDDO u_eddy_mix_tot(:, LL+1)=0. v_eddy_mix_tot(:, LL+1)=0. h_eddy_mix_tot(:, LL+1)=0. DO JW=1,NW u_eddy_mix_tot(:, LL+1) = u_eddy_mix_tot(:, LL+1) + u_eddy_mix_p(JW, :) v_eddy_mix_tot(:, LL+1) = v_eddy_mix_tot(:, LL+1) + v_eddy_mix_p(JW, :) h_eddy_mix_tot(:, LL+1) = h_eddy_mix_tot(:, LL+1) + h_eddy_mix_p(JW, :) ENDDO DO JW=1,NW ! d_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL+1) + d_eddy_mix(JW,:) d_eddy_mix_tot(:, LL+1) = d_eddy_mix_tot(:, LL+1) + d_eddy_mix(JW,:) ! u_eddy_mix_tot(:, LL+1) = u_eddy_mix_tot(:, LL+1)+ u_eddy_mix_p(JW, :) ! v_eddy_mix_tot(:, LL+1) = v_eddy_mix_tot(:, LL+1)+ v_eddy_mix_p(JW, :) ENDDO ! if (ll.ge.55) then ! print*, "this line is for tunning" !endif DO QQ=1,NQ DO JW=1,NW pq_eddy_mix_p(JW, :, QQ) = d_eddy_mix(JW, :)* (zq(:, LL + 1,QQ)- zq(:, LL, QQ)) & /(ZH(:, LL + 1)- ZH(:, LL)) & *SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) ! pq_eddy_mix_tot(:, LL+1,QQ) = pq_eddy_mix_tot(:, LL+1,QQ) & ! + pq_eddy_mix_p(JW, :, QQ) ENDDO ENDDO pq_eddy_mix_tot(:, LL+1,:)=0. DO JW=1,NW DO QQ=1, NQ pq_eddy_mix_tot(:, LL+1,QQ) = pq_eddy_mix_tot(:, LL+1,QQ) + pq_eddy_mix_p(JW, :, QQ) ENDDO endDO ! print*, "this line is for tunning" ENDDO !LL = LAUNCH, nlayer - 1 d_eddy_mix_tot(:, LAUNCH) = d_eddy_mix_tot(:, LAUNCH+1) d_eddy_mix_tot(:, nlayer + 1) = d_eddy_mix_tot(:, nlayer) d_eddy_mix_tot(:,:) = DTIME/DELTAT/REAL(NW) * d_eddy_mix_tot(:,:) & + (1.-DTIME/DELTAT) * de_eddymix_rto(:,:) call write_output('nonoro_d_mixing_tot','Total EP Flux along U in nonoro', 'm2s-1',d_eddy_mix_tot(:,2:nlayer+1)) ! u_eddy_mix_tot(:, :) = 0. ! v_eddy_mix_tot(:, :) = 0. ! pq_eddy_mix_tot(:, :, :) = 0. ! DO LL = 1, nlayer-1 !u_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL)*(UU(:, LL + 1) - UU(:, LL)) & ! /(ZH(:, LL + 1) - ZH(:, LL)) !*sign(1.,u_shear(:, LL)) !v_eddy_mix_tot(:, LL) = d_eddy_mix_tot(:, LL)*(VV(:, LL + 1) - VV(:, LL)) & ! ! /(ZH(:, LL + 1) - ZH(:, LL)) !*sign(1.,u_shear(:, LL)) ! DO QQ=1, NQ ! pq_eddy_mix_tot(:, LL,QQ) = d_eddy_mix_tot(:, LL) & ! * (zq(:, LL + 1,QQ)- zq(:, LL, QQ)) & ! /(ZH(:, LL + 1)- ZH(:, LL)) ! ENDDO !ENDDO !LL = LAUNCH, nlayer - 1 de_eddymix_rto(:,:) = d_eddy_mix_tot(:,:) !----------------------------------------------------------------------------------------------------------------------- ! 5. TENDENCY CALCULATIONS !----------------------------------------------------------------------------------------------------------------------- ! 5.1 Flow rectification at the top and in the low layers ! -------------------------------------------------------- ! Warning, this is the total on all GW u_eddy_mix_tot(:, nlayer + 1) = 0. v_eddy_mix_tot(:, nlayer + 1) = 0. h_eddy_mix_tot(:, nlayer + 1) = 0. pq_eddy_mix_tot(:, nlayer + 1,:)=0. ! Here, big change compared to FLott version: ! We compensate (u_epflux_tot(:, LAUNCH), ie total emitted upward flux ! over the layers max(1,LAUNCH-3) to LAUNCH-1 DO LL = 1, max(1,LAUNCH-3) u_eddy_mix_tot(:, LL) = 0. v_eddy_mix_tot(:, LL) = 0. h_eddy_mix_tot(:, LL) = 0. end DO DO QQ=1,NQ DO LL = 1, max(1,LAUNCH-3) pq_eddy_mix_tot(:, LL,QQ) = 0. end DO ENDDO DO LL = max(2,LAUNCH-2), LAUNCH-1 u_eddy_mix_tot(:, LL) = u_eddy_mix_tot(:, LL - 1) + u_eddy_mix_tot(:, LAUNCH) & * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) v_eddy_mix_tot(:, LL) = v_eddy_mix_tot(:, LL - 1) + v_eddy_mix_tot(:, LAUNCH) & * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) h_eddy_mix_tot(:, LL) = h_eddy_mix_tot(:, LL - 1) + h_eddy_mix_tot(:, LAUNCH) & * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) end DO DO QQ=1,NQ DO LL = max(2,LAUNCH-2), LAUNCH-1 pq_eddy_mix_tot(:, LL,QQ) = pq_eddy_mix_tot(:, LL-1,QQ)+pq_eddy_mix_tot(:, LAUNCH,QQ)& * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) end DO ENDDO !u_eddy_mix_tot(:,:) = DTIME/DELTAT/REAL(NW) * u_eddy_mix_tot(:,:) & ! + (1.-DTIME/DELTAT) * df_eddymix_flx(:,:) !do ii=1,ngrid ! if (u_eddy_mix_tot(ii,58).lt. -8000.) then ! print*, ii ! endif !enddo ! This way, the total flux from GW is zero, but there is a net transport ! (upward) that should be compensated by circulation ! and induce additional friction at the surface call write_output('nonoro_u_mixing_tot','Total EP Flux along U in nonoro', '',u_eddy_mix_tot(:,2:nlayer+1)) call write_output('nonoro_v_mixing_tot','Total EP Flux along V in nonoro', '',v_eddy_mix_tot(:,2:nlayer+1)) ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 !--------------------------------------------- DO LL = 1, nlayer !Notice here the d_u and d_v are tendency (For Mars) but not increment(Venus). d_u(:, LL) = (u_eddy_mix_tot(:, LL + 1) - u_eddy_mix_tot(:, LL)) & / (ZH(:, LL + 1) - ZH(:, LL)) !d_v(:, LL) = (v_eddy_mix_tot(:, LL + 1) - v_eddy_mix_tot(:, LL)) & ! / (ZH(:, LL + 1) - ZH(:, LL)) d_h(:, LL) = (h_eddy_mix_tot(:, LL + 1) - h_eddy_mix_tot(:, LL)) & / (ZH(:, LL + 1) - ZH(:, LL)) ENDDO DO QQ=1,NQ DO ll = 1, nlayer d_pq(:, ll, QQ) = (pq_eddy_mix_tot(:, LL + 1,QQ) - pq_eddy_mix_tot(:, LL, QQ)) & / (ZH(:, LL + 1) - ZH(:, LL)) end DO ENDDO !df_eddymix_flx(:,:) = u_eddy_mix_tot(:,:) !d_pq(:, :, :)=0. !d_t(:,:) = 0. !d_v(:,:) = 0. !zustr(:) = 0. !zvstr(:) = 0. ! call write_output('nonoro_d_u','nonoro_d_u', '',d_u(:,:)) ! call write_output('nonoro_d_v','nonoro_d_v', '',d_v(:,:)) ! 5.3 Update tendency of wind with the previous (and saved) values !----------------------------------------------------------------- d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) & + (1.-DTIME/DELTAT) * du_eddymix_gwd(:,:) !d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & ! + (1.-DTIME/DELTAT) * dv_eddymix_gwd(:,:) du_eddymix_gwd(:,:) = d_u(:,:) !dv_eddymix_gwd(:,:) = d_v(:,:) d_v(:,:)=0. call write_output('du_eddymix_gwd','Tendency on U due to nonoro GW', 'm.s-2',du_eddymix_gwd(:,:)) !call write_output('dv_eddymix_gwd','Tendency on V due to nonoro GW', 'm.s-2',dv_eddymix_gwd(:,:)) d_h(:,:) = DTIME/DELTAT/REAL(NW) * d_h(:,:) & + (1.-DTIME/DELTAT) * dh_eddymix_gwd(:,:) do ii=1,ngrid d_t(ii,:) = d_h(ii,:) * (PP(ii,:) / PH(ii,1))**rcp enddo dh_eddymix_gwd(:,:)=d_h(:,:) DO QQ=1,NQ d_pq(:, :, QQ) =DTIME/DELTAT/REAL(NW) * d_pq(:, :, QQ) & + (1.-DTIME/DELTAT) * dq_eddymix_gwd(:, :, QQ) endDO END SUBROUTINE NONORO_GWD_MIX ! ======================================================== ! Subroutines used to allocate/deallocate module variables ! ======================================================== SUBROUTINE ini_nonoro_gwd_mix(ngrid,nlayer,nq) IMPLICIT NONE INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers INTEGER, INTENT (in) :: nq ! number of atmospheric tracers allocate(du_eddymix_gwd(ngrid,nlayer)) allocate(dv_eddymix_gwd(ngrid,nlayer)) allocate(dh_eddymix_gwd(ngrid,nlayer)) allocate(dq_eddymix_gwd(ngrid,nlayer,nq)) allocate(de_eddymix_rto(ngrid,nlayer+1)) allocate(df_eddymix_flx(ngrid,nlayer+1)) !du_eddymix_gwd(:,:)=0 !dv_eddymix_gwd(:,:)=0 ! allocate(east_gwstress(ngrid,nlayer)) ! east_gwstress(:,:)=0 ! allocate(west_gwstress(ngrid,nlayer)) ! west_gwstress(:,:)=0 END SUBROUTINE ini_nonoro_gwd_mix ! ---------------------------------- SUBROUTINE end_nonoro_gwd_mix IMPLICIT NONE if (allocated(du_eddymix_gwd)) deallocate(du_eddymix_gwd) if (allocated(dv_eddymix_gwd)) deallocate(dv_eddymix_gwd) if (allocated(dh_eddymix_gwd)) deallocate(dh_eddymix_gwd) if (allocated(dq_eddymix_gwd)) deallocate(dq_eddymix_gwd) if (allocated(de_eddymix_rto)) deallocate(de_eddymix_rto) if (allocated(df_eddymix_flx)) deallocate(df_eddymix_flx) ! if (allocated(east_gwstress)) deallocate(east_gwstress) ! if (allocated(west_gwstress)) deallocate(west_gwstress) END SUBROUTINE end_nonoro_gwd_mix !----------------------------------- ! FUNCTION MAXLOCATION(gwd_normal,gwd_satura) ! implicit none ! INTEGER MAXLOCATION ! REAL gwd_normal,gwd_satura ! IF (gwd_normal .GT. gwd_satura ) THEN ! MAXLOCATION=1 ! ELSEIF (gwd_normal .LT.gwd_satura) THEN ! MAXLOCATION=2 ! ELSE ! MAXLOCATION=1 ! ENDIF ! return ! END FUNCTION END MODULE nonoro_gwd_mix_mod