source: trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F @ 867

      SUBROUTINE phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq,
     .           day_ini,time,
     .           tsurf,tsoil,emis,q2,qsurf,cloudfrac,totcloudfrac,hice)

      USE surfdat_h
      USE comgeomfi_h
      USE tracer_h

      implicit none

c======================================================================
c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
c  Adaptation à Mars : Yann Wanherdrick 
c Objet: Lecture de l etat initial pour la physique
c======================================================================
#include "netcdf.inc"
#include "dimensions.h"
#include "dimphys.h"
#include "planete.h"
#include "comcstfi.h"

      INTEGER ngrid
c======================================================================
      INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
      PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
!======================================================================
!  Arguments:
!  ---------
!  inputs:
      character*(*) fichnom ! "startfi.nc" file
      integer tab0
      integer Lmodif
      integer nsoil ! # of soil layers
      integer nq
      integer day_ini
      real time

!  outputs:
      real tsurf(ngrid,nbsrf) ! surface temperature
      real tsoil(ngrid,nsoil,nbsrf) ! soil temperature
      real emis(ngrid) ! surface emissivity
      real q2(ngrid, llm+1) ! 
      real qsurf(ngrid,nq) ! tracers on surface
!      real co2ice(ngrid) ! co2 ice cover
      real cloudfrac(ngrid,nlayermx)
      real hice(ngrid), totcloudfrac(ngrid)


!======================================================================
!  Local variables:

!      INTEGER radpas
!      REAL co2_ppm
!      REAL solaire

      real xmin,xmax ! to display min and max of a field
c
      INTEGER ig,iq,lmax
      INTEGER nid, nvarid
      INTEGER ierr, i, nsrf
!      integer isoil 
!      INTEGER length
!      PARAMETER (length=100)
      CHARACTER*7 str7
      CHARACTER*2 str2
      CHARACTER*1 yes
c
      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
      INTEGER nqold

! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
      logical :: oldtracernames=.false.
      integer :: count
      character(len=30) :: txt ! to store some text

!! added variables by AS to allow to read only slices of startfi
      real :: toto(ngrid)
      integer :: sta   !! subscript in starti where we start to read
      integer, dimension(2) :: sta2d
      integer, dimension(2) :: siz2d

c AS: get the cursor that is stored in dimphys.h
c ... this allows to read only a part of startfi horiz grid
      sta = cursor

c
c ALLOCATE ARRAYS IN surfdat_h
c
      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))

c
c Ouvrir le fichier contenant l etat initial:
c

      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
      IF (ierr.NE.NF_NOERR) THEN
        write(6,*)' Pb d''ouverture du fichier '//fichnom
        CALL ABORT
      ENDIF

! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
!                    qsurf02, ...)
      count=0
      do iq=1,nq
        txt= " "
        write(txt,'(a5,i2.2)')'qsurf',iq
        ierr=NF_INQ_VARID(nid,txt,nvarid)
        if (ierr.ne.NF_NOERR) then
          ! did not find old tracer name
          exit ! might as well stop here
        else
          ! found old tracer name
          count=count+1
        endif
      enddo
      if (count.eq.nq) then
        write(*,*) "phyetat0:tracers seem to follow old naming ",
     &             "convention (qsurf01,qsurf02,...)"
        write(*,*) "   => will work for now ... "
        write(*,*) "      but you should run newstart to rename them"
        oldtracernames=.true.
      endif

c modifications possibles des variables de tab_cntrl
      write(*,*) 'TABFI de phyeta0',Lmodif,tab0
      call tabfi (ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
     .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)

c
c Lecture des latitudes (coordonnees):
c
      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_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
#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_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
         CALL abort
      ENDIF
c
c Lecture des aires des mailles:
c
      ierr = NF_INQ_VARID (nid, "area", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <area> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <area>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(area)
      xmax = MAXVAL(area)
      PRINT*,'Aires des mailles <area>:', xmin, xmax
c
c Lecture du geopotentiel au sol:
c
      ierr = NF_INQ_VARID (nid, "phisfi", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <phisfi> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,phisfi)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,phisfi)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <phisfi>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(phisfi)
      xmax = MAXVAL(phisfi)
      PRINT*,'Geopotentiel au sol <phisfi>:', xmin, xmax
c
c Lecture de l''albedo du sol nu:
c
      ierr = NF_INQ_VARID (nid, "albedodat", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <albedodat> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,albedodat)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,albedodat)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <albedodat>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(albedodat)
      xmax = MAXVAL(albedodat)
      PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
c
c ZMEA
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_VARA_DOUBLE(nid,nvarid,sta,ngrid,zmea)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,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, ngrid
         xmin = MIN(zmea(i),xmin)
         xmax = MAX(zmea(i),xmax)
      ENDDO
      PRINT*,'<zmea>:', xmin, xmax
c
c ZSTD
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_VARA_DOUBLE(nid,nvarid,sta,ngrid,zstd)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,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, ngrid
         xmin = MIN(zstd(i),xmin)
         xmax = MAX(zstd(i),xmax)
      ENDDO
      PRINT*,'<zstd>:', xmin, xmax
c
c ZSIG
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_VARA_DOUBLE(nid,nvarid,sta,ngrid,zsig)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,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, ngrid
         xmin = MIN(zsig(i),xmin)
         xmax = MAX(zsig(i),xmax)
      ENDDO
      PRINT*,'<zsig>:', xmin, xmax
c
c ZGAM
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_VARA_DOUBLE(nid,nvarid,sta,ngrid,zgam)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,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, ngrid
         xmin = MIN(zgam(i),xmin)
         xmax = MAX(zgam(i),xmax)
      ENDDO
      PRINT*,'<zgam>:', xmin, xmax
c
c ZTHE
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_VARA_DOUBLE(nid,nvarid,sta,ngrid,zthe)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,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, ngrid
         xmin = MIN(zthe(i),xmin)
         xmax = MAX(zthe(i),xmax)
      ENDDO
      PRINT*,'<zthe>:', xmin, xmax
c
c CO2 ice cover
c
! Ehouarn: from now on, there is no "co2ice" standalone field; it is supposed
! to be with stored in qsurf(i_co2_ice)
!      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
!      IF (ierr.NE.NF_NOERR) THEN
!         PRINT*, 'phyetat0: Le champ <co2ice> est absent'
!         CALL abort
!      ENDIF
!#ifdef NC_DOUBLE
!      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2ice)
!#else
!      ierr = NF_GET_VAR_REAL(nid, nvarid, co2ice)
!#endif
!      IF (ierr.NE.NF_NOERR) THEN
!         PRINT*, 'phyetat0: Lecture echouee pour <co2ice>'
!         CALL abort
!      ENDIF
!      xmin = 1.0E+20
!      xmax = -1.0E+20
!      xmin = MINVAL(co2ice)
!      xmax = MAXVAL(co2ice)
!      PRINT*,'CO2 ice cover <co2ice>:', xmin, xmax
c
c Lecture des temperatures du sol:
c
      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <tsurf> est absent'
         PRINT*, '          Mais je vais essayer de lire TS**'
         IF (nbsrf.GT.99) THEN
            PRINT*, "Trop de sous-mailles"
            CALL abort
         ENDIF
         DO nsrf = 1, nbsrf
           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, tsurf(1,nsrf))
#else
           ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(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
           xmin = MINVAL(tsurf)
           xmax = MAXVAL(tsurf)
           PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
         ENDDO
      ELSE
         PRINT*, 'phyetat0: Le champ <tsurf> est present'
         PRINT*, '          J ignore donc les autres temperatures TS**'
#ifdef NC_DOUBLE
         ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, tsurf)
#else
         ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, tsurf)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
            CALL abort
         ENDIF
         xmin = 1.0E+20
         xmax = -1.0E+20
         xmin = MINVAL(tsurf)
         xmax = MAXVAL(tsurf)
         PRINT*,'Temperature du sol <tsurf>', xmin, xmax
         IF (nbsrf >= 2) THEN
            DO nsrf = 2, nbsrf
               DO i = 1, ngrid
                  tsurf(i,nsrf) = tsurf(i,1)
               ENDDO
            ENDDO
         ENDIF
      ENDIF
c
c Lecture des temperatures du sol profond:
c
!      IF (nsoil.GT.99 .OR. nbsrf.GT.99) THEN
!         PRINT*, "Trop de couches ou sous-mailles"
!         CALL abort
!      ENDIF
!      DO nsrf = 1, nbsrf
!         DO isoil=1, nsoil
!            WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
!            ierr = NF_INQ_VARID (nid, 'tsoil', nvarid)
!            IF (ierr.NE.NF_NOERR) THEN
!               PRINT*, "phyetat0: Le champ <tsoil> est absent"
!               PRINT*, "          Il prend donc la valeur de surface"
!               DO i=1, ngrid
!                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
!               ENDDO
!            ELSE
!#ifdef NC_DOUBLE
!              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,1,nsrf))
!#else
!              ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,1,nsrf))
!#endif
!               IF (ierr.NE.NF_NOERR) THEN
!                  PRINT*, "Lecture echouee pour <tsoil>"
!                  CALL abort
!               ENDIF
!            ENDIF
!         ENDDO
!      ENDDO
!      xmin = 1.0E+20
!      xmax = -1.0E+20
!      xmin = MINVAL(tsoil)
!      xmax = MAXVAL(tsoil)
!      PRINT*,'Temperatures du sol profond <tsoil>', xmin, xmax
c
c Surface emissivity
c
      ierr = NF_INQ_VARID (nid, "emis", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <emis> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta,ngrid, emis)
#else
      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, emis)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <emis>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(emis)
      xmax = MAXVAL(emis)
      PRINT*,'Surface emissivity <emis>:', xmin, xmax

c
c Cloud fraction (added by BC 2010)
c
      ierr = NF_INQ_VARID (nid, "cloudfrac", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <cloudfrac> est absent'
      cloudfrac(:,:)=0.5
!         CALL abort
      ENDIF
      sta2d = (/sta,1/)
      siz2d = (/ngrid, nlayermx/)
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,cloudfrac)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,cloudfrac)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <cloudfrac>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(cloudfrac)
      xmax = MAXVAL(cloudfrac)
      PRINT*,'Cloud fraction <cloudfrac>:', xmin, xmax


c
c Total cloud fraction (added by BC 2010)
c
      ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid)
!      ierr = NF_INQ_VARID (nid, "totalfrac", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <totcloudfrac> est absent'
      totcloudfrac(:)=0.5
!         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, totcloudfrac)
#else
      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, totcloudfrac)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>'
         !CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(totcloudfrac)
      xmax = MAXVAL(totcloudfrac)
      PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax



c
c Height of oceanic ice (added by BC 2010)
c
      ierr = NF_INQ_VARID (nid, "hice", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <hice> est absent'
         CALL abort
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, hice)
#else
      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, hice)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <hice>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(hice)
      xmax = MAXVAL(hice)
      PRINT*,'Height of oceanic ice <hice>:', xmin, xmax


c
c pbl wind variance
c
      ierr = NF_INQ_VARID (nid, "q2", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <q2> est absent'
         CALL abort
      ENDIF
      sta2d = (/sta,1/)
      siz2d = (/ngrid, llm+1/)
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,q2)
#else
      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,q2)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Lecture echouee pour <q2>'
         CALL abort
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      xmin = MINVAL(q2)
      xmax = MAXVAL(q2)
      PRINT*,'pbl wind variance <q2>:', xmin, xmax
c
c tracer on surface
c

      IF(nq.GE.1) THEN
         nqold=nq
         DO iq=1,nq
!            str7(1:5)='qsurf'
!            WRITE(str7(6:7),'(i2.2)') iq
!            ierr = NF_INQ_VARID (nid,str7,nvarid)
           IF (oldtracernames) THEN
             txt=" "
             write(txt,'(a5,i2.2)')'qsurf',iq
           ELSE
             txt=noms(iq)
!             if (txt.eq."h2o_vap") then
               ! There is no surface tracer for h2o_vap;
               ! "h2o_ice" should be loaded instead
!               txt="h2o_ice"
!               write(*,*) 'phyetat0: loading surface tracer',
!     &                     ' h2o_ice instead of h2o_vap'
!             endif
           ENDIF ! of IF (oldtracernames) THEN
           ierr=NF_INQ_VARID(nid,txt,nvarid)
           IF (ierr.NE.NF_NOERR) THEN
             write(*,*) 'PHYETAT0: WARNING : surface tracer',trim(txt),
     &                  ' not found in file'
             write(*,*) trim(txt), ' set to 0'
             do ig=1,ngrid
               qsurf(ig,iq)=0.
             end do
             nqold=min(iq-1,nqold)
           ELSE
           toto(:) = 0.
#ifdef NC_DOUBLE
           ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,toto)
#else
           ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,toto)
#endif
             IF (ierr.NE.NF_NOERR) THEN
               PRINT*, 'phyetat0: Lecture echouee pour <',trim(txt),'>'
               CALL abort
             ENDIF
             qsurf(1:ngrid,iq) = toto(1:ngrid)
           ENDIF
           xmin = 1.0E+20
           xmax = -1.0E+20
           xmin = MINVAL(qsurf(1:ngrid,iq))
           xmax = MAXVAL(qsurf(1:ngrid,iq))
           PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
         ENDDO
         if ((nqold.lt.nq).and.(nqold.ge.1)) then
c        case when new tracer are added in addition to old ones
             write(*,*)'qsurf 1 to ', nqold,'were already present'
             write(*,*)'qsurf ', nqold+1,' to ', nq,'are new'
             write(*,*)' and initialized to zero'
             qsurf(:,nqold+1:nq)=0.0
!            yes=' '
!            do while ((yes.ne.'y').and.(yes.ne.'n'))
!             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
!             write(*,*) 'to #',nq-nqold+1,'->', nq,'   (y or n) ?'
!             read(*,fmt='(a)') yes
!            end do
!            if (yes.eq.'y') then
!              write(*,*) 'OK, let s reindex qsurf'
!                 do ig=1,ngrid
!                    do iq=nq,nq-nqold+1,-1
!                       qsurf(ig,iq)=qsurf(ig,iq-nq+nqold)
!                    end do
!                    do iq=nq-nqold,1,-1
!                       qsurf(ig,iq)= 0.
!                    end do
!                 end do
!            end if
         end if ! of if ((nqold.lt.nq).and.(nqold.ge.1))
      ENDIF ! of IF(nq.GE.1)

! Call to soil_settings, in order to read soil temperatures,
! as well as thermal inertia and volumetric heat capacity

      PRINT*,'SOIL INIT'
      call soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
c
c Fermer le fichier:
c
      ierr = NF_CLOSE(nid)
c

      RETURN
      END