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

CONTAINS

      SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, 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
!---------------------------------------------------------------------------------



     use comcstfi_mod, only: g, pi, cpp, r
     USE ioipsl_getin_p_mod, ONLY : getin_p
     use assert_m, only : assert
     use vertical_layers_mod, only : presnivs
     use callkeys_mod, only : epflux_max, & ! Max EP flux value at launching altitude (previous name: RUWMAX)
                              gwd_convective_source
     use inifis_mod
#ifdef CPP_XIOS
     use xios_output_mod, only: send_xios_field
#endif
 implicit none

! 0. DECLARATIONS:

     ! 0.1 INPUTS

     INTEGER, intent(in):: ngrid, nlayer
     REAL, intent(in):: dtime              ! Time step of the Physics
     REAL, intent(in):: zmax_therm(ngrid)  ! Altitude of max velocity thermals (m)
     REAL, intent(in):: pp(ngrid, nlayer)  ! Pressure at full levels
     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
     REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature
     REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind
     REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind
     
     ! 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
     REAL :: bvlow(ngrid)      ! Qu est-ce ? 
     REAL :: dz

     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, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber
     REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber
     !REAL, parameter:: cmax = 30.   ! Max horizontal absolute phase velocity
     REAL, parameter:: cmax = 100.  ! Test for Saturn: Max horizontal absolute phase velocity
     REAL, parameter:: cmin = 1.    ! Min horizontal absolute phase velocity
     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 = 0.4      ! Control the launching altitude
     REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.)

     REAL prec(ngrid)     ! precipitations
     REAL prmax           ! Max value of precipitation

     ! 0.3.2 Parameters of waves dissipations
     REAL, parameter:: sat   = 1.     ! saturation parameter
     REAL, parameter:: rdiss = 1.     ! coefficient of dissipation
     REAL, parameter:: zoisec = 1.e-6 ! 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
     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 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-6  ! 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 = 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(:,:)
!    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 

        bv(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1)) &
                    * r**2 / cpp / H0**2              &
                    + (tt(:, ll) - tt(:, ll - 1))     &
                    / (zh(:, ll) - zh(:, ll - 1))     &
                    * r / H0 
        bv(:, ll) = sqrt(max(bvsec, bv(:,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)

! 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.)

           ! angle (0 or pi so far)
           zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &
                        + 1.) * pi / 2.
           ! horizontal wavenumber amplitude
           zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2
           ! horizontal phase speed
           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
           ! 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_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):
                      * 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)                                          &
                      * sat**2 * kmin**2 / 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) = G * (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) = G * (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) !/ 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

