PROGRAM testphys1d
      IMPLICIT NONE

c=======================================================================
c   subject:
c   --------
c   PROGRAM useful to run physical part of the venusian GCM in a 1D column
c       
c Can be compiled with a command like (e.g. for 50 layers)
c  "makegcm -p venus -d 50 testphys1d"

c A REVOIR POUR VENUS....
c It requires the files "testphys1d.def" "callphys.def"
c      and a file describing the sigma layers (e.g. "z2sig.def")
c
c   author: Frederic Hourdin, R.Fournier,F.Forget (original Mars version)
c   ------- Sebastien Lebonnois (Venus version)
c   
c=======================================================================

c VENUS: revoir tous les include...... par rapport a calfis.F
c        voir aussi les arguments de calfis: clesphy0
c        voir aussi dans les include de physiq.F les variables qui 
c             doivent etre initialisees...
#include "dimensions.h"
#include "dimsoil.h"
#include "comcstfi.h"
#include "control.h"
#include "comvert.h"
#include "netcdf.inc"
#include "comg1d.h"
#include "logic.h"
#include "clesphys.h"
#include "iniprint.h"

c --------------------------------------------------------------
c  Declarations
c --------------------------------------------------------------
c
      INTEGER unit           ! unite de lecture de "testphys1d.def"
      INTEGER unitstart      ! unite d'ecriture de "startphy.nc"
      INTEGER nlayer,nlevel,nsoil,ndt
      INTEGER ilayer,ilevel,isoil,idt,iq
      LOGICAl firstcall,lastcall
c
      INTEGER day0          ! date initial (sol ; =0 a Ls=0)
      REAL day              ! date durant le run
      REAL time             ! time (0<time<1 ; time=0.5 a midi)
      REAL play(llm)   ! Pressure at the middle of the layers (Pa)
      REAL plev(llm+1) ! intermediate pressure levels (pa)
      REAL psurf,tsurf(1)      
      REAL u(llm),v(llm)  ! zonal, meridional wind
      REAL gru,grv   ! prescribed "geostrophic" background wind
      REAL temp(llm)   ! temperature at the middle of the layers
      REAL q(llm,nqtot) ! tracer mixing ratio (e.g. kg/kg)
      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
      REAL zlay(llm)   ! altitude estimee dans les couches (km)
      REAL long(1),lati(1),area(1)
      REAL cufi(1),cvfi(1)
      REAL phisfi(1),albedo(1)
      REAL solsw(1),sollwdown(1),dlw(1),radsol(1)

c    Physical and dynamical tandencies (e.g.  m.s-2, K/s, Pa/s)
      REAL du(llm),dv(llm),dtemp(llm)
      REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
      REAL dpsurf    
      REAL dq(llm,nqtot)
      REAL dqdyn(llm,nqtot)

c   Various intermediate variables
      REAL zls
      REAL phi(llm),s(llm),aps(llm),bps(llm)
      REAL pk(llm),pks, w(llm)
      INTEGER l, ierr, aslun,radpas
      REAL tmp1(0:llm),tmp2(0:llm)
      INTEGER        longcles
      PARAMETER    ( longcles = 20 )
      REAL clesphy0( longcles      )
                        

      character*2 str2

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

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

      lunout = 6

c ------------------------------------------------------
c  Constantes prescrites ICI
c ------------------------------------------------------

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

c     Constante de la Planete Venus
c     -----------------------------
      rad=6051300.               ! rayon de Venus (m)  ~6051300 m
      daysec= 1.0087e7           ! duree du sol (s)  ~1.e7 s
      omeg=4.*asin(1.)/19.4141e6 ! vitesse de rotation (rad.s-1)
      g= 8.87                    ! gravite (m.s-2) ~8.87
      mugaz=43.44                ! Masse molaire de l'atm (g.mol-1) ~43.44
! ADAPTATION GCM POUR CP(T)
! VENUS: Cp(T) = cpp*(T/T0)^nu 
! avec T0=460. et nu=0.35
!     cpp=1.0e3
      cpp=9.0e2      ! version constante
      r= 8.314511E+0 *1000.E+0/mugaz
      rcp= r/cpp

c ------------------------------------------------------
c  Lecture des parametres dans "testphys1d.def" 
c ------------------------------------------------------

c   Opening parameters file "testphys1d.def"
c   ---------------------------------------
      unit =97
      OPEN(unit,file='testphys1d.def',status='old',form='formatted'
     .     ,iostat=ierr)

      IF(ierr.ne.0) THEN
        write(*,*) 'Problem to open "testphys1d.def'
        write(*,*) 'Is it there ?'
        stop
      END IF

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
      PRINT *,'date de depart ?'
      READ(unit,*) day0
      day=float(day0)
      PRINT *,day0
c  Heure de demarrage
      PRINT *,'heure de debut de simulation (entre 0 et 24) ?'
      READ(unit,*) time
      time=time/24.E+0

c  Discretisation (Definition de la grille et des pas de temps)
c  --------------
c
      nlayer=llm
      nlevel=nlayer+1
      nsoil=nsoilmx
      PRINT *,'nombre de pas de temps par jour ?'
      READ(unit,*) day_step
      print*,day_step

      PRINT *,'nombre d appel au rayonnement par jour ?'
      READ(unit,*) nbapp_rad
      print*,nbapp_rad

      PRINT *,'nombre de pas physique entre chaque sortie ?'
      READ(unit,*) nbphypers
      print*,nbphypers

      PRINT *,'nombre de jours simules ?'
      READ(unit,*) ndt
      print*,ndt

      ndt=ndt*day_step     
      dtphys=daysec/day_step  

c Pression de surface sur la planete
c ------------------------------------
c
      PRINT *,'pression au sol'
      READ(unit,*) psurf
      PRINT *,psurf
c Pression de reference  ! voir dyn3d/etat0_venus
c     pa     =  5.e4 
      pa     =  1.e6 
      preff  = 9.2e6 ! 92 bars
c     preff  = psurf
 
c  latitude/longitude 
c  -------------------
      PRINT *,'latitude en degres ?'
      READ(unit,*) lati(1)
      PRINT *,lati(1)
      long(1)=0.E+0

c  Initialisation albedo 
c  ----------------------
c
      PRINT *,'Albedo du sol nu ?'
      READ(unit,*) albedo(1)
      PRINT *,albedo(1)

c   Initialisation speciales "physiq"
c   ---------------------------------
c   la surface de chaque maille est inutile en 1D --->
      area(1)=1.E+0
c de meme ?
      cufi(1)=1.E+0
      cvfi(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

      CALL iniphysiq(1,llm,daysec,day0,dtphys,
     .            lati,long,area,cufi,cvfi,rad,g,r,cpp)

c   Initialisation pour prendre en compte les vents en 1-D
c   ------------------------------------------------------
 
c    vent geostrophique
      PRINT *,'composante vers l est du vent geostrophique (U) ?'
      READ(unit,*) gru
      PRINT *,'composante vers le nord du vent geostrophique (V) ?'
      READ(unit,*) grv

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

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

c    Vertical Coordinates  (hybrids)
c    """"""""""""""""""""
      CALL  disvert
      
c     Calcul au milieu des couches : Vient de la version Mars
c     WARNING : le choix de placer le milieu des couches au niveau de
c     pression intermédiaire est arbitraire et pourrait etre modifié.
c     C'est fait de la meme facon dans disvert

      DO l = 1, llm
       aps(l) =  0.5 *( ap(l) +ap(l+1))
       bps(l) =  0.5 *( bp(l) +bp(l+1))
      ENDDO

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

      DO ilayer=1,nlayer
        play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
        pk(ilayer)  =cpp*(play(ilayer)/preff)**rcp
c       write(120,*) ilayer,plev(ilayer),play(ilayer)
      ENDDO
c     write(120,*) nlevel,plev(nlevel)
c     stop
      
      pks=cpp*(psurf/preff)**rcp

c  profil de temperature et altitude au premier appel
c  --------------------------------------------------

c modif par rapport a Mars: 
c   on envoie dz/T=-log(play/psurf)*r/g dans profile
      tmp1(0)=0.0
      tmp1(1)= -log(play(1)/psurf)*r/g
      DO ilayer=2,nlayer
        tmp1(ilayer)=-log(play(ilayer)/play(ilayer-1))*r/g
      ENDDO
      call profile(unit,nlayer+1,tmp1,tmp2)
c file testphys1d.def closed in profile!

      print*,"               Pression        Altitude     Temperature"
      ilayer=1
      tsurf(1)=tmp2(0)
       temp(1)=tmp2(1)
       zlay(1)=tmp2(1)*tmp1(1)
      print*,"           0",tsurf(1)
      print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
      DO ilayer=2,nlayer
        temp(ilayer)=tmp2(ilayer)
        zlay(ilayer)=zlay(ilayer-1)+tmp2(ilayer)*tmp1(ilayer)
        print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
      ENDDO
      
c     temperature du sous-sol
c     ~~~~~~~~~~~~~~~~~~~~~~~
      DO isoil=1,nsoil
         tsoil(isoil)=tsurf(1)
      ENDDO

c    Initialisation des traceurs
c    ---------------------------

      DO iq=1,nqtot
        DO ilayer=1,nlayer
           q(ilayer,iq) = 0.
        ENDDO
      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.

c  Ecriture de "startphy.nc"
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"...

      clesphy0(1) = 0.
      clesphy0(2) = nbapp_rad
      clesphy0(3) = 0. ! cycle diurne
      clesphy0(4) = 1. ! soil_model
      clesphy0(5) = 1. ! new_oliq
      clesphy0(6) = 0.
      clesphy0(7) = 0.
      clesphy0(8) = 0.

      solsw(1)    = 0.
      sollwdown(1)= 0.
      dlw(1)      = 0.
      radsol(1)   = 0.
     
      radpas      = NINT(day_step/clesphy0(2))
      soil_model  = .true.
      new_oliq    = .true.

      call phyredem("startphy.nc  ",dtphys,radpas,
     .              lati,long,
     .              tsurf,tsoil,albedo,
     .              solsw,sollwdown,dlw,radsol,
     .              temp,q)

c=======================================================================
c  BOUCLE TEMPORELLE DU MODELE 1D 
c=======================================================================
c
      firstcall=.true.
      lastcall=.false.

      DO idt=1,ndt
c        IF (idt.eq.ndt) lastcall=.true.
        IF (idt.eq.ndt-day_step-1) then       !test
         lastcall=.true.
c toujours nulle dans le cas de Venus, pour l'instant...
         zls = 0.0
c        write(103,*) 'Ls=',zls*180./pi
c        write(103,*) 'Lat=', lati(1)
c        write(103,*) 'RunEnd - Atmos. Temp. File'
c        write(103,*) 'RunEnd - Atmos. Temp. File'
c        write(104,*) 'Ls=',zls*180./pi
c        write(104,*) 'Lat=', lati(1)
c        write(104,*) 'RunEnd - Atmos. Temp. File'
        ENDIF

c    calcul du geopotentiel 
c     ~~~~~~~~~~~~~~~~~~~~~
! ADAPTATION GCM POUR CP(T)
      DO ilayer=1,nlayer
        s(ilayer)=(play(ilayer)/psurf)**rcp
      ENDDO
      phi(1)=cpp*temp(1)*(1.E+0-s(1))
      DO ilayer=2,nlayer
         phi(ilayer)=phi(ilayer-1)+
     &     cpp*(temp(ilayer-1)/s(ilayer-1)+temp(ilayer)/s(ilayer))*0.5
     &        *(s(ilayer-1)-s(ilayer))

      ENDDO

c       appel de la physique
c       --------------------

      CALL physiq (1,llm,nqtot,
     ,     firstcall,lastcall,
     ,     day,time,dtphys,
     ,     plev,play,pk,phi,phisfi,
     ,     presnivs,clesphy0,
     ,     u,v,temp,q,  
     ,     w,
C - sorties
     s     du,dv,dtemp,dq,dpsurf)

c     print*,"DT APRES PHYSIQ=",day,time
c     print*,dtemp
c     print*,temp
c     print*," "
c     stop

c       evolution du vent : modele 1D
c       -----------------------------
 
c       la physique calcule les derivees temporelles de u et v.
c       Pas de coriolis 
          DO ilayer=1,nlayer
             du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
             dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
          ENDDO
c      
c         call writeg1d(1,nlayer,du,"du","du")
c         call writeg1d(1,nlayer,dv,"dv","dv")
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         call writeg1d(1,nlayer,u,"u","Vent zonal")
c         call writeg1d(1,nlayer,v,"v","Vent meridien")

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

           call writeg1d(1,nlayer,play,"press","Pressure")
           call writeg1d(1,nlayer,phi,"phi","Geopot")

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

        CALL endg1d(1,nlayer,zlay/1000.,ndt)
        print*,"Temperature finale:"
        print*,temp
        
c stabilite
      DO ilayer=1,nlayer
        zlay(ilayer) = phi(ilayer)/8870.  !en km
      ENDDO
      DO ilayer=2,nlayer
        tmp1(ilayer) =
     .    (temp(ilayer)-temp(ilayer-1))/(zlay(ilayer)-zlay(ilayer-1)) 
     .   + 8870./cpp
      ENDDO

      OPEN(11,file='profile.new')
      write (11,*) tsurf
      DO ilayer=1,nlayer
        write (11,*) zlay(ilayer),temp(ilayer),tmp1(ilayer)
      ENDDO

c    ========================================================
      END
 
c***********************************************************************
c***********************************************************************

#include "../dyn3d/disvert.F"
#include "../dyn3d/abort_gcm.F"
#include "../dyn3d/dump2d.F"