PROGRAM gcm1d
      IMPLICIT NONE

#include "dimensions.h"
#include "dimphy.h"
#include "YOMCST.h"
#include "comg1d.h"
#include "clesphys.h"
#include "control.h"

c    Arguments :
c    -----------

      integer ngrid,nlayer,longcles,nqmax
      parameter (longcles=20,nqmax=3)
cfleur on  ne se sert pas de nqmax mais de nqmx defini dnas dimensions.h depend de la 
c      compilation
      integer radpas
      integer nday,it,iflag_con,unit,i,itap,n_cooling
      integer iphys_ver,iadv_tvl,i_cvg,i_hum

      real fnday
      real h,kappa,z1,z2,sbid,sigbid
      real plev(klev+1),play(klev), psol
      real temp(klev),u(klev),v(klev),tsurf,q(klev,nqmx),w(klev+1)
      real ug(klev),vg(klev)
      real du(klev),dv(klev),dt(klev),dpsrf,dq(klev,nqmx)
      real du_dyn(klev),dv_dyn(klev)
     :    ,dt_dyn(klev),dq_dyn(klev,nqmax)
      real phis,presnivs(klev),clesphy0(longcles)
      real time,timestep,ecritphy,day,tho
      real co2_ppm,solaire
      real rlat,rlon,tsol,radsol,psol_f,tsol_f,qsol_f,qsol
      real rugmer,rugsrel,snow,agesno,deltat,zmea,zstd,zsig,sn
      real zgam,zthe,zpic,zval
      real ema_sig(klev),ema_w(klev),ncst_cbmf,ema_cbmfz,ema_pcb
     : ,zz_f(klev),vu_f(klev),vv_f(klev),t_f(klev),q_f(klev,nqmax)

      real temp0(klev),q0(klev,3)

      real dt_cooling(klev),dq_cooling(klev)
      real d_t_cool(klev),d_q_cool(klev)
      real d_t_adv(klev),d_q_adv(klev,nqmax)
      real d_t_cvg(klev),d_q_cvg(klev)
      real ht(100),hq(100),hw(100)

      real phy_nat(360)
      real phy_alb(360)
      real phy_sst(360)
      real phy_bil(360)
      real phy_rug(360)
      real phy_ice(360)

      logical cycle_diurne,soil_model,new_oliq,ok_orodr
     :       ,ok_orolf,ok_limitvrai

      logical firstcall,lastcall,flag_cool
      logical itsourcecont
c      itsourcecont permet de choisir entre les sources pour plus de 5 traceurs
      character*80 ans,file_forctl, file_fordat, file_start,file
      character*2 cnbl

c-----------------------------------------------------------------------

      COMMON/comvert/
     * s(llm),sig(llm+1),ds(llm),dsig(llm),dsig1(llm),sdsig(llm)
     . ,sig_s(llm)

      REAL s,sig,ds,dsig,dsig1,sdsig,sig_s

c-----------------------------------------------------------------------
c INCLUDE 'temps.h':

      COMMON/temps/itaufin,dtd,
     s  day_ini,day_end,anne_ini

      INTEGER  itaufin
      INTEGER*4 day_ini,day_end,anne_ini
      REAL     dtd

c-----------------------------------------------------------------------
c   dynamical tendencies

      INTEGER l,ierr,aslun,nlevel,iq,ll

      REAL longitude,latitude
      REAL zlay(klev),phi(klev)
      REAL paire 

      DATA latitude,longitude/0.,0./

c-----------------------------------------------------------------------
c    Initialisations  des constantes
c    -------------------------------

c constantes
c ----------

      time=0.
      tho=3600. 
      it=0

      call suphec

c parametres lus dans execution_1D:
c ---------------------------------

      read(*,*) fnday
      print*, 'fnday',fnday
      read(*,*) ecritphy
      print*, 'ecritphy',ecritphy
      read(*,*) timestep
      print*, 'timestep',timestep
      read(*,*) iflag_con 
      print*, 'iflag_con',iflag_con
      read(*,*) flag_cool 
      print*, 'flag_cool',flag_cool
      read(*,'(a)') cnbl
      print*, 'cnbl',cnbl
      read(*,*) radpas   
      print*, 'radpas',radpas
      read(*,*) anne_ini ! sans importance mais il faut qqchose
      print*, 'anne_ini',anne_ini
      
c si ecritphy=0: on ecrit la physique a chaque pas de temps:

      if (ecritphy.eq.0.) then
       ecritphy = timestep/rday
      endif
      write(*,*) 'ECRITPHY = ',ecritphy

      firstcall=.true.
      lastcall=.false.
      day=1.
      itsourcecont=.true.

      open(78,file='transil.dat',form='unformatted',
     s    access='direct',recl=4*llm*llm)

      ngrid=1
      IF (ngrid.NE.klon) THEN
         PRINT*,'STOP in inifis'
         PRINT*,'Probleme de dimenesions :'
         PRINT*,'ngrid     = ',ngrid
         PRINT*,'klon   = ',klon
         STOP
      ENDIF
      nlayer=klev
      nlevel=nlayer+1

c   -------------------------------------------------------------------
c   Initialisation de la discretisation verticale
c   ---------------------------------------------
c sb      CALL disvert(llm,rkappa,sig,dsig,s,ds,dsig1,sdsig)

c   -------------------------------------------------------------------
c   Profils initiaux.
c   -----------------

      unit = 80
      open(unit,file='start'//cnbl//'.data',form='formatted')
      read(unit,*) 
      read(unit,*) 
      read(unit,*) (play(l),l=1,nlayer)
      read(unit,*) (plev(l),l=1,nlevel)
      read(unit,*) (temp(l),l=1,nlayer)
      read(unit,*) (q(l,1),l=1,nlayer)
      read(unit,*) (q(l,2),l=1,nlayer)
      read(unit,*) (ema_sig(l),l=1,nlayer)
      read(unit,*) (ema_w(l),l=1,nlayer) ! reinitialise + loin???
      read(unit,*) ncst_cbmf, ema_cbmfz, ema_pcb
      read(unit,*) (zz_f(l),l=1,nlayer)
      read(unit,*) (vu_f(l),l=1,nlayer)
      read(unit,*) (vv_f(l),l=1,nlayer)
      close(unit)
      write(*,*) 'Lecture fichier start.data ok'

      DO l = 1, nlayer
         u(l) = vu_f(l)
         v(l) = vv_f(l)
c##################
cfleur ATTENTION q(3) ce nest plus la glace mais le radon si on active phytrac
c##################3
         q(l,3) = 0. ! on initialise la glace a zero
         PRINT*,plev(l),play(l)
c on recopie le profil initial dans les champs de forcing:
         t_f(l) = temp(l)
         q_f(l,1) = q(l,1)
         q_f(l,2) = q(l,2)
         q_f(l,3) = q(l,3)
      ENDDO


c Lecture/creation des conditions aux limites:
c --------------------------------------------

      unit = 81
      open(unit,file='startphy.data',form='formatted')
      read(unit,*)
      read(unit,*) co2_ppm
      read(unit,*) solaire
      read(unit,*) rlat
      read(unit,*) tsol
      read(unit,*) radsol
      close(unit)
      write(*,*) 'Lecture fichier startphy.data ok'

      unit = 82
      open(unit,file='condsurf.data',form='formatted')
      read(unit,*) psol_f
      read(unit,*) tsol_f
      read(unit,*) qsol_f
      close(unit)
      write(*,*) 'Lecture fichier condsurf.data ok'

      if (play(1).lt.10000.) then
        do l = 1, nlayer
         play(l) = play(l)*100.
        enddo
      endif
      if (plev(1).lt.10000.) then
        do l = 1, nlevel
         plev(l) = plev(l)*100.
        enddo
      endif

      if (abs(psol_f-plev(1)) .gt. 1.) then
         print *,' Incompatibilite entre psol et profil'
     :   ,' de pression'
         print *,' psol = ',psol_f,' plev(1) = ',plev(1)
         stop
      endif

      tsol    = tsol_f ! temp au sol prise dans condsurf.data
      psol    = psol_f ! pression au sol prise dans condsurf.data
      qsol    = qsol_f
      rugmer  = 0.0001 ! valeur de cchlim.data
      rugsrel = 0.0    ! (rugsrel = rugoro)
      snow    = 0.0
      sn=0.
      agesno  = 50.0
      rlon    = 0.0
      deltat  = 0.0    ! ne sert que pour les slab_ocean
      phis = 0.
      zmea = 0.
      zstd = 0.
      zsig = 0.
      zgam = 0.
      zthe = 0.
      zpic = 0.
      zval = 0.

      do i=1,360
      phy_sst(i) = tsol
      phy_nat(i) = 0.0  ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
      phy_alb(i) = 0.15 ! albedo land only (old value condsurf_jyg=0.3)
      phy_bil(i) = 1.0  ! ne sert que pour les slab_ocean
      phy_rug(i) = 0.1  ! longueur rugosite utilisee sur land only 
      phy_ice(i) = 0.0  ! fraction de glace (?)
      enddo

      call writelim
     s   (phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice)



c controles du run (en 3D: lus dans run.def) :

      cycle_diurne = .FALSE.
      soil_model   = .FALSE.
      new_oliq     = .FALSE.
      ok_orodr     = .FALSE.
      ok_orolf     = .FALSE.
      ok_limitvrai = .FALSE.

      do i = 1, longcles
       clesphy0(i) = 0.
      enddo

      nbapp_rad = NINT(86400./radpas/timestep)

      clesphy0(1) = FLOAT( iflag_con )
      clesphy0(2) = FLOAT( nbapp_rad )
      if (cycle_diurne ) clesphy0(3) = 1.
      if (  soil_model ) clesphy0(4) = 1.
      if (    new_oliq ) clesphy0(5) = 1.
      if (    ok_orodr ) clesphy0(6) = 1.
      if (    ok_orolf ) clesphy0(7) = 1.
      if ( ok_limitvrai) clesphy0(8) = 1.


c -------------------------------------------------------------------
c Discretisation verticale:
c -------------------------

c calcul sig,dsig,s,ds,dsig1,sdsig a partir des play, plev lus:

      do l = 1, nlevel
       sig(l)= plev(l)/plev(1)
      enddo

      do l = 1, nlayer
       s(l)= sig(l)**RKAPPA
      enddo

      do l = 2, nlayer
       ds(l) = s(l-1) - s(l)
      enddo
      ds(1)  = 1. - s(1)

      do l = 1, nlayer
       dsig(l)  = sig(l)-sig(l+1)
       sdsig(l) = s(l) * dsig(l)
       dsig1(l) = 1./dsig(l)
       presnivs(l) = play(l)/100./100.
      enddo

      print*,'Diagnostique de la discretisation verticale'

      print*
      h=7.
      kappa = RKAPPA
      print*,'comparaison de sig(l) et (s(l)+s(l+1))/2)**(1/K)'
      do 14 l=1,llm-1
         sigbid=(0.5*(s(l)+s(l+1)))**(1./kappa)
         print*,'sig(',l+1,')  = ',sig(l+1),
     :           '    valeur approchee :',sigbid,'   ',dsig(l)
14    continue
      print*
      print*,'comparaison de s(l) et (sig(l)+sig(l+1))/2)**K'
      do 15 l=1,llm
         sbid=(0.5*(sig(l+1)+sig(l)))**kappa
         print*,'  s(',l,')  = ',s(l),
     :           '    valeur approchee :',sbid
15    continue

      print*,'Altitude approchee z,dz'
      print*
      z1=0.
      print*,'   l       Z      DZ      Ztop   dsig'
      DO 18 l=1,llm-1
         z2=-h*log(sig(l+1))
         write(*,'(i5,3x,4f8.4)') l,-h*log(s(l))/kappa,z2-z1,z2
     :    ,dsig(l)
         write(14,'(3x,i5,1f10.4)') l,-h*log(s(l))/kappa
         z1=z2
18    CONTINUE
      write(*,'(i5,3x,3f8.4)') l,-h*log(s(llm))/kappa
      write(14,'(3x,i5,1f10.4)') l,-h*log(s(llm))/kappa

c      print*,'Correspondance approx pression-altitude:'
c      z1 = 0.
c      do l = 1, llm
c       z1 = z1 + 287.*temp(l)/9.81/play(l)*(plev(l)-plev(l+1))
c       write(*,*) l,play(l),z1
c      enddo

c -------------------------------------------------------------------
c Ecriture de l'etat initial:
c ---------------------------

c     call ini_fis(timestep,radpas,iflag_con,anne_ini
c    :            ,co2_ppm,solaire,rlat,rlon,tsol,deltat
c    :            ,qsol,snow,radsol,rugmer,agesno,zmea,zstd,zsig,zgam
c    :            ,zthe,zpic,zval,rugsrel
c    :            ,phy_sst,phy_nat,phy_alb,phy_bil,phy_rug,phy_ice)

      call physdem(rlon, rlat, timestep,radpas,co2_ppm,
     .                   solaire,tsol, qsol,
     .                   sn, radsol, deltat, rugmer,
     .                   agesno, zmea, zstd, zsig,
     .                   zgam, zthe, zpic, zval,
     .                   rugsrel)



c -------------------------------------------------------------------
c Options de la simulation:
c -------------------------

      unit = 83
      open(unit,file='version.data',form='formatted')
      read(unit,*) iphys_ver,iadv_tvl,i_cvg,i_hum
      close(unit)
      write(*,*) 'Lecture fichier version.data ok'

      write(*,*) 'physiq # ',iphys_ver,' (0=std; 1=forced glob cis;
     : 2=forced loc cis)'
      write(*,*) 'advection switch : ',iadv_tvl,' (0=no adv ;
     : 1=forced adv)'
      write(*,*) 'convergence switch : ',i_cvg,' (0=no conv;
     : 1=forced conv)'



c Eventuellement, preparation de la "bulle froide":
c -------------------------------------------------

      IF (flag_cool) THEN
      unit = 84
      open(unit,file='cool_buble.data',form='formatted')
      read(unit,*) n_cooling
      do l = 1, nlayer
      read(unit,*) dt_cooling(l),dq_cooling(l)
      enddo
      close(unit)
      write(*,*) 'Lecture fichier cool_buble.data ok'
      ENDIF


c Eventuellement, preparation du forcage par la convergence:
c ----------------------------------------------------------

      IF (i_cvg .EQ. 1) THEN
      file_forctl = 'forcing.ctl'
      file_fordat = 'forcing.dat'
      file_start  = 'start'//cnbl//'.data'

      call copie(klev,play,psol,file_forctl)
      call get_uvd2(itap,file_forctl,file_fordat,file_start
     :             ,ht,hq,hw)
      ENDIF


c-----------------------------------------------------------------------
c  initialisation pour GRADS-1D
c  ----------------------------

      g1d_nlayer=nlayer
      g1d_nomfich='grads1d.dat'
      g1d_unitfich=30
      g1d_nomctl='grads1d.ctl'
      g1d_unitctl=31
      g1d_premier=.true.

      file='sort'
      call inigrads(1,1
     s  ,0.,1.,-2.,2.,1,0.,-2.,2.,1.
     s  ,llm,presnivs,1000.
     s  ,1800.,file,'Diagconvect')
      print*,'Fin de Initialisation de wrgras'


c=======================================================================
c   DEBUT DE L'INTEGRATION TEMPORELLE:
c   ==================================

      itap = 1

1     continue

c  calcul du geopotentiel:
c  -----------------------

      phi(1)=RD*temp(1)*
     :(plev(1)-play(1))/(.5*(plev(1)+play(1)))
      do l = 1, nlayer-1
         phi(l+1)=phi(l)+RD*(temp(l)+temp(l+1))*
     :   (play(l)-play(l+1))/(play(l)+play(l+1))
      enddo

c        PRINT*,'altitude km  et T (K)'
c        WRITE(6,'(2f10.2)') (.001*phi(l)/g,temp(l),l=1,nlayer)


c pour etre coherent avec LMDZ.3, on passe 2 arguments 
c en plus a physiq: phis et aire (pas utilises en 1D):
      paire = 1. ! aire de la maille


c  eventuellement, induction de la convection par chauffage du
c  bas et refroidissement du haut de la couche limite:
c  ------------------------------------------------------------

      IF (flag_cool) THEN

      call cool_pool(itap
     e           ,n_cooling,dt_cooling,dq_cooling
     s           ,d_t_cool,d_q_cool)

      do l = 1, nlayer
       temp(l) = temp(l) + d_t_cool(l)
       q(l,1)  = q(l,1)  + d_q_cool(l)
      enddo

      ELSE

      do l = 1, nlayer 
       d_t_cool(l) = 0.
       d_q_cool(l) = 0.
      enddo

      ENDIF


c  eventuellement, ajouter une tendance relative a une
c  translation du domaine (ex: pour suivre une ligne de grains):
c  -------------------------------------------------------------

      IF (iadv_tvl .EQ. 1) THEN

      call advect_tvl(timestep,temp,q,vu_f,vv_f,t_f,q_f
     :               ,d_t_adv,d_q_adv)

      do l = 1, nlayer
       temp(l) = temp(l) + d_t_adv(l)
       q(l,1)  = q(l,1)  + d_q_adv(l,1)
       q(l,2)  = q(l,2)  + d_q_adv(l,2)
       q(l,3)  = q(l,3)  + d_q_adv(l,3)
      enddo

      ELSE

      do l = 1, nlayer 
       d_t_adv(l) = 0.
       d_q_adv(l,1) = 0.
       d_q_adv(l,2) = 0.
       d_q_adv(l,3) = 0.
      enddo

      ENDIF


c  eventuellement, ajouter une tendance relative a la
c  convergence de grande echelle:
c  -------------------------------------------------------------

      IF (i_cvg .EQ. 1) THEN

      call get_uvd(itap,timestep,tsol,qsol,file_fordat,ht,hq,hw)

      do l = 1, nlayer

       temp0(l) = temp(l) ! memoire de "l'avant dynamique"
       q0(l,1) = q(l,1)
       q0(l,2) = q(l,2)
       q0(l,3) = q(l,3)

       d_t_cvg(l) = ht(l) * timestep
       d_q_cvg(l) = hq(l) * timestep
       temp(l) = temp(l) + d_t_cvg(l)
       q(l,1)  = q(l,1)  + d_q_cvg(l)
c       write(*,*)'d_t_cvg,d_q_cvg:',d_t_cvg(l),d_q_cvg(l) 
       if (q(l,1).lt.0.) then ! sb
       print*,'OVAP negative dans main!'
       print*,'itap,l,q,d_q_cvg:',itap,l,q(l,1),d_q_cvg(l),hq(l)
       q(l,1)  = MAX(q(l,1),1.e-10) ! evite les humidites negatives
       endif
      enddo

      if (itap.ge.30) then
        print*,'hq(l),ht(l),dq,dt,q:'
     :  ,hq(16),ht(16),d_q_cvg(16),d_t_cvg(16),q(16,1)
      endif

      ELSE

      do l = 1, nlayer
       d_t_cvg(l) = 0.
       d_q_cvg(l) = 0.
      enddo

      ENDIF

      do l = 1, nlayer
       dt_dyn(l) = d_t_cvg(l) / timestep
       dq_dyn(l,1) = d_q_cvg(l) / timestep
       dq_dyn(l,2) = 0.0
       dq_dyn(l,3) = 0.0
      enddo


c   calcul des tendances physiques:
c   -------------------------------

      if (itsourcecont) then
        do iq=5,nqmx
          ll=min(iq-4,18)
          if(firstcall) then 
c          q(ll,iq)=5/((plev(ll)-plev(ll+1))/rg)
c           q(ll,iq)=5 
c          q(ll,iq)=5*play(ll)/((plev(ll)-plev(ll+1))/rg)/(RD*temp(ll))
            q(ll,iq)=5*play(ll)/(RD*temp(ll))
          endif
        print*,'q ll iq',q(ll,iq),ll,iq
        enddo
      else
        do iq=5,nqmx
          ll=min(iq-4,18)
          q(ll,iq)=q(ll,iq)+1
             do l=1,llm
               q(l,iq)=q(l,iq)*(1-timestep/tho)
        print*,'q ll iq',q(ll,iq),ll,iq
             enddo   
        enddo
 
      endif

c      CALL physiq(ecritphy,ngrid,nlayer,nqmx,
c     :            firstcall,lastcall,
c     :            day,day,time,timestep,
c     :            plev,play,phi,phis,paire,
c     :            u,v,temp,q,
c     :            du_dyn, dv_dyn,  dt_dyn, dq_dyn,
c     :            w,
c     :            du,dv,dt,dq,dpsrf)

c     CALL physiq (ecritphy,ngrid,nlayer,nqmx,firstcall,lastcall,
c    ,   day,day,time,timestep,plev,play,phi,phis,paire,
c     ,             presnivs,clesphy0,u,v,temp,q,  
c pour calculer proprement la tendance de cape liee a la
c dynamique, il faut entrer aussi temp0 et q0:
c    ,             presnivs,clesphy0,u,v,temp,q,temp0,q0,  
c    ,             du_dyn,dv_dyn,dt_dyn,dq_dyn,w,
c - sorties
c    s             du,dv,dt,dq,dpsrf)

      print*,'PAS DE TEMPS ',timestep
      print*,'ATTENTION!!! Il faudra passer temp0 et q0'
      CALL physiq(ngrid,nlayer,nqmx,
     :            firstcall,lastcall,
     :            day,day,time,timestep,
     :            plev,play,phi,phis,paire,presnivs,clesphy0,
     :            u,v,temp,q,
     :            w,
     :            du,dv,dt,dq,dpsrf)




      firstcall=.false.


c   Ajout des tendances
c   -------------------
      DO l=1,nlayer
         u(l)=u(l)+timestep*du(l)
         v(l)=v(l)+timestep*dv(l)

         temp(l)=temp(l)+timestep*dt(l)
      ENDDO

      do iq=1,nqmx
         do l=1,nlayer
            q(l,iq)=q(l,iq)+timestep*dq(l,iq)
         enddo
      enddo
c      write(78,rec=itap) ((q(l,iq),l=1,llm),iq=5,22)
      write(78,rec=itap) ((q(l,iq),iq=5,22),l=1,llm)
      itap = itap + 1


         time=time+timestep/rday
         it=it+1
         if(mod(it,2000).eq.0) print*,'TIME=',time
         if(time.gt.fnday) then
            CALL endg1d(1,nlayer,play,int(time/ecritphy),timestep)
            stop
         endif

         GOTO 1

      END

c=======================================================================
c=======================================================================
c    FIN DU PROGRAMMES
c    CI-DESSOUS, QUELQUES SOUS-PROGRAMMES UTILS
c=======================================================================
c=======================================================================

#include "1DUTILS.h"
#include "1Dconv.h"