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(:, :) ! Eastward stress profile REAL, allocatable, save :: west_gwstress(:, :) ! Westward stress profile !$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) ! D.BARDET 01/2020 - reproductibility of ! the launching altitude ! calculation ! - wave characteristic ! calculation using MOD ! - adding east_gwstress ! and west_gwstress variables ! ! ADAPTED FOR GENERIC D.BARDET 01/2020 ! 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(level). Bugs in BV,d_u,d_v fixed. !--------------------------------------------------------------------------------- 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 #ifdef CPP_XIOS use xios_output_mod, only: send_xios_field #endif implicit none ! 0. DECLARATIONS: ! 0.1 INPUTS INTEGER, intent(in):: ngrid ! number of atmospheric columns INTEGER, intent(in):: nlayer ! number of atmospheric columns REAL, intent(in):: dtime ! Time step of the Physics 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) ! Temperature at full levels(K) REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels(m/s) REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind 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 wind at full levels REAL :: bvlow(ngrid) ! initial N at each grid (not used) REAL :: dz ! depth of the GW source if gwd_convective_source=.true. INTEGER ii, jj, ll ! 0.3.0 Time scale of the like cycle of the waves parametrized 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 REAL kstar ! Control value to constrain the min horizontal wavenumber by the horizontal resolution REAL, save :: kmax ! Max horizontal wavenumber=N/u, u(=30) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax=62.8 km !$OMP THREADPRIVATE(kmax) !REAL, parameter:: kmax = 4.e-5 ! Max horizontal wavenumber=N/u, u(=70) zonal wind velocity at launch,lambda min,lambda=2*pi/kmax REAL, parameter:: kmin = 2.e-6 ! Min horizontal wavenumber=1/sqrt(DxDy) Dx and Dy horizontal grid,lambda max,lambda=2*pi/kmax=314 km REAL, save :: cmax ! Max horizontal absolute phase velocity=zonal wind at launch !$OMP THREADPRIVATE(cmax) !REAL, parameter:: cmax = 100. ! Test for Saturn: Max horizontal absolute phase velocity !REAL, parameter:: cmax = 70. ! Test for Saturn: Max horizontal absolute phase velocity REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity(not used) REAL cpha ! 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, parameter:: xlaunch = 0.2 ! Control the launching altitude by pressure REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer in m (approx.) REAL prec(ngrid) ! precipitations REAL prmax ! Max value of precipitation ! 0.3.2 Parameters of waves dissipations REAL, save :: sat ! saturation parameter(tunable) !$OMP THREADPRIVATE(sat) REAL, save :: rdiss ! coefficient of dissipation(tunable) !$OMP THREADPRIVATE(rdiss) REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic freguency ! 0.3.3 Background flow at 1/2 levels and vertical coordinate REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere REAL, save:: H0 ! characteristic height of the atmosphere !$OMP THREADPRIVATE(H0) !REAL, parameter:: pr = 250 ! Reference pressure [Pa] !REAL, parameter:: tr = 220. ! Reference temperature [K] ! TEST FOR SATURN STRATOSPHERE REAL, parameter:: pr = 110 ! Reference pressure [Pa] REAL, parameter:: tr = 130. ! Reference temperature [K] REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H) 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-12 ! 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. at 1/2 levels REAL, parameter:: bvsec = 1.e-6 ! Security to avoid negative BV (!!IMPORTANT: tunable, depends on which Planet) REAL href(nlayer + 1) ! Reference altitude for launching altitude REAL ran_num_1, ran_num_2, ran_num_3 LOGICAL, save :: firstcall = .true. !$OMP THREADPRIVATE(firstcall) LOGICAL, save :: gwd_convective_source = .false. !$OMP THREADPRIVATE(gwd_convective_source) !-------------------------------------------------------------------------------------------------- ! 1. INITIALISATION !-------------------------------------------------------------------------------------------------- IF (firstcall) THEN write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" epflux_max = 0.0 ! Default value !! call getin_p("nonoro_gwd_epflux_max", epflux_max) write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max sat = 1. ! default gravity waves saturation value !! call getin_p("nonoro_gwd_sat", sat) write(*,*) "nonoro_gwd_ran: sat=", sat cmax = 30. ! default gravity waves phase velocity value !! call getin_p("nonoro_gwd_cmax", cmax) write(*,*) "nonoro_gwd_ran: cmax=", cmax rdiss=1 ! 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 ! 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 gwd_convective_source = .false. ! Compute subroutine's 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 Inv of pressure, temperature 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) ! Launching altitude for reproductible case 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 end DO ! Log-pressure vertical coordinate !DO ll = 1, nlayer + 1 ! zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec)) !end DO 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)) ! 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 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) !----------------------------------------------------------------------------------------------------------------------- ! 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 ran_num_1 = mod(tt(ii, jw) * 10., 1.) ran_num_2 = mod(tt(ii, jw) * 100., 1.) ! Angle (0 or pi so far) zp(jw, ii) = (sign(1., 0.5 - ran_num_1) + 1.) * pi / 2. ! Horizontal wavenumber amplitude ! TN+GG April/June 2020 (rev2381) ! "Individual waves are not supposed to occupy the entire domain, but only a faction of it" Lott+2012 !zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 kstar = pi / sqrt(cell_area(ii)) zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * 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 = cpha + 2. * cmax * & (ran_num_3 - 0.5) * & sqrt(3.) / sqrt(na * 1.) end DO IF (cpha < 0.) THEN cpha = - 1. * cpha 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 * 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 !DO jw = 1, nw ! print*,'epflux_max just after waves charac. ramdon:', epflux_max ! 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 altutide: intr_freq_p(jw, :) = zo(jw, :) & - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) & - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) end DO IF (gwd_convective_source) THEN DO jw = 1, nw ! VERSION WITH CONVECTIVE SOURCE ! designed for Earth ! Vertical velocity at launch level, value to ensure the ! imposed mmt flux factor related to the convective forcing: ! precipitations. ! tanh limitation to values above prmax: ! WWP(JW, :) = epflux_0(JW, :) & ! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2 ! Here, we neglected the kinetic energy providing of the thermodynamic ! phase change ! ! Factor related to the characteristics of the waves: wwp(jw, :) = wwp(jw, :) * zk(jw, :)**3 / kmin / bvlow(:) & / MAX(ABS(intr_freq_p(jw, :)), zoisec)**3 ! Moderation by the depth of the source (dz here): wwp(jw, :) = wwp(jw, :) & * exp(- bvlow(:)**2 / max(abs(intr_freq_p(jw, :)), zoisec)**2 & * zk(jw, :)**2 * dz**2) ! Put the stress in the right direction: u_epflux_p(jw, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & * bv(:, launch) * cos(zp(jw, :)) * wwp(jw, :)**2 v_epflux_p(JW, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & * bv(:, launch) * sin(zp(jw, :)) * wwp(jw, :)**2 end DO ELSE ! 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 ENDIF ! 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 ! Warning! 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.0* 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, ll)) / float(nw) west_gwstress(:, ll) = west_gwstress(:, ll) & + min(0., u_epflux_p(jw, ll)) / 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_(:, 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 ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 !--------------------------------------------- DO LL = 1, nlayer !Notice here the d_u and d_v are tendency (i.e. -1/rho*dEP/dz For Mars) but not increment(Venus). d_u(:, LL) = - (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) & / (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL))) d_v(:, LL) = - (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) & / (rho(:,ll) * (ZH(:, LL + 1) - ZH(:, LL))) ENDDO d_t(:,:) = 0. ! 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(:,:) ! Cosmetic: evaluation of the cumulated stress ZUSTR(:) = 0. ZVSTR(:) = 0. DO LL = 1, nlayer ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL)) ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL)) ENDDO #ifdef CPP_XIOS call send_xios_field("du_nonoro", d_u) call send_xios_field("dv_nonoro", d_v) call send_xios_field("east_gwstress", east_gwstress) call send_xios_field("west_gwstress", west_gwstress) #endif END SUBROUTINE nonoro_gwd_ran ! =================================================================== ! Subroutines used to write variables of memory in start files ! =================================================================== 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)) allocate(west_gwstress(ngrid, nlayer)) 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