Index: /LMDZ4/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 686)
+++ /LMDZ4/trunk/libf/phylmd/physiq.F	(revision 687)
@@ -42,7 +42,6 @@
 #define histmth
 #define histins
-#define histISCCP
-#define histREGDYN
-#define histmthNMC
+c #define histISCCP
+c #define histmthNMC
 c======================================================================
 c    modif   ( P. Le Van ,  12/10/98 )
@@ -70,5 +69,8 @@
 c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
 c omega---input-R-vitesse verticale en Pa/s
-c
+cIM comgeomphy.h BEG
+c cuphy----input-R-resolution des mailles en x (m)
+c cvphy----input-R-resolution des mailles en y (m)
+cIM comgeomphy.h END
 c d_u-----output-R-tendance physique de "u" (m/s/s)
 c d_v-----output-R-tendance physique de "v" (m/s/s)
@@ -125,10 +127,11 @@
 cIM "slab" ocean
       REAL tslab(klon)    !Temperature du slab-ocean
-      SAVE tslab
+cIM      SAVE tslab
       REAL seaice(klon)   !glace de mer (kg/m2) 
-      SAVE seaice
+cIM      SAVE seaice
       REAL fluxo(klon)    !flux turbulents ocean-glace de mer 
       REAL fluxg(klon)    !flux turbulents ocean-atmosphere
-c
+cIM
+      REAL amn, amx
 c======================================================================
 c Clef controlant l'activation du cycle diurne:
@@ -250,5 +253,4 @@
 #include "raddim.h"
 c
-cIM 080304   REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
       REAL swdn0(klon,klevp1), swdn(klon,klevp1)
       REAL swup0(klon,klevp1), swup(klon,klevp1)
@@ -380,13 +382,37 @@
       INTEGER imin_debut, nbpti
       INTEGER jmin_debut, nbptj 
-c
-      REAL nbsunlit(nregISCtot,klon)  !nbsunlit : moyenne de sunlit
-      INTEGER ncol, seed(klon)
-
+cIM parametres ISCCP BEG
+      INTEGER nbapp_isccp,isccppas
+      INTEGER n, napisccp
+c     PARAMETER(napisccp=3)
+      PARAMETER(napisccp=1)
+      INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
+      DATA ifreq_isccp/3/
+      SAVE ifreq_isccp
+      CHARACTER*5 typinout(napisccp)
+      DATA typinout/'i3od'/
+cIM verif boxptop BEG
+      CHARACTER*1 verticaxe(napisccp)
+      DATA verticaxe/'1'/ 
+cIM verif boxptop END
+      INTEGER nvlev(napisccp)
+c     INTEGER nvlev
+      REAL t1, aa
+      REAL seed_re(klon,napisccp)
+      INTEGER seed_old(klon,napisccp)
+      SAVE seed_old
+      INTEGER iphy(iim,jjmp1)
+cIM parametres ISCCP END
+c
 c ncol = nb. de sous-colonnes pour chaque maille du GCM 
-c     PARAMETER(ncol=100)
-c     PARAMETER(ncol=625)
-c     PARAMETER(ncol=10)
-      PARAMETER(ncol=25)
+c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM 
+      INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
+      REAL nbsunlit(nregISCtot,klon,napisccp)  !nbsunlit : moyenne de sunlit
+      PARAMETER(ncolmx=1500)
+c
+cIM verif boxptop BEG
+      REAL vertlev(ncolmx,napisccp)
+cIM verif boxptop END
+c
       REAL tautab(0:255)
       INTEGER invtau(-20:45000)
@@ -413,15 +439,7 @@
       REAL dem_sH2B(klon,klev)
       REAL dem_cH2B(klon,klev)
-
-c output from ISCCP simulator
-      REAL fq_isccp(klon,7,7)
-      REAL totalcldarea(klon) 
-      REAL meanptop(klon)
-      REAL meantaucld(klon)
-      REAL boxtau(klon,ncol)
-      REAL boxptop(klon,ncol) 
-c
-      INTEGER l, kmax, lmax
-      PARAMETER(kmax=8, lmax=8)
+c
+      INTEGER kmax, lmax, lmax3
+      PARAMETER(kmax=8, lmax=8, lmax3=3)
       INTEGER kmaxm1, lmaxm1
       PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
@@ -430,8 +448,23 @@
      .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)
-      PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
+cIM 051005     PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
+      PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
       REAL o500(klon)
 c
@@ -440,11 +473,14 @@
       INTEGER nreg, nbregdyn
       PARAMETER(nbregdyn=5)
-      REAL histoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
-      REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbregdyn)
-      REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbregdyn) 
-      SAVE nhistoWt
+cIM 051005 BEG
+c     REAL histoW(iwmax,nbregdyn,napisccp)
+c     REAL nhistoW(iwmax,nbregdyn,napisccp)
+c     REAL histoWn(iwmax,nbregdyn)
+c     REAL nhistoWn(iwmax,nbregdyn)
+cIM 090905 END
 
       INTEGER linv
       INTEGER pct_ocean(klon,nbregdyn)
+      SAVE pct_ocean
  
 c sorties ISCCP
@@ -455,6 +491,4 @@
 c     save ok_isccp, ecrit_isccp, nid_isccp        
       save nid_isccp        
-cIM 090704 BEG
-      INTEGER nbapp_isccp,isccppas
 
 #undef histISCCP
@@ -467,28 +501,19 @@
 cIM 190504 #endif
 
-c sorties statistiques regime dynamique
-c     logical ok_regdyn
-c     real ecrit_regdyn
-      integer nid_regdyn
-c     save ok_regdyn, ecrit_regdyn, nid_regdyn
-      save nid_regdyn
-
-#undef histREGDYN
-#define histREGDYN
-cIM 190504 #ifdef histREGDYN
-c     data ok_regdyn,ecrit_regdyn/.true.,0.125/
-c     data ok_regdyn,ecrit_regdyn/.true.,1./
-cIM 190504    data ok_regdyn/.true./
-cIM 190504 #else
-cIM 190504   data ok_regdyn/.false./
-cIM 190504 #endif 
-
       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./
-      DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
+cIM bad 151205     DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
+      DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
 
 c cldtopres pression au sommet des nuages
-      REAL cldtopres(lmaxm1)
-      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
+      REAL cldtopres(lmaxm1), cldtopres3(lmax3)
+cIM 151205 erreur     DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
+      DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
+      DATA cldtopres3/440., 680., 1000./
+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) 
+cIM 051005 END
 
       INTEGER komega, nhoriRD 
@@ -503,6 +528,6 @@
 c cnameisccp
       CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
-      DATA cnameisccp/'pc< 50hPa, tau< 0.3', 
-     .                'pc= 50-180hPa, tau< 0.3',
+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',
@@ -510,5 +535,5 @@
      .                'pc= 560-680hPa, tau< 0.3',
      .                'pc= 680-800hPa, tau< 0.3',
-     .                'pc< 50hPa, tau= 0.3-1.3',
+     .                'pc= 800-1000hPa, tau< 0.3',
      .                'pc= 50-180hPa, tau= 0.3-1.3',
      .                'pc= 180-310hPa, tau= 0.3-1.3',
@@ -517,5 +542,5 @@
      .                'pc= 560-680hPa, tau= 0.3-1.3',
      .                'pc= 680-800hPa, tau= 0.3-1.3',
-     .                'pc< 50hPa, tau= 1.3-3.6',
+     .                'pc= 800-1000hPa, tau= 0.3-1.3',
      .                'pc= 50-180hPa, tau= 1.3-3.6',
      .                'pc= 180-310hPa, tau= 1.3-3.6',
@@ -524,5 +549,5 @@
      .                'pc= 560-680hPa, tau= 1.3-3.6',
      .                'pc= 680-800hPa, tau= 1.3-3.6',
-     .                'pc< 50hPa, tau= 3.6-9.4',
+     .                'pc= 800-1000hPa, tau= 1.3-3.6',
      .                'pc= 50-180hPa, tau= 3.6-9.4',
      .                'pc= 180-310hPa, tau= 3.6-9.4',
@@ -531,5 +556,5 @@
      .                'pc= 560-680hPa, tau= 3.6-9.4',
      .                'pc= 680-800hPa, tau= 3.6-9.4',
-     .                'pc< 50hPa, tau= 9.4-23',
+     .                'pc= 800-1000hPa, tau= 3.6-9.4',
      .                'pc= 50-180hPa, tau= 9.4-23',
      .                'pc= 180-310hPa, tau= 9.4-23',
@@ -538,5 +563,5 @@
      .                'pc= 560-680hPa, tau= 9.4-23',
      .                'pc= 680-800hPa, tau= 9.4-23',
-     .                'pc< 50hPa, tau= 23-60',
+     .                'pc= 800-1000hPa, tau= 9.4-23',
      .                'pc= 50-180hPa, tau= 23-60',
      .                'pc= 180-310hPa, tau= 23-60',
@@ -545,5 +570,5 @@
      .                'pc= 560-680hPa, tau= 23-60',
      .                'pc= 680-800hPa, tau= 23-60',
-     .                'pc< 50hPa, tau> 60.',
+     .                'pc= 800-1000hPa, tau= 23-60',
      .                'pc= 50-180hPa, tau> 60.',
      .                'pc= 180-310hPa, tau> 60.',
@@ -551,5 +576,6 @@
      .                'pc= 440-560hPa, tau> 60.',
      .                'pc= 560-680hPa, tau> 60.',
-     .                'pc= 680-800hPa, tau> 60.'/
+     .                'pc= 680-800hPa, tau> 60.',
+     .                'pc= 800-1000hPa, tau> 60.'/
 c
 c     REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
@@ -562,14 +588,12 @@
 c
       logical ok_hf
-cIM200505     integer ecrit_hf
-cIM200505    integer ecrit_hf2mth
-cIM200505    save ecrit_hf2mth
 c
       integer nid_hf, nid_hf3d
-cIM200505     save ok_hf, ecrit_hf, nid_hf, nid_hf3d
       save ok_hf, nid_hf, nid_hf3d
 
 c  QUESTION : noms de variables ?
 
+#undef histhf
+#define histhf
 #ifdef histhf
 cIM 130904   data ok_hf,ecrit_hf/.true.,0.25/
@@ -623,4 +647,7 @@
       REAL ftsol(klon,nbsrf)
       SAVE ftsol                  ! temperature du sol
+cIM
+      REAL newsst(klon) !temperature de l'ocean
+      SAVE newsst
 c
       REAL ftsoil(klon,nsoilmx,nbsrf)
@@ -778,4 +805,6 @@
 cym
       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
@@ -810,4 +839,7 @@
       REAL wo(klon,klev)
       SAVE wo                     ! ozone
+cIM sorties
+      REAL un_jour
+      PARAMETER(un_jour=86400.)
 c======================================================================
 c
@@ -924,5 +956,5 @@
       REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
       real zqsat(klon,klev)
-      INTEGER i, k, iq, ig, j, nsrf, ll, iiq
+      INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq
       REAL t_coup
       PARAMETER (t_coup=234.0)
@@ -1050,6 +1082,4 @@
 c
 c======================================================================
-cIM200505     INTEGER ecrit_mth
-cIM200505     SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
 c
 cIM cf. AM 081204 BEG
@@ -1068,13 +1098,4 @@
 
 c
-cIM200505     INTEGER ecrit_day
-cIM200505     SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
-c
-cIM200505     INTEGER ecrit_ins
-cIM200505     SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
-c
-cIM200505     INTEGER ecrit_reg
-cIM200505     SAVE ecrit_reg   ! frequence d'ecriture
-c
       integer itau_w   ! pas de temps ecriture = itap + itau_phy
 c
@@ -1090,4 +1111,11 @@
 
       REAL zx_rh(klon,klev)
+cIM RH a 2m (la surface)
+      REAL rh2m(klon), qsat2m(klon)
+      REAL zx_rh2m(klon,nbsrf), zx_qsat2m(klon,nbsrf)
+      REAL zx_qs1(klon,nbsrf), zx_t1(klon,nbsrf), zdelta1(klon,nbsrf)
+      REAL zcor1(klon,nbsrf)
+      REAL tpot(klon), tpote(klon)
+      REAL Lheat
 
       INTEGER        length
@@ -1096,4 +1124,6 @@
 c
       INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
+cIM
+      INTEGER ndex2d1(iwmax)
 c
 cIM AMIP2 BEG
@@ -1127,5 +1157,7 @@
 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
 cIM 280405 BEG
@@ -1140,9 +1172,11 @@
 cIM 280405 END
 c
-      INTEGER nhori, nvert, nvert1
-c     REAL zstok
-      REAL zsto, zout, zsto1, zsto2
-c     REAL zstoave, zstoin 
-      REAL zstophy, zstorad, zstohf, zstoday, zstomth
+      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
+
       real zjulian
       save zjulian
@@ -1176,9 +1210,5 @@
       REAL      zero_v(klon)
       CHARACTER*15 ztit
-      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
-      SAVE      ip_ebil
-      DATA      ip_ebil/0/
-      INTEGER   if_ebil ! level for energy conserv. dignostics
-      SAVE      if_ebil
+
 c+jld ec_conser
       REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
@@ -1245,6 +1275,6 @@
          SAVE clwcon0
          SAVE paire_ter
-         SAVE nhistoW
-         SAVE histoW
+c        SAVE nhistoW
+c        SAVE histoW
 c SAVE anne 20/09/2005
          SAVE pblh
@@ -1267,7 +1297,15 @@
 #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======================================================================
       modname = 'physiq'
-      IF (if_ebil.ge.1) THEN
+cIM
+      IF (ip_ebil_phy.ge.1) THEN
         DO i=1,klon
           zero_v(i)=0.
@@ -1319,6 +1357,6 @@
          clwcon(:,:) = 0.0
          paire_ter(:) = 0.0
-         nhistoW(:,:,:,:) = 0.0
-         histoW(:,:,:,:) = 0.0
+c        nhistoW(:,:,:,:) = 0.0
+c        histoW(:,:,:,:) = 0.0
 ! fin anne
 ! Anne 12/09/2005
@@ -1339,5 +1377,6 @@
          wfbils(:,:)=0
 cym      
-         IF (if_ebil.ge.1) d_h_vcol_phy=0.
+cIM      
+         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
 c
 c appel a la lecture du run.def physique
@@ -1345,5 +1384,6 @@
          call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
      .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
-     .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
+cIM  .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
+     .                  iflag_cldcon,ratqsbas,ratqshaut,
      .                  ok_ade, ok_aie, 
      .                  bl95_b0, bl95_b1,
@@ -1363,5 +1403,6 @@
      .       ocean, tslab,seaice,
      .       fqsurf,qsol,fsnow,
-     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
+cIM 220306  .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
+     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollw,
      .       dlw,radsol,frugs,agesno,clesphy0,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
@@ -1425,4 +1466,8 @@
 
          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
@@ -1500,14 +1545,18 @@
 cIM200505    .                   ecrit_reg
 cIM200505        ENDIF
-c
-cIM 230505 BEG
-         ecrit_ins = NINT(ecrit_ins/dtime)
-         ecrit_hf = NINT(ecrit_hf/dtime)
-c        ecrit_hf2mth = 4*30
-         ecrit_day = NINT(ecrit_day/dtime)
-         ecrit_mth = NINT(ecrit_mth/dtime)
-         ecrit_tra = NINT(86400.*ecrit_tra/dtime)
-         ecrit_reg = NINT(ecrit_reg/dtime)
-cIM 230505 END
+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, ecrit_tra en secondes, chaque pas de temps de la physique
+         ecrit_ins = dtime
+         ecrit_tra = 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
+         ecrit_day = ecrit_day * un_jour
+         ecrit_mth = ecrit_mth * un_jour
+         ecrit_reg = ecrit_reg * un_jour
+cIM 030306 END
 c
 c Initialiser le couplage si necessaire
@@ -1548,5 +1597,4 @@
 c#include "ini_bilKP_ins.h"
 c#include "ini_bilKP_ave.h"
-#include "ini_histday_seri.h"
 #endif
 
@@ -1563,17 +1611,13 @@
 #endif
 
+#undef histmthNMC
+#define histmthNMC
 #ifdef histmthNMC
 #include "ini_histmthNMC.h"
 #endif
 
-#ifdef histREGDYN
-#include "ini_histREGDYN.h"
-#endif
-
-c#undef histmthNMC
-c#define histmthNMC
-#ifdef histmthNMC
-#include "ini_histmthNMC.h"
-#endif
+#include "ini_histday_seri.h"
+
+#include "ini_paramLMDZ_phy.h"
 
 #endif
@@ -1684,8 +1728,8 @@
         ENDDO
       ENDDO
-C
-      IF (if_ebil.ge.1) THEN 
+cIM
+      IF (ip_ebil_phy.ge.1) THEN 
         ztit='after dynamic'
-        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
+        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)
@@ -1694,5 +1738,5 @@
 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
+        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
@@ -1730,5 +1774,31 @@
 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(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
@@ -1764,11 +1834,11 @@
       ENDDO
       ENDDO
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after reevap'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
+        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
+         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
@@ -1839,4 +1909,15 @@
       fder = dlw
 
+      IF (check) THEN
+       amn=MIN(tslab(1),1000.)
+       amx=MAX(tslab(1),-1000.)
+       DO i=2, klon
+        amn=MIN(tslab(i),amn)
+        amx=MAX(tslab(i),amx)
+       ENDDO
+c
+       PRINT*,' debut avant clqh min max tslab',amn,amx
+      ENDIF !(check) THEN
+c
       CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
      e            t_seri,q_seri,u_seri,v_seri,
@@ -1857,8 +1938,6 @@
      s            dsens, devap,
      s            ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
-cIM cf. AM 081204 BEG
      s            pblh,capCL,oliqCL,cteiCL,pblT,
      s            therm,trmb1,trmb2,trmb3,plcl,
-cIM cf. AM 081204 END
      s            fqcalving, ffonte, run_off_lic_0,
 cIM "slab" ocean
@@ -1903,11 +1982,11 @@
       ENDDO
       ENDDO
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after clmain'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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
+         call diagphy(airephy,ztit,ip_ebil_phy
      e      , zero_v, zero_v, zero_v, zero_v, sens
      e      , evap  , zero_v, zero_v, ztsol
@@ -1957,4 +2036,7 @@
             wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
      $         + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
+cIM
+            wfbilo(i,nsrf) = ( fevap(i,nsrf) - 
+     $      (rain_fall(i) + snow_fall(i)) ) * pctsrf(i,nsrf)
             zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
             zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
@@ -1983,4 +2065,14 @@
       ENDDO
 
+      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 d_ts min max ftsol',itap,amn,amx
+      ENDIF !(check) THEN
 c
 c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
@@ -1988,28 +2080,25 @@
       DO nsrf = 1, nbsrf
         DO i = 1, klon
-          IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
-cccIM
-          IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
-cIM cf JLD ??
-          IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) 
-     .    fqcalving(i,nsrf) = zxfqcalving(i)
-cIM cf. AM 081204 BEG
-          IF (pctsrf(i,nsrf) .LT. epsfra) pblh(i,nsrf)=s_pblh(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) plcl(i,nsrf)=s_lcl(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) capCL(i,nsrf)=s_capCL(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) oliqCL(i,nsrf)=s_oliqCL(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) cteiCL(i,nsrf)=s_cteiCL(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) pblT(i,nsrf)=s_pblT(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) therm(i,nsrf)=s_therm(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) trmb1(i,nsrf)=s_trmb1(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) trmb2(i,nsrf)=s_trmb2(i)
-          IF (pctsrf(i,nsrf) .LT. epsfra) trmb3(i,nsrf)=s_trmb3(i)
+          IF (pctsrf(i,nsrf) .LT. epsfra.OR.t2m(i,nsrf).EQ.0.) THEN
+           ftsol(i,nsrf) = zxtsol(i)
+           t2m(i,nsrf) = zt2m(i)
+           q2m(i,nsrf) = zq2m(i)
+           u10m(i,nsrf) = zu10m(i)
+           v10m(i,nsrf) = zv10m(i)
+           ffonte(i,nsrf) = zxffonte(i)
+           fqcalving(i,nsrf) = zxfqcalving(i)
+           pblh(i,nsrf)=s_pblh(i)
+           plcl(i,nsrf)=s_lcl(i)
+           capCL(i,nsrf)=s_capCL(i)
+           oliqCL(i,nsrf)=s_oliqCL(i) 
+           cteiCL(i,nsrf)=s_cteiCL(i)
+           pblT(i,nsrf)=s_pblT(i)
+           therm(i,nsrf)=s_therm(i)
+           trmb1(i,nsrf)=s_trmb1(i)
+           trmb2(i,nsrf)=s_trmb2(i)
+           trmb3(i,nsrf)=s_trmb3(i)
+          ENDIF
         ENDDO
       ENDDO
-c
 c
 c Calculer la derive du flux infrarouge
@@ -2174,11 +2263,11 @@
         ENDDO
       ENDDO
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after convect'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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
+         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
@@ -2264,8 +2353,8 @@
 c
 c===================================================================
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after dry_adjust'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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)
@@ -2362,11 +2451,11 @@
          WRITE(lunout,*)"Precip=", zx_t
       ENDIF
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after fisrt'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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
+        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
@@ -2471,8 +2560,8 @@
          snow_fall(i) = snow_con(i) + snow_lsc(i)
       ENDDO
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit="after diagcld"
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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)
@@ -2501,4 +2590,56 @@
       ENDDO
       ENDDO
+c
+cIM Calculer l'humidite relative a 2m (rh2m) pour diagnostique
+cIM ajout dependance type surface
+      DO i = 1, klon
+       rh2m(i)=0.
+       qsat2m(i)=0.
+      DO nsrf=1, nbsrf
+         zx_t1(i,nsrf) = t2m(i,nsrf)
+         IF (thermcep) THEN
+            zdelta1(i,nsrf) = MAX(0.,SIGN(1.,rtt-zx_t1(i,nsrf)))
+            zx_qs1(i,nsrf)  = r2es * 
+     $      FOEEW(zx_t1(i,nsrf),zdelta1(i,nsrf))/paprs(i,1)
+            zx_qs1(i,nsrf)  = MIN(0.5,zx_qs1(i,nsrf))
+            zcor1(i,nsrf)   = 1./(1.-retv*zx_qs1(i,nsrf))
+            zx_qs1(i,nsrf)  = zx_qs1(i,nsrf)*zcor1(i,nsrf)
+         ELSE
+c
+           IF (zx_t.LT.RTT) THEN
+              zx_qs = qsats(zx_t)/paprs(i,1)
+           ELSE
+              zx_qs = qsatl(zx_t)/paprs(i,1)
+           ENDIF
+         ENDIF
+       zx_rh2m(i,nsrf) = q2m(i,nsrf)/zx_qs1(i,nsrf)
+       zx_qsat2m(i,nsrf)=zx_qs1(i,nsrf)
+       rh2m(i) = rh2m(i)+zx_rh2m(i,nsrf)*pctsrf(i,nsrf)
+       qsat2m(i)=qsat2m(i)+zx_qsat2m(i,nsrf)*pctsrf(i,nsrf)
+      ENDDO !nsrf
+      ENDDO
+c
+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
+c
 cjq - introduce the aerosol direct and first indirect radiative forcings
 cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
@@ -2643,11 +2784,11 @@
       ENDDO
       ENDDO
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after rad'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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
+        call diagphy(airephy,ztit,ip_ebil_phy
      e      , topsw, toplw, solsw, sollw, zero_v
      e      , zero_v, zero_v, zero_v, ztsol
@@ -2788,8 +2929,8 @@
      C               aam, torsfc)
 cIM cf. FLott END
-c
-      IF (if_ebil.ge.2) THEN 
+cIM
+      IF (ip_ebil_phy.ge.2) THEN 
         ztit='after orography'
-        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
+        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)
@@ -2883,5 +3024,6 @@
      s                   ve, vq, ue, uq)
 c
-cIM diag. bilKP
+cIM global posePB BEG
+      IF(1.EQ.0) THEN
 c
       CALL transp_lay (paprs,zxtsol,
@@ -2889,4 +3031,6 @@
      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
@@ -2902,7 +3046,8 @@
       END DO 
 c-jld ec_conser
-      IF (if_ebil.ge.1) THEN 
+cIM
+      IF (ip_ebil_phy.ge.1) THEN 
         ztit='after physic'
-        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
+        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)
@@ -2911,5 +3056,5 @@
 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
+        call diagphy(airephy,ztit,ip_ebil_phy
      e      , topsw, toplw, solsw, sollw, sens
      e      , evap, rain_fall, snow_fall, ztsol
@@ -3001,6 +3146,6 @@
 c
 cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
-c#include "write_bilKP_ins.h"
-c#include "write_bilKP_ave.h"
+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:
@@ -3025,5 +3170,4 @@
 #ifdef histday
 #include "write_histday.h"
-#include "write_histday_seri.h"
 #endif
 
@@ -3036,16 +3180,15 @@
 #endif
 
-#ifdef histREGDYN
-#include "write_histREGDYN.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
@@ -3066,5 +3209,6 @@
      .      fqsurf, qsol,
      .      fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
-     .      solsw, sollwdown,dlw,
+cIM  .      solsw, sollwdown,dlw,
+     .      solsw, sollw,dlw,
      .      radsol,frugs,agesno,
      .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
