!
! $Header$
!
c
c
      SUBROUTINE phyetat0 (fichnom,dtime,co2_ppm_etat0,solaire_etat0,
     .           rlat_p, rlon_p, pctsrf_p, tsol_p,
     .           ocean_in, ok_veget_in,
     .           alb1_p, alb2_p,
     .           rain_fall_p, snow_fall_p,solsw_p, sollw_p,
     .           radsol_p,clesphy0,
     .           zmea_p,zstd_p,zsig_p,zgam_p,zthe_p,zpic_p,zval_p,
     .           rugsrel_p,tabcntr0,
     .           t_ancien_p,q_ancien_p,ancien_ok_p, rnebcon_p, ratqs_p,
     .           clwcon_p)

      USE dimphy
      USE mod_grid_phy_lmdz
      USE mod_phys_lmdz_para
      USE iophy
      USE ocean_slab_mod,   ONLY : ocean_slab_init
      USE ocean_cpl_mod,    ONLY : ocean_cpl_init
      USE ocean_forced_mod, ONLY : ocean_forced_init
      USE fonte_neige_mod,  ONLY : fonte_neige_init
      USE pbl_surface_mod,  ONLY : pbl_surface_init
      USE surface_data,     ONLY : ocean, ok_veget

      IMPLICIT none
c======================================================================
c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
c Objet: Lecture de l'etat initial pour la physique
c======================================================================
#include "dimensions.h"
#include "netcdf.inc"
#include "indicesol.h"
#include "dimsoil.h"
#include "clesphys.h"
#include "temps.h"
#include "thermcell.h"
#include "compbl.h"
c======================================================================
      CHARACTER*(*) fichnom
      REAL dtime
      INTEGER radpas
      REAL rlat_p(klon), rlon_p(klon)
      REAL co2_ppm_etat0
      REAL solaire_etat0
      REAL tsol_p(klon,nbsrf)
      REAL tsoil_p(klon,nsoilmx,nbsrf)
      REAL tslab_p(klon), seaice_p(klon)
      REAL qsurf_p(klon,nbsrf)
      REAL qsol_p(klon)
      REAL snow_p(klon,nbsrf)
      REAL alb1_p(klon,nbsrf)   ! albedo in visible SW interval
      REAL alb2_p(klon,nbsrf)   ! albedo in near IR interval
      REAL evap_p(klon,nbsrf)
      REAL radsol_p(klon)
      REAL rain_fall_p(klon)
      REAL snow_fall_p(klon)
      REAL sollw_p(klon)
      real solsw_p(klon)
      real fder_p(klon)
      REAL frugs_p(klon,nbsrf)
      REAL agesno_p(klon,nbsrf)
      REAL zmea_p(klon)
      REAL zstd_p(klon)
      REAL zsig_p(klon)
      REAL zgam_p(klon)
      REAL zthe_p(klon)
      REAL zpic_p(klon)
      REAL zval_p(klon)
      REAL rugsrel_p(klon)
      REAL pctsrf_p(klon, nbsrf)
      REAL run_off_lic_0_p(klon)
      
      REAL t_ancien_p(klon,klev), q_ancien_p(klon,klev)
      real rnebcon_p(klon,klev),clwcon_p(klon,klev)
      real ratqs_p(klon,klev)
      LOGICAL,SAVE ::  ancien_ok
      LOGICAL ::  ancien_ok_p
      
      REAL zmasq_glo(klon_glo)
      REAL rlat(klon_glo), rlon(klon_glo)
      REAL tsol(klon_glo,nbsrf)
      REAL tsoil(klon_glo,nsoilmx,nbsrf)
cIM "slab" ocean
      REAL tslab(klon_glo), seaice(klon_glo)
      REAL qsurf(klon_glo,nbsrf)
      REAL qsol(klon_glo)
      REAL snow(klon_glo,nbsrf)
      REAL alb1(klon_glo,nbsrf)
      REAL alb2(klon_glo,nbsrf)
      REAL evap(klon_glo,nbsrf)
      REAL radsol(klon_glo)
      REAL rain_fall(klon_glo)
      REAL snow_fall(klon_glo)
      REAL sollw(klon_glo)
      real solsw(klon_glo)
      real fder(klon_glo)
      REAL frugs(klon_glo,nbsrf)
      REAL agesno(klon_glo,nbsrf)
      REAL zmea(klon_glo)
      REAL zstd(klon_glo)
      REAL zsig(klon_glo)
      REAL zgam(klon_glo)
      REAL zthe(klon_glo)
      REAL zpic(klon_glo)
      REAL zval(klon_glo)
      REAL rugsrel(klon_glo)
      REAL pctsrf(klon_glo, nbsrf)
      REAL fractint(klon_glo)
      REAL run_off_lic_0(klon_glo)
      REAL t_ancien(klon_glo,klev)
      REAL q_ancien(klon_glo,klev)
      real rnebcon(klon_glo,klev)
      real clwcon(klon_glo,klev)
      real ratqs(klon_glo,klev)

      CHARACTER*6 ocean_in
      LOGICAL ok_veget_in

      INTEGER        longcles
      PARAMETER    ( longcles = 20 )
      REAL clesphy0( longcles )
c
      REAL xmin, xmax
c
      INTEGER nid, nvarid
      INTEGER ierr, i, nsrf, isoil ,k
      INTEGER length
      PARAMETER (length=100)
      REAL tab_cntrl(length), tabcntr0(length)
      REAL,SAVE :: tab_cntrl_omp(length)
      CHARACTER*7 str7
      CHARACTER*2 str2

c FH1D
c     real iolat(jjm+1)
      real iolat(jjm+1-1/iim)
c
c Ouvrir le fichier contenant l'etat initial:
c

c$OMP MASTER
      print *,'MASTER -x , omp_rank=',omp_rank
c$OMP END MASTER

c$OMP MASTER
      IF (is_mpi_root) THEN
        print*,'fichnom ',fichnom
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
        IF (ierr.NE.NF_NOERR) THEN
          write(6,*)' Pb d''ouverture du fichier '//fichnom
          write(6,*)' ierr = ', ierr
          CALL ABORT
        ENDIF
      ENDIF
c$OMP END MASTER
c
c Lecture des parametres de controle:
c
c$OMP MASTER
      IF (is_mpi_root) THEN
      
      ierr = NF_INQ_VARID (nid, "controle", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <controle> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <controle>'
         CALL abort
      ENDIF
      ENDIF

c$OMP END MASTER
       
       CALL bcast(tab_cntrl)
       
c
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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
c
         tab_cntrl(1)=dtime
         tab_cntrl(2)=radpas
         co2_ppm_etat0      = tab_cntrl(3)
         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)

       

         IF( clesphy0(1).NE.tab_cntrl( 5 ) )  THEN
             clesphy0(1)=tab_cntrl( 5 )
         ENDIF

         IF( clesphy0(2).NE.tab_cntrl( 6 ) )  THEN
             clesphy0(2)=tab_cntrl( 6 )
         ENDIF

         IF( clesphy0(3).NE.tab_cntrl( 7 ) )  THEN
             clesphy0(3)=tab_cntrl( 7 )
         ENDIF

         IF( clesphy0(4).NE.tab_cntrl( 8 ) )  THEN
             clesphy0(4)=tab_cntrl( 8 )
         ENDIF

         IF( clesphy0(5).NE.tab_cntrl( 9 ) )  THEN
             clesphy0(5)=tab_cntrl( 9 )
         ENDIF

         IF( clesphy0(6).NE.tab_cntrl( 10 ) )  THEN
             clesphy0(6)=tab_cntrl( 10 )
         ENDIF

         IF( clesphy0(7).NE.tab_cntrl( 11 ) )  THEN
             clesphy0(7)=tab_cntrl( 11 )
         ENDIF

         IF( clesphy0(8).NE.tab_cntrl( 12 ) )  THEN
             clesphy0(8)=tab_cntrl( 12 )
         ENDIF


c
c Lecture des latitudes (coordonnees):
c
c$OMP MASTER
      IF (is_mpi_root) THEN
      
      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <latitude> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlat)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, rlat)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
         CALL abort
      ENDIF

c
c Lecture des longitudes (coordonnees):
c
      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <longitude> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rlon)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, rlon)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
         CALL abort
      ENDIF
C
C
C Lecture du masque terre mer
C

      ierr = NF_INQ_VARID (nid, "masque", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmasq_glo)
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, zmasq_glo)
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <masque>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <masque> est absent'
          PRINT*, 'fichier startphy non compatible avec phyetat0'
C      CALL abort
      ENDIF

       
C Lecture des fractions pour chaque sous-surface
C
C initialisation des sous-surfaces
C
      pctsrf = 0.
C
C fraction de terre
C

      ierr = NF_INQ_VARID (nid, "FTER", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid,        
     .                             pctsrf(1 : klon_glo,is_ter))
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, 
     .                           pctsrf(1 : klon_glo,is_ter))
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <FTER>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <FTER> est absent'
c@$$         CALL abort
      ENDIF

C
C fraction de glace de terre
C
      ierr = NF_INQ_VARID (nid, "FLIC", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, 
     .                             pctsrf(1 : klon_glo,is_lic))
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, 
     .                           pctsrf(1 : klon_glo,is_lic))
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <FLIC>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <FLIC> est absent'
c@$$         CALL abort
      ENDIF
C
C fraction d'ocean
C
      ierr = NF_INQ_VARID (nid, "FOCE", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, 
     .                             pctsrf(1 : klon_glo,is_oce))
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, 
     .                           pctsrf(1 : klon_glo,is_oce))
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <FOCE>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <FOCE> est absent'
c@$$         CALL abort
      ENDIF

C
C fraction glace de mer
C
      ierr = NF_INQ_VARID (nid, "FSIC", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, 
     .                             pctsrf(1 : klon_glo,is_sic))
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, 
     .                           pctsrf(1 : klon_glo, is_sic))
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <FSIC>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <FSIC> est absent'
c@$$         CALL abort
      ENDIF

C
C  Verification de l'adequation entre le masque et les sous-surfaces
C
      fractint( 1 : klon_glo) = pctsrf(1 : klon_glo, is_ter) 
     $    + pctsrf(1 : klon_glo, is_lic)
      DO i = 1 , klon_glo
        IF ( abs(fractint(i) - zmasq_glo(i) ) .GT. EPSFRA ) THEN
            WRITE(*,*) 'phyetat0: attention fraction terre pas ', 
     $          'coherente ', i, zmasq_glo(i), pctsrf(i, is_ter)
     $          ,pctsrf(i, is_lic)
        ENDIF 
      END DO 
      fractint (1 : klon_glo) =  pctsrf(1 : klon_glo, is_oce) 
     $    + pctsrf(1 : klon_glo, is_sic)
      DO i = 1 , klon_glo
        IF ( abs( fractint(i) - (1. - zmasq_glo(i))) .GT. EPSFRA ) THEN
            WRITE(*,*) 'phyetat0 attention fraction ocean pas ', 
     $          'coherente ', i, zmasq_glo(i) , pctsrf(i, is_oce)
     $          ,pctsrf(i, is_sic)
        ENDIF 
      END DO 

C
c Lecture des temperatures du sol:
c

      ierr = NF_INQ_VARID (nid, "TS", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "TS"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <TS"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsol(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, tsol(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <TS"//str2//">"
             CALL abort
           ENDIF

           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              xmin = MIN(tsol(i,nsrf),xmin)
              xmax = MAX(tsol(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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsol(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, tsol(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <TS>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            xmin = MIN(tsol(i,1),xmin)
            xmax = MAX(tsol(i,1),xmax)
         ENDDO
         PRINT*,'Temperature du sol <TS>', xmin, xmax
         DO nsrf = 2, nbsrf
         DO i = 1, klon_glo
            tsol(i,nsrf) = tsol(i,1)
         ENDDO
         ENDDO
      ENDIF

c
c Lecture des temperatures du sol profond:
c
      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
      ENDIF
      WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
      ierr = NF_INQ_VARID (nid, 'Tsoil'//str7, nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <Tsoil"//str7//"> est absent"
         PRINT*, "          Il prend donc la valeur de surface"
         DO i=1, klon_glo
             tsoil(i,isoil,nsrf)=tsol(i,nsrf)
         ENDDO
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,isoil,nsrf))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,isoil,nsrf))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "Lecture echouee pour <Tsoil"//str7//">"
            CALL abort
         ENDIF
      ENDIF
      ENDDO
      ENDDO
c
cIM "slab" ocean 
c
c Lecture de tslab (pour slab ocean seulement):      
c
      IF (ocean_in .eq. 'slab  ') then
        ierr = NF_INQ_VARID (nid, "TSLAB", nvarid)
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, "phyetat0: Le champ <TSLAB> est absent"
          CALL abort
        ENDIF
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tslab)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, tslab)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, "phyetat0: Lecture echouee pour <TSLAB>"
          CALL abort
        ENDIF
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon_glo
          xmin = MIN(tslab(i),xmin)
          xmax = MAX(tslab(i),xmax)
        ENDDO
        PRINT*,'Min, Max tslab (utilise si OCEAN=slab )', xmin, xmax
c
c Lecture de seaice (pour slab ocean seulement):
c
        ierr = NF_INQ_VARID (nid, "SEAICE", nvarid)
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, "phyetat0: Le champ <SEAICE> est absent"
          CALL abort
        ENDIF
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, seaice)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, seaice)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, "phyetat0: Lecture echouee pour <SEAICE>"
          CALL abort
        ENDIF
        xmin = 1.0E+20
        xmax = -1.0E+20
        DO i = 1, klon_glo
          xmin = MIN(seaice(i),xmin)
          xmax = MAX(seaice(i),xmax)
        ENDDO
        PRINT*,'Masse de la glace de mer (utilise si OCEAN=slab)', 
     $  xmin, xmax
      ELSE
        tslab = 0.
        seaice = 0.
      ENDIF
c
c Lecture de l'humidite de l'air juste au dessus du sol:
c
      ierr = NF_INQ_VARID (nid, "QS", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "QS"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <QS"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, qsurf(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, qsurf(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <QS"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, qsurf(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, qsurf(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <QS>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            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_glo
            qsurf(i,nsrf) = qsurf(i,1)
         ENDDO
         ENDDO
      ENDIF
C
C Eau dans le sol (pour le modele de sol "bucket")
C
      ierr = NF_INQ_VARID (nid, "QSOL", nvarid)
      IF (ierr .EQ.  NF_NOERR) THEN
#ifdef NC_DOUBLE
          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, qsol)
#else
          ierr = NF_GET_VAR_REAL(nid, nvarid, qsol)
#endif
          IF (ierr.NE.NF_NOERR) THEN
              PRINT*, 'phyetat0: Lecture echouee pour <QSOL>'
              CALL abort
          ENDIF
      else
          PRINT*, 'phyetat0: Le champ <QSOL> est absent'
          PRINT*, '          Valeur par defaut nulle'
          qsol(:)=0.
c@$$         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
        xmin = MIN(qsol(i),xmin)
        xmax = MAX(qsol(i),xmax)
      ENDDO
      PRINT*,'Eau dans le sol (mm) <QSOL>', xmin, xmax
c
c Lecture de neige au sol:
c
      ierr = NF_INQ_VARID (nid, "SNOW", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "SNOW"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <SNOW"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, snow(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, snow(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <SNOW"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, snow(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, snow(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <SNOW>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            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_glo
            snow(i,nsrf) = snow(i,1)
         ENDDO
         ENDDO
      ENDIF
c
c Lecture de albedo de l'interval visible au sol:
c
      ierr = NF_INQ_VARID (nid, "ALBE", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "ALBE"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <ALBE"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb1(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, alb1(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <ALBE"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              xmin = MIN(alb1(i,nsrf),xmin)
              xmax = MAX(alb1(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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb1(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, alb1(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <ALBE>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            xmin = MIN(alb1(i,1),xmin)
            xmax = MAX(alb1(i,1),xmax)
         ENDDO
         PRINT*,'Neige du sol <ALBE>', xmin, xmax
         DO nsrf = 2, nbsrf
         DO i = 1, klon_glo
            alb1(i,nsrf) = alb1(i,1)
         ENDDO
         ENDDO
      ENDIF

c
c Lecture de albedo au sol dans l'interval proche infra-rouge:
c
      ierr = NF_INQ_VARID (nid, "ALBLW", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ALBLW> est absent'
c        PRINT*, '          Mais je vais essayer de lire ALBLW**'
         PRINT*, '          Mais je vais prendre ALBE**'
         DO nsrf = 1, nbsrf
           DO i = 1, klon_glo
             alb2(i,nsrf) = alb1(i,nsrf)
           ENDDO
         ENDDO
      ELSE
         PRINT*, 'phyetat0: Le champ <ALBLW> est present'
         PRINT*, '          J ignore donc les autres ALBLW**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, alb2(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, alb2(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <ALBLW>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            xmin = MIN(alb2(i,1),xmin)
            xmax = MAX(alb2(i,1),xmax)
         ENDDO
         PRINT*,'Neige du sol <ALBLW>', xmin, xmax
         DO nsrf = 2, nbsrf
         DO i = 1, klon_glo
            alb2(i,nsrf) = alb2(i,1)
         ENDDO
         ENDDO
      ENDIF
c
c Lecture de evaporation:  
c
      ierr = NF_INQ_VARID (nid, "EVAP", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "EVAP"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <EVAP"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <EVAP"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, evap(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, evap(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <EVAP>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            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_glo
            evap(i,nsrf) = evap(i,1)
         ENDDO
         ENDDO
      ENDIF
c
c Lecture precipitation liquide:
c
      ierr = NF_INQ_VARID (nid, "rain_f", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <rain_f> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rain_fall)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, rain_fall)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <rain_f>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(rain_fall(i),xmin)
         xmax = MAX(rain_fall(i),xmax)
      ENDDO
      PRINT*,'Precipitation liquide rain_f:', xmin, xmax
c
c Lecture precipitation solide:
c
      ierr = NF_INQ_VARID (nid, "snow_f", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <snow_f> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, snow_fall)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, snow_fall)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <snow_f>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(snow_fall(i),xmin)
         xmax = MAX(snow_fall(i),xmax)
      ENDDO
      PRINT*,'Precipitation solide snow_f:', xmin, xmax
c
c Lecture rayonnement solaire au sol:
c
      ierr = NF_INQ_VARID (nid, "solsw", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <solsw> est absent'
         PRINT*, 'mis a zero'
         solsw = 0.
      ELSE
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, solsw)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, solsw)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <solsw>'
          CALL abort
        ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(solsw(i),xmin)
         xmax = MAX(solsw(i),xmax)
      ENDDO
      PRINT*,'Rayonnement solaire au sol solsw:', xmin, xmax
c
c Lecture rayonnement IF au sol:
c
      ierr = NF_INQ_VARID (nid, "sollw", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <sollw> est absent'
         PRINT*, 'mis a zero'
         sollw = 0.
      ELSE
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, sollw)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, sollw)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <sollw>'
          CALL abort
        ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(sollw(i),xmin)
         xmax = MAX(sollw(i),xmax)
      ENDDO
      PRINT*,'Rayonnement IF au sol sollw:', xmin, xmax
      
      ENDIF  ! is_mpi_root
c$OMP END MASTER


c$OMP MASTER
      IF (is_mpi_root) THEN
c
c Lecture derive des flux:
c
      ierr = NF_INQ_VARID (nid, "fder", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <fder> est absent'
         PRINT*, 'mis a zero'
         fder = 0.
      ELSE
#ifdef NC_DOUBLE
        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, fder)
#else
        ierr = NF_GET_VAR_REAL(nid, nvarid, fder)
#endif
        IF (ierr.NE.NF_NOERR) THEN
          PRINT*, 'phyetat0: Lecture echouee pour <fder>'
          CALL abort
        ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(fder(i),xmin)
         xmax = MAX(fder(i),xmax)
      ENDDO
      PRINT*,'Derive des flux fder:', xmin, xmax

c
c Lecture du rayonnement net au sol:
c
      ierr = NF_INQ_VARID (nid, "RADS", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <RADS> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, radsol)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, radsol)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <RADS>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(radsol(i),xmin)
         xmax = MAX(radsol(i),xmax)
      ENDDO
      PRINT*,'Rayonnement net au sol radsol:', xmin, xmax
c
c Lecture de la longueur de rugosite 
c
c
      ierr = NF_INQ_VARID (nid, "RUG", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "RUG"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <RUG"//str2//"> est absent"
              CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, frugs(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, frugs(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <RUG"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, frugs(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, frugs(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <RUG>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            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_glo
            frugs(i,nsrf) = frugs(i,1)
         ENDDO
         ENDDO
      ENDIF

c
c Lecture de l'age de la neige:
c
      ierr = NF_INQ_VARID (nid, "AGESNO", nvarid)
      IF (ierr.NE.NF_NOERR) 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
           ENDIF
           WRITE(str2,'(i2.2)') nsrf
           ierr = NF_INQ_VARID (nid, "AGESNO"//str2, nvarid)
           IF (ierr.NE.NF_NOERR) THEN
              PRINT*, "phyetat0: Le champ <AGESNO"//str2//"> est absent"
              agesno = 50.0
           ENDIF
#ifdef NC_DOUBLE
           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, agesno(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, agesno(1,nsrf))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, "phyetat0: Lecture echouee pour <AGESNO"//str2//">"
             CALL abort
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           DO i = 1, klon_glo
              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**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, agesno(1,1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, agesno(1,1))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <AGESNO>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon_glo
            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_glo
            agesno(i,nsrf) = agesno(i,1)
         ENDDO
         ENDDO
      ENDIF

c
      ierr = NF_INQ_VARID (nid, "ZMEA", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmea)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zmea)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZMEA>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zmea(i),xmin)
         xmax = MAX(zmea(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zmea:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "ZSTD", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zstd)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zstd)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZSTD>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zstd(i),xmin)
         xmax = MAX(zstd(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zstd:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "ZSIG", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zsig)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zsig)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZSIG>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zsig(i),xmin)
         xmax = MAX(zsig(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zsig:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "ZGAM", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zgam)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zgam)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZGAM>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zgam(i),xmin)
         xmax = MAX(zgam(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zgam:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "ZTHE", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zthe)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zthe)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZTHE>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zthe(i),xmin)
         xmax = MAX(zthe(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zthe:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "ZPIC", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZPIC> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zpic)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zpic)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZPIC>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zpic(i),xmin)
         xmax = MAX(zpic(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zpic:', xmin, xmax
c
      ierr = NF_INQ_VARID (nid, "ZVAL", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ZVAL> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zval)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, zval)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <ZVAL>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(zval(i),xmin)
         xmax = MAX(zval(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax
c
c
      ierr = NF_INQ_VARID (nid, "RUGSREL", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <RUGSREL> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rugsrel)
#else
      ierr = NF_GET_VAR_REAL(nid, nvarid, rugsrel)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <RUGSREL>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon_glo
         xmin = MIN(rugsrel(i),xmin)
         xmax = MAX(rugsrel(i),xmax)
      ENDDO
      PRINT*,'Rugosite relief (ecart-type) rugsrel:', xmin, xmax
c
c
      ancien_ok = .TRUE.
c
      ierr = NF_INQ_VARID (nid, "TANCIEN", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <TANCIEN> est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         ancien_ok = .FALSE.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, t_ancien)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, t_ancien)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <TANCIEN>"
            CALL abort
         ENDIF
      ENDIF
c
      ierr = NF_INQ_VARID (nid, "QANCIEN", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <QANCIEN> est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         ancien_ok = .FALSE.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q_ancien)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, q_ancien)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <QANCIEN>"
            CALL abort
         ENDIF
      ENDIF
c
      clwcon=0.
      ierr = NF_INQ_VARID (nid, "CLWCON", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ CLWCON est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         clwcon = 0.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, clwcon)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, clwcon)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <CLWCON>"
            CALL abort
         ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINval(clwcon)
      xmax = MAXval(clwcon)
      PRINT*,'Eau liquide convective (ecart-type) clwcon:', xmin, xmax
c
      rnebcon=0.
      ierr = NF_INQ_VARID (nid, "RNEBCON", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ RNEBCON est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         rnebcon = 0.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, rnebcon)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, rnebcon)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <RNEBCON>"
            CALL abort
         ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINval(rnebcon)
      xmax = MAXval(rnebcon)
      PRINT*,'Nebulosite convective (ecart-type) rnebcon:', xmin, xmax

c
c Lecture ratqs
c
       ratqs=0.
      ierr = NF_INQ_VARID (nid, "RATQS", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <RATQS> est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         ratqs = 0.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, ratqs)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, ratqs)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <RATQS>"
            CALL abort
         ENDIF
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINval(ratqs)
      xmax = MAXval(ratqs)
      PRINT*,'(ecart-type) ratqs:', xmin, xmax
c
c Lecture run_off_lic_0
c
      ierr = NF_INQ_VARID (nid, "RUNOFFLIC0", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <RUNOFFLIC0> est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         run_off_lic_0 = 0.
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, run_off_lic_0)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, run_off_lic_0)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <RUNOFFLIC0>"
            CALL abort
         ENDIF
      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


c Fermer le fichier:
c
      ierr = NF_CLOSE(nid)
      ENDIF ! is_mpi_root
c
c$OMP END MASTER

c$OMP MASTER
cym  en attendant mieux
        iolat(1)=rlat(1)

!FH1D   
!iolat(jjm+1)=rlat(klon_glo)
	iolat(jjm+1-1/iim)=rlat(klon_glo)
        if (iim.gt.1) then
	do i=2,jjm
	  iolat(i)=rlat(2+(i-2)*iim)
	enddo
        endif

	CALL bcast_mpi(iolat)
	CALL bcast_mpi(rlon)

!FH1D
!       call init_iophy(iolat,rlon(2:iim+1))
        call init_iophy(iolat,rlon(2-1/iim:iim+1-1/iim))
       
c$OMP END MASTER
      	
      call Scatter( rlat,rlat_p)
      call Scatter( rlon,rlon_p)
      call Scatter( tsol,tsol_p)
      call Scatter( tsoil,tsoil_p)
      call Scatter( tslab,tslab_p)
      call Scatter( seaice,seaice_p)
      call Scatter( qsurf,qsurf_p)
      call Scatter( qsol,qsol_p)
      call Scatter( snow,snow_p)
      call Scatter( alb1,alb1_p)
      call Scatter( alb2,alb2_p)
      call Scatter( evap,evap_p)
      call Scatter( radsol,radsol_p)
      call Scatter( rain_fall,rain_fall_p)
      call Scatter( snow_fall,snow_fall_p)
      call Scatter( sollw,sollw_p)
      call Scatter( solsw,solsw_p)
      call Scatter( fder,fder_p)
      call Scatter( frugs,frugs_p)
      call Scatter( agesno,agesno_p)
      call Scatter( zmea,zmea_p)
      call Scatter( zstd,zstd_p)
      call Scatter( zsig,zsig_p)
      call Scatter( zgam,zgam_p)
      call Scatter( zthe,zthe_p)
      call Scatter( zpic,zpic_p)
      call Scatter( zval,zval_p)
      call Scatter( rugsrel,rugsrel_p)
      call Scatter( pctsrf,pctsrf_p)
      call Scatter( run_off_lic_0,run_off_lic_0_p)
      call Scatter( t_ancien,t_ancien_p)
      call Scatter( q_ancien,q_ancien_p)
      call Scatter( rnebcon,rnebcon_p)
      call Scatter( clwcon,clwcon_p)
      call Scatter( ratqs,ratqs_p)
      call Scatter( zmasq_glo,zmasq)

c
c Initilalize variables in module surface_data
c
      ok_veget = ok_veget_in
      ocean    = ocean_in
c
c Initialize module pbl_surface_mod 
c
      CALL pbl_surface_init(qsol_p, fder_p, snow_p, qsurf_p,
     $     evap_p, frugs_p, agesno_p, tsoil_p)

c Initialize ocean module according to ocean type
      IF ( ocean == 'slab' ) THEN
c        initilalize module ocean_slab_init
         CALL ocean_slab_init(dtime, tslab_p, seaice_p, pctsrf_p)
      ELSEIF ( ocean == 'couple' ) THEN
c        initilalize module ocean_cpl_init
         CALL ocean_cpl_init(dtime, rlon_p, rlat_p)
      ELSE
c        initilalize module ocean_forced_init
         CALL ocean_forced_init
      ENDIF
c
c Initilialize module fonte_neige_mod      
c
      CALL fonte_neige_init(run_off_lic_0_p)


      RETURN
      END