!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phyetat0.F,v 1.2 2004/06/22 11:45:33 lmdzadmin Exp $
!
c
c
      SUBROUTINE phyetat0 (fichnom,dtime,
     .            rlat,rlon, tsol,tsoil,
     .           albe, solsw, sollw,
     .           fder,radsol,resch4,
     .           tabcntr0,
     .           t_ancien,ancien_ok)
c======================================================================
c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
c Objet: Lecture de l'etat initial pour la physique
c======================================================================
      use dimphy
      IMPLICIT none
#include "dimensions.h"
#include "netcdf.inc"
#include "dimsoil.h"
#include "clesphys.h"
#include "temps.h"
c======================================================================
      CHARACTER*(*) fichnom
      REAL dtime
      INTEGER radpas,chimpas
      REAL rlat(klon), rlon(klon)   ! in degrees
      REAL tsol(klon)
      REAL tsoil(klon,nsoilmx)
      REAL albe(klon)
cIM BEG alblw
      REAL alblw(klon)
cIM END alblw
      REAL radsol(klon)
      REAL sollw(klon)
      real solsw(klon)
      real fder(klon)

      REAL t_ancien(klon,klev)
      LOGICAL ancien_ok

      REAL resch4(klon)
      INTEGER ig0

      REAL xmin, xmax
c
      INTEGER nid, nvarid
      INTEGER ierr, i, nsrf, isoil 
      INTEGER length
      PARAMETER (length=100)
      REAL tab_cntrl(length), tabcntr0(length)
      CHARACTER*2 str2
c
c Ouvrir le fichier contenant l'etat initial:
c
      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
c
c Lecture des parametres de controle:
c
      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
      ELSE
c
         DO i = 1, length
           tabcntr0( i ) = tab_cntrl( i )
         ENDDO
c

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

      ENDIF

      itau_phy = tab_cntrl(15)

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_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 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*, "phyetat0: Lecture echouee pour <TS>"
         CALL abort
      ELSE
         PRINT*, 'phyetat0: Le champ <TS> est present'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsol(1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, tsol(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
            xmin = MIN(tsol(i),xmin)
            xmax = MAX(tsol(i),xmax)
         ENDDO
         PRINT*,'Temperature du sol <TS>', xmin, xmax
      ENDIF
c
c Lecture des temperatures du sol profond:
c
      DO isoil=1, nsoilmx
      IF (isoil.GT.99) THEN
         PRINT*, "Trop de couches"
         CALL abort
      ENDIF
      WRITE(str2,'(i2.2)') isoil
      ierr = NF_INQ_VARID (nid, 'Tsoil'//str2, nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <Tsoil"//str2//"> est absent"
         PRINT*, "          Il prend donc la valeur de surface"
         DO i=1, klon
             tsoil(i,isoil)=tsol(i)
         ENDDO
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,isoil))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,isoil))
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "Lecture echouee pour <Tsoil"//str2//">"
            CALL abort
         ENDIF
      ENDIF
      ENDDO
c
c Lecture de albedo au sol:
c
      ierr = NF_INQ_VARID (nid, "ALBE", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'phyetat0: Le champ <ALBE> est absent'
         PRINT*, "phyetat0: Lecture echouee pour <ALBE>"
         CALL abort
      ELSE
         PRINT*, 'phyetat0: Le champ <ALBE> est present'
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albe(1))
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid, albe(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
            xmin = MIN(albe(i),xmin)
            xmax = MAX(albe(i),xmax)
         ENDDO
         PRINT*,'Albedo du sol <ALBE>', xmin, xmax
      ENDIF

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
         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
         xmin = MIN(sollw(i),xmin)
         xmax = MAX(sollw(i),xmax)
      ENDDO
      PRINT*,'Rayonnement IF au sol sollw:', xmin, xmax

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
         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
         xmin = MIN(radsol(i),xmin)
         xmax = MAX(radsol(i),xmax)
      ENDDO
      PRINT*,'Rayonnement net au sol radsol:', xmin, xmax
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     Par defaut on cree 2 bandes de methane au pole Nord et au pole Sud
c     (entre 75 et 85 degres de latitude) de 2 metres.
c     Les poles sont sec !
      resch4(1) = 0.    ! pole nord = 1 point
      DO ig0=2,klon
          if ((rlat(ig0).ge.75..and.rlat(ig0).le.85.).or. 
     &        (rlat(ig0).ge.-85.and.rlat(ig0).le.-75.)) then
            resch4(ig0) = 2.
          else
            resch4(ig0) = 0.
          endif
      ENDDO
      resch4(klon) = 0.   ! pole sud = 1 point

      ierr = NF_INQ_VARID (nid, "RESCH4", nvarid)
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, "phyetat0: Le champ <RESCH4> est absent"
         PRINT*, "Pas de reservoir de methane mais je continue..."
         PRINT*, "Pour info, je met 2 metres de methane sur 2 bandes"
         PRINT*, "comprises entre 75 et 85 degres de latitude dans  "
         PRINT*, "chaque hemisphere."
      ELSE
#ifdef NC_DOUBLE
         ierr = NF_GET_VAR_DOUBLE(nid, nvarid, resch4)
#else
         ierr = NF_GET_VAR_REAL(nid, nvarid,resch4)
#endif
         IF (ierr.NE.NF_NOERR) THEN
            PRINT*, "phyetat0: Lecture echouee pour <reservoir>"
            CALL abort
         ENDIF
      ENDIF
c
c Fermer le fichier:
c
      ierr = NF_CLOSE(nid)
c
      RETURN
      END