MODULE lmdz_wake

! $Id: lmdz_wake.F90 4588 2023-06-28 22:54:46Z fhourdin $

CONTAINS

                tenv0, qe0, omgb, &

                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &

                sigd_con, Cin, &

                dth, hw, wape, fip, gfl, &

                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &

                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &

                d_deltat_gw, &                                                      ! tendencies

                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)             ! tendencies

  ! **************************************************************

  ! *

  ! WAKE                                                        *

  ! retour a un Pupper fixe                                *

  ! *

  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *

  ! modified by :   ROEHRIG Romain        01/29/2007            *

  ! **************************************************************

  USE lmdz_wake_ini , ONLY : wake_ini

  USE lmdz_wake_ini , ONLY : prt_level,epsim1,RG,RD

  USE lmdz_wake_ini , ONLY : stark, wdens_ref, coefgw, alpk, wk_pupper

  USE lmdz_wake_ini , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl

  USE lmdz_wake_ini , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin

  USE lmdz_wake_ini , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc

  USE lmdz_wake_ini , ONLY : iflag_wk_profile

  IMPLICIT NONE

  ! ============================================================================

  ! But : Decrire le comportement des poches froides apparaissant dans les

  ! grands systemes convectifs, et fournir l'energie disponible pour

  ! le declenchement de nouvelles colonnes convectives.

  ! State variables :

  ! deltatw    : temperature difference between wake and off-wake regions

  ! deltaqw    : specific humidity difference between wake and off-wake regions

  ! sigmaw     : fractional area covered by wakes.

  ! wdens      : number of wakes per unit area

  ! Variable de sortie :

  ! wape : WAke Potential Energy

  ! fip  : Front Incident Power (W/m2) - ALP

  ! gfl  : Gust Front Length per unit area (m-1)

  ! dtls : large scale temperature tendency due to wake

  ! dqls : large scale humidity tendency due to wake

  ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)

  ! dp_omgb : vertical gradient of large scale omega

  ! awdens  : densite de poches actives

  ! wdens   : densite de poches

  ! omgbdth: flux of Delta_Theta transported by LS omega

  ! dtKE   : differential heating (wake - unpertubed)

  ! dqKE   : differential moistening (wake - unpertubed)

  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)

  ! dp_deltomg  : vertical gradient of omg (s-1)

  ! wkspread  : spreading term in d_t_wake and d_q_wake

  ! deltatw     : updated temperature difference (T_w-T_u).

  ! deltaqw     : updated humidity difference (q_w-q_u).

  ! sigmaw      : updated wake fractional area.

  ! d_deltat_gw : delta T tendency due to GW

  ! Variables d'entree :

  ! aire : aire de la maille

  ! tenv0  : temperature dans l'environnement  (K)

  ! qe0  : humidite dans l'environnement     (kg/kg)

  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)

  ! dtdwn: source de chaleur due aux descentes (K/s)

  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)

  ! dta  : source de chaleur due courants satures et detrain  (K/s)

  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)

  ! wgen : number of wakes generated per unit area and per sec (/m^2/s)

  ! amdwn: flux de masse total des descentes, par unite de

  !        surface de la maille (kg/m2/s)

  ! amup : flux de masse total des ascendances, par unite de

  !        surface de la maille (kg/m2/s)

  ! sigd_con:

  ! Cin  : convective inhibition

  ! p    : pressions aux milieux des couches (Pa)

  ! ph   : pressions aux interfaces (Pa)

  ! pi  : (p/p_0)**kapa (adim)

  ! dtime: increment temporel (s)

  ! Variables internes :

  ! rhow : masse volumique de la poche froide

  ! rho  : environment density at P levels

  ! rhoh : environment density at Ph levels

  ! tenv   : environment temperature | may change within

  ! qe   : environment humidity    | sub-time-stepping

  ! the  : environment potential temperature

  ! thu  : potential temperature in undisturbed area

  ! tu   :  temperature  in undisturbed area

  ! qu   : humidity in undisturbed area

  ! dp_omgb: vertical gradient og LS omega

  ! omgbw  : wake average vertical omega

  ! dp_omgbw: vertical gradient of omgbw

  ! omgbdq : flux of Delta_q transported by LS omega

  ! dth  : potential temperature diff. wake-undist.

  ! th1  : first pot. temp. for vertical advection (=thu)

  ! th2  : second pot. temp. for vertical advection (=thw)

  ! q1   : first humidity for vertical advection

  ! q2   : second humidity for vertical advection

  ! d_deltatw   : terme de redistribution pour deltatw

  ! d_deltaqw   : terme de redistribution pour deltaqw

  ! deltatw0   : deltatw initial

  ! deltaqw0   : deltaqw initial

  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)

  ! amflux : horizontal mass flux through wake boundary

  ! wdens_ref: initial number of wakes per unit area (3D) or per

  ! unit length (2D), at the beginning of each time step

  ! Tgw    : 1 sur la periode de onde de gravite

  ! Cgw    : vitesse de propagation de onde de gravite

  ! LL     : distance entre 2 poches

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

  ! Declaration de variables

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

  ! Arguments en entree

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

  INTEGER, INTENT(IN) :: klon,klev

  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi

  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb

  REAL,                             INTENT(IN)          :: dtime

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tenv0, qe0

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup

  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa

  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen

  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con

  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin

  !

  ! Input/Output

  ! State variables

  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw

  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw

  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens

  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens

  ! Sorties

  ! --------

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: wkspread    !  unused (jyg)

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg

  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar

  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw

  ! Tendencies of state variables (2 is appended to the names of fields which are the cumul of fields

  !                                 computed at each sub-timestep; e.g. d_wdens2 is the cumul of d_wdens)

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw

  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2

  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2

  ! Variables internes

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

  ! Variables a fixer

  REAL                                                  :: delta_t_min

  INTEGER                                               :: nsub

  REAL                                                  :: dtimesub

  REAL                                                  :: wdens0

  ! IM 080208

  ! Variables de sauvegarde

!!  REAL, DIMENSION (klon)                                :: sigmaw1

  ! Variables liees a la dynamique de population 1

  ! Variables liees a la dynamique de population 2

!!  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col

  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act

  REAL                                                  :: drdt_pos

  REAL                                                  :: tau_wk_inv_min

  ! Some components of the tendencies of state variables

  ! Variables pour les GW

  ! Variables liees au calcul de hw

  ! Variables liees au test de la forme triangulaire du profil de Delta_theta

  ! Sub-timestep tendencies and related variables

  ! Autres variables internes

  INTEGER                                               ::isubstep, k, i, igout

  REAL                                                  :: sigmaw_targ

  REAL                                                  :: wdens_targ

  REAL                                                  :: d_sigmaw_targ

  REAL                                                  :: d_wdens_targ

  REAL                                                  :: rrd1, rrd2

  ! cc nrlmd

!!  REAL, DIMENSION (klon)                                :: nat_rate

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

  ! Initialisations

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

  ! ALON = 3.e5

  ! alon = 1.E6

  !  Provisionnal; to be suppressed when f_shear is parameterized

  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)

  ! coefgw : Coefficient pour les ondes de gravite

  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)

  ! wdens : Densite surfacique de poche froide

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

  ! cc nrlmd      coefgw=10

  ! coefgw=1

  ! wdens0 = 1.0/(alon**2)

  ! cc nrlmd      wdens = 1.0/(alon**2)

  ! cc nrlmd      stark = 0.50

  ! CRtest

  ! cc nrlmd      alpk=0.1

  ! alpk = 1.0

  ! alpk = 0.5

  ! alpk = 0.05

!print *,'XXXX dtime input ', dtime

  ! Initialisation de toutes des densites a wdens_ref.

  ! Les densites peuvent evoluer si les poches debordent

  ! (voir au tout debut de la boucle sur les substeps)

  !jyg<

  !!  wdens(:) = wdens_ref

   ENDDO

  !>jyg

 ENDIF  ! (iflag_wk_pop_dyn == 0)

  ! print*,'stark',stark

  ! print*,'alpk',alpk

  ! print*,'wdens',wdens

  ! print*,'coefgw',coefgw

  ! cc

  ! Minimum value for |T_wake - T_undist|. Used for wake top definition

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

  delta_t_min = 0.2

  ! 1. - Save initial values, initialize tendencies, initialize output fields

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

!jyg<

!!  DO k = 1, klev

!!    DO i = 1, klon

!!      ppi(i, k) = pi(i, k)

!!      deltatw0(i, k) = deltatw(i, k)

!!      deltaqw0(i, k) = deltaqw(i, k)

!!      tenv(i, k) = tenv0(i, k)

!!      qe(i, k) = qe0(i, k)

!!      dtls(i, k) = 0.

!!      dqls(i, k) = 0.

!!      d_deltat_gw(i, k) = 0.

!!      d_tenv(i, k) = 0.

!!      d_qe(i, k) = 0.

!!      d_deltatw(i, k) = 0.

!!      d_deltaqw(i, k) = 0.

!!      ! IM 060508 beg

!!      d_deltatw2(i, k) = 0.

!!      d_deltaqw2(i, k) = 0.

!!      ! IM 060508 end

!!    END DO

!!  END DO

      ENDIF

!!  DO i = 1, klon

!!   sigmaw_in(i) = sigmaw(i)

!!  END DO

!>jyg

!

!!    wdens(:) = wdens(:) + wgen(:)*dtime

!!    d_wdens2(:) = wgen(:)*dtime

!!  ELSE

  ENDIF  ! (iflag_wk_pop_dyn >= 1)

  ! sigmaw1=sigmaw

  ! IF (sigd_con.GT.sigmaw1) THEN

  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con

  ! ENDIF

!

    END DO

!

       END DO

    ENDIF ! (iflag_wk_pop_dyn == 2)

  ELSE

    END DO

  ENDIF  ! (iflag_wk_pop_dyn >= 1)

!

    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))

!jyg<

!!    sigmaw(i) = amax1(sigmaw(i), sigmad)

!!    sigmaw(i) = amin1(sigmaw(i), 0.99)

!  print *,'XXXX1 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)

!>jyg

  END DO

!

!<jyg

!

!  Vertical advection local variables

!>jyg

!

  ENDIF

  ! 2. - Prognostic part

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

  ! 2.1 - Undisturbed area and Wake integrals

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

  END DO

  ! Distance between wakes

  END DO

  ! Potential temperatures and humidity

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

      ! write(*,*)'wake 1',i,k,RD,tenv(i,k)

      ! write(*,*)'wake 2',rho(i,k)

        ! write(*,*)'wake 3',i,k,rd,tenv(i,k)

        ! write(*,*)'wake 4',i,k,rd,tenv(i,k)

      ELSE

        ! write(*,*)'wake 5',rd,(tenv(i,k)+tenv(i,k-1))

        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)

      END IF

      ! write(*,*)'wake 7',ppi(i,k)

      ! write(*,*)'wake 8',(RD*(tenv(i,k)+deltatw(i,k)))

    END DO

  END DO

      ELSE

        n2(i, k) = amax1(0., -RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &

      END IF

    END DO

  END DO

  END DO

  ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)

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

    END DO

  END DO

  END DO

      rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &

    END DO

  END DO

  ! Choose an integration bound well above wake top

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

  ! Determine Wake top pressure (Ptop) from buoyancy integral

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

  ! -1/ Pressure of the level where dth becomes less than delta_t_min.

  END DO

      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)

        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &

      END IF

    END DO

  END DO

  ! -2/ dth integral

  END DO

      END IF

    END DO

  END DO

  ! -3/ height of triangle with area= sum_dth and base = dthmin

  END DO

  ! -4/ now, get Ptop

  END DO

      END IF

    END DO

  END DO

  ENDIF

  ! -5/ Determination de ktop et kupper

    END DO

  END DO

  !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper

  ! -6/ Correct ktop and ptop

  END DO

      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &

        ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &

      END IF

    END DO

  END DO

  END DO

    END DO

  END DO

  ENDIF

  ! -5/ Set deltatw & deltaqw to 0 above kupper

      END IF

    END DO

  END DO

  ! Vertical gradient of LS omega

      END IF

    END DO

  END DO

  ! Integrals (and wake top level number)

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

  ! Initialize sum_thvu to 1st level virt. pot. temp.

  END DO

      END IF

    END DO

  END DO

  END DO

  ! 2.1 - WAPE and mean forcing computation

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

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

  ! Means

    ! av_thve = sum_thve/hw0

    wape(i) = -RG*hw0(i)*(av_dth(i)+ &

  END DO

  ! 2.2 Prognostic variable update

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

  ! Filter out bad wakes

      END IF

    END DO

  END DO

!jyg<

!!      sigmaw(i) = amax1(sigmad, sigd_con(i))

!  print *,'XXXX2 d_sigmaw2(i), sigmaw(i) ', d_sigmaw2(i), sigmaw(i)

!>jyg

    ELSE

    END IF

  END DO

  ! Check qx and qw positivity

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

    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)),  &

  END DO

      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &

      END IF

    END DO

  END DO

  END DO

  ENDIF

  ! C -----------------------------------------------------------------

  ! Sub-time-stepping

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

  nsub = 10

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

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

    !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper

    ! wk_adv is the logical flag enabling wake evolution in the time advance

    ! loop

    END DO

    ENDIF

    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement

    ! negatif de ktop a kupper --------

    ! cc           On calcule pour cela une densite wdens0 pour laquelle on

    ! aurait un entrainement nul ---

    !jyg<

    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou

    ! les poches sont insuffisantes pour accueillir tout le flux de masse

    ! des descentes unsaturees. Nous faisons alors l'hypothese que la

    ! convection profonde cree directement de nouvelles poches, sans passer

    ! par les thermiques. La nouvelle valeur de wdens est alors imposee.

      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',

      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)

           omg(i, kupper(i)+1)=-RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &

        ELSE

        ENDIF

        wdens0 = (sigmaw(i)/(4.*3.14))* &

        ENDIF

          ENDIF

        END IF

      END IF

    END DO

!jyg<

!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)