MODULE lmdz_gwd_ini
! -----------------------------------------------------
! This module initializes all the parameters needed
! for gravity wave drag (gwd) parameterizations
!------------------------------------------------------

   IMPLICIT NONE

   REAL :: RCPD, RD, RG, RPI, RA, ROMEGA, RKAPPA, RLVTT
   !$OMP THREADPRIVATE(RCPD, RD, RG, RPI, RA, ROMEGA, RKAPPA, RLVTT)

   INTEGER, SAVE :: lunout, prt_level            ! Logical unit number and level for standard output
   !$OMP THREADPRIVATE(lunout,prt_level)

   ! for orographic gravity wave effects
   LOGICAL :: ok_strato
   !$OMP THREADPRIVATE(ok_strato)

   LOGICAL, SAVE, PROTECTED :: ok_orodr = .true.  ! activate orographic gravity wave drag
   !$OMP THREADPRIVATE(ok_orodr)

   LOGICAL, SAVE, PROTECTED :: ok_orolf = .true.  ! activate orographic gravity wave lift
   !$OMP THREADPRIVATE(ok_orolf)

   REAL, SAVE, PROTECTED :: zpmm_orodr_t = 100.  ! threshold on zpeak-zmean [m] to activate sso gwd
   !$OMP THREADPRIVATE(zpmm_orodr_t)

   REAL, SAVE, PROTECTED :: zstd_orodr_t = 10.    ! threshold on zstd [m] to activate sso gwd
   !$OMP THREADPRIVATE(zstd_orodr_t)

   REAL, SAVE, PROTECTED :: zpmm_orolf_t = 100.  ! threshold on zpeak-zmean [m] to activate sso lift
   !$OMP THREADPRIVATE(zpmm_orolf_t)

   REAL, SAVE, PROTECTED :: nm_oro_t = -1      ! number of subgrid-scale mountains above which sso drag and lift are active
   !$OMP THREADPRIVATE(nm_oro_t)

   INTEGER, SAVE, PROTECTED :: NKTOPG ! Security value for blocked flow level
   !$OMP THREADPRIVATE(NKTOPG)

   INTEGER, SAVE, PROTECTED :: NSTRA ! An estimate to qualify the upper levels of the model where one wants to impose stress
   !$OMP THREADPRIVATE(NSTRA)

   REAL, SAVE, PROTECTED :: GFRCRIT = 1. ! Critical Non-dimensional mountain Height (HNC in (1), LOTT 1999)
   !$OMP THREADPRIVATE(GFRCRIT)

   REAL, SAVE, PROTECTED :: GKWAKE = 0.5 ! Bluff-body drag coefficient for low level wake (Cd in (2), LOTT 1999)
   !$OMP THREADPRIVATE(GKWAKE)

   REAL, SAVE, PROTECTED :: GRCRIT ! Critical Richardson Number (Ric, End of first column p791 of LOTT 1999)
   !$OMP THREADPRIVATE(GRCRIT)

   REAL, SAVE, PROTECTED :: GVCRIT
   !$OMP THREADPRIVATE(GVCRIT)

   REAL, SAVE, PROTECTED :: GKDRAG ! Gravity wave drag coefficient (G in (3), LOTT 1999)
   !$OMP THREADPRIVATE(GKDRAG)

   REAL, SAVE, PROTECTED :: GKLIFT ! Mountain Lift coefficient (Cl in (4), LOTT 1999)
   !$OMP THREADPRIVATE(GKLIFT)

   REAL, SAVE, PROTECTED :: GHMAX ! Not used
   !$OMP THREADPRIVATE(GHMAX)

   REAL, SAVE, PROTECTED :: GRAHILO ! Set-up the trapped waves fraction (Beta, End of first column,  LOTT 1999)
   !$OMP THREADPRIVATE(GRAHILO)

   REAL, SAVE, PROTECTED :: GSIGCR ! Security value for blocked flow depth
   !$OMP THREADPRIVATE(GSIGCR)

   REAL, SAVE, PROTECTED :: GSSEC ! Security min value for low-level B-V frequency
   !$OMP THREADPRIVATE(GSSEC)

   REAL, SAVE, PROTECTED :: GTSEC ! Security min value for anisotropy and GW stress.
   !$OMP THREADPRIVATE(GTSEC)

   REAL, SAVE, PROTECTED :: GVSEC ! Security min value for ulow
   !$OMP THREADPRIVATE(GVSEC)

   ! for non-orographic gravity wave drag
   LOGICAL, SAVE, PROTECTED :: ok_gwd_rando = .false.   ! activate non-orographic stochastic gravity wave drag params
   !$OMP THREADPRIVATE(ok_gwd_rando)

   REAL, SAVE, PROTECTED :: GWD_RANDO_RUWMAX = 2.    ! Maximum Eliassen-Palm flux at launch level, in "FLOTT_GWD_rando"
   !$OMP THREADPRIVATE(GWD_RANDO_RUWMAX)

   REAL, SAVE, PROTECTED :: GWD_RANDO_SAT = 0.25     ! saturation parameter in "FLOTT_GWD_rando"  S_c in equation (12) of Lott (JGR, vol 118, page 8897, 2013)
   !$OMP THREADPRIVATE(GWD_RANDO_SAT)

   REAL, SAVE, PROTECTED ::  GWD_FRONT_RUWMAX = 2.5  ! Same as GWD_RANDO params but for fronal GWs
   !$OMP THREADPRIVATE(GWD_FRONT_RUWMAX)

   REAL, SAVE, PROTECTED ::  GWD_FRONT_SAT = 0.60  ! Same as GWD_RANDO params but for fronal GWs
   !$OMP THREADPRIVATE(GWD_FRONT_SAT)

   INTEGER, SAVE, PROTECTED :: NK = 2
   !$OMP THREADPRIVATE(NK)

   INTEGER, SAVE, PROTECTED :: NP = 2
   !$OMP THREADPRIVATE(NP)

   INTEGER, SAVE, PROTECTED :: NO = 2
   !$OMP THREADPRIVATE(NO)

   INTEGER, SAVE, PROTECTED :: NA = 5 ! number of realizations to get the phase speed
   !$OMP THREADPRIVATE(NA)

   INTEGER, SAVE, PROTECTED :: NW
   !$OMP THREADPRIVATE(NW)

   ! for TKE production by SSO drag
   INTEGER, SAVE, PROTECTED :: addtkeoro = 0     ! activate TKE production by sso drag
   !$OMP THREADPRIVATE(addtkeoro)

   LOGICAL, SAVE, PROTECTED :: smallscales_tkeoro = .false.  ! considers all sso scales for TKE production by sso drag if addtkeoro >=2
   !$OMP THREADPRIVATE(smallscales_tkeoro)

   REAL, SAVE, PROTECTED :: alphatkeoro = 1.    ! tuning parameter for TKE production by sso drag
   !$OMP THREADPRIVATE(alphatkeoro)

   ! Others

   LOGICAL, SAVE :: firstcall = .true.
   !$OMP THREADPRIVATE(firstcall)

   LOGICAL, SAVE :: gwd_reproductibilite_mpiomp = .true.
   !$OMP THREADPRIVATE(gwd_reproductibilite_mpiomp)

CONTAINS

!===================================================================================================================================================
   SUBROUTINE gwd_ini(klon, klev, pplay, paprs, lunout_in, prt_level_in, &
                      RCPD_in, RD_in, RG_in, RPI_in, ROMEGA_in, RA_in, &
                      RKAPPA_in, RLVTT_in, ok_strato_in)

      USE ioipsl_getin_p_mod, ONLY: getin_p

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: klon, klev, lunout_in, prt_level_in
      REAL, DIMENSION(klon, klev) :: pplay
      REAL, DIMENSION(klon, klev + 1) :: paprs
      REAL, INTENT(IN) :: RCPD_in, RD_in, RG_in, RPI_in, ROMEGA_in
      REAL, INTENT(IN) :: RLVTT_IN, RA_in, RKAPPA_in
      LOGICAL, INTENT(IN) :: ok_strato_in
      CHARACTER(len=80) :: abort_message

      lunout = lunout_in
      prt_level = prt_level_in
      RCPD = RCPD_in
      RD = RD_in
      RG = RG_in
      RPI = RPI_in
      ROMEGA = ROMEGA_in
      RA = RA_in
      RKAPPA = RKAPPA_in
      RLVTT = RLVTT_in
      ok_strato = ok_strato_in

! default value of orographic gravity wave drag depend on ok_strato
      IF (ok_strato) THEN
         gkdrag = 0.1875
         grahilo = 0.1
         grcrit = 1.
         gklift = 0.25
      ELSE
         gkdrag = 0.2
         grahilo = 1.
         grcrit = 0.01
         gklift = 0.50
      END IF

! getin from .def files

      CALL getin_p('ok_orodr', ok_orodr)
      CALL getin_p('ok_orolf', ok_orolf)
      CALL getin_p('ok_gwd_rando', ok_gwd_rando)
      CALL getin_p('zpmm_orodr_t', zpmm_orodr_t)
      CALL getin_p('zpmm_orolf_t', zpmm_orolf_t)
      CALL getin_p('zstd_orodr_t', zstd_orodr_t)
      CALL getin_p('nm_oro_t', nm_oro_t)
      CALL getin_p('addtkeoro', addtkeoro)
      CALL getin_p('alphatkeoro', alphatkeoro)
      CALL getin_p('smallscales_tkeoro', smallscales_tkeoro)
      CALL getin_p('gwd_rando_ruwmax', gwd_rando_ruwmax)
      CALL getin_p('gwd_rando_sat', gwd_rando_sat)
      CALL getin_p('gwd_front_ruwmax', gwd_front_ruwmax)
      CALL getin_p('gwd_front_sat', gwd_front_sat)
      CALL getin_p('sso_gkdrag', gkdrag)
      CALL getin_p('sso_grahil', grahilo)
      CALL getin_p('sso_grcrit', grcrit)
      CALL getin_p('sso_gfrcri', gfrcrit)
      CALL getin_p('sso_gkwake', gkwake)
      CALL getin_p('sso_gklift', gklift)

! write in used.def files
      WRITE (lunout, *) 'gwd_ini, ok_orodr:', ok_orodr
      WRITE (lunout, *) 'gwd_ini, ok_orolf:', ok_orolf
      WRITE (lunout, *) 'gwd_ini, ok_gwd_rando:', ok_gwd_rando
      WRITE (lunout, *) 'gwd_ini, zpmm_orodr_t:', zpmm_orodr_t
      WRITE (lunout, *) 'gwd_ini, zstd_orodr_t:', zstd_orodr_t
      WRITE (lunout, *) 'gwd_ini, zpmm_orolf_t:', zpmm_orolf_t
      WRITE (lunout, *) 'gwd_ini, nm_oro_t:', nm_oro_t
      WRITE (lunout, *) 'gwd_ini, addtkeoro:', addtkeoro
      WRITE (lunout, *) 'gwd_ini, alphatkeoro:', alphatkeoro
      WRITE (lunout, *) 'gwd_ini, smallscales_tkeoro:', smallscales_tkeoro
      WRITE (lunout, *) 'gwd_ini, gwd_rando_ruwmax:', gwd_rando_ruwmax
      WRITE (lunout, *) 'gwd_ini, gwd_rando_sat:', gwd_rando_sat
      WRITE (lunout, *) 'gwd_ini, gwd_front_ruwmax:', gwd_front_ruwmax
      WRITE (lunout, *) 'gwd_ini, gwd_front_sat:', gwd_front_sat
      WRITE (lunout, *) 'gwd_ini, gklift:', gklift
      WRITE (lunout, *) 'gwd_ini, gkwake:', gkwake
      WRITE (lunout, *) 'gwd_ini, gfrcrit:', gfrcrit
      WRITE (lunout, *) 'gwd_ini, grcrit:', grcrit
      WRITE (lunout, *) 'gwd_ini, grahilo:', grahilo

! few checks:

      IF (klon == 1 .AND. ok_gwd_rando) THEN
         print*, 'stochastic non-orographic gravity wave drag param does not work in 1D'
         print*, 'I set ok_gwd_rando to false'
         ok_gwd_rando=.false.
      END IF

! initialisation of ogwd parameters

      IF (ok_strato) THEN
         CALL sugwd_strato(klon, klev, paprs, pplay)
      ELSE
         CALL sugwd(klon, klev, paprs, pplay)
      END IF

! initialisation of non-orographic gwd parameterisation
      NW = NK*NP*NO
      IF (ok_gwd_rando) THEN
         CALL gwd_rando_first(klev)
      END IF

   END SUBROUTINE gwd_ini

!===================================================================================================================================================
   SUBROUTINE sugwd_strato(nlon, nlev, paprs, pplay)

      ! **** *SUGWD* INITIALIZE PARAMETERS CONTROLLING GRAVITY WAVE DRAG

      ! PURPOSE.
      ! --------
      ! INITIALIZE PARAMETERS THAT CONTROLS THE
      ! GRAVITY WAVE DRAG PARAMETRIZATION.
      ! VERY IMPORTANT:
      ! ______________
      ! THIS ROUTINE SET_UP THE "TUNABLE PARAMETERS" OF THE
      ! VARIOUS SSO SCHEMES

      ! **   INTERFACE.
      ! ----------
      ! CALL *SUGWD* FROM *SUPHEC*
      ! -----        ------

      ! EXPLICIT ARGUMENTS :
      ! --------------------
      ! PAPRS,PPLAY : Pressure at semi and full model levels
      ! NLEV        : number of model levels
      ! NLON        : number of points treated in the physics

      ! IMPLICIT ARGUMENTS :
      ! --------------------
      ! -GFRCRIT-R:  Critical Non-dimensional mountain Height
      ! (HNC in (1),    LOTT 1999)
      ! -GKWAKE--R:  Bluff-body drag coefficient for low level wake
      ! (Cd in (2),     LOTT 1999)
      ! -GRCRIT--R:  Critical Richardson Number
      ! (Ric, End of first column p791 of LOTT 1999)
      ! -GKDRAG--R:  Gravity wave drag coefficient
      ! (G in (3),      LOTT 1999)
      ! -GKLIFT--R:  Mountain Lift coefficient
      ! (Cl in (4),     LOTT 1999)
      ! -GHMAX---R:  Not used
      ! -GRAHILO-R:  Set-up the trapped waves fraction
      ! (Beta , End of first column,  LOTT 1999)

      ! -GSIGCR--R:  Security value for blocked flow depth
      ! -NKTOPG--I:  Security value for blocked flow level
      ! -nstra----I:  An estimate to qualify the upper levels of
      ! the model where one wants to impose strees
      ! profiles
      ! -GSSECC--R:  Security min value for low-level B-V frequency
      ! -GTSEC---R:  Security min value for anisotropy and GW stress.
      ! -GVSEC---R:  Security min value for ulow

      ! METHOD.
      ! -------
      ! SEE DOCUMENTATION

      ! EXTERNALS.
      ! ----------
      ! NONE

      ! REFERENCE.
      ! ----------
      ! Lott, 1999: Alleviation of stationary biases in a GCM through...
      ! Monthly Weather Review, 127, pp 788-801.

      ! AUTHOR.
      ! -------
      ! FRANCOIS LOTT        *LMD*

      ! MODIFICATIONS.
      ! --------------
      ! ORIGINAL : 90-01-01 (MARTIN MILLER, ECMWF)
      ! LAST:  99-07-09     (FRANCOIS LOTT,LMD)
      ! ------------------------------------------------------------------

      USE dimphy
      USE mod_phys_lmdz_para
      USE mod_grid_phy_lmdz
      USE geometry_mod

      IMPLICIT NONE

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

      ! ARGUMENTS
      INTEGER nlon, nlev
      REAL paprs(nlon, nlev + 1)
      REAL pplay(nlon, nlev)

      INTEGER jk
      REAL zpr, ztop, zsigt, zpm1r
      INTEGER :: cell, ij, nstra_tmp, nktopg_tmp
      REAL :: current_dist, dist_min, dist_min_glo

      ! *       1.    SET THE VALUES OF THE PARAMETERS
      ! --------------------------------

      ghmax = 10000.
      zpr = 100000.
      ZTOP = 0.00005
      zsigt = 0.94

!ym Take the point at equator close to (0,0) coordinates.
      dist_min = 360
      dist_min_glo = 360.
      cell = -1
      DO ij = 1, klon
         current_dist = sqrt(longitude_deg(ij)**2 + latitude_deg(ij)**2)
         current_dist = current_dist*(1 + (1e-10*ind_cell_glo(ij))/klon_glo) ! For point unicity
         IF (dist_min > current_dist) THEN
            dist_min = current_dist
            cell = ij
         END IF
      END DO

      CALL reduce_min(dist_min, dist_min_glo)
      CALL bcast(dist_min_glo)
      IF (dist_min /= dist_min_glo) cell = -1
!ym in future find the point at equator close to (0,0) coordinates.
      PRINT *, 'SUGWD distmin dist_min_glo cell=', dist_min, dist_min_glo, cell

      nktopg_tmp = nktopg
      nstra_tmp = nstra

      IF (cell /= -1) THEN

         DO jk = 1, nlev
            !zpm1r = pplay(cell+1, jk)/paprs(cell+1, 1)
            zpm1r = pplay(cell, jk)/paprs(cell, 1)
            IF (zpm1r >= zsigt) THEN
               nktopg_tmp = jk
            END IF
            IF (zpm1r >= ztop) THEN
               nstra_tmp = jk
            END IF
         END DO
      ELSE
         nktopg_tmp = 0
         nstra_tmp = 0
      END IF

      CALL reduce_sum(nktopg_tmp, nktopg)
      CALL bcast(nktopg)
      CALL reduce_sum(nstra_tmp, nstra)
      CALL bcast(nstra)

      ! inversion car dans orodrag on compte les niveaux a l'envers
      nktopg = nlev - nktopg + 1
      nstra = nlev - nstra
      PRINT *, ' DANS SUGWD nktopg=', nktopg
      PRINT *, ' DANS SUGWD nstra=', nstra
      if (nstra == 0) call abort_physic("sugwd_strato", "no level in stratosphere", 1)

      gsigcr = 0.80 ! Top of low level flow
      gvcrit = 0.1

      WRITE (UNIT=6, FMT='('' *** SSO essential constants ***'')')
      WRITE (UNIT=6, FMT='('' *** SPECIFIED IN SUGWD ***'')')
      WRITE (UNIT=6, FMT='('' Gravity wave ct '',E13.7,'' '')') gkdrag
      WRITE (UNIT=6, FMT='('' Trapped/total wave dag '',E13.7,'' '')') grahilo
      WRITE (UNIT=6, FMT='('' Critical Richardson   = '',E13.7,'' '')') grcrit
      WRITE (UNIT=6, FMT='('' Critical Froude'',e13.7)') gfrcrit
      WRITE (UNIT=6, FMT='('' Low level Wake bluff cte'',e13.7)') gkwake
      WRITE (UNIT=6, FMT='('' Low level lift  cte'',e13.7)') gklift

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

      ! *       2.    SET VALUES OF SECURITY PARAMETERS
      ! ---------------------------------

      gvsec = 0.10
      gssec = 0.0001
      gtsec = 0.00001

      RETURN
   END SUBROUTINE sugwd_strato

!===================================================================================================================================================
   SUBROUTINE sugwd(nlon, nlev, paprs, pplay)

      USE dimphy
      USE mod_phys_lmdz_para
      USE mod_grid_phy_lmdz
      ! USE parallel

      ! **** *SUGWD* INITIALIZE PARAMETERS CONTROLLING GRAVITY WAVE DRAG

      ! PURPOSE.
      ! --------
      ! INITIALIZE PARAMETERS THAT CONTROLS THE
      ! GRAVITY WAVE DRAG PARAMETRIZATION.

      ! **   INTERFACE.
      ! ----------
      ! CALL *SUGWD* FROM *SUPHEC*
      ! -----        ------

      ! EXPLICIT ARGUMENTS :
      ! --------------------
      ! PSIG        : VERTICAL COORDINATE TABLE
      ! NLEV        : NUMBER OF MODEL LEVELS

      ! IMPLICIT ARGUMENTS :
      ! --------------------
      ! METHOD.
      ! -------
      ! SEE DOCUMENTATION

      ! EXTERNALS.
      ! ----------
      ! NONE

      ! REFERENCE.
      ! ----------
      ! ECMWF Research Department documentation of the IFS

      ! AUTHOR.
      ! -------
      ! MARTIN MILLER             *ECMWF*

      ! MODIFICATIONS.
      ! --------------
      ! ORIGINAL : 90-01-01
      ! ------------------------------------------------------------------
      IMPLICIT NONE

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

      INTEGER nlon, nlev, jk
      REAL paprs(nlon, nlev + 1)
      REAL pplay(nlon, nlev)
      REAL zpr, zstra, zsigt, zpm1r
      REAL :: pplay_glo(klon_glo, nlev)
      REAL :: paprs_glo(klon_glo, nlev + 1)

      ! *       1.    SET THE VALUES OF THE PARAMETERS
      ! --------------------------------

      ghmax = 10000.
      zpr = 100000.
      zstra = 0.1
      zsigt = 0.94

      CALL gather(pplay, pplay_glo)
      CALL bcast(pplay_glo)
      CALL gather(paprs, paprs_glo)
      CALL bcast(paprs_glo)

      DO jk = 1, nlev
         zpm1r = pplay_glo((klon_glo/2) + 1, jk)/paprs_glo((klon_glo/2) + 1, 1)
         IF (zpm1r >= zsigt) THEN
            nktopg = jk
         END IF
         zpm1r = pplay_glo((klon_glo/2) + 1, jk)/paprs_glo((klon_glo/2) + 1, 1)
         IF (zpm1r >= zstra) THEN
            nstra = jk
         END IF
      END DO

      ! inversion car dans orodrag on compte les niveaux a l'envers
      nktopg = nlev - nktopg + 1
      nstra = nlev - nstra
      PRINT *, ' DANS SUGWD nktopg=', nktopg
      PRINT *, ' DANS SUGWD nstra=', nstra

      gsigcr = 0.80
      gvcrit = 0.0

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

      ! *       2.    SET VALUES OF SECURITY PARAMETERS
      ! ---------------------------------

      gvsec = 0.10
      gssec = 1.E-12

      gtsec = 1.E-07

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

      RETURN
   END SUBROUTINE sugwd
!===================================================================================================================================================

   SUBROUTINE gwd_rando_first(klev)

      USE ioipsl_getin_p_mod, ONLY: getin_p

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: klev
      CHARACTER(LEN=20), PARAMETER :: modname = 'gwd_rando_first'
      CHARACTER(LEN=80) :: abort_message

      IF (firstcall) THEN
         ! Key to solve a non-reproductability issue. The aim is to test to get back to previous version
         ! to remove asap
         CALL getin_p('gwd_reproductibilite_mpiomp', gwd_reproductibilite_mpiomp)

         IF (NW + 3*NA >= KLEV) THEN
            abort_message = 'NW+3*NA>=KLEV Problem to generate waves associated with precip.'
            CALL abort_physic(modname, abort_message, 1)
         END IF

         IF (NW + 4*(NA - 1) + NA >= KLEV) THEN
             abort_message = 'NW+3*NA>=KLEV Problem to generate waves associated with fronts'
             CALL abort_physic(modname, abort_message, 1)
         END IF

         firstcall = .false.
      END IF
   END SUBROUTINE gwd_rando_first

!===================================================================================================================================================
END MODULE lmdz_gwd_ini
