MODULE etat0_database_mod
   IMPLICIT NONE

CONTAINS

   SUBROUTINE init_etat0
      USE xios_mod
      USE omp_para
      USE getin_mod, ONLY: getin
      IMPLICIT NONE
      CHARACTER(LEN=256) :: etat0_database_type

      ! @getin_doc etat0_database_type realm=initial type=String default=`legacy`
      !  When `etat0=="database"`, selects type of external file(s) from which initial state is read. Valid values :
      !  * `legacy` : reads from XIOS group `read_files_legacy`, see `xml/context_dynamico.xml`
      !  * `ERA5_forcing` : reads from XIOS group `read_files_ERA5_forcing`, see `xml/context_dynamico.xml`

      etat0_database_type = "legacy"
      CALL getin("etat0_database_type", etat0_database_type)

      IF (is_omp_master) THEN
         CALL xios_set_fieldgroup_attr("read_fields", enabled=.TRUE.)
         CALL xios_set_filegroup_attr("read_files", enabled=.TRUE.)
      END IF
   END SUBROUTINE init_etat0

   SUBROUTINE etat0(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
      USE icosa
      USE restart_mod
      USE wind_from_lonlat_mod
      USE write_field_mod
      USE time_mod
      USE transfert_mod
      USE xios_mod
      USE write_field_mod
      USE vertical_remap_mod
      USE theta2theta_rhodz_mod
      USE qsat_mod
      USE pression_mod, ONLY: pression
      USE omp_para
      TYPE(t_field), POINTER :: f_ps(:)
      TYPE(t_field), POINTER :: f_phis(:)
      TYPE(t_field), POINTER :: f_theta_rhodz(:)
      TYPE(t_field), POINTER :: f_u(:)
      TYPE(t_field), POINTER :: f_q(:)

      TYPE(t_field), POINTER, SAVE :: f_ulon_reg(:)
      TYPE(t_field), POINTER, SAVE :: f_ulat_reg(:)
      TYPE(t_field), POINTER, SAVE :: f_temp_reg(:)
      TYPE(t_field), POINTER, SAVE :: f_q_reg(:)

      TYPE(t_field), POINTER, SAVE :: f_ts(:)
      TYPE(t_field), POINTER, SAVE :: f_z(:)
      TYPE(t_field), POINTER, SAVE :: f_ulon(:)
      TYPE(t_field), POINTER, SAVE :: f_ulat(:)
      TYPE(t_field), POINTER, SAVE :: f_temp(:)
      TYPE(t_field), POINTER, SAVE :: f_q1(:)
      TYPE(t_field), POINTER, SAVE :: f_qsat(:)
      TYPE(t_field), POINTER, SAVE :: f_p(:)
      INTEGER :: nb_level
      REAL, ALLOCATABLE:: levels(:)
      INTEGER :: ind

      CALL xios_read_field("relief_db", f_phis)

      CALL writeField("relief_out", f_phis, once=.TRUE.)

      IF (is_omp_level_master) THEN
         DO ind = 1, ndomain
            IF (.NOT. assigned_domain(ind)) CYCLE
            f_phis(ind)%rval2d(:) = f_phis(ind)%rval2d(:)*g
         END DO
      END IF
!$OMP BARRIER

      IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", n_glo=nb_level)
      CALL bcast_omp(nb_level)
      ALLOCATE (levels(nb_level))

      IF (is_omp_master) CALL xios_get_axis_attr("lev_ecdyn", value=levels)
      CALL bcast_omp(levels)

      levels = levels*100  ! hectoPascal -> Pascal

      CALL allocate_field(f_ts, field_t, type_real, name="ts")
      CALL allocate_field(f_z, field_t, type_real, name="z")
      CALL allocate_field(f_ulon_reg, field_t, type_real, nb_level)
      CALL allocate_field(f_ulat_reg, field_t, type_real, nb_level)
      CALL allocate_field(f_temp_reg, field_t, type_real, nb_level)
      CALL allocate_field(f_q_reg, field_t, type_real, nb_level)

      CALL allocate_field(f_q1, field_t, type_real, llm)
      CALL allocate_field(f_qsat, field_t, type_real, llm)
      CALL allocate_field(f_p, field_t, type_real, llm + 1)
      CALL allocate_field(f_temp, field_t, type_real, llm)
      CALL allocate_field(f_ulon, field_t, type_real, llm)
      CALL allocate_field(f_ulat, field_t, type_real, llm)

      CALL xios_read_field("z_db", f_z)
      CALL xios_read_field("ps_db", f_ps)
      CALL xios_read_field("ts_db", f_ts)
      CALL writeField("ps_out", f_ps)

!$OMP BARRIER

!    CALL writeField("phis_out",f_phis,once=.TRUE.)
!    CALL writeField("ts_out",f_ts,once=.TRUE.)

! make correction to ps due to relief at higher resolution
! difference with LMDZ : tsol is taken from ECDYN.NC and not from ECPHY
      IF (is_omp_level_master) THEN
         DO ind = 1, ndomain
            IF (.NOT. assigned_domain(ind)) CYCLE
            f_ps(ind)%rval2d(:) = f_ps(ind)%rval2d(:)*(1.+(f_z(ind)%rval2d(:) - f_phis(ind)%rval2d(:))/287.0/f_ts(ind)%rval2d(:))
         END DO
      END IF
!$OMP BARRIER
      CALL transfert_request(f_ps, req_i0)
      CALL writeField("ps_out", f_ps)

      CALL xios_read_field("temp_db", f_temp_reg)
      CALL vertical_remap(levels, f_temp_reg, f_ps, f_temp)
      CALL transfert_request(f_temp, req_i0)
      CALL temperature2theta_rhodz(f_ps, f_temp, f_theta_rhodz)

      CALL xios_read_field("u_db", f_ulon_reg)
      CALL vertical_remap(levels, f_ulon_reg, f_ps, f_ulon)
      CALL xios_read_field("v_db", f_ulat_reg)
      CALL vertical_remap(levels, f_ulat_reg, f_ps, f_ulat)
      CALL transfert_request(f_ulat, req_i0)
      CALL transfert_request(f_ulon, req_i0)
      CALL ulonlat2un(f_ulon, f_ulat, f_u)

      CALL xios_read_field("q_db", f_q_reg)
      CALL vertical_remap(levels, f_q_reg, f_ps, f_q1)

      CALL pression(f_ps, f_p)
! difference with LMDZ : for qsat, pressure at mid layer is computed as a mean value pmid=(p(l)+p(l+1))/2
      CALL qsat(f_temp, f_p, f_qsat)
      CALL transfert_request(f_qsat, req_i0)

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         f_q(ind)%rval4d(:, :, :) = 1e-6
         f_q(ind)%rval4d(:, :, 1) = f_q1(ind)%rval3d(:, :)*f_qsat(ind)%rval3d(:, :)*0.01
         WHERE (f_q(ind)%rval4d(:, :, 1) < 0) f_q(ind)%rval4d(:, :, 1) = 0
      END DO

      CALL writeField("tempdb_out", f_temp_reg)
      CALL writeField("temp_out", f_temp)

      CALL deallocate_field(f_ts)
      CALL deallocate_field(f_z)
      CALL deallocate_field(f_ulon_reg)
      CALL deallocate_field(f_ulat_reg)
      CALL deallocate_field(f_temp_reg)
      CALL deallocate_field(f_q_reg)

      CALL deallocate_field(f_q1)
      CALL deallocate_field(f_qsat)
      CALL deallocate_field(f_p)
      CALL deallocate_field(f_temp)
      CALL deallocate_field(f_ulon)
      CALL deallocate_field(f_ulat)

   END SUBROUTINE etat0

END MODULE etat0_database_mod
