c $Header$
c
      SUBROUTINE physiq (nlon,
     .                   nlev,
     .                   nqmax,
     .                   debut,
     .                   lafin,
     .                   rjourvrai,
     .                   rjour_ecri,
     .                   gmtime,
     .                   pdtphys,
     .                   paprs,
     .                   pplay,
     .                   pphi,
     .                   pphis,
     .                   paire,
     .                   presnivs,
     .                   clesphy0,
     .                   u,
     .                   v,
     .                   t,
     .                   qx,
     .                   omega,
#ifdef INCA_CH4
     .                   flxmass_w,
#endif
     .                   d_u, 
     .                   d_v, 
     .                   d_t, 
     .                   d_qx, 
     .                   d_ps)
      USE ioipsl
#ifdef INCA
      USE chemshut 
      USE obs_pos
#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    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 paire---input-R-aire de chaque maille
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
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 "temps.h"
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 npas, nexca, itimestep
      logical rnpb
#ifdef INCA
      parameter(rnpb=.false.)
#else
      parameter(rnpb=.true.)
#endif
      PARAMETER (npas=1440)
      PARAMETER (nexca=48)
      PARAMETER (itimestep=1800)
      EXTERNAL fromcpl, intocpl, inicma
      REAL cpl_sst(iim,jjmp1), cpl_sic(iim,jjmp1)
      REAL cpl_alb_sst(iim,jjmp1), cpl_alb_sic(iim,jjmp1)
c======================================================================
c ok_ocean indique l'utilisation du modele oceanique "slab ocean",
c il faut bien sur s'assurer que le bilan energetique de reference
c a la surface de l'ocean est bien present dans le fichier des
c conditions aux limites, ainsi que l'indicateur du sol ne contient
c pas de glace oceanique (pas de valeur 3).
c
      LOGICAL ok_ocean
      PARAMETER (ok_ocean=.FALSE.)
      REAL cyang  ! capacite thermique de l'ocean superficiel
      PARAMETER (cyang=30.0 * 4.228e+06)
      REAL cbing  ! capacite thermique de la glace oceanique
      PARAMETER (cbing=1.0 * 4.228e+06)
      REAL cthermiq
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.)
      REAL soilcap(klon,nbsrf), soilflux(klon,nbsrf)
      SAVE soilcap, soilflux
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
      PARAMETER (ok_journe=.FALSE.)
c
      LOGICAL ok_mensuel ! sortir le fichier mensuel
      PARAMETER (ok_mensuel=.TRUE.)
c
      LOGICAL ok_instan ! sortir le fichier instantane
      PARAMETER (ok_instan=.FALSE.)
c
      LOGICAL ok_region ! sortir le fichier regional
      PARAMETER (ok_region=.FALSE.)
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
      INTEGER nvm           ! nombre de vegetations
      PARAMETER (nvm=8)
      REAL veget(klon,nvm)  ! couverture vegetale
      SAVE veget
c
c Variables argument:
c
      INTEGER nlon
      INTEGER nlev
      INTEGER nqmax
      REAL rjourvrai, rjour_ecri
      REAL gmtime
      REAL pdtphys
      LOGICAL debut, lafin
      REAL paprs(klon,klev+1)
      REAL pplay(klon,klev)
      REAL pphi(klon,klev)
      REAL pphis(klon)
      REAL paire(klon)
      REAL presnivs(klev)
      REAL znivsig(klev)

      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_u_dyn(klon,klev)
      REAL d_v_dyn(klon,klev)
      REAL d_t_dyn(klon,klev)
      REAL d_q_dyn(klon,klev)

      REAL omega(klon,klev)

      REAL flxmass_w(klon,klev)

      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)

      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
c
      REAL rlat(klon)
      SAVE rlat                   ! latitude pour chaque point
c
      REAL rlon(klon)
      SAVE rlon                   ! longitude 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
      SAVE co2_ppm                ! concentration du CO2
c
      REAL solaire
      SAVE solaire                ! constante solaire
c
      REAL ftsol(klon,nbsrf)
      SAVE ftsol                  ! temperature du sol
c
      REAL ftsoil(klon,nsoilmx,nbsrf)
      SAVE ftsoil                 ! temperature dans le sol
c
      REAL deltat(klon)
      SAVE deltat                 ! ecart avec la SST de reference
c
      REAL fqsol(klon,nbsrf)
      SAVE fqsol                  ! humidite du sol
c
      REAL fsnow(klon,nbsrf)
      SAVE fsnow                  ! epaisseur neigeuse
c
      REAL rugmer(klon)
      SAVE rugmer                 ! longeur de rugosite sur mer (m)
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
      REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
c
      REAL zuthe(klon),zvthe(klon)
      SAVE zuthe
      SAVE zvthe
      INTEGER igwd,igwdim,idx(klon),itest(klon)
c
      REAL agesno(klon)
      SAVE agesno                 ! age de la neige
c
      REAL alb_neig(klon)
      SAVE alb_neig               ! albedo de la neige
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 ema_work1(klon, klev), ema_work2(klon, klev)
      SAVE ema_work1, ema_work2
      REAL wdn(klon), tdn(klon), qdn(klon)
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
      LOGICAL offline           ! Controle du stockage ds "physique"
      PARAMETER (offline=.false.)
      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
      REAL          :: calday
#endif

cAA
      REAL rain_fall(klon) ! pluie
      REAL snow_fall(klon) ! neige
      REAL evap(klon), devap(klon) ! evaporation et sa derivee
      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
      REAL bils(klon) ! bilan de chaleur au sol
      REAL fder(klon) ! Derive de flux (sensible et latente) 
      REAL ruis(klon) ! ruissellement
      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
      REAL zxrugs(klon) ! longueur de rugosite
c
c Conditions aux limites
c
      INTEGER julien
      INTEGER idayvrai
      SAVE idayvrai
c
      INTEGER lmt_pas
      SAVE lmt_pas                ! frequence de mise a jour
      REAL pctsrf(klon,nbsrf)
      SAVE pctsrf                 ! sous-fraction du sol
      REAL lmt_sst(klon)
      SAVE lmt_sst                ! temperature de la surface ocean
      REAL lmt_bils(klon)
      SAVE lmt_bils               ! bilan de chaleur au sol
      REAL lmt_alb(klon)
      SAVE lmt_alb                ! temperature de la surface ocean
      REAL lmt_rug(klon)
      SAVE lmt_rug                ! longueur de rugosite
      REAL alb_eau(klon)
      SAVE alb_eau                ! albedo sur l'ocean
      REAL albsol(klon)
      SAVE albsol                 ! 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 albsno    ! calculer albedo sur neige
      EXTERNAL ajsec     ! ajustement sec
      EXTERNAL clmain    ! couche limite 
      EXTERNAL condsurf  ! lire les conditions aux limites
      EXTERNAL conlmd    ! convection (schema LMD)
cKE43
      EXTERNAL conema  ! 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
cAA
      EXTERNAL hgardfou  ! verifier les temperatures
      EXTERNAL hydrol    ! hydrologie du sol
      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
c
c Variables locales
c
      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
      REAL fluxq(klon,klev)   ! flux turbulent d'humidite
      REAL fluxt(klon,klev)   ! flux turbulent de chaleur
      REAL fluxu(klon,klev)   ! flux turbulent de vitesse u
      REAL fluxv(klon,klev)   ! flux turbulent de vitesse v
c
      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 topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
      REAL albpla(klon)
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
      SAVE  topsw0,toplw0,solsw0,sollw0, heat0, cool0
      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 zx_alb_lic, zx_alb_oce, zx_alb_ter, zx_alb_sic
      REAL zxtsol(klon), zxqsol(klon), zxsnow(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
      INTEGER i, k, iq, ig, j, nsrf, ll
      REAL t_coup
      PARAMETER (t_coup=234.0)
c
      REAL zphi(klon,klev)
      REAL zx_tmp_x(iim), zx_tmp_yjjmp1
      REAL zx_relief(iim,jjmp1)
      REAL zx_aire(iim,jjmp1)
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
      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_t_eva(klon,klev),d_q_eva(klon,klev)
      REAL rneb(klon,klev)
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)
      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 ratqs(klon,klev)
      integer flag_ratqs
      real zpt_conv(klon,klev)

c
c Variables liees a l'ecriture de la bande histoire physique
c
      INTEGER ecrit_mth
      SAVE ecrit_mth   ! frequence d'ecriture (fichier mensuel)
c
      INTEGER ecrit_day
      SAVE ecrit_day   ! frequence d'ecriture (fichier journalier)
c
      INTEGER ecrit_ins
      SAVE ecrit_ins   ! frequence d'ecriture (fichier instantane)
c
      INTEGER ecrit_reg
      SAVE ecrit_reg   ! frequence d'ecriture
c
      REAL oas_sols(klon), z_sols(iim,jjmp1)
      SAVE oas_sols
      REAL oas_nsol(klon), z_nsol(iim,jjmp1)
      SAVE oas_nsol
      REAL oas_rain(klon), z_rain(iim,jjmp1)
      SAVE oas_rain
      REAL oas_snow(klon), z_snow(iim,jjmp1)
      SAVE oas_snow
      REAL oas_evap(klon), z_evap(iim,jjmp1)
      SAVE oas_evap
      REAL oas_ruis(klon), z_ruis(iim,jjmp1)
      SAVE oas_ruis
      REAL oas_tsol(klon), z_tsol(iim,jjmp1)
      SAVE oas_tsol
      REAL oas_fder(klon), z_fder(iim,jjmp1)
      SAVE oas_fder
      REAL oas_albe(klon), z_albe(iim,jjmp1)
      SAVE oas_albe
      REAL oas_taux(klon), z_taux(iim,jjmp1)
      SAVE oas_taux
      REAL oas_tauy(klon), z_tauy(iim,jjmp1)
      SAVE oas_tauy
      REAL oas_ruisoce(klon), z_ruisoce(iim,jjmp1)
      SAVE oas_ruisoce
      REAL oas_ruisriv(klon), z_ruisriv(iim,jjmp1)
      SAVE oas_ruisriv
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)
      REAL u_seri(klon,klev), v_seri(klon,klev)
c
      REAL tr_seri(klon,klev,nbtr)
      REAL d_tr(klon,klev,nbtr)
      REAL source_tr(klon,nbtr)

      REAL zx_rh(klon,klev)
      REAL dtimeday,dtimecri,dtimexp9,fecri_pas,fecri86400,fecritday

      INTEGER        length
      PARAMETER    ( length = 100 )
      REAL tabcntr0( length       )
c
      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
      REAL zx_tmp_fi2d(klon)
      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
      SAVE nid_day, nid_mth, nid_ins
c
      INTEGER nhori, nvert
      REAL zsto, zout, zjulian
      integer idayref

      character*20 modname
      character*80 abort_message
      logical ok_sync

c
c Declaration des constantes et des fonctions thermodynamiques
c
#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"
c======================================================================
      modname = 'physiq'
      ok_sync=.TRUE.
      IF (nqmax .LT. 2) THEN
         PRINT*, 'eaux vapeur et liquide sont indispensables'
         CALL ABORT
      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
 
         IF (ok_oasis) THEN
            PRINT*, "Attentions! les parametres suivants sont fixes:"
            PRINT *,'***********************************************'
            PRINT*, "npas, nexca, itimestep=", npas, nexca, itimestep
            PRINT*, "Changer-les manuellement s il le faut"
            PRINT *,'***********************************************'
            CALL inicma( npas, nexca, itimestep)
         ENDIF
c
         IF (ok_ocean) THEN
            PRINT*, '************************'
            PRINT*, 'SLAB OCEAN est active, prenez precautions !'
            PRINT*, '************************'
         ENDIF
c
         DO k = 2, nvm          ! pas de vegetation
            DO i = 1, klon
               veget(i,k) = 0.0
            ENDDO
         ENDDO
         DO i = 1, klon
            veget(i,1) = 1.0    ! il n'y a que du desert
         ENDDO
         PRINT*, 'Pas de vegetation; desert partout'
c
c Initialiser les compteurs:
c

         itap    = 0
         itaprad = 0
c
         CALL phyetat0 ("startphy.nc",dtime,co2_ppm,solaire,
     .         rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
     .        radsol,rugmer,agesno,clesphy0,
     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
     .       t_ancien, q_ancien, ancien_ok )

c
         radpas = NINT( 86400./dtime/nbapp_rad)

c
         CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
     ,                    ok_instan, ok_region )
c
         IF (ABS(dtime-pdtphys).GT.0.001) THEN
            PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
            abort_message=' See above '
            call abort_gcm(modname,abort_message,1)
         ENDIF
         IF (nlon .NE. klon) THEN
            PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
            abort_message=' See above '
            call abort_gcm(modname,abort_message,1)
         ENDIF
         IF (nlev .NE. klev) THEN
            PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
            abort_message=' See above '
            call abort_gcm(modname,abort_message,1)
         ENDIF
c
         IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN 
           PRINT*, 'Nbre d appels au rayonnement insuffisant'
           PRINT*, "Au minimum 4 appels par jour si cycle diurne"
           abort_message=' See above '
           call abort_gcm(modname,abort_message,1)
         ENDIF
         PRINT*, "Clef pour la convection, iflag_con=", iflag_con
c
cKE43
c Initialisation pour la convection de K.E. (sb):
         IF (iflag_con.GE.3) THEN

         PRINT*, "*** Convection de Kerry Emanuel 4.3  "
         PRINT*, "On va utiliser le melange convectif des traceurs qui"
         PRINT*, "est calcule dans convect4.3"
         PRINT*, " !!! penser aux logical flags de phytrac"

          DO i = 1, klon
           ema_cbmf(i) = 0.
           ema_pcb(i)  = 0.
           ema_pct(i)  = 0.
           ema_workcbmf(i) = 0.
          ENDDO
         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
         IF (soil_model) THEN
            DO nsrf = 1, nbsrf
            CALL soil(dtime, nsrf, fsnow(1,nsrf),
     .                ftsol(1,nsrf), ftsoil(1,1,nsrf),
     .                soilcap(1,nsrf), soilflux(1,nsrf))
            ENDDO
         ENDIF
c
         lmt_pas = NINT(86400./dtime * 1.0)   ! tous les jours
         PRINT*,'La frequence de lecture surface est de ', lmt_pas
c
         ecrit_mth = NINT(86400./dtime *ecritphy)  ! tous les ecritphy jours
         IF (ok_mensuel) THEN
         PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
         ENDIF
         ecrit_day = NINT(86400./dtime *1.0)  ! tous les jours
         IF (ok_journe) THEN
         PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
         ENDIF
ccc         ecrit_ins = NINT(86400./dtime *0.5)  ! 2 fois par jour
         ecrit_ins = NINT(86400./dtime *0.25)  ! tous les jours
         IF (ok_instan) THEN
         PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
         ENDIF
         ecrit_reg = NINT(86400./dtime *0.25)  ! 4 fois par jour
         IF (ok_region) THEN
         PRINT*, 'La frequence de sortie region est de ', ecrit_reg
         ENDIF
c
c
      IF (ok_journe) THEN
c
         idayref = day_ini
         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
c
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
         DO i = 1, iim
            zx_lon(i,1) = rlon(i+1)
            zx_lon(i,jjmp1) = rlon(i+1)
         ENDDO
         DO ll=1,klev
            znivsig(ll)=float(ll)
         ENDDO
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
         CALL histbeg("histday", iim,zx_lon, jjmp1,zx_lat,
     .                 1,iim,1,jjmp1, 0, zjulian, dtime, 
     .                 nhori, nid_day)
         CALL histvert(nid_day, "presnivs", "Vertical levels", "mb",
     .                 klev, presnivs, nvert)
c        call histvert(nid_day, 'sig_s', 'Niveaux sigma','-',
c    .              klev, znivsig, nvert)
c
         zsto = dtime
         zout = dtime * FLOAT(ecrit_day)
c
         CALL histdef(nid_day, "phis", "Surface geop. height", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "once", zsto,zout)
c
         CALL histdef(nid_day, "aire", "Grid area", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "once", zsto,zout)
c
c Champs 2D:
c
         CALL histdef(nid_day, "tsol", "Surface Temperature", "K",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "psol", "Surface Pressure", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "rain", "Precipitation", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "snow", "Snow fall", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "evap", "Evaporation", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "tops", "Solar rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "topl", "IR rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "sols", "Solar rad. at surf.", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "soll", "IR rad. at surface", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "bils", "Surf. total heat flux", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "sens", "Sensible heat flux", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "fder", "Heat flux derivation", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "frtu", "Zonal wind stress", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "frtv", "Meridional wind stress", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "ruis", "Runoff", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "sicf", "Sea-ice fraction", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "cldl", "Low-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "cldm", "Mid-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "cldh", "High-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "cldt", "Total cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "cldq", "Cloud liquid water path", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
c Champs 3D:
c
         CALL histdef(nid_day, "temp", "Air temperature", "K",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "ovap", "Specific humidity", "Kg/Kg",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "geop", "Geopotential height", "m",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "vitu", "Zonal wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "vitv", "Meridional wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "vitw", "Vertical wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_day, "pres", "Air pressure", "Pa",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histend(nid_day)
c
         ndex2d = 0
         ndex3d = 0
c
      ENDIF ! fin de test sur ok_journe
c
      IF (ok_mensuel) THEN
c
         idayref = day_ini
         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
c
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
         DO i = 1, iim
            zx_lon(i,1) = rlon(i+1)
            zx_lon(i,jjmp1) = rlon(i+1)
         ENDDO
         DO ll=1,klev
            znivsig(ll)=float(ll)
         ENDDO
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
         CALL histbeg("histmth", iim,zx_lon, jjmp1,zx_lat,
     .                 1,iim,1,jjmp1, 0, zjulian, dtime, 
     .                 nhori, nid_mth)
         CALL histvert(nid_mth, "presnivs", "Vertical levels", "mb",
     .                 klev, presnivs, nvert)
c        call histvert(nid_mth, 'sig_s', 'Niveaux sigma','-',
c    .              klev, znivsig, nvert)
c
         zsto = dtime
         zout = dtime * FLOAT(ecrit_mth)
c
         CALL histdef(nid_mth, "phis", "Surface geop. height", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "once",  zsto,zout)
c
         CALL histdef(nid_mth, "aire", "Grid area", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "once",  zsto,zout)
c
c Champs 2D:
c
         CALL histdef(nid_mth, "tsol", "Surface Temperature", "K",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "psol", "Surface Pressure", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "qsol", "Surface humidity", "mm",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "rain", "Precipitation", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "plul", "Large-scale Precip.", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "pluc", "Convective Precip.", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "snow", "Snow fall", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ages", "Snow age", "day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "evap", "Evaporation", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "tops", "Solar rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "topl", "IR rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "sols", "Solar rad. at surf.", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "soll", "IR rad. at surface", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "tops0", "Solar rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "topl0", "IR rad. at TOA", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "sols0", "Solar rad. at surf.", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "soll0", "IR rad. at surface", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "bils", "Surf. total heat flux", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "fder", "Heat flux derivation", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "frtu", "Zonal wind stress", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "frtv", "Meridional wind stress", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ruis", "Runoff", "mm/day",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "sicf", "Sea-ice fraction", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "albs", "Surface albedo", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cdrm", "Momentum drag coef.", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cdrh", "Heat drag coef.", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cldl", "Low-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cldm", "Mid-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cldh", "High-level cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cldt", "Total cloudiness", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "cldq", "Cloud liquid water path", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ue", "Zonal energy transport", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ve", "Merid energy transport", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "uq", "Zonal humidity transport", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "vq", "Merid humidity transport", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32, 
     .                "ave(X)", zsto,zout)
cKE43
      IF (iflag_con .GE. 3) THEN ! sb
c
         CALL histdef(nid_mth, "cape", "Conv avlbl pot ener", "J/Kg",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "pbase", "Cld base pressure", "hPa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ptop", "Cld top pressure", "hPa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "fbase", "Cld base mass flux", "Kg/m2/s",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "ave(X)", zsto,zout)
c
c
      ENDIF
c34EK
c
c Champs 3D:
c
         CALL histdef(nid_mth, "temp", "Air temperature", "K",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "ovap", "Specific humidity", "Kg/Kg",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "geop", "Geopotential height", "m",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "vitu", "Zonal wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "vitv", "Meridional wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "vitw", "Vertical wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "pres", "Air pressure", "Pa",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "rneb", "Cloud fraction", "-",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "rhum", "Relative humidity", "-",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "oliq", "Liquid water content", "kg/kg",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtdyn", "Dynamics dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dqdyn", "Dynamics dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtcon", "Convection dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dqcon", "Convection dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtlsc", "Condensation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dqlsc", "Condensation dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtvdf", "Boundary-layer dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dqvdf", "Boundary-layer dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dteva", "Reevaporation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dqeva", "Reevaporation dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)

         CALL histdef(nid_mth, "ptconv", "POINTS CONVECTIFS"," ",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)

         CALL histdef(nid_mth, "ratqs", "RATQS"," ",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)

c
         CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)

         CALL histdef(nid_mth, "dqajs", "Dry adjust. dQ", "Kg/Kg/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtswr", "SW radiation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtsw0", "SW radiation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtlwr", "LW radiation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dtlw0", "LW radiation dT", "K/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "duvdf", "Boundary-layer dU", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dvvdf", "Boundary-layer dV", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         IF (ok_orodr) THEN
         CALL histdef(nid_mth, "duoro", "Orography dU", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dvoro", "Orography dV", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         ENDIF
C
         IF (ok_orolf) THEN
         CALL histdef(nid_mth, "dulif", "Orography dU", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dvlif", "Orography dV", "m/s2",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
         ENDIF
C
         CALL histdef(nid_mth, "ozone", "Ozone concentration", "-",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         if (nqmax.GE.3) THEN
         DO iq=1,nqmax-2
         IF (iq.LE.99) THEN
         WRITE(str2,'(i2.2)') iq
         CALL histdef(nid_mth, "trac"//str2, "Tracer No."//str2, "-",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
         ELSE
         PRINT*, "Trop de traceurs"
         CALL abort
         ENDIF
         ENDDO
         ENDIF
c
cKE43
      IF (iflag_con.GE.3) THEN ! (sb)
c
         CALL histdef(nid_mth, "upwd", "saturated updraft", "Kg/m2/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dnwd", "saturated downdraft","Kg/m2/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth, "dnwd0", "unsat. downdraft", "Kg/m2/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
         CALL histdef(nid_mth,"Ma","undilute adiab updraft","Kg/m2/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
c
c
      ENDIF
c34EK
         CALL histend(nid_mth)
c
         ndex2d = 0
         ndex3d = 0
c
      ENDIF ! fin de test sur ok_mensuel
c
c
      IF (ok_instan) THEN
c
         idayref = day_ini
         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
c
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
         DO i = 1, iim
            zx_lon(i,1) = rlon(i+1)
            zx_lon(i,jjmp1) = rlon(i+1)
         ENDDO
         DO ll=1,klev
            znivsig(ll)=float(ll)
         ENDDO
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
         CALL histbeg("histins", iim,zx_lon, jjmp1,zx_lat,
     .                 1,iim,1,jjmp1, 0, zjulian, dtime,
     .                 nhori, nid_ins)
         CALL histvert(nid_ins, "presnivs", "Vertical levels", "mb",
     .                 klev, presnivs, nvert)
c        call histvert(nid_ins, 'sig_s', 'Niveaux sigma','-',
c    .              klev, znivsig, nvert)
c
         zsto = dtime * ecrit_ins
         zout = dtime * ecrit_ins
C
         CALL histdef(nid_ins, "phis", "Surface geop. height", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "once", zsto,zout)
c
         CALL histdef(nid_ins, "aire", "Grid area", "-",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "once", zsto,zout)
c
c Champs 2D:
c
         CALL histdef(nid_ins, "psol", "Surface Pressure", "Pa",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histdef(nid_ins, "topl", "OLR", "W/m2",
     .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     .                "inst(X)", zsto,zout)
c
c Champs 3D:
c
         CALL histdef(nid_ins, "temp", "Temperature", "K",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histdef(nid_ins, "vitu", "Zonal wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histdef(nid_ins, "vitv", "Merid wind", "m/s",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histdef(nid_ins, "geop", "Geopotential height", "m",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histdef(nid_ins, "pres", "Air pressure", "Pa",
     .                iim,jjmp1,nhori, klev,1,klev,nvert, 32,
     .                "inst(X)", zsto,zout)
c
         CALL histend(nid_ins)
c
         ndex2d = 0
         ndex3d = 0
c
      ENDIF
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
         IF (ok_oasis) THEN
         DO i = 1, klon
           oas_sols(i) = 0.0
           oas_nsol(i) = 0.0
           oas_rain(i) = 0.0
           oas_snow(i) = 0.0
           oas_evap(i) = 0.0
           oas_ruis(i) = 0.0
           oas_tsol(i) = 0.0
           oas_fder(i) = 0.0
           oas_albe(i) = 0.0
           oas_taux(i) = 0.0
           oas_tauy(i) = 0.0
         ENDDO
         ENDIF
c

#ifdef INCA
           calday = FLOAT(MOD(NINT(xjour),360)) + gmtime
           print *, 'initial time ', xjour, calday
#ifdef INCAINFO
           PRINT *, 'Appel CHEMINI ...'
#endif
           CALL chemini( rpi,
     $                   rg,
     $                   ra,
     $                   paire,
     $                   rlat,
     $                   rlon,
     $                   calday,
     $                   tracnam,
     $                   natsnam,
     $                   mxoutflds,
     $                   outinst,
     $                   outtimav,
     $                   klon,
     $                   nqmax)
#ifdef INCAINFO
           PRINT *, 'OK.'
#endif
#endif

      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)
      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
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)
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
         idayvrai = NINT(xjour)
         PRINT *,' PHYS cond  julien ',julien,idayvrai,ok_limitvrai
         CALL condsurf(julien,idayvrai, pctsrf ,
     .                  lmt_sst,lmt_alb,lmt_rug,lmt_bils  )
         CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
      ENDIF
cccccccccc
      IF (ok_oasis .AND. MOD(itap-1,nexca).EQ.0) THEN
C
         CALL fromcpl(itap,jjmp1*iim,
     .        cpl_sst,cpl_sic,cpl_alb_sst,cpl_alb_sic)
         DO i = 1, iim-1 ! un seul point pour le pole nord
            cpl_sst(i,1) = cpl_sst(iim,1)
            cpl_sic(i,1) = cpl_sic(iim,1)
            cpl_alb_sst(i,1) = cpl_alb_sst(iim,1)
            cpl_alb_sic(i,1) = cpl_alb_sic(iim,1)
         ENDDO
         DO i = 2, iim ! un seul point pour le pole sud
            cpl_sst(i,jjmp1) = cpl_sst(1,jjmp1)
            cpl_sic(i,jjmp1) = cpl_sic(1,jjmp1)
            cpl_alb_sst(i,jjmp1) = cpl_alb_sst(1,jjmp1)
            cpl_alb_sic(i,jjmp1) = cpl_alb_sic(1,jjmp1)
         ENDDO
c
         ig = 1
         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
     .       pctsrf(ig,is_sic).GT.epsfra) THEN
            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
     .                        - (cpl_sic(1,1)-pctsrf(ig,is_sic))
            pctsrf(ig,is_sic) = cpl_sic(1,1)
            lmt_sst(ig) = cpl_sst(1,1)
         ENDIF
         DO j = 2, jjm
         DO i = 1, iim
         ig = ig + 1
         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
     .       pctsrf(ig,is_sic).GT.epsfra) THEN
           pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
     .                       - (cpl_sic(i,j)-pctsrf(ig,is_sic))
           pctsrf(ig,is_sic) = cpl_sic(i,j)
           lmt_sst(ig) = cpl_sst(i,j)
         ENDIF
         ENDDO
         ENDDO
         ig = ig + 1
         IF (pctsrf(ig,is_oce).GT.epsfra .OR.
     .       pctsrf(ig,is_sic).GT.epsfra) THEN
            pctsrf(ig,is_oce) = pctsrf(ig,is_oce)
     .                        - (cpl_sic(1,jjmp1)-pctsrf(ig,is_sic))
            pctsrf(ig,is_sic) = cpl_sic(1,jjmp1)
            lmt_sst(ig) = cpl_sst(1,jjmp1)
         ENDIF
c
      ENDIF ! ok_oasis
cccccccccc
c
 
      IF (ok_ocean) THEN
         DO i = 1, klon
            ftsol(i,is_oce) = lmt_sst(i) + deltat(i)
         ENDDO

      ELSE
         DO i = 1, klon
            ftsol(i,is_oce) = lmt_sst(i)
         ENDDO

      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))
         zlsdcp=RLSTT/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
c Appeler la diffusion verticale (programme de couche limite)
c
      DO i = 1, klon
         frugs(i,is_ter) = SQRT(lmt_rug(i)**2+rugoro(i)**2)
         frugs(i,is_lic) = rugoro(i)
         frugs(i,is_oce) = rugmer(i)
         frugs(i,is_sic) = 0.001
         zxrugs(i) = 0.0
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
      ENDDO
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
      ENDDO
      ENDDO
c
      CALL clmain(dtime,pctsrf,
     e            t_seri,q_seri,u_seri,v_seri,soil_model,
     e            ftsol,soilcap,soilflux,paprs,pplay,radsol,
     e            fsnow,fqsol,
     e            rlat, frugs,
     s            d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
     s            fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,rugmer,
     s            dsens, devap,
     s            ycoefh,yu1,yv1) 
c
      DO i = 1, klon
         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
         evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
         fder(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
c Incrementer la temperature du sol
c
      DO i = 1, klon
         zxtsol(i) = 0.0
         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
         ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
         zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
      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)
        ENDDO
      ENDDO

c
c Appeler le modele du sol
c
      IF (soil_model) THEN
         DO nsrf = 1, nbsrf
         CALL soil(dtime, nsrf, fsnow(1,nsrf),
     .             ftsol(1,nsrf), ftsoil(1,1,nsrf),
     .             soilcap(1,nsrf), soilflux(1,nsrf))
         ENDDO
      ENDIF
c
c Calculer la derive du flux infrarouge
c
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         fder(i) = fder(i) - 4.0*RSIGMA*zxtsol(i)**3 *
     .                       (ftsol(i,nsrf)-zxtsol(i))
     .                      *pctsrf(i,nsrf)
      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,paire)
         PRINT*, "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, fluxq(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)
      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:
          if (nqmax .GE. 4) then
              ntra = nbtr
          else
              ntra = 1 
          endif
          if (iflag_con.eq.4) then ! vectorise
          CALL conemav (dtime,paprs,pplay,t_seri,q_seri,
     .        u_seri,v_seri,tr_seri,nbtr,
     .        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,
c    .        Ma,cape,tvp,(/(nint(rflag(i)),i=1,size(rflag))/),
     .        Ma,cape,tvp,iflagctrl,
     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr)

          else

          CALL conema (dtime,paprs,pplay,t_seri,q_seri,
     .        u_seri,v_seri,tr_seri,nbtr,
     .        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)
          endif


          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
          PRINT*, "iflag_con non-prevu", iflag_con
          CALL abort
      ENDIF
      CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
     .              d_u_con, d_v_con)
      DO k = 1, klev
        DO i = 1, klon
          t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
          q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
          u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
          v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
        ENDDO
      ENDDO
      IF (check) THEN
          za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
          PRINT*, "aprescon=", za
          zx_t = 0.0
          za = 0.0
          DO i = 1, klon
            za = za + paire(i)/FLOAT(klon)
            zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
          ENDDO
          zx_t = zx_t/za*dtime
          PRINT*, "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
      IF (nqmax.GT.2) THEN !--melange convectif de traceurs
c
          IF (iflag_con .NE. 2 .AND.  iflag_con .LT. 3 ) THEN 
              PRINT*, 'Pour l instant, seul conflx fonctionne ',
     $            'avec traceurs', iflag_con
              PRINT*,' Mettre iflag_con', 
     $            ' = 2  ou 4 dans run.def et repasser'
              CALL abort
              ENDIF 
c
      ENDIF !--nqmax.GT.2
c
c Appeler l'ajustement sec
c
      CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
      DO k = 1, klev
      DO i = 1, klon
         t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
         q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
      ENDDO
      ENDDO

c   RATQS
      if (iflag_con.eq.2) then
          flag_ratqs=0
      else
          flag_ratqs=1
      endif
      call calcratqs (flag_ratqs ,
     I            paprs,pplay,q_seri,d_t_con,d_t_ajs
     O           ,ratqs,zpt_conv)
c
c Appeler le processus de condensation a grande echelle
c et le processus de precipitation
c
      CALL fisrtilp_tr(dtime,paprs,pplay,
     .           t_seri, q_seri,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)
      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,paire)
         PRINT*, "apresilp=", za
         zx_t = 0.0
         za = 0.0
         DO i = 1, klon
            za = za + paire(i)/FLOAT(klon)
            zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
        ENDDO
         zx_t = zx_t/za*dtime
         PRINT*, "Precip=", zx_t
      ENDIF
c
c Nuages diagnostiques:
c
      IF (iflag_con.EQ.2) THEN ! seulement pour Tiedtke
      CALL diagcld1(paprs,pplay,
     .             rain_con,snow_con,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
      ENDIF
c
c Nuages stratus artificiels:
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
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
      ENDDO
      ENDDO

#ifdef INCA
           calday = FLOAT(julien) + gmtime

#ifdef INCA_AER
      call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
     &   prfl,psfl,pctsrf,paire)
#endif

#ifdef INCAINFO
           PRINT *, 'Appel CHEMHOOK_BEGIN ...'
#endif
           CALL chemhook_begin (calday,
     $                          pctsrf(1,3),
     $                          rlat,
     $                          rlon,
     $                          paire,
     $                          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, 
     $                          tr_seri)
#ifdef INCAINFO
           PRINT *, 'OK.'
#endif
#endif

c
c Calculer les parametres optiques des nuages et quelques
c parametres pour diagnostiques:
c
      CALL nuage (paprs, pplay,
     .            t_seri, cldliq, cldfra, cldtau, cldemi,
     .            cldh, cldl, cldm, cldt, cldq)
c
c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
c
      IF (MOD(itaprad,radpas).EQ.0) THEN
      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)
c        CALL zenith(zlongi,gmtime,rlat,rlon,rmu0,fract) !va disparaitre
        CALL alboc_cd(rmu0,alb_eau)
      ELSE
        CALL angle(zlongi,rlat,fract,rmu0)
        CALL alboc(FLOAT(julien),rlat,alb_eau)
      ENDIF
      CALL albsno(veget,agesno,alb_neig)
      DO i = 1, klon
         zx_alb_oce = alb_eau(i)
         IF (pctsrf(i,is_oce).GT.epsfra .AND. ftsol(i,is_oce).LT.271.35)
     .   zx_alb_oce = 0.6 ! pour slab_ocean
         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_lic)/(fsnow(i,is_lic)+10.0)))
         zx_alb_lic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_ter)/(fsnow(i,is_ter)+10.0)))
         zx_alb_ter = alb_neig(i)*zfra + lmt_alb(i)*(1.0-zfra)
         zfra = MAX(0.0,MIN(1.0,fsnow(i,is_sic)/(fsnow(i,is_sic)+10.0)))
         zx_alb_sic = alb_neig(i)*zfra + 0.6*(1.0-zfra)
         albsol(i) = zx_alb_oce * pctsrf(i,is_oce)
     .             + zx_alb_lic * pctsrf(i,is_lic)
     .             + zx_alb_ter * pctsrf(i,is_ter)
     .             + zx_alb_sic * pctsrf(i,is_sic)
      ENDDO
      CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
     e            (dist, rmu0, fract, co2_ppm, solaire,
     e             paprs, pplay,zxtsol,albsol, t_seri,q_seri,wo,
     e             cldfra, cldemi, cldtau,
     s             heat,heat0,cool,cool0,radsol,albpla,
     s             topsw,toplw,solsw,sollw,
     s             topsw0,toplw0,solsw0,sollw0)
      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
c Calculer l'hydrologie de la surface
c
      CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, evap,
     .            agesno, ftsol,fqsol,fsnow, ruis)
c
      DO i = 1, klon
         zxqsol(i) = 0.0
         zxsnow(i) = 0.0
      ENDDO
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         zxqsol(i) = zxqsol(i) + fqsol(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
      DO nsrf = 1, nbsrf
      DO i = 1, klon
         IF (pctsrf(i,nsrf).LT.epsfra) THEN
            fqsol(i,nsrf) = zxqsol(i)
            fsnow(i,nsrf) = zxsnow(i)
         ENDIF
      ENDDO
      ENDDO
c
c Calculer le bilan du sol et la derive de temperature (couplage)
c
      DO i = 1, klon
         bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
      ENDDO
      IF (ok_ocean) THEN
         DO i = 1, klon
            cthermiq = cyang
            IF (ftsol(i,is_oce).LT. 271.35) cthermiq = cbing
            IF (pctsrf(i,is_oce).GT.epsfra) deltat(i) = deltat(i) +
     .                          (bils(i)-lmt_bils(i))/cthermiq * dtime
            IF (deltat(i).GT.15.0 ) deltat(i) = 15.0
            IF (deltat(i).LT.-15.0) deltat(i) = -15.0
         ENDDO
      ENDIF
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
        igwdim=MAX(1,igwd)
c
        CALL drag_noro(klon,klev,dtime,paprs,pplay,
     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
     e                   igwd,igwdim,idx,itest,
     e                   t_seri, u_seri, v_seri,
     s                   zulow, zvlow, zustr, zvstr,
     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
        igwdim=MAX(1,igwd)
c
        CALL lift_noro(klon,klev,dtime,paprs,pplay,
     e                   rlat,zmea,zstd, zsig, zgam, zthe,zpic,zval,
     e                   igwd,igwdim,idx,itest,
     e                   t_seri, u_seri, v_seri,
     s                   zulow, zvlow, zustr, zvstr,
     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
cAA
cAA Installation de l'interface online-offline pour traceurs
cAA
c====================================================================
c   Calcul  des tendances traceurs
c====================================================================
C Pascale : il faut quand meme apeller phytrac car il gere les sorties
cKE43       des traceurs => il faut donc mettre des flags a .false.
      IF (iflag_con.GE.3) THEN
c           on ajoute les tendances calculees par KE43
        DO iq=1, nqmax-2 ! Sandrine a -3 ???
        DO k = 1, nlev
        DO i = 1, klon
          tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
        ENDDO
        ENDDO
        WRITE(iqn,'(i2.2)') iq
        CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
        ENDDO
CMAF modif pour garder info du nombre de traceurs auxquels
C la physique s'applique
      ELSE
CMAF modif pour garder info du nombre de traceurs auxquels
C la physique s'applique
C
      call phytrac (rnpb,
     I                   itap, julien, gmtime, 
     I                   debut,lafin,
     I                   nqmax-2,
     I                   nlon,nlev,dtime,
     I                   u,v,t,paprs,pplay,
     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
     I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
     I                   frac_impa, frac_nucl,
     I                   rlon,paire,pphis,pphi,
     I                   albsol, 
     I                   qx(1,1,1), rhcl, 
     I                   cldfra, rneb, diafra, cldliq, itop_con,
     I                   ibas_con,
     I                   pmflxr,pmflxs,prfl,psfl,flxmass_w,
     O                   tr_seri)
      ENDIF 

      IF (offline) THEN

	 call phystokenc (
     I                   nlon,nlev,pdtphys,rlon,rlat,
     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
     I                   ycoefh,yu1,yv1,ftsol,pctsrf,
     I                   frac_impa, frac_nucl,
     I                   pphis,paire,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
c Accumuler les variables a stocker dans les fichiers histoire:
c
      IF (ok_oasis) THEN ! couplage oasis
      DO i = 1, klon
        oas_sols(i) = oas_sols(i) + solsw(i)          / FLOAT(nexca)
        oas_nsol(i) = oas_nsol(i) + (bils(i)-solsw(i))/ FLOAT(nexca)
        oas_rain(i) = oas_rain(i) + rain_fall(i)      / FLOAT(nexca)
        oas_snow(i) = oas_snow(i) + snow_fall(i)      / FLOAT(nexca)
        oas_evap(i) = oas_evap(i) + evap(i)           / FLOAT(nexca)
        oas_tsol(i) = oas_tsol(i) + zxtsol(i)         / FLOAT(nexca)
        oas_fder(i) = oas_fder(i) + fder(i)           / FLOAT(nexca)
        oas_albe(i) = oas_albe(i) + albsol(i)         / FLOAT(nexca)
        oas_taux(i) = oas_taux(i) + fluxu(i,1)        / FLOAT(nexca)
        oas_tauy(i) = oas_tauy(i) + fluxv(i,1)        / FLOAT(nexca)
        oas_ruis(i) = oas_ruis(i) + ruis(i)       /FLOAT(nexca)/dtime
      ENDDO
      ENDIF
c
c
      IF (ok_journe) THEN
c
      ndex2d = 0
      ndex3d = 0
c
c Champs 2D:
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
         CALL histwrite(nid_day,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
         CALL histwrite(nid_day,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)
C
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
      CALL histwrite(nid_day,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = paprs(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_day,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_day,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
      CALL histwrite(nid_day,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
      CALL histwrite(nid_day,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
      CALL histwrite(nid_day,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
      CALL histwrite(nid_day,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
      CALL histwrite(nid_day,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
      CALL histwrite(nid_day,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
      CALL histwrite(nid_day,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
      CALL histwrite(nid_day,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
      CALL histwrite(nid_day,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
      CALL histwrite(nid_day,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = fluxu(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_day,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = fluxv(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_day,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_day,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
      CALL histwrite(nid_day,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
      CALL histwrite(nid_day,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
      CALL histwrite(nid_day,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
      CALL histwrite(nid_day,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
      CALL histwrite(nid_day,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
c Champs 3D:
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
      CALL histwrite(nid_day,"temp",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
      CALL histwrite(nid_day,"ovap",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
      CALL histwrite(nid_day,"geop",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
      CALL histwrite(nid_day,"vitu",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
      CALL histwrite(nid_day,"vitv",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
      CALL histwrite(nid_day,"vitw",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
      CALL histwrite(nid_day,"pres",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      if (ok_sync) then
        call histsync(nid_day)
      endif
      ENDIF
C
      IF (ok_mensuel) THEN
c
      ndex2d = 0
      ndex3d = 0
c
c Champs 2D:
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
         CALL histwrite(nid_mth,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
C
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
         CALL histwrite(nid_mth,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)

      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
      CALL histwrite(nid_mth,"tsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = paprs(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxqsol,zx_tmp_2d)
      CALL histwrite(nid_mth,"qsol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"rain",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"plul",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"pluc",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, snow_fall,zx_tmp_2d)
      CALL histwrite(nid_mth,"snow",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, agesno,zx_tmp_2d)
      CALL histwrite(nid_mth,"ages",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, evap,zx_tmp_2d)
      CALL histwrite(nid_mth,"evap",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw,zx_tmp_2d)
      CALL histwrite(nid_mth,"tops",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
      CALL histwrite(nid_mth,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw,zx_tmp_2d)
      CALL histwrite(nid_mth,"sols",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw,zx_tmp_2d)
      CALL histwrite(nid_mth,"soll",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, topsw0,zx_tmp_2d)
      CALL histwrite(nid_mth,"tops0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw0,zx_tmp_2d)
      CALL histwrite(nid_mth,"topl0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, solsw0,zx_tmp_2d)
      CALL histwrite(nid_mth,"sols0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sollw0,zx_tmp_2d)
      CALL histwrite(nid_mth,"soll0",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, bils,zx_tmp_2d)
      CALL histwrite(nid_mth,"bils",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, sens,zx_tmp_2d)
      CALL histwrite(nid_mth,"sens",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, fder,zx_tmp_2d)
      CALL histwrite(nid_mth,"fder",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ruis,zx_tmp_2d)
      CALL histwrite(nid_mth,"ruis",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = fluxu(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"frtu",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = fluxv(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"frtv",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      DO i = 1, klon
         zx_tmp_fi2d(i) = pctsrf(i,is_sic)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_mth,"sicf",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, albsol,zx_tmp_2d)
      CALL histwrite(nid_mth,"albs",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragm,zx_tmp_2d)
      CALL histwrite(nid_mth,"cdrm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cdragh,zx_tmp_2d)
      CALL histwrite(nid_mth,"cdrh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldl,zx_tmp_2d)
      CALL histwrite(nid_mth,"cldl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldm,zx_tmp_2d)
      CALL histwrite(nid_mth,"cldm",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldh,zx_tmp_2d)
      CALL histwrite(nid_mth,"cldh",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldt,zx_tmp_2d)
      CALL histwrite(nid_mth,"cldt",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cldq,zx_tmp_2d)
      CALL histwrite(nid_mth,"cldq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ue,zx_tmp_2d)
      CALL histwrite(nid_mth,"ue",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, ve,zx_tmp_2d)
      CALL histwrite(nid_mth,"ve",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, uq,zx_tmp_2d)
      CALL histwrite(nid_mth,"uq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, vq,zx_tmp_2d)
      CALL histwrite(nid_mth,"vq",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
cKE43
      IF (iflag_con .GE. 3) THEN ! sb
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, cape,zx_tmp_2d)
      CALL histwrite(nid_mth,"cape",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1,pbase,zx_tmp_2d)
      CALL histwrite(nid_mth,"pbase",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_pct,zx_tmp_2d)
      CALL histwrite(nid_mth,"ptop",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1,ema_cbmf,zx_tmp_2d)
      CALL histwrite(nid_mth,"fbase",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
c
      ENDIF
c34EK
c
c Champs 3D:
C
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
      CALL histwrite(nid_mth,"temp",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,ivap), zx_tmp_3d)
      CALL histwrite(nid_mth,"ovap",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
      CALL histwrite(nid_mth,"geop",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
      CALL histwrite(nid_mth,"vitu",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
      CALL histwrite(nid_mth,"vitv",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, omega, zx_tmp_3d)
      CALL histwrite(nid_mth,"vitw",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
      CALL histwrite(nid_mth,"pres",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldfra, zx_tmp_3d)
      CALL histwrite(nid_mth,"rneb",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zx_rh, zx_tmp_3d)
      CALL histwrite(nid_mth,"rhum",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cldliq, zx_tmp_3d)
      CALL histwrite(nid_mth,"oliq",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_dyn, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtdyn",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_dyn, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqdyn",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_con, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtcon",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_con, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqcon",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_lsc, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtlsc",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_lsc, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqlsc",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_vdf, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtvdf",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_vdf, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqvdf",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_eva, zx_tmp_3d)
      CALL histwrite(nid_mth,"dteva",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_eva, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqeva",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zpt_conv, zx_tmp_3d)
      CALL histwrite(nid_mth,"ptconv",itap,zx_tmp_3d,
     .                                   iim*(jjmp1)*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, ratqs, zx_tmp_3d)
      CALL histwrite(nid_mth,"ratqs",itap,zx_tmp_3d,
     .                                   iim*(jjmp1)*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_t_ajs, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtajs",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_q_ajs, zx_tmp_3d)
      CALL histwrite(nid_mth,"dqajs",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtswr",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, heat0, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtsw0",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtlwr",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, cool0, zx_tmp_3d)
      CALL histwrite(nid_mth,"dtlw0",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_vdf, zx_tmp_3d)
      CALL histwrite(nid_mth,"duvdf",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_vdf, zx_tmp_3d)
      CALL histwrite(nid_mth,"dvvdf",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      IF (ok_orodr) THEN
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_oro, zx_tmp_3d)
      CALL histwrite(nid_mth,"duoro",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_oro, zx_tmp_3d)
      CALL histwrite(nid_mth,"dvoro",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      ENDIF
C
      IF (ok_orolf) THEN
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_u_lif, zx_tmp_3d)
      CALL histwrite(nid_mth,"dulif",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, d_v_lif, zx_tmp_3d)
      CALL histwrite(nid_mth,"dvlif",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
      ENDIF
C
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, wo, zx_tmp_3d)
      CALL histwrite(nid_mth,"ozone",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      IF (nqmax.GE.3) THEN
      DO iq=1,nqmax-2
      IF (iq.LE.99) THEN
         CALL gr_fi_ecrit(klev,klon,iim,jjmp1, qx(1,1,iq+2), zx_tmp_3d)
         WRITE(str2,'(i2.2)') iq
         CALL histwrite(nid_mth,"trac"//str2,itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
      ELSE
         PRINT*, "Trop de traceurs"
         CALL abort
      ENDIF
      ENDDO
      ENDIF
cKE43
      IF (iflag_con.GE.3) THEN ! (sb)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, upwd, zx_tmp_3d)
      CALL histwrite(nid_mth,"upwd",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd, zx_tmp_3d)
      CALL histwrite(nid_mth,"dnwd",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, dnwd0, zx_tmp_3d)
      CALL histwrite(nid_mth,"dnwd0",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, Ma, zx_tmp_3d)
      CALL histwrite(nid_mth,"Ma",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
c
      ENDIF
c34EK
c
      if (ok_sync) then
        call histsync(nid_mth)
      endif
      ENDIF
c
      IF (ok_instan) THEN
c
      ndex2d = 0
      ndex3d = 0
c
c Champs 2D:
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,pphis,zx_tmp_2d)
         CALL histwrite(nid_ins,"phis",i,zx_tmp_2d,iim*jjmp1,ndex2d)
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjmp1,paire,zx_tmp_2d)
         CALL histwrite(nid_ins,"aire",i,zx_tmp_2d,iim*jjmp1,ndex2d)

      DO i = 1, klon
         zx_tmp_fi2d(i) = paprs(i,1)
      ENDDO
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zx_tmp_fi2d,zx_tmp_2d)
      CALL histwrite(nid_ins,"psol",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
      CALL gr_fi_ecrit(1, klon,iim,jjmp1, toplw,zx_tmp_2d)
      CALL histwrite(nid_ins,"topl",itap,zx_tmp_2d,iim*jjmp1,ndex2d)
c
c Champs 3D:
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, t_seri, zx_tmp_3d)
      CALL histwrite(nid_ins,"temp",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, u_seri, zx_tmp_3d)
      CALL histwrite(nid_ins,"vitu",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, v_seri, zx_tmp_3d)
      CALL histwrite(nid_ins,"vitv",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, zphi, zx_tmp_3d)
      CALL histwrite(nid_ins,"geop",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, pplay, zx_tmp_3d)
      CALL histwrite(nid_ins,"pres",itap,zx_tmp_3d,
     .                                   iim*jjmp1*klev,ndex3d)
c
      if (ok_sync) then
        call histsync(nid_ins)
      endif
      ENDIF
c
      IF (ok_oasis .AND. mod(itap,nexca).EQ.0) THEN
c
c Je ne traite pas le ruissellement, pour l'instant (qui m'aidera ?)
         DO i = 1, klon
            oas_ruisoce(i) = 0.0
            oas_ruisriv(i) = 0.0
         ENDDO
c
         ig = 1
         DO i = 1, iim
            z_sols(i,1) = oas_sols(ig)
            z_nsol(i,1) = oas_nsol(ig)
            z_rain(i,1) = oas_rain(ig)
            z_snow(i,1) = oas_snow(ig)
            z_evap(i,1) = oas_evap(ig)
            z_ruisoce(i,1) = oas_ruisoce(ig)
            z_ruisriv(i,1) = oas_ruisriv(ig)
            z_tsol(i,1) = oas_tsol(ig)
            z_fder(i,1) = oas_fder(ig)
            z_albe(i,1) = oas_albe(ig)
            z_taux(i,1) = oas_taux(ig)
            z_tauy(i,1) = oas_tauy(ig)
         ENDDO
         DO j = 2, jjm
         DO i = 1, iim
            ig = ig + 1
            z_sols(i,j) = oas_sols(ig)
            z_nsol(i,j) = oas_nsol(ig)
            z_rain(i,j) = oas_rain(ig)
            z_snow(i,j) = oas_snow(ig)
            z_evap(i,j) = oas_evap(ig)
            z_ruisoce(i,j) = oas_ruisoce(ig)
            z_ruisriv(i,j) = oas_ruisriv(ig)
            z_tsol(i,j) = oas_tsol(ig)
            z_fder(i,j) = oas_fder(ig)
            z_albe(i,j) = oas_albe(ig)
            z_taux(i,j) = oas_taux(ig)
            z_tauy(i,j) = oas_tauy(ig)
         ENDDO
         ENDDO
         ig = ig + 1
         DO i = 1, iim
            z_sols(i,jjmp1)    = oas_sols(ig)
            z_nsol(i,jjmp1)    = oas_nsol(ig)
            z_rain(i,jjmp1)    = oas_rain(ig)
            z_snow(i,jjmp1)    = oas_snow(ig)
            z_evap(i,jjmp1)    = oas_evap(ig)
            z_ruisoce(i,jjmp1) = oas_ruisoce(ig)
            z_ruisriv(i,jjmp1) = oas_ruisriv(ig)
            z_tsol(i,jjmp1)    = oas_tsol(ig)
            z_fder(i,jjmp1)    = oas_fder(ig)
            z_albe(i,jjmp1)    = oas_albe(ig)
            z_taux(i,jjmp1)    = oas_taux(ig)
            z_tauy(i,jjmp1)    = oas_tauy(ig)
         ENDDO
c
c Passer les champs au coupleur:
c
         CALL intocpl(itap,jjmp1*iim,
     .                   z_sols, z_nsol,
     .                   z_rain, z_snow, z_evap,
     .                   z_ruisoce, z_ruisriv,
     .                   z_tsol, z_fder, z_albe,
     .                   z_taux, z_tauy)
         DO i = 1, klon
           oas_sols(i) = 0.0
           oas_nsol(i) = 0.0
           oas_rain(i) = 0.0
           oas_snow(i) = 0.0
           oas_evap(i) = 0.0
           oas_ruis(i) = 0.0
           oas_tsol(i) = 0.0
           oas_fder(i) = 0.0
           oas_albe(i) = 0.0
           oas_taux(i) = 0.0
           oas_tauy(i) = 0.0
         ENDDO
      ENDIF
c
c Ecrire la bande regionale (binaire grads)
      IF (ok_region .AND. mod(itap,ecrit_reg).eq.0) THEN
         CALL ecriregs(84,zxtsol)
         CALL ecriregs(84,paprs(1,1))
         CALL ecriregs(84,topsw)
         CALL ecriregs(84,toplw)
         CALL ecriregs(84,solsw)
         CALL ecriregs(84,sollw)
         CALL ecriregs(84,rain_fall)
         CALL ecriregs(84,snow_fall)
         CALL ecriregs(84,evap)
         CALL ecriregs(84,sens)
         CALL ecriregs(84,bils)
         CALL ecriregs(84,pctsrf(1,is_sic))
         CALL ecriregs(84,fluxu(1,1))
         CALL ecriregs(84,fluxv(1,1))
         CALL ecriregs(84,ue)
         CALL ecriregs(84,ve)
         CALL ecriregs(84,uq)
         CALL ecriregs(84,vq)
c
         CALL ecrirega(84,u_seri)
         CALL ecrirega(84,v_seri)
         CALL ecrirega(84,omega)
         CALL ecrirega(84,t_seri)
         CALL ecrirega(84,zphi)
         CALL ecrirega(84,q_seri)
         CALL ecrirega(84,cldfra)
         CALL ecrirega(84,cldliq)
         CALL ecrirega(84,pplay)


cc         CALL ecrirega(84,d_t_dyn)
cc         CALL ecrirega(84,d_q_dyn)
cc         CALL ecrirega(84,heat)
cc         CALL ecrirega(84,cool)
cc         CALL ecrirega(84,d_t_con)
cc         CALL ecrirega(84,d_q_con)
cc         CALL ecrirega(84,d_t_lsc)
cc         CALL ecrirega(84,d_q_lsc)
      ENDIF

#ifdef INCA
#ifdef INCAINFO
           PRINT *, 'Appel CHEMHOOK_END ...'
#endif
           CALL chemhook_end (calday,
     $                        dtime,
     $                        pplay,
     $                        t_seri,
     $                        tr_seri,
     $                        nbtr,
     $                        paprs,
     $                        anne_ini,
     $                        day_ini,
     $                        xjour)
#ifdef INCAINFO
           PRINT *, 'OK.'
#endif
#endif

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
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====================================================================
c Si c'est la fin, il faut conserver l'etat de redemarrage
c====================================================================
c
      IF (lafin) THEN
ccc         IF (ok_oasis) CALL quitcpl
         CALL phyredem ("restartphy.nc",dtime,radpas,co2_ppm,solaire,
     .        rlat,rlon,ftsol,ftsoil,deltat,fqsol,fsnow,
     .        radsol,rugmer,agesno,
     .        zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
     .        t_ancien, q_ancien)
      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, j, 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