!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
!
c
      SUBROUTINE physiq (nlon,nlev,nqmax,
     .            debut,lafin,rjourvrai,gmtime,pdtphys,
     .            paprs,pplay,ppk,pphi,pphis,presnivs,
     .            u,v,t,qx,
     .            omega,
     .            d_u, d_v, d_t, d_qx, d_ps)

c======================================================================
c
c Modifications pour la physique de Titan
c  adaptation a partir de celle de Venus
c     S. Lebonnois (LMD/CNRS) Mai 2008
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    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
c debut-----input-L-variable logique indiquant le premier passage
c lafin-----input-L-variable logique indiquant le dernier passage
c rjourvrai-input-R-NBjours
c gmtime----input-R-temps universel dans la journee (fraction de jour)
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 ppk  -----input-R-fonction d'Exner au milieu de couche
c pphi------input-R-geopotentiel de chaque couche (g z) (reference sol)
c pphis-----input-R-geopotentiel du sol
c presnivs--input_R_pressions approximat. des milieux couches ( en PA)
c u---------input-R-vitesse dans la direction X (de O a E) en m/s
c v---------input-R-vitesse Y (de S a N) en m/s
c t---------input-R-temperature (K)
c qx--------input-R-mass mixing ratio traceurs (kg/kg) 
c d_t_dyn---input-R-tendance dynamique pour "t" (K/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======================================================================
      USE ioipsl
!      USE histcom ! not needed; histcom is included in ioipsl
      USE infotrac
      USE control_mod
      use dimphy
      USE comgeomphy
      IMPLICIT none
c======================================================================
c   CLEFS CPP POUR LES IO
c   =====================
#define histmth
#define histday
#define histins
c======================================================================
#include "dimensions.h"
      integer jjmp1
      parameter (jjmp1=jjm+1-1/jjm)
#include "dimsoil.h"
#include "clesphys.h"
#include "temps.h"
#include "iniprint.h"
#include "logic.h"
#include "comorbit.h"
#include "microtab.h"
#include "diagmuphy.h"
#include "tabcontrol.h"
#include "itemps.h"
c======================================================================
      LOGICAL ok_mensuel ! sortir le fichier mensuel
      save ok_mensuel
c      PARAMETER (ok_mensuel=.true.)
c
      LOGICAL ok_journe ! sortir le fichier journalier
      save ok_journe
c      PARAMETER (ok_journe=.true.)
c
      LOGICAL ok_instan ! sortir le fichier instantane
      save ok_instan
c      PARAMETER (ok_instan=.true.)
c
c======================================================================
c
c Variables argument:
c
      INTEGER nlon
      INTEGER nlev
      INTEGER nqmax
      REAL rjourvrai
      REAL gmtime
      REAL pdtphys
      LOGICAL debut, lafin
      REAL paprs(klon,klev+1)
      REAL pplay(klon,klev)
      REAL pphi(klon,klev)
      REAL pphis(klon)
      REAL presnivs(klev)

! ADAPTATION GCM POUR CP(T)
      REAL ppk(klon,klev)

      REAL u(klon,klev)
      REAL v(klon,klev)
      REAL t(klon,klev)
      REAL qx(klon,klev,nqmax)

      REAL,save,allocatable :: t_ancien(:,:)
      REAL,save,allocatable :: u_ancien(:,:)
      LOGICAL ancien_ok
      SAVE ancien_ok

      REAL d_u_dyn(klon,klev)
      REAL d_t_dyn(klon,klev)

      REAL omega(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)

      REAL,save,allocatable :: swnet(:,:)
      REAL,save,allocatable :: lwnet(:,:)
c
c Variables propres a la physique
c
      REAL,save,allocatable :: radsol(:) ! bilan radiatif au sol calcule par code radiatif
      REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche) 
      INTEGER,save :: itap        ! compteur pour la physique
      REAL,save,allocatable :: ftsol(:)    ! temperature du sol
      REAL,save,allocatable :: ftsoil(:,:) ! temperature dans le sol
      REAL,save,allocatable :: falbe(:)    ! albedo 
      REAL delp(klon,klev)        ! epaisseur d'une couche
      
CMODDEB FLOTT
c
c  Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
c
      REAL,save,allocatable :: zmea(:)   ! orographie moyenne
      REAL,save,allocatable :: zstd(:)   ! deviation standard de l'OESM
      REAL,save,allocatable :: zsig(:)   ! pente de l'OESM
      REAL,save,allocatable :: zgam(:)   ! anisotropie de l'OESM
      REAL,save,allocatable :: zthe(:)   ! orientation de l'OESM
      REAL,save,allocatable :: zpic(:)   ! Maximum de l'OESM
      REAL,save,allocatable :: zval(:)   ! Minimum de l'OESM
      REAL,save,allocatable :: rugoro(:) ! longueur de rugosite de l'OESM

      INTEGER igwd,idx(klon),itest(klon)
c
c  Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro

      REAL zulow(klon),zvlow(klon)
      REAL zustrdr(klon), zvstrdr(klon)
      REAL zustrli(klon), zvstrli(klon)
      REAL zustrhi(klon), zvstrhi(klon)

c Pour calcul GW drag oro et nonoro: CALCUL de N2:
      real zdzlev(klon,klev)
      real ztlev(klon,klev),zpklev(klon,klev)
      real ztetalay(klon,klev),ztetalev(klon,klev)
      real zdtetalev(klon,klev)
      real zn2(klon,klev) ! BV^2 at plev

c Pour les bilans de moment angulaire, 
      integer bilansmc
c Pour le transport de ballons
      integer ballons
c j'ai aussi besoin
c du stress de couche limite a la surface:

      REAL zustrcl(klon),zvstrcl(klon) 

c et du stress total c de la physique:

      REAL zustrph(klon),zvstrph(klon) 
c
      REAL,save,allocatable :: zuthe(:),zvthe(:)

c Variables locales:
c
      REAL cdragh(klon) ! drag coefficient pour T and Q
      REAL cdragm(klon) ! drag coefficient pour vent
c
cAA  Pour TRACEURS 
cAA
      REAL,save,allocatable :: source(:,:)
      integer nmicro
      save    nmicro

      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
      character*8 nom
      REAL qaer(klon,klev,nqmax)

      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
      REAL,save,allocatable :: dlw(:)  ! derivee infra rouge
      REAL,save,allocatable :: fder(:) ! Derive de flux (sensible et latente) 
      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

c======================================================================
c
c Declaration des procedures appelees
c
      EXTERNAL ajsec     ! ajustement sec
      EXTERNAL clmain    ! couche limite 
      EXTERNAL hgardfou  ! verifier les temperatures
      EXTERNAL orbite    ! calculer l'orbite 
      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
c      EXTERNAL transp    ! transport total de l'eau et de l'energie
      EXTERNAL abort_gcm
      EXTERNAL printflag
      EXTERNAL zenang
      EXTERNAL diagetpq
      EXTERNAL conf_phys
      EXTERNAL diagphy
      EXTERNAL mucorr
      EXTERNAL phytrac
c
c Variables locales
c
CXXX PB 
      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 flux_dyn(klon,klev)  ! flux de chaleur produit par la dynamique
      REAL flux_ajs(klon,klev)  ! flux de chaleur ajustement sec
      REAL flux_ec(klon,klev)   ! flux de chaleur Ec
c
c Le rayonnement n'est pas calcule tous les pas, il faut donc
c                      sauvegarder les sorties du rayonnement
      REAL,save,allocatable :: heat(:,:)    ! chauffage solaire
      REAL,save,allocatable :: cool(:,:)    ! refroidissement infrarouge
      REAL,save,allocatable :: dtrad(:,:)   ! K s-1 
      REAL,save,allocatable :: topsw(:), toplw(:)
      REAL,save,allocatable :: solsw(:), sollw(:)
      REAL,save,allocatable :: sollwdown(:) ! downward LW flux at surface
      REAL tmpout(klon,klev)  ! K s-1 

      REAL    dtimerad
      INTEGER itaprad
      SAVE itaprad,dtimerad
      REAL zdtime
c

c CHIMIE

      REAL    dtimechim
      INTEGER itapchim,appel_chim
      SAVE itapchim,dtimechim

c ORBITE

      REAL dist, rmu0(klon), fract(klon), pdecli
      REAL zday
      REAL zls,zlsdeg,zlsm1
c
      INTEGER i, k, iq, ig, j, ll, l
c
      REAL zphi(klon,klev)
      REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
c
c Variables du changement
c
c ajs: ajustement sec
c vdf: couche limite (Vertical DiFfusion)
c mph: microphysique
c kim: chimie
      REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
c
      REAL d_ts(klon)
c
      REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
      REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
c
      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)

CMOD LOTT: Tendances Orography Sous-maille
      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)
C          Tendances Ondes de G non oro (runs strato).
      REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
      REAL d_t_hin(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 itau_w   ! pas de temps ecriture = itap + itau_phy

c Variables locales pour effectuer les appels en serie
c
      REAL t_seri(klon,klev)
      REAL u_seri(klon,klev), v_seri(klon,klev)
c
      REAL tr_seri(klon,klev,nqmax)
c
c pour ioipsl
      INTEGER ndex2d(iim*jjmp1),ndex3d(iim*jjmp1*klev)
      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
      REAL zx_tmp_2d(iim,jjmp1),zx_tmp_3d(iim,jjmp1,klev)
      REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)

      INTEGER nid_day, nid_mth, nid_ins
      SAVE nid_day, nid_mth, nid_ins
c
      INTEGER nhori, nvert, idayref
      REAL zsto, zout, zsto1, zsto2, zero
      parameter (zero=0.0e0)
      real zjulian
      save zjulian

      CHARACTER*2  str2
      character*20 modname
      character*80 abort_message
      logical ok_sync

      character*30 nom_fichier
      character*10 varname
      character*40 vartitle
      character*20 varunits
C     Variables liees au bilan d'energie et d'enthalpi
      REAL ztsol(klon)
      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
      REAL      d_h_vcol_phy
      REAL      fs_bound, fq_bound
      SAVE      d_h_vcol_phy
      REAL      zero_v(klon),zero_v2(klon,klev)
      CHARACTER*15 ztit
      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
      SAVE      ip_ebil
      DATA      ip_ebil/2/
      INTEGER   if_ebil ! level for energy conserv. dignostics
      SAVE      if_ebil
c+jld ec_conser
      REAL d_t_ec(klon,klev)    ! tendance du a la conversion Ec -> E thermique
c-jld ec_conser

c TEST VENUS...
      REAL mang(klon,klev)    ! moment cinetique
      REAL mangtot            ! moment cinetique total

      CHARACTER*2 str1

c Temporaire avant de trouver mieux :
c Recuperation des TAU du TR
      REAL t_tauhvd(klon,klev),t_khvd(klon,klev)
      REAL t_tcld(klon,klev),t_kcld(klon,klev)
      REAL t_kcvd(klon,klev)
c  ASTUCE POUR EVITER klon... EN ATTENDANT MIEUX
      INTEGER   ngrid
      PARAMETER (ngrid=(jjm-1)*iim+2)  ! = klon
      INTEGER NSPECV,NSPECI,NLAYER
      PARAMETER (NSPECV=24,NSPECI=46,NLAYER=llm)
      REAL TAUHID(ngrid,NLAYER,NSPECI)
     &               ,TAUCID(ngrid,NLAYER,NSPECI)
     &               ,TAUGID(ngrid,NLAYER,NSPECI)
     &               ,TAUHVD(ngrid,NLAYER,NSPECV)
     &               ,TAUCVD(ngrid,NLAYER,NSPECV)
     &               ,TAUGVD(ngrid,NLAYER,NSPECV)

      COMMON /TAUD/   TAUHID,TAUCID,TAUGID,TAUHVD,TAUCVD,TAUGVD
      COMMON /PLANT/ CSUBP,F0PI
      REAL CSUBP,F0PI

* common relatifs au nuages
      real rmcbar(ngrid,NLAYER),xfbar(ngrid,NLAYER,4) 
      integer ncount(ngrid,NLAYER)
      COMMON/rnuabar/ncount,rmcbar,xfbar

       REAL ch4(klon,jjm+1),dch4(jjm+1)
       INTEGER ig0
       integer ich4
       common/ch4ind/ich4

c     flux de chaleur latente d'evaporation CH4
      REAL fclat(klon)
c     reservoir de surface
      REAL,save,allocatable :: reservoir(:)

c Declaration des constantes et des fonctions thermodynamiques
c
#include "YOMCST.h"

c======================================================================
c INITIALISATIONS
c======================================================================

      modname = 'physiq'
      ok_sync=.TRUE.

      bilansmc = 0
      ballons  = 0

      IF (if_ebil.ge.1) THEN
        DO i=1,klon
          zero_v(i)=0.
        END DO 
        DO i=1,klon
         DO j=1,klev
          zero_v2(i,j)=0.
         END DO 
        END DO 
      END IF 
      
c PREMIER APPEL SEULEMENT
c========================
      IF (debut) THEN
         allocate(t_ancien(klon,klev),u_ancien(klon,klev))
         allocate(swnet(klon,klevp1),lwnet(klon,klevp1))
         allocate(radsol(klon),ftsol(klon),falbe(klon))
         allocate(rlev(klon,klevp1),ftsoil(klon,nsoilmx))
         allocate(zmea(klon),zstd(klon),zsig(klon),zgam(klon))
         allocate(zthe(klon),zpic(klon),zval(klon),rugoro(klon))
         allocate(zuthe(klon),zvthe(klon),dlw(klon),fder(klon))
         allocate(heat(klon,klev),cool(klon,klev))
         allocate(dtrad(klon,klev),topsw(klon),toplw(klon))
         allocate(solsw(klon),sollw(klon),sollwdown(klon))
         allocate(source(klon,nqmax))
         allocate(reservoir(klon))

         CALL suphec ! initialiser constantes et parametres phys.

         IF (if_ebil.ge.1) d_h_vcol_phy=0.
c
c appel a la lecture du physiq.def 
c
         call conf_phys(ok_mensuel,ok_journe,ok_instan,if_ebil)

c
c Initialiser les compteurs:
c
         itap        = 0
         itaprad     = 0
         itapchim    = 0
         ncount(:,:) = 0

c         
c Lecture startphy.nc :
c
c REMETTRE TOUS LES PARAMETRES POUR OROGW...  A FAIRE POUR TITAN
         CALL phyetat0 ("startphy.nc",
     .       rlatd,rlond,ftsol,ftsoil,
     .       falbe, solsw, sollw,
     .       dlw,radsol,reservoir,
c     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
     .       t_ancien, ancien_ok)

c dtime est defini dans tabcontrol.h et lu dans startphy 
c pdtphys est calcule a partir des nouvelles conditions:
c Reinitialisation du pas de temps physique quand changement
         IF (ABS(dtime-pdtphys).GT.0.001) THEN
            WRITE(lunout,*) 'Pas physique a change',dtime,
     .                        pdtphys
c           abort_message='Pas physique n est pas correct '
c           call abort_gcm(modname,abort_message,1)
c----------------
c pour initialiser convenablement le time_counter, il faut tenir compte
c du changement de dtime en changeant itau_phy (point de depart)
            itau_phy = NINT(itau_phy*dtime/pdtphys)
c----------------
            dtime=pdtphys
         ENDIF

         radpas  = NINT( RDAY/pdtphys/nbapp_rad)
         chimpas =   radpas*nbapp_rad/nbapp_chim

         CALL printflag( ok_mensuel,ok_journe, ok_instan )

c
c Initialiser les pas de temps:
c
      dtimerad = dtime*FLOAT(radpas)  ! pas de temps du rayonnement (s)
c      PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
            
      dtimechim = dtime*FLOAT(chimpas)  ! pas de temps de la chimie (s)
c      PRINT*,'dtimechim,dtime,chimpas',dtimechim,dtime,chimpas


c INITIALISATION ORBITE

         CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)

c---------
c FLOTT
c       IF (ok_orodr) THEN
c         DO i=1,klon
c         rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
c         ENDDO
c         CALL SUGWD(klon,klev,paprs,pplay)
c         DO i=1,klon
c         zuthe(i)=0.
c         zvthe(i)=0.
c         if(zstd(i).gt.10.)then
c           zuthe(i)=(1.-zgam(i))*cos(zthe(i))
c           zvthe(i)=(1.-zgam(i))*sin(zthe(i))
c         endif
c         ENDDO
c       ENDIF

      if (bilansmc.eq.1) then
C  OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
C  DU BILAN DE MOMENT ANGULAIRE.
      open(27,file='aaam_bud.out',form='formatted')
      open(28,file='fields_2d.out',form='formatted')
      write(*,*)'Ouverture de aaam_bud.out (FL Vous parle)'
      write(*,*)'Ouverture de fields_2d.out (FL Vous parle)'
      endif !bilansmc

c--------------SLEBONNOIS
C  OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
C  DES BALLONS
      if (ballons.eq.1) then
      open(30,file='ballons-lat.out',form='formatted')
      open(31,file='ballons-lon.out',form='formatted')
      open(32,file='ballons-u.out',form='formatted')
      open(33,file='ballons-v.out',form='formatted')
      open(34,file='ballons-alt.out',form='formatted')
      write(*,*)'Ouverture des ballons*.out'
      endif !ballons
c-------------

c---------
C TRACEURS
C source dans couche limite
          source = 0.0 ! pas de source, pour l'instant
C
c Si microphysique offline, pas besoin d'avoir de traceurs microphysiques
c car on lit les profils verticaux des qaer dans une look-up table pour 
c le rayonnement. 

c  calcul de nmicro
c !!!! Les traceurs microphysiques doivent etre toujours en premiers!!

      nmicro = 0
      do iq=1,nqmax
         nom = tname(iq)
c        print*,iq,"nom=",nom,"tname=",tname(iq)
         print*,iq,"nom=",nom
         if (nom(1:1).eq."q") then
	   nmicro = nmicro+1
	 endif
      enddo
      print*,"nmicro=",nmicro

c
c Verifications:
c
         IF ((nmicro.eq.0).and.(microfi.eq.1)) THEN
           abort_message="MICROPHYSIQUE ONLINE, MAIS NMICRO=0..."
           call abort_gcm(modname,abort_message,1) 
         ENDIF
         IF (microfi.lt.1.and.clouds.eq.1) THEN
          write(lunout,*)"microfi.lt.1.and.clouds.eq.1"
          abort_message = 
     &    "Impossible de faire des nuages sans microphysique..."
          call abort_gcm(modname,abort_message,1)
         ENDIF
         IF (nlon .NE. klon) THEN
            WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, 
     .                      klon
            abort_message='nlon et klon ne sont pas coherents'
            call abort_gcm(modname,abort_message,1)
         ENDIF
         IF (nlev .NE. klev) THEN
            WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
     .                       klev
            abort_message='nlev et klev ne sont pas coherents'
            call abort_gcm(modname,abort_message,1)
         ENDIF
c
         IF (dtime*FLOAT(radpas).GT.(RDAY*0.25).AND.cycle_diurne)
     $    THEN 
           WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
           WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
           abort_message='Nbre d appels au rayonnement insuffisant'
           call abort_gcm(modname,abort_message,1)
         ENDIF
c
         WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
     .                   iflag_ajs
c
         ecrit_mth = NINT(RDAY/dtime) *nday  ! tous les nday jours
         IF (ok_mensuel) THEN
         WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
     .                   ecrit_mth
         ENDIF

         ecrit_day = NINT(RDAY/dtime *1.0)  ! tous les jours
         IF (ok_journe) THEN
         WRITE(lunout,*)'La frequence de sortie journaliere est de ',
     .                   ecrit_day
         ENDIF

         ecrit_ins = NINT(RDAY/dtime*ecritphy)  ! Fraction de jour reglable
         IF (ok_instan) THEN
         WRITE(lunout,*)'La frequence de sortie instant. est de ', 
     .                   ecrit_ins
         ENDIF

c Initialisation des sorties 
c===========================

#ifdef CPP_IOIPSL

#ifdef histmth
#include "ini_histmth.h"
#endif

#ifdef histday
#include "ini_histday.h"
#endif

#ifdef histins
#include "ini_histins.h"
#endif

#endif

c
c Initialiser les valeurs de u pour calculs tendances
c (pour T, c'est fait dans phyetat0)
c
      DO k = 1, klev
      DO i = 1, klon
         u_ancien(i,k) = u(i,k)
      ENDDO
      ENDDO

      rmcbar  = 0.
      xfbar   = 0.
         
      ENDIF ! debut
c====================================================================
c======================================================================

c   Creer un reservoir de surface infini 
c
      reservoir(:) = 2.

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)
      ENDDO
      ENDDO
      DO iq = 1, nqmax
      DO  k = 1, klev
      DO  i = 1, klon
         tr_seri(i,k,iq) = qx(i,k,iq)
      ENDDO
      ENDDO
      ENDDO
C
      DO i = 1, klon
          ztsol(i) = ftsol(i) 
      ENDDO
C
      IF (if_ebil.ge.1) THEN 
        ztit='after dynamic'
        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
C     Comme les tendances de la physique sont ajoute dans la dynamique,
C     on devrait avoir que la variation d'entalpie par la dynamique
C     est egale a la variation de la physique au pas de temps precedent.
C     Donc la somme de ces 2 variations devrait etre nulle.
        call diagphy(airephy,ztit,ip_ebil
     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     e      , zero_v, zero_v, zero_v, ztsol
     e      , d_h_vcol+d_h_vcol_phy, d_qt, 0.
     s      , fs_bound, fq_bound )
      END IF 

c====================================================================
c Diagnostiquer la tendance dynamique
c
      IF (ancien_ok) THEN
         DO k = 1, klev
         DO i = 1, klon
            d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
            d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
         ENDDO
         ENDDO

! ADAPTATION GCM POUR CP(T)
         do i=1,klon
          flux_dyn(i,1) = 0.0
          do j=2,klev
            flux_dyn(i,j) = flux_dyn(i,j-1)
     . +cpdet(t_seri(i,j-1))/RG*d_t_dyn(i,j-1)*(paprs(i,j-1)-paprs(i,j))
          enddo
         enddo
         
      ELSE
         DO k = 1, klev
         DO i = 1, klon
            d_u_dyn(i,k) = 0.0
            d_t_dyn(i,k) = 0.0
         ENDDO
         ENDDO
         ancien_ok = .TRUE.
      ENDIF
c====================================================================
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   calcul du geopotentiel aux niveaux intercouches
c   ponderation des altitudes au niveau des couches en dp/p

      DO l=1,klev
         DO i=1,klon
c           zzlay(i,l)=zphi(i,l)/RG
c SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
            zzlay(i,l)=RG*RA*RA/(RG*RA-zphi(i,l))-RA
         ENDDO
      ENDDO
      DO i=1,klon
c        zzlev(i,1)=0.
c CORRECTION 13/01/2011  
c (correspond a la position de la surface en ce point vs RA)
         zzlev(i,1)=pphis(i)/RG
      ENDDO
      DO l=2,klev
         DO i=1,klon
            z1=(pplay(i,l-1)+paprs(i,l))/(pplay(i,l-1)-paprs(i,l))
            z2=(paprs(i,l)+pplay(i,l))/(paprs(i,l)-pplay(i,l))
            zzlev(i,l)=(z1*zzlay(i,l-1)+z2*zzlay(i,l))/(z1+z2)
         ENDDO
      ENDDO
      DO i=1,klon
         zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
      ENDDO

c- - - - - - - - - - - - - - - -
c DIAGNOSTIQUE GRILLE VERTICALE
c- - - - - - - - - - - - - - - -
c     print*,"DIAGNOSTIQUE GRILLE VERTICALE"
c     i=klon/2
c     print*,"Niveau  Pression  Altitude    (lev puis lay)"
c     do l=1,klev
c      print*,l,paprs(i,l),zzlev(i,l)
c      print*,l,pplay(i,l),zzlay(i,l)
c     enddo
c     print*,klev+1,paprs(i,klev+1),zzlev(i,klev+1)
c     stop

c====================================================================
c
c Verifier les temperatures
c
      CALL hgardfou(t_seri,ftsol,'debutphy')
c====================================================================
c
c Incrementer le compteur de la physique
c
      itap   = itap + 1

c====================================================================
c
c Epaisseurs couches

      DO k = 1, klev
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO




c====================================================================
c ORBITE ET ECLAIREMENT
c====================================================================


c Pour TITAN:
c  calcul de la longitude solaire
          CALL solarlong(rjourvrai+gmtime,zls)
          zlsdeg = zls*180./RPI      ! zls est en radians !!
          print*,'Ls',zlsdeg

      CALL orbite(zls,dist,pdecli) 
      IF (debut) zlsm1=zls

c dans zenang, Ls en degres ; dans mucorr, Ls en radians
      IF (cycle_diurne) THEN
        zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
        CALL zenang(zlsdeg,gmtime,zdtime,rlatd,rlond,rmu0,fract)
      ELSE
        call mucorr(klon,zls,rlatd,rmu0,fract)
      ENDIF

c====================================================================
c COUCHE LIMITE 
c====================================================================

c-------------------------------
c TEST: on ne tient pas compte des calculs de clmain mais on force
c l'equilibre radiatif du sol
      if (1.eq.0) then
              if (debut) then
                print*,"ATTENTION, CLMAIN SHUNTEE..."
              endif

      DO i = 1, klon
         sens(i) = 0.0e0 ! flux de chaleur sensible au sol
         fder(i) = 0.0e0
         dlw(i)  = 0.0e0
      ENDDO

c Incrementer la temperature du sol
c
      DO i = 1, klon
         d_ts(i)  = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
         ftsol(i) = ftsol(i) + d_ts(i)
         do j=1,nsoilmx
           ftsoil(i,j)=ftsol(i)
         enddo
      ENDDO

c-------------------------------
      else
c-------------------------------

      fder = dlw

c     print*,"radsol avant clmain=",radsol(klon/2)
c     print*,"solsw avant clmain=",solsw(klon/2)
c     print*,"sollw avant clmain=",sollw(klon/2)

c  CLMAIN

! ADAPTATION GCM POUR CP(T)
      CALL clmain(dtime,itap,
     e            t_seri,u_seri,v_seri,
     e            rmu0, 
     e            ftsol,
     $            ftsoil,
     $            paprs,pplay,ppk,radsol,falbe,
     e            solsw, sollw, sollwdown, fder,
     e            rlond, rlatd, cuphy, cvphy, 
     e            debut, lafin,
     s            d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
     s            fluxt,fluxu,fluxv,cdragh,cdragm,
     s            dsens,
     s            ycoefh,yu1,yv1) 

c     print*,"radsol apres clmain=",radsol(klon/2)
c     print*,"solsw apres clmain=",solsw(klon/2)
c     print*,"sollw apres clmain=",sollw(klon/2)

CXXX Incrementation des flux
      DO i = 1, klon
         sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
         fder(i) = dlw(i) + dsens(i) 
      ENDDO
CXXX

      DO k = 1, klev
      DO i = 1, klon
         t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
         d_t_vdf(i,k)= d_t_vdf(i,k)/dtime          ! K/s
         u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
         d_u_vdf(i,k)= d_u_vdf(i,k)/dtime          ! (m/s)/s
         v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
         d_v_vdf(i,k)= d_v_vdf(i,k)/dtime          ! (m/s)/s
      ENDDO
      ENDDO

c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime

C TRACEURS

      d_tr_vdf = 0.
      if (iflag_trac.eq.1) then
      DO iq=1, nqmax
          CALL cltrac(dtime,ycoefh,t_seri,
     s               tr_seri(1,1,iq), source,
     e               paprs, pplay, delp,
     s               d_tr_vdf(1,1,iq))

          tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
          d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime          ! /s
      ENDDO
      endif

      IF (if_ebil.ge.2) THEN 
        ztit='after clmain'
        CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
         call diagphy(airephy,ztit,ip_ebil
     e      , zero_v, zero_v, zero_v, zero_v, sens
     e      , zero_v, zero_v, zero_v, ztsol
     e      , d_h_vcol, d_qt, d_ec
     s      , fs_bound, fq_bound )
      END IF 
C
c
c Incrementer la temperature du sol
c
c     print*,'Tsol avant clmain:',ftsol(klon/2)
      DO i = 1, klon
         ftsol(i) = ftsol(i) + d_ts(i)
      ENDDO
c     print*,'DTsol apres clmain:',d_ts(klon/2)
c     print*,'Tsol apres clmain:',ftsol(klon/2)

c Calculer la derive du flux infrarouge
c
      DO i = 1, klon
            dlw(i) = - 4.0*emis*RSIGMA*ftsol(i)**3 
      ENDDO

c-------------------------------
      endif  ! fin du TEST

c
c Appeler l'ajustement sec
c
c===================================================================
c Convection seche 
c===================================================================
c
      d_t_ajs(:,:)=0.
      d_u_ajs(:,:)=0.
      d_v_ajs(:,:)=0.
      d_tr_ajs(:,:,:)=0.
c
      IF(prt_level>9)WRITE(lunout,*)
     .    'AVANT LA CONVECTION SECHE , iflag_ajs='
     s   ,iflag_ajs

      if(iflag_ajs.eq.0) then
c  Rien
c  ====
         IF(prt_level>9)WRITE(lunout,*)'pas de convection'

      else if(iflag_ajs.eq.1) then

c  Ajustement sec
c  ==============
         IF(prt_level>9)WRITE(lunout,*)'ajsec'

! ADAPTATION GCM POUR CP(T)
         CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
     .              tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)

! ADAPTATION GCM POUR CP(T)
         do i=1,klon
          flux_ajs(i,1) = 0.0
          do j=2,klev
            flux_ajs(i,j) = flux_ajs(i,j-1)
     .        + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
     .                                 *delp(i,j-1)
          enddo
         enddo
         
         t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
         d_t_ajs(:,:)= d_t_ajs(:,:)/dtime          ! K/s
         u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
         d_u_ajs(:,:)= d_u_ajs(:,:)/dtime          ! (m/s)/s
         v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
         d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
      if (iflag_trac.eq.1) then
         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
         d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime          ! /s
      endif
      
c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime

      endif
c
      IF (if_ebil.ge.2) THEN 
        ztit='after dry_adjust'
        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
        call diagphy(airephy,ztit,ip_ebil
     e      , zero_v, zero_v, zero_v, zero_v, sens
     e      , zero_v, zero_v, zero_v, ztsol
     e      , d_h_vcol, d_qt, d_ec
     s      , fs_bound, fq_bound )
      END IF 

c====================================================================
c   MICROPHYSIQUE ET CHIMIE
c====================================================================

      d_tr_mph(:,:,:)=0.
      d_tr_kim(:,:,:)=0.
      
c on recupere tr_seri inchange, d_tr_micro, d_tr_chim, tous les trois sur nqmax
c on recupere aussi qaer pour le mettre dans les sorties
c  si microfi=1, sortie de qaer(1:nmicro)
c  si nmicro != nqmax et si chimi, sortie de tr_seri(nmicro+1:nqmax) 

c faire un test comme pour rayonnement, avec chimi en plus comme flag, 
c pour voir si chimie appelee -> bouleen, qui passe dans phytrac.
c faut aussi le pas de temps chimique: dtimechim, a passer..

      appel_chim = 0
      IF (MOD(itapchim,chimpas).EQ.0) THEN
c             print*,'CHIMIE ', 
c    $             ' (itapchim=',itapchim,'/chimpas=',chimpas,')'
       appel_chim = 1
       itapchim = 0
      ENDIF
      itapchim = itapchim + 1

      if (iflag_trac.eq.1) then
c         call begintime(tt0)
         call phytrac (debut,lafin,
     .                   nqmax,nmicro,dtime,appel_chim,dtimechim,
     .                   paprs,pplay,delp,t,rmu0,fract,pdecli,zls,
     .                   yu1,yv1,zzlev,zzlay,ftsol,
     .                   tr_seri,qaer,d_tr_mph,d_tr_kim,
     .                   fclat,reservoir)

c         call endtime(tt0,tt1)
c         ttphytra=ttphytra+tt1

c ----- ICI on ajuste radsol en tenant compte du flux de chaleur latente 
c       d'evaporation du reservoir.
c       NOTE : c'est pas tres elegant mais ca permet d'eviter d'aller
c              toucher a clmain.
        if (clouds.eq.1) then
          radsol(:) = radsol(:)+fclat(:)    !test pas de flx de chaleur latente
        endif

        if (microfi.ge.1) then
         tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
     .                        + d_tr_mph(:,:,1:nmicro)*dtime
	endif
c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
c       retourne des traceurs nuls et il y a parfois des valeurs negatives qui trainent.
c       Pour ne diffuser le probleme, on force les valeurs negatives a ZERO.
        DO iq=1,nmicro
          DO i=1,klon
            DO l=1,klev
              if (tr_seri(i,l,iq).lt.0.) then
                 tr_seri(i,l,iq) = 0.
              endif
            ENDDO
          ENDDO
        ENDDO

c condensation:
c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
        if ((clouds.eq.1.or.(chimi)).and.nqmax.gt.nmicro) then
          tr_seri(:,:,nmicro+1:nqmax) = tr_seri(:,:,nmicro+1:nqmax)
     .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
        endif
        if ((chimi).and.(nqmax.gt.nmicro)) then
c chimie:
         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
	endif

      endif

c       ch4=0.
c       do l=1,llm
c         ch4(1,l) = tr_seri(1,l,ich4)
c         do j=2,jjm
c           ig0=1+(j-2)*iim
c           do i=1,iim
c             ch4(j,l)= ch4(j,l)  + tr_seri(ig0+i,l,ich4)/iim
c           enddo
c         enddo
c         ch4(jjm+1,l) = tr_seri(klon,l,ich4)
c       enddo
c       do j=1,jjm+1
c         write(501,*) j,ch4(j,1)
c       enddo
c       do l=1,llm
c         write(502,'(I3,49(ES24.17,1X))') l, 
c     &   (ch4(j,l),j=1,jjm+1)
c       enddo
c       write(501,*) ""
c       write(502,*) ""

c------------------
c test condensation
c     do i=1,nqmax 
c       if(tname(i).eq."HCN") then
c          print*,"HCN="
c          do k=1,klev
c           print*,k,tr_seri(klon/2,k,i),d_tr_mph(klon/2,k,i)*dtime
c    v      ,d_tr_kim(klon/2,k,i)*dtime
c          enddo
c          stop
c       endif
c     enddo
c------------------

c====================================================================
c RAYONNEMENT
c====================================================================

      IF (MOD(itaprad,radpas).EQ.0) THEN
c             print*,'RAYONNEMENT ', 
c    $             ' (itaprad=',itaprad,'/radpas=',radpas,')'

c ATTENTION, (klon/2) ne marche pas toujours............
c     print*,"radsol avant radlwsw=",radsol(klon/2)
c     print*,"solsw avant radlwsw=",solsw(klon/2)
c     print*,"sollw avant radlwsw=",sollw(klon/2)
c     print*,"avant radlwsw"

c   ----------------
c   Calcul du rayon moyen des gouttes et des fractions volumique pour le TR
c  ----------------
      IF (clouds.eq.1) THEN
        DO i=1,klon
          DO j=1,klev
            rmcbar(i,j)=rmcbar(i,j)/MAX(FLOAT(ncount(i,j)),1.)
            xfbar(i,j,:)=xfbar(i,j,:)/MAX(FLOAT(ncount(i,j)),1.)
          ENDDO 
        ENDDO
      ENDIF
     
c      call begintime(tt0)
      CALL radlwsw
     e            (dist, rmu0, fract, falbe, zzlev,
     e             paprs, pplay,ftsol, t_seri, nqmax, nmicro,
     c             tr_seri, qaer,
     s             heat,cool,radsol,
     s             topsw,toplw,solsw,sollw,
     s             sollwdown,
     s             lwnet, swnet)
c      call endtime(tt0,tt1)
c      ttrad=ttrad+tt1

c     print*,"apres radlwsw"
c     mise a zero du rayon moyen des gouttes et des fractions volumique 
      IF (clouds.eq.1) THEN
        rmcbar(:,:)  = 0.
        xfbar(:,:,:) = 0.
        ncount(:,:)  = 0
      ENDIF

c     print*,"radsol apres radlwsw=",radsol(klon/2)
c     print*,"solsw apres radlwsw=",solsw(klon/2)
c     print*,"sollw apres radlwsw=",sollw(klon/2)
      itaprad = 0
      DO k = 1, klev
       DO i = 1, klon
         dtrad(i,k) = heat(i,k)-cool(i,k)     !K/s
       ENDDO
      ENDDO
c       print*,"heat (K/s) =",heat(klon/2,:)
c       print*,"cool (K/s) =",cool(klon/2,:)
c       print*,"dtrad1 (K/s) =",dtrad(1,:)
c       print*,"dtrad2 (K/s) =",dtrad(klon/2,:)
c       print*,"dtrad3 (K/s) =",dtrad(klon,:)
	    
      ENDIF
      itaprad = itaprad + 1
c====================================================================
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) + dtrad(i,k) * dtime
      ENDDO
      ENDDO
 
      IF (if_ebil.ge.2) THEN 
        ztit='after rad'
        CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
        call diagphy(airephy,ztit,ip_ebil
     e      , topsw, toplw, solsw, sollw, zero_v
     e      , zero_v, zero_v, zero_v, ztsol
     e      , d_h_vcol, d_qt, d_ec
     s      , fs_bound, fq_bound )
      END IF 
c

c====================================================================
c+jld ec_conser
      DO k = 1, klev
      DO i = 1, klon
        d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
     $      *(u(i,k)**2+v(i,k)**2-u_seri(i,k)**2-v_seri(i,k)**2)
        t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
        d_t_ec(i,k) = d_t_ec(i,k)/dtime
       END DO 
      END DO 
         do i=1,klon
          flux_ec(i,1) = 0.0
          do j=2,klev
            flux_ec(i,j) = flux_ec(i,j-1)
     . +cpdet(t_seri(i,j-1))/RG*d_t_ec(i,j-1)*delp(i,j-1)
          enddo
         enddo
         
c-jld ec_conser
c====================================================================

      IF (if_ebil.ge.1) THEN 
        ztit='after physic'
        CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
     e      , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
C     Comme les tendances de la physique sont ajoute dans la dynamique,
C     on devrait avoir que la variation d'entalpie par la dynamique
C     est egale a la variation de la physique au pas de temps precedent.
C     Donc la somme de ces 2 variations devrait etre nulle.
        call diagphy(airephy,ztit,ip_ebil
     e      , topsw, toplw, solsw, sollw, sens
     e      , zero_v, zero_v, zero_v, ztsol
     e      , d_h_vcol, d_qt, d_ec
     s      , fs_bound, fq_bound )
C
      d_h_vcol_phy=d_h_vcol
C
      END IF 
C
c====================================================================
c   Calcul  des gravity waves  FLOTT
c====================================================================
c
c      if (ok_orodr.or.ok_gw_nonoro) then
cc  CALCUL DE N2
c       do i=1,klon
c        do k=2,klev
c	  ztlev(i,k)  = (t_seri(i,k)+t_seri(i,k-1))/2.
c	  zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
c	enddo
c       enddo
c       call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
c       call t2tpot(klon*klev,t_seri,ztetalay,ppk)
c       do i=1,klon
c        do k=2,klev
c	  zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
c	  zdzlev(i,k)    = (zphi(i,k)-zphi(i,k-1))/RG
c          zn2(i,k) = RG*zdtetalev(i,k)/(ztetalev(i,k)*zdzlev(i,k))
c          zn2(i,k) = max(zn2(i,k),1.e-12)  ! securite
c	enddo
c       enddo
c
c      endif
c      
cc ----------------------------ORODRAG
c      IF (ok_orodr) THEN
cc
cc  selection des points pour lesquels le shema est actif:
c        igwd=0
c        DO i=1,klon
c        itest(i)=0
cc        IF ((zstd(i).gt.10.0)) THEN
c        IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
c          itest(i)=1
c          igwd=igwd+1
c          idx(igwd)=i
c        ENDIF
c        ENDDO
cc        igwdim=MAX(1,igwd)
cc
cc A ADAPTER POUR TITAN !!!
c        CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
c     e                   zmea,zstd, zsig, zgam, zthe,zpic,zval,
c     e                   igwd,idx,itest,
c     e                   t_seri, u_seri, v_seri,
c     s                   zulow, zvlow, zustrdr, zvstrdr,
c     s                   d_t_oro, d_u_oro, d_v_oro)
c
cc  ajout des tendances
c           t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
c           d_t_oro(:,:)= d_t_oro(:,:)/dtime          ! K/s
c           u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
c           d_u_oro(:,:)= d_u_oro(:,:)/dtime          ! (m/s)/s
c           v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
c           d_v_oro(:,:)= d_v_oro(:,:)/dtime          ! (m/s)/s
cc    
c      ELSE
c         d_t_oro = 0.
c         d_u_oro = 0.
c         d_v_oro = 0.
c	 zustrdr = 0.
c	 zvstrdr = 0.
cc
c      ENDIF ! fin de test sur ok_orodr
cc
cc ----------------------------OROLIFT
c      IF (ok_orolf) THEN
cc
cc  selection des points pour lesquels le shema est actif:
c        igwd=0
c        DO i=1,klon
c        itest(i)=0
c        IF ((zpic(i)-zmea(i)).GT.100.) THEN
c          itest(i)=1
c          igwd=igwd+1
c          idx(igwd)=i
c        ENDIF
c        ENDDO
cc        igwdim=MAX(1,igwd)
cc
cc A ADAPTER POUR VENUS!!!
cc            CALL lift_noro(klon,klev,dtime,paprs,pplay,
cc     e                   rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
cc     e                   igwd,idx,itest,
cc     e                   t_seri, u_seri, v_seri,
cc     s                   zulow, zvlow, zustrli, zvstrli,
cc     s                   d_t_lif, d_u_lif, d_v_lif               )
c
cc
cc  ajout des tendances
c           t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
c           d_t_lif(:,:)= d_t_lif(:,:)/dtime          ! K/s
c           u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
c           d_u_lif(:,:)= d_u_lif(:,:)/dtime          ! (m/s)/s
c           v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
c           d_v_lif(:,:)= d_v_lif(:,:)/dtime          ! (m/s)/s
cc
c      ELSE
c         d_t_lif = 0.
c         d_u_lif = 0.
c         d_v_lif = 0.
c         zustrli = 0.
c         zvstrli = 0.
cc
c      ENDIF ! fin de test sur ok_orolf
c
cc ---------------------------- NON-ORO GRAVITY WAVES
c       IF(ok_gw_nonoro) then
c
c      call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
c     e               t_seri, u_seri, v_seri,
c     o               zustrhi,zvstrhi,
c     o               d_t_hin, d_u_hin, d_v_hin)
c
cc  ajout des tendances
c
c         t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
c         d_t_hin(:,:)= d_t_hin(:,:)/dtime          ! K/s
c         u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
c         d_u_hin(:,:)= d_u_hin(:,:)/dtime          ! (m/s)/s
c         v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
c         d_v_hin(:,:)= d_v_hin(:,:)/dtime          ! (m/s)/s
c
c      ELSE
c         d_t_hin = 0.
c         d_u_hin = 0.
c         d_v_hin = 0.
c         zustrhi = 0.
c         zvstrhi = 0.
c
c      ENDIF ! fin de test sur ok_gw_nonoro
c
c====================================================================
c Transport de ballons 
c====================================================================
      if (ballons.eq.1) then
         CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
c    C               t,pplay,u,v,pphi)   ! alt above surface (smoothed for GCM)
     C               t,pplay,u,v,zphi)   ! alt above planet average radius
      endif !ballons

c====================================================================
c Bilan de mmt angulaire
c====================================================================
      if (bilansmc.eq.1) then
CMODDEB FLOTT
C  CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE

      DO i = 1, klon
        zustrph(i)=0.
        zvstrph(i)=0.
        zustrcl(i)=0.
        zvstrcl(i)=0.
      ENDDO
      DO k = 1, klev
      DO i = 1, klon
       zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
     c            (paprs(i,k)-paprs(i,k+1))/rg
       zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
     c            (paprs(i,k)-paprs(i,k+1))/rg
       zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
     c            (paprs(i,k)-paprs(i,k+1))/rg
       zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
     c            (paprs(i,k)-paprs(i,k+1))/rg
      ENDDO
      ENDDO

      CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
     C               ra,rg,romega,
     C               rlatd,rlond,pphis,
     C               zustrdr,zustrli,zustrcl,
     C               zvstrdr,zvstrli,zvstrcl,
     C               paprs,u,v)
                     
CCMODFIN FLOTT
      endif !bilansmc

c=======================================================================
c   SORTIES
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
      ENDDO
      ENDDO
c     print*,"vnatphy=",v(705,:)
c     print*,"unatphy=",u(705,:)
c
      DO iq = 1, nqmax
      DO  k = 1, klev
      DO  i = 1, klon
         d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
      ENDDO
      ENDDO
      ENDDO
      
c------------------------
c Calcul moment cinetique
c------------------------
c TEST 
c     mangtot = 0.0
c     DO k = 1, klev
c     DO i = 1, klon
c       mang(i,k) = RA*cos(rlatd(i)*RPI/180.)
c    .     *(u_seri(i,k)+RA*cos(rlatd(i)*RPI/180.)*ROMEGA)
c    .     *airephy(i)*delp(i,k)/RG
c       mangtot=mangtot+mang(i,k)
c     ENDDO
c     ENDDO
c     print*,"Moment cinetique total = ",mangtot
c
c------------------------
c
c Sauvegarder les valeurs de t et u a la fin de la physique:
c
      DO k = 1, klev
      DO i = 1, klon
         u_ancien(i,k) = u_seri(i,k)
         t_ancien(i,k) = t_seri(i,k)
      ENDDO
      ENDDO
c
c=============================================================
c   Ecriture des sorties
c=============================================================
      
#ifdef CPP_IOIPSL

#ifdef histmth
#include "write_histmth.h"
#endif

#ifdef histday
#include "write_histday.h"
#endif

#ifdef histins
#include "write_histins.h"
#endif
              
#endif

c====================================================================
c Si c'est la fin, il faut conserver l'etat de redemarrage
c====================================================================
c
      IF (lafin) THEN
         itau_phy = itau_phy + itap
         lsinit   = zlsdeg
c REMETTRE TOUS LES PARAMETRES POUR OROGW... A FAIRE POUR TITAN
         CALL phyredem ("restartphy.nc",
     .      rlatd, rlond, ftsol, ftsoil,
     .      falbe,
     .      solsw, sollw,dlw,
     .      radsol,reservoir,
c     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,
     .      t_ancien)
      
c--------------FLOTT
CMODEB LOTT
C  FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
C  DU BILAN DE MOMENT ANGULAIRE.
      if (bilansmc.eq.1) then
        write(*,*)'Fermeture de aaam_bud.out (FL Vous parle)'
        close(27)                                     
        close(28)                                     
      endif !bilansmc
CMODFIN
c-------------
c--------------SLEBONNOIS
C  FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
C  DES BALLONS
      if (ballons.eq.1) then
        write(*,*)'Fermeture des ballons*.out'
        close(30)                                     
        close(31)                                     
        close(32)                                     
        close(33)                                     
        close(34)                                     
      endif !ballons
c-------------
      ENDIF
      

      RETURN
      END



***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************

      SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
      IMPLICIT none
c
c Tranformer une variable de la grille physique a
c la grille d'ecriture
c
      INTEGER nfield,nlon,iim,jjmp1, jjm
      REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
c
      INTEGER i, n, ig
c
      jjm = jjmp1 - 1
      DO n = 1, nfield
         DO i=1,iim
            ecrit(i,n) = fi(1,n)
            ecrit(i+jjm*iim,n) = fi(nlon,n)
         ENDDO
         DO ig = 1, nlon - 2
           ecrit(iim+ig,n) = fi(1+ig,n)
         ENDDO
      ENDDO
      RETURN
      END
