PROGRAM testphys1d_triton

! to use  'getin'
      USE ioipsl_getincom 
      use planet_h          
      use comgeomfi_h
      use comsoil_h
      IMPLICIT NONE

!==================================================================
!     
!     Purpose
!     -------
!     Run the physics package of the universal model in a 1D column.
!     
!     It can be compiled with a command like (e.g. for 25 layers):
!                                  "makegcm -p pluto -d 25 testphys1d"
!     It requires the files "callphys.def", "z2sig.def",
!     "traceur.def", and "run1d.def" with a line "INCLUDEDEF=callphys.def"
!
!     Authors
!     -------
!     Frederic Hourdin
!     R. Fournier
!     F. Forget
!     F. Montmessin (water ice added)

!==================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "surfdat.h"
!#include "comsoil.h"
#include "comdiurn.h"
#include "callkeys.h"
#include "comcstfi.h"
#include "comsaison.h"
#include "control.h"
#include "comvert.h"
#include "netcdf.inc"
#include "comg1d.h"
#include "fisice.h"
#include "logic.h"
#include "advtrac.h"

c --------------------------------------------------------------
c  Declarations
c --------------------------------------------------------------
c
      INTEGER nlayer,nlevel,nsoil,ndt
      INTEGER ilayer,ilevel,isoil,idt,iq
      LOGICAl firstcall,lastcall
c
      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
      INTEGER lecttsoil     ! lecture of tsoil from proftsoil
      INTEGER lecthaze      ! lecture of haze from profhaze
      REAL day           ! date durant le run
      REAL AU           ! astronomical unit AU=149 million km
      REAL time             ! time (0<time<1 ; time=0.5 a midi)
      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
      REAL psurf,tsurf      
      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
      REAL gru,grv   ! prescribed "geostrophic" background wind
      REAL temp(nlayermx)   ! temperature at the middle of the layers
      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
      integer :: i_n2=0  ! tracer index of n2 ice
      integer :: i_ch4_ice=0  ! tracer index of Ch4 ice
      integer :: i_ch4_gas=0  ! tracer index of Ch4 gas
      integer :: i_co_ice=0  ! tracer index of Co ice
      integer :: i_co_gas=0  ! tracer index of Co gas
      integer :: i_prec_haze=0  ! tracer index of haze precursor
      integer :: i_haze=0  ! tracer index of haze
      integer :: i_haze10=0  ! tracer index of haze
      integer :: i_haze30=0  ! tracer index of haze
      integer :: i_haze50=0  ! tracer index of haze
      integer :: i_haze100=0  ! tracer index of haze
      REAL emis             ! surface layer
      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
      REAL  Lsperi,excentric

c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
      REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
      REAL dpsurf    
      REAL dq(nlayermx,nqmx)

c   Various intermediate variables
      REAL zls
      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
      REAL pks, ptif, w(nlayermx)
      INTEGER ierr
      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
      Logical  tracerdyn
      integer :: nq=1 ! number of tracers

      character (len=7) :: str7

      logical saveprofile

!     added by RW for zlay computation
      real Hscale, rho, dz

      real ch4ref ! % ch4
      real coref ! % co
      real hazeref ! haze kg/kg
      real prechazeref ! prec haze kg/kg


c=======================================================================

c=======================================================================
c INITIALISATION
c=======================================================================

      saveprofile=.true.

c ------------------------------------------------------
c  Default values for constants (corresponding to Triton)
c ------------------------------------------------------

      pi=2.E+0*asin(1.E+0)

c     Constante de Triton
c     ----------------------------
      rad=1353400.               ! rayon de Triton (m)
      daysec=507773            ! duree du sol (s) 
      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)
      g=0.78                 ! gravite (m.s-2) 
      mugaz=28.                ! Masse molaire de l'atm (g.mol-1) ~28
      rcp=.2853                ! = r/cp  ~0.2853
      r= 8.314511E+0 *1000.E+0/mugaz
      cpp= r/rcp
      year_day = 10247         ! duree de l'annee (sols) 
      periheli = 4444.45       ! dist.min. soleil-Triton (Mkm) 
      aphelie = 4545.67         ! dist.max. soleil-Triton (Mkm) 
      peri_day =  0      ! date du perihelie (sols depuis printemps)
      obliquit=28.32         ! Obliquite de la planete (deg)         
      excentric= 0.      ! Excentricite 
 
c     Parametres Couche limite et Turbulence 
c     --------------------------------------
      z0 =  1.e-2                ! surface roughness (m) ~0.01 
 
c     ----------------------------------------------------
      hybrid=.true.

!     Constants for stokes.F90
      molrad=1.93e-10 ! N2 TB16!!!
      visc=6.67e-6     ! N2 TB16!!!

c ------------------------------------------------------
c  Load parameters from "run.def" 
c ------------------------------------------------------
      open(90,file='run.def',status='old',form='formatted',
     &     iostat=ierr)
      if (ierr.ne.0) then
        write(*,*) 'Cannot find required file "run.def"'
        write(*,*) '  (which should contain some input parameters'
        write(*,*) '   along with the following line:'
        write(*,*) '   INCLUDEDEF=callphys.def'
        write(*,*) '   )'
        write(*,*) ' ... might as well stop here ...'
        stop
      else
        close(90)
      endif

! check if we are going to run with or without tracers
      write(*,*) "Run with or without tracer transport ?"
      tracer=.false. ! default value
      call getin("tracer",tracer)
      write(*,*) " tracer = ",tracer

! while we're at it, check if there is a 'traceur.def' file
! and preocess it, if necessary. Otherwise initialize tracer names
      if (tracer) then
      ! load tracer names from file 'traceur.def'
        open(90,file='traceur.def',status='old',form='formatted',
     &       iostat=ierr)
        if (ierr.eq.0) then
          write(*,*) "testphys1d: Reading file traceur.def"
          ! read number of tracers:
          read(90,*,iostat=ierr) nq
          if (ierr.ne.0) then
            write(*,*) "testphys1d: error reading number of tracers"
            write(*,*) "   (first line of traceur.def) "
            stop
          else
            ! check that the number of tracers is indeed nqmx
            if (nq.ne.nqmx) then
              write(*,*) "testphys1d: error, wrong number of tracers:"
              write(*,*) "nq=",nq," whereas nqmx=",nqmx
              stop
            endif
          endif
        
          ! initialize advection schemes to Van-Leer for all tracers
          do iq=1,nq
            iadv(iq)=3 ! Van-Leer
          enddo
        
          do iq=1,nq
          ! minimal version, just read in the tracer names, 1 per line
            read(90,*,iostat=ierr) tnom(iq)
            write(*,*) 'Ini 1d reading traceur.def: ', tnom(iq),iq
            if (ierr.ne.0) then
              write(*,*) 'testphys1d: error reading tracer names...'
              stop
            endif
          enddo !of do iq=1,nq

! TB check for n2_ice / ch4 tracers:
         write(*,*) 'Tracers: tnom=', tnom(:)
         do iq=1,nq
           if (tnom(iq)=="n2") then
             i_n2=iq
           elseif (tnom(iq)=="ch4_ice") then
             i_ch4_ice=iq
           elseif (tnom(iq)=="ch4_gas") then
             i_ch4_gas=iq
           elseif (tnom(iq)=="prec_haze") then
             i_prec_haze=iq
           elseif (tnom(iq)=="haze") then
             i_haze=iq
           elseif (tnom(iq)=="haze10") then
             i_haze10=iq
           elseif (tnom(iq)=="haze30") then
             i_haze30=iq
           elseif (tnom(iq)=="haze50") then
             i_haze50=iq
           elseif (tnom(iq)=="haze100") then
             i_haze100=iq
           elseif (tnom(iq)=="co_ice") then
             i_co_ice=iq
           elseif (tnom(iq)=="co_gas") then
             i_co_gas=iq
           endif
         enddo

        else  ! ierr=0
          write(*,*) 'Cannot find required file "traceur.def"'
          write(*,*) ' If you want to run with tracers, I need it'
          write(*,*) ' ... might as well stop here ...'
          stop
        endif
        close(90)

      else
      ! we still need to set (dummy) tracer names for physdem1
        nq=nqmx
        do iq=1,nq
          write(str7,'(a1,i2.2)')'q',iq
          tnom(iq)=str7
        enddo
        ! actually, we'll need at least one "n2_ice" tracer
        ! (for surface N2 ice)
        write(*,*) "No tracer ! tracer=false"
        tnom(1)="n2"
        i_n2=1
      endif ! of if (tracer)

c  Date et heure locale du debut du run
c  ------------------------------------
c  Date (en sols depuis le solstice de printemps) du debut du run
      day0 = 0 ! default value for day0
      write(*,*) 'Initial date (in sols ; =0 at Ls=0)?'
      call getin("day0",day0)
      day=float(day0)
      write(*,*) " day0 = ",day0
c  Heure de demarrage
      time=0 ! default value for time
      write(*,*)'Initial local time (in hours, between 0 and 24)?'
      call getin("time",time)
      write(*,*)" time = ",time
      time=time/24.E+0 ! convert time (hours) to fraction of sol

c  Discretisation (Definition de la grille et des pas de temps)
c  --------------
c
      nlayer=nlayermx
      nlevel=nlayer+1
      nsoil=nsoilmx
      PRINT *,'Dims nlevel=',nlevel,' nsoil=',nsoil

      PRINT *,'Length of day (s) ?'
      call getin("daysec",daysec)
      write(*,*) " daysec = ",daysec
      omeg=4.*asin(1.)/(daysec)  ! vitesse de rotation (rad.s-1)

      PRINT *,'Length of year (sol) ?'
      call getin("year_day",year_day)
      write(*,*) " year_day= ",year_day

      day_step=48 ! default value for day_step
      PRINT *,'Number of time steps per sol ?'
      call getin("day_step",day_step)
      write(*,*) " day_step = ",day_step

      ndt=10 ! default value for ndt
      PRINT *,'Number of sols to run ?'
      call getin("ndt",ndt)
      write(*,*) " ndt = ",ndt

      ndt=ndt*day_step     
      dtphys=daysec/day_step  

c Pression de surface sur la planete
c ------------------------------------
      psurf=1. ! default value for psurf
      PRINT *,'Surface pressure (Pa) ?'
      call getin("psurf",psurf)
      write(*,*) " psurf = ",psurf
c Pression de reference
      preff=1. ! these values are not needed in 1D  
      pa=0.25*preff
c Gravity of planet
c -----------------
      g=0.78 ! default value for g
      PRINT *,'Gravity ?'
      call getin("g",g)
      write(*,*) " g = ",g

c Molar mass of atmosphere
c ------------------------
      PRINT *,'Molar mass of atmosphere ?'
      call getin("mugaz",mugaz)
      write(*,*) " mugaz = ",mugaz

c Specific heat capacity of atmosphere
c ------------------------------------
      PRINT *,'Specific heat capacity of atmosphere ?'
      call getin("cpp",cpp)
      write(*,*) " cpp = ",cpp

      r= 8.314511E+0 *1000.E+0/mugaz
      rcp=r/cpp
 
c  latitude/longitude
c  ------------------
      lati(1)=0 ! default value for lati(1)
      PRINT *,'latitude (in degrees) ?'
      call getin("latitude",lati(1))
      write(*,*) " latitude = ",lati(1)
      lati(1)=lati(1)*pi/180.E+0
      long(1)=0.E+0
      long(1)=long(1)*pi/180.E+0

c  periastron/apoastron
c  --------------------
      PRINT *,'periastron (in AU) ?'
      call getin("periheli",periheli)
      write(*,*) "periastron = ",periheli
      AU=149.597870700 ! km
      periheli=periheli*AU ! AU to Mkm
      PRINT *,'apoastron (in AU) ?'
      call getin("aphelie",aphelie)
      write(*,*) "apoastron = ",aphelie
      aphelie=aphelie*AU ! AU to Mkm

      excentric=(1-periheli*periheli/(aphelie*aphelie) )**0.5

c  obliquity
c  ---------
      PRINT *,'obliquity (in deg) ?'
      call getin("obliquit",obliquit)
      write(*,*) "obliquity = ",obliquit

c  Initialisation albedo /emis / inertie du sol
c  --------------------------------------
      albedodat(1)=0.2 ! default value for albedodat
      PRINT *,'Albedo of bare ground ?'
      call getin("albedo",albedodat(1))
      write(*,*) " albedo = ",albedodat(1)

      emis=0.8 ! default value for emissivity
      PRINT *,'Emissivity of bare ground ?'
      call getin("emis",emis)
      write(*,*) " emis = ",emis

      inertiedat(1,1)=400 ! default value for inertiedat
      PRINT *,'Soil thermal inertia (SI) ?'
      call getin("inertia",inertiedat(1,1))
      write(*,*) " inertia = ",inertiedat(1,1)
c
c  pour le schema d'ondes de gravite
c  ---------------------------------
      zmea(1)=0.E+0
      zstd(1)=0.E+0
      zsig(1)=0.E+0
      zgam(1)=0.E+0
      zthe(1)=0.E+0

c    Initialisation des traceurs
c    ---------------------------
      DO iq=1,nqmx
        if (iq.eq.i_n2) then
           DO ilayer=1,nlayer
              q(ilayer,iq) = 1 
           ENDDO
           write(*,*) 'ini 1D traceur ',iq,' : q_n2 = ', q(:,iq)
        else if (iq.eq.i_ch4_gas) then
          ch4ref=0.5 ! default value for vmr ch4 
          PRINT *,'vmr CH4 (%) ?'
          call getin("ch4ref",ch4ref)
          write(*,*) " CH4 ref (%) = ",ch4ref
          DO ilayer=1,nlayer
             q(ilayer,iq) = 0.01*ch4ref*16./28.
          ENDDO
          write(*,*) 'ini 1D traceur ',iq,' : q_ch4_gas = ', q(:,iq)
        else if (iq.eq.i_co_gas) then
          coref=0.05 ! default value for vmr ch4 
          PRINT *,'vmr CO (%) ?'
          call getin("coref",coref)
          write(*,*) " CO ref (%) = ",coref
          DO ilayer=1,nlayer
             q(ilayer,iq) = 0.01*coref*28./28.
          ENDDO
          write(*,*) 'ini 1D traceur ',iq,' : q_co_gas = ', q(:,iq)
        else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30).or.
     &     (iq.eq.i_haze).or.(iq.eq.i_haze50).or.(iq.eq.i_haze100)) then
          hazeref=0. ! default value for haze mmr
          PRINT *,'qhaze (kg/kg) ?'
          call getin("hazeref",hazeref)
          write(*,*) " haze ref (kg/kg) = ",hazeref
          DO ilayer=1,nlayer
             q(ilayer,iq) = hazeref
          ENDDO
          write(*,*) 'ini 1D traceur ',iq,' : q_haze = ', q(:,iq)
        else if (iq.eq.i_prec_haze) then
          prechazeref=0. ! default value for vmr ch4 
          PRINT *,'qprechaze (kg/kg) ?'
          call getin("prechazeref",prechazeref)
          write(*,*) " prec haze ref (kg/kg) = ",prechazeref
          DO ilayer=1,nlayer
             q(ilayer,iq) = prechazeref
          ENDDO
          write(*,*) 'ini 1D traceur ',iq,' : q_prec_haze = ', q(:,iq)
          
        else
           DO ilayer=1,nlayer
              q(ilayer,iq) = 0.
           ENDDO
        endif
      ENDDO

c     TB: clean surface
      DO iq=1,nqmx
        qsurf(iq) = 0.
      ENDDO

c   Initialisation speciales "physiq"
c   ---------------------------------
c   la surface de chaque maille est inutile en 1D --->
      area(1)=1.E+0

c   le geopotentiel au sol est inutile en 1D car tout est controle
c   par la pression de surface --->
      phisfi(1)=0.E+0

c   "inifis" reproduit un certain nombre d'initialisations deja faites
c  + lecture des clefs de callphys.def
c  + calcul de la frequence d'appel au rayonnement
c  + calcul des sinus et cosinus des longitude latitude

!   Possible issue with dtphys in input and include!!!
      CALL inifis(1,llm,day0,daysec,dtphys,
     .            lati,long,area,rad,g,r,cpp)
c   Initialisation pour prendre en compte les vents en 1-D
c   ------------------------------------------------------
      ptif=2.E+0*omeg*sinlat(1)

c    vent geostrophique
      gru=10. ! default value for gru
      PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
      call getin("u",gru)
      write(*,*) " u = ",gru
      grv=0. !default value for grv
      PRINT *,'meridional northward component of the geostrophic',
     &' wind (m/s) ?'
      call getin("v",grv)
      write(*,*) " v = ",grv

c     Initialisation des vents  au premier pas de temps
      DO ilayer=1,nlayer
         u(ilayer)=gru
         v(ilayer)=grv
      ENDDO

c     energie cinetique turbulente
      DO ilevel=1,nlevel
         q2(ilevel)=0.E+0
      ENDDO

c     Surface 
c     -------------------------------------------
      if(i_n2.gt.0)then
         qsurf(i_n2)=0 ! default value for N2ice
         print*,'Initial N2 ice on the surface (kg.m-2)'
         call getin("n2ice",qsurf(i_n2))
         write(*,*) " n2ice = ",qsurf(i_n2)
      endif

c  calcul des pressions et altitudes en utilisant les niveaux sigma
c  ----------------------------------------------------------------

c    Vertical Coordinates
c    """"""""""""""""""""
      hybrid=.false.
      PRINT *,'Hybrid coordinates ?'
      call getin("hybrid",hybrid)
      write(*,*) " hybrid = ", hybrid

      CALL disvert

      ! we want only the scale height from z2sig, in order to compute zlay
      open(91,file="z2sig.def",form='formatted')
      read(91,*) Hscale
      close(91)

      DO ilevel=1,nlevel
        plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
      ENDDO

      DO ilayer=1,nlayer
        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
      ENDDO

      DO ilayer=1,nlayer
!        zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
!     &   /g
        zlay(ilayer)=-1000.0*Hscale*log(play(ilayer)/plev(1))
      ENDDO

c   Orbital parameters
c   ------------------
      PRINT *,'Solar longitude of perihelion Lsperi '
      call getin("Lsperi",Lsperi)
      write(*,*) " Lsperi = ", Lsperi
      
      Lsperi = Lsperi * pi/180.   ! Ls of perihelion

c     Calcul du sol correspondant a Lsperi
      call call_dayperi(Lsperi,excentric,peri_day,year_day)
      PRINT *,'date of perihelion (sol)',peri_day

c   Profil de temperature au premier appel
c   --------------------------------------
      pks=psurf**rcp

c   Altitude en km dans profile: on divise zlay par 1000
      tmp1(0)=0.E+0
      DO ilayer=1,nlayer
        tmp1(ilayer)=zlay(ilayer)/1000.E+0
      ENDDO
      call profile(nlayer+1,tmp1,tmp2)

      tsurf=tmp2(0)
      DO ilayer=1,nlayer
        temp(ilayer)=tmp2(ilayer)
      ENDDO

!   Initialize soil properties and temperature
!   ------------------------------------------
      volcapa=1.e6 ! volumetric heat capacity

      lecttsoil=0 ! default value for lecttsoil
      call getin("lecttsoil",lecttsoil)
      if (lecttsoil==1) then
         OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
         DO isoil=1,nsoil
            READ (14,*) tsoil(isoil)
            inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
         ENDDO
         GOTO 201
101      STOP'fichier proftsoil inexistant'
201      CONTINUE
         CLOSE(14)

      else 
        DO isoil=1,nsoil
         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
         tsoil(isoil)=tsurf  ! soil temperature
        ENDDO
      endif

!   Initialize depths
!   -----------------
      do isoil=0,nsoil-1
        mlayer(isoil)=2.e-4*(2.**(isoil-0.5)) ! mid-layer depth
      enddo
      do isoil=1,nsoil
        layer(isoil)=2.e-4*(2.**(isoil-1)) ! layer depth
      enddo

c   Initialisation pour les sorties GRADS dans "g1d.dat" et "g1d.ctl"
c   ----------------------------------------------------------------
c   (effectuee avec "writeg1d" a partir de "physiq.F"
c   ou tout autre subroutine physique, y compris celle ci).

      g1d_nlayer=nlayer
      g1d_nomfich='g1d.dat'
      g1d_unitfich=40
      g1d_nomctl='g1d.ctl'
      g1d_unitctl=41
      g1d_premier=.true.
      g2d_premier=.true.

!     Initialize haze profile
!     ------------------------------------------
      if (haze) then

        lecthaze=0 ! default value for lecthaze
        call getin("lecthaze",lecthaze)
        if (lecthaze==1) then
         OPEN(15,file='profhaze',status='old',form='formatted',err=301)
         DO iq=1,nq
          if (iq.eq.i_haze) then
            DO ilayer=1,nlayer
              READ (15,*) q(ilayer,iq)
            ENDDO
          endif
         ENDDO
         GOTO 401
301      STOP'fichier profhaze inexistant'
401      CONTINUE
         CLOSE(15)
        endif
      endif

c   Ecriture de "startfi"
c   --------------------
c   (Ce fichier sera aussitot relu au premier
c   appel de "physiq", mais il est necessaire pour passer
c   les variables purement physiques a "physiq"...

      call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
     .              dtphys,float(day0),
     .              time,tsurf,tsoil,emis,q2,qsurf,
     .              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)

c=======================================================================
c=======================================================================
c=======================================================================
c  BOUCLE TEMPORELLE DU MODELE 1D 
c=======================================================================
c=======================================================================
c=======================================================================

      firstcall=.true.
      lastcall=.false.

      DO idt=1,ndt
        IF (idt.eq.ndt-day_step-1) then       !test
         lastcall=.true.
         call solarlong(day*1.0,zls)
         write(103,*) 'Ls=',zls*180./pi
         write(103,*) 'Lat=', lati(1)*180./pi
         write(103,*) 'RunEnd - Atmos. Temp. File'
         write(103,*) 'RunEnd - Atmos. Temp. File'
         write(104,*) 'Ls=',zls*180./pi
         write(104,*) 'Lat=', lati(1)
         write(104,*) 'RunEnd - Atmos. Temp. File'
        ENDIF

c    calcul du geopotentiel 
c     ~~~~~~~~~~~~~~~~~~~~~
        DO ilayer=1,nlayer
            s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
            h(ilayer)=cpp*temp(ilayer)/(pks*s(ilayer))
        ENDDO

        phi(1)=pks*h(1)*(1.E+0-s(1))
        DO ilayer=2,nlayer
          phi(ilayer)=phi(ilayer-1)+
     &               pks*(h(ilayer-1)+h(ilayer))*.5E+0
     &                  *(s(ilayer-1)-s(ilayer))
        ENDDO

c       appel de la physique
c       --------------------
        CALL physiq (1,llm,nqmx,
     ,     firstcall,lastcall,
     ,     day,time,dtphys,
     ,     plev,play,phi,
     ,     u, v,temp, q,  
     ,     w,
C - sorties
     s     du, dv, dtemp, dq,dpsurf,tracerdyn)


c       evolution du vent : modele 1D
c       -----------------------------
 
c       la physique calcule les derivees temporelles de u et v.
c       on y rajoute betement un effet Coriolis.
c
c       DO ilayer=1,nlayer
c          du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
c          dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
c       ENDDO

c       Pour certain test : pas de coriolis a l'equateur
c       if(lati(1).eq.0.) then
        DO ilayer=1,nlayer
             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
        ENDDO
c       end if
c      
c       Calcul du temps au pas de temps suivant
c       ---------------------------------------
        firstcall=.false.
        time=time+dtphys/daysec
        IF (time.gt.1.E+0) then
            time=time-1.E+0
            day=day+1
        ENDIF

c       calcul des vitesses et temperature au pas de temps suivant
c       ----------------------------------------------------------

        DO ilayer=1,nlayer
           u(ilayer)=u(ilayer)+dtphys*du(ilayer)
           v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
           temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
        ENDDO

c       calcul des pressions au pas de temps suivant
c       ----------------------------------------------------------

        psurf=psurf+dtphys*dpsurf   ! evolution de la pression de surface
        DO ilevel=1,nlevel
          plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
        ENDDO
        DO ilayer=1,nlayer
          play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
        ENDDO

c       calcul traceur au pas de temps suivant
c       --------------------------------------
        DO iq = 1, nqmx
          DO ilayer=1,nlayer
             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
          ENDDO
        END DO

      ENDDO   ! fin de la boucle temporelle

c    ========================================================
c    GESTION DES SORTIE
c    ========================================================

c    fermeture pour conclure les sorties format grads dans "g1d.dat"
c    et "g1d.ctl" (effectuee avec "writeg1d" a partir de "physiq.F"
c    ou tout autre subroutine physique, y compris celle ci).

! save atm temperature profile
      if(saveprofile)then
         OPEN(12,file='profile.out',form='formatted')
         write(12,*) tsurf
         DO ilayer=1,nlayermx
	    write(12,*) temp(ilayer)
         ENDDO
         CLOSE(12)
      endif
! save haze profile
      if (haze.and.lecthaze.eq.1) then 
         OPEN(16,file='profhaze.out',form='formatted')
         DO iq=1,nq
          if (iq.eq.i_haze) then
            DO ilayer=1,nlayer
              write(16,*) q(ilayer,iq)
            ENDDO
          endif
         ENDDO
         CLOSE(16)
      endif

      ! here we compute altitude CORRECTLY!!!

!      print*,'zlay=',zlay/1000.
!      rho = play(1)/(R*tsurf)
!      zlay(1)=(plev(1)-plev(2))/(g*rho)
!      do ilayer=2,nlayer
!         rho  = play(ilayer)/(R*temp(ilayer))
!         rho  = play(ilayer)/(R*230.0)
!         dz = (plev(ilayer-1)-plev(ilayer))/(g*rho)
!         zlay(ilayer)=zlay(ilayer-1)+dz
!      enddo
!      print*,'zlay=',zlay/1000.
!      stop

c        CALL endg1d(1,nlayer,zphi/(g*1000.),ndt)
      do ilayer=1,nlayer
        zlay(ilayer)=-18.0e3*log(play(ilayer)/psurf)
      enddo
      CALL endg1d(1,nlayer,zlay/1000.,ndt)

c    ========================================================
      END
 
c***********************************************************************
c***********************************************************************
c     Subroutines Bidons utilise seulement en 3D, mais
c     necessaire a la compilation de testphys1d en 1-D

      subroutine gr_fi_dyn
      RETURN
      END
 
c***********************************************************************
c***********************************************************************

#include "../dyn3d/disvert.F"
#include "call_dayperi.F"