SUBROUTINE etat0_netcdf
    
      USE startvar
      USE ioipsl
      !
      IMPLICIT NONE
      !
#include "dimensions.h"
#include "paramet.h"
      !
      !
      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 
     .KLON=KFDIA-KIDIA+1,KLEV=llm
      !
#include "comgeom2.h"
#include "comvert.h"
#include "comconst.h"
      !
      REAL :: latfi(klon), lonfi(klon)
      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
     . psol(iip1, jjp1), phis(iip1, jjp1)
      REAL :: p3d(iip1, jjp1, llm+1)
      REAL :: uvent(iip1, jjp1, llm)
      REAL :: vvent(iip1, jjm, llm)
      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
      REAL :: tsol(klon), qsol(klon), sn(klon), radsol(klon)
      REAL :: deltat(klon), rugmer(klon), agesno(klon)
      REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
      REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1)
      REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1)
      REAL :: qd(iip1, jjp1, llm)
      !
      CHARACTER*80 :: varname
      !
      INTEGER :: i,j, ig, l
      REAL :: xpi
      !
      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
      REAL :: workvar(iip1,jjp1,llm)
      !
      REAL ::  prefkap, unskap
      !
      REAL :: q_sat
      EXTERNAL q_sat
      real :: time_step,t_ops,t_wrt

#include "comdissnew.h"
#include "control.h"
#include "serre.h"
#include "clesph0.h"

      INTEGER  ::        longcles
      PARAMETER      ( longcles  = 20 )
      REAL :: clesphy0 ( longcles       )
      REAL :: p(iip1,jjp1,llm)
      INTEGER :: itau, iday
      REAL :: masse(iip1,jjp1,llm)
      REAL :: xpn,xps,xppn(iim),xpps(iim)
      real :: time
      REAL :: phi(ip1jmp1,llm)
      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
      REAL :: w(ip1jmp1,llm)
      REAL ::phystep,co2_ppm,solaire
      INTEGER :: radpas

      CHARACTER*80 :: visu_file
      INTEGER :: visuid
      !
      !   Constantes 
      !
      pi     = 4. * ATAN(1.)
      rad    = 6371229.
      omeg   = 4.* ASIN(1.)/(24.*3600.)
      g      = 9.8
      daysec = 86400.
      kappa  = 0.2857143
      cpp    = 1004.70885
      !
      preff     = 101325.
      unskap = 1./kappa
      !
      jmp1    = jjm + 1
      !
      !    Construct a grid
      !

      CALL defrun_new(99,.TRUE.,clesphy0)

      dtvr   = daysec/FLOAT(day_step)
      print*,'dtvr',dtvr

      CALL inicons0()
      CALL inigeom()
      !
      CALL inifilr()
      !
      latfi(1) = ASIN(1.0)
      DO j = 2, jjm
        DO i = 1, iim
          latfi((j-2)*iim+1+i)=  rlatu(j)
        ENDDO
      ENDDO
      latfi(klon) = - ASIN(1.0)
      !
      lonfi(1) = 0.0
      DO j = 2, jjm
        DO i = 1, iim
          lonfi((j-2)*iim+1+i) =  rlonv(i)
        ENDDO
      ENDDO
      lonfi(klon) = 0.0
      !
      xpi = 2.0 * ASIN(1.0)
      DO ig = 1, klon
        latfi(ig) = latfi(ig) * 180.0 / xpi
        lonfi(ig) = lonfi(ig) * 180.0 / xpi
      ENDDO
      !
      varname = 'relief'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      orog(:,:) = 0.0
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0)
      !
      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
      WRITE(*,'(49I1)') INT(orog(:,:))
      !
      varname = 'rugosite'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      rugo(:,:) = 0.0
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0)
      !
      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 
      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
      !
      varname = 'masque'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      masque(:,:) = 0.0
       CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0)
      !
      WRITE(*,*) 'MASQUE construit : Masque'
      WRITE(*,'(49I1)') INT(masque(:,:))
      !
      !
      varname = 'psol'
      psol(:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0)
      !
      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 
      !  anyway.
      !
      WRITE(*,*) 'PSOL :', psol(10,20)
      WRITE(*,*) ap(:), bp(:)
      CALL pression(ip1jmp1, ap, bp, psol, p3d)
      WRITE(*,*) 'P3D :', p3d(10,20,:)
      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
      WRITE(*,*) 'PK:', pk(10,20,:)
      !
      !
      !
      prefkap =  preff  ** kappa
      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
      DO l = 1, llm
        DO j=1,jjp1
          DO i =1, iip1
            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
           ENDDO
        ENDDO
      ENDDO
      !
      WRITE(*,*) 'PLS :', pls(10,20,:)
      !
      varname = 'surfgeo'
      phis(:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0)
      !
      varname = 'u'
      uvent(:,:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
     . workvar, uvent, 0.0)
      !  
      varname = 'v'
      vvent(:,:,:) = 0.0
      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
     . workvar, vvent, 0.0)
      !
      varname = 't'
      t3d(:,:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
     . workvar, t3d, 0.0)
      !
      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
     .                          maxval(t3d(:,:,:))
      varname = 'tpot'
      tpot(:,:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
     . pk, tpot, 0.0)
      !
      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
     .                          maxval(t3d(:,:,:))
      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
     .                          maxval(pls(:,:,:))
      DO l = 1, llm
        DO j=1,jjp1
          DO i =1, iip1-1
           qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
          ENDDO
          qsat(iip1,j,l) = qsat(1,j,l)
        ENDDO
      ENDDO
      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
     .                           maxval(qsat(:,:,:))
      !
      WRITE(*,*) 'QSAT :', qsat(10,20,:)
      !
      varname = 'q'
      qd(:,:,:) = 0.0
      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
     .                           maxval(qsat(:,:,:))
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
     . qsat, qd, 0.0)
      q3d(:,:,:,1) = qd(:,:,:)
      !
      varname = 'tsol'
      ! This line needs to be replaced by a call to restget to get the values in the restart file
      tsol(:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0)
      !
      WRITE(*,*) 'TSOL construit :'
      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
      !
      varname = 'qsol'
      qsol(:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0)
      !
      varname = 'snow'
      sn(:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0)
      !
      varname = 'rads'
      radsol(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0)
      !
      varname = 'deltat'
      deltat(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0)
      !
      varname = 'rugmer'
      rugmer(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0)
      !
      varname = 'agsno'
      agesno(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0)

      varname = 'zmea'
      zmea(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0)
      varname = 'zstd'
      zstd(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0)
      varname = 'zsig'
      zsig(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0)
      varname = 'zgam'
      zgam(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0)
      varname = 'zthe'
      zthe(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0)
      varname = 'zpic'
      zpic(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0)
      varname = 'zval'
      zval(:) = 0.0
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0)
      rugsrel(:) = 0.0
C Calcul intermediaire
c 
      CALL massdair( p3d, masse  )
c

      print *,' ALPHAX ',alphax

      DO  l = 1, llm
        DO  i    = 1, iim
          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
        ENDDO
          xpn      = SUM(xppn)/apoln
          xps      = SUM(xpps)/apols
        DO i   = 1, iip1
          masse(   i   ,   1     ,  l )   = xpn
          masse(   i   ,   jjp1  ,  l )   = xps
        ENDDO
      ENDDO
      q3d(iip1,:,:,:) = q3d(1,:,:,:)
      phis(iip1,:) = phis(1,:)

C Ecriture


      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
     *                tetagdiv, tetagrot , tetatemp              )
      print*,'sortie inidissip'
      itau = 0
      iday = dayref +itau/day_step
      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
c     
      IF(time.GT.1)  THEN
       time = time - 1
       iday = iday + 1
      ENDIF
      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
      print*,'sortie geopot'
      
      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
     *                phi,w, pbaru,pbarv,time+iday-dayref   )
       print*,'sortie caldyn0'     
      CALL dynredem0("start.nc",dayref,anneeref,phis,nqmx)
      print*,'sortie dynredem0'
      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
     .                            psol)
      print*,'sortie dynredem1' 
C
C Ecriture etat initial physique
C
      phystep   = dtvr * FLOAT(iphysiq)
      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
      co2_ppm   = 330.0
      solaire   = 1370.0

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

C     Sortie Visu pour les champs dynamiques
      print*,'sortie visu'
      time_step = 1.
      t_ops = 2.
      t_wrt = 2.
      itau = 2.
      visu_file='Etat0_visu.nc'
      CALL initdynav(visu_file,dayref,anneeref,time_step,
     .              t_ops, t_wrt, nqmx, visuid)
      CALL writedynav(visuid, nqmx, itau,vvent ,
     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
      print*,'entree histclo'
      CALL histclo
      RETURN
      !
      END SUBROUTINE etat0_netcdf