PROGRAM modif_etat0
      
      USE startvar
      USE ioipsl
      !
      IMPLICIT NONE
      !
#include "dimensions.h"
#include "paramet.h"
      !
      !
#include "dimphy.h"
      !
#include "comgeom2.h"
#include "comvert.h"
#include "comconst.h"
#include "indicesol.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 :: tsol(klon,nbsrf), qsol(klon,nbsrf), sn(klon,nbsrf)
      REAL :: evap(klon,nbsrf), albe(klon,nbsrf)
      real :: rain_fall(klon), snow_fall(klon)
      real :: sollw(klon), solsw(klon)
      REAL :: radsol(klon),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 :: pctsrf(klon, nbsrf)
      REAL :: tsoil(klon, nsoilmx, nbsrf)
      REAL :: t_ancien(klon, klev), q_ancien(klon, klev)
      REAL :: tabcntrl0(100)
      LOGICAL :: ancien_ok
      ! declarations pour lecture glace de mer
      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
      INTEGER :: itau(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
      !
#include "comdissnew.h"
#include "control.h"
#include "serre.h"
#include "clesph0.h"

      INTEGER  ::        longcles
      PARAMETER      ( longcles  = 20 )
      REAL :: clesphy0 ( longcles       )
      REAL ::phystep,co2_ppm,solaire
      REAL :: unskap
      INTEGER :: radpas

      CHARACTER*80 :: visu_file
      INTEGER :: visuid
      !
      ! indices de boucle 
      INTEGER :: ji
      !   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()
      !
      !
      !! 1. READ orography 
      !
      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(*,*) masque(:,:)
      WRITE(*,*) 'orographie'
      WRITE(*,*) orog(:,:)
      !
      !! 2. create masque and READ land iice fraction
      !!
      ! lit bande startphy a modifier
      !
     
      CALL phyetat0("startphy.nc", phystep, co2_ppm, solaire, latfi,  
     $    lonfi, pctsrf, tsol, tsoil, deltat, qsol, sn, 
     $    albe, evap, rain_fall, snow_fall, solsw, sollw,
     $    radsol, rugmer, agesno, clesphy0, 
     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, 
     $    rugsrel, tabcntrl0, t_ancien, q_ancien, ancien_ok)
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 = '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, itau, 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
C Ecriture etat initial physique
C
      WRITE(*,*) 'zmasq avant phyredem'
      WRITE(*,*) zmasq

      call phyredem("restartphy.nc",phystep,radpas, co2_ppm, solaire,
     $    latfi, lonfi, pctsrf, tsol, tsoil, deltat, qsol, sn, 
     $    albe, evap, rain_fall, snow_fall, solsw, sollw, 
     $    radsol, rugmer,  agesno, 
     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, 
     $    t_ancien, q_ancien)

      CALL histclo
      !
      END