Index: LMDZ4/branches/V3_test/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/V3_test/libf/phylmd/physiq.F	(revision 703)
+++ LMDZ4/branches/V3_test/libf/phylmd/physiq.F	(revision 704)
@@ -3,4 +3,6 @@
 !
 c
+c#define IO_DEBUG
+
       SUBROUTINE physiq (nlon,nlev,nqmax,
      .            debut,lafin,rjourvrai,gmtime,pdtphys,
@@ -17,7 +19,12 @@
 
       USE ioipsl
-      USE histcom
+      USE comgeomphy
+      USE write_field_phy
+      USE dimphy
+      USE iophy
+      USE misc_mod, mydebug=>debug
+      USE vampir
 #ifdef INCA
-      USE chemshut
+cym      USE chemshut
       USE species_names
 #ifdef INCA_CH4
@@ -84,5 +91,7 @@
       integer jjmp1
       parameter (jjmp1=jjm+1-1/jjm)
-#include "dimphy.h"
+      integer iip1
+      parameter (iip1=iim+1)
+cym#include "dimphy.h"
 #include "regdim.h"
 #include "indicesol.h"
@@ -92,5 +101,5 @@
 #include "logic.h"
 #include "temps.h"
-#include "comgeomphy.h"
+cym#include "comgeomphy.h"
 #include "advtrac.h"
 #include "iniprint.h"
@@ -110,4 +119,5 @@
 #include "oasis.h"
       INTEGER,SAVE :: npas, nexca
+c$OMP THREADPRIVATE(npas, nexca)
       logical rnpb
 #ifdef INCA
@@ -119,9 +129,10 @@
       character*6 ocean
       SAVE ocean
-
+c$OMP THREADPRIVATE(ocean)
 c      parameter (ocean = 'force ')
 c     parameter (ocean = 'couple')
       logical ok_ocean
       SAVE ok_ocean
+c$OMP THREADPRIVATE(ok_ocean)
 c
 cIM "slab" ocean
@@ -144,4 +155,5 @@
       logical ok_veget
       save ok_veget
+c$OMP THREADPRIVATE(ok_veget)
 c     parameter (ok_veget = .true.)
 c      parameter (ok_veget = .false.)
@@ -163,12 +175,15 @@
       LOGICAL ok_journe ! sortir le fichier journalier
       save ok_journe
+c$OMP THREADPRIVATE(ok_journe)
 c      PARAMETER (ok_journe=.true.)
 c
       LOGICAL ok_mensuel ! sortir le fichier mensuel
       save ok_mensuel
+c$OMP THREADPRIVATE(ok_mensuel)
 c      PARAMETER (ok_mensuel=.true.)
 c
       LOGICAL ok_instan ! sortir le fichier instantane
       save ok_instan
+c$OMP THREADPRIVATE(ok_instan)
 c      PARAMETER (ok_instan=.true.)
 c
@@ -179,6 +194,7 @@
       REAL fm_therm(klon,klev+1)
       REAL entr_therm(klon,klev)
-      real q2(klon,klev+1,nbsrf)
-      save q2
+      real,allocatable,save :: q2(:,:,:)
+c$OMP THREADPRIVATE(q2)
+cym      save q2
 c======================================================================
 c
@@ -215,9 +231,10 @@
       REAL qx(klon,klev,nqmax)
 
-      REAL t_ancien(klon,klev), q_ancien(klon,klev)
-      SAVE t_ancien, q_ancien
+      REAL,allocatable,save :: t_ancien(:,:), q_ancien(:,:)
+c$OMP THREADPRIVATE(t_ancien, q_ancien)
+cym      SAVE 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)
@@ -241,6 +258,8 @@
       CHARACTER*3 ctetaSTD(nbteta)
       DATA ctetaSTD/'350','380','405'/
+c$OMP THREADPRIVATE(ctetaSTD)
       REAL rtetaSTD(nbteta)
       DATA rtetaSTD/350., 380., 405./
+c$OMP THREADPRIVATE(rtetaSTD)     
 c
       REAL PVteta(klon,nbteta)
@@ -249,26 +268,32 @@
 cMI Amip2 PV a theta constante
 
-      INTEGER klevp1, klevm1
-      PARAMETER(klevp1=klev+1,klevm1=klev-1)
-#include "raddim.h"
-c
-      REAL swdn0(klon,klevp1), swdn(klon,klevp1)
-      REAL swup0(klon,klevp1), swup(klon,klevp1)
-      SAVE swdn0 , swdn, swup0, swup
-c
-      REAL SWdn200clr(klon), SWdn200(klon)
-      REAL SWup200clr(klon), SWup200(klon)
-      SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
-c
-      REAL lwdn0(klon,klevp1), lwdn(klon,klevp1)
-      REAL lwup0(klon,klevp1), lwup(klon,klevp1)
-      SAVE lwdn0 , lwdn, lwup0, lwup 
-c
-      REAL LWdn200clr(klon), LWdn200(klon)
-      REAL LWup200clr(klon), LWup200(klon)
-      SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
-c
-      REAL LWdnTOA(klon), LWdnTOAclr(klon)
-      SAVE LWdnTOA, LWdnTOAclr
+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)
+cym      SAVE swdn0 , swdn, swup0, swup
+c
+      REAL,allocatable,save :: SWdn200clr(:), SWdn200(:)
+      REAL,allocatable,save :: SWup200clr(:), SWup200(:)
+c$OMP THREADPRIVATE(SWdn200clr, SWdn200, SWup200clr, SWup200)
+cym      SAVE SWdn200clr, SWdn200, SWup200clr, SWup200
+c
+      REAL,allocatable,save :: lwdn0(:,:), lwdn(:,:)
+      REAL,allocatable,save :: lwup0(:,:), lwup(:,:)
+c$OMP THREADPRIVATE(lwdn0 , lwdn, lwup0, lwup)
+cym      SAVE lwdn0 , lwdn, lwup0, lwup 
+c
+      REAL,allocatable,save :: LWdn200clr(:), LWdn200(:)
+      REAL,allocatable,save :: LWup200clr(:), LWup200(:)
+c$OMP THREADPRIVATE(LWdn200clr, LWdn200, LWup200clr, LWup200)
+cym      SAVE LWdn200clr, LWdn200, LWup200clr, LWup200
+c
+      REAL,allocatable,save :: LWdnTOA(:), LWdnTOAclr(:)
+c$OMP THREADPRIVATE(LWdnTOA, LWdnTOAclr)
+cym      SAVE LWdnTOA, LWdnTOAclr
 c
 cIM Amip2
@@ -281,8 +306,10 @@
      .60000., 50000., 40000., 30000., 25000., 20000.,
      .15000., 10000., 7000., 5000., 3000., 2000., 1000./
+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  '/
+c$OMP THREADPRIVATE(clevSTD)
 c
       CHARACTER*3 bb2
@@ -298,15 +325,18 @@
       PARAMETER(nout=3) !nout=1 : day; =2 : mth; =3 : NMC
 c
-      REAL tsumSTD(klon,nlevSTD,nout)
-      REAL usumSTD(klon,nlevSTD,nout), vsumSTD(klon,nlevSTD,nout)
-      REAL wsumSTD(klon,nlevSTD,nout), phisumSTD(klon,nlevSTD,nout)
-      REAL qsumSTD(klon,nlevSTD,nout), rhsumSTD(klon,nlevSTD,nout)
-c
-      SAVE tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD, 
-     .     qsumSTD, rhsumSTD
+      REAL,SAVE,ALLOCATABLE :: tsumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: usumSTD(:,:,:), vsumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: wsumSTD(:,:,:), phisumSTD(:,:,:)
+      REAL,SAVE,ALLOCATABLE :: qsumSTD(:,:,:), rhsumSTD(:,:,:)
+c
+cym      SAVE tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD, 
+cym     .     qsumSTD, rhsumSTD
+c$OMP THREADPRIVATE(tsumSTD, usumSTD, vsumSTD, wsumSTD, phisumSTD) 
+c$OMP THREADPRIVATE(qsumSTD, rhsumSTD)
 c
       logical oknondef(klon,nlevSTD,nout)
-      real tnondef(klon,nlevSTD,nout) 
-      save tnondef
+      real,SAVE,ALLOCATABLE :: tnondef(:,:,:) 
+c$OMP THREADPRIVATE(tnondef)
+cym      save tnondef
 c
 c les produits uvSTD, vqSTD, .., T2STD sont calcules
@@ -319,8 +349,8 @@
       real wqSTD(klon,nlevSTD)
 c
-      real uvsumSTD(klon,nlevSTD,nout)
-      real vqsumSTD(klon,nlevSTD,nout)
-      real vTsumSTD(klon,nlevSTD,nout)
-      real wqsumSTD(klon,nlevSTD,nout)
+      real,save,allocatable :: uvsumSTD(:,:,:)
+      real,save,allocatable :: vqsumSTD(:,:,:)
+      real,save,allocatable :: vTsumSTD(:,:,:)
+      real,save,allocatable :: wqsumSTD(:,:,:)
 c
       real vphiSTD(klon,nlevSTD)
@@ -330,12 +360,15 @@
       real T2STD(klon,nlevSTD)
 c
-      real vphisumSTD(klon,nlevSTD,nout)
-      real wTsumSTD(klon,nlevSTD,nout)
-      real u2sumSTD(klon,nlevSTD,nout)
-      real v2sumSTD(klon,nlevSTD,nout)
-      real T2sumSTD(klon,nlevSTD,nout)
-c
-      SAVE uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD
-      SAVE vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD
+      real,save,allocatable :: vphisumSTD(:,:,:)
+      real,save,allocatable :: wTsumSTD(:,:,:)
+      real,save,allocatable :: u2sumSTD(:,:,:)
+      real,save,allocatable :: v2sumSTD(:,:,:)
+      real,save,allocatable :: T2sumSTD(:,:,:)
+c
+cym      SAVE uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD
+cym      SAVE vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD
+c$OMP THREADPRIVATE(uvsumSTD, vqsumSTD, vTsumSTD, wqsumSTD)
+c$OMP THREADPRIVATE(vphisumSTD, wTsumSTD, u2sumSTD, v2sumSTD, T2sumSTD)
+
 cMI Amip2
 c
@@ -355,5 +388,5 @@
       REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
 
-      INTEGER kp1
+      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)
@@ -369,6 +402,6 @@
 cv3.4
       INTEGER debug, debugcol
-      INTEGER npoints
-      PARAMETER(npoints=klon) 
+cym      INTEGER npoints
+cym      PARAMETER(npoints=klon) 
 c
       INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
@@ -390,9 +423,12 @@
       DATA ifreq_isccp/3/
       SAVE ifreq_isccp
+c$OMP THREAPRIVATE(ifreq_isccp)
       CHARACTER*5 typinout(napisccp)
       DATA typinout/'i3od'/
+c$OMP THREAPRIVATE(typinout)
 cIM verif boxptop BEG
       CHARACTER*1 verticaxe(napisccp)
       DATA verticaxe/'1'/ 
+c$OMP THREAPRIVATE(verticaxe)
 cIM verif boxptop END
       INTEGER nvlev(napisccp)
@@ -400,7 +436,9 @@
       REAL t1, aa
       REAL seed_re(klon,napisccp)
-      INTEGER seed_old(klon,napisccp)
-      SAVE seed_old
-      INTEGER iphy(iim,jjmp1)
+      INTEGER,ALLOCATABLE,SAVE :: seed_old(:,:)
+cym      SAVE seed_old
+c$OMP THREADPRIVATE(seed_old)
+cym !!!! A voir plus tard 
+cym      INTEGER iphy(iim,jjmp1)
 cIM parametres ISCCP END
 c
@@ -481,6 +519,7 @@
 
       INTEGER linv
-      INTEGER pct_ocean(klon,nbregdyn)
-      SAVE pct_ocean
+      INTEGER,ALLOCATABLE,SAVE :: pct_ocean(:,:)
+c$OMP THREADPRIVATE(pct_ocean)
+cym      SAVE pct_ocean
  
 c sorties ISCCP
@@ -491,5 +530,5 @@
 c     save ok_isccp, ecrit_isccp, nid_isccp        
       save nid_isccp        
-
+c$OMP THREADPRIVATE(ok_isccp, ecrit_isccp, nid_isccp,nid_isccp)
 #undef histISCCP
 #define histISCCP
@@ -505,5 +544,5 @@
 cIM bad 151205     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
       DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
-
+c$OMP THREADPRIVATE(zx_tau,zx_pc)
 c cldtopres pression au sommet des nuages
       REAL cldtopres(lmaxm1), cldtopres3(lmax3)
@@ -511,4 +550,5 @@
       DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
       DATA cldtopres3/440., 680., 1000./
+c$OMP THREADPRIVATE(cldtopres,cldtopres3)
 cIM 051005 BEG
       REAL tmp_his1_3d(iwmax,kmaxm1,lmaxm1,nbregdyn,napisccp) 
@@ -525,4 +565,5 @@
       CHARACTER *3 pclev(lmaxm1)
       DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
+c$OMP THREADPRIVATE(taulev,pclev)
 c
 c cnameisccp
@@ -578,4 +619,5 @@
      .                'pc= 680-800hPa, tau> 60.',
      .                'pc= 800-1000hPa, tau> 60.'/
+c$OMP THREADPRIVATE(cnameisccp)
 c
 c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
@@ -591,9 +633,9 @@
       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 ?
 
 #undef histhf
-#define histhf
+c#define histhf
 #ifdef histhf
 cIM 130904   data ok_hf,ecrit_hf/.true.,0.25/
@@ -602,5 +644,4 @@
       data ok_hf/.false./
 #endif
-
       INTEGER        longcles
       PARAMETER    ( longcles = 20 )
@@ -611,4 +652,5 @@
       REAL xjour
       SAVE xjour
+c$OMP THREADPRIVATE(xjour)
 c
 c
@@ -617,19 +659,25 @@
       REAL dtime
       SAVE dtime                  ! pas temporel de la physique
+c$OMP THREADPRIVATE(dtime)
 c
       INTEGER radpas
       SAVE radpas                 ! frequence d'appel rayonnement
-c
-      REAL radsol(klon)
-      SAVE radsol               ! bilan radiatif au sol calcule par code radiatif
-c
-      REAL rlat(klon)
-      SAVE rlat                   ! latitude pour chaque point
-c
-      REAL rlon(klon)
-      SAVE rlon                   ! longitude pour chaque point
-c
-      REAL rlonPOS(klon)
-      SAVE rlonPOS                ! longitudes > 0. pour chaque point
+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
@@ -638,4 +686,5 @@
       INTEGER itap
       SAVE itap                   ! compteur pour la physique
+c$OMP THREADPRIVATE(itap)
 c
       REAL co2_ppm_etat0
@@ -645,34 +694,46 @@
       real slp(klon) ! sea level pressure
 
-      REAL ftsol(klon,nbsrf)
-      SAVE ftsol                  ! temperature du sol
+      REAL,allocatable,save :: ftsol(:,:)
+c$OMP THREADPRIVATE(ftsol)
+cym      SAVE ftsol                  ! temperature du sol
+
 cIM
-      REAL newsst(klon) !temperature de l'ocean
-      SAVE newsst
-c
-      REAL ftsoil(klon,nsoilmx,nbsrf)
-      SAVE ftsoil                 ! temperature dans le sol
-c
-      REAL fevap(klon,nbsrf)
-      SAVE fevap                 ! evaporation
-      REAL fluxlat(klon,nbsrf)
-      SAVE fluxlat
-c
-      REAL deltat(klon)
-      SAVE deltat                 ! ecart avec la SST de reference
-c
-      REAL fqsurf(klon,nbsrf)
-      SAVE fqsurf                 ! humidite de l'air au contact de la surface
-c
-      REAL qsol(klon)
-      SAVE qsol                  ! hauteur d'eau dans le sol
-c
-      REAL fsnow(klon,nbsrf)
-      SAVE fsnow                  ! epaisseur neigeuse
-c
-      REAL falbe(klon,nbsrf)
-      SAVE falbe                  ! albedo par type de surface
-      REAL falblw(klon,nbsrf)
-      SAVE falblw                 ! albedo par type de surface
+      REAL,SAVE,ALLOCATABLE :: newsst(:) !temperature de l'ocean
+c$OMP THREADPRIVATE(newsst)
+cym     SAVE newsst
+c
+      REAL,allocatable,save :: ftsoil(:,:,:)
+c$OMP THREADPRIVATE(ftsoil)
+cym      SAVE ftsoil                 ! temperature dans le sol
+c
+      REAL,allocatable,save :: fevap(:,:)
+c$OMP THREADPRIVATE(fevap)
+cym      SAVE fevap                 ! evaporation
+      REAL,allocatable,save :: fluxlat(:,:)
+c$OMP THREADPRIVATE(fluxlat)
+cym      SAVE fluxlat
+c
+      REAL,allocatable,save :: deltat(:)
+c$OMP THREADPRIVATE(deltat)
+cym      SAVE deltat                 ! ecart avec la SST de reference
+c
+      REAL,allocatable,save :: fqsurf(:,:)
+c$OMP THREADPRIVATE(fqsurf)
+cym      SAVE fqsurf                 ! humidite de l'air au contact de la surface
+c
+      REAL,allocatable,save :: qsol(:)
+c$OMP THREADPRIVATE(qsol)
+cym      SAVE qsol                  ! hauteur d'eau dans le sol
+c
+      REAL,allocatable,save :: fsnow(:,:)
+c$OMP THREADPRIVATE(fsnow)
+cym      SAVE fsnow                  ! epaisseur neigeuse
+c
+      REAL,allocatable,save :: falbe(:,:)
+c$OMP THREADPRIVATE(falbe)
+cym      SAVE falbe                  ! albedo par type de surface
+      REAL,allocatable,save :: falblw(:,:)
+c$OMP THREADPRIVATE(falblw)      
+cym      SAVE falblw                 ! albedo par type de surface
 
 c
@@ -680,73 +741,94 @@
 c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
 c
-      REAL zmea(klon)
-      SAVE zmea                   ! orographie moyenne
-c
-      REAL zstd(klon)
-      SAVE zstd                   ! deviation standard de l'OESM
-c
-      REAL zsig(klon)
-      SAVE zsig                   ! pente de l'OESM
-c
-      REAL zgam(klon)
-      save zgam                   ! anisotropie de l'OESM
-c
-      REAL zthe(klon)
-      SAVE zthe                   ! orientation de l'OESM
-c
-      REAL zpic(klon)
-      SAVE zpic                   ! Maximum de l'OESM
-c
-      REAL zval(klon)
-      SAVE zval                   ! Minimum de l'OESM
-c
-      REAL rugoro(klon)
-      SAVE rugoro                 ! longueur de rugosite de l'OESM
+      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 zuthe(klon),zvthe(klon)
-      SAVE zuthe
-      SAVE zvthe
+      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)
-      SAVE agesno                 ! age de la neige
-c
-      REAL alb_neig(klon)
-      SAVE alb_neig               ! albedo de la neige
-c
-      REAL run_off_lic_0(klon)
-      SAVE run_off_lic_0
+      REAL,allocatable,save :: agesno(:,:)
+c$OMP THREADPRIVATE(agesno)
+cym      SAVE agesno                 ! age de la neige
+c
+      REAL,allocatable,save :: alb_neig(:)
+c$OMP THREADPRIVATE(alb_neig)
+cym      SAVE alb_neig               ! albedo de la neige
+c
+      REAL,allocatable,save :: run_off_lic_0(:)
+c$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 ema_workcbmf(klon)   ! cloud base mass flux
-      SAVE ema_workcbmf
-
-      REAL ema_cbmf(klon)       ! cloud base mass flux
-      SAVE ema_cbmf
-
-      REAL ema_pcb(klon)        ! cloud base pressure
-      SAVE ema_pcb
-
-      REAL ema_pct(klon)        ! cloud top pressure
-      SAVE ema_pct
+      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
-
-      REAL Ma(klon,klev)        ! undilute upward mass flux
-      SAVE Ma
-      REAL qcondc(klon,klev)    ! in-cld water content from convect
-      SAVE qcondc 
-      REAL ema_work1(klon, klev), ema_work2(klon, klev)
-      SAVE ema_work1, ema_work2
+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 wd(klon) ! sb
-      SAVE wd       ! sb
+      REAL,allocatable,save :: wd(:) ! sb
+c$OMP THREADPRIVATE(wd)
+cym      SAVE wd       ! sb
 
 c Variables locales pour la couche limite (al1):
@@ -766,19 +848,26 @@
       REAL yu1(klon)            ! vents dans la premiere couche U
       REAL yv1(klon)            ! vents dans la premiere couche V
-      REAL ffonte(klon,nbsrf)    !Flux thermique utilise pour fondre la neige
-      REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface 
+      REAL,SAVE,ALLOCATABLE :: ffonte(:,:)    !Flux thermique utilise pour fondre la neige
+c$OMP THREAPRIVATE(ffonte)
+      REAL,SAVE,ALLOCATABLE :: fqcalving(:,:) !Flux d'eau "perdu" par la surface 
+c$OMP THREAPRIVATE(fqcalving)
+      REAL,SAVE,ALLOCATABLE :: fqfonte(:,:)  !Quantite d'eau de fonte des glaciers
+c$OMP THREAPRIVATE(fqcalving)
 c                               !et necessaire pour limiter la
 c                               !hauteur de neige, en kg/m2/s
-      REAL zxffonte(klon), zxfqcalving(klon)
-
-c$$$      LOGICAL offline           ! Controle du stockage ds "physique"
-c$$$      PARAMETER (offline=.false.)
-c$$$      INTEGER physid
-      REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
-      save pfrac_impa
-      REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
-      save pfrac_nucl
-      REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
-      save pfrac_1nucl
+      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)
@@ -789,18 +878,24 @@
 
 cAA
-      REAL rain_fall(klon) ! pluie
-      REAL snow_fall(klon) ! neige
-      save snow_fall, rain_fall
+      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 total_rain(klon), nday_rain(klon)
-      save nday_rain
-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 dlw(klon)    ! derivee infra rouge
+      REAL,allocatable,save :: dlw(:)    ! derivee infra rouge
+c$OMP THREADPRIVATE(dlw)
 cym
-      SAVE dlw
+cym      SAVE dlw
 cym
       REAL bils(klon) ! bilan de chaleur au sol
@@ -809,6 +904,7 @@
       REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
 C                             ! type de sous-surface et pondere par la fraction
-      REAL fder(klon) ! Derive de flux (sensible et latente) 
-      save fder
+      REAL,allocatable,save :: fder(:) ! Derive de flux (sensible et latente)
+c$OMP THREADPRIVATE(fder) 
+cym      save fder
       REAL ve(klon) ! integr. verticale du transport meri. de l'energie
       REAL vq(klon) ! integr. verticale du transport meri. de l'eau
@@ -816,6 +912,7 @@
       REAL uq(klon) ! integr. verticale du transport zonal de l'eau
 c
-      REAL frugs(klon,nbsrf) ! longueur de rugosite
-      save frugs
+      REAL,allocatable,save :: frugs(:,:) ! longueur de rugosite
+c$OMP THREADPRIVATE(frugs) 
+cym      save frugs
       REAL zxrugs(klon) ! longueur de rugosite
 c
@@ -826,17 +923,27 @@
       INTEGER lmt_pas
       SAVE lmt_pas                ! frequence de mise a jour
-      REAL pctsrf(klon,nbsrf)
+c$OMP THREADPRIVATE(lmt_pas) 
+      REAL,allocatable,save :: pctsrf(:,:)
+c$OMP THREADPRIVATE(pctsrf) 
 cIM
       REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
-      REAL paire_ter(klon)        !surfaces terre 
-c
-      SAVE pctsrf                 ! sous-fraction du sol
-      REAL albsol(klon)
-      SAVE albsol                 ! albedo du sol total
-      REAL albsollw(klon)
-      SAVE albsollw                 ! albedo du sol total
-
-      REAL wo(klon,klev)
-      SAVE wo                     ! ozone
+
+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
@@ -885,11 +992,9 @@
 c Variables locales
 c
-      real clwcon(klon,klev),rnebcon(klon,klev)
-      real clwcon0(klon,klev),rnebcon0(klon,klev)
-cIM cf. AM 081204 BEG
-      real clwcon0th(klon,klev),rnebcon0th(klon,klev)
-cIM cf. AM 081204 END
-      save rnebcon, clwcon
-
+      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
@@ -911,27 +1016,39 @@
       REAL zxfluxv(klon, klev)
 CXXX
-      REAL heat(klon,klev)    ! chauffage solaire
-      REAL heat0(klon,klev)   ! chauffage solaire ciel clair
-      REAL cool(klon,klev)    ! refroidissement infrarouge
-      REAL cool0(klon,klev)   ! refroidissement infrarouge ciel clair
-      REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
-      real sollwdown(klon)    ! downward LW flux at surface
+      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 sollwdownclr(klon)    ! downward CS LW flux at surface 
-      real toplwdown(klon)       ! downward CS LW flux at TOA
-      real toplwdownclr(klon)    ! downward CS LW flux at TOA
+      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 topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
-      REAL albpla(klon)
+      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
-      SAVE  heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
-      SAVE  sollwdownclr, toplwdown, toplwdownclr
-      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
+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)
@@ -961,19 +1078,30 @@
 c
       REAL zphi(klon,klev)
-      REAL zx_relief(iim,jjmp1)
-      REAL zx_aire(iim,jjmp1)
+cym A voir plus tard !!
+cym      REAL zx_relief(iim,jjmp1)
+cym      REAL zx_aire(iim,jjmp1)
 c
 cIM cf. AM Variables locales pour la CLA (hbtm2)
 c
-      REAL pblh(klon, nbsrf)           ! Hauteur de couche limite
-      REAL plcl(klon, nbsrf)           ! Niveau de condensation de la CLA
-      REAL capCL(klon, nbsrf)          ! CAPE de couche limite
-      REAL oliqCL(klon, nbsrf)          ! eau_liqu integree de couche limite
-      REAL cteiCL(klon, nbsrf)          ! cloud top instab. crit. couche limite
-      REAL pblt(klon, nbsrf)          ! T a la Hauteur de couche limite
-      REAL therm(klon, nbsrf)
-      REAL trmb1(klon, nbsrf)          ! deep_cape
-      REAL trmb2(klon, nbsrf)          ! inhibition 
-      REAL trmb3(klon, nbsrf)          ! Point Omega
+      REAL,SAVE,ALLOCATABLE :: pblh(:, :)           ! Hauteur de couche limite
+c$OMP THREAPRIVATE(pblh)
+      REAL,SAVE,ALLOCATABLE :: plcl(:, :)           ! Niveau de condensation de la CLA
+c$OMP THREAPRIVATE(plcl)
+      REAL,SAVE,ALLOCATABLE :: capCL(:, :)          ! CAPE de couche limite
+c$OMP THREAPRIVATE(capCL)
+      REAL,SAVE,ALLOCATABLE :: oliqCL(:, :)          ! eau_liqu integree de couche limite
+c$OMP THREAPRIVATE(oliqCL)
+      REAL,SAVE,ALLOCATABLE :: cteiCL(:, :)          ! cloud top instab. crit. couche limite
+c$OMP THREAPRIVATE(cteiCL)
+      REAL,SAVE,ALLOCATABLE :: pblt(:, :)          ! T a la Hauteur de couche limite
+c$OMP THREAPRIVATE(pblt)
+      REAL,SAVE,ALLOCATABLE :: therm(:, :)
+c$OMP TREADPRIVATE(therm)
+      REAL,SAVE,ALLOCATABLE :: trmb1(:, :)          ! deep_cape
+c$OMP TREADPRIVATE(trmb1)
+      REAL,SAVE,ALLOCATABLE :: trmb2(:, :)          ! inhibition 
+c$OMP TREADPRIVATE(trmb2)
+      REAL,SAVE,ALLOCATABLE :: trmb3(:, :)          ! Point Omega
+c$OMP TREADPRIVATE(trmb3)
 c Grdeurs de sorties
       REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
@@ -988,12 +1116,15 @@
       REAL dnwd0(klon,klev)     ! unsaturated downdraft mass flux
       REAL tvp(klon,klev)       ! virtual temp of lifted parcel
-      REAL cape(klon)           ! CAPE
-      SAVE cape
+      REAL,allocatable,save :: cape(:)           ! CAPE
+c$OMP THREADPRIVATE(cape)
+cym      SAVE cape
       CHARACTER*40 capemaxcels  !max(CAPE)
 
-      REAL pbase(klon)          ! cloud base pressure
-      SAVE pbase
-      REAL bbase(klon)          ! cloud base buoyancy
-      SAVE bbase
+      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
@@ -1018,5 +1149,6 @@
 c vdf: couche limite (Vertical DiFfusion)
       REAL d_t_con(klon,klev),d_q_con(klon,klev)
-      REAL d_u_con(klon,klev),d_v_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)
@@ -1028,6 +1160,7 @@
 *********************************************************
 *     declarations
-      real zqasc(klon,klev)
-      save zqasc
+      real,save,allocatable :: zqasc(:,:)
+c$OMP THREADPRIVATE(zqasc)
+cym      save zqasc
       
 *********************************************************
@@ -1041,10 +1174,15 @@
       REAL prfl(klon,klev+1), psfl(klon,klev+1)
 c
-      INTEGER ibas_con(klon), itop_con(klon)
+      INTEGER,allocatable,save :: ibas_con(:), itop_con(:)
+c$OMP THREADPRIVATE(ibas_con,itop_con)
 cym
-      SAVE ibas_con,itop_con
+cym      SAVE ibas_con,itop_con
 cym
-      REAL rain_con(klon), rain_lsc(klon)
-      REAL snow_con(klon), snow_lsc(klon)
+      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
@@ -1058,7 +1196,11 @@
       REAL d_u_oli(klon,klev), d_v_oli(klon,klev) !tendances dues a oro et lif 
 
-      REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
+      REAL,allocatable,save :: ratqs(:,:)
+c$OMP THREADPRIVATE(ratqs)
+      REAL ratqss(klon,klev),ratqsc(klon,klev)
       real ratqsbas,ratqshaut
-      save ratqsbas,ratqshaut, ratqs
+cym      save ratqsbas,ratqshaut, ratqs
+      save ratqsbas,ratqshaut
+c$OMP THREADPRIVATE(ratqsbas,ratqshaut)
       real zpt_conv(klon,klev)
 
@@ -1068,10 +1210,12 @@
       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
@@ -1087,4 +1231,5 @@
       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
@@ -1123,4 +1268,5 @@
       REAL tabcntr0( length       )
 c
+
       INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
 cIM
@@ -1138,4 +1284,5 @@
       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)
@@ -1147,5 +1294,6 @@
 cIM 
       REAL airetot, pi
-      REAL zm_wo(jjmp1, klev)
+cym A voir plus tard
+cym      REAL zm_wo(jjmp1, klev)
 cIM AMIP2 END
 c
@@ -1153,5 +1301,6 @@
       REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
 #ifdef histmthNMC
-      REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
+cym   A voir plus tard !!!!
+cym      REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
       REAL zx_tmp_fiNC(klon,nlevSTD) 
 #endif
@@ -1164,8 +1313,10 @@
       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, 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.
@@ -1182,7 +1333,9 @@
       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
@@ -1208,11 +1361,20 @@
       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
@@ -1220,6 +1382,8 @@
 c-jld ec_conser
 cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
-      REAL t2m(klon,nbsrf), q2m(klon,nbsrf)   !temperature, humidite a 2m
-      REAL u10m(klon,nbsrf), v10m(klon,nbsrf) !vents a 10m
+      REAL,SAVE,ALLOCATABLE :: t2m(:,:), q2m(:,:)   !temperature, humidite a 2m
+c$OMP THREADPRIVATE(t2m,q2m)
+      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
@@ -1228,6 +1392,7 @@
 cjq   Aerosol effects (Johannes Quaas, 27/11/2003)
       REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
-      REAL sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
-      SAVE sulfate_pi
+      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
@@ -1239,11 +1404,15 @@
 
       ! Aerosol optical properties
-      REAL tau_ae(klon,klev,2), piz_ae(klon,klev,2)
-      REAL cg_ae(klon,klev,2)
-
-      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
+      REAL,SAVE,ALLOCATABLE :: tau_ae(:,:,:), piz_ae(:,:,:)
+c$OMP THREADPRIVATE(tau_ae,piz_ae)
+      REAL,SAVE,ALLOCATABLE :: cg_ae(:,:,:)
+c$OMP THREADPRIVATE(cg_ae)
+
+      REAL,SAVE,ALLOCATABLE :: topswad(:), solswad(:) ! Aerosol direct effect.
+c$OMP THREADPRIVATE(topswad,solswad)
       ! ok_ade=T -ADE=topswad-topsw
 
-      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
+      REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) ! Aerosol indirect effect.
+c$OMP THREADPRIVATE(topswai(,solswai)
       ! ok_aie=T ->
       !        ok_ade=T -AIE=topswai-topswad
@@ -1257,45 +1426,12 @@
 cym 
       SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
-cym
-c Anne
-         SAVE u10m
-         SAVE v10m
-         SAVE t2m
-         SAVE q2m
-         SAVE ffonte
-         SAVE fqcalving
-         SAVE piz_ae
-         SAVE tau_ae
-         SAVE cg_ae
-         SAVE rain_con
-         SAVE snow_con
-         SAVE topswai
-         SAVE topswad
-         SAVE solswai
-         SAVE solswad
-         SAVE d_u_con
-         SAVE d_v_con
-         SAVE rnebcon0
-         SAVE clwcon0
-         SAVE paire_ter
-c        SAVE nhistoW
-c        SAVE histoW
-c SAVE anne 20/09/2005
-         SAVE pblh
-         SAVE plcl
-         SAVE capCL
-         SAVE oliqCL
-         SAVE cteiCL
-         SAVE pblt
-         SAVE therm
-         SAVE trmb1
-         SAVE trmb2
-         SAVE trmb3
-
-c fin Anne
-cjq-end
+c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
+
 c
 c Declaration des constantes et des fonctions thermodynamiques
 c
+      REAL Field_tmp(klon2,klevp1)
+      LOGICAL,SAVE :: first=.true.
+c$OMP THREADPRIVATE(first)
 #include "YOMCST.h"
 #include "YOETHF.h"
@@ -1309,5 +1445,158 @@
 c
 c======================================================================
-      modname = 'physiq'
+
+cym => necessaire pour iflag_con != 2    
+      pmfd(:,:) = 0.
+      pen_u(:,:) = 0.
+      pen_d(:,:) = 0.
+      pde_d(:,:) = 0.
+      pde_u(:,:) = 0.
+      aam=0.
+      torsfc=0.
+cym => pour le couple ocean => revoir dans clmain/intersurf
+      fluxg(:)=0.
+      fluxo(:)=0.      
+
+      if (first) then 
+      
+      allocate( t_ancien(klon,klev), q_ancien(klon,klev))
+      allocate( q2(klon,klev+1,nbsrf))
+      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( ftsoil(klon,nsoilmx,nbsrf))
+      allocate( fevap(klon,nbsrf))
+      allocate( fluxlat(klon,nbsrf))
+      allocate( deltat(klon))
+      allocate( fqsurf(klon,nbsrf))
+      allocate( qsol(klon))
+      allocate( fsnow(klon,nbsrf))
+      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( agesno(klon,nbsrf))
+      allocate( alb_neig(klon))
+      allocate( run_off_lic_0(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( dlw(klon)   )
+      allocate( fder(klon) )
+      allocate( frugs(klon,nbsrf) )
+      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( therm(klon, nbsrf))
+      allocate( rain_con(klon))
+      allocate( pblt(klon, nbsrf))
+      allocate( t2m(klon,nbsrf), q2m(klon,nbsrf) )  
+      allocate( u10m(klon,nbsrf), v10m(klon,nbsrf)) 
+      allocate( topswad(klon), solswad(klon))
+      allocate( topswai(klon), solswai(klon) )
+      allocate( ffonte(klon,nbsrf))
+      allocate( fqcalving(klon,nbsrf))
+      allocate( fqfonte(klon,nbsrf)) 
+      allocate( pblh(klon, nbsrf))
+      allocate( plcl(klon, nbsrf))  
+      allocate( capCL(klon, nbsrf))   
+      allocate( oliqCL(klon, nbsrf))       
+      allocate( cteiCL(klon, nbsrf))       
+      allocate( trmb1(klon, nbsrf))   
+      allocate( trmb2(klon, nbsrf))     
+      allocate( trmb3(klon, nbsrf)) 
+      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.
+	run_off_lic_0(:)=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
@@ -1341,4 +1630,5 @@
          ffonte(:,:)=0.
          fqcalving(:,:)=0.
+	 fqfonte(:,:)=0.
          piz_ae(:,:,:)=0.
          tau_ae(:,:,:)=0.
@@ -1414,4 +1704,13 @@
      .       run_off_lic_0)
 
+       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 : pb sous surface au point ', i, 
+     $           pctsrf(i, 1 : nbsrf)
+         ENDIF 
+      ENDDO
+ 
 c   ATTENTION : il faudra a terme relire q2 dans l'etat initial
          q2(:,:,:)=1.e-8
@@ -1644,4 +1943,6 @@
 c
 #ifdef INCA
+           call VTe(VTphysiq)
+	   call VTb(VTinca)
            iii = MOD(NINT(xjour),360)
            calday = FLOAT(iii) + gmtime
@@ -1666,4 +1967,6 @@
            WRITE(lunout,*) 'OK.'
 #endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
 #endif
 c
@@ -1903,6 +2206,6 @@
        DO nsrf = 1, nbsrf
        DO i = 1, klon
-c$$$        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
-c$$$        fsollw(i,nsrf) = sollw(i)
+c@$$        fsollw(i,nsrf) = sollwdown(i) - RSIGMA*ftsol(i,nsrf)**4
+c@$$        fsollw(i,nsrf) = sollw(i)
          fsollw(i,nsrf) = sollw(i)
      $      + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf))
@@ -1910,7 +2213,24 @@
        ENDDO
        ENDDO
-
+      
+cYM   !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+cYM         Attention verrue
+cYM    ---> A supprimer plus tard       
+cYM         pour etre integre dans
+cYM         ORCHIDEE       
+      DO i = 1, klon
+        sollwdown(i)=sollw(i)+RSIGMA*ztsol(i)**4
+      ENDDO
+cYM  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
+      
       fder = dlw
 
+      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 (check) THEN
        amn=MIN(tslab(1),1000.)
@@ -1944,5 +2264,5 @@
      s            pblh,capCL,oliqCL,cteiCL,pblT,
      s            therm,trmb1,trmb2,trmb3,plcl,
-     s            fqcalving, ffonte, run_off_lic_0,
+     s            fqcalving, fqfonte,ffonte, run_off_lic_0,
 cIM "slab" ocean
      s            fluxo, fluxg, tslab, seaice) 
@@ -1986,4 +2306,13 @@
       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 
@@ -2013,4 +2342,5 @@
          zxffonte(i) = 0.0
          zxfqcalving(i) = 0.0
+	 zxfqfonte(i) = 0.0
 cIM cf. AM 081204 BEG
 c
@@ -2054,4 +2384,6 @@
             zxfqcalving(i) = zxfqcalving(i) + 
      .                      fqcalving(i,nsrf)*pctsrf(i,nsrf)
+            zxfqfonte(i) = zxfqfonte(i) + 
+     .                      fqfonte(i,nsrf)*pctsrf(i,nsrf)
 cIM cf. AM 081204 BEG
             s_pblh(i) = s_pblh(i) + pblh(i,nsrf)*pctsrf(i,nsrf)
@@ -2092,4 +2424,5 @@
            ffonte(i,nsrf) = zxffonte(i)
            fqcalving(i,nsrf) = zxfqcalving(i)
+           fqfonte(i,nsrf) = zxfqfonte(i)
            pblh(i,nsrf)=s_pblh(i)
            plcl(i,nsrf)=s_lcl(i)
@@ -2206,4 +2539,10 @@
           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
@@ -2267,4 +2606,12 @@
         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 
@@ -2467,4 +2814,12 @@
      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-------------------------------------------------------------------
@@ -2515,5 +2870,5 @@
       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 précédent diminué d'un facteur
+c  convection et du calcul du pas de temps pr��ent diminu�d'un facteur
 c  facttemps
 c      facttemps=pdtphys/1.e4
@@ -2665,4 +3020,6 @@
 
 #ifdef INCA
+      call VTe(VTphysiq)
+      call VTb(VTinca)
            calday = FLOAT(julien) + gmtime
 
@@ -2722,4 +3079,6 @@
            WRITE(lunout,*)'OK.'
 #endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
 #endif
 c     
@@ -2760,4 +3119,12 @@
      .               + 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
+      
       CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
      e            (dist, rmu0, fract, 
@@ -2788,4 +3155,12 @@
       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 
@@ -2875,4 +3250,11 @@
       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
@@ -2910,4 +3292,11 @@
 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.
@@ -2925,11 +3314,14 @@
 cIM calcul composantes axiales du moment angulaire et couple des montagnes
 c
-      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)
+      IF (monocpu) 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
@@ -3096,4 +3488,6 @@
 c
 #ifdef INCA
+      call VTe(VTphysiq)
+      call VTb(VTinca)
 #ifdef INCAINFO
            WRITE(lunout,*)'Appel CHEMHOOK_END ...'
@@ -3106,7 +3500,5 @@
      $                        nbtr,
      $                        paprs,
-#ifdef INCA_CH4
      $                        q_seri,
-#endif
      $                        annee_ref,
      $                        day_ini,
@@ -3115,6 +3507,5 @@
      $                        pphi,
      $                        pphis,
-     $                        zx_rh,
-     $                        qx(1,1,1))
+     $                        zx_rh)
 #else
      $                        xjour)
@@ -3123,4 +3514,6 @@
            WRITE(lunout,*)'OK.'
 #endif
+      call VTe(VTinca)
+      call VTb(VTphysiq)
 #endif
 
@@ -3129,4 +3522,11 @@
 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
