source: LMDZ5/trunk/libf/phylmd/phyetat0.F90 @ 2188

! $Id: phyetat0.F90 2188 2015-01-30 16:25:12Z 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
  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
       du_gwd_rando, dv_gwd_rando, entr_therm, f0, falb1, falb2, fm_therm, &
       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, 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 qsol(klon)
  REAL snow(klon, nbsrf)
  REAL evap(klon, nbsrf)
  real fder(klon)
  REAL frugs(klon, nbsrf)
  REAL agesno(klon, nbsrf)
  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
  REAL tab_cntrl(length), tabcntr0(length)
  CHARACTER*7 str7
  CHARACTER*2 str2
  LOGICAL :: found

  ! 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 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 albedo de l'interval visible au sol:

  CALL get_field("ALBE", falb1(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <ALBE> est absent'
     PRINT*, '          Mais je vais essayer de lire ALBE**'
     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("ALBE"//str2, falb1(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(falb1(i, nsrf), xmin)
           xmax = MAX(falb1(i, nsrf), xmax)
        ENDDO
        PRINT*, 'Albedo du sol ALBE**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <ALBE> est present'
     PRINT*, '          J ignore donc les autres ALBE**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(falb1(i, 1), xmin)
        xmax = MAX(falb1(i, 1), xmax)
     ENDDO
     PRINT*, 'Neige du sol <ALBE>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           falb1(i, nsrf) = falb1(i, 1)
        ENDDO
     ENDDO
  ENDIF

  ! Lecture de albedo au sol dans l'interval proche infra-rouge:

  CALL get_field("ALBLW", falb2(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <ALBLW> est absent'
     PRINT*, '          Mais je vais prendre ALBE**'
     DO nsrf = 1, nbsrf
        DO i = 1, klon
           falb2(i, nsrf) = falb1(i, nsrf)
        ENDDO
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <ALBLW> est present'
     PRINT*, '          J ignore donc les autres ALBLW**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(falb2(i, 1), xmin)
        xmax = MAX(falb2(i, 1), xmax)
     ENDDO
     PRINT*, 'Neige du sol <ALBLW>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           falb2(i, nsrf) = falb2(i, 1)
        ENDDO
     ENDDO
  ENDIF

  ! Lecture de evaporation:  

  CALL get_field("EVAP", evap(:, 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, evap(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(evap(i, nsrf), xmin)
           xmax = MAX(evap(i, nsrf), xmax)
        ENDDO
        PRINT*, 'evap 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(evap(i, 1), xmin)
        xmax = MAX(evap(i, 1), xmax)
     ENDDO
     PRINT*, 'Evap du sol <EVAP>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           evap(i, nsrf) = evap(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.
     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 

  CALL get_field("RUG", frugs(:, 1), found)
  IF (.NOT. found) THEN
     PRINT*, 'phyetat0: Le champ <RUG> est absent'
     PRINT*, '          Mais je vais essayer de lire RUG**'
     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("RUG"//str2, frugs(:, nsrf))
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(frugs(i, nsrf), xmin)
           xmax = MAX(frugs(i, nsrf), xmax)
        ENDDO
        PRINT*, 'rugosite du sol RUG**:', nsrf, xmin, xmax
     ENDDO
  ELSE
     PRINT*, 'phyetat0: Le champ <RUG> est present'
     PRINT*, '          J ignore donc les autres RUG**'
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO i = 1, klon
        xmin = MIN(frugs(i, 1), xmin)
        xmax = MAX(frugs(i, 1), xmax)
     ENDDO
     PRINT*, 'rugosite <RUG>', xmin, xmax
     DO nsrf = 2, nbsrf
        DO i = 1, klon
           frugs(i, nsrf) = frugs(i, 1)
        ENDDO
     ENDDO
  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
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon
           xmin = MIN(agesno(i, nsrf), xmin)
           xmax = MAX(agesno(i, nsrf), xmax)
        ENDDO
        PRINT*, 'Age de la neige AGESNO**:', nsrf, xmin, xmax
     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)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zmea(i), xmin)
     xmax = MAX(zmea(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zmea:', xmin, xmax

  CALL get_field("ZSTD", zstd)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zstd(i), xmin)
     xmax = MAX(zstd(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zstd:', xmin, xmax

  CALL get_field("ZSIG", zsig)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zsig(i), xmin)
     xmax = MAX(zsig(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zsig:', xmin, xmax

  CALL get_field("ZGAM", zgam)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zgam(i), xmin)
     xmax = MAX(zgam(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zgam:', xmin, xmax

  CALL get_field("ZTHE", zthe)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zthe(i), xmin)
     xmax = MAX(zthe(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zthe:', xmin, xmax

  CALL get_field("ZPIC", zpic)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zpic(i), xmin)
     xmax = MAX(zpic(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zpic:', xmin, xmax

  CALL get_field("ZVAL", zval)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(zval(i), xmin)
     xmax = MAX(zval(i), xmax)
  ENDDO
  PRINT*, 'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax

  CALL get_field("RUGSREL", rugoro)
  xmin = 1.0E+20
  xmax = -1.0E+20
  DO i = 1, klon
     xmin = MIN(rugoro(i), xmin)
     xmax = MAX(rugoro(i), xmax)
  ENDDO
  PRINT*, 'Rugosite relief (ecart-type) rugsrel:', xmin, xmax

  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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(clwcon)
  xmax = MAXval(clwcon)
  PRINT*, 'Eau liquide convective (ecart-type) clwcon:', xmin, xmax

  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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(rnebcon)
  xmax = MAXval(rnebcon)
  PRINT*, 'Nebulosite convective (ecart-type) rnebcon:', xmin, xmax

  ! 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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(ratqs)
  xmax = MAXval(ratqs)
  PRINT*, '(ecart-type) ratqs:', xmin, xmax

  ! 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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(run_off_lic_0)
  xmax = MAXval(run_off_lic_0)
  PRINT*, '(ecart-type) run_off_lic_0:', xmin, xmax

  ! 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
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO k = 1, klev+1
           DO i = 1, klon
              xmin = MIN(pbl_tke(i, k, nsrf), xmin)
              xmax = MAX(pbl_tke(i, k, nsrf), xmax)
           ENDDO
        ENDDO
        PRINT*, 'Turbulent kinetic energyl TKE**:', nsrf, xmin, xmax
     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
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO k = 1, klev+1
        DO i = 1, klon
          xmin = MIN(wake_delta_pbl_tke(i,k,nsrf),xmin)
          xmax = MAX(wake_delta_pbl_tke(i,k,nsrf),xmax)
        ENDDO
      ENDDO
      PRINT*,'TKE difference (w)-(x) DELTATKE**:', nsrf, xmin, xmax
    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
        xmin = 1.0E+20
        xmax = -1.0E+20
         DO i = 1, klon
            xmin = MIN(delta_tsurf(i, nsrf), xmin)
            xmax = MAX(delta_tsurf(i, nsrf), xmax)
         ENDDO
        PRINT*, 'delta_tsurf:', xmin, xmax
     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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(zmax0)
  xmax = MAXval(zmax0)
  PRINT*, '(ecart-type) zmax0:', xmin, xmax

  !           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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(f0)
  xmax = MAXval(f0)
  PRINT*, '(ecart-type) f0:', xmin, xmax

  ! 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
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO k = 1, klev
        DO i = 1, klon
           xmin = MIN(sig1(i, k), xmin)
           xmax = MAX(sig1(i, k), xmax)
        ENDDO
     ENDDO
     PRINT*, 'sig1:', xmin, xmax
  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
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO k = 1, klev
        DO i = 1, klon
           xmin = MIN(w01(i, k), xmin)
           xmax = MAX(w01(i, k), xmax)
        ENDDO
     ENDDO
     PRINT*, 'w01:', xmin, xmax
  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
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO k = 1, klev
        DO i = 1, klon
           xmin = MIN(wake_deltat(i, k), xmin)
           xmax = MAX(wake_deltat(i, k), xmax)
        ENDDO
     ENDDO
     PRINT*, 'wake_deltat:', xmin, xmax
  ENDIF

  ! wake_deltaq

  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
     xmin = 1.0E+20
     xmax = -1.0E+20
     DO k = 1, klev
        DO i = 1, klon
           xmin = MIN(wake_deltaq(i, k), xmin)
           xmax = MAX(wake_deltaq(i, k), xmax)
        ENDDO
     ENDDO
     PRINT*, 'wake_deltaq:', xmin, xmax
  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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(wake_s)
  xmax = MAXval(wake_s)
  PRINT*, '(ecart-type) wake_s:', xmin, xmax

  ! 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
  xmin = 1.0E+20
  xmax = -1.0E+20
  xmin = MINval(wake_cstar)
  xmax = MAXval(wake_cstar)
  PRINT*, '(ecart-type) wake_cstar:', xmin, xmax

  ! 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
      ALLOCATE(tslab(klon,nslay), stat=ierr)
      IF (ierr /= 0) CALL abort_gcm &
         ('phyetat0', 'pb allocation tslab', 1)
      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)=ftsol(:,is_oce)
          END DO
      END IF 
      print*, "calling slab_init"
      CALL ocean_slab_init(dtime, pctsrf)
  END IF ! Slab        

  ! on ferme le fichier
  CALL close_startphy

  ! Initialize module pbl_surface_mod 

  CALL pbl_surface_init(qsol, fder, snow, qsurf, &
       evap, frugs, agesno, 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