Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F	(revision 1470)
+++ 	(revision )
@@ -1,194 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-
-c    Auteur :  P. Le Van .
-c
-      IMPLICIT NONE
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-#include "logic.h"
-c
-c=======================================================================
-c
-c
-c    s = sigma ** kappa   :  coordonnee  verticale
-c    dsig(l)            : epaisseur de la couche l ds la coord.  s
-c    sig(l)             : sigma a l'interface des couches l et l-1
-c    ds(l)              : distance entre les couches l et l-1 en coord.s
-c
-c=======================================================================
-c
-      REAL pa,preff
-      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
-      REAL presnivs(llm)
-c
-c   declarations:
-c   -------------
-c
-      REAL sig(llm+1),dsig(llm)
-       real zzz(1:llm+1)
-       real dzz(1:llm)
-      real zk,zkm1,dzk1,dzk2,k0,k1
-c
-      INTEGER l
-      REAL snorm,dsigmin
-      REAL alpha,beta,gama,delta,deltaz,h
-      INTEGER np,ierr
-      REAL pi,x
-
-      REAL SSUM
-c
-c-----------------------------------------------------------------------
-c
-      pi=2.*ASIN(1.)
-
-      OPEN(99,file='sigma.def',status='old',form='formatted',
-     s   iostat=ierr)
-
-c-----------------------------------------------------------------------
-c   cas 1 on lit les options dans sigma.def:
-c   ----------------------------------------
-
-      IF (ierr.eq.0) THEN
-
-      READ(99,*) h           ! hauteur d'echelle 8.
-      READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
-      READ(99,*) beta        ! facteur d'acroissement en haut 1.3
-      READ(99,*) k0          ! nombre de couches dans la transition surf
-      READ(99,*) k1          ! nombre de couches dans la transition haute
-      CLOSE(99)
-      alpha=deltaz/(llm*h)
-      write(lunout,*)'h,alpha,k0,k1,beta'
-
-c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
-
-      alpha=deltaz/tanh(1./k0)*2.
-      zkm1=0.
-      sig(1)=1.
-      do l=1,llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
-     + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
-
-        dzk1=alpha*tanh(l/k0)
-        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
-        write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
-        zkm1=zk
-      enddo
-
-      sig(llm+1)=0.
-
-c
-       DO 2  l = 1, llm
-       dsig(l) = sig(l)-sig(l+1)
-   2   CONTINUE
-c
-
-      ELSE
-c-----------------------------------------------------------------------
-c   cas 2 ancienne discretisation (LMD5...):
-c   ----------------------------------------
-
-      WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
-
-      if (ok_strato) then
-         if (llm==39) then
-            dsigmin=0.3
-         else if (llm==50) then
-            dsigmin=1.
-         else
-            WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
-            dsigmin=1.
-         endif
-         WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
-      endif
-
-      h=7.
-      snorm  = 0.
-      DO l = 1, llm
-         x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
-
-         IF (ok_strato) THEN
-           dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
-     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
-         ELSE
-           dsig(l) = 1.0 + 7.0 * SIN(x)**2
-         ENDIF
-
-         snorm = snorm + dsig(l)
-      ENDDO
-      snorm = 1./snorm
-      DO l = 1, llm
-         dsig(l) = dsig(l)*snorm
-      ENDDO
-      sig(llm+1) = 0.
-      DO l = llm, 1, -1
-         sig(l) = sig(l+1) + dsig(l)
-      ENDDO
-
-      ENDIF
-
-
-      DO l=1,llm
-        nivsigs(l) = REAL(l)
-      ENDDO
-
-      DO l=1,llmp1
-        nivsig(l)= REAL(l)
-      ENDDO
-
-c
-c    ....  Calculs  de ap(l) et de bp(l)  ....
-c    .........................................
-c
-c
-c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
-c
-
-      bp(llmp1) =   0.
-
-      DO l = 1, llm
-cc
-ccc    ap(l) = 0.
-ccc    bp(l) = sig(l)
-
-      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-      ap(l) = pa * ( sig(l) - bp(l) )
-c
-      ENDDO
-
-      bp(1)=1.
-      ap(1)=0.
-
-      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
-
-      write(lunout,*)' BP '
-      write(lunout,*)  bp
-      write(lunout,*)' AP '
-      write(lunout,*)  ap
-
-      write(lunout,*)
-     .'Niveaux de pressions approximatifs aux centres des'
-      write(lunout,*)'couches calcules pour une pression de surface =',
-     .                 preff
-      write(lunout,*)
-     .     'et altitudes equivalentes pour une hauteur d echelle de'
-      write(lunout,*)'8km'
-      DO l = 1, llm
-       dpres(l) = bp(l) - bp(l+1)
-       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
-       write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
-     .        log(preff/presnivs(l))*8.
-     .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
-     .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
-      ENDDO
-
-      write(lunout,*)' PRESNIVS '
-      write(lunout,*)presnivs
-
-      RETURN
-      END
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90	(revision 1472)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90	(revision 1472)
@@ -0,0 +1,142 @@
+! $Id$
+
+SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+
+  ! Auteur : P. Le Van
+
+  IMPLICIT NONE
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  ! s = sigma ** kappa : coordonnee verticale
+  ! dsig(l) : epaisseur de la couche l ds la coord. s
+  ! sig(l) : sigma a l'interface des couches l et l-1
+  ! ds(l) : distance entre les couches l et l-1 en coord.s
+
+  REAL pa, preff
+  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
+  REAL presnivs(llm)
+
+  REAL sig(llm+1), dsig(llm)
+  real zk, zkm1, dzk1, dzk2, k0, k1
+
+  INTEGER l
+  REAL dsigmin
+  REAL alpha, beta, deltaz, h
+  INTEGER iostat
+  REAL pi, x
+
+  !-----------------------------------------------------------------------
+
+  pi = 2 * ASIN(1.)
+
+  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
+
+  IF (iostat == 0) THEN
+     ! cas 1 on lit les options dans sigma.def:
+     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
+     READ(99, *) beta ! facteur d'acroissement en haut 1.3
+     READ(99, *) k0 ! nombre de couches dans la transition surf
+     READ(99, *) k1 ! nombre de couches dans la transition haute
+     CLOSE(99)
+     alpha=deltaz/(llm*h)
+     write(lunout, *)'h, alpha, k0, k1, beta'
+
+     alpha=deltaz/tanh(1./k0)*2.
+     zkm1=0.
+     sig(1)=1.
+     do l=1, llm
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
+             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
+        zk=-h*log(sig(l+1))
+
+        dzk1=alpha*tanh(l/k0)
+        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
+        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
+        zkm1=zk
+     enddo
+
+     sig(llm+1)=0.
+
+     DO l = 1, llm
+        dsig(l) = sig(l)-sig(l+1)
+     end DO
+  ELSE
+     if (ok_strato) then
+        if (llm==39) then
+           dsigmin=0.3
+        else if (llm==50) then
+           dsigmin=1.
+        else
+           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
+           dsigmin=1.
+        endif
+        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
+     endif
+
+     DO l = 1, llm
+        x = pi * (l - 0.5) / (llm + 1)
+
+        IF (ok_strato) THEN
+           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
+        ELSE
+           dsig(l) = 1.0 + 7.0 * SIN(x)**2
+        ENDIF
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig(l) = sig(l+1) + dsig(l)
+     ENDDO
+  ENDIF
+
+  DO l=1, llm
+     nivsigs(l) = REAL(l)
+  ENDDO
+
+  DO l=1, llmp1
+     nivsig(l)= REAL(l)
+  ENDDO
+
+  ! .... Calculs de ap(l) et de bp(l) ....
+  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
+
+  bp(llmp1) = 0.
+
+  DO l = 1, llm
+     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+     ap(l) = pa * ( sig(l) - bp(l) )
+  ENDDO
+
+  bp(1)=1.
+  ap(1)=0.
+
+  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
+
+  write(lunout, *)' BP '
+  write(lunout, *) bp
+  write(lunout, *)' AP '
+  write(lunout, *) ap
+
+  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
+  write(lunout, *)'couches calcules pour une pression de surface =', preff
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
+  write(lunout, *)'8km'
+  DO l = 1, llm
+     dpres(l) = bp(l) - bp(l+1)
+     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
+          log(preff/presnivs(l))*8. &
+          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
+  ENDDO
+
+  write(lunout, *) 'PRESNIVS '
+  write(lunout, *) presnivs
+
+END SUBROUTINE disvert
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F	(revision 1470)
+++ 	(revision )
@@ -1,285 +1,0 @@
-!
-! $Id$
-!
-c
-c
-      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
-
-      USE filtreg_mod
-      USE infotrac, ONLY : nqtot
-      USE control_mod, ONLY: day_step,planet_type
-#ifdef CPP_IOIPSL
-      USE IOIPSL
-#else
-! if not using IOIPSL, we still need to use (a local version of) getin
-      USE ioipsl_getincom
-#endif
-      USE Write_Field
- 
-
-c%W%    %G%
-c=======================================================================
-c
-c   Author:    Frederic Hourdin      original: 15/01/93
-c   -------
-c
-c   Subject:
-c   ------
-c
-c   Method:
-c   --------
-c
-c   Interface:
-c   ----------
-c
-c      Input:
-c      ------
-c
-c      Output:
-c      -------
-c
-c=======================================================================
-      IMPLICIT NONE
-c-----------------------------------------------------------------------
-c   Declararations:
-c   ---------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "academic.h"
-#include "ener.h"
-#include "temps.h"
-#include "iniprint.h"
-#include "logic.h"
-
-c   Arguments:
-c   ----------
-
-      real time_0
-
-c   variables dynamiques
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
-      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
-      REAL ps(ip1jmp1)                       ! pression  au sol
-      REAL masse(ip1jmp1,llm)                ! masse d'air
-      REAL phis(ip1jmp1)                     ! geopotentiel au sol
-
-c   Local:
-c   ------
-
-      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
-      real tetajl(jjp1,llm)
-      INTEGER i,j,l,lsup,ij
-
-      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
-      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
-      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
-      LOGICAL ok_pv                ! Polar Vortex
-      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
-      
-      real zz,ran1
-      integer idum
-
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
-
-c-----------------------------------------------------------------------
-! 1. Initializations for Earth-like case
-! --------------------------------------
-c
-        ! initialize planet radius, rotation rate,...
-        call conf_planete
-
-        time_0=0.
-        day_ref=1
-        annee_ref=0
-
-        im         = iim
-        jm         = jjm
-        day_ini    = 1
-        dtvr    = daysec/REAL(day_step)
-        zdtvr=dtvr
-        etot0      = 0.
-        ptot0      = 0.
-        ztot0      = 0.
-        stot0      = 0.
-        ang0       = 0.
-
-        if (llm.eq.1) then
-          ! specific initializations for the shallow water case
-          kappa=1
-        endif
-        
-        CALL iniconst
-        CALL inigeom
-        CALL inifilr
-
-        if (llm.eq.1) then
-          ! initialize fields for the shallow water case, if required
-          if (.not.read_start) then
-            phis(:)=0.
-            q(:,:,:)=0
-            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
-          endif
-        endif
-
-        if (iflag_phys.eq.2) then
-          ! initializations for the academic case
-          
-!         if (planet_type=="earth") then
-
-          ! 1. local parameters
-          ! by convention, winter is in the southern hemisphere
-          ! Geostrophic wind or no wind?
-          ok_geost=.TRUE.
-          CALL getin('ok_geost',ok_geost)
-          ! Constants for Newtonian relaxation and friction
-          k_f=1.                !friction 
-          CALL getin('k_j',k_f)
-          k_f=1./(daysec*k_f)
-          k_c_s=4.  !cooling surface
-          CALL getin('k_c_s',k_c_s)
-          k_c_s=1./(daysec*k_c_s)
-          k_c_a=40. !cooling free atm
-          CALL getin('k_c_a',k_c_a)
-          k_c_a=1./(daysec*k_c_a)
-          ! Constants for Teta equilibrium profile
-          teta0=315.     ! mean Teta (S.H. 315K)
-          CALL getin('teta0',teta0)
-          ttp=200.       ! Tropopause temperature (S.H. 200K)
-          CALL getin('ttp',ttp)
-          eps=0.         ! Deviation to N-S symmetry(~0-20K)
-          CALL getin('eps',eps)
-          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
-          CALL getin('delt_y',delt_y)
-          delt_z=10.     ! Vertical Gradient (S.H. 10K)
-          CALL getin('delt_z',delt_z)
-          ! Polar vortex
-          ok_pv=.false.
-          CALL getin('ok_pv',ok_pv)
-          phi_pv=-50.            ! Latitude of edge of vortex
-          CALL getin('phi_pv',phi_pv)
-          phi_pv=phi_pv*pi/180.
-          dphi_pv=5.             ! Width of the edge
-          CALL getin('dphi_pv',dphi_pv)
-          dphi_pv=dphi_pv*pi/180.
-          gam_pv=4.              ! -dT/dz vortex (in K/km)
-          CALL getin('gam_pv',gam_pv)
-          
-          ! 2. Initialize fields towards which to relax
-          ! Friction
-          knewt_g=k_c_a
-          DO l=1,llm
-           zsig=presnivs(l)/preff
-           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
-           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
-          ENDDO
-          DO j=1,jjp1
-            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
-          ENDDO
-          
-          ! Potential temperature 
-          DO l=1,llm
-           zsig=presnivs(l)/preff
-           tetastrat=ttp*zsig**(-kappa)
-           tetapv=tetastrat
-           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
-               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
-           ENDIF
-           DO j=1,jjp1
-             ! Troposphere
-             ddsin=sin(rlatu(j))
-             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
-     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
-             ! Profil stratospherique isotherme (+vortex)
-             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
-             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
-             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
-           ENDDO
-          ENDDO ! of DO l=1,llm
- 
-!          CALL writefield('theta_eq',tetajl)
-
-          do l=1,llm
-            do j=1,jjp1
-              do i=1,iip1
-                 ij=(j-1)*iip1+i
-                 tetarappel(ij,l)=tetajl(j,l)
-              enddo
-            enddo
-          enddo
-
-
-!         else
-!          write(lunout,*)"iniacademic: planet types other than earth",
-!     &                   " not implemented (yet)."
-!          stop
-!         endif ! of if (planet_type=="earth")
-
-          ! 3. Initialize fields (if necessary)
-          IF (.NOT. read_start) THEN
-            ! surface pressure
-            ps(:)=preff
-            ! ground geopotential
-            phis(:)=0.
-            
-            CALL pression ( ip1jmp1, ap, bp, ps, p       )
-            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-            CALL massdair(p,masse)
-
-            ! bulk initialization of temperature
-            teta(:,:)=tetarappel(:,:)
-            
-            ! geopotential
-            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-            
-            ! winds
-            if (ok_geost) then
-              call ugeostr(phi,ucov)
-            else
-              ucov(:,:)=0.
-            endif
-            vcov(:,:)=0.
-            
-            ! bulk initialization of tracers
-            if (planet_type=="earth") then
-              ! Earth: first two tracers will be water
-              do i=1,nqtot
-                if (i.eq.1) q(:,:,i)=1.e-10
-                if (i.eq.2) q(:,:,i)=1.e-15
-                if (i.gt.2) q(:,:,i)=0.
-              enddo
-            else
-              q(:,:,:)=0
-            endif ! of if (planet_type=="earth")
-
-            ! add random perturbation to temperature
-            idum  = -1
-            zz = ran1(idum)
-            idum  = 0
-            do l=1,llm
-              do ij=iip2,ip1jm
-                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
-              enddo
-            enddo
-
-            ! maintain periodicity in longitude
-            do l=1,llm
-              do ij=1,ip1jmp1,iip1
-                teta(ij+iim,l)=teta(ij,l)
-              enddo
-            enddo
-
-          ENDIF ! of IF (.NOT. read_start)
-        endif ! of if (iflag_phys.eq.2)
-        
-      END
-c-----------------------------------------------------------------------
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90	(revision 1472)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90	(revision 1472)
@@ -0,0 +1,252 @@
+!
+! $Id$
+!
+SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+
+  USE filtreg_mod
+  USE infotrac, ONLY : nqtot
+  USE control_mod, ONLY: day_step,planet_type
+#ifdef CPP_IOIPSL
+  USE IOIPSL
+#else
+  ! if not using IOIPSL, we still need to use (a local version of) getin
+  USE ioipsl_getincom
+#endif
+  USE Write_Field
+
+  !   Author:    Frederic Hourdin      original: 15/01/93
+
+  IMPLICIT NONE
+
+  !   Declararations:
+  !   ---------------
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "comgeom.h"
+  include "academic.h"
+  include "ener.h"
+  include "temps.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  !   Arguments:
+  !   ----------
+
+  real time_0
+
+  !   variables dynamiques
+  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
+  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
+  REAL ps(ip1jmp1)                       ! pression  au sol
+  REAL masse(ip1jmp1,llm)                ! masse d'air
+  REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+  !   Local:
+  !   ------
+
+  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+  REAL pks(ip1jmp1)                      ! exner au  sol
+  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
+  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
+  REAL phi(ip1jmp1,llm)                  ! geopotentiel
+  REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
+  real tetajl(jjp1,llm)
+  INTEGER i,j,l,lsup,ij
+
+  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+  LOGICAL ok_pv                ! Polar Vortex
+  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+
+  real zz,ran1
+  integer idum
+
+  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+
+  !-----------------------------------------------------------------------
+  ! 1. Initializations for Earth-like case
+  ! --------------------------------------
+  !
+  ! initialize planet radius, rotation rate,...
+  call conf_planete
+
+  time_0=0.
+  day_ref=1
+  annee_ref=0
+
+  im         = iim
+  jm         = jjm
+  day_ini    = 1
+  dtvr    = daysec/REAL(day_step)
+  zdtvr=dtvr
+  etot0      = 0.
+  ptot0      = 0.
+  ztot0      = 0.
+  stot0      = 0.
+  ang0       = 0.
+
+  if (llm == 1) then
+     ! specific initializations for the shallow water case
+     kappa=1
+  endif
+
+  CALL iniconst
+  CALL inigeom
+  CALL inifilr
+
+  if (llm == 1) then
+     ! initialize fields for the shallow water case, if required
+     if (.not.read_start) then
+        phis(:)=0.
+        q(:,:,:)=0
+        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+     endif
+  endif
+
+  academic_case: if (iflag_phys == 2) then
+     ! initializations
+
+     ! 1. local parameters
+     ! by convention, winter is in the southern hemisphere
+     ! Geostrophic wind or no wind?
+     ok_geost=.TRUE.
+     CALL getin('ok_geost',ok_geost)
+     ! Constants for Newtonian relaxation and friction
+     k_f=1.                !friction 
+     CALL getin('k_j',k_f)
+     k_f=1./(daysec*k_f)
+     k_c_s=4.  !cooling surface
+     CALL getin('k_c_s',k_c_s)
+     k_c_s=1./(daysec*k_c_s)
+     k_c_a=40. !cooling free atm
+     CALL getin('k_c_a',k_c_a)
+     k_c_a=1./(daysec*k_c_a)
+     ! Constants for Teta equilibrium profile
+     teta0=315.     ! mean Teta (S.H. 315K)
+     CALL getin('teta0',teta0)
+     ttp=200.       ! Tropopause temperature (S.H. 200K)
+     CALL getin('ttp',ttp)
+     eps=0.         ! Deviation to N-S symmetry(~0-20K)
+     CALL getin('eps',eps)
+     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+     CALL getin('delt_y',delt_y)
+     delt_z=10.     ! Vertical Gradient (S.H. 10K)
+     CALL getin('delt_z',delt_z)
+     ! Polar vortex
+     ok_pv=.false.
+     CALL getin('ok_pv',ok_pv)
+     phi_pv=-50.            ! Latitude of edge of vortex
+     CALL getin('phi_pv',phi_pv)
+     phi_pv=phi_pv*pi/180.
+     dphi_pv=5.             ! Width of the edge
+     CALL getin('dphi_pv',dphi_pv)
+     dphi_pv=dphi_pv*pi/180.
+     gam_pv=4.              ! -dT/dz vortex (in K/km)
+     CALL getin('gam_pv',gam_pv)
+
+     ! 2. Initialize fields towards which to relax
+     ! Friction
+     knewt_g=k_c_a
+     DO l=1,llm
+        zsig=presnivs(l)/preff
+        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+     ENDDO
+     DO j=1,jjp1
+        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+     ENDDO
+
+     ! Potential temperature 
+     DO l=1,llm
+        zsig=presnivs(l)/preff
+        tetastrat=ttp*zsig**(-kappa)
+        tetapv=tetastrat
+        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+        ENDIF
+        DO j=1,jjp1
+           ! Troposphere
+           ddsin=sin(rlatu(j))
+           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
+                -delt_z*(1.-ddsin*ddsin)*log(zsig)
+           ! Profil stratospherique isotherme (+vortex)
+           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+           tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+        ENDDO
+     ENDDO
+
+     !          CALL writefield('theta_eq',tetajl)
+
+     do l=1,llm
+        do j=1,jjp1
+           do i=1,iip1
+              ij=(j-1)*iip1+i
+              tetarappel(ij,l)=tetajl(j,l)
+           enddo
+        enddo
+     enddo
+
+     ! 3. Initialize fields (if necessary)
+     IF (.NOT. read_start) THEN
+        ! surface pressure
+        ps(:)=preff
+        ! ground geopotential
+        phis(:)=0.
+
+        CALL pression ( ip1jmp1, ap, bp, ps, p       )
+        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL massdair(p,masse)
+
+        ! bulk initialization of temperature
+        teta(:,:)=tetarappel(:,:)
+
+        ! geopotential
+        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+
+        ! winds
+        if (ok_geost) then
+           call ugeostr(phi,ucov)
+        else
+           ucov(:,:)=0.
+        endif
+        vcov(:,:)=0.
+
+        ! bulk initialization of tracers
+        if (planet_type=="earth") then
+           ! Earth: first two tracers will be water
+           do i=1,nqtot
+              if (i == 1) q(:,:,i)=1.e-10
+              if (i == 2) q(:,:,i)=1.e-15
+              if (i.gt.2) q(:,:,i)=0.
+           enddo
+        else
+           q(:,:,:)=0
+        endif ! of if (planet_type=="earth")
+
+        ! add random perturbation to temperature
+        idum  = -1
+        zz = ran1(idum)
+        idum  = 0
+        do l=1,llm
+           do ij=iip2,ip1jm
+              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+           enddo
+        enddo
+
+        ! maintain periodicity in longitude
+        do l=1,llm
+           do ij=1,ip1jmp1,iip1
+              teta(ij+iim,l)=teta(ij,l)
+           enddo
+        enddo
+
+     ENDIF ! of IF (.NOT. read_start)
+  endif academic_case
+
+END SUBROUTINE iniacademic
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F	(revision 1470)
+++ 	(revision )
@@ -1,194 +1,0 @@
-!
-! $Id$
-!
-      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
-
-c    Auteur :  P. Le Van .
-c
-      IMPLICIT NONE
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
-#include "logic.h"
-c
-c=======================================================================
-c
-c
-c    s = sigma ** kappa   :  coordonnee  verticale
-c    dsig(l)            : epaisseur de la couche l ds la coord.  s
-c    sig(l)             : sigma a l'interface des couches l et l-1
-c    ds(l)              : distance entre les couches l et l-1 en coord.s
-c
-c=======================================================================
-c
-      REAL pa,preff
-      REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
-      REAL presnivs(llm)
-c
-c   declarations:
-c   -------------
-c
-      REAL sig(llm+1),dsig(llm)
-       real zzz(1:llm+1)
-       real dzz(1:llm)
-      real zk,zkm1,dzk1,dzk2,k0,k1
-c
-      INTEGER l
-      REAL snorm,dsigmin
-      REAL alpha,beta,gama,delta,deltaz,h
-      INTEGER np,ierr
-      REAL pi,x
-
-      REAL SSUM
-c
-c-----------------------------------------------------------------------
-c
-      pi=2.*ASIN(1.)
-
-      OPEN(99,file='sigma.def',status='old',form='formatted',
-     s   iostat=ierr)
-
-c-----------------------------------------------------------------------
-c   cas 1 on lit les options dans sigma.def:
-c   ----------------------------------------
-
-      IF (ierr.eq.0) THEN
-
-      READ(99,*) h           ! hauteur d'echelle 8.
-      READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
-      READ(99,*) beta        ! facteur d'acroissement en haut 1.3
-      READ(99,*) k0          ! nombre de couches dans la transition surf
-      READ(99,*) k1          ! nombre de couches dans la transition haute
-      CLOSE(99)
-      alpha=deltaz/(llm*h)
-      write(lunout,*)'h,alpha,k0,k1,beta'
-
-c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
-
-      alpha=deltaz/tanh(1./k0)*2.
-      zkm1=0.
-      sig(1)=1.
-      do l=1,llm
-        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
-     + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
-        zk=-h*log(sig(l+1))
-
-        dzk1=alpha*tanh(l/k0)
-        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
-        write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
-        zkm1=zk
-      enddo
-
-      sig(llm+1)=0.
-
-c
-       DO 2  l = 1, llm
-       dsig(l) = sig(l)-sig(l+1)
-   2   CONTINUE
-c
-
-      ELSE
-c-----------------------------------------------------------------------
-c   cas 2 ancienne discretisation (LMD5...):
-c   ----------------------------------------
-
-      WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
-
-      if (ok_strato) then
-         if (llm==39) then
-            dsigmin=0.3
-         else if (llm==50) then
-            dsigmin=1.
-         else
-            WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
-            dsigmin=1.
-         endif
-         WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
-      endif
-
-      h=7.
-      snorm  = 0.
-      DO l = 1, llm
-         x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
-
-         IF (ok_strato) THEN
-           dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
-     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
-         ELSE
-           dsig(l) = 1.0 + 7.0 * SIN(x)**2
-         ENDIF
-
-         snorm = snorm + dsig(l)
-      ENDDO
-      snorm = 1./snorm
-      DO l = 1, llm
-         dsig(l) = dsig(l)*snorm
-      ENDDO
-      sig(llm+1) = 0.
-      DO l = llm, 1, -1
-         sig(l) = sig(l+1) + dsig(l)
-      ENDDO
-
-      ENDIF
-
-
-      DO l=1,llm
-        nivsigs(l) = REAL(l)
-      ENDDO
-
-      DO l=1,llmp1
-        nivsig(l)= REAL(l)
-      ENDDO
-
-c
-c    ....  Calculs  de ap(l) et de bp(l)  ....
-c    .........................................
-c
-c
-c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
-c
-
-      bp(llmp1) =   0.
-
-      DO l = 1, llm
-cc
-ccc    ap(l) = 0.
-ccc    bp(l) = sig(l)
-
-      bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
-      ap(l) = pa * ( sig(l) - bp(l) )
-c
-      ENDDO
-
-      bp(1)=1.
-      ap(1)=0.
-
-      ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
-
-      write(lunout,*)' BP '
-      write(lunout,*)  bp
-      write(lunout,*)' AP '
-      write(lunout,*)  ap
-
-      write(lunout,*)
-     .'Niveaux de pressions approximatifs aux centres des'
-      write(lunout,*)'couches calcules pour une pression de surface =',
-     .                 preff
-      write(lunout,*)
-     .     'et altitudes equivalentes pour une hauteur d echelle de'
-      write(lunout,*)'8km'
-      DO l = 1, llm
-       dpres(l) = bp(l) - bp(l+1)
-       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
-       write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
-     .        log(preff/presnivs(l))*8.
-     .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
-     .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
-      ENDDO
-
-      write(lunout,*)' PRESNIVS '
-      write(lunout,*)presnivs
-
-      RETURN
-      END
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90	(revision 1472)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90	(revision 1472)
@@ -0,0 +1,142 @@
+! $Id$
+
+SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
+
+  ! Auteur : P. Le Van
+
+  IMPLICIT NONE
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  ! s = sigma ** kappa : coordonnee verticale
+  ! dsig(l) : epaisseur de la couche l ds la coord. s
+  ! sig(l) : sigma a l'interface des couches l et l-1
+  ! ds(l) : distance entre les couches l et l-1 en coord.s
+
+  REAL pa, preff
+  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
+  REAL presnivs(llm)
+
+  REAL sig(llm+1), dsig(llm)
+  real zk, zkm1, dzk1, dzk2, k0, k1
+
+  INTEGER l
+  REAL dsigmin
+  REAL alpha, beta, deltaz, h
+  INTEGER iostat
+  REAL pi, x
+
+  !-----------------------------------------------------------------------
+
+  pi = 2 * ASIN(1.)
+
+  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
+
+  IF (iostat == 0) THEN
+     ! cas 1 on lit les options dans sigma.def:
+     READ(99, *) h ! hauteur d'echelle 8.
+     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
+     READ(99, *) beta ! facteur d'acroissement en haut 1.3
+     READ(99, *) k0 ! nombre de couches dans la transition surf
+     READ(99, *) k1 ! nombre de couches dans la transition haute
+     CLOSE(99)
+     alpha=deltaz/(llm*h)
+     write(lunout, *)'h, alpha, k0, k1, beta'
+
+     alpha=deltaz/tanh(1./k0)*2.
+     zkm1=0.
+     sig(1)=1.
+     do l=1, llm
+        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
+             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
+        zk=-h*log(sig(l+1))
+
+        dzk1=alpha*tanh(l/k0)
+        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
+        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
+        zkm1=zk
+     enddo
+
+     sig(llm+1)=0.
+
+     DO l = 1, llm
+        dsig(l) = sig(l)-sig(l+1)
+     end DO
+  ELSE
+     if (ok_strato) then
+        if (llm==39) then
+           dsigmin=0.3
+        else if (llm==50) then
+           dsigmin=1.
+        else
+           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
+           dsigmin=1.
+        endif
+        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
+     endif
+
+     DO l = 1, llm
+        x = pi * (l - 0.5) / (llm + 1)
+
+        IF (ok_strato) THEN
+           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
+                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
+        ELSE
+           dsig(l) = 1.0 + 7.0 * SIN(x)**2
+        ENDIF
+     ENDDO
+     dsig = dsig / sum(dsig)
+     sig(llm+1) = 0.
+     DO l = llm, 1, -1
+        sig(l) = sig(l+1) + dsig(l)
+     ENDDO
+  ENDIF
+
+  DO l=1, llm
+     nivsigs(l) = REAL(l)
+  ENDDO
+
+  DO l=1, llmp1
+     nivsig(l)= REAL(l)
+  ENDDO
+
+  ! .... Calculs de ap(l) et de bp(l) ....
+  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
+
+  bp(llmp1) = 0.
+
+  DO l = 1, llm
+     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
+     ap(l) = pa * ( sig(l) - bp(l) )
+  ENDDO
+
+  bp(1)=1.
+  ap(1)=0.
+
+  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
+
+  write(lunout, *)' BP '
+  write(lunout, *) bp
+  write(lunout, *)' AP '
+  write(lunout, *) ap
+
+  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
+  write(lunout, *)'couches calcules pour une pression de surface =', preff
+  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
+  write(lunout, *)'8km'
+  DO l = 1, llm
+     dpres(l) = bp(l) - bp(l+1)
+     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
+          log(preff/presnivs(l))*8. &
+          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
+          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
+  ENDDO
+
+  write(lunout, *) 'PRESNIVS '
+  write(lunout, *) presnivs
+
+END SUBROUTINE disvert
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F	(revision 1470)
+++ 	(revision )
@@ -1,285 +1,0 @@
-!
-! $Id$
-!
-c
-c
-      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
-
-      USE filtreg_mod
-      USE infotrac, ONLY : nqtot
-      USE control_mod, ONLY: day_step,planet_type
-#ifdef CPP_IOIPSL
-      USE IOIPSL
-#else
-! if not using IOIPSL, we still need to use (a local version of) getin
-      USE ioipsl_getincom
-#endif
-      USE Write_Field
- 
-
-c%W%    %G%
-c=======================================================================
-c
-c   Author:    Frederic Hourdin      original: 15/01/93
-c   -------
-c
-c   Subject:
-c   ------
-c
-c   Method:
-c   --------
-c
-c   Interface:
-c   ----------
-c
-c      Input:
-c      ------
-c
-c      Output:
-c      -------
-c
-c=======================================================================
-      IMPLICIT NONE
-c-----------------------------------------------------------------------
-c   Declararations:
-c   ---------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "academic.h"
-#include "ener.h"
-#include "temps.h"
-#include "iniprint.h"
-#include "logic.h"
-
-c   Arguments:
-c   ----------
-
-      real time_0
-
-c   variables dynamiques
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
-      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
-      REAL ps(ip1jmp1)                       ! pression  au sol
-      REAL masse(ip1jmp1,llm)                ! masse d'air
-      REAL phis(ip1jmp1)                     ! geopotentiel au sol
-
-c   Local:
-c   ------
-
-      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
-      real tetajl(jjp1,llm)
-      INTEGER i,j,l,lsup,ij
-
-      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
-      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
-      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
-      LOGICAL ok_pv                ! Polar Vortex
-      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
-      
-      real zz,ran1
-      integer idum
-
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
-
-c-----------------------------------------------------------------------
-! 1. Initializations for Earth-like case
-! --------------------------------------
-c
-        ! initialize planet radius, rotation rate,...
-        call conf_planete
-
-        time_0=0.
-        day_ref=1
-        annee_ref=0
-
-        im         = iim
-        jm         = jjm
-        day_ini    = 1
-        dtvr    = daysec/REAL(day_step)
-        zdtvr=dtvr
-        etot0      = 0.
-        ptot0      = 0.
-        ztot0      = 0.
-        stot0      = 0.
-        ang0       = 0.
-
-        if (llm.eq.1) then
-          ! specific initializations for the shallow water case
-          kappa=1
-        endif
-        
-        CALL iniconst
-        CALL inigeom
-        CALL inifilr
-
-        if (llm.eq.1) then
-          ! initialize fields for the shallow water case, if required
-          if (.not.read_start) then
-            phis(:)=0.
-            q(:,:,:)=0
-            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
-          endif
-        endif
-
-        if (iflag_phys.eq.2) then
-          ! initializations for the academic case
-          
-!         if (planet_type=="earth") then
-
-          ! 1. local parameters
-          ! by convention, winter is in the southern hemisphere
-          ! Geostrophic wind or no wind?
-          ok_geost=.TRUE.
-          CALL getin('ok_geost',ok_geost)
-          ! Constants for Newtonian relaxation and friction
-          k_f=1.                !friction 
-          CALL getin('k_j',k_f)
-          k_f=1./(daysec*k_f)
-          k_c_s=4.  !cooling surface
-          CALL getin('k_c_s',k_c_s)
-          k_c_s=1./(daysec*k_c_s)
-          k_c_a=40. !cooling free atm
-          CALL getin('k_c_a',k_c_a)
-          k_c_a=1./(daysec*k_c_a)
-          ! Constants for Teta equilibrium profile
-          teta0=315.     ! mean Teta (S.H. 315K)
-          CALL getin('teta0',teta0)
-          ttp=200.       ! Tropopause temperature (S.H. 200K)
-          CALL getin('ttp',ttp)
-          eps=0.         ! Deviation to N-S symmetry(~0-20K)
-          CALL getin('eps',eps)
-          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
-          CALL getin('delt_y',delt_y)
-          delt_z=10.     ! Vertical Gradient (S.H. 10K)
-          CALL getin('delt_z',delt_z)
-          ! Polar vortex
-          ok_pv=.false.
-          CALL getin('ok_pv',ok_pv)
-          phi_pv=-50.            ! Latitude of edge of vortex
-          CALL getin('phi_pv',phi_pv)
-          phi_pv=phi_pv*pi/180.
-          dphi_pv=5.             ! Width of the edge
-          CALL getin('dphi_pv',dphi_pv)
-          dphi_pv=dphi_pv*pi/180.
-          gam_pv=4.              ! -dT/dz vortex (in K/km)
-          CALL getin('gam_pv',gam_pv)
-          
-          ! 2. Initialize fields towards which to relax
-          ! Friction
-          knewt_g=k_c_a
-          DO l=1,llm
-           zsig=presnivs(l)/preff
-           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
-           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
-          ENDDO
-          DO j=1,jjp1
-            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
-          ENDDO
-          
-          ! Potential temperature 
-          DO l=1,llm
-           zsig=presnivs(l)/preff
-           tetastrat=ttp*zsig**(-kappa)
-           tetapv=tetastrat
-           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
-               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
-           ENDIF
-           DO j=1,jjp1
-             ! Troposphere
-             ddsin=sin(rlatu(j))
-             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
-     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
-             ! Profil stratospherique isotherme (+vortex)
-             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
-             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
-             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
-           ENDDO
-          ENDDO ! of DO l=1,llm
- 
-!          CALL writefield('theta_eq',tetajl)
-
-          do l=1,llm
-            do j=1,jjp1
-              do i=1,iip1
-                 ij=(j-1)*iip1+i
-                 tetarappel(ij,l)=tetajl(j,l)
-              enddo
-            enddo
-          enddo
-
-
-!         else
-!          write(lunout,*)"iniacademic: planet types other than earth",
-!     &                   " not implemented (yet)."
-!          stop
-!         endif ! of if (planet_type=="earth")
-
-          ! 3. Initialize fields (if necessary)
-          IF (.NOT. read_start) THEN
-            ! surface pressure
-            ps(:)=preff
-            ! ground geopotential
-            phis(:)=0.
-            
-            CALL pression ( ip1jmp1, ap, bp, ps, p       )
-            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-            CALL massdair(p,masse)
-
-            ! bulk initialization of temperature
-            teta(:,:)=tetarappel(:,:)
-            
-            ! geopotential
-            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-            
-            ! winds
-            if (ok_geost) then
-              call ugeostr(phi,ucov)
-            else
-              ucov(:,:)=0.
-            endif
-            vcov(:,:)=0.
-            
-            ! bulk initialization of tracers
-            if (planet_type=="earth") then
-              ! Earth: first two tracers will be water
-              do i=1,nqtot
-                if (i.eq.1) q(:,:,i)=1.e-10
-                if (i.eq.2) q(:,:,i)=1.e-15
-                if (i.gt.2) q(:,:,i)=0.
-              enddo
-            else
-              q(:,:,:)=0
-            endif ! of if (planet_type=="earth")
-
-            ! add random perturbation to temperature
-            idum  = -1
-            zz = ran1(idum)
-            idum  = 0
-            do l=1,llm
-              do ij=iip2,ip1jm
-                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
-              enddo
-            enddo
-
-            ! maintain periodicity in longitude
-            do l=1,llm
-              do ij=1,ip1jmp1,iip1
-                teta(ij+iim,l)=teta(ij,l)
-              enddo
-            enddo
-
-          ENDIF ! of IF (.NOT. read_start)
-        endif ! of if (iflag_phys.eq.2)
-        
-      END
-c-----------------------------------------------------------------------
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90	(revision 1472)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90	(revision 1472)
@@ -0,0 +1,252 @@
+!
+! $Id$
+!
+SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+
+  USE filtreg_mod
+  USE infotrac, ONLY : nqtot
+  USE control_mod, ONLY: day_step,planet_type
+#ifdef CPP_IOIPSL
+  USE IOIPSL
+#else
+  ! if not using IOIPSL, we still need to use (a local version of) getin
+  USE ioipsl_getincom
+#endif
+  USE Write_Field
+
+  !   Author:    Frederic Hourdin      original: 15/01/93
+
+  IMPLICIT NONE
+
+  !   Declararations:
+  !   ---------------
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "comgeom.h"
+  include "academic.h"
+  include "ener.h"
+  include "temps.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  !   Arguments:
+  !   ----------
+
+  real time_0
+
+  !   variables dynamiques
+  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
+  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
+  REAL ps(ip1jmp1)                       ! pression  au sol
+  REAL masse(ip1jmp1,llm)                ! masse d'air
+  REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+  !   Local:
+  !   ------
+
+  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+  REAL pks(ip1jmp1)                      ! exner au  sol
+  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
+  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
+  REAL phi(ip1jmp1,llm)                  ! geopotentiel
+  REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
+  real tetajl(jjp1,llm)
+  INTEGER i,j,l,lsup,ij
+
+  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+  LOGICAL ok_pv                ! Polar Vortex
+  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+
+  real zz,ran1
+  integer idum
+
+  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+
+  !-----------------------------------------------------------------------
+  ! 1. Initializations for Earth-like case
+  ! --------------------------------------
+  !
+  ! initialize planet radius, rotation rate,...
+  call conf_planete
+
+  time_0=0.
+  day_ref=1
+  annee_ref=0
+
+  im         = iim
+  jm         = jjm
+  day_ini    = 1
+  dtvr    = daysec/REAL(day_step)
+  zdtvr=dtvr
+  etot0      = 0.
+  ptot0      = 0.
+  ztot0      = 0.
+  stot0      = 0.
+  ang0       = 0.
+
+  if (llm == 1) then
+     ! specific initializations for the shallow water case
+     kappa=1
+  endif
+
+  CALL iniconst
+  CALL inigeom
+  CALL inifilr
+
+  if (llm == 1) then
+     ! initialize fields for the shallow water case, if required
+     if (.not.read_start) then
+        phis(:)=0.
+        q(:,:,:)=0
+        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+     endif
+  endif
+
+  academic_case: if (iflag_phys == 2) then
+     ! initializations
+
+     ! 1. local parameters
+     ! by convention, winter is in the southern hemisphere
+     ! Geostrophic wind or no wind?
+     ok_geost=.TRUE.
+     CALL getin('ok_geost',ok_geost)
+     ! Constants for Newtonian relaxation and friction
+     k_f=1.                !friction 
+     CALL getin('k_j',k_f)
+     k_f=1./(daysec*k_f)
+     k_c_s=4.  !cooling surface
+     CALL getin('k_c_s',k_c_s)
+     k_c_s=1./(daysec*k_c_s)
+     k_c_a=40. !cooling free atm
+     CALL getin('k_c_a',k_c_a)
+     k_c_a=1./(daysec*k_c_a)
+     ! Constants for Teta equilibrium profile
+     teta0=315.     ! mean Teta (S.H. 315K)
+     CALL getin('teta0',teta0)
+     ttp=200.       ! Tropopause temperature (S.H. 200K)
+     CALL getin('ttp',ttp)
+     eps=0.         ! Deviation to N-S symmetry(~0-20K)
+     CALL getin('eps',eps)
+     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+     CALL getin('delt_y',delt_y)
+     delt_z=10.     ! Vertical Gradient (S.H. 10K)
+     CALL getin('delt_z',delt_z)
+     ! Polar vortex
+     ok_pv=.false.
+     CALL getin('ok_pv',ok_pv)
+     phi_pv=-50.            ! Latitude of edge of vortex
+     CALL getin('phi_pv',phi_pv)
+     phi_pv=phi_pv*pi/180.
+     dphi_pv=5.             ! Width of the edge
+     CALL getin('dphi_pv',dphi_pv)
+     dphi_pv=dphi_pv*pi/180.
+     gam_pv=4.              ! -dT/dz vortex (in K/km)
+     CALL getin('gam_pv',gam_pv)
+
+     ! 2. Initialize fields towards which to relax
+     ! Friction
+     knewt_g=k_c_a
+     DO l=1,llm
+        zsig=presnivs(l)/preff
+        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+     ENDDO
+     DO j=1,jjp1
+        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+     ENDDO
+
+     ! Potential temperature 
+     DO l=1,llm
+        zsig=presnivs(l)/preff
+        tetastrat=ttp*zsig**(-kappa)
+        tetapv=tetastrat
+        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+        ENDIF
+        DO j=1,jjp1
+           ! Troposphere
+           ddsin=sin(rlatu(j))
+           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
+                -delt_z*(1.-ddsin*ddsin)*log(zsig)
+           ! Profil stratospherique isotherme (+vortex)
+           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+           tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+        ENDDO
+     ENDDO
+
+     !          CALL writefield('theta_eq',tetajl)
+
+     do l=1,llm
+        do j=1,jjp1
+           do i=1,iip1
+              ij=(j-1)*iip1+i
+              tetarappel(ij,l)=tetajl(j,l)
+           enddo
+        enddo
+     enddo
+
+     ! 3. Initialize fields (if necessary)
+     IF (.NOT. read_start) THEN
+        ! surface pressure
+        ps(:)=preff
+        ! ground geopotential
+        phis(:)=0.
+
+        CALL pression ( ip1jmp1, ap, bp, ps, p       )
+        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL massdair(p,masse)
+
+        ! bulk initialization of temperature
+        teta(:,:)=tetarappel(:,:)
+
+        ! geopotential
+        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+
+        ! winds
+        if (ok_geost) then
+           call ugeostr(phi,ucov)
+        else
+           ucov(:,:)=0.
+        endif
+        vcov(:,:)=0.
+
+        ! bulk initialization of tracers
+        if (planet_type=="earth") then
+           ! Earth: first two tracers will be water
+           do i=1,nqtot
+              if (i == 1) q(:,:,i)=1.e-10
+              if (i == 2) q(:,:,i)=1.e-15
+              if (i.gt.2) q(:,:,i)=0.
+           enddo
+        else
+           q(:,:,:)=0
+        endif ! of if (planet_type=="earth")
+
+        ! add random perturbation to temperature
+        idum  = -1
+        zz = ran1(idum)
+        idum  = 0
+        do l=1,llm
+           do ij=iip2,ip1jm
+              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+           enddo
+        enddo
+
+        ! maintain periodicity in longitude
+        do l=1,llm
+           do ij=1,ip1jmp1,iip1
+              teta(ij+iim,l)=teta(ij,l)
+           enddo
+        enddo
+
+     ENDIF ! of IF (.NOT. read_start)
+  endif academic_case
+
+END SUBROUTINE iniacademic
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F	(revision 1470)
+++ 	(revision )
@@ -1,612 +1,0 @@
-!
-! $Id$
-!
-c
-      SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
-     s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
-     s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
-     s                   frac_impa, frac_nucl,
-     s                   prfl, psfl, rhcl, zqta, fraca,
-     s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
-
-c
-      USE dimphy
-      IMPLICIT none
-c======================================================================
-c Auteur(s): Z.X. Li (LMD/CNRS)
-c Date: le 20 mars 1995
-c Objet: condensation et precipitation stratiforme.
-c        schema de nuage
-c======================================================================
-c======================================================================
-cym#include "dimensions.h"
-cym#include "dimphy.h"
-#include "YOMCST.h"
-#include "tracstoke.h"
-#include "fisrtilp.h"
-c
-c Arguments:
-c
-      REAL dtime ! intervalle du temps (s)
-      REAL paprs(klon,klev+1) ! pression a inter-couche
-      REAL pplay(klon,klev) ! pression au milieu de couche
-      REAL t(klon,klev) ! temperature (K)
-      REAL q(klon,klev) ! humidite specifique (kg/kg)
-      REAL d_t(klon,klev) ! incrementation de la temperature (K)
-      REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
-      REAL d_ql(klon,klev) ! incrementation de l'eau liquide
-      REAL rneb(klon,klev) ! fraction nuageuse
-      REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
-      REAL rhcl(klon,klev) ! humidite relative en ciel clair
-      REAL rain(klon) ! pluies (mm/s)
-      REAL snow(klon) ! neige (mm/s)
-      REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
-      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
-      REAL ztv(klon,klev)
-      REAL zqta(klon,klev),fraca(klon,klev) 
-      REAL sigma1(klon,klev),sigma2(klon,klev)
-      REAL qltot(klon,klev),ctot(klon,klev)
-      REAL zpspsk(klon,klev),ztla(klon,klev)
-      REAL zthl(klon,klev)
-
-      logical lognormale(klon)
-
-cAA
-c Coeffients de fraction lessivee : pour OFF-LINE
-c
-      REAL pfrac_nucl(klon,klev)
-      REAL pfrac_1nucl(klon,klev)
-      REAL pfrac_impa(klon,klev)
-c
-c Fraction d'aerosols lessivee par impaction et par nucleation
-c POur ON-LINE
-c
-      REAL frac_impa(klon,klev)
-      REAL frac_nucl(klon,klev)
-      real zct      ,zcl
-cAA
-c
-c Options du programme:
-c
-      REAL seuil_neb ! un nuage existe vraiment au-dela
-      PARAMETER (seuil_neb=0.001)
-
-      INTEGER ninter ! sous-intervals pour la precipitation
-      INTEGER ncoreczq  
-      INTEGER iflag_cldcon
-      PARAMETER (ninter=5)
-      LOGICAL evap_prec ! evaporation de la pluie
-      PARAMETER (evap_prec=.TRUE.)
-      REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
-      logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
-
-      real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
-      real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
-      real erf   
-      REAL qcloud(klon)
-c
-      LOGICAL cpartiel ! condensation partielle
-      PARAMETER (cpartiel=.TRUE.)
-      REAL t_coup
-      PARAMETER (t_coup=234.0)
-c
-c Variables locales:
-c
-      INTEGER i, k, n, kk
-      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
-      REAL zrfl(klon), zrfln(klon), zqev, zqevt
-      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
-      REAL ztglace, zt(klon)
-      INTEGER nexpo ! exponentiel pour glace/eau
-      REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
-      REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
-c
-      LOGICAL appel1er
-      SAVE appel1er
-c$OMP THREADPRIVATE(appel1er)
-c
-c---------------------------------------------------------------
-c
-cAA Variables traceurs:
-cAA  Provisoire !!! Parametres alpha du lessivage
-cAA  A priori on a 4 scavenging # possibles
-c
-      REAL a_tr_sca(4)
-      save a_tr_sca
-c$OMP THREADPRIVATE(a_tr_sca)
-c
-c Variables intermediaires
-c
-      REAL zalpha_tr
-      REAL zfrac_lessi
-      REAL zprec_cond(klon)
-cAA
-      REAL zmair, zcpair, zcpeau
-C     Pour la conversion eau-neige
-      REAL zlh_solid(klon), zm_solid
-cIM 
-cym      INTEGER klevm1
-c---------------------------------------------------------------
-c
-c Fonctions en ligne:
-c
-      REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
-      REAL zzz
-#include "YOETHF.h"
-#include "FCTTRE.h"
-      fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
-      fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
-c
-      DATA appel1er /.TRUE./
-cym
-      zdelq=0.0
-      
-      print*,'NUAGES4 A. JAM'
-      IF (appel1er) THEN
-c
-         PRINT*, 'fisrtilp, ninter:', ninter
-         PRINT*, 'fisrtilp, evap_prec:', evap_prec
-         PRINT*, 'fisrtilp, cpartiel:', cpartiel
-         IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
-          PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
-          PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
-c         CALL abort
-         ENDIF
-         appel1er = .FALSE.
-c
-cAA initialiation provisoire
-       a_tr_sca(1) = -0.5
-       a_tr_sca(2) = -0.5
-       a_tr_sca(3) = -0.5
-       a_tr_sca(4) = -0.5
-c
-cAA Initialisation a 1 des coefs des fractions lessivees 
-c
-!cdir collapse
-      DO k = 1, klev
-       DO i = 1, klon
-          pfrac_nucl(i,k)=1.
-          pfrac_1nucl(i,k)=1.
-          pfrac_impa(i,k)=1.
-       ENDDO 
-      ENDDO 
-
-      ENDIF          !  test sur appel1er
-c
-cMAf Initialisation a 0 de zoliq
-c      DO i = 1, klon
-c         zoliq(i)=0.
-c      ENDDO 
-c Determiner les nuages froids par leur temperature
-c  nexpo regle la raideur de la transition eau liquide / eau glace.
-c
-      ztglace = RTT - 15.0
-      nexpo = 6
-ccc      nexpo = 1
-c
-c Initialiser les sorties:
-c
-!cdir collapse
-      DO k = 1, klev+1
-      DO i = 1, klon
-         prfl(i,k) = 0.0
-         psfl(i,k) = 0.0
-      ENDDO
-      ENDDO
-
-!cdir collapse
-      DO k = 1, klev
-      DO i = 1, klon
-         d_t(i,k) = 0.0
-         d_q(i,k) = 0.0
-         d_ql(i,k) = 0.0
-         rneb(i,k) = 0.0
-         radliq(i,k) = 0.0
-         frac_nucl(i,k) = 1. 
-         frac_impa(i,k) = 1. 
-      ENDDO
-      ENDDO
-      DO i = 1, klon
-         rain(i) = 0.0
-         snow(i) = 0.0
-         zoliq(i)=0.
-c     ENDDO
-c
-c Initialiser le flux de precipitation a zero
-c
-c     DO i = 1, klon
-         zrfl(i) = 0.0
-         zneb(i) = seuil_neb
-      ENDDO
-c
-c
-cAA Pour plus de securite 
-
-      zalpha_tr   = 0.
-      zfrac_lessi = 0.
-
-cAA----------------------------------------------------------
-c
-      ncoreczq=0
-c Boucle verticale (du haut vers le bas)
-c
-cIM : klevm1
-cym      klevm1=klev-1
-      DO 9999 k = klev, 1, -1
-c
-cAA----------------------------------------------------------
-c
-      DO i = 1, klon
-         zt(i)=t(i,k)
-         zq(i)=q(i,k)
-      ENDDO
-c
-c Calculer la varition de temp. de l'air du a la chaleur sensible
-C transporter par la pluie.
-C Il resterait a rajouter cet effet de la chaleur sensible sur les
-C flux de surface, du a la diff. de temp. entre le 1er niveau et la
-C surface.
-C
-      IF(k.LE.klevm1) THEN         
-         DO i = 1, klon
-cIM
-            zmair=(paprs(i,k)-paprs(i,k+1))/RG
-            zcpair=RCPD*(1.0+RVTMP2*zq(i))
-            zcpeau=RCPD*RVTMP2
-            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
-     $           + zmair*zcpair*zt(i) )
-     $           / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
-C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
-         ENDDO
-      ENDIF
-c
-c
-c Calculer l'evaporation de la precipitation
-c
-
-
-      IF (evap_prec) THEN
-      DO i = 1, klon
-      IF (zrfl(i) .GT.0.) THEN
-         IF (thermcep) THEN
-           zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
-           zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-           zqs(i)=MIN(0.5,zqs(i))
-           zcor=1./(1.-RETV*zqs(i))
-           zqs(i)=zqs(i)*zcor
-         ELSE
-           IF (zt(i) .LT. t_coup) THEN
-              zqs(i) = qsats(zt(i)) / pplay(i,k)
-           ELSE
-              zqs(i) = qsatl(zt(i)) / pplay(i,k)
-           ENDIF
-         ENDIF
-         zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
-         zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
-     .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
-         zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
-     .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
-         zqev = MIN (zqev, zqevt)
-         zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
-     .                            /RG/dtime
-
-c pour la glace, on rï¿½vapore toute la prï¿½ip dans la couche du dessous
-c la glace venant de la couche du dessus est simplement dans la couche
-c du dessous.
-
-         IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
-
-         zq(i) = zq(i) - (zrfln(i)-zrfl(i))
-     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-         zt(i) = zt(i) + (zrfln(i)-zrfl(i))
-     .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
-     .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-         zrfl(i) = zrfln(i)
-      ENDIF
-      ENDDO
-      ENDIF
-c
-c Calculer Qs et L/Cp*dQs/dT:
-c
-      IF (thermcep) THEN
-         DO i = 1, klon
-           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
-           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
-           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
-           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
-           zqs(i) = MIN(0.5,zqs(i))
-           zcor = 1./(1.-RETV*zqs(i))
-           zqs(i) = zqs(i)*zcor
-           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
-         ENDDO
-      ELSE
-         DO i = 1, klon
-            IF (zt(i).LT.t_coup) THEN
-               zqs(i) = qsats(zt(i))/pplay(i,k)
-               zdqs(i) = dqsats(zt(i),zqs(i))
-            ELSE
-               zqs(i) = qsatl(zt(i))/pplay(i,k)
-               zdqs(i) = dqsatl(zt(i),zqs(i))
-            ENDIF
-         ENDDO
-      ENDIF
-c
-c Determiner la condensation partielle et calculer la quantite
-c de l'eau condensee:
-c
-
-      IF (cpartiel) THEN
-
-c        print*,'Dans partiel k=',k
-c
-c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
-c   nuageuse a partir des PDF de Sandrine Bony.
-c   rneb  : fraction nuageuse
-c   zqn   : eau totale dans le nuage
-c   zcond : eau condensee moyenne dans la maille.
-c           on prend en compte le rï¿½hauffement qui diminue la partie condensee
-c
-c   Version avec les raqts
-
-         if (iflag_pdf.eq.0) then
-
-           do i=1,klon
-            zdelq = min(ratqs(i,k),0.99) * zq(i)
-            rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
-            zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
-           enddo
-
-         else
-c
-c   Version avec les nouvelles PDFs.
-           do i=1,klon
-              if(zq(i).lt.1.e-15) then
-                ncoreczq=ncoreczq+1
-                zq(i)=1.e-15
-              endif
-           enddo 
-
-              if (iflag_cldcon>=5) then
-
-                 call cloudth(klon,klev,k,ztv,
-     .           zq,zqta,fraca,
-     .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
-     .           ratqs,zqs,t)
-
-                 do i=1,klon
-                 rneb(i,k)=ctot(i,k)
-                 zqn(i)=qcloud(i)
-                 enddo
-
-              endif
-
-! Pour iflag_cldcon<=4, on prend toujours la lognormale
-! Dans le cas iflag_cldcon=5, on prend systÃ©matiquement la bi-gaussienne
-! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques
-
-            lognormale(:)=
-     .      iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10)
-            do i=1,klon
-            if (lognormale(i)) then
-            zpdf_sig(i)=ratqs(i,k)*zq(i)
-            zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
-            zpdf_delta(i)=log(zq(i)/zqs(i))
-            zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
-            zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
-            zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
-            zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
-            zpdf_e1(i)=1.-erf(zpdf_e1(i))
-            zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
-            zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
-            zpdf_e2(i)=1.-erf(zpdf_e2(i))
-            endif
-            enddo
-
-            do i=1,klon
-            if (lognormale(i)) then
-            if (zpdf_e1(i).lt.1.e-10) then
-               rneb(i,k)=0.
-               zqn(i)=zqs(i)
-            else
-               rneb(i,k)=0.5*zpdf_e1(i)
-               zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
-            endif
-            endif
-            
-           enddo 
-
-
-        endif ! iflag_pdf
-
-        DO i=1,klon
-           IF (rneb(i,k) .LE. 0.0) THEN
-              zqn(i) = 0.0
-              rneb(i,k) = 0.0
-              zcond(i) = 0.0
-              rhcl(i,k)=zq(i)/zqs(i)
-           ELSE IF (rneb(i,k) .GE. 1.0) THEN
-              zqn(i) = zq(i)
-              rneb(i,k) = 1.0                  
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
-              rhcl(i,k)=1.0
-           ELSE
-              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-           ENDIF
-        ENDDO
-!         do i=1,klon
-!            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
-!            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
-!            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
-!c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
-!c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
-!c  la convection.
-!c  ATTENTION !!! Il va falloir verifier tout ca.
-!            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
-!c           print*,'ZDQS ',zdqs(i)
-!c--Olivier
-!            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
-!            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
-!            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
-!c--fin
-!           ENDDO
-      ELSE
-         DO i = 1, klon
-            IF (zq(i).GT.zqs(i)) THEN
-               rneb(i,k) = 1.0
-            ELSE
-               rneb(i,k) = 0.0
-            ENDIF
-            zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
-         ENDDO
-      ENDIF
-c
-      DO i = 1, klon
-         zq(i) = zq(i) - zcond(i)
-c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
-         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
-      ENDDO
-c
-c Partager l'eau condensee en precipitation et eau liquide nuageuse
-c
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         zoliq(i) = zcond(i)
-         zrho(i) = pplay(i,k) / zt(i) / RD
-         zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
-         zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
-         zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
-         zfice(i) = zfice(i)**nexpo
-         zneb(i) = MAX(rneb(i,k), seuil_neb)
-         radliq(i,k) = zoliq(i)/REAL(ninter+1)
-      ENDIF
-      ENDDO
-c
-      DO n = 1, ninter
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-
-         IF (zneb(i).EQ.seuil_neb) THEN
-             ztot = 0.0
-         ELSE
-c  quantite d'eau a eliminer: zchau
-c  meme chose pour la glace: zfroi
-             if (ptconv(i,k)) then
-                zcl   =cld_lc_con
-                zct   =1./cld_tau_con
-                zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
-     .              *fallvc(zrhol(i)) * zfice(i)
-             else
-                zcl   =cld_lc_lsc
-                zct   =1./cld_tau_lsc
-                zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
-     .              *fallvs(zrhol(i)) * zfice(i)
-             endif
-             zchau    = zct   *dtime/REAL(ninter) * zoliq(i)
-     .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
-             ztot    = zchau    + zfroi
-             ztot    = MAX(ztot   ,0.0)
-         ENDIF
-         ztot    = MIN(ztot,zoliq(i))
-         zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
-         radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
-      ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i = 1, klon
-      IF (rneb(i,k).GT.0.0) THEN
-         d_ql(i,k) = zoliq(i)
-         zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
-     .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
-      ENDIF
-      IF (zt(i).LT.RTT) THEN
-        psfl(i,k)=zrfl(i)
-      ELSE
-        prfl(i,k)=zrfl(i)
-      ENDIF
-      ENDDO
-c
-c Calculer les tendances de q et de t:
-c
-      DO i = 1, klon
-         d_q(i,k) = zq(i) - q(i,k)
-         d_t(i,k) = zt(i) - t(i,k)
-      ENDDO
-c
-cAA--------------- Calcul du lessivage stratiforme  -------------
-
-      DO i = 1,klon
-c
-         zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
-     .                * (paprs(i,k)-paprs(i,k+1))/RG
-         IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-cAA lessivage nucleation LMD5 dans la couche elle-meme
-            if (t(i,k) .GE. ztglace) THEN
-               zalpha_tr = a_tr_sca(3)
-            else
-               zalpha_tr = a_tr_sca(4)
-            endif
-            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
-c
-c nucleation avec un facteur -1 au lieu de -0.5
-            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
-            pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
-         ENDIF
-c
-      ENDDO      ! boucle sur i
-c
-cAA Lessivage par impaction dans les couches en-dessous
-      DO kk = k-1, 1, -1
-        DO i = 1, klon
-          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
-            if (t(i,kk) .GE. ztglace) THEN
-              zalpha_tr = a_tr_sca(1)
-            else
-              zalpha_tr = a_tr_sca(2)
-            endif
-            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
-            pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
-            frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
-          ENDIF
-        ENDDO
-      ENDDO
-c
-cAA----------------------------------------------------------
-c                     FIN DE BOUCLE SUR K   
- 9999 CONTINUE
-c
-cAA-----------------------------------------------------------
-c
-c Pluie ou neige au sol selon la temperature de la 1ere couche
-c
-      DO i = 1, klon
-      IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
-         snow(i) = zrfl(i)
-         zlh_solid(i) = RLSTT-RLVTT
-      ELSE
-         rain(i) = zrfl(i)
-         zlh_solid(i) = 0.
-      ENDIF
-      ENDDO
-C
-C For energy conservation : when snow is present, the solification
-c latent heat is considered.
-      DO k = 1, klev
-        DO i = 1, klon
-          zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
-          zmair=(paprs(i,k)-paprs(i,k+1))/RG
-          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
-          d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
-        END DO 
-      END DO
-c
-
-      if (ncoreczq>0) then
-         print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
-      endif
-      RETURN
-      END
Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90	(revision 1472)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90	(revision 1472)
@@ -0,0 +1,619 @@
+!
+! $Id$
+!
+!
+SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
+     d_t, d_q, d_ql, rneb, radliq, rain, snow, &
+     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
+     frac_impa, frac_nucl, &
+     prfl, psfl, rhcl, zqta, fraca, &
+     ztv, zpspsk, ztla, zthl, iflag_cldcon)
+
+  !
+  USE dimphy
+  IMPLICIT none
+  !======================================================================
+  ! Auteur(s): Z.X. Li (LMD/CNRS)
+  ! Date: le 20 mars 1995
+  ! Objet: condensation et precipitation stratiforme.
+  !        schema de nuage
+  !======================================================================
+  !======================================================================
+  !ym include "dimensions.h"
+  !ym include "dimphy.h"
+  include "YOMCST.h"
+  include "tracstoke.h"
+  include "fisrtilp.h"
+  !
+  ! Arguments:
+  !
+  REAL dtime ! intervalle du temps (s)
+  REAL paprs(klon,klev+1) ! pression a inter-couche
+  REAL pplay(klon,klev) ! pression au milieu de couche
+  REAL t(klon,klev) ! temperature (K)
+  REAL q(klon,klev) ! humidite specifique (kg/kg)
+  REAL d_t(klon,klev) ! incrementation de la temperature (K)
+  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
+  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
+  REAL rneb(klon,klev) ! fraction nuageuse
+  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
+  REAL rhcl(klon,klev) ! humidite relative en ciel clair
+  REAL rain(klon) ! pluies (mm/s)
+  REAL snow(klon) ! neige (mm/s)
+  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
+  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
+  REAL ztv(klon,klev)
+  REAL zqta(klon,klev),fraca(klon,klev) 
+  REAL sigma1(klon,klev),sigma2(klon,klev)
+  REAL qltot(klon,klev),ctot(klon,klev)
+  REAL zpspsk(klon,klev),ztla(klon,klev)
+  REAL zthl(klon,klev)
+
+  logical lognormale(klon)
+
+  !AA
+  ! Coeffients de fraction lessivee : pour OFF-LINE
+  !
+  REAL pfrac_nucl(klon,klev)
+  REAL pfrac_1nucl(klon,klev)
+  REAL pfrac_impa(klon,klev)
+  !
+  ! Fraction d'aerosols lessivee par impaction et par nucleation
+  ! POur ON-LINE
+  !
+  REAL frac_impa(klon,klev)
+  REAL frac_nucl(klon,klev)
+  real zct      ,zcl
+  !AA
+  !
+  ! Options du programme:
+  !
+  REAL seuil_neb ! un nuage existe vraiment au-dela
+  PARAMETER (seuil_neb=0.001)
+
+  INTEGER ninter ! sous-intervals pour la precipitation
+  INTEGER ncoreczq  
+  INTEGER iflag_cldcon
+  PARAMETER (ninter=5)
+  LOGICAL evap_prec ! evaporation de la pluie
+  PARAMETER (evap_prec=.TRUE.)
+  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
+  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
+
+  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
+  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
+  real erf   
+  REAL qcloud(klon)
+  !
+  LOGICAL cpartiel ! condensation partielle
+  PARAMETER (cpartiel=.TRUE.)
+  REAL t_coup
+  PARAMETER (t_coup=234.0)
+  !
+  ! Variables locales:
+  !
+  INTEGER i, k, n, kk
+  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
+  REAL zrfl(klon), zrfln(klon), zqev, zqevt
+  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
+  REAL ztglace, zt(klon)
+  INTEGER nexpo ! exponentiel pour glace/eau
+  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
+  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
+  !
+  LOGICAL appel1er
+  SAVE appel1er
+  !$OMP THREADPRIVATE(appel1er)
+  !
+  !---------------------------------------------------------------
+  !
+  !AA Variables traceurs:
+  !AA  Provisoire !!! Parametres alpha du lessivage
+  !AA  A priori on a 4 scavenging # possibles
+  !
+  REAL a_tr_sca(4)
+  save a_tr_sca
+  !$OMP THREADPRIVATE(a_tr_sca)
+  !
+  ! Variables intermediaires
+  !
+  REAL zalpha_tr
+  REAL zfrac_lessi
+  REAL zprec_cond(klon)
+  !AA
+  REAL zmair, zcpair, zcpeau
+  !     Pour la conversion eau-neige
+  REAL zlh_solid(klon), zm_solid
+  !IM 
+  !ym      INTEGER klevm1
+  !---------------------------------------------------------------
+  !
+  ! Fonctions en ligne:
+  !
+  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
+  REAL zzz
+  include "YOETHF.h"
+  include "FCTTRE.h"
+  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
+  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
+  !
+  DATA appel1er /.TRUE./
+  !ym
+  zdelq=0.0
+
+  print*,'NUAGES4 A. JAM'
+  IF (appel1er) THEN
+     !
+     PRINT*, 'fisrtilp, ninter:', ninter
+     PRINT*, 'fisrtilp, evap_prec:', evap_prec
+     PRINT*, 'fisrtilp, cpartiel:', cpartiel
+     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
+        PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
+        PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
+        !         CALL abort
+     ENDIF
+     appel1er = .FALSE.
+     !
+     !AA initialiation provisoire
+     a_tr_sca(1) = -0.5
+     a_tr_sca(2) = -0.5
+     a_tr_sca(3) = -0.5
+     a_tr_sca(4) = -0.5
+     !
+     !AA Initialisation a 1 des coefs des fractions lessivees 
+     !
+     !cdir collapse
+     DO k = 1, klev
+        DO i = 1, klon
+           pfrac_nucl(i,k)=1.
+           pfrac_1nucl(i,k)=1.
+           pfrac_impa(i,k)=1.
+        ENDDO
+     ENDDO
+
+  ENDIF          !  test sur appel1er
+  !
+  !MAf Initialisation a 0 de zoliq
+  !      DO i = 1, klon
+  !         zoliq(i)=0.
+  !      ENDDO 
+  ! Determiner les nuages froids par leur temperature
+  !  nexpo regle la raideur de la transition eau liquide / eau glace.
+  !
+  ztglace = RTT - 15.0
+  nexpo = 6
+  !cc      nexpo = 1
+  !
+  ! Initialiser les sorties:
+  !
+  !cdir collapse
+  DO k = 1, klev+1
+     DO i = 1, klon
+        prfl(i,k) = 0.0
+        psfl(i,k) = 0.0
+     ENDDO
+  ENDDO
+
+  !cdir collapse
+  DO k = 1, klev
+     DO i = 1, klon
+        d_t(i,k) = 0.0
+        d_q(i,k) = 0.0
+        d_ql(i,k) = 0.0
+        rneb(i,k) = 0.0
+        radliq(i,k) = 0.0
+        frac_nucl(i,k) = 1. 
+        frac_impa(i,k) = 1. 
+     ENDDO
+  ENDDO
+  DO i = 1, klon
+     rain(i) = 0.0
+     snow(i) = 0.0
+     zoliq(i)=0.
+     !     ENDDO
+     !
+     ! Initialiser le flux de precipitation a zero
+     !
+     !     DO i = 1, klon
+     zrfl(i) = 0.0
+     zneb(i) = seuil_neb
+  ENDDO
+  !
+  !
+  !AA Pour plus de securite 
+
+  zalpha_tr   = 0.
+  zfrac_lessi = 0.
+
+  !AA----------------------------------------------------------
+  !
+  ncoreczq=0
+  ! Boucle verticale (du haut vers le bas)
+  !
+  !IM : klevm1
+  !ym      klevm1=klev-1
+  DO k = klev, 1, -1
+     !
+     !AA----------------------------------------------------------
+     !
+     DO i = 1, klon
+        zt(i)=t(i,k)
+        zq(i)=q(i,k)
+     ENDDO
+     !
+     ! Calculer la varition de temp. de l'air du a la chaleur sensible
+     ! transporter par la pluie.
+     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
+     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
+     ! surface.
+     !
+     IF(k.LE.klevm1) THEN         
+        DO i = 1, klon
+           !IM
+           zmair=(paprs(i,k)-paprs(i,k+1))/RG
+           zcpair=RCPD*(1.0+RVTMP2*zq(i))
+           zcpeau=RCPD*RVTMP2
+           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
+                + zmair*zcpair*zt(i) ) &
+                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
+           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
+        ENDDO
+     ENDIF
+     !
+     !
+     ! Calculer l'evaporation de la precipitation
+     !
+
+
+     IF (evap_prec) THEN
+        DO i = 1, klon
+           IF (zrfl(i) .GT.0.) THEN
+              IF (thermcep) THEN
+                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
+                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+                 zqs(i)=MIN(0.5,zqs(i))
+                 zcor=1./(1.-RETV*zqs(i))
+                 zqs(i)=zqs(i)*zcor
+              ELSE
+                 IF (zt(i) .LT. t_coup) THEN
+                    zqs(i) = qsats(zt(i)) / pplay(i,k)
+                 ELSE
+                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
+                 ENDIF
+              ENDIF
+              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
+              zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
+                   * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+              zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
+                   * RG*dtime/(paprs(i,k)-paprs(i,k+1))
+              zqev = MIN (zqev, zqevt)
+              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
+                   /RG/dtime
+
+              ! pour la glace, on ré-évapore toute la précip dans la
+              ! couche du dessous
+              ! la glace venant de la couche du dessus est simplement
+              ! dans la couche du dessous.
+
+              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
+
+              zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
+                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
+                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
+                   * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+              zrfl(i) = zrfln(i)
+           ENDIF
+        ENDDO
+     ENDIF
+     !
+     ! Calculer Qs et L/Cp*dQs/dT:
+     !
+     IF (thermcep) THEN
+        DO i = 1, klon
+           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
+           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
+           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
+           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
+           zqs(i) = MIN(0.5,zqs(i))
+           zcor = 1./(1.-RETV*zqs(i))
+           zqs(i) = zqs(i)*zcor
+           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
+        ENDDO
+     ELSE
+        DO i = 1, klon
+           IF (zt(i).LT.t_coup) THEN
+              zqs(i) = qsats(zt(i))/pplay(i,k)
+              zdqs(i) = dqsats(zt(i),zqs(i))
+           ELSE
+              zqs(i) = qsatl(zt(i))/pplay(i,k)
+              zdqs(i) = dqsatl(zt(i),zqs(i))
+           ENDIF
+        ENDDO
+     ENDIF
+     !
+     ! Determiner la condensation partielle et calculer la quantite
+     ! de l'eau condensee:
+     !
+
+     IF (cpartiel) THEN
+
+        !        print*,'Dans partiel k=',k
+        !
+        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
+        !   nuageuse a partir des PDF de Sandrine Bony.
+        !   rneb  : fraction nuageuse
+        !   zqn   : eau totale dans le nuage
+        !   zcond : eau condensee moyenne dans la maille.
+        !  on prend en compte le réchauffement qui diminue la partie
+        ! condensee
+        !
+        !   Version avec les raqts
+
+        if (iflag_pdf.eq.0) then
+
+           do i=1,klon
+              zdelq = min(ratqs(i,k),0.99) * zq(i)
+              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
+              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
+           enddo
+
+        else
+           !
+           !   Version avec les nouvelles PDFs.
+           do i=1,klon
+              if(zq(i).lt.1.e-15) then
+                 ncoreczq=ncoreczq+1
+                 zq(i)=1.e-15
+              endif
+           enddo
+
+           if (iflag_cldcon>=5) then
+
+              call cloudth(klon,klev,k,ztv, &
+                   zq,zqta,fraca, &
+                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
+                   ratqs,zqs,t)
+
+              do i=1,klon
+                 rneb(i,k)=ctot(i,k)
+                 zqn(i)=qcloud(i)
+              enddo
+
+           endif
+
+           if (iflag_cldcon <= 4) then
+              lognormale = .true.
+           elseif (iflag_cldcon == 6) then
+              ! lognormale en l'absence des thermiques
+              lognormale = fraca(:,k) < 1e-10
+           else
+              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
+              ! bi-gaussienne
+              lognormale = .false.
+           end if
+
+           do i=1,klon
+              if (lognormale(i)) then
+                 zpdf_sig(i)=ratqs(i,k)*zq(i)
+                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
+                 zpdf_delta(i)=log(zq(i)/zqs(i))
+                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
+                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
+                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
+                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
+                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
+                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
+                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
+                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
+              endif
+           enddo
+
+           do i=1,klon
+              if (lognormale(i)) then
+                 if (zpdf_e1(i).lt.1.e-10) then
+                    rneb(i,k)=0.
+                    zqn(i)=zqs(i)
+                 else
+                    rneb(i,k)=0.5*zpdf_e1(i)
+                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+                 endif
+              endif
+
+           enddo
+
+
+        endif ! iflag_pdf
+
+        DO i=1,klon
+           IF (rneb(i,k) .LE. 0.0) THEN
+              zqn(i) = 0.0
+              rneb(i,k) = 0.0
+              zcond(i) = 0.0
+              rhcl(i,k)=zq(i)/zqs(i)
+           ELSE IF (rneb(i,k) .GE. 1.0) THEN
+              zqn(i) = zq(i)
+              rneb(i,k) = 1.0                  
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))
+              rhcl(i,k)=1.0
+           ELSE
+              zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+           ENDIF
+        ENDDO
+        !         do i=1,klon
+        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
+        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
+        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
+        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
+        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
+        !c  la convection.
+        !c  ATTENTION !!! Il va falloir verifier tout ca.
+        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
+        !c           print*,'ZDQS ',zdqs(i)
+        !c--Olivier
+        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
+        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
+        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
+        !c--fin
+        !           ENDDO
+     ELSE
+        DO i = 1, klon
+           IF (zq(i).GT.zqs(i)) THEN
+              rneb(i,k) = 1.0
+           ELSE
+              rneb(i,k) = 0.0
+           ENDIF
+           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
+        ENDDO
+     ENDIF
+     !
+     DO i = 1, klon
+        zq(i) = zq(i) - zcond(i)
+        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
+        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
+     ENDDO
+     !
+     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
+     !
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           zoliq(i) = zcond(i)
+           zrho(i) = pplay(i,k) / zt(i) / RD
+           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
+           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
+           zfice(i) = zfice(i)**nexpo
+           zneb(i) = MAX(rneb(i,k), seuil_neb)
+           radliq(i,k) = zoliq(i)/REAL(ninter+1)
+        ENDIF
+     ENDDO
+     !
+     DO n = 1, ninter
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0) THEN
+              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+
+              IF (zneb(i).EQ.seuil_neb) THEN
+                 ztot = 0.0
+              ELSE
+                 !  quantite d'eau a eliminer: zchau
+                 !  meme chose pour la glace: zfroi
+                 if (ptconv(i,k)) then
+                    zcl   =cld_lc_con
+                    zct   =1./cld_tau_con
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvc(zrhol(i)) * zfice(i)
+                 else
+                    zcl   =cld_lc_lsc
+                    zct   =1./cld_tau_lsc
+                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
+                         *fallvs(zrhol(i)) * zfice(i)
+                 endif
+                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
+                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
+                 ztot    = zchau    + zfroi
+                 ztot    = MAX(ztot   ,0.0)
+              ENDIF
+              ztot    = MIN(ztot,zoliq(i))
+              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
+              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+           ENDIF
+        ENDDO
+     ENDDO
+     !
+     DO i = 1, klon
+        IF (rneb(i,k).GT.0.0) THEN
+           d_ql(i,k) = zoliq(i)
+           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
+                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+        ENDIF
+        IF (zt(i).LT.RTT) THEN
+           psfl(i,k)=zrfl(i)
+        ELSE
+           prfl(i,k)=zrfl(i)
+        ENDIF
+     ENDDO
+     !
+     ! Calculer les tendances de q et de t:
+     !
+     DO i = 1, klon
+        d_q(i,k) = zq(i) - q(i,k)
+        d_t(i,k) = zt(i) - t(i,k)
+     ENDDO
+     !
+     !AA--------------- Calcul du lessivage stratiforme  -------------
+
+     DO i = 1,klon
+        !
+        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
+             * (paprs(i,k)-paprs(i,k+1))/RG
+        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+           !AA lessivage nucleation LMD5 dans la couche elle-meme
+           if (t(i,k) .GE. ztglace) THEN
+              zalpha_tr = a_tr_sca(3)
+           else
+              zalpha_tr = a_tr_sca(4)
+           endif
+           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 
+           !
+           ! nucleation avec un facteur -1 au lieu de -0.5
+           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
+           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
+        ENDIF
+        !
+     ENDDO      ! boucle sur i
+     !
+     !AA Lessivage par impaction dans les couches en-dessous
+     DO kk = k-1, 1, -1
+        DO i = 1, klon
+           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
+              if (t(i,kk) .GE. ztglace) THEN
+                 zalpha_tr = a_tr_sca(1)
+              else
+                 zalpha_tr = a_tr_sca(2)
+              endif
+              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
+              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
+              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
+           ENDIF
+        ENDDO
+     ENDDO
+     !
+     !AA----------------------------------------------------------
+     !                     FIN DE BOUCLE SUR K   
+  end DO
+  !
+  !AA-----------------------------------------------------------
+  !
+  ! Pluie ou neige au sol selon la temperature de la 1ere couche
+  !
+  DO i = 1, klon
+     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
+        snow(i) = zrfl(i)
+        zlh_solid(i) = RLSTT-RLVTT
+     ELSE
+        rain(i) = zrfl(i)
+        zlh_solid(i) = 0.
+     ENDIF
+  ENDDO
+  !
+  ! For energy conservation : when snow is present, the solification
+  ! latent heat is considered.
+  DO k = 1, klev
+     DO i = 1, klon
+        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
+        zmair=(paprs(i,k)-paprs(i,k+1))/RG
+        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
+        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
+     END DO
+  END DO
+  !
+
+  if (ncoreczq>0) then
+     print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
+  endif
+
+END SUBROUTINE fisrtilp
