SUBROUTINE etat0_netcdf
    
      USE startvar
      USE ioipsl
      !
      IMPLICIT NONE
      !
#include "dimensions.h"
#include "paramet.h"
      !
      !
c      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 
c     .KLON=KFDIA-KIDIA+1,KLEV=llm
      !
#include "comgeom2.h"
#include "comvert.h"
#include "comconst.h"
#include "indicesol.h"
#include "dimphy.h"
#include "dimsoil.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)
      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
      REAL :: tsoil(klon,nsoilmx,nbsrf) 
      REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
      REAL :: solsw(klon), sollw(klon)
      REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon),rugmer(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)
      REAL :: pctsrf(klon, nbsrf)
      REAL :: t_ancien(klon,klev), q_ancien(klon,klev)      !
      ! declarations pour lecture glace de mer
      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
      INTEGER :: itaul(1), fid
      REAL :: lev(1), date, dt
      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
      REAL :: flic_tmp(iip1, jjp1)
      REAL :: champint(iim, jjp1)
      !

      CHARACTER*80 :: varname
      !
      INTEGER :: i,j, ig, l, ji
      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(:,:))
      !
      !


C
C on initialise les sous surfaces
C
      pctsrf=0.
      !cree le masque a partir du fichier relief
      varname = 'zmasq'
      zmasq(:) = 0.
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
      WHERE (zmasq(1 : klon) .LE. EPSFRA)
          zmasq(1 : klon) = 0.
      END WHERE 
      WRITE(*,*)zmasq




      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'
      q3d(:,:,:,:) = 0.0
      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
C lecture du fichier glace de terre pour fixer la fraction de terre 
C et de glace de terre
C
      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
     $    , fid)
      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(dlon_lic(iml_lic), stat=iret)
      ALLOCATE(dlat_lic(jml_lic), stat=iret)
      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
     $    , 1, 1, fraclic)
      CALL flinclo(fid)
C
C interpolation sur la grille T du modele
C
      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 
     $    iml_lic, jml_lic
c
C sil les coordonnees sont en degres, on les transforme
C
      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
      ENDIF 
      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 
          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
      ENDIF 

      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 
C
      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
     $    ,iim, jjp1,
     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
C
C passage sur la grille physique
C
      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
     $    pctsrf(1:klon, is_lic))
C adequation avec le maque terre/mer
      WHERE (pctsrf(1 : klon, is_lic) .LE. EPSFRA ) 
          pctsrf(1 : klon, is_lic) = 0. 
      END WHERE
      WHERE (zmasq( 1 : klon) .LE. EPSFRA) 
          pctsrf(1 : klon, is_lic) = 0.
      END WHERE 
      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
      DO ji = 1, klon
        IF (zmasq(ji) .GT. EPSFRA) THEN 
            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
                pctsrf(ji, is_lic) = zmasq(ji)
                pctsrf(ji, is_ter) = 0.
            ELSE 
                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
            ENDIF 
        ENDIF 
      END DO 
C
C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
C
      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
          pctsrf(1 : klon, is_oce) = 0.
      END WHERE 
C
C verif que somme des sous surface = 1
C
      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) 
     $    .GT. EPSFRA)
      IF (ji .NE. 0) THEN
          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
      ENDIF 





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

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

c
c Initialisation 
c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
c
      tsolsrf(:,is_ter) = tsol
      tsolsrf(:,is_lic) = tsol
      tsolsrf(:,is_oce) = tsol
      tsolsrf(:,is_sic) = tsol
      snsrf(:,is_ter) = sn
      snsrf(:,is_lic) = sn
      snsrf(:,is_oce) = sn
      snsrf(:,is_sic) = sn
      albe(:,is_ter) = 0.08
      albe(:,is_lic) = 0.6
      albe(:,is_oce) = 0.5
      albe(:,is_sic) = 0.6
      evap(:,:) = 0.
      qsolsrf(:,is_ter) = qsol
      qsolsrf(:,is_lic) = qsol
      qsolsrf(:,is_oce) = 150.
      qsolsrf(:,is_sic) = 150.
      do i = 1, nbsrf
        do j = 1, nsoilmx
          tsoil(:,j,i) = tsol
        enddo
      enddo
      rain_fall = 0.; snow_fall = 0.
      solsw = 165.
      sollw = -53.
      t_ancien = 273.15
      q_ancien = 0.
      agesno = 0.
      deltat = 0.
      frugs(:,is_oce) = rugmer
      frugs(:,is_ter) = rugmer
      frugs(:,is_lic) = rugmer
      frugs(:,is_sic) = rugmer

      call physdem("startphy.nc",phystep,radpas, co2_ppm, solaire,
     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf, 
     $    albe, evap, rain_fall, snow_fall, solsw, sollw, 
     $    radsol, frugs,  agesno, 
     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, 
     $    t_ancien, q_ancien)

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