module phyetat0_mod

implicit none

contains

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

  use dimphy, only: klon
  use iophy, only: init_iophy_new
  use phys_state_var_mod, only: ftsol, ftsoil, age, dlw, falbe, fder, q2, &
                                radsol, sollw, sollwdown, solsw, &
                                ancien_ok, t_ancien, &
                                zmea, zgam, zpic, zsig, zstd, zthe, zval
  use iostart, only: get_var, get_field, open_startphy, close_startphy
  use geometry_mod, only: longitude_deg, latitude_deg
  use time_phylmdz_mod, only: itau_phy, raz_date, pdtphys
  use ioipsl_getin_p_mod, only: getin_p
  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
                                east_gwstress, west_gwstress

implicit none
!======================================================================
include "dimsoil.h"
include "clesphys.h"
include "tabcontrol.h"
!======================================================================

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

! global variables read in startphy.nc file

! open physics initial state file:
if (startphy_file) then
  call open_startphy(fichnom)
endif

!
! Load control parameters:
!
IF (startphy_file) THEN
  CALL get_var("controle",tab_cntrl,found)
  IF (.not.found) THEN
    write(*,*) modname//': Field <controle> is missing'
    call abort_physic(modname,"Missing <controle>",1)
  ENDIF
       
  DO i = 1, length
    tabcntr0( i ) = tab_cntrl( i )
  ENDDO

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

  itau_phy = tab_cntrl(15)

! Warning, if raz_date is active :
! itau_phy must be re-set to zero after phyetat0 !
  IF (raz_date.eq.1) THEN
    itau_phy=0
  ENDIF

ELSE
  tabcntr0(:)=1 ! dummy initialization
  ! Initialize parameter or get values from def files
  dtime=pdtphys
  radpas=1
  itau_phy=0
ENDIF ! of IF (startphy_file)

IF (startphy_file) THEN
  ! read latitudes and make a sanity check (because already known from dyn)
  call get_field("latitude",lat_startphy,found)
  IF (.not.found) THEN
    write(*,*) modname//':: Field <latitude> is missing'
    call abort_physic(modname,"Missing <latitude>",1)
  ENDIF
  DO i=1,klon
    IF (ABS(lat_startphy(i)-latitude_deg(i))>=0.01) THEN
      WRITE(*,*) modname//": Warning! Latitude discrepancy wrt startphy file:",&
                 " i=",i," lat_startphy(i)=",lat_startphy(i),&
                 " latitude_deg(i)=",latitude_deg(i)
      call abort_physic(modname,"<latitude> values discrepency",1)
    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
    write(*,*)'phyetat0: Field <longitude> is missing'
    call abort_physic(modname,"Missing <longitude>",1)
  ENDIF
  DO i=1,klon
    IF (ABS(lon_startphy(i)-longitude_deg(i))>=0.01) THEN
      WRITE(*,*) modname//": Warning! Longitude discrepancy wrt startphy file:",&
                 " i=",i," lon_startphy(i)=",lon_startphy(i),&
                 " longitude_deg(i)=",longitude_deg(i)
      call abort_physic(modname,"<latitude> values discrepency",1)
    ENDIF
  ENDDO
ENDIF ! of IF (startphy_file)

! read in other variables here ...

IF (startphy_file) THEN
  ! Load surface temperature:
  CALL get_field("TS",ftsol(:),found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <TS> is missing"
    call abort_physic(modname,"Missing <TS>",1)
  ELSE
    WRITE(*,*) modname//": Field <TS> is present"
    WRITE(*,*) 'Surface temperature <TS>', minval(ftsol), maxval(ftsol)
  ENDIF
ELSE
  ! Dummy initialization, but in fact this is later handled in physiq
  ftsol(:)=0
ENDIF ! of IF (startphy_file)

IF (startphy_file) THEN
  ! Load sub-surface temperatures:
  DO isoil=1, nsoilmx
    IF (isoil.GT.99) THEN
       WRITE(*,*) "Too many layers!"
       call abort_physic(modname,"Too many subsurface layers",1)
    ENDIF
    WRITE(str2,'(i2.2)') isoil
    CALL get_field('Tsoil'//str2,ftsoil(:,isoil),found)
    IF (.not.found) THEN
      WRITE(*,*) modname//": Field <Tsoil"//str2//"> is missing"
      WRITE(*,*) "           it will be nitialized with surface value"
      DO i=1, klon
             ftsoil(i,isoil)=ftsol(i)
      ENDDO
    ENDIF
  ENDDO
ELSE
  ! Dummy initialization, but in fact this is later handled in physiq
  ftsoil(:,:)=0
ENDIF ! of IF (startphy_file)

IF (startphy_file) THEN
  ! Load surface albedo:
  CALL get_field("ALBE", falbe,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ALBE> is missing"
    call abort_physic(modname,"Missing <ALBE>",1)
  ENDIF
ELSE
  ! Dummy initialization: read value from def file 
  surface_albedo=0.5 ! default
  CALL getin_p("surface_albedo",surface_albedo)
  falbe(:)=surface_albedo
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Surface albedo <ALBE>', minval(falbe), maxval(falbe)

IF (startphy_file) THEN
  ! Solar flux to the surface:
  CALL get_field("solsw",solsw,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <solsw> is missing"
    WRITE(*,*) "set to zero"
    solsw = 0.
  ENDIF
ELSE
  ! Dummy initialization
  solsw(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Solar flux to the surface solsw:', minval(solsw), maxval(solsw)

IF (startphy_file) THEN
  ! IR flux to the surface:
  CALL get_field("sollw",sollw,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <sollw> is missing"
    WRITE(*,*) "set to zero"
    sollw = 0.
  ENDIF
ELSE
  ! Dummy initialization
  sollw(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'IR flux to the surface sollw:', minval(sollw), maxval(solsw)

IF (startphy_file) THEN
  ! Derivative of the sensible and latent fluxes:
  CALL get_field("fder",fder,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <fder> is missing"
    WRITE(*,*) "set to zero"
    fder = 0.
  ENDIF
ELSE
  ! Dummy initialization
  fder(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Derivative of the flux fder:', minval(fder), maxval(fder)

IF (startphy_file) THEN
  ! Derivative of the IR flux:
  CALL get_field("dlw",dlw,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <dlw> is missing"
    WRITE(*,*) "set to zero"
    dlw = 0.
  ENDIF
ELSE
  ! Dummy initialization
  dlw(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Derivative of the IR flux dlw:', minval(dlw), maxval(dlw)

IF (startphy_file) THEN
  ! Downward IR flux at the surface:
  CALL get_field("sollwdown",sollwdown,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <sollwdown> is missing"
    WRITE(*,*) "set to zero"
    sollwdown = 0.
  ENDIF
ELSE
  ! Dummy initialization
  sollwdown(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Downward IR flux at the surface sollwdown:', minval(sollwdown), maxval(sollwdown)

IF (startphy_file) THEN
  ! Net flux at the surface:
  CALL get_field("RADS",radsol,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <RADS> is missing"
    call abort_physic(modname,"Missing <RADS>",1)
  ENDIF
ELSE
  ! Dummy initialization
  radsol(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Net flux at the surface radsol:', minval(radsol), maxval(radsol)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZMEA",zmea,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZMEA> is missing"
    WRITE(*,*) "set to zero"
    zmea=0.
  ENDIF
ELSE
  zmea(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zmea:', minval(zmea), maxval(zmea)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZSTD",zstd,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZSTD> is missing"
    WRITE(*,*) "set to zero"
    zstd=0.
  ENDIF
ELSE
  zstd(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zstd:', minval(zstd), maxval(zstd)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZSIG",zsig,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZSIG> is missing"
    WRITE(*,*) "set to zero"
    zsig=0.
  ENDIF
ELSE
  zsig(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zsig:', minval(zsig), maxval(zsig)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZGAM",zgam,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZGAM> is missing"
    WRITE(*,*) "set to zero"
    zgam=0.
  ENDIF
ELSE
  zgam(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zgam:', minval(zgam), maxval(zgam)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZTHE",zthe,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZTHE> is missing"
    WRITE(*,*) "set to zero"
    zthe=0.
  ENDIF
ELSE
  zthe(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zthe:', minval(zthe), maxval(zthe)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZPIC",zpic,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZPIC> is missing"
    WRITE(*,*) "set to zero"
    zpic=0.
  ENDIF
ELSE
  zpic(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zpic:', minval(zpic), maxval(zpic)

IF (startphy_file) THEN
  ! Load sub-grid scale orography parameters:
  CALL get_field("ZVAL",zval,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <ZVAL> is missing"
    WRITE(*,*) "set to zero"
    zval=0.
  ENDIF
ELSE
  zval(:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Subgrid scale orography zval:', minval(zval), maxval(zval)

IF (startphy_file) THEN
  ! Lecture de TANCIEN:
  ancien_ok = .TRUE.

  CALL get_field("TANCIEN",t_ancien,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <TANCIEN> is missing"
    WRITE(*,*) "Slightly biaised start. But let's keep going."
    ancien_ok = .FALSE.
  ENDIF
ELSE
  ancien_ok=.false.
ENDIF

IF (startphy_file) THEN
  CALL get_field("age",age,found)
  IF (.not.found) THEN
    PRINT*, "phyetat0: Age of air is missing"
    PRINT*, "Reinitialising age of air to 0"
    age(:,:) = 0.
  ENDIF
ELSE
  age(:,:) = 0.
ENDIF

IF (startphy_file) THEN
  ! Load Q2 the TKE at interlayer:
  CALL get_field("Q2",q2,found)
  IF (.not.found) THEN
    WRITE(*,*) modname//": Field <Q2> is missing"
    WRITE(*,*) "set to zero"
    q2(:,:)=0.
  ENDIF
ELSE
  ! Dummy initialization
  q2(:,:)=0
ENDIF ! of IF (startphy_file)
WRITE(*,*) 'Turbulent Kinetic Energy', minval(q2), maxval(q2)

! Non-orographic gravity waves 
if (startphy_file) then
   call get_field("du_nonoro_gwd",du_nonoro_gwd,found)
   if (.not.found) then
      write(*,*) "phyetat0: <du_nonoro_gwd> not in file"
      du_nonoro_gwd(:,:)=0.
   endif
endif ! of if (startphy_file)
if (startphy_file) then
   call get_field("dv_nonoro_gwd",dv_nonoro_gwd,found)
   if (.not.found) then
      write(*,*) "phyetat0: <dv_nonoro_gwd> not in file"
      dv_nonoro_gwd(:,:)=0.
   endif
endif ! of if (startphy_file)
if (startphy_file) then
   call get_field("east_gwstress",east_gwstress,found)
   if (.not.found) then
      write(*,*) "phyetat0: <east_gwstress> not in file"
      east_gwstress(:,:)=0.
   endif
endif ! of if (startphy_file)
if (startphy_file) then
   call get_field("west_gwstress",west_gwstress,found)
   if (.not.found) then
      write(*,*) "phyetat0: <west_gwstress> not in file"
      west_gwstress(:,:)=0.
   endif
endif ! of if (startphy_file)

! close file
IF (startphy_file) call close_startphy

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

end subroutine phyetat0

end module phyetat0_mod