source: trunk/LMDZ.COMMON/libf/dyn3d_common/dynetat0.F @ 1300

!
! $Id $
!
      SUBROUTINE dynetat0(fichnom,vcov,ucov,
     .                    teta,q,masse,ps,phis,time0)

      USE infotrac, only: tname, nqtot
      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror,
     &                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,
     &                  nf90_inquire_dimension,nf90_close

      use control_mod, only : planet_type, timestart

      IMPLICIT NONE

c=======================================================================
c
c   Auteur:  P. Le Van / L.Fairhead
c   -------
c
c   objet:
c   ------
c
c   Lecture de l'etat initial
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

#include "dimensions.h"
#include "paramet.h"
#include "temps.h"
#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"
#include "ener.h"
#include "netcdf.inc"
#include "description.h"
#include "serre.h"
#include "logic.h"
#include "iniprint.h"

c   Arguments:
c   ----------

      CHARACTER(len=*),INTENT(IN) :: fichnom
      REAL,INTENT(OUT) :: vcov(iip1,jjm,llm)
      REAL,INTENT(OUT) :: ucov(iip1,jjp1,llm)
      REAL,INTENT(OUT) :: teta(iip1,jjp1,llm)
      REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot)
      REAL,INTENT(OUT) :: masse(iip1,jjp1,llm)
      REAL,INTENT(OUT) :: ps(iip1,jjp1)
      REAL,INTENT(OUT) :: phis(iip1,jjp1)
      REAL,INTENT(OUT) :: time0

c   Variables 
c
      INTEGER length,iq
      PARAMETER (length = 100)
      REAL tab_cntrl(length) ! tableau des parametres du run
      INTEGER ierr, nid, nvarid

      character(len=12) :: start_file_type="earth" ! default start file type
      INTEGER idecal


      REAL,ALLOCATABLE :: time(:) ! times stored in start
      INTEGER timelen ! number of times stored in the file
      INTEGER indextime ! index of selected time
      !REAL  hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1

      INTEGER edges(4),corner(4)
      integer :: i

c-----------------------------------------------------------------------

c  Ouverture NetCDF du fichier etat initial

      ierr=nf90_open(fichnom,NF90_NOWRITE,nid)
      IF (ierr.NE.nf90_noerr) THEN
        write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc'
        write(lunout,*)trim(nf90_strerror(ierr))
        CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

c
      ierr = nf90_inq_varid (nid, "controle", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <controle> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, tab_cntrl)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echoue pour <controle>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      !!! AS: idecal is a hack to be able to read planeto starts...
      !!!     .... while keeping everything OK for LMDZ EARTH
      if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then
          write(lunout,*)'dynetat0 : Planeto-like start file'
          start_file_type="planeto"
          idecal = 4
          annee_ref  = 2000
      else
          write(lunout,*)'dynetat0 : Earth-like start file'
          idecal = 5
          annee_ref  = tab_cntrl(5)
      endif


      im         = tab_cntrl(1)
      jm         = tab_cntrl(2)
      lllm       = tab_cntrl(3)
      if (start_file_type.eq."earth") then
        day_ref    = tab_cntrl(4)
      else
        day_ini    = tab_cntrl(4)
        day_ref=0
      endif
      rad        = tab_cntrl(idecal+1)
      omeg       = tab_cntrl(idecal+2)
      g          = tab_cntrl(idecal+3)
      cpp        = tab_cntrl(idecal+4)
      kappa      = tab_cntrl(idecal+5)
      daysec     = tab_cntrl(idecal+6)
      dtvr       = tab_cntrl(idecal+7)
      etot0      = tab_cntrl(idecal+8)
      ptot0      = tab_cntrl(idecal+9)
      ztot0      = tab_cntrl(idecal+10)
      stot0      = tab_cntrl(idecal+11)
      ang0       = tab_cntrl(idecal+12)
      pa         = tab_cntrl(idecal+13)
      preff      = tab_cntrl(idecal+14)
c
      clon       = tab_cntrl(idecal+15)
      clat       = tab_cntrl(idecal+16)
      grossismx  = tab_cntrl(idecal+17)
      grossismy  = tab_cntrl(idecal+18)
c
      IF ( tab_cntrl(idecal+19).EQ.1. )  THEN
        fxyhypb  = . TRUE .
c        dzoomx   = tab_cntrl(25)
c        dzoomy   = tab_cntrl(26)
c        taux     = tab_cntrl(28)
c        tauy     = tab_cntrl(29)
      ELSE
        fxyhypb = . FALSE .
        ysinus  = . FALSE .
        IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = . TRUE. 
      ENDIF

      if (planet_type=="mars") then ! so far this is only for Mars
        hour_ini = tab_cntrl(29)
      else
        hour_ini=0
      endif

      if (start_file_type.eq."earth") then
        day_ini = tab_cntrl(30)
        itau_dyn = tab_cntrl(31)
        start_time = tab_cntrl(32)
      else
        day_ini=tab_cntrl(4)
        itau_dyn=0
        start_time=0
      endif
c   .................................................................
c
c
      write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa ',
     &               rad,omeg,g,cpp,kappa

      IF(   im.ne.iim           )  THEN
          PRINT 1,im,iim
          STOP
      ELSE  IF( jm.ne.jjm       )  THEN
          PRINT 2,jm,jjm
          STOP
      ELSE  IF( lllm.ne.llm     )  THEN
          PRINT 3,lllm,llm
          STOP
      ENDIF

      ierr=nf90_inq_varid(nid, "rlonu", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <rlonu> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, rlonu)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "rlatu", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <rlatu> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, rlatu)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "rlonv", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <rlonv> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, rlonv)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "rlatv", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <rlatv> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, rlatv)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour rlatv"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "cu", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <cu> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, cu)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <cu>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "cv", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <cv> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, cv)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <cv>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "aire", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <aire> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, aire)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <aire>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr = nf90_inq_varid (nid, "phisinit", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <phisinit> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr = nf90_get_var(nid, nvarid, phis)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

! read time axis
      ierr = nf90_inq_varid (nid, "temps", nvarid)
      IF (ierr .NE. nf90_noerr) THEN
        write(lunout,*)"dynetat0: Le champ <temps> est absent"
        write(lunout,*)"dynetat0: J essaie <Time>"
        ierr = nf90_inq_varid (nid, "Time", nvarid)
        IF (ierr .NE. nf90_noerr) THEN
           write(lunout,*)"dynetat0: Le champ <Time> est absent"
           write(lunout,*)trim(nf90_strerror(ierr))
           CALL ABORT_gcm("dynetat0", "", 1)
        ENDIF
        ! Get the length of the "Time" dimension
        ierr = nf90_inq_dimid(nid,"Time",nvarid)
        ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
        allocate(time(timelen))
        ! Then look for the "Time" variable
        ierr  =nf90_inq_varid(nid,"Time",nvarid)
        ierr = nf90_get_var(nid, nvarid, time)
        IF (ierr .NE. nf90_noerr) THEN
           write(lunout,*)"dynetat0: Lecture echouee <Time>"
           write(lunout,*)trim(nf90_strerror(ierr))
           CALL ABORT_gcm("dynetat0", "", 1)
        ENDIF
      ELSE   
        ! Get the length of the "temps" dimension
        ierr = nf90_inq_dimid(nid,"temps",nvarid)
        ierr = nf90_inquire_dimension(nid,nvarid,len=timelen)
        allocate(time(timelen))
        ! Then look for the "temps" variable
        ierr = nf90_inq_varid (nid, "temps", nvarid)
        ierr = nf90_get_var(nid, nvarid, time)
        IF (ierr .NE. nf90_noerr) THEN
           write(lunout,*)"dynetat0: Lecture echouee <temps>"
           write(lunout,*)trim(nf90_strerror(ierr))
           CALL ABORT_gcm("dynetat0", "", 1)
        ENDIF
      ENDIF

! select the desired time
      IF (timestart .lt. 0) THEN  ! default: we use the last time value
        indextime = timelen
      ELSE  ! else we look for the desired value in the time axis
       indextime = 0
        DO i=1,timelen
          IF (abs(time(i) - timestart) .lt. 0.01) THEN
             indextime = i
             EXIT
          ENDIF
        ENDDO
        IF (indextime .eq. 0) THEN
          write(lunout,*)"Time", timestart," is not in "
     &                                      //trim(fichnom)//"!!"
          write(lunout,*)"Stored times are:"
          DO i=1,timelen
             PRINT*, time(i)
          ENDDO
          CALL ABORT_gcm("dynetat0", "", 1)
        ENDIF
      ENDIF

      if (planet_type=="mars") then
        ! In start the absolute date is day_ini + hour_ini + time
        ! For now on, in the GCM dynamics, it is day_ini + time0
        time0 = time(indextime) + hour_ini
        day_ini = day_ini + INT(time0)
        time0 = time0 - INT(time0) ! time0 devient le nouveau hour_ini
        hour_ini = time0
      else
        time0 = time(indextime) 
      endif
      
      PRINT*, "dynetat0: Selected time ",time(indextime),
     .        " at index ",indextime
      
      DEALLOCATE(time)

! read vcov
      corner(1)=1
      corner(2)=1
      corner(3)=1
      corner(4)=indextime
      edges(1)=iip1
      edges(2)=jjm
      edges(3)=llm
      edges(4)=1
      ierr=nf90_inq_varid(nid,"vcov",nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <vcov> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr=nf90_get_var(nid,nvarid,vcov,corner,edges)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <vcov>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

! read ucov
      corner(1)=1
      corner(2)=1
      corner(3)=1
      corner(4)=indextime
      edges(1)=iip1
      edges(2)=jjp1
      edges(3)=llm
      edges(4)=1
      ierr=nf90_inq_varid(nid,"ucov",nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <ucov> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr=nf90_get_var(nid,nvarid,ucov,corner,edges)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <ucov>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
 
! read teta (same corner/edges as ucov)
      ierr=nf90_inq_varid(nid,"teta",nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <teta> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr=nf90_get_var(nid,nvarid,teta,corner,edges)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <teta>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

! read tracers (same corner/edges as ucov)
      IF(nqtot.GE.1) THEN
      DO iq=1,nqtot
        ierr= nf90_inq_varid(nid,tname(iq),nvarid)
        IF (ierr .NE. nf90_noerr) THEN
           write(lunout,*)"dynetat0: Le traceur <"//trim(tname(iq))//
     &                    "> est absent"
           write(lunout,*)"          Il est donc initialise a zero"
           q(:,:,:,iq)=0.
        ELSE
           ierr=nf90_get_var(nid,nvarid,q(:,:,:,iq),corner,edges)
          IF (ierr .NE. nf90_noerr) THEN
            write(lunout,*)"dynetat0: Lecture echouee pour "
     &                                //trim(tname(iq))
            write(lunout,*)trim(nf90_strerror(ierr))
            CALL ABORT_gcm("dynetat0", "", 1)
          ENDIF
        ENDIF
      ENDDO
      ENDIF

!read masse (same corner/edges as ucov)
      ierr=nf90_inq_varid(nid,"masse",nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <masse> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr=nf90_get_var(nid,nvarid,masse,corner,edges)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <masse>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

! read ps
      corner(1)=1
      corner(2)=1
      corner(3)=indextime
      edges(1)=iip1
      edges(2)=jjp1
      edges(3)=1
      ierr=nf90_inq_varid(nid,"ps",nvarid)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Le champ <ps> est absent"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF
      ierr=nf90_get_var(nid,nvarid,ps,corner,edges)
      IF (ierr .NE. nf90_noerr) THEN
         write(lunout,*)"dynetat0: Lecture echouee pour <ps>"
         write(lunout,*)trim(nf90_strerror(ierr))
         CALL ABORT_gcm("dynetat0", "", 1)
      ENDIF

      ierr=nf90_close(nid)

      if (planet_type/="mars") then
       day_ini=day_ini+INT(time0) ! obsolete stuff ; 0<time<1 anyways
       time0=time0-INT(time0)
      endif

  1   FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem
     *arrage est differente de la valeur parametree iim =',i4//)
   2  FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem
     *arrage est differente de la valeur parametree jjm =',i4//)
   3  FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema
     *rrage est differente de la valeur parametree llm =',i4//)
   4  FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema
     *rrage est differente de la valeur  dtinteg =',i4//)

      RETURN
      END