Index: LMDZ6/trunk/libf/phylmd/calwake.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calwake.F90	(revision 4084)
+++ LMDZ6/trunk/libf/phylmd/calwake.F90	(revision 4085)
@@ -208,14 +208,16 @@
   END DO
 
-  CALL wake(znatsurf, p, ph, pi, dtime, &
+
+
+  CALL wake(klon,klev,znatsurf, p, ph, pi, dtime, &
     te, qe, omgbe, &
     dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
     sigd0, Cin, &
-    dtw, dqw, sigmaw, awdens, wdens, &                                   ! state variables
+    dtw, dqw, sigmaw, awdens, wdens, &                      ! state variables
     dth, hw, wape, fip, gfl, &
     dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, &
     dtke, dqke, omg, dp_deltomg, spread, cstar, &
     d_deltat_gw, &
-    d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens)                     ! tendencies
+    d_deltatw, d_deltaqw, d_sigmaw, d_awdens, d_wdens)      ! tendencies
 
 !
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4084)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4085)
@@ -78,4 +78,5 @@
     USE write_field_phy
     USE lscp_mod, ONLY : lscp
+    USE wake_ini_mod, ONLY : wake_ini
 
     !USE cmp_seri_mod
@@ -1725,8 +1726,12 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! Nouvelle initialisation pour le rayonnement RRTM
-       !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        CALL iniradia(klon,klev,paprs(1,1:klev+1))
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+       CALL wake_ini(rg,rd,rv,prt_level)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
        !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -2145,4 +2150,5 @@
       sollwdown(:)= sollwdown(:) + betalwoff *(-1.*ZFLDN0(:,1) - &
                     sollwdown(:))
+
 
 
Index: LMDZ6/trunk/libf/phylmd/wake.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4084)
+++ LMDZ6/trunk/libf/phylmd/wake.F90	(revision 4085)
@@ -2,6 +2,7 @@
 ! $Id$
 
-SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
-                te0, qe0, omgb, &
+
+SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, &
+                tenv0, qe0, omgb, &
                 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
                 sigd_con, Cin, &
@@ -9,5 +10,5 @@
                 dth, hw, wape, fip, gfl, &
                 dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
-                dtke, dqke, omg, dp_deltomg, spread, cstar, &
+                dtke, dqke, omg, dp_deltomg, wkspread, cstar, &
                 d_deltat_gw, &
                 d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2)     ! tendencies
@@ -23,8 +24,13 @@
   ! **************************************************************
 
-  USE ioipsl_getin_p_mod, ONLY : getin_p
-  USE dimphy
-  use mod_phys_lmdz_para
-  USE print_control_mod, ONLY: prt_level
+
+  USE wake_ini_mod , ONLY : wake_ini
+  USE wake_ini_mod , ONLY : prt_level,epsim1,RG,RD
+  USE wake_ini_mod , ONLY : stark, wdens_ref, coefgw, alpk, wdens_ref, stark, coefgw, alpk
+  USE wake_ini_mod , ONLY : crep_upper, crep_sol, tau_cv, rzero, aa0, flag_wk_check_trgl
+  USE wake_ini_mod , ONLY : iflag_wk_act, iflag_wk_check_trgl, iflag_wk_pop_dyn, wdensmin
+  USE wake_ini_mod , ONLY : sigmad, hwmin, wapecut, cstart, sigmaw_max, dens_rate, epsilon_loc
+
+
   IMPLICIT NONE
   ! ============================================================================
@@ -57,5 +63,5 @@
   ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
   ! dp_deltomg  : vertical gradient of omg (s-1)
-  ! spread  : spreading term in d_t_wake and d_q_wake
+  ! 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).
@@ -66,5 +72,5 @@
 
   ! aire : aire de la maille
-  ! te0  : temperature dans l'environnement  (K)
+  ! tenv0  : temperature dans l'environnement  (K)
   ! qe0  : humidite dans l'environnement     (kg/kg)
   ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
@@ -90,5 +96,5 @@
   ! rho  : environment density at P levels
   ! rhoh : environment density at Ph levels
-  ! te   : environment temperature | may change within
+  ! tenv   : environment temperature | may change within
   ! qe   : environment humidity    | sub-time-stepping
   ! the  : environment potential temperature
@@ -121,10 +127,9 @@
   ! -------------------------------------------------------------------------
 
-  include "YOMCST.h"
-  include "cvthermo.h"
 
   ! Arguments en entree
   ! --------------------
 
+  INTEGER, INTENT(IN) :: klon,klev
   INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
@@ -132,5 +137,5 @@
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
   REAL,                             INTENT(IN)          :: dtime
-  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
+  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: tenv0, qe0
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
   REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
@@ -155,5 +160,5 @@
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
   REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
-  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
+  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
@@ -169,39 +174,8 @@
 
   ! Variables à fixer
-  INTEGER, SAVE                                         :: igout
-  !$OMP THREADPRIVATE(igout)
-  LOGICAL, SAVE                                         :: first = .TRUE.
-  !$OMP THREADPRIVATE(first)
-!jyg<
-!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
-  REAL, SAVE, DIMENSION(2)                              :: wdens_ref
-  REAL, SAVE                                            :: stark, coefgw, alpk
-!>jyg
-  REAL, SAVE                                            :: crep_upper, crep_sol  
-  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
-
-  REAL, SAVE                                            :: tau_cv
-  !$OMP THREADPRIVATE(tau_cv)
-  REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
-  !$OMP THREADPRIVATE(rzero, aa0)
-
-  LOGICAL, SAVE                                         :: flag_wk_check_trgl
-  !$OMP THREADPRIVATE(flag_wk_check_trgl)
-  INTEGER, SAVE                                         :: iflag_wk_check_trgl
-  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
-  INTEGER, SAVE                                         :: iflag_wk_pop_dyn
-  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
 
   REAL                                                  :: delta_t_min
   INTEGER                                               :: nsub
   REAL                                                  :: dtimesub
-  REAL, SAVE                                            :: wdensmin
-  !$OMP THREADPRIVATE(wdensmin)
-  REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
-  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
-  REAL, SAVE                                            :: sigmaw_max
-  !$OMP THREADPRIVATE(sigmaw_max)  
-  REAL, SAVE                                            :: dens_rate
-  !$OMP THREADPRIVATE(dens_rate)
   REAL                                                  :: wdens0
   ! IM 080208
@@ -211,5 +185,5 @@
   REAL, DIMENSION (klon, klev)                          :: deltatw0
   REAL, DIMENSION (klon, klev)                          :: deltaqw0
-  REAL, DIMENSION (klon, klev)                          :: te, qe
+  REAL, DIMENSION (klon, klev)                          :: tenv, qe
 !!  REAL, DIMENSION (klon)                                :: sigmaw1
 
@@ -222,6 +196,4 @@
   REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
   LOGICAL, DIMENSION (klon)                             :: kill_wake
-  INTEGER, SAVE                                         :: iflag_wk_act
-  !$OMP THREADPRIVATE(iflag_wk_act)
   REAL                                                  :: drdt_pos
   REAL                                                  :: tau_wk_inv_min
@@ -246,15 +218,12 @@
   ! Sub-timestep tendencies and related variables
   REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
-  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
+  REAL, DIMENSION (klon, klev)                          :: d_tenv, d_qe
   REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
   REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
   REAL, DIMENSION (klon)                                :: q0_min, q1_min
   LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
-  REAL, SAVE                                            :: epsilon
-  !$OMP THREADPRIVATE(epsilon)
-  DATA epsilon/1.E-15/
 
   ! Autres variables internes
-  INTEGER                                               ::isubstep, k, i
+  INTEGER                                               ::isubstep, k, i, igout
 
   REAL                                                  :: wdens_targ
@@ -307,28 +276,8 @@
   REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
 
+
   ! -------------------------------------------------------------------------
   ! Initialisations
   ! -------------------------------------------------------------------------
-
-  ! print*, 'wake initialisations'
-
-  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
-  ! -------------------------------------------------------------------------
-
-!!  DATA wapecut, sigmad, hwmin/5., .02, 10./
-!!  DATA wapecut, sigmad, hwmin/1., .02, 10./
-  DATA sigmad, hwmin/.02, 10./
-!!  DATA wdensmin/1.e-12/
-  DATA wdensmin/1.e-14/
-  ! cc nrlmd
-  DATA sigmaw_max/0.4/
-  DATA dens_rate/0.1/
-  ! cc
-  ! Longueur de maille (en m)
-  ! -------------------------------------------------------------------------
-
-  ! ALON = 3.e5
-  ! alon = 1.E6
-
   !  Provisionnal; to be suppressed when f_shear is parameterized
   f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
@@ -353,76 +302,6 @@
   ! alpk = 0.05
 
- if (first) then
-
-  igout = klon/2+1/klon
-
-  crep_upper = 0.9
-  crep_sol = 1.0
-
-  ! Get wapecut from parameter file
-  wapecut = 1.
-  CALL getin_p('wapecut', wapecut)
-
-  ! cc nrlmd Lecture du fichier wake_param.data
-  stark=0.33
-  CALL getin_p('stark',stark)
-  cstart = stark*sqrt(2.*wapecut)
-
-  alpk=0.25
-  CALL getin_p('alpk',alpk)
-!jyg<
-!!  wdens_ref=8.E-12
-!!  CALL getin_p('wdens_ref',wdens_ref)
-  wdens_ref(1)=8.E-12
-  wdens_ref(2)=8.E-12
-  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
-  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
-!>jyg
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!------------------------------------------------------------------------ 
-
-  iflag_wk_pop_dyn = 0
-  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 
-                                                    ! and wdens prognostic
-  iflag_wk_act = 0
-  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
-                                            ! 1: act(:)=1.
-                                            ! 2: act(:)=f(Wape)
-
-  rzero = 5000.
-  CALL getin_p('rzero_wk', rzero)
-  aa0 = 3.14*rzero*rzero
-!
-  tau_cv = 4000.
-  CALL getin_p('tau_cv', tau_cv)
-
-!------------------------------------------------------------------------ 
-
-  coefgw=4.
-  CALL getin_p('coefgw',coefgw)
-
-  WRITE(*,*) 'stark=', stark
-  WRITE(*,*) 'alpk=', alpk
-!jyg<
-!!  WRITE(*,*) 'wdens_ref=', wdens_ref
-  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
-  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
-!>jyg
-  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
-  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
-  WRITE(*,*) 'coefgw=', coefgw 
-
-  flag_wk_check_trgl=.false.
-  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
-  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
-  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
-  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
-  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
-  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
-
-  first=.false.
- endif
+
+ igout = klon/2+1/klon
 
  IF (iflag_wk_pop_dyn == 0) THEN
@@ -457,10 +336,10 @@
 !!      deltatw0(i, k) = deltatw(i, k)
 !!      deltaqw0(i, k) = deltaqw(i, k)
-!!      te(i, k) = te0(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_te(i, k) = 0.
+!!      d_tenv(i, k) = 0.
 !!      d_qe(i, k) = 0.
 !!      d_deltatw(i, k) = 0.
@@ -475,10 +354,10 @@
       deltatw0(:,:) = deltatw(:,:)
       deltaqw0(:,:) = deltaqw(:,:)
-      te(:,:) = te0(:,:)
+      tenv(:,:) = tenv0(:,:)
       qe(:,:) = qe0(:,:)
       dtls(:,:) = 0.
       dqls(:,:) = 0.
       d_deltat_gw(:,:) = 0.
-      d_te(:,:) = 0.
+      d_tenv(:,:) = 0.
       d_qe(:,:) = 0.
       d_deltatw(:,:) = 0.
@@ -547,5 +426,5 @@
 dtke(:,:) = 0.
 dqke(:,:) = 0.
-spread(:,:) = 0.
+wkspread(:,:) = 0.
 omgbdth(:,:) = 0.
 omg(:,:) = 0.
@@ -618,25 +497,25 @@
   DO k = 1, klev
     DO i = 1, klon
-      ! write(*,*)'wake 1',i,k,rd,te(i,k)
-      rho(i, k) = p(i, k)/(rd*te(i,k))
+      ! write(*,*)'wake 1',i,k,RD,tenv(i,k)
+      rho(i, k) = p(i, k)/(RD*tenv(i,k))
       ! write(*,*)'wake 2',rho(i,k)
       IF (k==1) THEN
-        ! write(*,*)'wake 3',i,k,rd,te(i,k)
-        rhoh(i, k) = ph(i, k)/(rd*te(i,k))
-        ! write(*,*)'wake 4',i,k,rd,te(i,k)
+        ! write(*,*)'wake 3',i,k,rd,tenv(i,k)
+        rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
+        ! write(*,*)'wake 4',i,k,rd,tenv(i,k)
         zhh(i, k) = 0
       ELSE
-        ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
-        rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
+        ! write(*,*)'wake 5',rd,(tenv(i,k)+tenv(i,k-1))
+        rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(i,k-1)))
         ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
-        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
+        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
       END IF
       ! write(*,*)'wake 7',ppi(i,k)
-      the(i, k) = te(i, k)/ppi(i, k)
-      thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
-      tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
+      the(i, k) = tenv(i, k)/ppi(i, k)
+      thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
+      tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
       qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
-      ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
-      rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
+      ! write(*,*)'wake 8',(RD*(tenv(i,k)+deltatw(i,k)))
+      rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
       dth(i, k) = deltatw(i, k)/ppi(i, k)
     END DO
@@ -648,5 +527,5 @@
         n2(i, k) = 0
       ELSE
-        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
+        n2(i, k) = amax1(0., -RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
                              (p(i,k+1)-p(i,k-1)))
       END IF
@@ -676,6 +555,6 @@
 
   DO i = 1, klon
-    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
-    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
+    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
+    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.
     rhow_moyen(i, 1) = rhow(i, 1)
   END DO
@@ -683,5 +562,5 @@
   DO k = 2, klev
     DO i = 1, klon
-      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
+      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG) + 1.
       epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
       rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
@@ -735,5 +614,5 @@
   DO k = 1, klev
     DO i = 1, klon
-      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
+      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
       IF (dz(i)>0) THEN
         z(i) = z(i) + dz(i)
@@ -760,8 +639,8 @@
   DO k = 1, klev
     DO i = 1, klon
-      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
+      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw0(i)-z(i))
       IF (dz(i)>0) THEN
         z(i) = z(i) + dz(i)
-        ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
+        ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
       END IF
     END DO
@@ -856,5 +735,5 @@
   DO k = 1, klev
     DO i = 1, klon
-      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
+      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
       IF (dz(i)>0) THEN
         z(i) = z(i) + dz(i)
@@ -896,5 +775,5 @@
     av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
 
-    wape(i) = -rg*hw0(i)*(av_dth(i)+ &
+    wape(i) = -RG*hw0(i)*(av_dth(i)+ &
         epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
 
@@ -1002,6 +881,6 @@
       ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
       IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
-        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
-                               rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
+        omg(i, kupper(i)+1) = -RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
+                               RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
         wdens0 = (sigmaw(i)/(4.*3.14))* &
           ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
@@ -1154,5 +1033,5 @@
       DO i = 1, klon
         IF (wk_adv(i) .AND. k<=ktop(i)) THEN
-          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
+          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*RG)
           z(i) = z(i) + dz(i)
           dp_deltomg(i, k) = dp_deltomg(i, 1)
@@ -1164,5 +1043,5 @@
     DO i = 1, klon
       IF (wk_adv(i)) THEN
-        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
+        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*RG)
         ztop(i) = z(i) + dztop(i)
         omgtop(i) = dp_deltomg(i, 1)*ztop(i)
@@ -1182,5 +1061,5 @@
     DO i = 1, klon
       IF (wk_adv(i)) THEN
-        omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
+        omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i)
         dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
       END IF
@@ -1190,5 +1069,5 @@
       DO i = 1, klon
         IF (wk_adv(i) .AND. k<=ktop(i)) THEN
-          omg(i, k) = -rho(i, k)*rg*omg(i, k)
+          omg(i, k) = -rho(i, k)*RG*omg(i, k)
           dp_deltomg(i, k) = dp_deltomg(i, 1)
         END IF
@@ -1200,6 +1079,6 @@
     DO i = 1, klon
       IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
-        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
-          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
+        omg(i, kupper(i)+1) = -RG*amdwn(i, kupper(i)+1)/sigmaw(i) + &
+          RG*amup(i, kupper(i)+1)/(1.-sigmaw(i))
         dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
           (ptop(i)-pupper(i))
@@ -1368,5 +1247,5 @@
           ! C
           ! -----------------------------------------------------------------
-          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
+          d_tenv(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
                                   rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
                                  (ph(i,k)-ph(i,k+1)) &
@@ -1380,5 +1259,5 @@
                                  (ph(i,k)-ph(i,k+1)) )
         ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
-          d_te(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
+          d_tenv(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
 
           d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
@@ -1469,7 +1348,7 @@
                          sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
 !>jyg
-            spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
-
-          ! cc        spread(i,k) =
+            wkspread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
+
+          ! cc        wkspread(i,k) =
           ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
           ! cc     $  sigmaw(i)
@@ -1488,10 +1367,10 @@
           ! Sans GW
 
-          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
+          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-wkspread(k)*deltatw(k))
 
           ! GW formule 1
 
           ! deltatw(k) = deltatw(k)+dtimesub*
-          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
+          ! $         (ff+dtKE(k) - wkspread(k)*deltatw(k)-Tgw(k)*deltatw(k))
 
           ! GW formule 2
@@ -1530,5 +1409,5 @@
     ! Scale tendencies so that water vapour remains positive in w and x.
 
-    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
+    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
       d_deltaqw, sigmaw, d_sigmaw, alpha)
 
@@ -1542,5 +1421,5 @@
       DO i = 1, klon
         IF (wk_adv(i) .AND. k<=kupper(i)) THEN
-          d_te(i, k) = alpha(i)*d_te(i, k)
+          d_tenv(i, k) = alpha(i)*d_tenv(i, k)
           d_qe(i, k) = alpha(i)*d_qe(i, k)
           d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
@@ -1562,5 +1441,5 @@
       DO i = 1, klon
         IF (wk_adv(i) .AND. k<=kupper(i)) THEN
-          dtls(i, k) = dtls(i, k) + d_te(i, k)
+          dtls(i, k) = dtls(i, k) + d_tenv(i, k)
           dqls(i, k) = dqls(i, k) + d_qe(i, k)
           ! cc nrlmd
@@ -1574,7 +1453,7 @@
       DO i = 1, klon
         IF (wk_adv(i) .AND. k<=kupper(i)) THEN
-          te(i, k) = te0(i, k) + dtls(i, k)
+          tenv(i, k) = tenv0(i, k) + dtls(i, k)
           qe(i, k) = qe0(i, k) + dqls(i, k)
-          the(i, k) = te(i, k)/ppi(i, k)
+          the(i, k) = tenv(i, k)/ppi(i, k)
           deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
           deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
@@ -1658,5 +1537,5 @@
       DO i = 1, klon
         IF (wk_adv(i)) THEN
-          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
+          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG)
           IF (dz(i)>0) THEN
             z(i) = z(i) + dz(i)
@@ -1689,8 +1568,8 @@
       DO i = 1, klon
         IF (wk_adv(i)) THEN
-          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
+          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw(i)-z(i))
           IF (dz(i)>0) THEN
             z(i) = z(i) + dz(i)
-            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
+            ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
             ktop(i) = k
           END IF
@@ -1789,5 +1668,5 @@
       DO i = 1, klon
         IF (wk_adv(i)) THEN !!! nrlmd
-          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
+          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
           IF (dz(i)>0) THEN
             z(i) = z(i) + dz(i)
@@ -1832,5 +1711,5 @@
         av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
 
-        wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
+        wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
                               av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
       END IF
@@ -1927,17 +1806,17 @@
       IF (ok_qx_qw(i)) THEN
         ! cc
-        rho(i, k) = p(i, k)/(rd*te(i,k))
+        rho(i, k) = p(i, k)/(RD*tenv(i,k))
         IF (k==1) THEN
-          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
+          rhoh(i, k) = ph(i, k)/(RD*tenv(i,k))
           zhh(i, k) = 0
         ELSE
-          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
-          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
-        END IF
-        the(i, k) = te(i, k)/ppi(i, k)
-        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
-        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
+          rhoh(i, k) = ph(i, k)*2./(RD*(tenv(i,k)+tenv(i,k-1)))
+          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1)
+        END IF
+        the(i, k) = tenv(i, k)/ppi(i, k)
+        thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
+        tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i)
         qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
-        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
+        rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k)))
         dth(i, k) = deltatw(i, k)/ppi(i, k)
       END IF
@@ -1967,6 +1846,6 @@
       IF (ok_qx_qw(i)) THEN
         ! cc
-        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
-        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg)
+        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
+        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*RG)
         IF (dz(i)>0) THEN
           z(i) = z(i) + dz(i)
@@ -2017,5 +1896,5 @@
       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
 
-      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
+      wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
                              av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
     END IF
@@ -2224,5 +2103,5 @@
     END DO
   END DO
-!jyg<
+
   DO i = 1, klon
     d_sigmaw2(i) = d_sigmaw2(i)/dtime
@@ -2230,12 +2109,9 @@
     d_wdens2(i) = d_wdens2(i)/dtime
   ENDDO
-!>jyg
-
-
-
-  RETURN
+
+ RETURN
 END SUBROUTINE wake
 
-SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
+SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qe, d_qe, deltaqw, &
     d_deltaqw, sigmaw, d_sigmaw, alpha)
   ! ------------------------------------------------------
@@ -2258,6 +2134,5 @@
   REAL alpha1(nlon)
   REAL x, a, b, c, discrim
-  REAL epsilon
-  ! DATA epsilon/1.e-15/
+  REAL epsilon_loc
   INTEGER i,k
 
@@ -2280,5 +2155,5 @@
         b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
           deltaqw(i, k)*d_sigmaw(i)
-        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
+        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc
         discrim = b*b - 4.*a*c
         ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
@@ -2308,6 +2183,2 @@
   RETURN
 END SUBROUTINE wake_vec_modulation
-
-
-
-
Index: LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90	(revision 4085)
+++ LMDZ6/trunk/libf/phylmd/wake_ini_mod.F90	(revision 4085)
@@ -0,0 +1,207 @@
+MODULE wake_ini_mod
+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
+
+  ! -------------------------------------------------------------------------
+  ! Déclaration de variables
+  ! -------------------------------------------------------------------------
+
+  ! Variables à fixer
+!jyg<
+!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
+  INTEGER, SAVE, PROTECTED                                         :: prt_level
+  REAL, SAVE, PROTECTED, DIMENSION(2)                              :: wdens_ref
+  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk
+!>jyg
+  REAL, SAVE, PROTECTED                                            :: crep_upper, crep_sol  
+  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
+
+  REAL, SAVE, PROTECTED                                            :: tau_cv
+  !$OMP THREADPRIVATE(tau_cv)
+  REAL, SAVE, PROTECTED                                            :: rzero, aa0 ! minimal wake radius and area
+  !$OMP THREADPRIVATE(rzero, aa0)
+
+  LOGICAL, SAVE, PROTECTED                                         :: flag_wk_check_trgl
+  !$OMP THREADPRIVATE(flag_wk_check_trgl)
+  INTEGER, SAVE, PROTECTED                                         :: iflag_wk_act
+  !$OMP THREADPRIVATE(iflag_wk_act)
+
+  INTEGER, SAVE, PROTECTED                                         :: iflag_wk_check_trgl
+  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
+  INTEGER, SAVE, PROTECTED                                         :: iflag_wk_pop_dyn
+  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
+
+  REAL, SAVE, PROTECTED                                            :: wdensmin
+  !$OMP THREADPRIVATE(wdensmin)
+  REAL, SAVE, PROTECTED                                            :: sigmad, hwmin, wapecut, cstart
+  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
+  REAL, SAVE, PROTECTED                                            :: sigmaw_max
+  !$OMP THREADPRIVATE(sigmaw_max)  
+  REAL, SAVE, PROTECTED                                            :: dens_rate
+  !$OMP THREADPRIVATE(dens_rate)
+  REAL, SAVE, PROTECTED                                            :: epsilon_loc
+  !$OMP THREADPRIVATE(epsilon_loc)
+  REAL, SAVE, PROTECTED                                            :: epsim1,RG,RD
+  !$OMP THREADPRIVATE(epsim1,RG,RD)
+
+
+
+CONTAINS
+
+  ! =========================================================================
+  SUBROUTINE wake_ini(rg_in,rd_in,rv_in,prt_lev)
+  ! =========================================================================
+
+  ! **************************************************************
+  ! *
+  ! WAKE                                                        *
+  ! retour a un Pupper fixe                                *
+  ! *
+  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
+  ! modified by :   ROEHRIG Romain        01/29/2007            *
+  ! **************************************************************
+
+  ! -------------------------------------------------------------------------
+  ! Initialisations
+  ! -------------------------------------------------------------------------
+
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  real eps
+  integer, intent(in) :: prt_lev
+  real, intent(in) :: rg_in,rd_in,rv_in
+
+  prt_level=prt_lev
+  epsilon_loc=1.E-15
+  wapecut=1. ! previously 5.
+  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
+  sigmad=0.02
+  hwmin=10.
+!!  DATA wdensmin/1.e-12/
+  wdensmin=1.e-14
+  ! cc nrlmd
+  sigmaw_max=0.4
+  dens_rate=0.1
+
+  eps = rd_in/rv_in
+  epsim1 = 1.0/eps - 1.0
+  RG=rg_in
+  RD=rd_in
+
+
+  ! cc
+  ! 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é 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
+
+
+
+  crep_upper = 0.9
+  crep_sol = 1.0
+
+  ! Get wapecut from parameter file
+  wapecut = 1.
+
+print*,'wapecut',wapecut
+  CALL getin_p('wapecut', wapecut)
+print*,'wapecut',wapecut
+
+  ! cc nrlmd Lecture du fichier wake_param.data
+
+
+  ! cc nrlmd Lecture du fichier wake_param.data
+  stark=0.33
+  CALL getin_p('stark',stark)
+  cstart = stark*sqrt(2.*wapecut)
+
+  alpk=0.25
+  CALL getin_p('alpk',alpk)
+!jyg<
+!!  wdens_ref=8.E-12
+!!  CALL getin_p('wdens_ref',wdens_ref)
+  wdens_ref(1)=8.E-12
+  wdens_ref(2)=8.E-12
+  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
+  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
+!>jyg
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!------------------------------------------------------------------------ 
+
+  iflag_wk_pop_dyn = 0
+  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed 
+                                                    ! and wdens prognostic
+  iflag_wk_act = 0
+  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
+                                            ! 1: act(:)=1.
+                                            ! 2: act(:)=f(Wape)
+
+  rzero = 5000.
+  CALL getin_p('rzero_wk', rzero)
+  aa0 = 3.14*rzero*rzero
+!
+  tau_cv = 4000.
+  CALL getin_p('tau_cv', tau_cv)
+
+!------------------------------------------------------------------------ 
+
+  coefgw=4.
+  CALL getin_p('coefgw',coefgw)
+
+  WRITE(*,*) 'stark=', stark
+  WRITE(*,*) 'alpk=', alpk
+!jyg<
+!!  WRITE(*,*) 'wdens_ref=', wdens_ref
+  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
+  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
+!>jyg
+  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
+  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
+  WRITE(*,*) 'coefgw=', coefgw 
+
+  flag_wk_check_trgl=.false.
+  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
+  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
+  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
+  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
+  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
+  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
+
+ RETURN
+
+END SUBROUTINE wake_ini
+
+END MODULE wake_ini_mod
