Index: /LMDZ.3.3/trunk/libf/dyn3d/modif_etat0.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/modif_etat0.F	(revision 96)
+++ /LMDZ.3.3/trunk/libf/dyn3d/modif_etat0.F	(revision 96)
@@ -0,0 +1,259 @@
+      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