Index: trunk/LMDZ.PLUTO/changelog.txt =================================================================== --- trunk/LMDZ.PLUTO/changelog.txt (revision 3934) +++ trunk/LMDZ.PLUTO/changelog.txt (revision 3935) @@ -1879,2 +1879,6 @@ flag to trigger outputing a diagsoil.nc file or not. +== 24/10/2025 == Jliu +Add non-orographic gravity waves and induced turbulence to the model. +The model is compileable with the two new shcemes, however, carefull +parameter selection is needed since the current parameters are for Mars. Index: trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3934) +++ trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3935) @@ -6,4 +6,6 @@ logical,save :: calladj,calltherm,n2cond,callsoil,calllott !$OMP THREADPRIVATE(calladj,calltherm,n2cond,callsoil,calllott) + logical,save :: calllott_nonoro +!$OMP THREADPRIVATE(calllott_nonoro) logical,save :: callconduct,callmolvis,callmoldiff !$OMP THREADPRIVATE(callconduct,callmolvis,callmoldiff) Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3934) +++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3935) @@ -23,5 +23,11 @@ mugaz, pi, avocado use planetwide_mod, only: planetwide_sumval + use nonoro_gwd_mix_mod, only: calljliu_gwimix use callkeys_mod + use nonoro_gwd_ran_mod, only: ini_nonoro_gwd_ran, & + end_nonoro_gwd_ran + use nonoro_gwd_mix_mod, only: ini_nonoro_gwd_mix, & + end_nonoro_gwd_mix + use surfdat_h use wstats_mod, only: callstats @@ -373,7 +379,23 @@ if (is_master) write(*,*) trim(rname)//": callsoil = ",callsoil + if (is_master) write(*,*)trim(rname)//& + ": call orographic gravity waves ?" calllott=.false. ! default value call getin_p("calllott",calllott) if (is_master) write(*,*)trim(rname)//": calllott = ",calllott + + if (is_master) write(*,*)trim(rname)//& + ": call non-orographic gravity waves ?" + calllott_nonoro=.false. ! default value + call getin_p("calllott_nonoro",calllott_nonoro) + if (is_master) write(*,*)trim(rname)//& + ": calllott_nonoro = ",calllott_nonoro + + if (is_master) write(*,*)trim(rname)//& + ": call jliu's non-orogrphic GW-induced turbulence" + calljliu_gwimix=.false. ! default value + call getin_p("calljliu_gwimix",calljliu_gwimix) + if (is_master) write(*,*)trim(rname)//& + ": calljliu_gwimix = ",calljliu_gwimix if (is_master) write(*,*)trim(rname)//& @@ -1500,4 +1522,12 @@ call ini_comsoil_h(ngrid) + ! allocate arrays in "nonoro_gwd_ran_mod" + call end_nonoro_gwd_ran + call ini_nonoro_gwd_ran(ngrid,nlayer) + + ! allocate arrays in "nonoro_gwd_mix_mod" + call end_nonoro_gwd_mix + call ini_nonoro_gwd_mix(ngrid,nlayer,nq) + END SUBROUTINE inifis Index: trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_mix_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_mix_mod.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_mix_mod.F90 (revision 3935) @@ -0,0 +1,920 @@ +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_mod, 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 + + implicit none + + 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 d_wave(NW,ngrid) ! coefficients consider the tracers' gradients + 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 zq_var(ngrid,nlayer,nq) + REAL dzq_var(ngrid,nlayer,nq) + REAL mdzq_var(nlayer,nq) + REAL zq_ave(nlayer,nq) + REAL dzq_ave(nlayer,nq) + REAL zq_ratio(nlayer,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(:,:) + +! tracer average and gradients for mixing + zq_ave(:,:)=0. + zq_ave(:,:) = SUM(zq(:,:,:), dim=1) / real(ngrid) + + zq_var(:,:,:)=0. + DO ii=1,ngrid + zq_var(ii,:,:)=zq(ii,:,:)-zq_ave(:,:) + endDO + + dzq_var(:,:,:)=0. + mdzq_var(:,:) =0. + dzq_ave(:,:) =0. + zq_ratio(:,:) =0. + DO LL=1,nlayer-1 + dzq_var(:, LL, :) = (zq_var(:, LL+1, :) - zq_var(:, LL, :))**2. + mdzq_var(LL, :) = SUM(dzq_var(:, LL, :), dim=1)/real(ngrid) + dzq_ave(LL, :) = (zq_ave(LL+1, :) - zq_ave(LL, :))**2. + zq_ratio(LL,:) = MAX(0., (mdzq_var(LL+1,:) + mdzq_var(LL,:)))& + /MAX(1E-15, (dzq_ave(LL+1,:) + dzq_ave(LL,:))) + ! do qq=1,nq + ! print*, 'ratio=', zq_ratio(LL,QQ) + ! enddo + endDO + dzq_var(:,nlayer,:) = dzq_var(:,nlayer-1, :) + mdzq_var(nlayer,:) = mdzq_var(nlayer-1, :) + dzq_ave(nlayer,:) = dzq_ave(nlayer-1, :) +! 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_wave(:, :) =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 + ! d_wave(JW,:) = d_eddy_mix(JW, :)*(mdzq_var(LL,QQ)/dzq_ave(LL,QQ))**2. + DO QQ=1,NQ + DO JW=1,NW + d_wave(JW, :) = d_eddy_mix(JW, :)*zq_ratio(LL,QQ) + pq_eddy_mix_p(JW, :, QQ) = d_wave(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 + + do QQ=1,NQ + dq_eddymix_gwd(:, :, QQ)=d_pq(:, :, 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 Index: trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_ran_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_ran_mod.F90 (revision 3935) +++ trunk/LMDZ.PLUTO/libf/phypluto/nonoro_gwd_ran_mod.F90 (revision 3935) @@ -0,0 +1,513 @@ +MODULE nonoro_gwd_ran_mod + +IMPLICIT NONE + +REAL,allocatable,save :: du_nonoro_gwd(:,:) ! Zonal wind tendency due to GWD +REAL,allocatable,save :: dv_nonoro_gwd(:,:) ! Meridional wind tendency due to GWD +REAL,ALLOCATABLE,SAVE :: east_gwstress(:,:) ! Profile of eastward stress +REAL,ALLOCATABLE,SAVE :: west_gwstress(:,:) ! Profile of westward stress + +!$OMP THREADPRIVATE(du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress) + +CONTAINS + + SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, cpnew, rnew, pp, & + zmax_therm, pt, pu, pv, pdt, pdu, pdv, & + zustr,zvstr,d_t, d_u, d_v) + + !-------------------------------------------------------------------------------- + ! Parametrization of the momentum flux deposition due to a discrete + ! number of gravity waves. + ! F. Lott + ! Version 14, Gaussian distribution of the source + ! LMDz model online version + ! ADAPTED FOR VENUS / F. LOTT + S. LEBONNOIS + ! Version adapted on 03/04/2013: + ! - input flux compensated in the deepest layers + ! + ! ADAPTED FOR MARS G.GILLI 02/2016 + ! Revision with F.Forget 06/2016 Variable EP-flux according to + ! PBL variation (max velocity thermals) + ! UPDATED D.BARDET 01/2020 - reproductibility of the + ! launching altitude calculation + ! - wave characteristic + ! calculation using MOD + ! - adding east_gwstress and + ! west_gwstress variables + ! UPDATED J.LIU 12/2021 The rho (density) at the specific + ! locations is introduced. The equation + ! of EP-flux is corrected by adding the + ! term of density at launch (source) + ! altitude. + ! + !--------------------------------------------------------------------------------- + + use comcstfi_mod, only: g, pi, r + 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 + + implicit none + + CHARACTER (LEN=20) :: modname='nonoro_gwd_ran' + 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 + 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) :: 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) :: 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 + + ! 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 :: BVLOW(ngrid) ! initial N at each grid (not used) + + INTEGER II, JJ, LL + + ! 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 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 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-5 ! Security to avoid negative BV + REAL HREF(nlayer + 1) ! Reference pressure for launching altitude + + + REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 + + + LOGICAL,SAVE :: firstcall = .true. +!$OMP THREADPRIVATE(firstcall) + +!----------------------------------------------------------------------------------------------------------------------- +! 1. INITIALISATIONS +!----------------------------------------------------------------------------------------------------------------------- + IF (firstcall) THEN + write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" + epflux_max = 5.E-4 ! Mars' value !! + call getin_p("nonoro_gwd_epflux_max", epflux_max) + write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max + sat = 1.5 ! default gravity waves saturation value !! + call getin_p("nonoro_gwd_sat", sat) + write(*,*) "nonoro_gwd_ran: sat=", sat +! cmax = 50. ! default gravity waves phase velocity value !! +! call getin_p("nonoro_gwd_cmax", cmax) +! write(*,*) "nonoro_gwd_ran: cmax=", cmax + rdiss=0.07 ! default coefficient of dissipation !! + call getin_p("nonoro_gwd_rdiss", rdiss) + write(*,*) "nonoro_gwd_ran: rdiss=", rdiss + kmax=1.E-4 ! default Max horizontal wavenumber !! + call getin_p("nonoro_gwd_kmax", kmax) + write(*,*) "nonoro_gwd_ran: kmax=", kmax + kmin=7.E-6 ! default Max horizontal wavenumber !! + call getin_p("nonoro_gwd_kmin", kmin) + write(*,*) "nonoro_gwd_ran: kmin=", kmin + xlaunch=0.6 ! default GW luanch altitude controller !! + call getin_p("nonoro_gwd_xlaunch", xlaunch) + write(*,*) "nonoro_gwd_ran: xlaunch=", xlaunch + ! Characteristic vertical scale height + H0 = r * tr / g + ! Control + if (deltat .LT. dtime) THEN + call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1) + endif + if (nlayer .LT. nw) THEN + call abort_physic("nonoro_gwd_ran","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(:,:) +! 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. + !---------------------------------------------------------------------------- + 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 + 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) + * ABS(intr_freq_p(JW, :))**3 / (2.*BV(:, LL+1)) & + * rho(:,launch)*exp(-zh(:, ll + 1) / H0) & + * SAT**2 *KMIN**2 / ZK(JW, :)**4) + end DO + ! END OF W(KB)ARNING + + ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress + DO JW = 1, NW + u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) * WWP(JW, :) + v_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :)) * WWP(JW, :) + end DO + + u_epflux_tot(:, LL + 1) = 0. + v_epflux_tot(:, LL + 1) = 0. + DO JW = 1, NW + u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :) + v_epflux_tot(:, LL + 1) = v_epflux_tot(:, LL + 1) + v_epflux_p(JW, :) + EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW) + WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW) + end DO + end DO ! DO LL = LAUNCH, nlayer - 1 + +!----------------------------------------------------------------------------------------------------------------------- +! 5. TENDENCY CALCULATIONS +!----------------------------------------------------------------------------------------------------------------------- + + ! 5.1 Flow rectification at the top and in the low layers + ! -------------------------------------------------------- + ! Warning, this is the total on all GW + u_epflux_tot(:, nlayer + 1) = 0. + v_epflux_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_epflux_tot(:, LL) = 0. + v_epflux_tot(:, LL) = 0. + end DO + DO LL = max(2,LAUNCH-2), LAUNCH-1 + u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * & + (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * & + (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + & + EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ & + (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + & + WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ & + (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) + end DO + ! 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_epflux_tot','Total EP Flux along U in nonoro', '',u_epflux_tot(:,2:nlayer+1)) + call write_output('nonoro_v_epflux_tot','Total EP Flux along V in nonoro', '',v_epflux_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) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) & + / (PH(:, LL + 1) - PH(:, LL)) + d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) & + / (PH(:, LL + 1) - PH(:, LL)) + ENDDO + d_t(:,:) = 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_nonoro_gwd(:,:) + d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & + + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:) + du_nonoro_gwd(:,:) = d_u(:,:) + dv_nonoro_gwd(:,:) = d_v(:,:) + + call write_output('du_nonoro_gwd','Tendency on U due to nonoro GW', 'm.s-2',du_nonoro_gwd(:,:)) + call write_output('dv_nonoro_gwd','Tendency on V due to nonoro GW', 'm.s-2',dv_nonoro_gwd(:,:)) + + ! Cosmetic: evaluation of the cumulated stress + ZUSTR(:) = 0. + ZVSTR(:) = 0. + DO LL = 1, nlayer + ZUSTR(:) = ZUSTR(:) + D_U(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME + ZVSTR(:) = ZVSTR(:) + D_V(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME + ENDDO + + END SUBROUTINE NONORO_GWD_RAN + + + +! ======================================================== +! Subroutines used to allocate/deallocate module variables +! ======================================================== + SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer) + + IMPLICIT NONE + + INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns + INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers + + allocate(du_nonoro_gwd(ngrid,nlayer)) + allocate(dv_nonoro_gwd(ngrid,nlayer)) + allocate(east_gwstress(ngrid,nlayer)) + east_gwstress(:,:)=0 + allocate(west_gwstress(ngrid,nlayer)) + west_gwstress(:,:)=0 + + END SUBROUTINE ini_nonoro_gwd_ran +! ---------------------------------- + SUBROUTINE end_nonoro_gwd_ran + + IMPLICIT NONE + + if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd) + if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd) + if (allocated(east_gwstress)) deallocate(east_gwstress) + if (allocated(west_gwstress)) deallocate(west_gwstress) + + END SUBROUTINE end_nonoro_gwd_ran + +END MODULE nonoro_gwd_ran_mod Index: trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90 (revision 3934) +++ trunk/LMDZ.PLUTO/libf/phypluto/phyetat0_mod.F90 (revision 3935) @@ -27,4 +27,7 @@ get_field, get_var, inquire_field, & inquire_dimension, inquire_dimension_length + use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd + use nonoro_gwd_mix_mod, only: du_eddymix_gwd, dv_eddymix_gwd, de_eddymix_rto, & + df_eddymix_flx, dh_eddymix_gwd, dq_eddymix_gwd use callkeys_mod, only: surfalbedo,surfemis, callsoil use mod_phys_lmdz_para, only : is_master Index: trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3934) +++ trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (revision 3935) @@ -47,8 +47,11 @@ use calldrag_noro_mod, only: calldrag_noro use time_phylmdz_mod, only: daysec + use nonoro_gwd_ran_mod, only: nonoro_gwd_ran + use nonoro_gwd_mix_mod, only: nonoro_gwd_mix, calljliu_gwimix + use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, & calllott, callrad, callsoil, nosurf, & callconduct,callmolvis,callmoldiff, & - corrk, & + corrk, calllott_nonoro, & diurnal, enertest, fat1au, & icetstep, intheat, iradia, kastprof, & @@ -67,4 +70,5 @@ global1d, szangle, & callmufi, evol1d + use check_fields_mod, only: check_physics_fields use conc_mod, only: rnew, cpnew, ini_conc_mod @@ -476,4 +480,6 @@ REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer) REAL d_t_hin(ngrid,nlayer) + REAL d_u_mix(ngrid,nlayer), d_v_mix(ngrid,nlayer) + REAL d_t_mix(ngrid,nlayer), zdq_mix(ngrid,nlayer,nq) ! Diagnostics 2D of gw_nonoro REAL zustrhi(ngrid), zvstrhi(ngrid) @@ -1384,4 +1390,50 @@ endif ! end of 'calladj' +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! IV.b Non-orographic gravity waves and induced turbulence +! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (calllott_nonoro) THEN + CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, & + cpnew,rnew, & + zplay, & + zmax_th, & + pt, pu, pv, & + pdt, pdu, pdv, & + zustrhi,zvstrhi, & + d_t_hin, d_u_hin, d_v_hin) + IF (calljliu_gwimix) THEN + CALL nonoro_gwd_mix(ngrid,nlayer,ptimestep, & + nq,cpnew, rnew, & + zplay, & + zmax_th, & + pt, pu, pv, pq, zh, & + !loss, chemical reaction loss rates + pdt, pdu, pdv, pdq, zdh, & + ! zustrhi,zvstrhi, + zdq_mix, d_t_mix, d_u_mix, d_v_mix) + ENDIF + + ! Update tendencies + pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer) + & + d_t_hin(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer) + & + d_u_hin(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer) + & + d_v_hin(1:ngrid,1:nlayer) + ! Update tendencies of gw mixing + IF (calljliu_gwimix) THEN + pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer) + & + d_t_mix(1:ngrid,1:nlayer) + pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer) + & + d_u_mix(1:ngrid,1:nlayer) + pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer) + & + d_v_mix(1:ngrid,1:nlayer) + pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq) + & + zdq_mix(1:ngrid,1:nlayer,1:nq) + ENDIF + + + ENDIF ! of IF (calllott_nonoro) + !----------------------------------------------- ! V. Nitrogen condensation-sublimation :