MODULE etat0_temperature_mod
  USE prec
  USE icosa, ONLY: llm,nqtot
  IMPLICIT NONE
  PRIVATE

  REAL(rstd), SAVE, ALLOCATABLE :: t_profile(:)
!$OMP THREADPRIVATE(t_profile)

  PUBLIC :: getin_etat0, compute_etat0

CONTAINS
  
  SUBROUTINE getin_etat0
    USE getin_mod, ONLY: getin
    USE mpipara, ONLY: is_mpi_root
    USE omp_para, ONLY: omp_in_parallel
    USE transfert_omp_mod, ONLY: bcast_omp
    USE free_unit_mod, ONLY : free_unit
    USE abort_mod
    INTEGER :: unit,ok
    INTEGER :: l
    CHARACTER(len=255) :: temperature_file
  
    temperature_file="profile.in" ! default file name
    ! but users may want to use some other file name
    CALL getin("temperature_profile_file",temperature_file)
    
    ALLOCATE(t_profile(llm))
 
 
 !$OMP MASTER
    unit=free_unit()
    OPEN(unit,file=temperature_file,status="old",action="read",iostat=ok)
    IF (ok/=0) THEN
      WRITE(*,*) "getin_etat0 error: input file ",trim(temperature_file)," not found!"
      CALL dynamico_abort( "Could not open temperature file." )
    ENDIF
    ! read in t_profile() line by line, starting from first atmospheric
    ! layer, up to model top
    DO l=1,llm
      READ(unit,fmt=*,iostat=ok) t_profile(l)
      IF (ok/=0) THEN
        WRITE(*,*) "getin_etat0 error: failed reading t_profile(l) for l=",l
        CALL dynamico_abort( "Could not read temperature file. " )
      ENDIF
    ENDDO
    
    CLOSE(unit)
    IF (is_mpi_root) THEN
      WRITE(*,*) "Using input temperature profile from file ",trim(temperature_file),":"
      DO l=1,llm
         WRITE(*,*) "  TEMP(l=",l,")=",t_profile(l)
      ENDDO
    ENDIF
!$OMP END MASTER

    IF (omp_in_parallel()) THEN
      CALL bcast_omp(t_profile)
    ENDIF
    
  END SUBROUTINE getin_etat0

  SUBROUTINE compute_etat0(ngrid, phis, ps, temp, ulon, ulat, q)
    USE earth_const, ONLY: preff
    INTEGER, INTENT(IN)    :: ngrid
    REAL(rstd),INTENT(OUT) :: phis(ngrid)
    REAL(rstd),INTENT(OUT) :: ps(ngrid)
    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm)
    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
    INTEGER :: l
    phis(:)=0
    ps(:)=preff
    DO l=1,llm
      temp(:,l)=t_profile(l)
    ENDDO
    ulon(:,:)=0
    ulat(:,:)=0
    q(:,:,:)=0

  END SUBROUTINE compute_etat0

END MODULE etat0_temperature_mod