Index: LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90	(revision 2435)
+++ LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90	(revision 2488)
@@ -7,4 +7,5 @@
 SUBROUTINE cv3_param(nd, k_upper, delt)
 
+  USE ioipsl_getin_p_mod, ONLY : getin_p
   use mod_phys_lmdz_para
   IMPLICIT NONE
@@ -39,5 +40,4 @@
   INTEGER, INTENT(IN)              :: k_upper
   REAL, INTENT(IN)                 :: delt ! timestep (seconds)
-
 
 ! Local variables
@@ -65,63 +65,54 @@
 
   IF (first) THEN
-
 ! -- "microphysical" parameters:
-    sigdz = 0.01
-    spfac = 0.15
-    pbcrit = 150.0
-    ptcrit = 500.0
 ! IM beg: ajout fis. reglage ep
-    flag_epkeorig = 1
-    elcrit = 0.0003
-    tlcrit = -55.0
+! CR+JYG: shedding coefficient (used when iflag_mix_adiab=1)
 ! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
 
     omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
-
 ! -- misc:
-
     dtovsh = -0.2 ! dT for overshoot
-    dpbase = -40. ! definition cloud base (400m above LCL)
 ! cc      dttrig = 5.   ! (loose) condition for triggering
     dttrig = 10. ! (loose) condition for triggering
-    flag_wb = 1
-    wbmax = 6. ! (m/s) adiab updraught speed at LFC (used in cv3p1_closure)
-
-! -- rate of approach to quasi-equilibrium:
-
     dtcrit = -2.0
-    tau = 8000.
-
 ! -- end of convection
-
-    tau_stop = 15000.
-    ok_convstop = .False.
-
-    ok_intermittent = .False.
-
 ! -- interface cloud parameterization:
-
     delta = 0.01 ! cld
-
 ! -- interface with boundary-layer (gust factor): (sb)
-
     betad = 10.0 ! original value (from convect 4.3)
 
-   !$OMP MASTER
-    OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', ERR=9999)
-    READ (99, *, END=9998) dpbase
-    READ (99, *, END=9998) pbcrit
-    READ (99, *, END=9998) ptcrit
-    READ (99, *, END=9998) sigdz
-    READ (99, *, END=9998) spfac
-    READ (99, *, END=9998) tau
-    READ (99, *, END=9998) flag_wb
-    READ (99, *, END=9998) wbmax
-    READ (99, *, END=9998) ok_convstop
-    READ (99, *, END=9998) tau_stop
-    READ (99, *, END=9998) ok_intermittent
-9998 CONTINUE
-    CLOSE (99)
-9999 CONTINUE
+! Var interm pour le getin
+     dpbase=-40.
+     CALL getin_p('dpbase',dpbase)
+     pbcrit=150.0
+     CALL getin_p('pbcrit',pbcrit)
+     ptcrit=500.0
+     CALL getin_p('ptcrit',ptcrit)
+     sigdz=0.01
+     CALL getin_p('sigdz',sigdz)
+     spfac=0.15
+     CALL getin_p('spfac',spfac)
+     tau=8000.
+     CALL getin_p('tau',tau)
+     flag_wb=1
+     CALL getin_p('flag_wb',flag_wb)
+     wbmax=6.
+     CALL getin_p('wbmax',wbmax)
+     ok_convstop=.False.
+     CALL getin_p('ok_convstop',ok_convstop)
+     tau_stop=15000.
+     CALL getin_p('tau_stop',tau_stop)
+     ok_intermittent=.False.
+     CALL getin_p('ok_intermittent',ok_intermittent)
+     coef_peel=0.25
+     CALL getin_p('coef_peel',coef_peel)
+
+     flag_epKEorig=1
+     CALL getin_p('flag_epKEorig',flag_epKEorig)
+     elcrit=0.0003
+     CALL getin_p('elcrit',elcrit)
+     tlcrit=-55.0
+     CALL getin_p('tlcrit',tlcrit)
+
     WRITE (*, *) 'dpbase=', dpbase
     WRITE (*, *) 'pbcrit=', pbcrit
@@ -130,42 +121,15 @@
     WRITE (*, *) 'spfac=', spfac
     WRITE (*, *) 'tau=', tau
-    WRITE (*, *) 'flag_wb =', flag_wb
-    WRITE (*, *) 'wbmax =', wbmax
-    WRITE (*, *) 'ok_convstop =', ok_convstop
-    WRITE (*, *) 'tau_stop =', tau_stop
-    WRITE (*, *) 'ok_intermittent =', ok_intermittent
-
-! IM Lecture du fichier ep_param.data
-    OPEN (79, FILE='ep_param.data', STATUS='old', FORM='formatted', ERR=7999)
-    READ (79, *, END=7998) flag_epkeorig
-    READ (79, *, END=7998) elcrit
-    READ (79, *, END=7998) tlcrit
-7998 CONTINUE
-    CLOSE (79)
-7999 CONTINUE
-    WRITE (*, *) 'flag_epKEorig', flag_epkeorig
+    WRITE (*, *) 'flag_wb=', flag_wb
+    WRITE (*, *) 'wbmax=', wbmax
+    WRITE (*, *) 'ok_convstop=', ok_convstop
+    WRITE (*, *) 'tau_stop=', tau_stop
+    WRITE (*, *) 'ok_intermittent=', ok_intermittent
+    WRITE (*, *) 'coef_peel=', coef_peel
+
+    WRITE (*, *) 'flag_epKEorig=', flag_epKEorig
     WRITE (*, *) 'elcrit=', elcrit
     WRITE (*, *) 'tlcrit=', tlcrit
-! IM end: ajout fis. reglage ep
-  !$OMP END MASTER
-
-   CALL bcast(dpbase)
-   CALL bcast(pbcrit)
-   CALL bcast(ptcrit)
-   CALL bcast(sigdz)
-   CALL bcast(spfac)
-   CALL bcast(tau)
-   CALL bcast(flag_wb)
-   CALL bcast(wbmax)
-   CALL bcast(ok_convstop)
-   CALL bcast(tau_stop)
-   CALL bcast(ok_intermittent)
-
-   CALL bcast(flag_epkeorig)
-   CALL bcast(elcrit)
-   CALL bcast(tlcrit)
-
     first = .FALSE.
-
   END IF ! (first)
 
@@ -4178,8 +4142,10 @@
                           ft, fq, fu, fv, ftra, &
                           Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
+                          epmax_diag, & ! epmax_cape
                           iflag1, &
                           precip1, sig1, w01, &
                           ft1, fq1, fu1, fv1, ftra1, &
-                          Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
+                          Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1, &
+                          epmax_diag1) ! epmax_cape
   IMPLICIT NONE
 
@@ -4198,4 +4164,5 @@
   REAL qcondc(nloc, nd)
   REAL wd(nloc), cape(nloc)
+  REAL epmax_diag(nloc)
 
 !outputs:
@@ -4209,4 +4176,5 @@
   REAL qcondc1(nloc, nd)
   REAL wd1(nloc), cape1(nloc)
+  REAL epmax_diag1(len) ! epmax_cape
 
 !local variables:
@@ -4218,4 +4186,5 @@
     wd1(idcum(i)) = wd(i)
     cape1(idcum(i)) = cape(i)
+    epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape
   END DO
 
@@ -4252,2 +4221,148 @@
   RETURN
 END SUBROUTINE cv3_uncompress
+
+
+        subroutine cv3_epmax_fn_cape(nloc,ncum,nd &
+                 , ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
+                 , pbase, p, ph, tv, buoy, sig, w0,iflag &
+                 , epmax_diag)
+        implicit none
+
+        ! On fait varier epmax en fn de la cape
+        ! Il faut donc recalculer ep, et hp qui a déjà été calculé et
+        ! qui en dépend
+        ! Toutes les autres variables fn de ep sont calculées plus bas.
+
+  include "cvthermo.h"
+  include "cv3param.h"  
+  include "conema3.h"
+  include "cvflag.h"
+
+! inputs:
+      INTEGER, INTENT (IN)                               :: ncum, nd, nloc
+      INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
+      REAL, DIMENSION (nloc), INTENT (IN)                :: hnk,pbase
+      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, lv, lf, tv, h
+      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: clw, buoy,frac
+      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig,w0
+      INTEGER, DIMENSION (nloc), INTENT (IN)             :: iflag(nloc)
+      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
+      REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
+! inouts:
+      REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: ep,hp  
+! outputs
+      REAL, DIMENSION (nloc), INTENT (OUT)           :: epmax_diag
+
+! local
+      integer i,k    
+!      real hp_bak(nloc,nd)
+!      real ep_bak(nloc,nd)
+      real m_loc(nloc,nd)
+      real sig_loc(nloc,nd)
+      real w0_loc(nloc,nd)
+      integer iflag_loc(nloc)
+      real cape(nloc)
+        
+        if (coef_epmax_cape.gt.1e-12) then
+
+        ! il faut calculer la cape: on fait un calcule simple car tant qu'on ne
+        ! connait pas ep, on ne connait pas les mélanges, ddfts etc... qui sont
+        ! necessaires au calcul de la cape dans la nouvelle physique
+        
+!        write(*,*) 'cv3_routines check 4303'
+        do i=1,ncum
+        do k=1,nd
+          sig_loc(i,k)=sig(i,k)
+          w0_loc(i,k)=w0(i,k)
+          iflag_loc(i)=iflag(i)
+!          ep_bak(i,k)=ep(i,k)
+        enddo ! do k=1,nd
+        enddo !do i=1,ncum
+
+!        write(*,*) 'cv3_routines check 4311'
+!        write(*,*) 'nl=',nl
+        CALL cv3_closure(nloc, ncum, nd, icb, inb, & ! na->nd
+          pbase, p, ph, tv, buoy, &
+          sig_loc, w0_loc, cape, m_loc,iflag_loc)
+
+!        write(*,*) 'cv3_routines check 4316'
+!        write(*,*) 'ep(1,:)=',ep(1,:)
+        do i=1,ncum
+           epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i))
+           epmax_diag(i)=amax1(epmax_diag(i),0.0)
+!           write(*,*) 'i,icb,inb,cape,epmax_diag=', &
+!                i,icb(i),inb(i),cape(i),epmax_diag(i)
+           do k=1,nl
+                ep(i,k)=ep(i,k)/epmax*epmax_diag(i)
+                ep(i,k)=amax1(ep(i,k),0.0)
+                ep(i,k)=amin1(ep(i,k),epmax_diag(i))
+           enddo
+        enddo
+ !       write(*,*) 'ep(1,:)=',ep(1,:)
+
+      !write(*,*) 'cv3_routines check 4326'
+! On recalcule hp:
+!      do k=1,nl
+!        do i=1,ncum
+!	  hp_bak(i,k)=hp(i,k)
+!	enddo
+!      enddo
+      do k=1,nl
+        do i=1,ncum
+	  hp(i,k)=h(i,k)
+	enddo
+      enddo
+
+  IF (cvflag_ice) THEN
+
+      do k=minorig+1,nl
+       do i=1,ncum
+        if((k.ge.icb(i)).and.(k.le.inb(i)))then
+          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* &
+                              ep(i, k)*clw(i, k)
+        endif
+       enddo
+      enddo !do k=minorig+1,n
+  ELSE !IF (cvflag_ice) THEN
+
+      DO k = minorig + 1, nl
+       DO i = 1, ncum
+        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
+          hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
+        endif
+       enddo
+      enddo !do k=minorig+1,n
+
+  ENDIF !IF (cvflag_ice) THEN     
+      !write(*,*) 'cv3_routines check 4345'
+!      do i=1,ncum  
+!       do k=1,nl
+!        if ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-1).or. &
+!            ((abs(hp_bak(i,k)-hp(i,k))/hp_bak(i,k).gt.1e-4).and. &
+!            (ep(i,k)-ep_bak(i,k).lt.1e-4))) then
+!           write(*,*) 'i,k=',i,k
+!           write(*,*) 'coef_epmax_cape=',coef_epmax_cape
+!           write(*,*) 'epmax_diag(i)=',epmax_diag(i)
+!           write(*,*) 'ep(i,k)=',ep(i,k)
+!           write(*,*) 'ep_bak(i,k)=',ep_bak(i,k)
+!           write(*,*) 'hp(i,k)=',hp(i,k)
+!           write(*,*) 'hp_bak(i,k)=',hp_bak(i,k)
+!           write(*,*) 'h(i,k)=',h(i,k)
+!           write(*,*) 'nk(i)=',nk(i)
+!           write(*,*) 'h(i,nk(i))=',h(i,nk(i))
+!           write(*,*) 'lv(i,k)=',lv(i,k)
+!           write(*,*) 't(i,k)=',t(i,k)
+!           write(*,*) 'clw(i,k)=',clw(i,k)
+!           write(*,*) 'cpd,cpv=',cpd,cpv
+!           stop
+!        endif
+!       enddo !do k=1,nl
+!      enddo !do i=1,ncum  
+      endif !if (coef_epmax_cape.gt.1e-12) then
+      !write(*,*) 'cv3_routines check 4367'
+
+      return
+      end subroutine cv3_epmax_fn_cape
+
+
+
