source: trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90 @ 1661

!
! $Id $
!
subroutine phyetat0(fichnom)
! Load initial state for the physics
! and do some resulting initializations

      USE dimphy
      USE mod_grid_phy_lmdz
      USE mod_phys_lmdz_para
      USE iophy
      USE phys_state_var_mod
      USE iostart
      use geometry_mod, only: longitude_deg, latitude_deg
      USE time_phylmdz_mod, only: itau_phy, raz_date

implicit none
!======================================================================
! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
! Objet: Lecture de l'etat initial pour la physique
!======================================================================
#include "netcdf.inc"
#include "dimsoil.h"
#include "clesphys.h"
#include "tabcontrol.h"
!======================================================================

character(len=*),intent(in) :: fichnom ! input file name
REAL    :: xmin, xmax
LOGICAL :: found
REAL    :: tab_cntrl(length)
integer :: i,isoil
CHARACTER(len=2) :: str2
REAL :: lon_startphy(klon), lat_startphy(klon)

! les variables globales lues dans le fichier restart

! open physics initial state file:
call open_startphy(fichnom)

!
! Lecture des parametres de controle:
!
      CALL get_var("controle",tab_cntrl,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <controle> est absent'
         CALL abort
      ENDIF
       
      DO i = 1, length
           tabcntr0( i ) = tab_cntrl( i )
      ENDDO


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

      itau_phy = tab_cntrl(15)

! Attention si raz_date est active :
! il faut remettre a zero itau_phy apres phyetat0 !
      IF (raz_date.eq.1) THEN
        itau_phy=0
      ENDIF

! read latitudes and make a sanity check (because already known from dyn)
call get_field("latitude",lat_startphy,found)
IF (.not.found) THEN
  PRINT*, 'phyetat0: Le champ <latitude> est absent'
  CALL abort
ENDIF
DO i=1,klon
  IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.01) THEN
    WRITE(*,*) "phyetat0: Warning! Latitude discrepancy wrt startphy file:",&
               " i=",i," lat_startphy(i)=",lat_startphy(i),&
               " latitude_deg(i)=",latitude_deg(i)
    CALL abort
  ENDIF
ENDDO

! read longitudes and make a sanity check (because already known from dyn)
call get_field("longitude",lon_startphy,found)
IF (.not.found) THEN
  PRINT*, 'phyetat0: Le champ <longitude> est absent'
  CALL abort
ENDIF
DO i=1,klon
  IF (ABS(lon_startphy(i)-longitude_deg(i))>=0.01) THEN
    WRITE(*,*) "phyetat0: Warning! Longitude discrepancy wrt startphy file:",&
               " i=",i," lon_startphy(i)=",lon_startphy(i),&
               " longitude_deg(i)=",longitude_deg(i)
    CALL abort
  ENDIF
ENDDO

! read in other variables here ...

! Lecture des temperatures du sol:

       CALL get_field("TS",ftsol(:),found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <TS> est absent'
         PRINT*, "phyetat0: Lecture echouee pour <TS>"
         CALL abort
      ELSE
         PRINT*, 'phyetat0: Le champ <TS> est present'
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon
            xmin = MIN(ftsol(i),xmin)
            xmax = MAX(ftsol(i),xmax)
         ENDDO
         PRINT*,'Temperature du sol <TS>', xmin, xmax
      ENDIF


! Lecture des temperatures du sol profond:

      DO isoil=1, nsoilmx
      IF (isoil.GT.99) THEN
         PRINT*, "Trop de couches"
         CALL abort
      ENDIF
      WRITE(str2,'(i2.2)') isoil
      CALL get_field('Tsoil'//str2,ftsoil(:,isoil),found)
      IF (.not.found) THEN
         PRINT*, "phyetat0: Le champ <Tsoil"//str2//"> est absent"
         PRINT*, "          Il prend donc la valeur de surface"
         DO i=1, klon
             ftsoil(i,isoil)=ftsol(i)
         ENDDO
      ENDIF
      ENDDO

! Lecture de albedo au sol:

      CALL get_field("ALBE", falbe,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ALBE> est absent'
         PRINT*, "phyetat0: Lecture echouee pour <ALBE>"
         CALL abort
      ELSE
         xmin = 1.0E+20
         xmax = -1.0E+20
         DO i = 1, klon
            xmin = MIN(falbe(i),xmin)
            xmax = MAX(falbe(i),xmax)
         ENDDO
         PRINT*,'Albedo du sol <ALBE>', xmin, xmax
      ENDIF

! Lecture rayonnement solaire au sol:

      CALL get_field("solsw",solsw,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <solsw> est absent'
         PRINT*, 'mis a zero'
         solsw = 0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(solsw(i),xmin)
         xmax = MAX(solsw(i),xmax)
      ENDDO
      PRINT*,'Rayonnement solaire au sol solsw:', xmin, xmax

! Lecture rayonnement IR au sol:

      CALL get_field("sollw",sollw,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <sollw> est absent'
         PRINT*, 'mis a zero'
         sollw = 0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(sollw(i),xmin)
         xmax = MAX(sollw(i),xmax)
      ENDDO
      PRINT*,'Rayonnement IR au sol sollw:', xmin, xmax

! Lecture derive des flux:

      CALL get_field("fder",fder,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <fder> est absent'
         PRINT*, 'mis a zero'
         fder = 0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(fder(i),xmin)
         xmax = MAX(fder(i),xmax)
      ENDDO
      PRINT*,'Derive des flux fder:', xmin, xmax

! Lecture derive flux IR:

      CALL get_field("dlw",dlw,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <dlw> est absent'
         PRINT*, 'mis a zero'
         dlw = 0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(dlw(i),xmin)
         xmax = MAX(dlw(i),xmax)
      ENDDO
      PRINT*,'Derive flux IR dlw:', xmin, xmax

! Lecture rayonnement IR vers le bas au sol:

      CALL get_field("sollwdown",sollwdown,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
         PRINT*, 'mis a zero'
         sollwdown = 0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(sollwdown(i),xmin)
         xmax = MAX(sollwdown(i),xmax)
      ENDDO
      PRINT*,'Flux IR vers le bas au sol sollwdown:', xmin, xmax

! Lecture du rayonnement net au sol:

      CALL get_field("RADS",radsol,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <RADS> est absent'
         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

! Lecture de l'orographie sous-maille:

      CALL get_field("ZMEA",zmea,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
         PRINT*, 'mis a zero'
         zmea=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zmea(i),xmin)
         xmax = MAX(zmea(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zmea:', xmin, xmax

      CALL get_field("ZSTD",zstd,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
         PRINT*, 'mis a zero'
         zstd=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zstd(i),xmin)
         xmax = MAX(zstd(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zstd:', xmin, xmax

      CALL get_field("ZSIG",zsig,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
         PRINT*, 'mis a zero'
         zsig=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zsig(i),xmin)
         xmax = MAX(zsig(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zsig:', xmin, xmax

      CALL get_field("ZGAM",zgam,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
         PRINT*, 'mis a zero'
         zgam=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zgam(i),xmin)
         xmax = MAX(zgam(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zgam:', xmin, xmax

      CALL get_field("ZTHE",zthe,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
         PRINT*, 'mis a zero'
         zthe=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zthe(i),xmin)
         xmax = MAX(zthe(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zthe:', xmin, xmax

      CALL get_field("ZPIC",zpic,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZPIC> est absent'
         PRINT*, 'mis a zero'
         zpic=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zpic(i),xmin)
         xmax = MAX(zpic(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zpic:', xmin, xmax

      CALL get_field("ZVAL",zval,found)
      IF (.not.found) THEN
         PRINT*, 'phyetat0: Le champ <ZVAL> est absent'
         PRINT*, 'mis a zero'
         zval=0.
      ENDIF
      xmin = 1.0E+20
      xmax = -1.0E+20
      DO i = 1, klon
         xmin = MIN(zval(i),xmin)
         xmax = MAX(zval(i),xmax)
      ENDDO
      PRINT*,'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax

! Lecture de TANCIEN:

      ancien_ok = .TRUE.

      CALL get_field("TANCIEN",t_ancien,found)
      IF (.not.found) THEN
         PRINT*, "phyetat0: Le champ <TANCIEN> est absent"
         PRINT*, "Depart legerement fausse. Mais je continue"
         ancien_ok = .FALSE.
      ENDIF

! close file
call close_startphy

! do some more initializations
call init_iophy_new(latitude_deg,longitude_deg)

end subroutine phyetat0