!

! $Id: acama_gwd_rando_m.F90 3977 2021-08-25 17:24:20Z fhourdin $

!

module ACAMA_GWD_rando_m

  implicit none

contains

    ! Parametrization of the momentum flux deposition due to a discrete

    ! number of gravity waves.

    ! Author: F. Lott, A. de la Camara

    ! July, 24th, 2014

    ! Gaussian distribution of the source, source is vorticity squared

    ! Reference: de la Camara and Lott (GRL, 2015, vol 42, 2071-2078 )

    ! Lott et al (JAS, 2010, vol 67, page 157-170)

    ! Lott et al (JAS, 2012, vol 69, page 2134-2151)

!  ONLINE:

    use dimphy, only: klon, klev

    use assert_m, only: assert

    USE ioipsl_getin_p_mod, ONLY : getin_p

    USE vertical_layers_mod, ONLY : presnivs

    include "YOMCST.h"

    include "clesphys.h"

!  OFFLINE:

!   include "dimensions.h"

!   include "dimphy.h"

!END DIFFERENCE

    include "YOEGWD.h"

    ! 0. DECLARATIONS:

    ! 0.1 INPUTS

    REAL, intent(in)::DTIME ! Time step of the Physics

    REAL, intent(in):: PP(:, :) ! (KLON, KLEV) Pressure at full levels

    REAL, intent(in):: ROT(:,:) ! Relative vorticity

    REAL, intent(in):: TT(:, :) ! (KLON, KLEV) Temp at full levels

    REAL, intent(in):: UU(:, :) ! (KLON, KLEV) Zonal wind at full levels

    REAL, intent(in):: VV(:, :) ! (KLON, KLEV) Merid wind at full levels

    REAL, intent(in):: PLAT(:) ! (KLON) LATITUDE

    ! 0.2 OUTPUTS

    REAL, intent(out):: zustr(:), zvstr(:) ! (KLON) Surface Stresses

    REAL, intent(inout):: d_u(:, :), d_v(:, :)

    REAL, intent(inout):: east_gwstress(:, :) !  Profile of eastward stress

    REAL, intent(inout):: west_gwstress(:, :) !  Profile of westward stress

    ! (KLON, KLEV) tendencies on winds

    ! O.3 INTERNAL ARRAYS

    INTEGER II, JJ, LL

    ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED

    REAL DELTAT

    ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS

    INTEGER, PARAMETER:: NK = 2, NP = 2, NO = 2, NW = NK * NP * NO

    INTEGER JK, JP, JO, JW

    INTEGER, PARAMETER:: NA = 5  !number of realizations to get the phase speed

    REAL KMIN, KMAX ! Min and Max horizontal wavenumbers

    REAL CMIN, CMAX ! Min and Max absolute ph. vel.

    REAL CPHA ! absolute PHASE VELOCITY frequency

    ! Waves Intr. freq. at the 1/2 lev surrounding the full level

    ! Wave EP-fluxes at the 2 semi levels surrounding the full level

    ! Fluxes X and Y for each waves at 1/2 Levels

    INTEGER LAUNCH, LTROP ! Launching altitude and tropo altitude

    REAL XLAUNCH ! Controle the launching altitude

    REAL XTROP ! SORT of Tropopause altitude

    REAL PRMAX ! Maximum value of PREC, and for which our linear formula

    ! for GWs parameterisation apply

    ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS

    REAL RDISS, ZOISEC ! COEFF DE DISSIPATION, SECURITY FOR INTRINSIC FREQ

    REAL CORSEC ! SECURITY FOR INTRINSIC CORIOLIS

    REAL RUWFRT,SATFRT

    ! 0.3.3 BACKGROUND FLOW AT 1/2 LEVELS AND VERTICAL COORDINATE

    REAL H0 ! Characteristic Height of the atmosphere

    REAL DZ ! Characteristic depth of the source!

    REAL PR, TR ! Reference Pressure and Temperature

    REAL PSEC ! Security to avoid division by 0 pressure

    REAL BVSEC ! Security to avoid negative BVF

    LOGICAL, SAVE :: gwd_reproductibilite_mpiomp=.true.

    LOGICAL, SAVE :: firstcall = .TRUE.

  !$OMP THREADPRIVATE(firstcall,gwd_reproductibilite_mpiomp)

    CHARACTER (LEN=20) :: modname='acama_gwd_rando_m'

    CHARACTER (LEN=80) :: abort_message

    ! Cle introduite pour resoudre un probleme de non reproductibilite

    ! Le but est de pouvoir tester de revenir a la version precedenete

    ! A eliminer rapidement

    ENDIF

!    CALL iophys_ini(dtime)

  ENDIF

    !-----------------------------------------------------------------

    ! 1. INITIALISATIONS

    ! 1.1 Basic parameter

    ! Are provided from elsewhere (latent heat of vaporization, dry

    ! gaz constant for air, gravity constant, heat capacity of dry air

    ! at constant pressure, earth rotation rate, pi).

    ! 1.2 Tuning parameters of V14

! Values for linear in rot (recommended):

!   RUWFRT=0.005 ! As RUWMAX but for frontal waves

!   SATFRT=1.00  ! As SAT    but for frontal waves

! Values when rot^2 is used

!    RUWFRT=0.02  ! As RUWMAX but for frontal waves

!    SATFRT=1.00  ! As SAT    but for frontal waves

!    CMAX = 30.   ! Characteristic phase speed

! Values when rot^2*EXP(-pi*sqrt(J)) is used

!   RUWFRT=2.5   ! As RUWMAX but for frontal waves ~ N0*F0/4*DZ

!   SATFRT=0.60   ! As SAT    but for frontal waves

    CMAX = 50.    ! Characteristic phase speed

! Phase speed test

!   RUWFRT=0.01

!   CMAX = 50.   ! Characteristic phase speed (TEST)

! Values when rot^2 and exp(-m^2*dz^2) are used

!   RUWFRT=0.03  ! As RUWMAX but for frontal waves

!   SATFRT=1.00  ! As SAT    but for frontal waves

! CRUCIAL PARAMETERS FOR THE WIND FILTERING

    XLAUNCH=0.95 ! Parameter that control launching altitude

    RDISS = 0.5  ! Diffusion parameter

    ! maximum of rain for which our theory applies (in kg/m^2/s)

    DZ = 1000. ! Characteristic depth of the source

    XTROP=0.2 ! Parameter that control tropopause altitude

    DELTAT=24.*3600. ! Time scale of the waves (first introduced in 9b)

!   DELTAT=DTIME     ! No AR-1 Accumulation, OR OFFLINE

    KMIN = 2.E-5

    ! minimum horizontal wavenumber (inverse of the subgrid scale resolution)

    KMAX = 1.E-3  ! Max horizontal wavenumber

    CMIN = 1.     ! Min phase velocity

    TR = 240. ! Reference Temperature

    PR = 101300. ! Reference pressure

    BVSEC = 5.E-3 ! Security to avoid negative BVF

    PSEC = 1.E-6 ! Security to avoid division by 0 pressure

    ZOISEC = 1.E-6 ! Security FOR 0 INTRINSIC FREQ

!  ONLINE

    call assert(klon == (/size(pp, 1), size(tt, 1), size(uu, 1), &

         size(vv, 1), size(rot,1), size(zustr), size(zvstr), size(d_u, 1), &

         size(d_v, 1), &

        size(east_gwstress,1), size(west_gwstress,1) /), &

    call assert(klev == (/size(pp, 2), size(tt, 2), size(uu, 2), &

         size(vv, 2), size(d_u, 2), size(d_v, 2), &

         size(east_gwstress,2), size(west_gwstress,2) /), &

!  END ONLINE

!       PRINT *, 'flott_gwd_rando: deltat < dtime!'

!       STOP 1

    ENDIF

!       PRINT *, 'flott_gwd_rando: you will have problem with random numbers'

!       STOP 1

    ENDIF

    ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS

    ! Pressure and Inv of pressure

    end DO

    ! Launching altitude

       ! Reprend la formule qui calcule PH en fonction de PP=play

       end DO

    ELSE

    ENDIF

    LAUNCH=0

    LTROP =0

    ENDDO

    ENDDO

    !LAUNCH=22 ; LTROP=33

!   print*,'LAUNCH=',LAUNCH,'LTROP=',LTROP

!   PRINT *,'LAUNCH IN ACAMARA:',LAUNCH

    ! Log pressure vert. coordinate

    end DO

    ! BV frequency

       ! BVSEC: BV Frequency (UH USED IS AS A TEMPORARY ARRAY DOWN TO WINDS)

       UH(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) &

            * RD**2 / RCPD / H0**2 + (TT(:, LL) &

    end DO

    BVLOW = 0.5 * (TT(:, LTROP )+ TT(:, LAUNCH)) &

         * RD**2 / RCPD / H0**2 + (TT(:, LTROP ) &

    ! SMOOTHING THE BV HELPS

    end DO

    ! WINDS

       UZ(:, LL) = ABS((SQRT(UU(:, LL)**2+VV(:, LL)**2) &

          - SQRT(UU(:,LL-1)**2+VV(:, LL-1)**2)) &

    end DO

   ! BAROTROPIC VORTICITY AND INTEGRATED CORIOLIS PARAMETER

        !ROTBA(:) = ROTBA(:) + (ROT(:,LL)+ROT(:,LL+1))/2./RG*(PP(:,LL)-PP(:,LL+1))

        ! Introducing the complete formula (exp of Richardson number):

        ROTBA(:) = ROTBA(:) + &

                !((ROT(:,LL)+ROT(:,LL+1))/2.)**2 &

                (CORIO(:)*TANH(ABS(ROT(:,LL)+ROT(:,LL+1))/2./CORIO(:)))**2 &

                /RG*(PP(:,LL)-PP(:,LL+1)) &

                * EXP(-RPI*BV(:,LL+1)/UZ(:,LL+1)) &

!               * DZ*BV(:,LL+1)/4./ABS(CORIO(:))

!ARRET

    ENDDO

    !   PRINT *,'MAX ROTBA:',MAXVAL(ROTBA)

    !   ROTBA(:)=(1.*ROTBA(:)  & ! Testing zone

    !           +0.15*CORIO(:)**2                &

    !           /(COS(PLAT(:)*RPI/180.)+0.02)  &

    !           )*DZ*0.01/0.0001/4. ! & ! Testing zone

    !   MODIF GWD4 AFTER 1985

    !            *(1.25+SIN(PLAT(:)*RPI/180.))/(1.05+SIN(PLAT(:)*RPI/180.))/1.25

    !          *1./(COS(PLAT(:)*RPI/180.)+0.02)

    !    CORIO(:) = MAX(ROMEGA*2.*ABS(SIN(PLAT(:)*RPI/180.)),ZOISEC)/RG*PP(:,1)

    ! 3 WAVES CHARACTERISTICS CHOSEN RANDOMLY AT THE LAUNCH ALTITUDE

    ! The mod functions of weird arguments are used to produce the

    ! waves characteristics in an almost stochastic way

    JW = 0

             ! Angle

                ! Angle (0 or PI so far)

                ! ZP(JW, II) = (SIGN(1., 0.5 - MOD(TT(II, JW) * 10., 1.)) + 1.) &

                !      * RPI / 2.

                ! Angle between 0 and pi

! TEST WITH POSITIVE WAVES ONLY (Part I/II)

!               ZP(JW, II) = 0.

                ! Horizontal wavenumber amplitude

                ! Horizontal phase speed

                CPHA = 0.

                    CPHA = CPHA + &

                END DO

! TEST WITH POSITIVE WAVES ONLY (Part II/II)

!               ZP(JW, II) = 0.

                ENDIF

                ! Absolute frequency is imposed

                ! Intrinsic frequency is imposed

                ZO(JW, II) = ZO(JW, II) &

                     + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) &

                ! Momentum flux at launch lev

                ! LAUNCHED RANDOM WAVES WITH LOG-NORMAL AMPLITUDE

                !  RIGHT IN THE SH (GWD4 after 1990)

                    RUW0(JW, II) = RUW0(JW,II) + &

                END DO

                RUW0(JW, II) = RUWFRT &

                          * EXP(RUW0(JW,II))/1250. &  ! 2 mpa at south pole

                ! RUW0(JW, II) = RUWFRT

             ENDDO

    end DO

    ! 4. COMPUTE THE FLUXES

    ! 4.0

    ! 4.1 Vertical velocity at launching altitude to ensure

    ! the correct value to the imposed fluxes.

       ! Evaluate intrinsic frequency at launching altitude:

       ZOP(JW, :) = ZO(JW, :) &

            - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) &

       ! VERSION WITH FRONTAL SOURCES

       ! Momentum flux at launch level imposed by vorticity sources

       ! tanh limitation for values above CORIO (inertial instability).

       ! WWP(JW, :) = RUW0(JW, :) &

       WWP(JW, :) = RUWFRT      &

       !     * (CORIO(:)*TANH(ROTBA(:)/CORIO(:)))**2 &

       !    * ABS((CORIO(:)*TANH(ROTBA(:)/CORIO(:)))*CORIO(:)) &

       !  CONSTANT FLUX

       !    * (CORIO(:)*CORIO(:)) &

       ! MODERATION BY THE DEPTH OF THE SOURCE (DZ HERE)

       !      *EXP(-BVLOW(:)**2/MAX(ABS(ZOP(JW, :)),ZOISEC)**2 &

       !      *ZK(JW, :)**2*DZ**2) &

       ! COMPLETE FORMULA:

            !* CORIO(:)**2*TANH(ROTBA(:)/CORIO(:)**2) &

            * ROTBA(:) &

       !  RESTORE DIMENSION OF A FLUX

       !     *RD*TR/PR

       !     *1. + RUW0(JW, :)

       ! Factor related to the characteristics of the waves: NONE

       ! Moderation by the depth of the source (dz here): NONE

       ! Put the stress in the right direction:

    end DO

    ! 4.2 Uniform values below the launching altitude

       end DO

    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, and vi) testing the breaking.

       ! Warning: all the physics is here (passage from one level

       ! to the next)

          ! Intrinsic Frequency

          ZOP(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LL + 1) &

          ! No breaking (Eq.6)

          ! Dissipation (Eq. 8)

          WWP(JW, :) = WWM(JW, :) * EXP(- 4. * RDISS * PR / (PH(:, LL + 1) &

               + PH(:, LL)) * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 &

               / MAX(ABS(ZOP(JW, :) + ZOM(JW, :)) / 2., ZOISEC)**4 &

          ! Critical levels (forced to zero if intrinsic frequency changes sign)

          ! Saturation (Eq. 12)

          WWP(JW, :) = min(WWP(JW, :), MAX(0., &

               SIGN(1., ZOP(JW, :) * ZOM(JW, :))) * ABS(ZOP(JW, :))**3 &

          !    / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * SATFRT**2 * KMIN**2 &

               / BV(:, LL + 1) * EXP(- ZH(:, LL + 1) / H0) * KMIN**2 &

!              *(SATFRT*(2.5+1.5*TANH((ZH(:,LL+1)/H0-8.)/2.)))**2 &

               *SATFRT**2       &

       end DO

       ! Evaluate EP-flux from Eq. 7 and give the right orientation to

       ! the stress

       end DO

       end DO

    end DO

    ! 5 CALCUL DES TENDANCES:

    ! 5.1 Rectification des flux au sommet et dans les basses couches

    end DO

    ! AR-1 RECURSIVE FORMULA (13) IN VERSION 4

       D_U(:, LL) = (1.-DTIME/DELTAT) * D_U(:, LL) + DTIME/DELTAT/REAL(NW) * &

            RG * (RUW(:, LL + 1) - RUW(:, LL)) &

!  NO AR1 FOR MERIDIONAL TENDENCIES

!      D_V(:, LL) = (1.-DTIME/DELTAT) * D_V(:, LL) + DTIME/DELTAT/REAL(NW) * &

       D_V(:, LL) =                                            1./REAL(NW) * &

            RG * (RVW(:, LL + 1) - RVW(:, LL)) &

    ENDDO

    ! Cosmetic: evaluation of the cumulated stress

!      ZVSTR = ZVSTR + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL))/DTIME

    ENDDO

! COSMETICS TO VISUALIZE ROTBA

end module ACAMA_GWD_rando_m