! $Id: phyetat0.F90 2243 2015-03-24 13:28:51Z fhourdin $

SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)

  USE dimphy, only: klon, zmasq, klev, nslay
  USE iophy, ONLY : init_iophy_new
  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
  USE fonte_neige_mod,  ONLY : fonte_neige_init
  USE pbl_surface_mod,  ONLY : pbl_surface_init
  USE surface_data,     ONLY : type_ocean, version_ocean
  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
       qsol, fevap, z0m, z0h, agesno, &
       du_gwd_rando, dv_gwd_rando, entr_therm, f0, fm_therm, &
       falb_dir, falb_dif, &
       ftsol, pbl_tke, pctsrf, q_ancien, radpas, radsol, rain_fall, ratqs, &
       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
       wake_s, zgam, &
       zmax0, zmea, zpic, zsig, &
       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl
  USE iostart, ONLY : close_startphy, get_field, get_var, open_startphy
  USE infotrac, only: nbtr, type_trac, tname, niadv
  USE traclmdz_mod,    ONLY : traclmdz_from_restart
  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init

  IMPLICIT none
  !======================================================================
  ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
  ! Objet: Lecture de l'etat initial pour la physique
  !======================================================================
  include "dimensions.h"
  include "netcdf.inc"
  include "dimsoil.h"
  include "clesphys.h"
  include "temps.h"
  include "thermcell.h"
  include "compbl.h"
  include "YOMCST.h"
  !======================================================================
  CHARACTER*(*) fichnom

  ! les variables globales lues dans le fichier restart

  REAL tsoil(klon, nsoilmx, nbsrf)
  REAL qsurf(klon, nbsrf)
  REAL snow(klon, nbsrf)
  real fder(klon)
  REAL run_off_lic_0(klon)
  REAL fractint(klon)
  REAL trs(klon, nbtr)
  REAL zts(klon)

  CHARACTER*6 ocean_in
  LOGICAL ok_veget_in

  INTEGER        longcles
  PARAMETER    ( longcles = 20 )
  REAL clesphy0( longcles )

  REAL xmin, xmax

  INTEGER nid, nvarid
  INTEGER ierr, i, nsrf, isoil , k
  INTEGER length
  PARAMETER (length=100)
  INTEGER it, iiq, isw
  REAL tab_cntrl(length), tabcntr0(length)
  CHARACTER*7 str7
  CHARACTER*2 str2
  LOGICAL :: found,phyetat0_get,phyetat0_srf

  ! FH1D
  !     real iolat(jjm+1)
  real iolat(jjm+1-1/(iim*jjm))

  ! Ouvrir le fichier contenant l'etat initial:

  CALL open_startphy(fichnom)

  ! Lecture des parametres de controle:

  CALL get_var("controle", tab_cntrl)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
  ! Les constantes de la physiques sont lues dans la physique seulement.
  ! Les egalites du type
  !             tab_cntrl( 5 )=clesphy0(1)
  ! sont remplacees par
  !             clesphy0(1)=tab_cntrl( 5 )
  ! On inverse aussi la logique.
  ! On remplit les tab_cntrl avec les parametres lus dans les .def
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  DO i = 1, length
     tabcntr0( i ) = tab_cntrl( i )
  ENDDO

  tab_cntrl(1)=dtime
  tab_cntrl(2)=radpas

  ! co2_ppm : value from the previous time step
  IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
     co2_ppm = tab_cntrl(3)
     RCO2    = co2_ppm * 1.0e-06  * 44.011/28.97 
     ! ELSE : keep value from .def
  END IF

  ! co2_ppm0 : initial value of atmospheric CO2 (from create_etat0_limit.e .def)
  co2_ppm0   = tab_cntrl(16)

  solaire_etat0      = tab_cntrl(4)
  tab_cntrl(5)=iflag_con
  tab_cntrl(6)=nbapp_rad

  if (cycle_diurne) tab_cntrl( 7) =1.
  if (soil_model) tab_cntrl( 8) =1.
  if (new_oliq) tab_cntrl( 9) =1.
  if (ok_orodr) tab_cntrl(10) =1.
  if (ok_orolf) tab_cntrl(11) =1.
  if (ok_limitvrai) tab_cntrl(12) =1.

  itau_phy = tab_cntrl(15)

  clesphy0(1)=tab_cntrl( 5 )
  clesphy0(2)=tab_cntrl( 6 )
  clesphy0(3)=tab_cntrl( 7 )
  clesphy0(4)=tab_cntrl( 8 )
  clesphy0(5)=tab_cntrl( 9 )
  clesphy0(6)=tab_cntrl( 10 )
  clesphy0(7)=tab_cntrl( 11 )
  clesphy0(8)=tab_cntrl( 12 )

  ! Lecture des latitudes (coordonnees):

  CALL get_field("latitude", rlat)

  ! Lecture des longitudes (coordonnees):

  CALL get_field("longitude", rlon)

  ! Lecture du masque terre mer

  CALL get_field("masque", zmasq, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <masque> est absent'
     PRINT *, 'fichier startphy non compatible avec phyetat0'
  ENDIF

  ! Lecture des fractions pour chaque sous-surface

  ! initialisation des sous-surfaces

  pctsrf = 0.

  ! fraction de terre

  CALL get_field("FTER", pctsrf(:, is_ter), found)
  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FTER> est absent'

  ! fraction de glace de terre

  CALL get_field("FLIC", pctsrf(:, is_lic), found)
  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FLIC> est absent'

  ! fraction d'ocean

  CALL get_field("FOCE", pctsrf(:, is_oce), found)
  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FOCE> est absent'

  ! fraction glace de mer

  CALL get_field("FSIC", pctsrf(:, is_sic), found)
  IF (.NOT. found) PRINT*, 'phyetat0: Le champ <FSIC> est absent'

  !  Verification de l'adequation entre le masque et les sous-surfaces

  fractint( 1 : klon) = pctsrf(1 : klon, is_ter)  &
       + pctsrf(1 : klon, is_lic)
  DO i = 1 , klon
     IF ( abs(fractint(i) - zmasq(i) ) .GT. EPSFRA ) THEN
        WRITE(*, *) 'phyetat0: attention fraction terre pas ',  &
             'coherente ', i, zmasq(i), pctsrf(i, is_ter) &
             , pctsrf(i, is_lic)
        WRITE(*, *) 'Je force la coherence zmasq=fractint'
        zmasq(i) = fractint(i)
     ENDIF
  END DO
  fractint (1 : klon) =  pctsrf(1 : klon, is_oce)  &
       + pctsrf(1 : klon, is_sic)
  DO i = 1 , klon
     IF ( abs( fractint(i) - (1. - zmasq(i))) .GT. EPSFRA ) THEN
        WRITE(*, *) 'phyetat0 attention fraction ocean pas ',  &
             'coherente ', i, zmasq(i) , pctsrf(i, is_oce) &
             , pctsrf(i, is_sic)
        WRITE(*, *) 'Je force la coherence zmasq=1.-fractint'
        zmasq(i) = 1. - fractint(i)
     ENDIF
  END DO

  ! Lecture des temperatures du sol:

  CALL get_field("TS", ftsol(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <TS> est absent'
     PRINT*, '          Mais je vais essayer de lire TS**'
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field("TS"//str2, ftsol(:, nsrf))

        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(ftsol(i, nsrf), xmin)
           xmax = MAX(ftsol(i, nsrf), xmax)
        ENDDO
        PRINT*, 'Temperature du sol TS**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <TS> est present'
     PRINT*, '          J ignore donc les autres temperatures TS**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(ftsol(i, 1), xmin)
        xmax = MAX(ftsol(i, 1), xmax)
     ENDDO
     PRINT*, 'Temperature du sol <TS>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           ftsol(i, nsrf) = ftsol(i, 1)
        ENDDO
     ENDDO
  ENDIF

!===================================================================
  ! Lecture des albedo difus et direct

  DO nsrf = 1, nbsrf
     DO isw=1, nsw
        IF (isw.GT.99 .AND. nsrf.GT.99) THEN
           PRINT*, "Trop de bandes SW ou sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str7, '(i2.2, "srf", i2.2)') isw, nsrf

        CALL get_field('A_dir_SW'//str7, falb_dir(:, isw, nsrf), found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: Le champ <A_dir_SW"//str7//"> est absent"
           PRINT*, "          Il prend donc la valeur de surface"
           DO i=1, klon
              falb_dir(i, isw, nsrf)=0.2
           ENDDO
        ENDIF
        CALL get_field('A_dif_SW'//str7, falb_dif(:, isw, nsrf), found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: Le champ <A_dif_SW"//str7//"> est absent"
           PRINT*, "          Il prend donc la valeur de surface"
           DO i=1, klon
              falb_dif(i, isw, nsrf)=0.2
           ENDDO
        ENDIF
     ENDDO
  ENDDO

!===================================================================
  ! Lecture des temperatures du sol profond:

  DO nsrf = 1, nbsrf
     DO isoil=1, nsoilmx
        IF (isoil.GT.99 .AND. nsrf.GT.99) THEN
           PRINT*, "Trop de couches ou sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str7, '(i2.2, "srf", i2.2)') isoil, nsrf

        CALL get_field('Tsoil'//str7, tsoil(:, isoil, nsrf), found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
           PRINT*, "          Il prend donc la valeur de surface"
           DO i=1, klon
              tsoil(i, isoil, nsrf)=ftsol(i, nsrf)
           ENDDO
        ENDIF
     ENDDO
  ENDDO

!===================================================================
  ! Lecture de l'humidite de l'air juste au dessus du sol:

  CALL get_field("QS", qsurf(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <QS> est absent'
     PRINT*, '          Mais je vais essayer de lire QS**'
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field("QS"//str2, qsurf(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(qsurf(i, nsrf), xmin)
           xmax = MAX(qsurf(i, nsrf), xmax)
        ENDDO
        PRINT*, 'Humidite pres du sol QS**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <QS> est present'
     PRINT*, '          J ignore donc les autres humidites QS**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(qsurf(i, 1), xmin)
        xmax = MAX(qsurf(i, 1), xmax)
     ENDDO
     PRINT*, 'Humidite pres du sol <QS>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           qsurf(i, nsrf) = qsurf(i, 1)
        ENDDO
     ENDDO
  ENDIF

  ! Eau dans le sol (pour le modele de sol "bucket")

  CALL get_field("QSOL", qsol, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <QSOL> est absent'
     PRINT*, '          Valeur par defaut nulle'
     qsol(:)=0.
  ENDIF

  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(qsol(i), xmin)
     xmax = MAX(qsol(i), xmax)
  ENDDO
  PRINT*, 'Eau dans le sol (mm) <QSOL>', xmin, xmax

  ! Lecture de neige au sol:

  CALL get_field("SNOW", snow(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <SNOW> est absent'
     PRINT*, '          Mais je vais essayer de lire SNOW**'
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field( "SNOW"//str2, snow(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(snow(i, nsrf), xmin)
           xmax = MAX(snow(i, nsrf), xmax)
        ENDDO
        PRINT*, 'Neige du sol SNOW**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <SNOW> est present'
     PRINT*, '          J ignore donc les autres neiges SNOW**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(snow(i, 1), xmin)
        xmax = MAX(snow(i, 1), xmax)
     ENDDO
     PRINT*, 'Neige du sol <SNOW>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           snow(i, nsrf) = snow(i, 1)
        ENDDO
     ENDDO
  ENDIF

  ! Lecture de evaporation:  

  CALL get_field("EVAP", fevap(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <EVAP> est absent'
     PRINT*, '          Mais je vais essayer de lire EVAP**'
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field("EVAP"//str2, fevap(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(fevap(i, nsrf), xmin)
           xmax = MAX(fevap(i, nsrf), xmax)
        ENDDO
        PRINT*, 'fevap du sol EVAP**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <EVAP> est present'
     PRINT*, '          J ignore donc les autres EVAP**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(fevap(i, 1), xmin)
        xmax = MAX(fevap(i, 1), xmax)
     ENDDO
     PRINT*, 'Evap du sol <EVAP>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           fevap(i, nsrf) = fevap(i, 1)
        ENDDO
     ENDDO
  ENDIF

  ! Lecture precipitation liquide:

  CALL get_field("rain_f", rain_fall)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(rain_fall(i), xmin)
     xmax = MAX(rain_fall(i), xmax)
  ENDDO
  PRINT*, 'Precipitation liquide rain_f:', xmin, xmax

  ! Lecture precipitation solide:

  CALL get_field("snow_f", snow_fall)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(snow_fall(i), xmin)
     xmax = MAX(snow_fall(i), xmax)
  ENDDO
  PRINT*, 'Precipitation solide snow_f:', xmin, xmax

  ! Lecture rayonnement solaire au sol:

  CALL get_field("solsw", solsw, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <solsw> est absent'
     PRINT*, 'mis a zero'
     solsw(:) = 0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(solsw(i), xmin)
     xmax = MAX(solsw(i), xmax)
  ENDDO
  PRINT*, 'Rayonnement solaire au sol solsw:', xmin, xmax

  ! Lecture rayonnement IF au sol:

  CALL get_field("sollw", sollw, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <sollw> est absent'
     PRINT*, 'mis a zero'
     sollw = 0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(sollw(i), xmin)
     xmax = MAX(sollw(i), xmax)
  ENDDO
  PRINT*, 'Rayonnement IF au sol sollw:', xmin, xmax

  CALL get_field("sollwdown", sollwdown, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
     PRINT*, 'mis a zero'
     sollwdown = 0.
     zts=0.
     do nsrf=1,nbsrf
        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
     enddo
     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
  ENDIF
!  print*,'TS SOLL',zts(klon/2),sollw(klon/2),sollwdown(klon/2)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(sollwdown(i), xmin)
     xmax = MAX(sollwdown(i), xmax)
  ENDDO
  PRINT*, 'Rayonnement IF au sol sollwdown:', xmin, xmax


  ! Lecture derive des flux:

  CALL get_field("fder", fder, found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <fder> est absent'
     PRINT*, 'mis a zero'
     fder = 0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(fder(i), xmin)
     xmax = MAX(fder(i), xmax)
  ENDDO
  PRINT*, 'Derive des flux fder:', xmin, xmax

  ! Lecture du rayonnement net au sol:

  CALL get_field("RADS", radsol)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(radsol(i), xmin)
     xmax = MAX(radsol(i), xmax)
  ENDDO
  PRINT*, 'Rayonnement net au sol radsol:', xmin, xmax

  ! Lecture de la longueur de rugosite 

IF (1==0) THEN ! A DERTRUIRE TOUT DE SUITE
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
! Retrocompatibilite. A nettoyer fin 2015
        CALL get_field("RUG"//str2, z0m(:, nsrf),found)
        IF (found) THEN
            z0h(:,nsrf)=z0m(:,nsrf)
            PRINT*,'Lecture de ',"RUG"//str2,' -> z0m/z0h (obsolete)'
        ELSE
            CALL get_field("Z0m"//str2, z0m(:, nsrf), found)
            IF (.NOT.found) Z0m=1.e-3 ! initialisation à 1mm au cas ou.
            CALL get_field("Z0h"//str2, z0h(:, nsrf), found)
            IF (.NOT.found) Z0h=1.e-3 ! initialisation à 1mm au cas ou.
        ENDIF
        PRINT*, 'rugosite Z0m',nsrf,minval(z0m(:, nsrf)),maxval(z0m(:, nsrf))
        PRINT*, 'rugosite Z0h',nsrf,minval(z0h(:, nsrf)),maxval(z0h(:, nsrf))

     ENDDO
ELSE
  PRINT*,'AVANT phyetat0_srf'
  found=phyetat0_srf(1,z0m,"RUG","Z0m ancien",0.001)
  PRINT*,'APRES phyetat0_srf'
  IF (found) THEN
     z0h(:,1:nbsrf)=z0m(:,1:nbsrf)
  ELSE
     found=phyetat0_srf(1,z0m,"Z0m","Roughness length, momentum ",0.001)
     found=phyetat0_srf(1,z0h,"Z0h","Roughness length, enthalpy ",0.001)
  ENDIF
ENDIF

  ! Lecture de l'age de la neige:

  CALL get_field("AGESNO", agesno(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <AGESNO> est absent'
     PRINT*, '          Mais je vais essayer de lire AGESNO**'
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field("AGESNO"//str2, agesno(:, nsrf), found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: Le champ <AGESNO"//str2//"> est absent"
           agesno = 50.0
        ENDIF
        PRINT*, 'agesno',nsrf,minval(agesno(:, nsrf)),maxval(agesno(:, nsrf))
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <AGESNO> est present'
     PRINT*, '          J ignore donc les autres AGESNO**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(agesno(i, 1), xmin)
        xmax = MAX(agesno(i, 1), xmax)
     ENDDO
     PRINT*, 'Age de la neige <AGESNO>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           agesno(i, nsrf) = agesno(i, 1)
        ENDDO
     ENDDO
  ENDIF

!  CALL get_field("ZMEA", zmea)
!  PRINT*, 'OROGRAPHIE SOUS-MAILLE zmea:',minval(zmea(:)),maxval(zmea(:))
  found=phyetat0_get(1,zmea,"ZMEA","mean orography",0.)

  CALL get_field("ZSTD", zstd)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zstd:',minval(zstd(:)),maxval(zstd(:))

  CALL get_field("ZSIG", zsig)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zsig:',minval(zsig(:)),maxval(zsig(:))

  CALL get_field("ZGAM", zgam)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zgam:',minval(zgam(:)),maxval(zgam(:))

  CALL get_field("ZTHE", zthe)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zthe:',minval(zthe(:)),maxval(zthe(:))

  CALL get_field("ZPIC", zpic)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zpic:',minval(zpic(:)),maxval(zpic(:))

  CALL get_field("ZVAL", zval)
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zval:',minval(zval(:)),maxval(zval(:))

  CALL get_field("RUGSREL", rugoro)
  PRINT*, 'Rugosite relief (ecart-type) rugsrel:',minval(rugoro(:)),maxval(rugoro(:))

  ancien_ok = .TRUE.

  CALL get_field("TANCIEN", t_ancien, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <TANCIEN> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ancien_ok = .FALSE.
  ENDIF

  CALL get_field("QANCIEN", q_ancien, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <QANCIEN> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ancien_ok = .FALSE.
  ENDIF

  CALL get_field("UANCIEN", u_ancien, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <UANCIEN> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ancien_ok = .FALSE.
  ENDIF

  CALL get_field("VANCIEN", v_ancien, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <VANCIEN> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ancien_ok = .FALSE.
  ENDIF

  clwcon=0.
  CALL get_field("CLWCON", clwcon, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ CLWCON est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
  ENDIF
  PRINT*,'Eau liquide convective (ecart-type) clwcon:',MINval(clwcon),MAXval(clwcon)


  rnebcon = 0.
  CALL get_field("RNEBCON", rnebcon, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ RNEBCON est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
  ENDIF
  PRINT*, 'Nebulosite convective (ecart-type) rnebcon:',MINval(rnebcon),MAXval(rnebcon)

  ! Lecture ratqs

  ratqs=0.
  CALL get_field("RATQS", ratqs, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <RATQS> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
  ENDIF
  PRINT*, '(ecart-type) ratqs:', MINval(ratqs),MAXval(ratqs)

  ! Lecture run_off_lic_0

  CALL get_field("RUNOFFLIC0", run_off_lic_0, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <RUNOFFLIC0> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     run_off_lic_0 = 0.
  ENDIF
  PRINT*, '(ecart-type) run_off_lic_0:', MINval(run_off_lic_0),MAXval(run_off_lic_0)

  ! Lecture de l'energie cinetique turbulente

  IF (iflag_pbl>1) then
     DO nsrf = 1, nbsrf
        IF (nsrf.GT.99) THEN
           PRINT*, "Trop de sous-mailles"
           call abort_gcm("phyetat0", "", 1)
        ENDIF
        WRITE(str2, '(i2.2)') nsrf
        CALL get_field("TKE"//str2, pbl_tke(:, 1:klev+1, nsrf), found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: <TKE"//str2//"> est absent"
           pbl_tke(:, :, nsrf)=1.e-8
        ENDIF
        PRINT*, 'Turbulent kinetic energyl TKE**:', nsrf, minval(pbl_tke(:,:,nsrf)),maxval(pbl_tke(:,:, nsrf))

     ENDDO
  ENDIF

! Lecture de l'ecart de TKE (w) - (x)
!
  IF (iflag_pbl>1 .AND. iflag_wake>=1  &
           .AND. iflag_pbl_split >=1 ) then
    DO nsrf = 1, nbsrf
      IF (nsrf.GT.99) THEN
        PRINT*, "Trop de sous-mailles"
        call abort_gcm("phyetat0", "", 1)
      ENDIF
      WRITE(str2,'(i2.2)') nsrf
      CALL get_field("DELTATKE"//str2, &
                    wake_delta_pbl_tke(:,1:klev+1,nsrf),found)
      IF (.NOT. found) THEN
        PRINT*, "phyetat0: <DELTATKE"//str2//"> est absent"
        wake_delta_pbl_tke(:,:,nsrf)=0.
      ENDIF
      PRINT*,'TKE difference (w)-(x) DELTATKE**:', nsrf,       &
      minval(wake_delta_pbl_tke(:,:,nsrf)),maxval(wake_delta_pbl_tke(:,:, nsrf))

    ENDDO

  ! delta_tsurf

    DO nsrf = 1, nbsrf
       IF (nsrf.GT.99) THEN
         PRINT*, "Trop de sous-mailles"
         call abort_gcm("phyetat0", "", 1)
       ENDIF
       WRITE(str2,'(i2.2)') nsrf
     CALL get_field("DELTA_TSURF"//str2, delta_tsurf(:,nsrf), found)
     IF (.NOT. found) THEN
        PRINT*, "phyetat0: Le champ <DELTA_TSURF"//str2//"> est absent"
        PRINT*, "Depart legerement fausse. Mais je continue"
        delta_tsurf(:,nsrf)=0.
     ELSE
        PRINT*, 'delta_tsurf:', nsrf,       &
      minval(delta_tsurf(:,nsrf)),maxval(delta_tsurf(:, nsrf))
     ENDIF
    ENDDO  ! nsrf = 1, nbsrf
  ENDIF   !(iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1 )

  ! zmax0
  CALL get_field("ZMAX0", zmax0, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <ZMAX0> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     zmax0=40.
  ENDIF
  PRINT*, '(ecart-type) zmax0:', MINval(zmax0),MAXval(zmax0)

  !           f0(ig)=1.e-5
  ! f0
  CALL get_field("F0", f0, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <f0> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     f0=1.e-5
  ENDIF
  PRINT*, '(ecart-type) f0:', MINval(f0),MAXval(f0)

  ! sig1 or ema_work1

  CALL get_field("sig1", sig1, found)
  IF (.NOT. found) CALL get_field("EMA_WORK1", sig1, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ sig1 est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     sig1=0.
  ELSE
     PRINT*, 'sig1:',minval(sig1(:,:)),maxval(sig1(:,:))
  ENDIF

  ! w01 or ema_work2

  CALL get_field("w01", w01, found)
  IF (.NOT. found) CALL get_field("EMA_WORK2", w01, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ w01 est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     w01=0.
  ELSE
     PRINT*, 'w01:', minval(w01(:,:)),maxval(w01(:,:))
  ENDIF

  ! wake_deltat

  CALL get_field("WAKE_DELTAT", wake_deltat, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <WAKE_DELTAT> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     wake_deltat=0.
  ELSE
     PRINT*, 'wake_deltat:',  minval(wake_deltat(:,:)),maxval(wake_deltat(:,:))
  ENDIF

  ! wake_deltaq

  found=phyetat0_get(klev,wake_deltaq,"WAKE_DELTAQ","Delta hum. wake/env",0.)
! CALL get_field("WAKE_DELTAQ", wake_deltaq, found)
! IF (.NOT. found) THEN
!    PRINT*, "phyetat0: Le champ <WAKE_DELTAQ> est absent"
!    PRINT*, "Depart legerement fausse. Mais je continue"
!    wake_deltaq=0.
! ELSE
!    PRINT*, 'wake_deltaq:',  minval(wake_deltaq(:,:)),maxval(wake_deltaq(:,:))
! ENDIF

  ! wake_s

  CALL get_field("WAKE_S", wake_s, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <WAKE_S> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     wake_s=0.
  ENDIF
  PRINT*, '(ecart-type) wake_s:', MINval(wake_s),MAXval(wake_s)

  ! wake_cstar

  CALL get_field("WAKE_CSTAR", wake_cstar, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <WAKE_CSTAR> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     wake_cstar=0.
  ENDIF
  PRINT*, '(ecart-type) wake_cstar:', MINval(wake_cstar),MAXval(wake_cstar)

  ! wake_pe

  CALL get_field("WAKE_PE", wake_pe, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <WAKE_PE> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     wake_pe=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(wake_pe)
  xmax = MAXval(wake_pe)
  PRINT*, '(ecart-type) wake_pe:', xmin, xmax

  ! wake_fip

  CALL get_field("WAKE_FIP", wake_fip, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <WAKE_FIP> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     wake_fip=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(wake_fip)
  xmax = MAXval(wake_fip)
  PRINT*, '(ecart-type) wake_fip:', xmin, xmax

  !  thermiques

  CALL get_field("FM_THERM", fm_therm, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <fm_therm> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     fm_therm=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(fm_therm)
  xmax = MAXval(fm_therm)
  PRINT*, '(ecart-type) fm_therm:', xmin, xmax

  CALL get_field("ENTR_THERM", entr_therm, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <entr_therm> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     entr_therm=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(entr_therm)
  xmax = MAXval(entr_therm)
  PRINT*, '(ecart-type) entr_therm:', xmin, xmax

  CALL get_field("DETR_THERM", detr_therm, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <detr_therm> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     detr_therm=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(detr_therm)
  xmax = MAXval(detr_therm)
  PRINT*, '(ecart-type) detr_therm:', xmin, xmax

  CALL get_field("ALE_BL", ale_bl, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <ale_bl> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ale_bl=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(ale_bl)
  xmax = MAXval(ale_bl)
  PRINT*, '(ecart-type) ale_bl:', xmin, xmax

  CALL get_field("ALE_BL_TRIG", ale_bl_trig, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <ale_bl_trig> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     ale_bl_trig=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(ale_bl_trig)
  xmax = MAXval(ale_bl_trig)
  PRINT*, '(ecart-type) ale_bl_trig:', xmin, xmax

  CALL get_field("ALP_BL", alp_bl, found)
  IF (.NOT. found) THEN
     PRINT*, "phyetat0: Le champ <alp_bl> est absent"
     PRINT*, "Depart legerement fausse. Mais je continue"
     alp_bl=0.
  ENDIF
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(alp_bl)
  xmax = MAXval(alp_bl)
  PRINT*, '(ecart-type) alp_bl:', xmin, xmax

  ! Read and send field trs to traclmdz

  IF (type_trac == 'lmdz') THEN
     DO it=1, nbtr
        iiq=niadv(it+2)
        CALL get_field("trs_"//tname(iiq), trs(:, it), found)
        IF (.NOT. found) THEN
           PRINT*,  &
                "phyetat0: Le champ <trs_"//tname(iiq)//"> est absent"
           PRINT*, "Depart legerement fausse. Mais je continue"
           trs(:, it) = 0.
        ENDIF
        xmin = 1.0E+20
        xmax = -1.0E+20
        xmin = MINval(trs(:, it))
        xmax = MAXval(trs(:, it))
        PRINT*, "(ecart-type) trs_"//tname(iiq)//" :", xmin, xmax

     END DO
     CALL traclmdz_from_restart(trs)

     IF (carbon_cycle_cpl) THEN
        ALLOCATE(co2_send(klon), stat=ierr)
        IF (ierr /= 0) CALL abort_gcm &
             ('phyetat0', 'pb allocation co2_send', 1)
        CALL get_field("co2_send", co2_send, found)
        IF (.NOT. found) THEN
           PRINT*, "phyetat0: Le champ <co2_send> est absent"
           PRINT*, "Initialisation uniforme a co2_ppm=", co2_ppm
           co2_send(:) = co2_ppm
        END IF
     END IF
  END IF

  if (ok_gwd_rando) then
     call get_field("du_gwd_rando", du_gwd_rando, found)
     if (.not. found) then
        print *, "du_gwd_rando not found, setting it to 0."
        du_gwd_rando = 0.
     end if

     call get_field("dv_gwd_rando", dv_gwd_rando, found)
     if (.not. found) then
        print *, "dv_gwd_rando not found, setting it to 0."
        dv_gwd_rando = 0.
     end if
  end if

  ! Initialize Slab variables
  IF ( type_ocean == 'slab' ) THEN
      print*, "calling slab_init"
      CALL ocean_slab_init(dtime, pctsrf)
      ! tslab
      CALL get_field("tslab", tslab, found)
      IF (.NOT. found) THEN 
          PRINT*, "phyetat0: Le champ <tslab> est absent"
          PRINT*, "Initialisation a tsol_oce"
          DO i=1,nslay
              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
          END DO
      END IF 
      ! Sea ice variables
      IF (version_ocean == 'sicINT') THEN
          CALL get_field("slab_tice", tice, found)
          IF (.NOT. found) THEN 
              PRINT*, "phyetat0: Le champ <tice> est absent"
              PRINT*, "Initialisation a tsol_sic"
                  tice(:)=ftsol(:,is_sic)
          END IF 
          CALL get_field("seaice", seaice, found)
          IF (.NOT. found) THEN
              PRINT*, "phyetat0: Le champ <seaice> est absent"
              PRINT*, "Initialisation a 0/1m suivant fraction glace"
              seaice(:)=0.
              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
                  seaice=917.
              END WHERE
          END IF
      END IF !sea ice INT
  END IF ! Slab        

  ! on ferme le fichier
  CALL close_startphy

  ! Initialize module pbl_surface_mod 

  CALL pbl_surface_init(fder, snow, qsurf, tsoil)

  ! Initialize module ocean_cpl_mod for the case of coupled ocean
  IF ( type_ocean == 'couple' ) THEN
     CALL ocean_cpl_init(dtime, rlon, rlat)
  ENDIF

  CALL init_iophy_new(rlat, rlon)

  ! Initilialize module fonte_neige_mod      
  CALL fonte_neige_init(run_off_lic_0)

END SUBROUTINE phyetat0

!===================================================================
FUNCTION phyetat0_get(nlev,field,name,descr,default)
!===================================================================
! Lecture d'un champ avec contrôle
! Function logique dont le resultat indique si la lecture
! s'est bien passée
! On donne une valeur par defaut dans le cas contraire
!===================================================================

USE iostart, ONLY : get_field
USE dimphy, only: klon 

IMPLICIT NONE
INCLUDE "iniprint.h"

LOGICAL phyetat0_get

! arguments
INTEGER,INTENT(IN) :: nlev
CHARACTER*(*),INTENT(IN) :: name,descr
REAL,INTENT(IN) :: default
REAL,DIMENSION(klon,nlev),INTENT(INOUT) :: field

! Local variables
LOGICAL found

   CALL get_field(name, field, found)
   IF (.NOT. found) THEN
     WRITE(lunout,*) "phyetat0: Le champ <",name,"> est absent"
     WRITE(lunout,*) "Depart legerement fausse. Mais je continue"
     field(:,:)=default
   ENDIF
   WRITE(lunout,*) name, descr, MINval(field),MAXval(field)
   phyetat0_get=found

RETURN
END FUNCTION phyetat0_get

!================================================================
FUNCTION phyetat0_srf(nlev,field,name,descr,default)
!===================================================================
! Lecture d'un champ par sous-surface avec contrôle
! Function logique dont le resultat indique si la lecture
! s'est bien passée
! On donne une valeur par defaut dans le cas contraire
!===================================================================

USE iostart, ONLY : get_field
USE dimphy, only: klon 
USE indice_sol_mod, only: nbsrf

IMPLICIT NONE
INCLUDE "iniprint.h"

LOGICAL phyetat0_srf
! arguments
INTEGER,INTENT(IN) :: nlev
CHARACTER*(*),INTENT(IN) :: name,descr
REAL,INTENT(IN) :: default
REAL,DIMENSION(klon,nlev,nbsrf),INTENT(INOUT) :: field

! Local variables
LOGICAL found,phyetat0_get
INTEGER nsrf
CHARACTER*2 str2
  
     IF (nbsrf.GT.99) THEN
        WRITE(lunout,*) "Trop de sous-mailles"
        call abort_gcm("phyetat0", "", 1)
     ENDIF

     DO nsrf = 1, nbsrf
        WRITE(str2, '(i2.2)') nsrf
        found= phyetat0_get(nlev,field(:,:, nsrf), &
        name//str2,descr//" srf:"//str2,default)
     ENDDO

     phyetat0_srf=found

RETURN
END FUNCTION phyetat0_srf