!
! $Header$
!
c
      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
cIM Amip2
     .            , dudyn
     .            , PVteta)

      USE ioipsl
      USE histcom
#ifdef INCA
      USE chemshut
      USE species_names
#ifdef INCA_CH4
!      USE obs_pos
#endif
#endif
      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 histins
#define histhf
#define histday
#define histmth
#define histISCCP
#define histmthNMC
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
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)
#include "dimphy.h"
#include "regdim.h"
#include "indicesol.h"
#include "dimsoil.h"
#include "clesphys.h"
#include "control.h"
#include "logic.h"
#include "temps.h"
#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======================================================================
c Parametres lies au coupleur OASIS:
#include "oasis.h"
      INTEGER,SAVE :: npas, nexca
      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      parameter (ocean = 'force ')
c     parameter (ocean = 'couple')
      logical ok_ocean
      SAVE ok_ocean
c
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
c
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     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      PARAMETER (ok_journe=.true.)
c
      LOGICAL ok_mensuel ! sortir le fichier mensuel
      save ok_mensuel
c      PARAMETER (ok_mensuel=.true.)
c
      LOGICAL ok_instan ! sortir le fichier instantane
      save ok_instan
c      PARAMETER (ok_instan=.true.)
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)
      real q2(klon,klev+1,nbsrf)
      save q2
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 zsurf(nbsrf)
cIM 
      INTEGER kinv
      real pir
cMI
      REAL u(klon,klev)
      REAL v(klon,klev)
      REAL t(klon,klev)
      REAL qx(klon,klev,nqmax)

      REAL t_ancien(klon,klev), q_ancien(klon,klev)
      SAVE t_ancien, q_ancien
      LOGICAL ancien_ok
      SAVE 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'/
      REAL rtetaSTD(nbteta)
      DATA rtetaSTD/350., 380., 405./
c
      REAL PVteta(klon,nbteta)
      REAL zx_tmp_3dte(iim,jjmp1,nbteta)
c
cMI Amip2 PV a theta constante

      INTEGER klevp1, klevm1
      PARAMETER(klevp1=klev+1,klevm1=klev-1)
#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)
      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
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./
      CHARACTER*4 clevSTD(nlevSTD)
      DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
     .'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
     .'70  ','50  ','30  ','20  ','10  '/
c
      CHARACTER*3 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 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
c
      logical oknondef(klon,nlevSTD,nout)
      real tnondef(klon,nlevSTD,nout) 
      save 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 uvsumSTD(klon,nlevSTD,nout)
      real vqsumSTD(klon,nlevSTD,nout)
      real vTsumSTD(klon,nlevSTD,nout)
      real wqsumSTD(klon,nlevSTD,nout)
c
      real vphiSTD(klon,nlevSTD)
      real wTSTD(klon,nlevSTD)
      real u2STD(klon,nlevSTD)
      real v2STD(klon,nlevSTD)
      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
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 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
      INTEGER npoints
      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 
c
      REAL nbsunlit(nregISCtot,klon)  !nbsunlit : moyenne de sunlit
      INTEGER ncol, seed(klon)

cIM ajout seed_re (REAL)et seed_old (integer)
      REAL seed_re(klon), aa
      INTEGER seed_old(klon)
      SAVE seed_old



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)
      REAL tautab(0:255)
      INTEGER invtau(-20:45000)
      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 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)
      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
      INTEGER iw, iwmax
      REAL wmin, pas_w
c     PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
      PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
      REAL o500(klon)
c
 
c sorties ISCCP

      integer nid_isccp
      save nid_isccp        
cIM 090704 BEG
      INTEGER nbapp_isccp,isccppas

      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./

c cldtopres pression au sommet des nuages
      REAL cldtopres(lmaxm1)
      DATA cldtopres/50., 180., 310., 440., 560., 680., 800./

      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'/
c
c cnameisccp
      CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
      DATA cnameisccp/'pc< 50hPa, tau< 0.3', 
     .                '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< 50hPa, tau= 0.3-1.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< 50hPa, tau= 1.3-3.6',
     .                '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< 50hPa, tau= 3.6-9.4',
     .                '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< 50hPa, tau= 9.4-23',
     .                '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< 50hPa, tau= 23-60',
     .                '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< 50hPa, tau> 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.'/
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  QUESTION : noms de variables ?

#ifdef histhf
      data ok_hf/.true./
#else
      data ok_hf/.false./
#endif

      INTEGER        longcles
      PARAMETER    ( longcles = 20 )
      REAL clesphy0( longcles      )
c
c Variables quasi-arguments
c
      REAL xjour
      SAVE xjour
c
c
c Variables propres a la physique
c
      REAL dtime
      SAVE dtime                  ! pas temporel de la physique
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
cc      INTEGER iflag_con
cc      SAVE iflag_con              ! indicateur de la convection
c
      INTEGER itap
      SAVE itap                   ! compteur pour la physique
c
      REAL co2_ppm_etat0
c
      REAL solaire_etat0
c
      real slp(klon) ! sea level pressure

      REAL ftsol(klon,nbsrf)
      SAVE ftsol                  ! temperature du sol
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

c
c
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
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
      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
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 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
      REAL wdn(klon), tdn(klon), qdn(klon)

      REAL wd(klon) ! sb
      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 ffonte(klon,nbsrf)    !Flux thermique utilise pour fondre la neige
      REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface 
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 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 rain_fall(klon) ! pluie
      REAL snow_fall(klon) ! neige
      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 evap(klon), devap(klon) ! evaporation et sa derivee
      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
      REAL dlw(klon)    ! derivee infra rouge
cym
      SAVE dlw
cym
      REAL bils(klon) ! bilan de chaleur au sol
      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 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) ! longueur de rugosite
      save frugs
      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
      REAL pctsrf(klon,nbsrf)
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
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 clmain    ! couche limite 
      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 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 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 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
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
cIM END
      REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
      REAL albpla(klon)
      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
c
      INTEGER itaprad
      SAVE 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)
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, iiq
      REAL t_coup
      PARAMETER (t_coup=234.0)
c
      REAL zphi(klon,klev)
      REAL zx_relief(iim,jjmp1)
      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
c Grdeurs 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 cape(klon)           ! CAPE
      SAVE cape
      CHARACTER*40 capemaxcels  !max(CAPE)

      REAL pbase(klon)          ! cloud base pressure
      SAVE pbase
      REAL bbase(klon)          ! cloud base buoyancy
      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 d_u_con(klon,klev),d_v_con(klon,klev)
      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 zqasc(klon,klev)
      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 ibas_con(klon), itop_con(klon)
cym
      SAVE ibas_con,itop_con
cym
      REAL rain_con(klon), rain_lsc(klon)
      REAL snow_con(klon), 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 ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
      real ratqsbas,ratqshaut
      save ratqsbas,ratqshaut, ratqs
      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
      save fact_cldcon,facttemps
      real facteur

      integer iflag_cldcon
      save 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      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)

      INTEGER        length
      PARAMETER    ( length = 100 )
      REAL tabcntr0( length       )
c
      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
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)
      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
      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 
      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
      SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
c
cIM 280405 BEG
      INTEGER nid_bilKPins, nid_bilKPave
      SAVE 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
c     REAL zstok
      REAL zsto, zout, zsto1, zsto2
c     REAL zstoave, zstoin 
      REAL zstophy, zstorad, zstohf, zstoday, zstomth
      real zjulian
      save 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
      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
      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
      REAL ZRCPD
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 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 sulfate_pi(klon, klev) ! SO4 aerosol concentration [ug/m3] (pre-industrial value)
      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 tau_ae(klon,klev,2), piz_ae(klon,klev,2)
      REAL cg_ae(klon,klev,2)

      REAL topswad(klon), solswad(klon) ! Aerosol direct effect.
      ! ok_ade=T -ADE=topswad-topsw

      REAL topswai(klon), solswai(klon) ! Aerosol indirect effect.
      ! 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
cym
cjq-end
c
c Declaration des constantes et des fonctions thermodynamiques
c
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
c======================================================================
      modname = 'physiq'
      IF (if_ebil.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.
         t2m(:,:)=0.
         q2m(:,:)=0.
         ffonte(:,:)=0.
         fqcalving(:,:)=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.
crv
c anne
         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 fin anne

cym
         wfbils(:,:)=0
cym      
         IF (if_ebil.ge.1) d_h_vcol_phy=0.
c
c appel a la lecture du run.def physique
c
         call conf_phys(pdtphys, ocean, ok_veget, ok_journe, ok_mensuel,
     .                  ok_instan, fact_cldcon, facttemps,ok_newmicro,
     .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
     .                  ok_ade, ok_aie, 
     .                  bl95_b0, bl95_b1,
     .                  iflag_thermals,nsplit_thermals)

c
c
c Initialiser les compteurs:
c

         frugs = 0.
         itap    = 0
         itaprad = 0
         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
     .       rlat,rlon,pctsrf, ftsol,ftsoil,
cIM "slab" ocean
     .       ocean, tslab,seaice,
     .       fqsurf,qsol,fsnow,
     .       falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
     .       dlw,radsol,frugs,agesno,clesphy0,
     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon, 
     .       run_off_lic_0)

c   ATTENTION : il faudra a terme relire q2 dans l'etat initial
         q2(:,:,:)=1.e-8
c
         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
         IF(ocean.NE.'force ') THEN
          ok_ocean=.TRUE.
         ENDIF
c
         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,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  "

          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
         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
         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
         WRITE(lunout,*)'La frequence de lecture surface est de ', 
     .                   lmt_pas
c
c Initialiser le couplage si necessaire
c
      npas = 0
      nexca = 0
      if (ocean == 'couple') then
        npas = itaufin/ iphysiq
        nexca = 86400 / dtime
        write(lunout,*)' ##### Ocean couple #####'
        write(lunout,*)' Valeurs des pas de temps'
        write(lunout,*)' npas = ', npas
        write(lunout,*)' nexca = ', nexca
      endif        
c
      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"
#include "ini_histday_seri.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

#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
           iii = MOD(NINT(xjour),360)
           calday = FLOAT(iii) + gmtime
           WRITE(lunout,*) 'initial time ', xjour, calday
#ifdef INCAINFO
           WRITE(lunout,*) 'Appel CHEMINI ...'
#endif
           CALL chemini( rpi,
     $                   rg,
     $                   ra,
     $                   airephy,
     $                   rlat,
     $                   rlon,
     $                   presnivs,
     $                   calday,
     $                   klon,
     $                   nqmax,
     $                   pdtphys,
     $                   annee_ref,
     $                   day_ini)
#ifdef INCAINFO
           WRITE(lunout,*) 'OK.'
#endif
#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
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
C
      IF (if_ebil.ge.1) THEN 
        ztit='after dynamic'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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
      CALL hgardfou(t_seri,ftsol,'debutphy')
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
c
      IF (if_ebil.ge.2) THEN 
        ztit='after reevap'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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
c Appeler la diffusion verticale (programme de couche limite)
c
      DO i = 1, klon
c       if (.not. ok_veget) then
c          frugs(i,is_ter) = SQRT(frugs(i,is_ter)**2+rugoro(i)**2)
c       endif 
c         frugs(i,is_lic) = rugoro(i)
c         frugs(i,is_oce) = rugmer(i)
c         frugs(i,is_sic) = 0.001
         zxrugs(i) = 0.0
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
c         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
        frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
      ENDDO
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
            zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
      ENDDO
      ENDDO
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
c
C     Calcul de l'abedo moyen par maille
      albsol(:)=0.
      albsollw(:)=0.
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
         albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
      ENDDO
      ENDDO
C
C     Repartition sous maille des flux LW et SW
C Modif OM+PASB+JLD
C Repartition du longwave par sous-surface linearisee
Cn

       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)
         fsollw(i,nsrf) = sollw(i)
     $      + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ftsol(i,nsrf))
         fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
       ENDDO
       ENDDO

      fder = dlw


      CALL clmain(dtime,itap,date0,pctsrf,pctsrf_new,
     e            t_seri,q_seri,u_seri,v_seri,
     e            julien, rmu0, co2_ppm, 
     e            ok_veget, ocean, npas, nexca, ftsol,
     $            soil_model,cdmmax, cdhmax,
     $            ksta, ksta_ter, ok_kzmin, ftsoil, qsol, 
cIM BAD    $            paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
     $            paprs,pplay, fsnow,fqsurf,fevap,falbe,falblw,
     $            fluxlat,
     e            rain_fall, snow_fall,
     e            fsolsw, fsollw, sollwdown, fder,
     e            rlon, rlat, cuphy, cvphy, frugs,
     e            debut, lafin, agesno,rugoro ,
     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
     s            q2,
     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
     s            fluxo, fluxg, tslab, seaice) 
c
CXXX PB
CXXX Incrementation des flux
CXXX

      zxfluxt=0.
      zxfluxq=0.
      zxfluxu=0.
      zxfluxv=0.
      DO nsrf = 1, nbsrf
        DO k = 1, klev
          DO i = 1, klon
            zxfluxt(i,k) = zxfluxt(i,k) + 
     $          fluxt(i,k,nsrf) * pctsrf( i, nsrf)
            zxfluxq(i,k) = zxfluxq(i,k) + 
     $          fluxq(i,k,nsrf) * pctsrf( i, nsrf)
            zxfluxu(i,k) = zxfluxu(i,k) + 
     $          fluxu(i,k,nsrf) * pctsrf( i, nsrf)
            zxfluxv(i,k) = zxfluxv(i,k) + 
     $          fluxv(i,k,nsrf) * pctsrf( i, nsrf)
          END DO 
        END DO 
      END DO 
      DO i = 1, klon
         sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
c         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
         evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
         fder(i) = dlw(i) + dsens(i) + devap(i)
      ENDDO


      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
c
      IF (if_ebil.ge.2) THEN 
        ztit='after clmain'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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
c Incrementer la temperature du sol
c
      DO i = 1, klon
         zxtsol(i) = 0.0
         zxfluxlat(i) = 0.0
c
         zt2m(i) = 0.0
         zq2m(i) = 0.0
         zu10m(i) = 0.0
         zv10m(i) = 0.0
cIM cf JLD ??
         zxffonte(i) = 0.0
         zxfqcalving(i) = 0.0
cIM cf. AM 081204 BEG
c
         s_pblh(i) = 0.0 
         s_lcl(i) = 0.0 
         s_capCL(i) = 0.0
         s_oliqCL(i) = 0.0
         s_cteiCL(i) = 0.0
         s_pblT(i) = 0.0
         s_therm(i) = 0.0
         s_trmb1(i) = 0.0
         s_trmb2(i) = 0.0
         s_trmb3(i) = 0.0
c
         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
      DO nsrf = 1, nbsrf
        DO i = 1, klon
c        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN 
            ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
cIM cf. JLD
            wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
     $         + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
            zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
cccIM
            zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
            zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
            zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
            zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
cIM cf JLD ??
            zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
            zxfqcalving(i) = zxfqcalving(i) + 
     .                      fqcalving(i,nsrf)*pctsrf(i,nsrf)
cIM cf. AM 081204 BEG
            s_pblh(i) = s_pblh(i) + pblh(i,nsrf)*pctsrf(i,nsrf)
            s_lcl(i) = s_lcl(i) + plcl(i,nsrf)*pctsrf(i,nsrf)
            s_capCL(i) = s_capCL(i) + capCL(i,nsrf) *pctsrf(i,nsrf)
            s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf) *pctsrf(i,nsrf)
            s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf) *pctsrf(i,nsrf)
            s_pblT(i) = s_pblT(i) + pblT(i,nsrf) *pctsrf(i,nsrf)
            s_therm(i) = s_therm(i) + therm(i,nsrf) *pctsrf(i,nsrf)
            s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) *pctsrf(i,nsrf)
            s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) *pctsrf(i,nsrf)
            s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) *pctsrf(i,nsrf)
c        ENDIF 
        ENDDO
      ENDDO

c
c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
c
      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)
        ENDDO
      ENDDO
c
c
c Calculer la derive du flux infrarouge
c
cXXX      DO nsrf = 1, nbsrf
      DO i = 1, klon
cXXX        IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN 
            dlw(i) = - 4.0*RSIGMA*zxtsol(i)**3 
cXXX     .          *(ftsol(i,nsrf)-zxtsol(i))
cXXX     .          *pctsrf(i,nsrf)
cXXX        ENDIF 
cXXX      ENDDO
      ENDDO
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, zxfluxq(1,1), 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

           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
c
      IF (if_ebil.ge.2) THEN 
        ztit='after convect'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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===================================================================
c
      IF (if_ebil.ge.2) THEN 
        ztit='after dry_adjust'
        CALL diagetpq(airephy,ztit,ip_ebil,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
c
      IF (if_ebil.ge.2) THEN 
        ztit='after fisrt'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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 
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 précédent diminué 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
cIM calcul nuages par le simulateur ISCCP
c
      IF (ok_isccp) THEN
#include "calcul_simulISCCP.h"
      ENDIF !ok_isccp

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
c
      IF (if_ebil.ge.2) THEN 
        ztit="after diagcld"
        CALL diagetpq(airephy,ztit,ip_ebil,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
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
         ! 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)
cym
      ELSE
        tau_ae(:,:,:)=0.0
	piz_ae(:,:,:)=0.0
	cg_ae(:,:,:)=0.0
cym     
      ENDIF

#ifdef INCA
           calday = FLOAT(julien) + gmtime

#ifdef INCA_AER
      call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
     &   prfl,psfl,pctsrf(1,3),airephy,xjour,rlat,rlon)
#endif

#ifdef INCAINFO
           WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
#endif

           CALL chemhook_begin (calday,
     $                          pctsrf(1,1),
     $                          rlat,
     $                          rlon,
     $                          airephy,
     $                          paprs,
     $                          pplay,
     $                          ycoefh,
     $                          pphi,
     $                          t_seri,
     $                          u,
     $                          v,
     $                          wo,
     $                          q_seri,
     $                          zxtsol,
     $                          zxsnow,
     $                          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
#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,
     e            sulfate, sulfate_pi,
     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
      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) ! ="=
      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 (if_ebil.ge.2) THEN 
        ztit='after rad'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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
      DO i = 1, klon
         zxqsurf(i) = 0.0
         zxsnow(i) = 0.0
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
         zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
      ENDDO
      ENDDO
c
c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
c
cXXX      DO nsrf = 1, nbsrf
cXXX      DO i = 1, klon
cXXX         IF (pctsrf(i,nsrf).LT.epsfra) THEN
cXXX            fqsurf(i,nsrf) = zxqsurf(i)
cXXX            fsnow(i,nsrf) = zxsnow(i)
cXXX         ENDIF
cXXX      ENDDO
cXXX      ENDDO
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 (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

      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
      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)
cIM cf. FLott END
c
      IF (if_ebil.ge.2) THEN 
        ztit='after orography'
        CALL diagetpq(airephy,ztit,ip_ebil,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,
#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 diag. bilKP
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
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
      IF (if_ebil.ge.1) THEN 
        ztit='after physic'
        CALL diagetpq(airephy,ztit,ip_ebil,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
     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
#ifdef INCAINFO
           WRITE(lunout,*)'Appel CHEMHOOK_END ...'
#endif
           CALL chemhook_end (calday,
     $                        dtime,
     $                        pplay,
     $                        t_seri,
     $                        tr_seri,
     $                        nbtr,
     $                        paprs,
#ifdef INCA_CH4
     $                        q_seri,
#endif
     $                        annee_ref,
     $                        day_ini,
#ifdef INCA_AER
     $                        xjour,
     $                        pphi,
     $                        pphis,
     $                        zx_rh,
     $                        qx(1,1,1))
#else
     $                        xjour)
#endif
#ifdef INCAINFO
           WRITE(lunout,*)'OK.'
#endif
#endif

c=============================================================
c
c Convertir les incrementations en tendances
c
      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
c#include "write_bilKP_ins.h"
c#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

#ifdef histhf
#include "write_histhf.h"
#endif

#ifdef histday
#include "write_histday.h"
#include "write_histday_seri.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

#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
ccc         IF (ok_oasis) CALL quitcpl
         CALL phyredem ("restartphy.nc",dtime,radpas,
     .      rlat, rlon, pctsrf, ftsol, ftsoil,
cIM "slab" ocean
     .      tslab, seaice, 
     .      fqsurf, qsol,
     .      fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
     .      solsw, sollwdown,dlw,
     .      radsol,frugs,agesno,
     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,run_off_lic_0)
      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
