Index: LMDZ4/tags/V3_4/libf/phylmd/hgardfou.F
===================================================================
--- LMDZ4/tags/V3_4/libf/phylmd/hgardfou.F	(revision 944)
+++ LMDZ4/tags/V3_4/libf/phylmd/hgardfou.F	(revision 944)
@@ -0,0 +1,133 @@
+!
+! $Header$
+!
+      SUBROUTINE hgardfou (rlon,rlat,pctsrf,t,tsol,text)
+      use dimphy
+      IMPLICIT none
+c======================================================================
+c Verifier la temperature
+c======================================================================
+#include "dimensions.h"
+cym#include "dimphy.h"
+#include "YOMCST.h"
+#include "indicesol.h"
+      REAL rlon(klon), rlat(klon), pctsrf(klon, nbsrf)
+      REAL t(klon,klev), tsol(klon,nbsrf)
+      CHARACTER*(*) text
+C
+      INTEGER i, k, nsrf
+      REAL zt(klon)
+      INTEGER jadrs(klon), jbad
+      LOGICAL ok
+c
+      LOGICAL firstcall
+      SAVE firstcall
+      DATA firstcall /.TRUE./
+      IF (firstcall) THEN
+         PRINT*, 'hgardfou garantit la temperature dans [100,370] K'
+         firstcall = .FALSE.
+c
+         DO i = 1, klon
+          print*,'i=',i,'rlon=',rlon(i),'rlat=',rlat(i)
+         ENDDO
+c
+      ENDIF
+c
+      ok = .TRUE.
+      DO k = 1, klev
+         DO i = 1, klon
+            zt(i) = t(i,k)
+         ENDDO
+#ifdef CRAY
+         CALL WHENFGT(klon, zt, 1, 370.0, jadrs, jbad)
+#else
+         jbad = 0
+         DO i = 1, klon
+         IF (zt(i).GT.370.0) THEN
+            jbad = jbad + 1
+            jadrs(jbad) = i
+         ENDIF
+         ENDDO
+#endif
+         IF (jbad .GT. 0) THEN
+           ok = .FALSE.
+           DO i = 1, jbad
+             PRINT *,'i,k,temperature,lon,lat,pourc ter,oce,lic,sic =',
+     $       jadrs(i),k,zt(jadrs(i)),rlon(jadrs(i)),rlat(jadrs(i)),
+     $       (pctsrf(jadrs(i),nsrf),nsrf=1,nbsrf)
+           ENDDO
+         ENDIF
+#ifdef CRAY
+         CALL WHENFLT(klon, zt, 1, 100.0, jadrs, jbad)
+#else
+         jbad = 0
+         DO i = 1, klon
+!         IF (zt(i).LT.100.0) THEN
+         IF (zt(i).LT.50.0) THEN
+            jbad = jbad + 1
+            jadrs(jbad) = i
+         ENDIF
+         ENDDO
+#endif
+         IF (jbad .GT. 0) THEN
+           ok = .FALSE.
+           DO i = 1, jbad
+             PRINT *,'i,k,temperature,lon,lat,pourc ter,oce,lic,sic =',
+     $       jadrs(i),k,zt(jadrs(i)),rlon(jadrs(i)),rlat(jadrs(i)),
+     $       (pctsrf(jadrs(i),nsrf),nsrf=1,nbsrf)
+           ENDDO
+         ENDIF
+      ENDDO
+c
+      DO nsrf = 1, nbsrf
+         DO i = 1, klon
+            zt(i) = tsol(i,nsrf)
+         ENDDO
+#ifdef CRAY
+         CALL WHENFGT(klon, zt, 1, 370.0, jadrs, jbad)
+#else
+         jbad = 0
+         DO i = 1, klon
+         IF (zt(i).GT.370.0) THEN
+            jbad = jbad + 1
+            jadrs(jbad) = i
+         ENDIF
+         ENDDO
+#endif
+         IF (jbad .GT. 0) THEN
+           ok = .FALSE.
+           DO i = 1, jbad
+            PRINT *,'i,nsrf,temperature,lon,lat,pourc ter,oce,lic,sic ='
+     $      ,jadrs(i),nsrf,zt(jadrs(i)),rlon(jadrs(i)),rlat(jadrs(i))
+     $      ,pctsrf(jadrs(i),nsrf)
+           ENDDO
+         ENDIF
+#ifdef CRAY
+         CALL WHENFLT(klon, zt, 1, 100.0, jadrs, jbad)
+#else
+         jbad = 0
+         DO i = 1, klon
+!         IF (zt(i).LT.100.0) THEN
+         IF (zt(i).LT.50.0) THEN
+            jbad = jbad + 1
+            jadrs(jbad) = i
+         ENDIF
+         ENDDO
+#endif
+         IF (jbad .GT. 0) THEN
+           ok = .FALSE.
+           DO i = 1, jbad
+            PRINT *,'i,nsrf,temperature,lon,lat,pourc ter,oce,lic,sic ='
+     $      ,jadrs(i),nsrf,zt(jadrs(i)),rlon(jadrs(i)),rlat(jadrs(i))
+     $      ,pctsrf(jadrs(i),nsrf)
+           ENDDO
+         ENDIF
+      ENDDO
+c
+      IF (.NOT. ok) THEN
+         PRINT*, 'hgardfou s arrete ', text
+         CALL abort
+      ENDIF
+
+      RETURN
+      END
Index: LMDZ4/tags/V3_4/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/tags/V3_4/libf/phylmd/physiq.F	(revision 944)
+++ LMDZ4/tags/V3_4/libf/phylmd/physiq.F	(revision 944)
@@ -0,0 +1,3395 @@
+!
+! $Header$
+!
+c
+c#define IO_DEBUG
+
+      SUBROUTINE physiq (nlon,nlev,nqmax,
+     .            debut,lafin,rjourvrai,gmtime,pdtphys,
+     .            paprs,pplay,pphi,pphis,presnivs,clesphy0,
+     .            u,v,t,qx,
+     .            omega,
+#ifdef INCA
+     .            flxmass_w,
+#endif
+     .            d_u, d_v, d_t, d_qx, d_ps
+     .            , dudyn
+     .            , PVteta)
+
+      USE ioipsl
+      USE comgeomphy
+      USE write_field_phy
+      USE dimphy
+      USE mod_grid_phy_lmdz
+      USE mod_phys_lmdz_para
+      USE iophy
+      USE misc_mod, mydebug=>debug
+      USE vampir
+      USE pbl_surface_mod, ONLY : pbl_surface
+
+
+      USE ocean_slab_mod, ONLY   : ocean_slab_get_vars
+      USE ocean_cpl_mod, ONLY    : ocean_cpl_get_vars
+      USE ocean_forced_mod, ONLY : ocean_forced_get_vars
+      USE fonte_neige_mod, ONLY  : fonte_neige_get_vars
+
+      IMPLICIT none
+c======================================================================
+c
+c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+c
+c Objet: Moniteur general de la physique du modele
+cAA      Modifications quant aux traceurs :
+cAA                  -  uniformisation des parametrisations ds phytrac
+cAA                  -  stockage des moyennes des champs necessaires
+cAA                     en mode traceur off-line 
+c======================================================================
+c   CLEFS CPP POUR LES IO
+c   =====================
+#define histhf
+#define histday
+#define histmth
+c#define histins
+#define histmthNMC
+#define histISCCP
+c======================================================================
+c    modif   ( P. Le Van ,  12/10/98 )
+c
+c  Arguments:
+c
+c nlon----input-I-nombre de points horizontaux
+c nlev----input-I-nombre de couches verticales
+c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
+c debut---input-L-variable logique indiquant le premier passage
+c lafin---input-L-variable logique indiquant le dernier passage
+c rjour---input-R-numero du jour de l'experience
+c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
+c pdtphys-input-R-pas d'integration pour la physique (seconde)
+c paprs---input-R-pression pour chaque inter-couche (en Pa)
+c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
+c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
+c pphis---input-R-geopotentiel du sol
+c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
+c u-------input-R-vitesse dans la direction X (de O a E) en m/s
+c v-------input-R-vitesse Y (de S a N) en m/s
+c t-------input-R-temperature (K)
+c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
+c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
+c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
+c omega---input-R-vitesse verticale en Pa/s
+c d_u-----output-R-tendance physique de "u" (m/s/s)
+c d_v-----output-R-tendance physique de "v" (m/s/s)
+c d_t-----output-R-tendance physique de "t" (K/s)
+c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
+c d_ps----output-R-tendance physique de la pression au sol
+cIM
+c PVteta--output-R-vorticite potentielle a des thetas constantes
+c======================================================================
+#include "dimensions.h"
+      integer jjmp1
+      parameter (jjmp1=jjm+1-1/jjm)
+      integer iip1
+      parameter (iip1=iim+1)
+
+#include "regdim.h"
+#include "indicesol.h"
+#include "dimsoil.h"
+#include "clesphys.h"
+#include "control.h"
+#include "logic.h"
+#include "temps.h"
+cym#include "comgeomphy.h"
+#include "advtrac.h"
+#include "iniprint.h"
+#include "thermcell.h"
+c======================================================================
+      LOGICAL ok_cvl  ! pour activer le nouveau driver pour convection KE
+      PARAMETER (ok_cvl=.TRUE.)
+      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
+      PARAMETER (ok_gust=.FALSE.)
+c======================================================================
+      LOGICAL check ! Verifier la conservation du modele en eau
+      PARAMETER (check=.FALSE.)
+      LOGICAL ok_stratus ! Ajouter artificiellement les stratus
+      PARAMETER (ok_stratus=.FALSE.)
+c======================================================================
+      logical rnpb
+#ifdef INCA
+      parameter(rnpb=.false.)
+#else
+      parameter(rnpb=.true.)
+#endif
+c      ocean = type de modele ocean a utiliser: force, slab, couple
+      character*6 ocean
+      SAVE ocean
+c$OMP THREADPRIVATE(ocean)
+
+cIM "slab" ocean
+      REAL tslab(klon)    !Temperature du slab-ocean
+      REAL seaice(klon)   !glace de mer (kg/m2) 
+      REAL fluxo(klon)    !flux turbulents ocean-glace de mer 
+      REAL fluxg(klon)    !flux turbulents ocean-atmosphere
+      REAL amn, amx
+c======================================================================
+c Clef controlant l'activation du cycle diurne:
+ccc      LOGICAL cycle_diurne
+ccc      PARAMETER (cycle_diurne=.FALSE.)
+c======================================================================
+c Modele thermique du sol, a activer pour le cycle diurne:
+ccc      LOGICAL soil_model
+ccc      PARAMETER (soil_model=.FALSE.)
+      logical ok_veget
+      save ok_veget
+c$OMP THREADPRIVATE(ok_veget)
+c     parameter (ok_veget = .true.)
+c      parameter (ok_veget = .false.)
+c======================================================================
+c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
+c le calcul du rayonnement est celle apres la precipitation des nuages.
+c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
+c la condensation et la precipitation. Cette cle augmente les impacts
+c radiatifs des nuages.
+ccc      LOGICAL new_oliq
+ccc      PARAMETER (new_oliq=.FALSE.)
+c======================================================================
+c Clefs controlant deux parametrisations de l'orographie:
+cc      LOGICAL ok_orodr
+ccc      PARAMETER (ok_orodr=.FALSE.)
+ccc      LOGICAL ok_orolf
+ccc      PARAMETER (ok_orolf=.FALSE.)
+c======================================================================
+      LOGICAL ok_journe ! sortir le fichier journalier
+      save ok_journe
+c$OMP THREADPRIVATE(ok_journe)
+c
+      LOGICAL ok_mensuel ! sortir le fichier mensuel
+      save ok_mensuel
+c$OMP THREADPRIVATE(ok_mensuel)
+c
+      LOGICAL ok_instan ! sortir le fichier instantane
+      save ok_instan
+c$OMP THREADPRIVATE(ok_instan)
+c
+      LOGICAL ok_region ! sortir le fichier regional
+      PARAMETER (ok_region=.FALSE.)
+c======================================================================
+c     pour phsystoke avec thermiques
+      REAL fm_therm(klon,klev+1)
+      REAL entr_therm(klon,klev)
+c======================================================================
+c
+      INTEGER ivap          ! indice de traceurs pour vapeur d'eau
+      PARAMETER (ivap=1)
+      INTEGER iliq          ! indice de traceurs pour eau liquide
+      PARAMETER (iliq=2)
+
+c
+c
+c Variables argument:
+c
+      INTEGER nlon
+      INTEGER nlev
+      INTEGER nqmax
+      REAL rjourvrai
+      REAL gmtime
+      REAL pdtphys
+      LOGICAL debut, lafin
+      REAL paprs(klon,klev+1)
+      REAL pplay(klon,klev)
+      REAL pphi(klon,klev)
+      REAL pphis(klon)
+      REAL presnivs(klev)
+      REAL znivsig(klev)
+      real pir
+
+      REAL u(klon,klev)
+      REAL v(klon,klev)
+      REAL t(klon,klev)
+      REAL qx(klon,klev,nqmax)
+
+      REAL,allocatable,save :: t_ancien(:,:), q_ancien(:,:)
+c$OMP THREADPRIVATE(t_ancien, q_ancien)
+      LOGICAL ancien_ok
+      SAVE ancien_ok
+c$OMP THREADPRIVATE(ancien_ok)
+      REAL d_t_dyn(klon,klev)
+      REAL d_q_dyn(klon,klev)
+
+      REAL omega(klon,klev)
+
+#ifdef INCA
+      REAL flxmass_w(klon,klev)
+#endif
+      REAL d_u(klon,klev)
+      REAL d_v(klon,klev)
+      REAL d_t(klon,klev)
+      REAL d_qx(klon,klev,nqmax)
+      REAL d_ps(klon)
+      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
+c
+cIM Amip2 PV a theta constante 
+c
+      INTEGER nbteta
+      PARAMETER(nbteta=3)
+      CHARACTER*3 ctetaSTD(nbteta)
+      DATA ctetaSTD/'350','380','405'/
+      SAVE ctetaSTD
+c$OMP THREADPRIVATE(ctetaSTD)
+      REAL rtetaSTD(nbteta)
+      DATA rtetaSTD/350., 380., 405./
+      SAVE rtetaSTD
+c$OMP THREADPRIVATE(rtetaSTD)     
+c
+      REAL PVteta(klon,nbteta)
+      REAL zx_tmp_3dte(iim,jjmp1,nbteta)
+c
+cMI Amip2 PV a theta constante
+
+cym      INTEGER klevp1, klevm1
+cym      PARAMETER(klevp1=klev+1,klevm1=klev-1)
+cym#include "raddim.h"
+c
+
+      REAL,allocatable,save :: swdn0(:,:), swdn(:,:)
+      REAL,allocatable,save :: swup0(:,:), swup(:,:)
+c$OMP THREADPRIVATE(swdn0 , swdn, swup0, swup)
+c
+      REAL,allocatable,save :: SWdn200clr(:), SWdn200(:)
+      REAL,allocatable,save :: SWup200clr(:), SWup200(:)
+c$OMP THREADPRIVATE(SWdn200clr, SWdn200, SWup200clr, SWup200)
+c
+      REAL,allocatable,save :: lwdn0(:,:), lwdn(:,:)
+      REAL,allocatable,save :: lwup0(:,:), lwup(:,:)
+c$OMP THREADPRIVATE(lwdn0 , lwdn, lwup0, lwup)
+c
+      REAL,allocatable,save :: LWdn200clr(:), LWdn200(:)
+      REAL,allocatable,save :: LWup200clr(:), LWup200(:)
+c$OMP THREADPRIVATE(LWdn200clr, LWdn200, LWup200clr, LWup200)
+c
+      REAL,allocatable,save :: LWdnTOA(:), LWdnTOAclr(:)
+c$OMP THREADPRIVATE(LWdnTOA, LWdnTOAclr)
+c
+cIM Amip2
+c variables a une pression donnee
+c
+      integer nlevSTD
+      PARAMETER(nlevSTD=17)
+      real rlevSTD(nlevSTD)
+      DATA rlevSTD/100000., 92500., 85000., 70000.,
+     .60000., 50000., 40000., 30000., 25000., 20000.,
+     .15000., 10000., 7000., 5000., 3000., 2000., 1000./
+      SAVE rlevstd
+c$OMP THREADPRIVATE(rlevSTD)
+      CHARACTER*4 clevSTD(nlevSTD)
+      DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
+     .'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
+     .'70  ','50  ','30  ','20  ','10  '/
+      SAVE clevSTD
+c$OMP THREADPRIVATE(clevSTD)
+c
+      CHARACTER*4 bb2
+      CHARACTER*2 bb3
+c
+      real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
+      real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
+      real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
+      real wlevSTD(klon,nlevSTD) 
+c
+c nout : niveau de output des variables a une pression donnee
+      INTEGER nout
+      PARAMETER(nout=3) !nout=1 : day; =2 : mth; =3 : NMC
+c
+      REAL,SAVE,ALLOCATABLE :: tsumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: usumSTD(:,:,:), vsumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: wsumSTD(:,:,:), phisumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: qsumSTD(:,:,:), rhsumSTD(:,:,:)
+c$OMP THREADPRIVATE(tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD) 
+c$OMP THREADPRIVATE(qsumSTD, rhsumSTD)
+c
+      logical oknondef(klon,nlevSTD,nout)
+      real,SAVE,ALLOCATABLE :: tnondef(:,:,:) 
+c$OMP THREADPRIVATE(tnondef)
+c
+c les produits uvSTD, vqSTD, .., T2STD sont calcules
+c a partir des valeurs instantannees toutes les 6 h
+c qui sont moyennees sur le mois
+c
+      real uvSTD(klon,nlevSTD)
+      real vqSTD(klon,nlevSTD)
+      real vTSTD(klon,nlevSTD)
+      real wqSTD(klon,nlevSTD)
+c
+      real,save,allocatable :: uvsumSTD(:,:,:)
+      real,save,allocatable :: vqsumSTD(:,:,:)
+      real,save,allocatable :: vTsumSTD(:,:,:)
+      real,save,allocatable :: wqsumSTD(:,:,:)
+c
+      real vphiSTD(klon,nlevSTD)
+      real wTSTD(klon,nlevSTD)
+      real u2STD(klon,nlevSTD)
+      real v2STD(klon,nlevSTD)
+      real T2STD(klon,nlevSTD)
+c
+      real,save,allocatable :: vphisumSTD(:,:,:)
+      real,save,allocatable :: wTsumSTD(:,:,:)
+      real,save,allocatable :: u2sumSTD(:,:,:)
+      real,save,allocatable :: v2sumSTD(:,:,:)
+      real,save,allocatable :: T2sumSTD(:,:,:)
+c
+c$OMP THREADPRIVATE(uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD)
+c$OMP THREADPRIVATE(vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD)
+
+cMI Amip2
+c
+#include "radepsi.h"
+#include "radopt.h"
+c
+c
+c prw: precipitable water
+      real prw(klon)
+
+      REAL convliq(klon,klev)  ! eau liquide nuageuse convective
+      REAL convfra(klon,klev)  ! fraction nuageuse convective
+
+      REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
+      REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
+      REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
+      REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
+
+      INTEGER linv, kp1
+c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
+c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
+      REAL flwp(klon), fiwp(klon)
+      REAL flwc(klon,klev), fiwc(klon,klev)
+      REAL flwp_c(klon), fiwp_c(klon)
+      REAL flwc_c(klon,klev), fiwc_c(klon,klev)
+      REAL flwp_s(klon), fiwp_s(klon)
+      REAL flwc_s(klon,klev), fiwc_s(klon,klev)
+
+cIM ISCCP simulator v3.4
+c dans clesphys.h top_height, overlap
+cv3.4
+      INTEGER debug, debugcol
+cym      INTEGER npoints
+cym      PARAMETER(npoints=klon) 
+c
+      INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
+      INTEGER nregISCtot
+      PARAMETER(nregISCtot=1) 
+c
+c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
+c y compris pour 1 point
+c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
+c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
+      INTEGER imin_debut, nbpti
+      INTEGER jmin_debut, nbptj 
+cIM parametres ISCCP BEG
+      INTEGER nbapp_isccp
+!     INTEGER nbapp_isccp,isccppas
+!     PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
+!                           !i.e. toutes les 3 heures 
+      INTEGER n, napisccp
+c     PARAMETER(napisccp=3)
+      PARAMETER(napisccp=1)
+      INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
+      DATA ifreq_isccp/3/
+      SAVE ifreq_isccp
+c$OMP THREADPRIVATE(ifreq_isccp)
+      CHARACTER*5 typinout(napisccp)
+      DATA typinout/'i3od'/
+      SAVE typinout
+c$OMP THREADPRIVATE(typinout)
+cIM verif boxptop BEG
+      CHARACTER*1 verticaxe(napisccp)
+      DATA verticaxe/'1'/ 
+      SAVE verticaxe
+c$OMP THREADPRIVATE(verticaxe)
+cIM verif boxptop END
+      INTEGER nvlev(napisccp)
+c     INTEGER nvlev
+      REAL t1, aa
+      REAL seed_re(klon,napisccp)
+      INTEGER,ALLOCATABLE,SAVE :: seed_old(:,:)
+c$OMP THREADPRIVATE(seed_old)
+cym !!!! A voir plus tard 
+cym      INTEGER iphy(iim,jjmp1)
+cIM parametres ISCCP END
+c
+c ncol = nb. de sous-colonnes pour chaque maille du GCM 
+c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM 
+c      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
+      INTEGER,SAVE :: ncol(napisccp)
+      INTEGER ncolmx, seed(klon,napisccp)
+      REAL nbsunlit(nregISCtot,klon,napisccp)  !nbsunlit : moyenne de sunlit
+c     PARAMETER(ncolmx=1500)
+      PARAMETER(ncolmx=300)
+c
+cIM verif boxptop BEG
+      REAL vertlev(ncolmx,napisccp)
+cIM verif boxptop END
+c
+      REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
+      INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
+c$OMP THREADPRIVATE(tautab,invtau)
+      REAL emsfc_lw
+      PARAMETER(emsfc_lw=0.99)
+c     REAL    ran0                      ! type for random number fuction
+c
+      REAL cldtot(klon,klev)
+c variables de haut en bas pour le simulateur ISCCP
+      REAL dtau_s(klon,klev) !tau nuages startiformes
+      REAL dtau_c(klon,klev) !tau nuages convectifs
+      REAL dem_s(klon,klev)  !emissivite nuages startiformes 
+      REAL dem_c(klon,klev)  !emissivite nuages convectifs
+c
+c variables de haut en bas pour le simulateur ISCCP
+      REAL pfull(klon,klev)
+      REAL phalf(klon,klev+1)
+      REAL qv(klon,klev)
+      REAL cc(klon,klev)
+      REAL conv(klon,klev)
+      REAL dtau_sH2B(klon,klev)
+      REAL dtau_cH2B(klon,klev)
+      REAL at(klon,klev)
+      REAL dem_sH2B(klon,klev)
+      REAL dem_cH2B(klon,klev)
+c
+      INTEGER kmax, lmax, lmax3
+      PARAMETER(kmax=8, lmax=8, lmax3=3)
+      INTEGER kmaxm1, lmaxm1
+      PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
+      INTEGER iimx7, jjmx7, jjmp1x7
+      PARAMETER(iimx7=iim*kmaxm1, jjmx7=jjm*lmaxm1, 
+     .jjmp1x7=jjmp1*lmaxm1)
+c
+c output from ISCCP simulator
+      REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
+      REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
+      REAL totalcldarea(klon,napisccp) 
+      REAL meanptop(klon,napisccp)
+      REAL meantaucld(klon,napisccp)
+      REAL boxtau(klon,ncolmx,napisccp)
+      REAL boxptop(klon,ncolmx,napisccp) 
+      REAL zx_tmp_fi3d_bx(klon,ncolmx)
+      REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
+c
+      REAL cld_fi3d(klon,lmax3)
+      REAL cld_3d(iim,jjmp1,lmax3)
+c
+      INTEGER iw, iwmax
+      REAL wmin, pas_w
+c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
+cIM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
+      PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
+      REAL o500(klon)
+c
+cIM: nbregdyn = nbre regions pour calculs statistiques sur output du ISCCP
+cIM: dynamiques  
+      INTEGER nreg, nbregdyn
+      PARAMETER(nbregdyn=5)
+
+      INTEGER,ALLOCATABLE,SAVE :: pct_ocean(:,:)
+c$OMP THREADPRIVATE(pct_ocean)
+cym      SAVE pct_ocean
+ 
+c sorties ISCCP
+
+      integer nid_isccp
+      save nid_isccp        
+c$OMP THREADPRIVATE(nid_isccp)
+
+      REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
+      DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
+      SAVE zx_tau
+      DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
+      SAVE zx_pc
+c$OMP THREADPRIVATE(zx_tau,zx_pc)
+c cldtopres pression au sommet des nuages
+      REAL cldtopres(lmaxm1), cldtopres3(lmax3)
+      DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
+      DATA cldtopres3/440., 680., 1000./
+      SAVE cldtopres,cldtopres3
+c$OMP THREADPRIVATE(cldtopres,cldtopres3)
+cIM 051005 BEG
+      REAL tmp_his1_3d(iwmax,kmaxm1,lmaxm1,nbregdyn,napisccp) 
+      REAL tmp_his2_3d(iwmax,kmaxm1,lmaxm1,nbregdyn,napisccp) 
+      REAL tmp_his3_3d(iwmax,kmaxm1,lmaxm1,nbregdyn,napisccp) 
+      INTEGER komega, nhoriRD 
+
+c taulev: numero du niveau de tau dans les sorties ISCCP
+      CHARACTER *4 taulev(kmaxm1)
+c     DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
+      DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
+      CHARACTER *3 pclev(lmaxm1)
+      DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
+      SAVE taulev,pclev
+c$OMP THREADPRIVATE(taulev,pclev)
+c
+c cnameisccp
+      CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
+cIM bad 151205     DATA cnameisccp/'pc< 50hPa, tau< 0.3', 
+      DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
+     .                'pc= 180-310hPa, tau< 0.3',
+     .                'pc= 310-440hPa, tau< 0.3',
+     .                'pc= 440-560hPa, tau< 0.3',
+     .                'pc= 560-680hPa, tau< 0.3',
+     .                'pc= 680-800hPa, tau< 0.3',
+     .                'pc= 800-1000hPa, tau< 0.3',
+     .                'pc= 50-180hPa, tau= 0.3-1.3',
+     .                'pc= 180-310hPa, tau= 0.3-1.3',
+     .                'pc= 310-440hPa, tau= 0.3-1.3',
+     .                'pc= 440-560hPa, tau= 0.3-1.3',
+     .                'pc= 560-680hPa, tau= 0.3-1.3',
+     .                'pc= 680-800hPa, tau= 0.3-1.3',
+     .                'pc= 800-1000hPa, tau= 0.3-1.3',
+     .                'pc= 50-180hPa, tau= 1.3-3.6',
+     .                'pc= 180-310hPa, tau= 1.3-3.6',
+     .                'pc= 310-440hPa, tau= 1.3-3.6',
+     .                'pc= 440-560hPa, tau= 1.3-3.6',
+     .                'pc= 560-680hPa, tau= 1.3-3.6',
+     .                'pc= 680-800hPa, tau= 1.3-3.6',
+     .                'pc= 800-1000hPa, tau= 1.3-3.6',
+     .                'pc= 50-180hPa, tau= 3.6-9.4',
+     .                'pc= 180-310hPa, tau= 3.6-9.4',
+     .                'pc= 310-440hPa, tau= 3.6-9.4',
+     .                'pc= 440-560hPa, tau= 3.6-9.4',
+     .                'pc= 560-680hPa, tau= 3.6-9.4',
+     .                'pc= 680-800hPa, tau= 3.6-9.4',
+     .                'pc= 800-1000hPa, tau= 3.6-9.4',
+     .                'pc= 50-180hPa, tau= 9.4-23',
+     .                'pc= 180-310hPa, tau= 9.4-23',
+     .                'pc= 310-440hPa, tau= 9.4-23',
+     .                'pc= 440-560hPa, tau= 9.4-23',
+     .                'pc= 560-680hPa, tau= 9.4-23',
+     .                'pc= 680-800hPa, tau= 9.4-23',
+     .                'pc= 800-1000hPa, tau= 9.4-23',
+     .                'pc= 50-180hPa, tau= 23-60',
+     .                'pc= 180-310hPa, tau= 23-60',
+     .                'pc= 310-440hPa, tau= 23-60',
+     .                'pc= 440-560hPa, tau= 23-60',
+     .                'pc= 560-680hPa, tau= 23-60',
+     .                'pc= 680-800hPa, tau= 23-60',
+     .                'pc= 800-1000hPa, tau= 23-60',
+     .                'pc= 50-180hPa, tau> 60.',
+     .                'pc= 180-310hPa, tau> 60.',
+     .                'pc= 310-440hPa, tau> 60.',
+     .                'pc= 440-560hPa, tau> 60.',
+     .                'pc= 560-680hPa, tau> 60.',
+     .                'pc= 680-800hPa, tau> 60.',
+     .                'pc= 800-1000hPa, tau> 60.'/
+       SAVE cnameisccp
+c$OMP THREADPRIVATE(cnameisccp)
+c
+c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
+c     INTEGER nhorix7
+cIM: region='3d' <==> sorties en global
+      CHARACTER*3 region
+      PARAMETER(region='3d')
+c
+cIM ISCCP simulator v3.4
+c
+      logical ok_hf
+c
+      integer nid_hf, nid_hf3d
+      save ok_hf, nid_hf, nid_hf3d
+c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
+c  QUESTION : noms de variables ?
+
+c#ifdef histhf
+c      data ok_hf/.true./
+c#else
+c      data ok_hf/.false./
+c#endif
+      INTEGER        longcles
+      PARAMETER    ( longcles = 20 )
+      REAL clesphy0( longcles      )
+c
+c Variables quasi-arguments
+c
+      REAL xjour
+      SAVE xjour
+c$OMP THREADPRIVATE(xjour)
+c
+c
+c Variables propres a la physique
+c
+      REAL dtime
+      SAVE dtime                  ! pas temporel de la physique
+c$OMP THREADPRIVATE(dtime)
+c
+      INTEGER radpas
+      SAVE radpas                 ! frequence d'appel rayonnement
+c$OMP THREADPRIVATE(radpas)
+c
+      REAL,allocatable,save :: radsol(:)
+c$OMP THREADPRIVATE(radsol)
+cym      SAVE radsol               ! bilan radiatif au sol calcule par code radiatif
+c
+      REAL,allocatable,save :: rlat(:)
+c$OMP THREADPRIVATE(rlat)
+cym      SAVE rlat                   ! latitude pour chaque point
+c
+      REAL,allocatable,save :: rlon(:)
+c$OMP THREADPRIVATE(rlon)
+cym      SAVE rlon                   ! longitude pour chaque point
+
+      REAL,SAVE,ALLOCATABLE :: rlonPOS(:)
+c$OMP THREADPRIVATE(rlonPOS)    
+cym      SAVE rlonPOS                ! longitudes > 0. pour chaque point
+c
+cc      INTEGER iflag_con
+cc      SAVE iflag_con              ! indicateur de la convection
+c
+      INTEGER itap
+      SAVE itap                   ! compteur pour la physique
+c$OMP THREADPRIVATE(itap)
+c
+      REAL co2_ppm_etat0
+c
+      REAL solaire_etat0
+c
+      real slp(klon) ! sea level pressure
+
+      REAL,allocatable,save :: ftsol(:,:)
+c$OMP THREADPRIVATE(ftsol)       ! temperature du sol
+
+cIM
+      REAL,SAVE,ALLOCATABLE :: newsst(:) !temperature de l'ocean
+c$OMP THREADPRIVATE(newsst)
+cym     SAVE newsst
+c
+      REAL fevap(klon,nbsrf)
+      REAL fluxlat(klon,nbsrf)
+c
+      REAL,allocatable,save :: deltat(:)
+c$OMP THREADPRIVATE(deltat)
+cym      SAVE deltat                 ! ecart avec la SST de reference
+c
+      REAL qsol(klon)
+c
+      REAL,allocatable,save :: falbe(:,:)
+c$OMP THREADPRIVATE(falbe)       ! albedo par type de surface
+c
+      REAL,allocatable,save :: falblw(:,:)
+c$OMP THREADPRIVATE(falblw)       ! albedo par type de surface
+
+c
+c
+c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
+c
+      REAL,allocatable,save :: zmea(:)
+c$OMP THREADPRIVATE(zmea)
+cym      SAVE zmea                   ! orographie moyenne
+c
+      REAL,allocatable,save :: zstd(:)
+c$OMP THREADPRIVATE(zstd)
+cym      SAVE zstd                   ! deviation standard de l'OESM
+c
+      REAL,allocatable,save :: zsig(:)
+c$OMP THREADPRIVATE(zsig)
+cym      SAVE zsig                   ! pente de l'OESM
+c
+      REAL,allocatable,save :: zgam(:)
+c$OMP THREADPRIVATE(zgam)
+cym      save zgam                   ! anisotropie de l'OESM
+c
+      REAL,allocatable,save :: zthe(:)
+c$OMP THREADPRIVATE(zthe)      
+cym      SAVE zthe                   ! orientation de l'OESM
+c
+      REAL,allocatable,save :: zpic(:)
+c$OMP THREADPRIVATE(zpic)
+cym      SAVE zpic                   ! Maximum de l'OESM
+c
+      REAL,allocatable,save :: zval(:)
+c$OMP THREADPRIVATE(zval)
+cym      SAVE zval                   ! Minimum de l'OESM
+c
+      REAL,allocatable,save :: rugoro(:)
+c$OMP THREADPRIVATE(rugoro)
+cym      SAVE rugoro                 ! longueur de rugosite de l'OESM
+c
+cIM 141004     REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
+      REAL zulow(klon),zvlow(klon)
+c
+      REAL,allocatable,save :: zuthe(:),zvthe(:)
+c$OMP THREADPRIVATE(zuthe,zvthe)
+cym      SAVE zuthe
+cym      SAVE zvthe
+      INTEGER igwd,idx(klon),itest(klon)
+c
+      REAL agesno(klon,nbsrf)
+c
+      REAL,allocatable,save :: alb_neig(:)
+c$OMP THREADPRIVATE(alb_neig)
+cym      SAVE alb_neig               ! albedo de la neige
+c
+c      REAL,allocatable,save :: run_off_lic_0(:)
+cc$OMP THREADPRIVATE(run_off_lic_0)
+cym      SAVE run_off_lic_0
+cKE43
+c Variables liees a la convection de K. Emanuel (sb):
+c
+      REAL,allocatable,save :: ema_workcbmf(:)   ! cloud base mass flux
+c$OMP THREADPRIVATE(ema_workcbmf)
+cym      SAVE ema_workcbmf
+
+      REAL,allocatable,save :: ema_cbmf(:)       ! cloud base mass flux
+c$OMP THREADPRIVATE(ema_cbmf)
+cym      SAVE ema_cbmf
+
+      REAL,allocatable,save :: ema_pcb(:)        ! cloud base pressure
+c$OMP THREADPRIVATE(ema_pcb)
+cym      SAVE ema_pcb
+
+      REAL,allocatable,save :: ema_pct(:)        ! cloud top pressure
+c$OMP THREADPRIVATE(ema_pct)
+cym      SAVE ema_pct
+
+      REAL bas, top             ! cloud base and top levels
+      SAVE bas
+      SAVE top
+c$OMP THREADPRIVATE(bas, top)
+
+      REAL,allocatable,save :: Ma(:,:)        ! undilute upward mass flux
+c$OMP THREADPRIVATE(Ma)
+cym      SAVE Ma
+      REAL,allocatable,save :: qcondc(:,:)    ! in-cld water content from convect
+c$OMP THREADPRIVATE(qcondc)
+cym      SAVE qcondc 
+      REAL,allocatable,save :: ema_work1(:, :), ema_work2(:, :)
+c$OMP THREADPRIVATE(ema_work1,ema_work2)
+cym      SAVE ema_work1, ema_work2
+      REAL wdn(klon), tdn(klon), qdn(klon)
+
+      REAL,allocatable,save :: wd(:) ! sb
+c$OMP THREADPRIVATE(wd)
+cym      SAVE wd       ! sb
+
+c Variables locales pour la couche limite (al1):
+c
+cAl1      REAL pblh(klon)           ! Hauteur de couche limite
+cAl1      SAVE pblh
+c34EK
+c
+c Variables locales:
+c
+      REAL cdragh(klon) ! drag coefficient pour T and Q
+      REAL cdragm(klon) ! drag coefficient pour vent
+cAA
+cAA  Pour phytrac 
+cAA
+      REAL ycoefh(klon,klev)    ! coef d'echange pour phytrac
+      REAL yu1(klon)            ! vents dans la premiere couche U
+      REAL yv1(klon)            ! vents dans la premiere couche V
+
+      REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
+
+c@$$      LOGICAL offline           ! Controle du stockage ds "physique"
+c@$$      PARAMETER (offline=.false.)
+c@$$      INTEGER physid
+      REAL,allocatable,save :: pfrac_impa(:,:)! Produits des coefs lessivage impaction
+c$OMP THREADPRIVATE(pfrac_impa)
+cym      save pfrac_impa
+      REAL,allocatable,save :: pfrac_nucl(:,:)! Produits des coefs lessivage nucleation
+c$OMP THREADPRIVATE(pfrac_nucl)
+cym      save pfrac_nucl
+      REAL,allocatable,save :: pfrac_1nucl(:,:)! Produits des coefs lessi nucl (alpha = 1)
+c$OMP THREADPRIVATE(pfrac_1nucl)
+cym      save pfrac_1nucl
+      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
+      REAL frac_nucl(klon,klev) ! idem (nucleation)
+#ifdef INCA
+      INTEGER       :: iii
+      REAL          :: calday
+#endif
+
+cAA
+      REAL,allocatable,save :: rain_fall(:) ! pluie
+c$OMP THREADPRIVATE(rain_fall)
+      REAL,allocatable,save :: snow_fall(:) ! neige
+c$OMP THREADPRIVATE(snow_fall)
+cym      save snow_fall, rain_fall
+
+cIM cf FH pour Tiedtke 080604
+      REAL rain_tiedtke(klon),snow_tiedtke(klon)
+c
+
+      REAL,allocatable,save :: total_rain(:), nday_rain(:)
+c$OMP THREADPRIVATE(total_rain,nday_rain)
+cym      save total_rain, nday_rain
+cIM 050204 END
+      REAL evap(klon), devap(klon) ! evaporation et sa derivee
+      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
+
+      REAL bils(klon) ! bilan de chaleur au sol
+      REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
+C                             ! type de sous-surface et pondere par la fraction
+      REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
+C                             ! type de sous-surface et pondere par la fraction
+      REAL fder(klon)         
+      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
+      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
+      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
+      REAL uq(klon) ! integr. verticale du transport zonal de l'eau
+c
+      REAL frugs(klon,nbsrf)
+      REAL zxrugs(klon) ! longueur de rugosite
+c
+c Conditions aux limites
+c
+      INTEGER julien
+c
+      INTEGER lmt_pas
+      SAVE lmt_pas                ! frequence de mise a jour
+c$OMP THREADPRIVATE(lmt_pas) 
+      REAL,allocatable,save :: pctsrf(:,:)
+c$OMP THREADPRIVATE(pctsrf) 
+cIM
+      REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
+
+cym      REAL paire_ter(klon)        !surfaces terre 
+      REAL,allocatable,save ::  paire_ter(:)        !surfaces terre 
+c$OMP THREADPRIVATE(paire_ter) 
+    
+cIM
+cym      SAVE pctsrf                 ! sous-fraction du sol
+      REAL,allocatable,save :: albsol(:)
+c$OMP THREADPRIVATE(albsol)
+cym      SAVE albsol                 ! albedo du sol total
+      REAL,allocatable,save :: albsollw(:)
+c$OMP THREADPRIVATE(albsollw)      
+cym      SAVE albsollw                 ! albedo du sol total
+
+      REAL,allocatable,save :: wo(:,:)
+c$OMP THREADPRIVATE(wo)
+cym     SAVE wo                     ! ozone
+
+cIM sorties
+      REAL un_jour
+      PARAMETER(un_jour=86400.)
+c======================================================================
+c
+c Declaration des procedures appelees
+c
+      EXTERNAL angle     ! calculer angle zenithal du soleil
+      EXTERNAL alboc     ! calculer l'albedo sur ocean
+      EXTERNAL ajsec     ! ajustement sec
+      EXTERNAL conlmd    ! convection (schema LMD)
+cKE43
+      EXTERNAL conema3  ! convect4.3
+      EXTERNAL fisrtilp  ! schema de condensation a grande echelle (pluie)
+cAA 
+      EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
+c                          ! stockage des coefficients necessaires au
+c                          ! lessivage OFF-LINE et ON-LINE
+      EXTERNAL hgardfou  ! verifier les temperatures
+      EXTERNAL nuage     ! calculer les proprietes radiatives
+      EXTERNAL o3cm      ! initialiser l'ozone
+      EXTERNAL orbite    ! calculer l'orbite terrestre
+      EXTERNAL ozonecm   ! prescrire l'ozone
+      EXTERNAL phyetat0  ! lire l'etat initial de la physique
+      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
+      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
+      EXTERNAL suphec    ! initialiser certaines constantes
+      EXTERNAL transp    ! transport total de l'eau et de l'energie
+      EXTERNAL ecribina  ! ecrire le fichier binaire global
+      EXTERNAL ecribins  ! ecrire le fichier binaire global
+      EXTERNAL ecrirega  ! ecrire le fichier binaire regional
+      EXTERNAL ecriregs  ! ecrire le fichier binaire regional
+cIM
+      EXTERNAL haut2bas  !variables de haut en bas
+      INTEGER lnblnk1
+      EXTERNAL lnblnk1   !enleve les blancs a la fin d'une variable de type
+                         !caracter
+      EXTERNAL ini_undefSTD  !initialise a 0 une variable a 1 niveau de pression
+      EXTERNAL undefSTD      !somme les valeurs definies d'1 var a 1 niveau de pression
+c     EXTERNAL moy_undefSTD  !moyenne d'1 var a 1 niveau de pression
+c     EXTERNAL moyglo_aire   !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
+c                            !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
+c
+c Variables locales
+c
+      real,allocatable,save :: clwcon(:,:),rnebcon(:,:)
+c$OMP THREADPRIVATE(clwcon,rnebcon)
+      real,allocatable,save :: clwcon0(:,:),rnebcon0(:,:)
+cym      save rnebcon, clwcon
+c$OMP THREADPRIVATE(clwcon0,rnebcon0)
+      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
+      REAL dialiq(klon,klev)  ! eau liquide nuageuse
+      REAL diafra(klon,klev)  ! fraction nuageuse
+      REAL cldliq(klon,klev)  ! eau liquide nuageuse
+      REAL cldfra(klon,klev)  ! fraction nuageuse
+      REAL cldtau(klon,klev)  ! epaisseur optique
+      REAL cldemi(klon,klev)  ! emissivite infrarouge
+c
+CXXX PB 
+      REAL fluxq(klon,klev, nbsrf)   ! flux turbulent d'humidite
+      REAL fluxt(klon,klev, nbsrf)   ! flux turbulent de chaleur
+      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
+      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
+c
+      REAL zxfluxt(klon, klev)
+      REAL zxfluxq(klon, klev)
+      REAL zxfluxu(klon, klev)
+      REAL zxfluxv(klon, klev)
+CXXX
+      REAL,allocatable,save :: heat(:,:)    ! chauffage solaire
+c$OMP THREADPRIVATE(heat)
+      REAL,allocatable,save :: heat0(:,:)   ! chauffage solaire ciel clair
+c$OMP THREADPRIVATE(heat0)
+      REAL,allocatable,save :: cool(:,:)    ! refroidissement infrarouge
+c$OMP THREADPRIVATE(cool)
+      REAL,allocatable,save :: cool0(:,:)   ! refroidissement infrarouge ciel clair
+c$OMP THREADPRIVATE(cool0)
+      REAL,allocatable,save :: topsw(:), toplw(:), solsw(:), sollw(:)
+c$OMP THREADPRIVATE(topsw,toplw,solsw,sollw)
+      real,allocatable,save :: sollwdown(:)    ! downward LW flux at surface
+c$OMP THREADPRIVATE(sollwdown)
+cIM BEG
+      real,allocatable,save :: sollwdownclr(:)    ! downward CS LW flux at surface 
+c$OMP THREADPRIVATE(sollwdownclr)
+      real,allocatable,save :: toplwdown(:)       ! downward CS LW flux at TOA
+c$OMP THREADPRIVATE(toplwdown)
+      real,allocatable,save :: toplwdownclr(:)    ! downward CS LW flux at TOA
+c$OMP THREADPRIVATE(toplwdownclr)
+cIM END
+      REAL,allocatable,save :: topsw0(:),toplw0(:),solsw0(:),sollw0(:)
+c$OMP THREADPRIVATE( topsw0,toplw0,solsw0,sollw0)
+      REAL,allocatable,save :: albpla(:)
+c$OMP THREADPRIVATE(albpla)
+      REAL fsollw(klon, nbsrf)   ! bilan flux IR pour chaque sous surface
+      REAL fsolsw(klon, nbsrf)   ! flux solaire absorb. pour chaque sous surface
+c Le rayonnement n'est pas calcule tous les pas, il faut donc
+c                      sauvegarder les sorties du rayonnement
+cym      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
+cym      SAVE  sollwdownclr, toplwdown, toplwdownclr
+cym      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
+c
+      INTEGER itaprad
+      SAVE itaprad
+c$OMP THREADPRIVATE(itaprad)
+c
+      REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
+      REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
+c
+      REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
+      REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
+c
+      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
+      REAL zxsnow_dummy(klon)
+c
+      REAL dist, rmu0(klon), fract(klon)
+      REAL zdtime, zlongi
+c
+      CHARACTER*2 str2
+      CHARACTER*2 iqn
+c
+      REAL qcheck
+      REAL z_avant(klon), z_apres(klon), z_factor(klon)
+      LOGICAL zx_ajustq
+c
+      REAL za, zb
+      REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
+      real zqsat(klon,klev)
+      INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq
+      REAL t_coup
+      PARAMETER (t_coup=234.0)
+c
+      REAL zphi(klon,klev)
+cym A voir plus tard !!
+cym      REAL zx_relief(iim,jjmp1)
+cym      REAL zx_aire(iim,jjmp1)
+c
+c Grandeurs de sorties
+      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
+      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
+      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
+      REAL s_trmb3(klon)
+cKE43
+c Variables locales pour la convection de K. Emanuel (sb):
+c
+      REAL upwd(klon,klev)      ! saturated updraft mass flux
+      REAL dnwd(klon,klev)      ! saturated downdraft mass flux
+      REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
+      REAL tvp(klon,klev)       ! virtual temp of lifted parcel
+      REAL,allocatable,save :: cape(:)           ! CAPE
+c$OMP THREADPRIVATE(cape)
+cym      SAVE cape
+      CHARACTER*40 capemaxcels  !max(CAPE)
+
+      REAL,allocatable,save :: pbase(:)          ! cloud base pressure
+c$OMP THREADPRIVATE(pbase)
+cym      SAVE pbase
+      REAL,allocatable,save :: bbase(:)          ! cloud base buoyancy
+c$OMP THREADPRIVATE(bbase)
+cym      SAVE bbase
+      REAL rflag(klon)          ! flag fonctionnement de convect
+      INTEGER iflagctrl(klon)          ! flag fonctionnement de convect
+c -- convect43:
+      INTEGER ntra              ! nb traceurs pour convect4.3
+      REAL pori_con(klon)    ! pressure at the origin level of lifted parcel
+      REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
+      REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
+      REAL dplcldt(klon), dplcldr(klon)
+c?     .     condm_con(klon,klev),conda_con(klon,klev),
+c?     .     mr_con(klon,klev),ep_con(klon,klev)
+c?     .    ,sadiab(klon,klev),wadiab(klon,klev)
+c --
+c34EK
+c
+c Variables du changement
+c
+c con: convection
+c lsc: condensation a grande echelle (Large-Scale-Condensation)
+c ajs: ajustement sec
+c eva: evaporation de l'eau liquide nuageuse
+c vdf: couche limite (Vertical DiFfusion)
+      REAL d_t_con(klon,klev),d_q_con(klon,klev)
+      REAL,SAVE,ALLOCATABLE :: d_u_con(:,:),d_v_con(:,:)
+c$OMP THREADPRIVATE(d_u_con,d_v_con)
+      REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
+      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
+      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
+      REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
+      REAL d_t_oli(klon,klev) !tendances dues a oro et lif 
+      REAL rneb(klon,klev)
+c
+*********************************************************
+*     declarations
+      real,save,allocatable :: zqasc(:,:)
+c$OMP THREADPRIVATE(zqasc)
+cym      save zqasc
+      
+*********************************************************
+cIM 081204 END
+c
+      REAL pmfu(klon,klev), pmfd(klon,klev)
+      REAL pen_u(klon,klev), pen_d(klon,klev)
+      REAL pde_u(klon,klev), pde_d(klon,klev)
+      INTEGER kcbot(klon), kctop(klon), kdtop(klon)
+      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
+      REAL prfl(klon,klev+1), psfl(klon,klev+1)
+c
+      INTEGER,allocatable,save :: ibas_con(:), itop_con(:)
+c$OMP THREADPRIVATE(ibas_con,itop_con)
+cym
+cym      SAVE ibas_con,itop_con
+cym
+      REAL,SAVE,ALLOCATABLE :: rain_con(:)
+c$OMP THREADPRIVATE(rain_con)
+      REAL rain_lsc(klon)
+      REAL,SAVE,ALLOCATABLE :: snow_con(:)
+c$OMP THREADPRIVATE(snow_con)
+      REAL snow_lsc(klon)
+      REAL d_ts(klon,nbsrf)
+c
+      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
+      REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
+c
+      REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
+      REAL d_t_oro(klon,klev)
+      REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
+      REAL d_t_lif(klon,klev)
+      REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif 
+
+      REAL,allocatable,save :: ratqs(:,:)
+c$OMP THREADPRIVATE(ratqs)
+      REAL ratqss(klon,klev),ratqsc(klon,klev)
+      real ratqsbas,ratqshaut
+cym      save ratqsbas,ratqshaut, ratqs
+      save ratqsbas,ratqshaut
+c$OMP THREADPRIVATE(ratqsbas,ratqshaut)
+      real zpt_conv(klon,klev)
+
+c Parametres lies au nouveau schema de nuages (SB, PDF)
+      real fact_cldcon
+      real facttemps
+      logical ok_newmicro
+      save ok_newmicro
+c$OMP THREADPRIVATE(ok_newmicro)
+      save fact_cldcon,facttemps
+c$OMP THREADPRIVATE(fact_cldcon,facttemps)
+      real facteur
+
+      integer iflag_cldcon
+      save iflag_cldcon
+c$OMP THREADPRIVATE(iflag_cldcon)
+      logical ptconv(klon,klev)
+cIM cf. AM 081204 BEG
+      logical ptconvth(klon,klev)
+cIM cf. AM 081204 END
+c
+c Variables liees a l'ecriture de la bande histoire physique
+c
+c======================================================================
+c
+cIM cf. AM 081204 BEG
+c   declarations pour sortir sur une sous-region
+      integer imin_ins,imax_ins,jmin_ins,jmax_ins
+      save imin_ins,imax_ins,jmin_ins,jmax_ins
+c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
+c      real lonmin_ins,lonmax_ins,latmin_ins
+c     s  ,latmax_ins
+c     data lonmin_ins,lonmax_ins,latmin_ins
+c    s  ,latmax_ins/
+c valeurs initiales     s   -5.,20.,41.,55./   
+c    s   100.,130.,-20.,20./
+c    s   -180.,180.,-90.,90./
+c======================================================================
+cIM cf. AM 081204 END
+
+c
+      integer itau_w   ! pas de temps ecriture = itap + itau_phy
+c
+c
+c Variables locales pour effectuer les appels en serie
+c
+      REAL t_seri(klon,klev), q_seri(klon,klev)
+      REAL ql_seri(klon,klev),qs_seri(klon,klev)
+      REAL u_seri(klon,klev), v_seri(klon,klev)
+c
+      REAL tr_seri(klon,klev,nbtr)
+      REAL d_tr(klon,klev,nbtr)
+
+      REAL zx_rh(klon,klev)
+cIM RH a 2m (la surface)
+      REAL rh2m(klon), qsat2m(klon)
+      REAL tpot(klon), tpote(klon)
+      REAL Lheat
+
+      INTEGER        length
+      PARAMETER    ( length = 100 )
+      REAL tabcntr0( length       )
+c
+      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
+cIM
+      INTEGER ndex2d1(iwmax)
+c
+cIM AMIP2 BEG
+      REAL moyglo, mountor
+cIM 141004 BEG
+      REAL zustrdr(klon), zvstrdr(klon)
+      REAL zustrli(klon), zvstrli(klon)
+      REAL zustrph(klon), zvstrph(klon)
+      REAL aam, torsfc
+cIM 141004 END
+cIM 190504 BEG
+      INTEGER ij, imp1jmp1
+      PARAMETER(imp1jmp1=(iim+1)*jjmp1)
+cym A voir plus tard
+      REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
+      REAL padyn(iim+1,jjmp1,klev+1)
+      REAL dudyn(iim+1,jjmp1,klev)
+      REAL rlatdyn(iim+1,jjmp1)
+cIM 190504 END
+      LOGICAL ok_msk
+      REAL msk(klon)
+cIM 
+      REAL airetot, pi
+cym A voir plus tard
+cym      REAL zm_wo(jjmp1, klev)
+cIM AMIP2 END
+c
+      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
+      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
+c#ifdef histmthNMC
+cym   A voir plus tard !!!!
+cym      REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
+      REAL zx_tmp_fiNC(klon,nlevSTD) 
+c#endif
+      REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D 
+      REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
+      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
+c
+      INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
+      INTEGER nid_ctesGCM
+      SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
+      SAVE nid_ctesGCM
+c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
+c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
+c
+cIM 280405 BEG
+      INTEGER nid_bilKPins, nid_bilKPave
+      SAVE nid_bilKPins, nid_bilKPave
+c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
+c
+      REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
+      REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
+      REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
+      REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
+c
+cIM 280405 END
+c
+      INTEGER nhori, nvert, nvert1, nvert3
+      REAL zsto, zsto1, zsto2
+      REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
+      REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
+      REAL zout_isccp(napisccp)
+      SAVE zcals, zcalh, zoutj, zout_isccp
+c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
+
+      real zjulian
+      save zjulian
+c$OMP THREADPRIVATE(zjulian)
+
+      character*20 modname
+      character*80 abort_message
+      logical ok_sync
+      real date0
+      integer idayref
+
+C essai writephys
+      integer fid_day, fid_mth, fid_ins
+      parameter (fid_ins = 1, fid_day = 2, fid_mth = 3) 
+      integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
+      parameter (prof2d_on = 1, prof3d_on = 2,
+     .           prof2d_av = 3, prof3d_av = 4)
+      character*30 nom_fichier
+      character*10 varname
+      character*40 vartitle
+      character*20 varunits
+C     Variables liees au bilan d'energie et d'enthalpi
+      REAL ztsol(klon)
+      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
+     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
+     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
+c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
+c$OMP+              h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
+      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
+      REAL      d_h_vcol_phy
+      REAL      fs_bound, fq_bound
+      SAVE      d_h_vcol_phy
+c$OMP THREADPRIVATE(d_h_vcol_phy)
+      REAL      zero_v(klon)
+      CHARACTER*15 ztit
+      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
+      SAVE      ip_ebil
+      DATA      ip_ebil/0/
+c$OMP THREADPRIVATE(ip_ebil)
+      INTEGER   if_ebil ! level for energy conserv. dignostics
+      SAVE      if_ebil
+c$OMP THREADPRIVATE(if_ebil)
+c+jld ec_conser
+      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
+      REAL ZRCPD
+c-jld ec_conser
+      REAL t2m(klon,nbsrf)  ! temperature a 2m
+      REAL q2m(klon,nbsrf)  ! humidite a 2m
+
+cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
+      REAL,SAVE,ALLOCATABLE :: u10m(:,:), v10m(:,:) !vents a 10m
+c$OMP THREADPRIVATE(u10m,v10m)
+      REAL zt2m(klon), zq2m(klon)             !temp., hum. 2m moyenne s/ 1 maille
+      REAL zu10m(klon), zv10m(klon)           !vents a 10m moyennes s/1 maille
+      CHARACTER*40 t2mincels, t2maxcels       !t2m min., t2m max
+      CHARACTER*40 tinst, tave, typeval
+cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
+      REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
+      REAL,allocatable,save :: sulfate_pi(:,:) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
+c$OMP THREADPRIVATE(sulfate_pi)
+cym      SAVE sulfate_pi
+
+      REAL cldtaupi(klon,klev)  ! Cloud optical thickness for pre-industrial (pi) aerosols
+
+      REAL re(klon, klev)       ! Cloud droplet effective radius
+      REAL fl(klon, klev)  ! denominator of re
+
+      REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
+
+      ! Aerosol optical properties
+      REAL,SAVE,ALLOCATABLE :: tau_ae(:,:,:), piz_ae(:,:,:)
+c$OMP THREADPRIVATE(tau_ae,piz_ae)
+      REAL,SAVE,ALLOCATABLE :: cg_ae(:,:,:)
+c$OMP THREADPRIVATE(cg_ae)
+
+#ifdef INCA_AER
+#ifdef CPP_COUPLE
+! Aerosol optical properties by INCA model 
+      REAL, SAVE, ALLOCATABLE  ::    tau_inca(:,:,:,:)
+      REAL, SAVE, ALLOCATABLE  ::    piz_inca(:,:,:,:)
+      REAL, SAVE, ALLOCATABLE  ::    cg_inca(:,:,:,:)
+      REAL, SAVE, ALLOCATABLE  ::    ccm(:,:,:) 
+      CHARACTER*4              ::    rfname(9) 
+
+      REAL,SAVE,ALLOCATABLE :: topswad_inca(:), solswad_inca(:) ! Aerosol direct effect.
+      REAL,SAVE,ALLOCATABLE :: topswad0_inca(:), solswad0_inca(:) ! Aerosol direct effect.
+      REAL,SAVE,ALLOCATABLE :: topswai_inca(:), solswai_inca(:) ! Aerosol indirect effect.
+      REAL,SAVE,ALLOCATABLE :: topsw_inca(:,:), solsw_inca(:,:) 
+      REAL,SAVE,ALLOCATABLE :: topsw0_inca(:,:), solsw0_inca(:,:)
+#endif
+#endif
+      REAL,SAVE,ALLOCATABLE :: topswad(:), solswad(:) ! Aerosol direct effect.
+c$OMP THREADPRIVATE(topswad,solswad)
+      ! ok_ade=T -ADE=topswad-topsw
+
+      REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) ! Aerosol indirect effect.
+c$OMP THREADPRIVATE(topswai,solswai)
+      ! ok_aie=T ->
+      !        ok_ade=T -AIE=topswai-topswad
+      !        ok_ade=F -AIE=topswai-topsw
+
+      REAL aerindex(klon)       ! POLDER aerosol index
+     
+      ! Parameters
+      LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
+      REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
+cym 
+      SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
+c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
+
+c
+c Declaration des constantes et des fonctions thermodynamiques
+c
+      LOGICAL,SAVE :: first=.true.
+c$OMP THREADPRIVATE(first)
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+cIM 100106 BEG : pouvoir sortir les ctes de la physique
+#include "conema3.h"
+#include "fisrtilp.h"
+#include "nuage.h"
+#include "compbl.h"
+cIM 100106 END : pouvoir sortir les ctes de la physique
+c
+c======================================================================
+
+cym => necessaire pour iflag_con != 2    
+      pmfd(:,:) = 0.
+      pen_u(:,:) = 0.
+      pen_d(:,:) = 0.
+      pde_d(:,:) = 0.
+      pde_u(:,:) = 0.
+      aam=0.
+      torsfc=0.
+
+      if (first) then 
+      
+      allocate( t_ancien(klon,klev), q_ancien(klon,klev))
+      allocate( swdn0(klon,klevp1), swdn(klon,klevp1))
+      allocate( swup0(klon,klevp1), swup(klon,klevp1))
+      allocate( SWdn200clr(klon), SWdn200(klon))
+      allocate( SWup200clr(klon), SWup200(klon))
+      allocate( lwdn0(klon,klevp1), lwdn(klon,klevp1))
+      allocate( lwup0(klon,klevp1), lwup(klon,klevp1))
+      allocate( LWdn200clr(klon), LWdn200(klon))
+      allocate( LWup200clr(klon), LWup200(klon))
+      allocate( LWdnTOA(klon), LWdnTOAclr(klon))
+      allocate( radsol(klon))
+      allocate( rlat(klon))
+      allocate( rlon(klon))
+      allocate( ftsol(klon,nbsrf))
+      allocate( deltat(klon))
+      allocate( falbe(klon,nbsrf))
+      allocate( falblw(klon,nbsrf))
+      allocate( zmea(klon))
+      allocate( zstd(klon))
+      allocate( zsig(klon))
+      allocate( zgam(klon))
+      allocate( zthe(klon))
+      allocate( zpic(klon))
+      allocate( zval(klon))
+      allocate( rugoro(klon))
+      allocate( zuthe(klon),zvthe(klon))
+      allocate( alb_neig(klon))
+      allocate( ema_workcbmf(klon))   
+      allocate( ema_cbmf(klon)) 
+      allocate( ema_pcb(klon))
+      allocate( ema_pct(klon))  
+      allocate( Ma(klon,klev) )
+      allocate( qcondc(klon,klev))  
+      allocate( ema_work1(klon, klev), ema_work2(klon, klev))
+      allocate( wd(klon) )
+      allocate( pfrac_impa(klon,klev))
+      allocate( pfrac_nucl(klon,klev))
+      allocate( pfrac_1nucl(klon,klev))
+      allocate( rain_fall(klon) )
+      allocate( snow_fall(klon) )
+      allocate( total_rain(klon), nday_rain(klon))
+      allocate( pctsrf(klon,nbsrf))
+      allocate( albsol(klon))
+      allocate( albsollw(klon))
+      allocate( wo(klon,klev))
+      allocate( clwcon(klon,klev),rnebcon(klon,klev))
+      allocate( heat(klon,klev)    )
+      allocate( heat0(klon,klev)  )
+      allocate( cool(klon,klev)    )
+      allocate( cool0(klon,klev)   )
+      allocate( topsw(klon), toplw(klon), solsw(klon), sollw(klon))
+      allocate( sollwdown(klon)    )
+      allocate( sollwdownclr(klon)  ) 
+      allocate( toplwdown(klon)      )
+      allocate( toplwdownclr(klon)   )
+      allocate( topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon))
+      allocate( albpla(klon))
+      allocate( cape(klon)   )        
+      allocate( pbase(klon)   )     
+      allocate( bbase(klon)    )     
+      allocate( ibas_con(klon), itop_con(klon))
+      allocate( ratqs(klon,klev))
+      allocate( sulfate_pi(klon, klev)) 
+      allocate( paire_ter(klon))
+      allocate(tsumSTD(klon,nlevSTD,nout))
+      allocate(usumSTD(klon,nlevSTD,nout))
+      allocate(vsumSTD(klon,nlevSTD,nout))
+      allocate(wsumSTD(klon,nlevSTD,nout))
+      allocate(phisumSTD(klon,nlevSTD,nout))
+      allocate(qsumSTD(klon,nlevSTD,nout))
+      allocate(rhsumSTD(klon,nlevSTD,nout))
+      allocate(uvsumSTD(klon,nlevSTD,nout))
+      allocate(vqsumSTD(klon,nlevSTD,nout))
+      allocate(vTsumSTD(klon,nlevSTD,nout))
+      allocate(wqsumSTD(klon,nlevSTD,nout))
+      allocate( vphisumSTD(klon,nlevSTD,nout))
+      allocate( wTsumSTD(klon,nlevSTD,nout))
+      allocate( u2sumSTD(klon,nlevSTD,nout))
+      allocate( v2sumSTD(klon,nlevSTD,nout))
+      allocate( T2sumSTD(klon,nlevSTD,nout))
+      allocate( seed_old(klon,napisccp))
+      allocate( pct_ocean(klon,nbregdyn))
+      allocate( rlonPOS(klon))
+      allocate( newsst(klon))
+      allocate( zqasc(klon,klev))
+      allocate( rain_con(klon))
+      allocate( u10m(klon,nbsrf), v10m(klon,nbsrf)) 
+      allocate( topswad(klon), solswad(klon))
+      allocate( topswai(klon), solswai(klon) )
+#ifdef INCA_AER
+#ifdef CPP_COUPLE
+      allocate( topswad_inca(klon), solswad_inca(klon))
+      allocate( topswad0_inca(klon), solswad0_inca(klon))
+      allocate( topswai_inca(klon), solswai_inca(klon))
+      allocate( topsw_inca(klon,9), solsw_inca(klon,9))
+      allocate( topsw0_inca(klon,9), solsw0_inca(klon,9))
+      allocate( tau_inca(klon,klev,9,2))
+      allocate( piz_inca(klon,klev,9,2))
+      allocate( cg_inca(klon,klev,9,2))
+      allocate( ccm(klon,klev,2) )
+#endif
+#endif
+      allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
+      allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
+      allocate( cg_ae(klon,klev,2))
+      allocate( snow_con(klon))
+      allocate( tnondef(klon,nlevSTD,nout)) 
+      allocate( d_u_con(klon,klev),d_v_con(klon,klev))            
+      
+      
+        paire_ter(:)=0.    
+        clwcon(:,:)=0.
+	rnebcon(:,:)=0.
+	ratqs(:,:)=0.
+	sollw(:)=0.
+	ema_work1(:,:)=0.
+	ema_work2(:,:)=0.
+cym Attention pbase pas initialise dans concvl !!!!
+        pbase(:)=0
+	first=.false.
+      endif 
+
+                  
+       modname = 'physiq'
+cIM
+      IF (ip_ebil_phy.ge.1) THEN
+        DO i=1,klon
+          zero_v(i)=0.
+        END DO 
+      END IF 
+      ok_sync=.TRUE.
+      IF (nqmax .LT. 2) THEN
+         abort_message = 'eaux vapeur et liquide sont indispensables'
+         CALL abort_gcm (modname,abort_message,1)
+      ENDIF
+      IF (debut) THEN
+         CALL suphec ! initialiser constantes et parametres phys.
+      ENDIF
+
+
+c======================================================================
+      xjour = rjourvrai
+c
+c Si c'est le debut, il faut initialiser plusieurs choses
+c          ********
+c
+       IF (debut) THEN
+C
+!rv
+         u10m(:,:)=0.
+         v10m(:,:)=0.
+         piz_ae(:,:,:)=0.
+         tau_ae(:,:,:)=0.
+         cg_ae(:,:,:)=0.
+         rain_con(:)=0.
+         snow_con(:)=0.
+         bl95_b0=0.
+         bl95_b1=0.
+         topswai(:)=0.
+         topswad(:)=0.
+         solswai(:)=0.
+         solswad(:)=0.
+#ifdef INCA_AER
+#ifdef CPP_COUPLE
+         tau_inca(:,:,:,:) = 0.
+         piz_inca(:,:,:,:) = 0.
+         cg_inca(:,:,:,:)  = 0.
+         ccm(:,:,:)        = 0.
+         topswai_inca(:)   = 0.
+         topswad_inca(:)   = 0.
+         topswad0_inca(:)  = 0.
+         topsw_inca(:,:)   = 0.
+         topsw0_inca(:,:)  = 0.
+         solswai_inca(:)   = 0.
+         solswad_inca(:)   = 0.
+         solswad0_inca(:)  = 0.
+         solsw_inca(:,:)   = 0.
+         solsw0_inca(:,:)  = 0.
+#endif
+#endif
+!rv
+!ACo
+         d_u_con(:,:) = 0.0
+         d_v_con(:,:) = 0.0
+         rnebcon0(:,:) = 0.0
+         clwcon0(:,:) = 0.0
+         rnebcon(:,:) = 0.0
+         clwcon(:,:) = 0.0
+         paire_ter(:) = 0.0
+c        nhistoW(:,:,:,:) = 0.0
+c        histoW(:,:,:,:) = 0.0
+
+cIM      
+         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
+c
+c appel a la lecture du run.def physique
+c
+         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
+     .                  ok_instan, ok_hf,
+     .                  fact_cldcon, facttemps,ok_newmicro,
+cIM  .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
+     .                  iflag_cldcon,ratqsbas,ratqshaut,
+     .                  ok_ade, ok_aie, 
+     .                  bl95_b0, bl95_b1,
+     .                  iflag_thermals,nsplit_thermals)
+
+c
+c
+c Initialiser les compteurs:
+c
+         itap    = 0
+         itaprad = 0
+
+         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
+     .       rlat,rlon,pctsrf, ftsol,
+     .       ocean, ok_veget,
+     .       falbe, falblw, rain_fall,snow_fall,
+     .       solsw, sollw,
+     .       radsol,clesphy0,
+     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
+     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon)
+
+       DO i=1,klon
+         IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) + 
+     $       pctsrf(i, is_oce) + pctsrf(i, is_sic)  - 1.) .GT. EPSFRA) 
+     $       THEN 
+            WRITE(*,*)
+     $   'physiq apres lecture de restart: pb sous surface au point ',
+     $   i, pctsrf(i, 1 : nbsrf)
+         ENDIF 
+      ENDDO
+
+         radpas = NINT( 86400./dtime/nbapp_rad)
+c
+C on remet le calendrier a zero
+c
+         IF (raz_date .eq. 1) THEN
+           itau_phy = 0
+         ENDIF
+
+cIM cf. AM 081204 BEG
+         PRINT*,'cycle_diurne3 =',cycle_diurne
+cIM cf. AM 081204 END
+c
+         CALL printflag( tabcntr0,radpas,ok_journe,
+     ,                    ok_instan, ok_region )
+c
+         IF (ABS(dtime-pdtphys).GT.0.001) THEN
+            WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
+     .                        pdtphys
+            abort_message='Pas physique n est pas correct '
+            call abort_gcm(modname,abort_message,1)
+         ENDIF
+         IF (nlon .NE. klon) THEN
+            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, 
+     .                      klon
+            abort_message='nlon et klon ne sont pas coherents'
+            call abort_gcm(modname,abort_message,1)
+         ENDIF
+         IF (nlev .NE. klev) THEN
+            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
+     .                       klev
+            abort_message='nlev et klev ne sont pas coherents'
+            call abort_gcm(modname,abort_message,1)
+         ENDIF
+c
+         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN 
+           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
+           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
+           abort_message='Nbre d appels au rayonnement insuffisant'
+           call abort_gcm(modname,abort_message,1)
+         ENDIF
+         WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
+         WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
+     .                   ok_cvl
+c
+cKE43
+c Initialisation pour la convection de K.E. (sb):
+         IF (iflag_con.GE.3) THEN
+
+         WRITE(lunout,*)"*** Convection de Kerry Emanuel 4.3  "
+         WRITE(lunout,*)
+     .      "On va utiliser le melange convectif des traceurs qui"
+         WRITE(lunout,*)"est calcule dans convect4.3"
+         WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
+
+          DO i = 1, klon
+           ema_cbmf(i) = 0.
+           ema_pcb(i)  = 0.
+           ema_pct(i)  = 0.
+           ema_workcbmf(i) = 0.
+          ENDDO
+cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
+          DO i = 1, klon
+           ibas_con(i) = 1
+           itop_con(i) = 1
+          ENDDO
+cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
+
+         ENDIF
+
+c34EK
+         IF (ok_orodr) THEN
+           DO i=1,klon
+             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
+           ENDDO
+
+cFRED           WRITE(lunout,*)'FRED rugoro ...'
+cFRED           WRITE(lunout,*) rugoro (416), rugoro (417), rugoro(418)
+cFRED     .                    , rugoro(419)
+cFRED           WRITE(lunout,*) rugoro (421), rugoro (422), rugoro(423)
+cFRED           WRITE(lunout,*) rugoro (424), rugoro (425), rugoro(426)
+cFRED           WRITE(lunout,*) rugoro (427), rugoro (428), rugoro(429)
+cFRED           WRITE(lunout,*) rugoro (430), rugoro (431), rugoro(567)
+cFRED     .                    , rugoro(568)
+cFRED           WRITE(lunout,*) rugoro (574), rugoro (575), rugoro(576)
+cFRED           WRITE(lunout,*) rugoro (577), rugoro (698), rugoro(699)
+cFRED           WRITE(lunout,*) rugoro (711), rugoro (712), rugoro(713)
+
+          WRITE(lunout,*)'FRED rugoro =0.'
+          rugoro=0.
+
+           CALL SUGWD(klon,klev,paprs,pplay)
+           DO i=1,klon
+             zuthe(i)=0.
+             zvthe(i)=0.
+             if(zstd(i).gt.10.)then
+               zuthe(i)=(1.-zgam(i))*cos(zthe(i))
+               zvthe(i)=(1.-zgam(i))*sin(zthe(i))
+             endif
+           ENDDO
+         ENDIF
+c
+c
+         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
+         WRITE(lunout,*)'La frequence de lecture surface est de ', 
+     .                   lmt_pas
+c
+cIM200505        ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
+c        IF (ok_mensuel) THEN
+c        WRITE(lunout,*)'La frequence de sortie mensuelle est de ', 
+c    .                   ecrit_mth
+c        ENDIF
+c        ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
+c        IF (ok_journe) THEN
+c        WRITE(lunout,*)'La frequence de sortie journaliere est de ',
+c    .                   ecrit_day
+c        ENDIF
+cIM 130904 BEG
+cIM 080205      ecrit_hf = 86400./dtime *0.25  ! toutes les 6h
+cIM 170305      
+c        ecrit_hf = 86400./dtime/12.  ! toutes les 2h
+cIM 230305      
+cIM200505        ecrit_hf = 86400./dtime *0.25  ! toutes les 6h
+c
+cIM200505        ecrit_hf2mth = ecrit_day/ecrit_hf*30
+c
+cIM200505        IF (ok_journe) THEN
+cIM200505        WRITE(lunout,*)'La frequence de sortie hf est de ',
+cIM200505    .                   ecrit_hf
+cIM200505        ENDIF
+cIM 130904 END
+ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
+ccc         ecrit_ins = NINT(86400./dtime *0.25)  ! 4 fois par jour
+c        ecrit_ins = NINT(86400./dtime/48.)  ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
+c        ecrit_ins = NINT(86400./dtime/12.)  ! toutes les deux heures
+cIM200505        ecrit_ins = NINT(86400./dtime/8.)  ! toutes les trois heures
+cIM200505        IF (ok_instan) THEN
+cIM200505        WRITE(lunout,*)'La frequence de sortie instant. est de ', 
+cIM200505    .                   ecrit_ins
+cIM200505        ENDIF
+cIM200505        ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
+cIM200505        IF (ok_region) THEN
+cIM200505        WRITE(lunout,*)'La frequence de sortie region est de ', 
+cIM200505    .                   ecrit_reg
+cIM200505        ENDIF
+cIM 030306 BEG
+cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
+cIM : ne pas modifier ecrit_hf2mth
+c
+         ecrit_hf2mth = 30*1/ecrit_hf 
+c ecrit_ins en secondes, chaque pas de temps de la physique
+         ecrit_ins = dtime
+cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
+         ecrit_hf = ecrit_hf * un_jour
+!IM
+         IF(ecrit_day.LE.1.) THEN
+          ecrit_day = ecrit_day * un_jour !en secondes
+         ENDIF
+!IM
+         ecrit_mth = ecrit_mth * un_jour
+         ecrit_reg = ecrit_reg * un_jour
+         ecrit_tra = ecrit_tra * un_jour
+         ecrit_ISCCP = ecrit_ISCCP * un_jour
+c
+         PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP',ecrit_hf,
+     .   ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP
+cIM 030306 END
+
+      capemaxcels = 't_max(X)'
+      t2mincels = 't_min(X)'
+      t2maxcels = 't_max(X)'
+      tinst = 'inst(X)'
+      tave = 'ave(X)'
+cIM cf. AM 081204 BEG
+      write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
+cIM cf. AM 081204 END
+c
+c=============================================================
+c   Initialisation des sorties
+c=============================================================
+
+#ifdef CPP_IOIPSL
+
+#ifdef histhf
+#include "ini_histhf.h"
+#endif
+
+#ifdef histday
+#include "ini_histday.h"
+cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
+c#include "ini_bilKP_ins.h"
+c#include "ini_bilKP_ave.h"
+#endif
+
+#ifdef histmth
+#include "ini_histmth.h"
+#endif
+
+#ifdef histins
+#include "ini_histins.h"
+#endif
+
+#ifdef histISCCP
+#include "ini_histISCCP.h"
+#endif
+
+#ifdef histmthNMC
+#include "ini_histmthNMC.h"
+#endif
+
+#include "ini_histday_seri.h"
+
+#include "ini_paramLMDZ_phy.h"
+
+#endif
+
+cXXXPB Positionner date0 pour initialisation de ORCHIDEE
+      date0 = zjulian
+C      date0 = day_ini
+      WRITE(*,*) 'physiq date0 : ',date0
+c
+c
+c
+c Prescrire l'ozone dans l'atmosphere
+c
+c
+cc         DO i = 1, klon
+cc         DO k = 1, klev
+cc            CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
+cc         ENDDO
+cc         ENDDO
+c
+#ifdef INCA
+           call VTe(VTphysiq)
+	   call VTb(VTinca)
+           iii = MOD(NINT(xjour),360)
+           calday = FLOAT(iii) + gmtime
+           WRITE(lunout,*) 'initial time ', xjour, calday
+#ifdef INCAINFO
+           WRITE(lunout,*) 'Appel CHEMINI ...'
+#endif
+           CALL chemini( 
+     $                   rg,
+     $                   ra,
+     $                   airephy,
+     $                   rlat,
+     $                   rlon,
+     $                   presnivs,
+     $                   calday,
+     $                   klon,
+     $                   nqmax,
+     $                   pdtphys,
+     $                   annee_ref,
+     $                   day_ini)
+#ifdef INCAINFO
+           WRITE(lunout,*) 'OK.'
+#endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
+#endif
+c
+      ENDIF
+c
+c   ****************     Fin  de   IF ( debut  )   ***************
+c
+c
+c Mettre a zero des variables de sortie (pour securite)
+c
+      DO i = 1, klon
+         d_ps(i) = 0.0
+      ENDDO
+      DO k = 1, klev
+      DO i = 1, klon
+         d_t(i,k) = 0.0
+         d_u(i,k) = 0.0
+         d_v(i,k) = 0.0
+      ENDDO
+      ENDDO
+      DO iq = 1, nqmax
+      DO k = 1, klev
+      DO i = 1, klon
+         d_qx(i,k,iq) = 0.0
+      ENDDO
+      ENDDO
+      ENDDO
+      da(:,:)=0.
+      mp(:,:)=0.
+      phi(:,:,:)=0.
+c
+c Ne pas affecter les valeurs entrees de u, v, h, et q
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         t_seri(i,k)  = t(i,k)
+         u_seri(i,k)  = u(i,k)
+         v_seri(i,k)  = v(i,k)
+         q_seri(i,k)  = qx(i,k,ivap)
+         ql_seri(i,k) = qx(i,k,iliq)
+         qs_seri(i,k) = 0.
+      ENDDO
+      ENDDO
+      IF (nqmax.GE.3) THEN
+      DO iq = 3, nqmax
+      DO  k = 1, klev
+      DO  i = 1, klon
+         tr_seri(i,k,iq-2) = qx(i,k,iq)
+      ENDDO
+      ENDDO
+      ENDDO
+      ELSE
+      DO k = 1, klev
+      DO i = 1, klon
+         tr_seri(i,k,1) = 0.0
+      ENDDO
+      ENDDO
+      ENDIF
+C
+      DO i = 1, klon
+        ztsol(i) = 0.
+      ENDDO
+      DO nsrf = 1, nbsrf
+        DO i = 1, klon
+          ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
+        ENDDO
+      ENDDO
+cIM
+      IF (ip_ebil_phy.ge.1) THEN 
+        ztit='after dynamic'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+C     Comme les tendances de la physique sont ajoute dans la dynamique,
+C     on devrait avoir que la variation d'entalpie par la dynamique
+C     est egale a la variation de la physique au pas de temps precedent.
+C     Donc la somme de ces 2 variations devrait etre nulle.
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
+     s      , fs_bound, fq_bound )
+      END IF 
+
+c Diagnostiquer la tendance dynamique
+c
+      IF (ancien_ok) THEN
+         DO k = 1, klev
+         DO i = 1, klon
+            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
+            d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
+         ENDDO
+         ENDDO
+      ELSE
+         DO k = 1, klev
+         DO i = 1, klon
+            d_t_dyn(i,k) = 0.0
+            d_q_dyn(i,k) = 0.0
+         ENDDO
+         ENDDO
+         ancien_ok = .TRUE.
+      ENDIF
+c
+c Ajouter le geopotentiel du sol:
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         zphi(i,k) = pphi(i,k) + pphis(i)
+      ENDDO
+      ENDDO
+c
+c Verifier les temperatures
+c
+cIM BEG
+      IF (check) THEN
+       amn=MIN(ftsol(1,is_ter),1000.)
+       amx=MAX(ftsol(1,is_ter),-1000.)
+       DO i=2, klon
+        amn=MIN(ftsol(i,is_ter),amn)
+        amx=MAX(ftsol(i,is_ter),amx)
+       ENDDO
+c
+       PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
+      ENDIF !(check) THEN
+cIM END
+c
+      CALL hgardfou(rlon,rlat,pctsrf,t_seri,ftsol,'debutphy')
+c
+cIM BEG
+      IF (check) THEN
+       amn=MIN(ftsol(1,is_ter),1000.)
+       amx=MAX(ftsol(1,is_ter),-1000.)
+       DO i=2, klon
+        amn=MIN(ftsol(i,is_ter),amn)
+        amx=MAX(ftsol(i,is_ter),amx)
+       ENDDO
+c
+       PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
+      ENDIF !(check) THEN
+cIM END
+c
+c Incrementer le compteur de la physique
+c
+      itap   = itap + 1
+      julien = MOD(NINT(xjour),360)
+      if (julien .eq. 0) julien = 360
+c
+c Mettre en action les conditions aux limites (albedo, sst, etc.).
+c Prescrire l'ozone et calculer l'albedo sur l'ocean.
+c
+      IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
+         WRITE(lunout,*)' PHYS cond  julien ',julien
+         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
+      ENDIF
+c
+c Re-evaporer l'eau liquide nuageuse
+c
+      DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
+      DO i = 1, klon
+         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
+c        zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
+         zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
+         zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
+         zb = MAX(0.0,ql_seri(i,k))
+         za = - MAX(0.0,ql_seri(i,k))
+     .                  * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
+         t_seri(i,k) = t_seri(i,k) + za
+         q_seri(i,k) = q_seri(i,k) + zb
+         ql_seri(i,k) = 0.0
+         d_t_eva(i,k) = za
+         d_q_eva(i,k) = zb
+      ENDDO
+      ENDDO
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after reevap'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+C
+      END IF 
+
+c
+C calculs necessaires au calcul de l'albedo dans l'interface
+c
+      CALL orbite(FLOAT(julien),zlongi,dist)
+      IF (cycle_diurne) THEN
+        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
+        CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
+      ELSE
+        rmu0 = -999.999
+      ENDIF
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c Appel au pbl_surface : Planetary Boudary Layer et Surface
+c Cela implique tous les interactions des sous-surfaces et la partie diffusion 
+c turbulent du couche limit. 
+c 
+c Certains varibales de sorties de pbl_surface sont utiliser que pour 
+c ecriture des fihiers hist_XXXX.nc, ces sont :
+c   qsol,      zq2m,      s_pblh,  s_lcl,
+c   s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
+c   s_therm,   s_trmb1,   s_trmb2, s_trmb3,
+c   zxrugs,    zu10m,     zv10m,   fder,
+c   zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+c   frugs,     agesno,    fsollw,  fsolsw,
+c   d_ts,      fevap,     fluxlat, t2m,
+c   wfbils,    wfbilo,    fluxt,   fluxu, fluxv,
+c
+c Certains ne sont pas utiliser du tout : 
+c   dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
+c
+      CALL pbl_surface( 
+     e     dtime,     date0,     itap,    julien,
+     e     debut,     lafin,
+     e     rlon,      rlat,      rugoro,  rmu0,     
+     e     rain_fall, snow_fall, solsw,   sollw,    
+     e     t_seri,    q_seri,    u_seri,  v_seri,   
+     e     pplay,     paprs,     pctsrf,            
+     +     ftsol,     falbe,     falblw,  u10m,   v10m,
+     s     sollwdown, cdragh,    cdragm,  yu1,    yv1,
+     s     albsol,    albsollw,  sens,    evap,  
+     s     zxtsol,    zxfluxlat, zt2m,    qsat2m, 
+     s     d_t_vdf,   d_q_vdf,   d_u_vdf, d_v_vdf,
+     s     ycoefh,    pctsrf_new,                
+     d     qsol,      zq2m,      s_pblh,  s_lcl,
+     d     s_capCL,   s_oliqCL,  s_cteiCL,s_pblT,
+     d     s_therm,   s_trmb1,   s_trmb2, s_trmb3,
+     d     zxrugs,    zu10m,     zv10m,   fder,
+     d     zxqsurf,   rh2m,      zxfluxu, zxfluxv,
+     d     frugs,     agesno,    fsollw,  fsolsw,
+     d     d_ts,      fevap,     fluxlat, t2m,
+     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
+     -     dsens,     devap,     zxsnow,
+     -     zxfluxt,   zxfluxq,   q2m,     fluxq )
+c
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      pctsrf(:,:) = pctsrf_new(:,:)
+      
+      DO k = 1, klev
+      DO i = 1, klon
+         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
+         q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
+         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
+         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
+      ENDDO
+      ENDDO
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after surface_main'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, sens
+     e      , evap  , zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+
+c
+c Appeler la convection (au choix)
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         conv_q(i,k) = d_q_dyn(i,k) 
+     .               + d_q_vdf(i,k)/dtime
+         conv_t(i,k) = d_t_dyn(i,k) 
+     .               + d_t_vdf(i,k)/dtime
+      ENDDO
+      ENDDO
+      IF (check) THEN
+         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+         WRITE(lunout,*) "avantcon=", za
+      ENDIF
+      zx_ajustq = .FALSE.
+      IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
+      IF (zx_ajustq) THEN
+         DO i = 1, klon
+            z_avant(i) = 0.0
+         ENDDO
+         DO k = 1, klev
+         DO i = 1, klon
+            z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
+     .                        *(paprs(i,k)-paprs(i,k+1))/RG
+         ENDDO
+         ENDDO
+      ENDIF
+      IF (iflag_con.EQ.1) THEN
+          stop'reactiver le call conlmd dans physiq.F'
+c     CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
+c    .             d_t_con, d_q_con,
+c    .             rain_con, snow_con, ibas_con, itop_con)
+      ELSE IF (iflag_con.EQ.2) THEN
+      CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
+     e            conv_t, conv_q, -evap, omega,
+     s            d_t_con, d_q_con, rain_con, snow_con,
+     s            pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
+     s            kcbot, kctop, kdtop, pmflxr, pmflxs)
+      WHERE (rain_con < 0.) rain_con = 0.
+      WHERE (snow_con < 0.) snow_con = 0.
+      DO i = 1, klon
+         ibas_con(i) = klev+1 - kcbot(i)
+         itop_con(i) = klev+1 - kctop(i)
+      ENDDO
+      ELSE IF (iflag_con.GE.3) THEN
+c nb of tracers for the KE convection:
+c MAF la partie traceurs est faite dans phytrac
+c on met ntra=1 pour limiter les appels mais on peut
+c supprimer les calculs / ftra.
+              ntra = 1
+c sb, oct02:
+c Schema de convection modularise et vectorise:
+c (driver commun aux versions 3 et 4)
+c
+          IF (ok_cvl) THEN ! new driver for convectL
+
+          CALL concvl (iflag_con,
+     .        dtime,paprs,pplay,t_seri,q_seri,
+     .        u_seri,v_seri,tr_seri,ntra,
+     .        ema_work1,ema_work2,
+     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
+     .        rain_con, snow_con, ibas_con, itop_con,
+     .        upwd,dnwd,dnwd0,
+     .        Ma,cape,tvp,iflagctrl,
+     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
+     .        pmflxr,pmflxs,
+     .        da,phi,mp)
+
+cIM cf. FH
+              clwcon0=qcondc
+              pmfu(:,:)=upwd(:,:)+dnwd(:,:)
+
+          ELSE ! ok_cvl
+c MAF conema3 ne contient pas les traceurs
+          CALL conema3 (dtime,
+     .        paprs,pplay,t_seri,q_seri,
+     .        u_seri,v_seri,tr_seri,ntra,
+     .        ema_work1,ema_work2,
+     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
+     .        rain_con, snow_con, ibas_con, itop_con,
+     .        upwd,dnwd,dnwd0,bas,top,
+     .        Ma,cape,tvp,rflag,
+     .        pbase
+     .        ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
+     .        ,clwcon0)
+
+          ENDIF ! ok_cvl
+
+c
+c Correction precip
+          rain_con = rain_con * cvl_corr
+          snow_con = snow_con * cvl_corr
+c
+
+           IF (.NOT. ok_gust) THEN
+           do i = 1, klon
+            wd(i)=0.0
+           enddo
+           ENDIF
+
+c =================================================================== c
+c Calcul des proprietes des nuages convectifs
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         zx_t = t_seri(i,k)
+         IF (thermcep) THEN
+            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
+            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
+            zx_qs  = MIN(0.5,zx_qs)
+            zcor   = 1./(1.-retv*zx_qs)
+            zx_qs  = zx_qs*zcor
+         ELSE
+           IF (zx_t.LT.t_coup) THEN
+              zx_qs = qsats(zx_t)/pplay(i,k)
+           ELSE
+              zx_qs = qsatl(zx_t)/pplay(i,k)
+           ENDIF
+         ENDIF
+         zqsat(i,k)=zx_qs
+      ENDDO
+      ENDDO
+
+c   calcul des proprietes des nuages convectifs
+             clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
+             call clouds_gno
+     s       (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
+
+c =================================================================== c
+
+          DO i = 1, klon
+            ema_pcb(i)  = pbase(i)
+          ENDDO
+          DO i = 1, klon
+            ema_pct(i)  = paprs(i,itop_con(i))
+          ENDDO
+          DO i = 1, klon
+            ema_cbmf(i) = ema_workcbmf(i)
+          ENDDO      
+      ELSE
+          WRITE(lunout,*) "iflag_con non-prevu", iflag_con
+          CALL abort
+      ENDIF
+
+c     CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
+c    .              d_u_con, d_v_con)
+
+      DO k = 1, klev
+        DO i = 1, klon
+         t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
+         q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
+         u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
+         v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
+        ENDDO
+      ENDDO
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after convect'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, rain_con, snow_con, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+C
+      IF (check) THEN
+          za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+          WRITE(lunout,*)"aprescon=", za
+          zx_t = 0.0
+          za = 0.0
+          DO i = 1, klon
+            za = za + airephy(i)/FLOAT(klon)
+            zx_t = zx_t + (rain_con(i)+
+     .                   snow_con(i))*airephy(i)/FLOAT(klon)
+          ENDDO
+          zx_t = zx_t/za*dtime
+          WRITE(lunout,*)"Precip=", zx_t
+      ENDIF
+      IF (zx_ajustq) THEN
+          DO i = 1, klon
+            z_apres(i) = 0.0
+          ENDDO
+          DO k = 1, klev
+            DO i = 1, klon
+              z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
+     .            *(paprs(i,k)-paprs(i,k+1))/RG
+            ENDDO
+          ENDDO
+          DO i = 1, klon
+            z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
+     .          /z_apres(i)
+          ENDDO
+          DO k = 1, klev
+            DO i = 1, klon
+              IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
+     .            z_factor(i).LT.(1.0-1.0E-08)) THEN
+                  q_seri(i,k) = q_seri(i,k) * z_factor(i)
+              ENDIF
+            ENDDO
+          ENDDO
+      ENDIF
+      zx_ajustq=.FALSE.
+c
+c===================================================================
+c Convection seche (thermiques ou ajustement)
+c===================================================================
+c
+      d_t_ajs(:,:)=0.
+      d_u_ajs(:,:)=0.
+      d_v_ajs(:,:)=0.
+      d_q_ajs(:,:)=0.
+      fm_therm(:,:)=0.
+      entr_therm(:,:)=0.
+c
+      IF(prt_level>9)WRITE(lunout,*)
+     .    'AVANT LA CONVECTION SECHE , iflag_thermals='
+     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+      if(iflag_thermals.lt.0) then
+c  Rien
+c  ====
+         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
+      else if(iflag_thermals.eq.0) then
+
+c  Ajustement sec
+c  ==============
+         IF(prt_level>9)WRITE(lunout,*)'ajsec'
+         CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
+         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
+         q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
+      else
+c  Thermiques
+c  ==========
+         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
+     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
+         call calltherm(pdtphys
+     s      ,pplay,paprs,pphi
+     s      ,u_seri,v_seri,t_seri,q_seri
+     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
+     s      ,fm_therm,entr_therm)
+      endif
+c
+c===================================================================
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after dry_adjust'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
+
+
+c-------------------------------------------------------------------------
+c  Caclul des ratqs
+c-------------------------------------------------------------------------
+
+c      print*,'calcul des ratqs'
+c   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
+c   ----------------
+c   on ecrase le tableau ratqsc calcule par clouds_gno
+      if (iflag_cldcon.eq.1) then
+         do k=1,klev
+         do i=1,klon
+            if(ptconv(i,k)) then
+              ratqsc(i,k)=ratqsbas
+     s        +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
+            else
+               ratqsc(i,k)=0.
+            endif
+         enddo
+         enddo
+      endif
+
+c   ratqs stables
+c   -------------
+      do k=1,klev
+cIM RAJOUT boucle do=i
+       do i=1, klon
+cIM         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
+cIM     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.) 
+         ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
+     s   min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.) 
+cIM      print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
+       enddo 
+      enddo
+
+
+c  ratqs final
+c  -----------
+      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
+c   les ratqs sont une conbinaison de ratqss et ratqsc
+c   ratqs final
+c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
+c   relaxation des ratqs
+c         facttemps=exp(-pdtphys/1.e4)
+         facteur=exp(-pdtphys*facttemps)
+         ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
+         ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
+c         print*,'calcul des ratqs fini'
+      else
+c   on ne prend que le ratqs stable pour fisrtilp
+         ratqs(:,:)=ratqss(:,:)
+      endif
+
+
+c
+c Appeler le processus de condensation a grande echelle
+c et le processus de precipitation
+c-------------------------------------------------------------------------
+      CALL fisrtilp(dtime,paprs,pplay,
+     .           t_seri, q_seri,ptconv,ratqs,
+     .           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
+     .           rain_lsc, snow_lsc,
+     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
+     .           frac_impa, frac_nucl,
+     .           prfl, psfl, rhcl)
+
+      WHERE (rain_lsc < 0) rain_lsc = 0.
+      WHERE (snow_lsc < 0) snow_lsc = 0.
+      DO k = 1, klev
+      DO i = 1, klon
+         t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
+         q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
+         ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
+         cldfra(i,k) = rneb(i,k)
+         IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
+      ENDDO
+      ENDDO
+      IF (check) THEN
+         za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
+         WRITE(lunout,*)"apresilp=", za
+         zx_t = 0.0
+         za = 0.0
+         DO i = 1, klon
+            za = za + airephy(i)/FLOAT(klon)
+            zx_t = zx_t + (rain_lsc(i)
+     .                  + snow_lsc(i))*airephy(i)/FLOAT(klon)
+        ENDDO
+         zx_t = zx_t/za*dtime
+         WRITE(lunout,*)"Precip=", zx_t
+      ENDIF
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after fisrt'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, rain_lsc, snow_lsc, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+c
+c-------------------------------------------------------------------
+c  PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
+c-------------------------------------------------------------------
+
+c 1. NUAGES CONVECTIFS
+c
+cIM cf FH
+c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
+       IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
+       snow_tiedtke=0.
+c     print*,'avant calcul de la pseudo precip '
+c     print*,'iflag_cldcon',iflag_cldcon
+       if (iflag_cldcon.eq.-1) then
+          rain_tiedtke=rain_con
+       else
+c       print*,'calcul de la pseudo precip '
+          rain_tiedtke=0.
+c         print*,'calcul de la pseudo precip 0'
+          do k=1,klev
+          do i=1,klon
+             if (d_q_con(i,k).lt.0.) then
+                rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
+     s         *(paprs(i,k)-paprs(i,k+1))/rg
+             endif
+          enddo
+          enddo
+       endif
+c
+c     call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
+c
+
+c Nuages diagnostiques pour Tiedtke
+      CALL diagcld1(paprs,pplay,
+cIM cf FH  .             rain_con,snow_con,ibas_con,itop_con,
+     .             rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
+     .             diafra,dialiq)
+      DO k = 1, klev
+      DO i = 1, klon
+      IF (diafra(i,k).GT.cldfra(i,k)) THEN
+         cldliq(i,k) = dialiq(i,k)
+         cldfra(i,k) = diafra(i,k)
+      ENDIF
+      ENDDO
+      ENDDO
+
+      ELSE IF (iflag_cldcon.eq.3) THEN
+c  On prend pour les nuages convectifs le max du calcul de la
+c  convection et du calcul du pas de temps precedent diminue d'un facteur
+c  facttemps
+c      facttemps=pdtphys/1.e4
+      facteur = pdtphys *facttemps
+      do k=1,klev
+         do i=1,klon
+            rnebcon(i,k)=rnebcon(i,k)*facteur
+            if (rnebcon0(i,k)*clwcon0(i,k).gt.rnebcon(i,k)*clwcon(i,k))
+     s      then
+                rnebcon(i,k)=rnebcon0(i,k)
+                clwcon(i,k)=clwcon0(i,k)
+            endif
+         enddo
+      enddo
+
+c
+cjq - introduce the aerosol direct and first indirect radiative forcings
+cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
+      IF (ok_ade.OR.ok_aie) THEN
+#if defined(CPP_COUPLE) && !defined(INCA_AER)
+         ! Get sulfate aerosol distribution
+         CALL readsulfate(rjourvrai, debut, sulfate)
+         CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
+
+         ! Calculate aerosol optical properties (Olivier Boucher)
+         CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
+     .        tau_ae, piz_ae, cg_ae, aerindex)
+#endif
+#if !defined(CPP_COUPLE)
+         ! Get sulfate aerosol distribution
+         CALL readsulfate(rjourvrai, debut, sulfate)
+         CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
+
+         ! Calculate aerosol optical properties (Olivier Boucher)
+         CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
+     .        tau_ae, piz_ae, cg_ae, aerindex)
+#endif
+cym
+      ELSE
+        tau_ae(:,:,:)=0.0
+	piz_ae(:,:,:)=0.0
+	cg_ae(:,:,:)=0.0
+cym     
+      ENDIF
+
+c
+cIM calcul nuages par le simulateur ISCCP
+c
+#ifdef histISCCP
+      IF (ok_isccp) THEN
+c SD initialisation
+         fq_isccp(:,:,:,:)= 0.
+         totalcldarea(:,:)= 0.
+         meanptop(:,:)    = 0.
+         meantaucld(:,:)  = 0.
+         nbsunlit(:,:,:)  = 0.
+c
+cIM appel simulateur toutes les  NINT(freq_ISCCP/dtime) heures
+       IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
+#include "calcul_simulISCCP.h"
+       ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime)))
+      ENDIF !ok_isccp
+#endif
+
+c   On prend la somme des fractions nuageuses et des contenus en eau
+      cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
+      cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
+
+      ENDIF
+
+c
+c 2. NUAGES STARTIFORMES
+c
+      IF (ok_stratus) THEN
+      CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
+      DO k = 1, klev
+      DO i = 1, klon
+      IF (diafra(i,k).GT.cldfra(i,k)) THEN
+         cldliq(i,k) = dialiq(i,k)
+         cldfra(i,k) = diafra(i,k)
+      ENDIF
+      ENDDO
+      ENDDO
+      ENDIF
+c
+c Precipitation totale
+c
+      DO i = 1, klon
+         rain_fall(i) = rain_con(i) + rain_lsc(i)
+         snow_fall(i) = snow_con(i) + snow_lsc(i)
+      ENDDO
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit="after diagcld"
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
+c
+c Calculer l'humidite relative pour diagnostique
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         zx_t = t_seri(i,k)
+         IF (thermcep) THEN
+            zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
+            zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
+            zx_qs  = MIN(0.5,zx_qs)
+            zcor   = 1./(1.-retv*zx_qs)
+            zx_qs  = zx_qs*zcor
+         ELSE
+           IF (zx_t.LT.t_coup) THEN
+              zx_qs = qsats(zx_t)/pplay(i,k)
+           ELSE
+              zx_qs = qsatl(zx_t)/pplay(i,k)
+           ENDIF
+         ENDIF
+         zx_rh(i,k) = q_seri(i,k)/zx_qs
+         zqsat(i,k)=zx_qs
+      ENDDO
+      ENDDO
+
+cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle 
+c   equivalente a 2m (tpote) pour diagnostique
+c
+      DO i = 1, klon
+       tpot(i)=zt2m(i)*(100000./paprs(i,1))**RKAPPA
+       IF (thermcep) THEN
+        IF(zt2m(i).LT.RTT) then
+	 Lheat=RLSTT
+	ELSE
+	 Lheat=RLVTT
+        ENDIF
+       ELSE
+        IF (zt2m(i).LT.RTT) THEN
+         Lheat=RLSTT
+        ELSE
+	 Lheat=RLVTT
+        ENDIF
+       ENDIF
+       tpote(i) = tpot(i)*     
+     . EXP((Lheat *qsat2m(i))/(RCPD*zt2m(i)))
+      ENDDO
+
+
+#ifdef INCA
+      call VTe(VTphysiq)
+      call VTb(VTinca)
+           calday = FLOAT(julien) + gmtime
+
+#ifdef INCA_AER
+      call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
+     &   ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
+#endif
+
+#ifdef INCAINFO
+           WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
+#endif
+
+           zxsnow_dummy(:) = 0.0
+
+           CALL chemhook_begin (calday,
+#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
+     $                          julien,
+     $                          gmtime,
+#endif
+     $                          pctsrf(1,1),
+     $                          rlat,
+     $                          rlon,
+     $                          airephy,
+     $                          paprs,
+     $                          pplay,
+     $                          ycoefh,
+     $                          pphi,
+     $                          t_seri,
+     $                          u,
+     $                          v,
+     $                          wo,
+     $                          q_seri,
+     $                          zxtsol,
+     $                          zxsnow_dummy,
+     $                          solsw,
+     $                          albsol,
+     $                          rain_fall,
+     $                          snow_fall,
+     $                          itop_con,
+     $                          ibas_con,
+     $                          cldfra,
+     $                          iim,
+     $                          jjm,
+#ifdef INCA_AER
+     $                          tr_seri,
+     $                          ftsol,
+     $                          paprs,
+     $                          cdragh,
+     $                          cdragm,
+     $                          pctsrf,
+     $				pdtphys,
+     $				itap)
+#else
+     $                          tr_seri)     
+#endif        
+
+
+#ifdef INCAINFO
+           WRITE(lunout,*)'OK.'
+#endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
+#endif
+c     
+c Calculer les parametres optiques des nuages et quelques
+c parametres pour diagnostiques:
+c
+      if (ok_newmicro) then
+      CALL newmicro (paprs, pplay,ok_newmicro,
+     .            t_seri, cldliq, cldfra, cldtau, cldemi,
+     .            cldh, cldl, cldm, cldt, cldq,
+     .            flwp, fiwp, flwc, fiwc,
+     e            ok_aie,
+#if defined(CPP_COUPLE) && defined(INCA_AER)  
+     e            ccm(:,:,1), ccm(:,:,2),
+#else
+     e            sulfate, sulfate_pi,
+#endif
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
+      else
+      CALL nuage (paprs, pplay,
+     .            t_seri, cldliq, cldfra, cldtau, cldemi,
+     .            cldh, cldl, cldm, cldt, cldq,
+     e            ok_aie,
+     e            sulfate, sulfate_pi,
+     e            bl95_b0, bl95_b1,
+     s            cldtaupi, re, fl)
+      
+      endif
+c
+c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
+c
+      IF (MOD(itaprad,radpas).EQ.0) THEN
+
+      DO i = 1, klon
+         albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
+     .             + falbe(i,is_lic) * pctsrf(i,is_lic)
+     .             + falbe(i,is_ter) * pctsrf(i,is_ter)
+     .             + falbe(i,is_sic) * pctsrf(i,is_sic)
+         albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
+     .               + falblw(i,is_lic) * pctsrf(i,is_lic)
+     .               + falblw(i,is_ter) * pctsrf(i,is_ter)
+     .               + falblw(i,is_sic) * pctsrf(i,is_sic)
+      ENDDO
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+      
+
+#if defined(CPP_COUPLE) && defined(INCA_AER)
+      CALL radlwsw_inca ! nouveau rayonnement (compatible Arpege-IFS)
+     e            (kdlon,kflev,dist, rmu0, fract, 
+     e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
+     e             wo,
+     e             cldfra, cldemi, cldtau,
+     s             heat,heat0,cool,cool0,radsol,albpla,
+     s             topsw,toplw,solsw,sollw,
+     s             sollwdown,
+     s             topsw0,toplw0,solsw0,sollw0,
+     s             lwdn0, lwdn, lwup0, lwup, 
+     s             swdn0, swdn, swup0, swup,
+     e             ok_ade, ok_aie, ! new for aerosol radiative effects
+     e             tau_inca, piz_inca, cg_inca, ! ="=
+     s             topswad_inca, solswad_inca, ! ="=
+     s             topswad0_inca, solswad0_inca, ! ="=
+     s             topsw_inca, topsw0_inca,
+     s             solsw_inca, solsw0_inca,
+     e             cldtaupi, ! ="=
+     s             topswai_inca, solswai_inca) ! ="=
+#else
+      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
+     e            (dist, rmu0, fract, 
+     e             paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
+     e             wo,
+     e             cldfra, cldemi, cldtau,
+     s             heat,heat0,cool,cool0,radsol,albpla,
+     s             topsw,toplw,solsw,sollw,
+     s             sollwdown,
+     s             topsw0,toplw0,solsw0,sollw0,
+     s             lwdn0, lwdn, lwup0, lwup, 
+     s             swdn0, swdn, swup0, swup,
+     e             ok_ade, ok_aie, ! new for aerosol radiative effects
+     e             tau_ae, piz_ae, cg_ae, ! ="=
+     s             topswad, solswad, ! ="=
+     e             cldtaupi, ! ="=
+     s             topswai, solswai) ! ="=
+#endif
+      itaprad = 0
+      ENDIF
+      itaprad = itaprad + 1
+c
+c Ajouter la tendance des rayonnements (tous les pas)
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         t_seri(i,k) = t_seri(i,k)
+     .               + (heat(i,k)-cool(i,k)) * dtime/86400.
+      ENDDO
+      ENDDO
+c
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+ 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after rad'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , topsw, toplw, solsw, sollw, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+c
+c
+c Calculer l'hydrologie de la surface
+c
+c      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
+c     .            agesno, ftsol,fqsurf,fsnow, ruis)
+c
+
+c
+c Calculer le bilan du sol et la derive de temperature (couplage)
+c
+      DO i = 1, klon
+c         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
+c a la demande de JLD
+         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
+      ENDDO
+c
+cmoddeblott(jan95)
+c Appeler le programme de parametrisation de l'orographie
+c a l'echelle sous-maille:
+c
+      IF (ok_orodr) THEN
+c
+c  selection des points pour lesquels le shema est actif:
+        igwd=0
+        DO i=1,klon
+        itest(i)=0
+c        IF ((zstd(i).gt.10.0)) THEN
+        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
+          itest(i)=1
+          igwd=igwd+1
+          idx(igwd)=i
+        ENDIF
+        ENDDO
+c        igwdim=MAX(1,igwd)
+c
+        CALL drag_noro(klon,klev,dtime,paprs,pplay,
+     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
+     e                   igwd,idx,itest,
+     e                   t_seri, u_seri, v_seri,
+cIM 141004    s                   zulow, zvlow, zustr, zvstr,
+     s                   zulow, zvlow, zustrdr, zvstrdr,
+     s                   d_t_oro, d_u_oro, d_v_oro)
+c
+c  ajout des tendances
+        DO k = 1, klev
+        DO i = 1, klon
+           t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
+           u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
+           v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
+        ENDDO
+        ENDDO
+c
+      ENDIF ! fin de test sur ok_orodr
+c
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+      
+      IF (ok_orolf) THEN
+c
+c  selection des points pour lesquels le shema est actif:
+        igwd=0
+        DO i=1,klon
+        itest(i)=0
+        IF ((zpic(i)-zmea(i)).GT.100.) THEN
+          itest(i)=1
+          igwd=igwd+1
+          idx(igwd)=i
+        ENDIF
+        ENDDO
+c        igwdim=MAX(1,igwd)
+c
+        CALL lift_noro(klon,klev,dtime,paprs,pplay,
+     e                   rlat,zmea,zstd,zpic,
+     e                   itest,
+     e                   t_seri, u_seri, v_seri,
+     s                   zulow, zvlow, zustrli, zvstrli,
+     s                   d_t_lif, d_u_lif, d_v_lif)
+c
+c  ajout des tendances
+        DO k = 1, klev
+        DO i = 1, klon
+           t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
+           u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
+           v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
+        ENDDO
+        ENDDO
+c
+      ENDIF ! fin de test sur ok_orolf
+c
+cIM cf. FLott BEG
+C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
+
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+      DO i = 1, klon
+        zustrph(i)=0.
+        zvstrph(i)=0.
+      ENDDO
+      DO k = 1, klev
+      DO i = 1, klon
+       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
+     c            (paprs(i,k)-paprs(i,k+1))/rg
+       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
+     c            (paprs(i,k)-paprs(i,k+1))/rg
+      ENDDO
+      ENDDO
+c
+cIM calcul composantes axiales du moment angulaire et couple des montagnes
+c
+      IF (is_sequential) THEN
+      
+        CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
+     C                 ra,rg,romega,
+     C                 rlat,rlon,pphis,
+     C                 zustrdr,zustrli,zustrph,
+     C                 zvstrdr,zvstrli,zvstrph,
+     C                 paprs,u,v,
+     C                 aam, torsfc)
+       ENDIF
+cIM cf. FLott END
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after orography'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+      END IF 
+c
+c
+cAA
+cAA Installation de l'interface online-offline pour traceurs
+cAA
+c====================================================================
+c   Calcul  des tendances traceurs
+c====================================================================
+C
+      call phytrac (     rnpb,
+     I                   itap, 
+     I                   julien, 
+     I                   gmtime,
+     I                   debut,
+     I                   lafin,
+     I                   nqmax-2,
+     I                   nlon,
+     I                   nlev,
+     I                   dtime,
+     I                   u,
+     I                   v,
+     I                   t,
+     I                   paprs,
+     I                   pplay,
+     I                   pmfu, 
+     I                   pmfd, 
+     I                   pen_u, 
+     I                   pde_u, 
+     I                   pen_d, 
+     I                   pde_d,
+     I                   ycoefh,
+     I                   fm_therm,
+     I                   entr_therm,
+     I                   yu1,
+     I                   yv1,
+     I                   ftsol,
+     I                   pctsrf,
+     I                   rlat,
+     I                   frac_impa, 
+     I                   frac_nucl,
+     I                   rlon,
+     I                   presnivs,
+     I                   pphis,
+     I                   pphi,
+     I                   albsol,
+     I                   qx(1,1,1), 
+     I                   rhcl,
+     I                   cldfra, 
+     I                   rneb, 
+     I                   diafra, 
+     I                   cldliq, 
+     I                   itop_con,
+     I                   ibas_con,
+     I                   pmflxr,
+     I                   pmflxs,
+     I                   prfl,
+     I                   psfl,
+     I                   da,
+     I                   phi,
+     I                   mp,
+     I                   upwd,
+     I                   dnwd,
+#ifdef INCA
+     I                   flxmass_w,
+#if defined(INCA_AER) && defined(CPP_COUPLE)
+     I                   tau_inca, 
+     I                   piz_inca, 
+     I                   cg_inca,
+     I                   ccm,
+     I                   rfname,
+#endif
+#endif
+     O                   tr_seri)
+
+      IF (offline) THEN
+
+         print*,'Attention on met a 0 les thermiques pour phystoke'
+	 call phystokenc (
+     I                   nlon,nlev,pdtphys,rlon,rlat,
+     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
+     I                   fm_therm,entr_therm,
+     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
+     I                   frac_impa, frac_nucl,
+     I                   pphis,airephy,dtime,itap)
+
+
+      ENDIF
+
+c
+c Calculer le transport de l'eau et de l'energie (diagnostique)
+c
+      CALL transp (paprs,zxtsol,
+     e                   t_seri, q_seri, u_seri, v_seri, zphi,
+     s                   ve, vq, ue, uq)
+c
+cIM global posePB BEG
+      IF(1.EQ.0) THEN
+c
+      CALL transp_lay (paprs,zxtsol,
+     e                   t_seri, q_seri, u_seri, v_seri, zphi,
+     s                   ve_lay, vq_lay, ue_lay, uq_lay)
+c
+      ENDIF !(1.EQ.0) THEN
+cIM global posePB END
+c Accumuler les variables a stocker dans les fichiers histoire:
+c
+c+jld ec_conser
+      DO k = 1, klev
+      DO i = 1, klon
+        ZRCPD = RCPD*(1.0+RVTMP2*q_seri(i,k))
+        d_t_ec(i,k)=0.5/ZRCPD
+     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
+        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
+        d_t_ec(i,k) = d_t_ec(i,k)/dtime
+       END DO 
+      END DO 
+c-jld ec_conser
+cIM
+      IF (ip_ebil_phy.ge.1) THEN 
+        ztit='after physic'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+C     Comme les tendances de la physique sont ajoute dans la dynamique,
+C     on devrait avoir que la variation d'entalpie par la dynamique
+C     est egale a la variation de la physique au pas de temps precedent.
+C     Donc la somme de ces 2 variations devrait etre nulle.
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , topsw, toplw, solsw, sollw, sens
+     e      , evap, rain_fall, snow_fall, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+C
+      d_h_vcol_phy=d_h_vcol
+C
+      END IF 
+C
+c=======================================================================
+c   SORTIES
+c=======================================================================
+
+cIM Interpolation sur les niveaux de pression du NMC
+c   -------------------------------------------------
+c
+#include "calcul_STDlev.h"
+c
+c slp sea level pressure
+      slp(:) = paprs(:,1)*exp(pphis(:)/(RD*t_seri(:,1)))
+c
+ccc prw = eau precipitable
+      DO i = 1, klon
+       prw(i) = 0.
+       DO k = 1, klev
+        prw(i) = prw(i) +
+     .           q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
+       ENDDO
+      ENDDO
+c
+cIM initialisation + calculs divers diag AMIP2
+c
+#include "calcul_divers.h"
+c
+#ifdef INCA
+      call VTe(VTphysiq)
+      call VTb(VTinca)
+#ifdef INCAINFO
+           WRITE(lunout,*)'Appel CHEMHOOK_END ...'
+#endif
+           CALL chemhook_end (calday,
+     $                        dtime,
+     $                        pplay,
+     $                        t_seri,
+     $                        tr_seri,
+     $                        nbtr,
+     $                        paprs,
+     $                        q_seri,
+     $                        annee_ref,
+     $                        day_ini,
+     $                        airephy,
+#ifdef INCA_AER
+     $                        xjour,
+     $                        pphi,
+     $                        pphis,
+     $                        zx_rh)
+#else
+     $                        xjour)
+#endif
+#ifdef INCAINFO
+           WRITE(lunout,*)'OK.'
+#endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
+#endif
+
+c=============================================================
+c
+c Convertir les incrementations en tendances
+c
+      if (mydebug) then
+        call writefield_phy('u_seri',u_seri,llm)
+        call writefield_phy('v_seri',v_seri,llm)
+        call writefield_phy('t_seri',t_seri,llm)
+	call writefield_phy('q_seri',q_seri,llm)
+      endif
+
+      DO k = 1, klev
+      DO i = 1, klon
+         d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
+         d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
+         d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
+         d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
+         d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
+      ENDDO
+      ENDDO
+c
+      IF (nqmax.GE.3) THEN
+      DO iq = 3, nqmax
+      DO  k = 1, klev
+      DO  i = 1, klon
+         d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
+      ENDDO
+      ENDDO
+      ENDDO
+      ENDIF
+c
+cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
+cIM global posePB#include "write_bilKP_ins.h"
+cIM global posePB#include "write_bilKP_ave.h"
+c
+c Sauvegarder les valeurs de t et q a la fin de la physique:
+c
+      DO k = 1, klev
+      DO i = 1, klon
+         t_ancien(i,k) = t_seri(i,k)
+         q_ancien(i,k) = q_seri(i,k)
+      ENDDO
+      ENDDO
+c
+c 22.03.04 BEG
+c=============================================================
+c   Ecriture des sorties
+c=============================================================
+#ifdef CPP_IOIPSL
+ 
+c Recupere des varibles calcule dans differents modules
+c pour ecriture dans histxxx.nc 
+
+      ! Get some variables from module mod_fonte_neige
+      CALL fonte_neige_get_vars(pctsrf, 
+     .     zxfqcalving, zxfqfonte, zxffonte)
+
+      IF (ocean == 'slab') THEN
+         ! Get some variables from module oceanslab
+         CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
+      ELSEIF (ocean == 'couple') THEN
+         ! Get some variables from module oceancpl
+         CALL ocean_cpl_get_vars(fluxo, fluxg)
+      ELSE 
+         ! Get some variables from module oceanforced
+         CALL ocean_forced_get_vars(fluxo, fluxg)
+      ENDIF
+
+#ifdef histhf
+#include "write_histhf.h"
+#endif
+
+#ifdef histday
+#include "write_histday.h"
+#endif
+
+#ifdef histmth
+#include "write_histmth.h"
+#endif
+
+#ifdef histins
+#include "write_histins.h"
+#endif
+
+#ifdef histISCCP
+#include "write_histISCCP.h"
+#endif
+
+#ifdef histmthNMC
+#include "write_histmthNMC.h"
+#endif
+
+#include "write_histday_seri.h"
+
+#include "write_paramLMDZ_phy.h"
+
+#endif
+
+c 22.03.04 END
+c
+c====================================================================
+c Si c'est la fin, il faut conserver l'etat de redemarrage
+c====================================================================
+c
+      
+
+      IF (lafin) THEN
+         itau_phy = itau_phy + itap
+         CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
+     .      rlat, rlon, pctsrf, ftsol,
+     .      falbe,falblw, rain_fall, 
+     .      snow_fall,
+     .      solsw, sollw,
+     .      radsol,
+     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
+     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
+      ENDIF
+      
+
+      RETURN
+      END
+      FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
+      IMPLICIT none
+c
+c Calculer et imprimer l'eau totale. A utiliser pour verifier
+c la conservation de l'eau
+c
+#include "YOMCST.h"
+      INTEGER klon,klev
+      REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
+      REAL aire(klon)
+      REAL qtotal, zx, qcheck
+      INTEGER i, k
+c
+      zx = 0.0
+      DO i = 1, klon
+         zx = zx + aire(i)
+      ENDDO
+      qtotal = 0.0
+      DO k = 1, klev
+      DO i = 1, klon
+         qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
+     .                     *(paprs(i,k)-paprs(i,k+1))/RG
+      ENDDO
+      ENDDO
+c
+      qcheck = qtotal/zx
+c
+      RETURN
+      END
+      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
+      IMPLICIT none
+c
+c Tranformer une variable de la grille physique a
+c la grille d'ecriture
+c
+      INTEGER nfield,nlon,iim,jjmp1, jjm
+      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
+c
+      INTEGER i, n, ig
+c
+      jjm = jjmp1 - 1
+      DO n = 1, nfield
+         DO i=1,iim
+            ecrit(i,n) = fi(1,n)
+            ecrit(i+jjm*iim,n) = fi(nlon,n)
+         ENDDO
+         DO ig = 1, nlon - 2
+           ecrit(iim+ig,n) = fi(1+ig,n)
+         ENDDO
+      ENDDO
+      RETURN
+      END
