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, pp, & pn2, presnivs, pt, pu, pv, & 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 ! READAPTED FOR VENUS S.LEBONNOIS 08/2021 !--------------------------------------------------------------------------------- USE ioipsl_getin_p_mod, ONLY : getin_p USE geometry_mod, ONLY: cell_area, latitude_deg !#ifdef CPP_XIOS ! use xios_output_mod, only: send_xios_field !#endif implicit none #include "YOMCST.h" ! 0. DECLARATIONS: ! 0.1 INPUTS INTEGER, intent(in):: ngrid, nlayer REAL, intent(in):: dtime ! Time step of the Physics REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels REAL, intent(in):: pn2(ngrid, nlayer) ! N2 (BV^2) at 1/2 levels REAL, intent(in):: presnivs(nlayer) ! Approximate pressure for reproductible launching altitude REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels ! 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 REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind ! 0.3 INTERNAL ARRAYS REAL :: tt(ngrid, nlayer) ! Temperature at full levels REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels 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 !!! TEST Kobs and BestFit Diogo REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber (lambda min 50 km) ! generic : kmax=N/u, u(=30) zonal wind velocity at launch REAL, parameter:: kmin = 1.e-5 ! Min horizontal wavenumber (lambda max 628 km) ! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid !---------------------------------------- ! VCD 1.1 tuning ! REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity !---------------------------------------- ! VCD 2.0 tuning REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity !---------------------------------------- ! generic : cmax=zonal wind at launch REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity !---------------------------------------- ! Nominal as Gilli2017 ! REAL, parameter:: epflux_max = 0.005 ! Max EP flux value at launching altitude (previous name: RUWMAX) !---------------------------------------- ! VCD 2.1 tuning REAL, parameter:: epflux_max = 0.001 ! Max EP flux value at launching altitude (previous name: RUWMAX) !---------------------------------------- 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) INTEGER launch ! Launching altitude REAL, parameter:: xlaunch = 5e-3 ! Value for top of cloud convective region ! 0.3.2 Parameters of waves dissipations !---------------------------------------- ! VCD 1.1 tuning ! REAL, parameter:: sat = 0.85 ! saturation parameter ! REAL, parameter:: rdiss = 0.1 ! coefficient of dissipation !---------------------------------------- ! VCD 2.0 tuning REAL, parameter:: sat = 0.6 ! saturation parameter REAL, parameter:: rdiss = 8.e-4 ! coefficient of dissipation !---------------------------------------- REAL, parameter:: zoisec = 1.e-8 ! 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 min_k(ngrid) ! min(kstar,kmin) REAL, save:: H0 ! characteristic height of the atmosphere REAL, save:: MDR ! characteristic mass density !---------------------------------------- ! VCD 1.1 tuning ! REAL, parameter:: pr = 5e5 ! Reference pressure [Pa] ! VENUS: cloud layer !---------------------------------------- ! VCD 2.0 tuning REAL, parameter:: pr = 5e4 ! Reference pressure [Pa] ! VENUS: cloud layer !---------------------------------------- REAL, parameter:: tr = 300. ! Reference temperature [K] ! VENUS: cloud layer 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-9 ! Security to avoid division by 0 pressure REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV REAL href(nlayer + 1) ! Reference altitude for launching altitude REAL ran_num_1, ran_num_2, ran_num_3 LOGICAL, save :: firstcall = .true. !-------------------------------------------------------------------------------------------------- ! 1. INITIALISATION !------------------ IF (firstcall) THEN write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" ! Characteristic vertical scale height H0 = RD * tr / RG ! Characteristic mass density at launch altitude MDR = pr / ( RD * tr ) ! 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 ! For VENUS, we use *_seri variables, that already integrate the different previous tendencies tt(:,:) = pt(:,:) uu(:,:) = pu(:,:) vv(:,:) = pv(:,:) ! 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 ! 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 ! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS bv(:, ll) = max(bvsec,sqrt(pn2(:,ll))) 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) ! TN+GG April/June 2020 ! "Individual waves are not supposed ! to occupy the entire domain, but ! only a faction of it" Lott+2012 ! minimum value of k between kmin and kstar DO ii = 1, ngrid kstar = RPI / sqrt(cell_area(ii)) min_k(ii) = max(kmin,kstar) end DO ! 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 DO ii = 1, ngrid ran_num_1 = mod(tt(ii, jw) * 10., 1.) ran_num_2 = mod(tt(ii, jw) * 100., 1.) !! OPTIONS GENERIC DIFF VENUS !! ! angle (random) - reference zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & + 1.) * RPI / 2. ! zp(jw, ii) = atan2(4.*vh(ii, launch),uh(ii, launch)) & ! + (sign(1., 0.5 - ran_num_1) + 1.) * RPI / 2. ! --------- TRY 00 ----------------------- !IF((abs(latitude_deg(ii)) .le. 55.)) THEN ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & ! * (10.) ! +/- 10 deg par rapport equateur !ELSE IF ((abs(latitude_deg(ii)) .le. 75.)) THEN ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & ! * (10. +(90.-10.) & ! * (latitude_deg(ii)-55.)/20.) !ELSE ! zp(jw, ii) = ran_num_1 * RPI !ENDIF !zp(jw, ii) = mod(zp(jw, ii),2.*RPI) ! ------ TRY 01------------------------------- !IF((abs(latitude_deg(ii)) .le. 55)) THEN ! zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & ! + 1.) * RPI / 2. !ELSE ! zp(jw, ii) = ran_num_1 * RPI !ENDIF ! ---- angle (0 or RPI) ----- !!zp(jw, ii) = RPI*mod(jw,2) ! horizontal wavenumber amplitude ! zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 zk(jw, ii) = min_k(ii) + (kmax - min_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. ! VCD 2.0 tuning 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 !cpha = cpha - uh(ii, launch) IF (cpha < 0.) THEN cpha = - 1. * cpha zp (jw, ii) = zp(jw, ii) + RPI ENDIF !IF (cpha < 1.) THEN ! cpha = 1. !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.) end DO 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 altutide: intr_freq_p(jw, :) = zo(jw, :) & - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) & - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) end DO ! 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 ! 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): * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) & * ((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) * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & * exp(-zh(:, ll + 1) / H0) & !! Correction for VCD 2.0 ! * sat**2 * kmin**2 / zk(jw, :)**4) * sat**2 * min_k(:)**2 * MDR / zk(jw, :)**4) end DO ! 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, ll)) / float(nw) end DO end DO ! DO LL = LAUNCH, nlayer - 1 ! print*,'u_epflux_tot just after launching:', u_epflux_tot ! print*,'v_epflux_tot just after launching:', v_epflux_tot ! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p) ! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p) ! 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 ! d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) d_u(:, LL) = RG * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) & / (PH(:, LL + 1) - PH(:, LL)) * DTIME ! d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) d_v(:, LL) = RG * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) & / (PH(:, LL + 1) - PH(:, LL)) * DTIME 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(:,:) ! print*,'u_epflux_tot just after tendency:', u_epflux_tot ! print*,'v_epflux_tot just after tendency:', v_epflux_tot ! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u) ! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v) ! Cosmetic: evaluation of the cumulated stress ZUSTR(:) = 0. ZVSTR(:) = 0. DO LL = 1, nlayer ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL)) ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (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