Changeset 2160 for LMDZ5/branches/testing/libf/phylmd/wake.F90

Timestamp:

Nov 28, 2014, 4:36:29 PM (10 years ago)

Author:

Laurent Fairhead

Message:

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phylmd/wake.F90 (modified) (8 diffs)

LMDZ5/branches/testing/libf/phylmd/wake.F90

@@ -19 +19 @@
 
   USE dimphy
+  use mod_phys_lmdz_para
   IMPLICIT NONE
   ! ============================================================================
@@ -151 +152 @@
   ! Variables à fixer
   REAL alon
-  REAL coefgw
-  REAL :: wdens_ref
-  REAL stark
-  REAL alpk
+  LOGICAL, SAVE :: first = .TRUE.
+  !$OMP THREADPRIVATE(first)
+  REAL, SAVE ::  stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol  
+  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
   REAL delta_t_min
   INTEGER nsub
@@ -232 +233 @@
 
   REAL, DIMENSION (klon, klev) :: crep
-  REAL crep_upper, crep_sol
 
   REAL, DIMENSION (klon, klev) :: ppi
@@ -280 +280 @@
   ! alpk = 0.05
 
+ if (first) then
   stark = 0.33
   alpk = 0.25
@@ -288 +289 @@
 
   ! cc nrlmd Lecture du fichier wake_param.data
+ !$OMP MASTER
   OPEN (99, FILE='wake_param.data', STATUS='old', FORM='formatted', ERR=9999)
   READ (99, *, END=9998) stark
@@ -296 +298 @@
   CLOSE (99)
 9999 CONTINUE
+ !$OMP END MASTER
+  CALL bcast(stark)
+  CALL bcast(alpk)
+  CALL bcast(wdens_ref)
+  CALL bcast(coefgw)
+
+  first=.false.
+ endif
 
   ! Initialisation de toutes des densites a wdens_ref.
@@ -1079 +1089 @@
           ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
 
-          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
-          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
+!jyg<
+!!
+!!---------------------------------------------------------------
+!! The change of delta_T due to PBL (vertical diffusion plus thermal plumes)
+!! is accounted for by the PBL and the Thermals schemes. It is now set to zero 
+!! within the Wake scheme.
+!!---------------------------------------------------------------
+          dtPBL(i,k) = 0.
+          dqPBL(i,k) = 0.
+! 
+!!          dtPBL(i,k)=wdtPBL(i,k) - udtPBL(i,k)
+!!          dqPBL(i,k)=wdqPBL(i,k) - udqPBL(i,k)
+! 
+!!          dtpbl(i, k) = (wdtpbl(i,k)/sigmaw(i)-udtpbl(i,k)/(1.-sigmaw(i)))
+!!          dqpbl(i, k) = (wdqpbl(i,k)/sigmaw(i)-udqpbl(i,k)/(1.-sigmaw(i)))
           ! print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
+!>jyg
+!
 
           ! cc nrlmd          Prise en compte du taux de mortalité
@@ -1795 +1820 @@
 END SUBROUTINE wake_vec_modulation
 
-SUBROUTINE wake_scal(p, ph, ppi, dtime, sigd_con, te0, qe0, omgb, dtdwn, &
-    dqdwn, amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, deltatw, &
-    deltaqw, dth, hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, &
-    dp_omgb, wdens, tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, &
-    spread, cstar, d_deltat_gw, d_deltatw2, d_deltaqw2)
-
-  ! **************************************************************
-  ! *
-  ! WAKE                                                        *
-  ! retour a un Pupper fixe                                *
-  ! *
-  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
-  ! modified by :   ROEHRIG Romain        01/29/2007            *
-  ! **************************************************************
-
-  USE dimphy
-  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.
-
-  ! Variables d'etat : deltatw    : ecart de temperature wake-undisturbed
-  ! area
-  ! deltaqw    : ecart d'humidite wake-undisturbed area
-  ! sigmaw     : fraction d'aire occupee par la poche.
-
-  ! 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   : hauteur de la poche
-  ! dp_omgb : vertical gradient of large scale omega
-  ! 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)
-  ! spread  : spreading term in dt_wake and dq_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
-  ! te0  : 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)
-  ! 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)
-  ! p    : pressions aux milieux des couches (Pa)
-  ! ph   : pressions aux interfaces (Pa)
-  ! ppi  : (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
-  ! te   : 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    : hw initial
-  ! sigmaw0: sigmaw initial
-  ! amflux : horizontal mass flux through wake boundary
-  ! wdens  : number of wakes per unit area (3D) or per
-  ! unit length (2D)
-  ! Tgw    : 1 sur la période de onde de gravité
-  ! Cgw    : vitesse de propagation de onde de gravité
-  ! LL     : distance entre 2 poches
-
-  ! -------------------------------------------------------------------------
-  ! Déclaration de variables
-  ! -------------------------------------------------------------------------
-
-  include "dimensions.h"
-  ! ccc      include "dimphy.h"
-  include "YOMCST.h"
-  include "cvthermo.h"
-  include "iniprint.h"
-
-  ! Arguments en entree
-  ! --------------------
-
-  REAL p(klev), ph(klev+1), ppi(klev)
-  REAL dtime
-  REAL te0(klev), qe0(klev)
-  REAL omgb(klev+1)
-  REAL dtdwn(klev), dqdwn(klev)
-  REAL wdtpbl(klev), wdqpbl(klev)
-  REAL udtpbl(klev), udqpbl(klev)
-  REAL amdwn(klev), amup(klev)
-  REAL dta(klev), dqa(klev)
-  REAL sigd_con
-
-  ! Sorties
-  ! --------
-
-  REAL deltatw(klev), deltaqw(klev), dth(klev)
-  REAL tu(klev), qu(klev)
-  REAL dtls(klev), dqls(klev)
-  REAL dtke(klev), dqke(klev)
-  REAL dtpbl(klev), dqpbl(klev)
-  REAL spread(klev)
-  REAL d_deltatgw(klev)
-  REAL d_deltatw2(klev), d_deltaqw2(klev)
-  REAL omgbdth(klev+1), omg(klev+1)
-  REAL dp_omgb(klev), dp_deltomg(klev)
-  REAL d_deltat_gw(klev)
-  REAL hw, sigmaw, wape, fip, gfl, cstar
-  INTEGER ktopw
-
-  ! Variables internes
-  ! -------------------
-
-  ! Variables à fixer
-  REAL alon
-  REAL coefgw
-  REAL wdens0, wdens
-  REAL stark
-  REAL alpk
-  REAL delta_t_min
-  REAL pupper
-  INTEGER nsub
-  REAL dtimesub
-  REAL sigmad, hwmin
-
-  ! Variables de sauvegarde
-  REAL deltatw0(klev)
-  REAL deltaqw0(klev)
-  REAL te(klev), qe(klev)
-  REAL sigmaw0, sigmaw1
-
-  ! Variables pour les GW
-  REAL ll
-  REAL n2(klev)
-  REAL cgw(klev)
-  REAL tgw(klev)
-
-  ! Variables liées au calcul de hw
-  REAL ptop_provis, ptop, ptop_new
-  REAL sum_dth
-  REAL dthmin
-  REAL z, dz, hw0
-  INTEGER ktop, kupper
-
-  ! Autres variables internes
-  INTEGER isubstep, k
-
-  REAL sum_thu, sum_tu, sum_qu, sum_thvu
-  REAL sum_dq, sum_rho
-  REAL sum_dtdwn, sum_dqdwn
-  REAL av_thu, av_tu, av_qu, av_thvu
-  REAL av_dth, av_dq, av_rho
-  REAL av_dtdwn, av_dqdwn
-
-  REAL rho(klev), rhoh(klev+1), rhow(klev)
-  REAL rhow_moyen(klev)
-  REAL zh(klev), zhh(klev+1)
-  REAL epaisseur1(klev), epaisseur2(klev)
-
-  REAL the(klev), thu(klev)
-
-  REAL d_deltatw(klev), d_deltaqw(klev)
-
-  REAL omgbw(klev+1), omgtop
-  REAL dp_omgbw(klev)
-  REAL ztop, dztop
-  REAL alpha_up(klev)
-
-  REAL rre1, rre2, rrd1, rrd2
-  REAL th1(klev), th2(klev), q1(klev), q2(klev)
-  REAL d_th1(klev), d_th2(klev), d_dth(klev)
-  REAL d_q1(klev), d_q2(klev), d_dq(klev)
-  REAL omgbdq(klev)
-
-  REAL ff, gg
-  REAL wape2, cstar2, heff
-
-  REAL crep(klev)
-  REAL crep_upper, crep_sol
-
-  ! -------------------------------------------------------------------------
-  ! Initialisations
-  ! -------------------------------------------------------------------------
-
-  ! print*, 'wake initialisations'
-
-  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
-  ! -------------------------------------------------------------------------
-
-  DATA sigmad, hwmin/.02, 10./
-
-  ! Longueur de maille (en m)
-  ! -------------------------------------------------------------------------
-
-  ! ALON = 3.e5
-  alon = 1.E6
-
-
-  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
-
-  ! coefgw : Coefficient pour les ondes de gravité
-  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
-  ! wdens : Densité de poche froide par maille
-  ! -------------------------------------------------------------------------
-
-  coefgw = 10
-  ! coefgw=1
-  ! wdens0 = 1.0/(alon**2)
-  wdens = 1.0/(alon**2)
-  stark = 0.50
-  ! CRtest
-  alpk = 0.1
-  ! alpk = 1.0
-  ! alpk = 0.5
-  ! alpk = 0.05
-  crep_upper = 0.9
-  crep_sol = 1.0
-
-
-  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
-  ! -------------------------------------------------------------------------
-
-  delta_t_min = 0.2
-
-
-  ! 1. - Save initial values and initialize tendencies
-  ! --------------------------------------------------
-
-  DO k = 1, klev
-    deltatw0(k) = deltatw(k)
-    deltaqw0(k) = deltaqw(k)
-    te(k) = te0(k)
-    qe(k) = qe0(k)
-    dtls(k) = 0.
-    dqls(k) = 0.
-    d_deltat_gw(k) = 0.
-    d_deltatw2(k) = 0.
-    d_deltaqw2(k) = 0.
-  END DO
-  ! sigmaw1=sigmaw
-  ! IF (sigd_con.GT.sigmaw1) THEN
-  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
-  ! ENDIF
-  sigmaw = max(sigmaw, sigd_con)
-  sigmaw = max(sigmaw, sigmad)
-  sigmaw = min(sigmaw, 0.99)
-  sigmaw0 = sigmaw
-  ! wdens=wdens0/(10.*sigmaw)
-  ! IF (sigd_con.GT.sigmaw1) THEN
-  ! print*, 'sigmaw1,sigd1', sigmaw, sigd_con
-  ! ENDIF
-
-  ! 2. - Prognostic part
-  ! =========================================================
-
-  ! print *, 'prognostic wake computation'
-
-
-  ! 2.1 - Undisturbed area and Wake integrals
-  ! ---------------------------------------------------------
-
-  z = 0.
-  ktop = 0
-  kupper = 0
-  sum_thu = 0.
-  sum_tu = 0.
-  sum_qu = 0.
-  sum_thvu = 0.
-  sum_dth = 0.
-  sum_dq = 0.
-  sum_rho = 0.
-  sum_dtdwn = 0.
-  sum_dqdwn = 0.
-
-  av_thu = 0.
-  av_tu = 0.
-  av_qu = 0.
-  av_thvu = 0.
-  av_dth = 0.
-  av_dq = 0.
-  av_rho = 0.
-  av_dtdwn = 0.
-  av_dqdwn = 0.
-
-  ! Potential temperatures and humidity
-  ! ----------------------------------------------------------
-
-  DO k = 1, klev
-    rho(k) = p(k)/(rd*te(k))
-    IF (k==1) THEN
-      rhoh(k) = ph(k)/(rd*te(k))
-      zhh(k) = 0
-    ELSE
-      rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
-      zhh(k) = (ph(k)-ph(k-1))/(-rhoh(k)*rg) + zhh(k-1)
-    END IF
-    the(k) = te(k)/ppi(k)
-    thu(k) = (te(k)-deltatw(k)*sigmaw)/ppi(k)
-    tu(k) = te(k) - deltatw(k)*sigmaw
-    qu(k) = qe(k) - deltaqw(k)*sigmaw
-    rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
-    dth(k) = deltatw(k)/ppi(k)
-    ll = (1-sqrt(sigmaw))/sqrt(wdens)
-  END DO
-
-  DO k = 1, klev - 1
-    IF (k==1) THEN
-      n2(k) = 0
-    ELSE
-      n2(k) = max(0., -rg**2/the(k)*rho(k)*(the(k+1)-the(k-1))/(p(k+ &
-        1)-p(k-1)))
-    END IF
-    zh(k) = (zhh(k)+zhh(k+1))/2
-
-    cgw(k) = sqrt(n2(k))*zh(k)
-    tgw(k) = coefgw*cgw(k)/ll
-  END DO
-
-  n2(klev) = 0
-  zh(klev) = 0
-  cgw(klev) = 0
-  tgw(klev) = 0
-
-  ! Calcul de la masse volumique moyenne de la colonne
-  ! -----------------------------------------------------------------
-
-  DO k = 1, klev
-    epaisseur1(k) = 0.
-    epaisseur2(k) = 0.
-  END DO
-
-  epaisseur1(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1.
-  epaisseur2(1) = -(ph(2)-ph(1))/(rho(1)*rg) + 1.
-  rhow_moyen(1) = rhow(1)
-
-  DO k = 2, klev
-    epaisseur1(k) = -(ph(k+1)-ph(k))/(rho(k)*rg) + 1.
-    epaisseur2(k) = epaisseur2(k-1) + epaisseur1(k)
-    rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+rhow(k)*epaisseur1(k))/ &
-      epaisseur2(k)
-  END DO
-
-
-  ! Choose an integration bound well above wake top
-  ! -----------------------------------------------------------------
-
-  ! Pupper = 50000.  ! melting level
-  pupper = 60000.
-  ! Pupper = 70000.
-
-
-  ! Determine Wake top pressure (Ptop) from buoyancy integral
-  ! -----------------------------------------------------------------
-
-  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
-
-  ptop_provis = ph(1)
-  DO k = 2, klev
-    IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
-      ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- &
-        1)+delta_t_min)*p(k))/(dth(k)-dth(k-1))
-      GO TO 25
-    END IF
-  END DO
-25 CONTINUE
-
-  ! -2/ dth integral
-
-  sum_dth = 0.
-  dthmin = -delta_t_min
-  z = 0.
-
-  DO k = 1, klev
-    dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
-    IF (dz<=0) GO TO 40
-    z = z + dz
-    sum_dth = sum_dth + dth(k)*dz
-    dthmin = min(dthmin, dth(k))
-  END DO
-40 CONTINUE
-
-  ! -3/ height of triangle with area= sum_dth and base = dthmin
-
-  hw0 = 2.*sum_dth/min(dthmin, -0.5)
-  hw0 = max(hwmin, hw0)
-
-  ! -4/ now, get Ptop
-
-  z = 0.
-  ptop = ph(1)
-
-  DO k = 1, klev
-    dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw0-z)
-    IF (dz<=0) GO TO 45
-    z = z + dz
-    ptop = ph(k) - rho(k)*rg*dz
-  END DO
-45 CONTINUE
-
-
-  ! -5/ Determination de ktop et kupper
-
-  DO k = klev, 1, -1
-    IF (ph(k+1)<ptop) ktop = k
-    IF (ph(k+1)<pupper) kupper = k
-  END DO
-
-  ! -6/ Correct ktop and ptop
-
-  ptop_new = ptop
-  DO k = ktop, 2, -1
-    IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
-      ptop_new = ((dth(k)+delta_t_min)*p(k-1)-(dth(k-1)+delta_t_min)*p(k))/ &
-        (dth(k)-dth(k-1))
-      GO TO 225
-    END IF
-  END DO
-225 CONTINUE
-
-  ptop = ptop_new
-
-  DO k = klev, 1, -1
-    IF (ph(k+1)<ptop) ktop = k
-  END DO
-
-  ! Set deltatw & deltaqw to 0 above kupper
-  ! -----------------------------------------------------------
-
-  DO k = kupper, klev
-    deltatw(k) = 0.
-    deltaqw(k) = 0.
-  END DO
-
-
-  ! Vertical gradient of LS omega
-  ! ------------------------------------------------------------
-
-  DO k = 1, kupper
-    dp_omgb(k) = (omgb(k+1)-omgb(k))/(ph(k+1)-ph(k))
-  END DO
-
-
-  ! Integrals (and wake top level number)
-  ! -----------------------------------------------------------
-
-  ! Initialize sum_thvu to 1st level virt. pot. temp.
-
-  z = 1.
-  dz = 1.
-  sum_thvu = thu(1)*(1.+eps*qu(1))*dz
-  sum_dth = 0.
-
-  DO k = 1, klev
-    dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
-    IF (dz<=0) GO TO 50
-    z = z + dz
-    sum_thu = sum_thu + thu(k)*dz
-    sum_tu = sum_tu + tu(k)*dz
-    sum_qu = sum_qu + qu(k)*dz
-    sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
-    sum_dth = sum_dth + dth(k)*dz
-    sum_dq = sum_dq + deltaqw(k)*dz
-    sum_rho = sum_rho + rhow(k)*dz
-    sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
-    sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
-  END DO
-50 CONTINUE
-
-  hw0 = z
-
-  ! 2.1 - WAPE and mean forcing computation
-  ! -------------------------------------------------------------
-
-  ! Means
-
-  av_thu = sum_thu/hw0
-  av_tu = sum_tu/hw0
-  av_qu = sum_qu/hw0
-  av_thvu = sum_thvu/hw0
-  ! av_thve = sum_thve/hw0
-  av_dth = sum_dth/hw0
-  av_dq = sum_dq/hw0
-  av_rho = sum_rho/hw0
-  av_dtdwn = sum_dtdwn/hw0
-  av_dqdwn = sum_dqdwn/hw0
-
-  wape = -rg*hw0*(av_dth+eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq))/ &
-    av_thvu
-
-  ! 2.2 Prognostic variable update
-  ! ------------------------------------------------------------
-
-  ! Filter out bad wakes
-
-  IF (wape<0.) THEN
-    IF (prt_level>=10) PRINT *, 'wape<0'
-    wape = 0.
-    hw = hwmin
-    sigmaw = max(sigmad, sigd_con)
-    fip = 0.
-    DO k = 1, klev
-      deltatw(k) = 0.
-      deltaqw(k) = 0.
-      dth(k) = 0.
-    END DO
-  ELSE
-    IF (prt_level>=10) PRINT *, 'wape>0'
-    cstar = stark*sqrt(2.*wape)
-  END IF
-
-  ! ------------------------------------------------------------------
-  ! Sub-time-stepping
-  ! ------------------------------------------------------------------
-
-  ! nsub=36
-  nsub = 10
-  dtimesub = dtime/nsub
-
-  ! ------------------------------------------------------------
-  DO isubstep = 1, nsub
-    ! ------------------------------------------------------------
-
-    ! print*,'---------------','substep=',isubstep,'-------------'
-
-    ! Evolution of sigmaw
-
-
-    gfl = 2.*sqrt(3.14*wdens*sigmaw)
-
-    sigmaw = sigmaw + gfl*cstar*dtimesub
-    sigmaw = min(sigmaw, 0.99) !!!!!!!!
-    ! wdens = wdens0/(10.*sigmaw)
-    ! sigmaw =max(sigmaw,sigd_con)
-    ! sigmaw =max(sigmaw,sigmad)
-
-    ! calcul de la difference de vitesse verticale poche - zone non perturbee
-
-    z = 0.
-    dp_deltomg(1:klev) = 0.
-    omg(1:klev+1) = 0.
-
-    omg(1) = 0.
-    dp_deltomg(1) = -(gfl*cstar)/(sigmaw*(1-sigmaw))
-
-    DO k = 2, ktop
-      dz = -(ph(k)-ph(k-1))/(rho(k-1)*rg)
-      z = z + dz
-      dp_deltomg(k) = dp_deltomg(1)
-      omg(k) = dp_deltomg(1)*z
-    END DO
-
-    dztop = -(ptop-ph(ktop))/(rho(ktop)*rg)
-    ztop = z + dztop
-    omgtop = dp_deltomg(1)*ztop
-
-
-    ! Conversion de la vitesse verticale de m/s a Pa/s
-
-    omgtop = -rho(ktop)*rg*omgtop
-    dp_deltomg(1) = omgtop/(ptop-ph(1))
-
-    DO k = 1, ktop
-      omg(k) = -rho(k)*rg*omg(k)
-      dp_deltomg(k) = dp_deltomg(1)
-    END DO
-
-    ! raccordement lineaire de omg de ptop a pupper
-
-    IF (kupper>ktop) THEN
-      omg(kupper+1) = -rg*amdwn(kupper+1)/sigmaw + rg*amup(kupper+1)/(1.- &
-        sigmaw)
-      dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(ptop-pupper)
-      DO k = ktop + 1, kupper
-        dp_deltomg(k) = dp_deltomg(kupper)
-        omg(k) = omgtop + (ph(k)-ptop)*dp_deltomg(kupper)
-      END DO
-    END IF
-
-    ! Compute wake average vertical velocity omgbw
-
-    DO k = 1, klev + 1
-      omgbw(k) = omgb(k) + (1.-sigmaw)*omg(k)
-    END DO
-
-    ! and its vertical gradient dp_omgbw
-
-    DO k = 1, klev
-      dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
-    END DO
-
-
-    ! Upstream coefficients for omgb velocity
-    ! --    (alpha_up(k) is the coefficient of the value at level k)
-    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
-
-    DO k = 1, klev
-      alpha_up(k) = 0.
-      IF (omgb(k)>0.) alpha_up(k) = 1.
-    END DO
-
-    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
-
-    rre1 = 1. - sigmaw
-    rre2 = sigmaw
-    rrd1 = -1.
-    rrd2 = 1.
-
-    ! Get [Th1,Th2], dth and [q1,q2]
-
-    DO k = 1, kupper + 1
-      dth(k) = deltatw(k)/ppi(k)
-      th1(k) = the(k) - sigmaw*dth(k) ! undisturbed area
-      th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake
-      q1(k) = qe(k) - sigmaw*deltaqw(k) ! undisturbed area
-      q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
-    END DO
-
-    d_th1(1) = 0.
-    d_th2(1) = 0.
-    d_dth(1) = 0.
-    d_q1(1) = 0.
-    d_q2(1) = 0.
-    d_dq(1) = 0.
-
-    DO k = 2, kupper + 1 !   loop on interfaces
-      d_th1(k) = th1(k-1) - th1(k)
-      d_th2(k) = th2(k-1) - th2(k)
-      d_dth(k) = dth(k-1) - dth(k)
-      d_q1(k) = q1(k-1) - q1(k)
-      d_q2(k) = q2(k-1) - q2(k)
-      d_dq(k) = deltaqw(k-1) - deltaqw(k)
-    END DO
-
-    omgbdth(1) = 0.
-    omgbdq(1) = 0.
-
-    DO k = 2, kupper + 1 !   loop on interfaces
-      omgbdth(k) = omgb(k)*(dth(k-1)-dth(k))
-      omgbdq(k) = omgb(k)*(deltaqw(k-1)-deltaqw(k))
-    END DO
-
-
-    ! -----------------------------------------------------------------
-    DO k = 1, kupper - 1
-      ! -----------------------------------------------------------------
-
-      ! Compute redistribution (advective) term
-
-      d_deltatw(k) = dtimesub/(ph(k)-ph(k+1))*(rrd1*omg(k)*sigmaw*d_th1(k)- &
-        rrd2*omg(k+1)*(1.-sigmaw)*d_th2(k+1)-(1.-alpha_up( &
-        k))*omgbdth(k)-alpha_up(k+1)*omgbdth(k+1))*ppi(k)
-      ! print*,'d_deltatw=',d_deltatw(k)
-
-      d_deltaqw(k) = dtimesub/(ph(k)-ph(k+1))*(rrd1*omg(k)*sigmaw*d_q1(k)- &
-        rrd2*omg(k+1)*(1.-sigmaw)*d_q2(k+1)-(1.-alpha_up( &
-        k))*omgbdq(k)-alpha_up(k+1)*omgbdq(k+1))
-      ! print*,'d_deltaqw=',d_deltaqw(k)
-
-      ! and increment large scale tendencies
-
-      dtls(k) = dtls(k) + dtimesub*((rre1*omg(k)*sigmaw*d_th1(k)-rre2*omg(k+ &
-        1)*(1.-sigmaw)*d_th2(k+1))/(ph(k)-ph(k+1))-sigmaw*(1.-sigmaw)*dth(k)* &
-        dp_deltomg(k))*ppi(k)
-      ! print*,'dtls=',dtls(k)
-
-      dqls(k) = dqls(k) + dtimesub*((rre1*omg(k)*sigmaw*d_q1(k)-rre2*omg(k+ &
-        1)*(1.-sigmaw)*d_q2(k+1))/(ph(k)-ph(k+1))-sigmaw*(1.-sigmaw)*deltaqw( &
-        k)*dp_deltomg(k))
-      ! print*,'dqls=',dqls(k)
-
-      ! -------------------------------------------------------------------
-    END DO
-    ! ------------------------------------------------------------------
-
-    ! Increment state variables
-
-    DO k = 1, kupper - 1
-
-      ! Coefficient de répartition
-
-      crep(k) = crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
-      crep(k) = crep(k) + crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
-
-
-      ! Reintroduce compensating subsidence term.
-
-      ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
-      ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
-      ! .                   /(1-sigmaw)
-      ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
-      ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
-      ! .                   /(1-sigmaw)
-
-      ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
-      ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
-      ! .                   /(1-sigmaw)
-      ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
-      ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
-      ! .                   /(1-sigmaw)
-
-      dtke(k) = (dtdwn(k)/sigmaw-dta(k)/(1.-sigmaw))
-      dqke(k) = (dqdwn(k)/sigmaw-dqa(k)/(1.-sigmaw))
-      ! print*,'dtKE=',dtKE(k)
-      ! print*,'dqKE=',dqKE(k)
-
-      dtpbl(k) = (wdtpbl(k)/sigmaw-udtpbl(k)/(1.-sigmaw))
-      dqpbl(k) = (wdqpbl(k)/sigmaw-udqpbl(k)/(1.-sigmaw))
-
-      spread(k) = (1.-sigmaw)*dp_deltomg(k) + gfl*cstar/sigmaw
-      ! print*,'spread=',spread(k)
-
-
-      ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
-      ! Jingmei
-
-      d_deltat_gw(k) = d_deltat_gw(k) - tgw(k)*deltatw(k)*dtimesub
-      ! print*,'d_delta_gw=',d_deltat_gw(k)
-      ff = d_deltatw(k)/dtimesub
-
-      ! Sans GW
-
-      ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
-
-      ! GW formule 1
-
-      ! deltatw(k) = deltatw(k)+dtimesub*
-      ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
-
-      ! GW formule 2
-
-      IF (dtimesub*tgw(k)<1.E-10) THEN
-        deltatw(k) = deltatw(k) + dtimesub*(ff+dtke(k)+dtpbl(k)-spread(k)* &
-          deltatw(k)-tgw(k)*deltatw(k))
-      ELSE
-        deltatw(k) = deltatw(k) + 1/tgw(k)*(1-exp(-dtimesub*tgw(k)))*(ff+dtke &
-          (k)+dtpbl(k)-spread(k)*deltatw(k)-tgw(k)*deltatw(k))
-      END IF
-
-      dth(k) = deltatw(k)/ppi(k)
-
-      gg = d_deltaqw(k)/dtimesub
-
-      deltaqw(k) = deltaqw(k) + dtimesub*(gg+dqke(k)+dqpbl(k)-spread(k)* &
-        deltaqw(k))
-
-      d_deltatw2(k) = d_deltatw2(k) + d_deltatw(k)
-      d_deltaqw2(k) = d_deltaqw2(k) + d_deltaqw(k)
-    END DO
-
-    ! And update large scale variables
-
-    DO k = 1, kupper
-      te(k) = te0(k) + dtls(k)
-      qe(k) = qe0(k) + dqls(k)
-      the(k) = te(k)/ppi(k)
-    END DO
-
-    ! Determine Ptop from buoyancy integral
-    ! ----------------------------------------------------------------------
-
-    ! -1/ Pressure of the level where dth changes sign.
-
-    ptop_provis = ph(1)
-
-    DO k = 2, klev
-      IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
-        ptop_provis = ((dth(k)+delta_t_min)*p(k-1)-(dth(k- &
-          1)+delta_t_min)*p(k))/(dth(k)-dth(k-1))
-        GO TO 65
-      END IF
-    END DO
-65  CONTINUE
-
-    ! -2/ dth integral
-
-    sum_dth = 0.
-    dthmin = -delta_t_min
-    z = 0.
-
-    DO k = 1, klev
-      dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
-      IF (dz<=0) GO TO 70
-      z = z + dz
-      sum_dth = sum_dth + dth(k)*dz
-      dthmin = min(dthmin, dth(k))
-    END DO
-70  CONTINUE
-
-    ! -3/ height of triangle with area= sum_dth and base = dthmin
-
-    hw = 2.*sum_dth/min(dthmin, -0.5)
-    hw = max(hwmin, hw)
-
-    ! -4/ now, get Ptop
-
-    ktop = 0
-    z = 0.
-
-    DO k = 1, klev
-      dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg), hw-z)
-      IF (dz<=0) GO TO 75
-      z = z + dz
-      ptop = ph(k) - rho(k)*rg*dz
-      ktop = k
-    END DO
-75  CONTINUE
-
-    ! -5/Correct ktop and ptop
-
-    ptop_new = ptop
-
-    DO k = ktop, 2, -1
-      IF (dth(k)>-delta_t_min .AND. dth(k-1)<-delta_t_min) THEN
-        ptop_new = ((dth(k)+delta_t_min)*p(k-1)-(dth(k-1)+delta_t_min)*p(k))/ &
-          (dth(k)-dth(k-1))
-        GO TO 275
-      END IF
-    END DO
-275 CONTINUE
-
-    ptop = ptop_new
-
-    DO k = klev, 1, -1
-      IF (ph(k+1)<ptop) ktop = k
-    END DO
-
-    ! -6/ Set deltatw & deltaqw to 0 above kupper
-
-    DO k = kupper, klev
-      deltatw(k) = 0.
-      deltaqw(k) = 0.
-    END DO
-
-    ! ------------------------------------------------------------------
-  END DO ! end sub-timestep loop
-  ! -----------------------------------------------------------------
-
-  ! Get back to tendencies per second
-
-  DO k = 1, kupper - 1
-    dtls(k) = dtls(k)/dtime
-    dqls(k) = dqls(k)/dtime
-    d_deltatw2(k) = d_deltatw2(k)/dtime
-    d_deltaqw2(k) = d_deltaqw2(k)/dtime
-    d_deltat_gw(k) = d_deltat_gw(k)/dtime
-  END DO
-
-  ! 2.1 - Undisturbed area and Wake integrals
-  ! ---------------------------------------------------------
-
-  z = 0.
-  sum_thu = 0.
-  sum_tu = 0.
-  sum_qu = 0.
-  sum_thvu = 0.
-  sum_dth = 0.
-  sum_dq = 0.
-  sum_rho = 0.
-  sum_dtdwn = 0.
-  sum_dqdwn = 0.
-
-  av_thu = 0.
-  av_tu = 0.
-  av_qu = 0.
-  av_thvu = 0.
-  av_dth = 0.
-  av_dq = 0.
-  av_rho = 0.
-  av_dtdwn = 0.
-  av_dqdwn = 0.
-
-  ! Potential temperatures and humidity
-  ! ----------------------------------------------------------
-
-  DO k = 1, klev
-    rho(k) = p(k)/(rd*te(k))
-    IF (k==1) THEN
-      rhoh(k) = ph(k)/(rd*te(k))
-      zhh(k) = 0
-    ELSE
-      rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
-      zhh(k) = (ph(k)-ph(k-1))/(-rhoh(k)*rg) + zhh(k-1)
-    END IF
-    the(k) = te(k)/ppi(k)
-    thu(k) = (te(k)-deltatw(k)*sigmaw)/ppi(k)
-    tu(k) = te(k) - deltatw(k)*sigmaw
-    qu(k) = qe(k) - deltaqw(k)*sigmaw
-    rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
-    dth(k) = deltatw(k)/ppi(k)
-
-  END DO
-
-  ! Integrals (and wake top level number)
-  ! -----------------------------------------------------------
-
-  ! Initialize sum_thvu to 1st level virt. pot. temp.
-
-  z = 1.
-  dz = 1.
-  sum_thvu = thu(1)*(1.+eps*qu(1))*dz
-  sum_dth = 0.
-
-  DO k = 1, klev
-    dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
-
-    IF (dz<=0) GO TO 51
-    z = z + dz
-    sum_thu = sum_thu + thu(k)*dz
-    sum_tu = sum_tu + tu(k)*dz
-    sum_qu = sum_qu + qu(k)*dz
-    sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
-    sum_dth = sum_dth + dth(k)*dz
-    sum_dq = sum_dq + deltaqw(k)*dz
-    sum_rho = sum_rho + rhow(k)*dz
-    sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
-    sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
-  END DO
-51 CONTINUE
-
-  hw0 = z
-
-  ! 2.1 - WAPE and mean forcing computation
-  ! -------------------------------------------------------------
-
-  ! Means
-
-  av_thu = sum_thu/hw0
-  av_tu = sum_tu/hw0
-  av_qu = sum_qu/hw0
-  av_thvu = sum_thvu/hw0
-  av_dth = sum_dth/hw0
-  av_dq = sum_dq/hw0
-  av_rho = sum_rho/hw0
-  av_dtdwn = sum_dtdwn/hw0
-  av_dqdwn = sum_dqdwn/hw0
-
-  wape2 = -rg*hw0*(av_dth+eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq))/ &
-    av_thvu
-
-
-  ! 2.2 Prognostic variable update
-  ! ------------------------------------------------------------
-
-  ! Filter out bad wakes
-
-  IF (wape2<0.) THEN
-    IF (prt_level>=10) PRINT *, 'wape2<0'
-    wape2 = 0.
-    hw = hwmin
-    sigmaw = max(sigmad, sigd_con)
-    fip = 0.
-    DO k = 1, klev
-      deltatw(k) = 0.
-      deltaqw(k) = 0.
-      dth(k) = 0.
-    END DO
-  ELSE
-    IF (prt_level>=10) PRINT *, 'wape2>0'
-    cstar2 = stark*sqrt(2.*wape2)
-
-  END IF
-
-  ktopw = ktop
-
-  IF (ktopw>0) THEN
-
-    ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
-    ! cc       heff = 600.
-    ! Utilisation de la hauteur hw
-    ! c       heff = 0.7*hw
-    heff = hw
-
-    fip = 0.5*rho(ktopw)*cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
-    fip = alpk*fip
-    ! jyg2
-  ELSE
-    fip = 0.
-  END IF
-
-
-  ! Limitation de sigmaw
-
-  ! sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
-  ! alors il disparait en se mélangeant à la partie undisturbed
-
-  ! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
-  IF ((sigmaw>0.9) .OR. ((wape>=wape2) .AND. (wape2<= &
-      1.0)) .OR. (ktopw<=2)) THEN
-    ! IM cf NR/JYG 251108    .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
-    ! IF (sigmaw.GT.0.9) THEN
-    DO k = 1, klev
-      dtls(k) = 0.
-      dqls(k) = 0.
-      deltatw(k) = 0.
-      deltaqw(k) = 0.
-    END DO
-    wape = 0.
-    hw = hwmin
-    sigmaw = sigmad
-    fip = 0.
-  END IF
-
-  RETURN
-END SUBROUTINE wake_scal
-
-
-
+