MODULE etat0_dcmip2016_baroclinic_wave_mod
   USE icosa
   USE caldyn_vars_mod
   IMPLICIT NONE
   PRIVATE

   INTEGER, SAVE :: testcase
   !$OMP THREADPRIVATE(testcase)

   INTEGER :: perturbation
   !$OMP THREADPRIVATE(perturbation)

   PUBLIC getin_etat0, compute_etat0

CONTAINS

   SUBROUTINE getin_etat0
      USE mpipara, ONLY: is_mpi_root
      USE tracer_icosa_mod
      IMPLICIT NONE
      LOGICAL :: is_moist
      CHARACTER(LEN=255) :: str_perturbation

      IF (nqtot < 5) THEN
         IF (is_mpi_root) THEN
            PRINT *, "nqtot must be at least 5 for test case dcmip2016_baroclinic_wave"
         END IF
         STOP
      END IF

      str_perturbation = "exponential"
      CALL getin("dcmip2016_baroclinic_wave_perturbation", str_perturbation)
      IF (TRIM(str_perturbation) == "exponential") THEN
         perturbation = 0
      ELSE IF (TRIM(str_perturbation) == "stream") THEN
         perturbation = 1
      END IF

   END SUBROUTINE getin_etat0

   SUBROUTINE compute_etat0(ngrid, lon, lat, phis, ps, vol, ulon, ulat, q)
      USE icosa
      USE disvert_mod
      USE omp_para
      USE dcmip2016_baroclinic_wave_mod, ONLY: baroclinic_wave_test
      USE earth_const
      USE terminator, ONLY: initial_value_Terminator
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ngrid
      REAL(rstd), INTENT(IN) :: lon(ngrid)
      REAL(rstd), INTENT(IN) :: lat(ngrid)
      REAL(rstd), INTENT(OUT) :: phis(ngrid)
      REAL(rstd), INTENT(OUT) :: ps(ngrid)
      REAL(rstd), INTENT(OUT) :: vol(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, PARAMETER :: moist = 1, deep = 0
      INTEGER :: zcoords
      REAL :: p, z
      REAL :: temp, rho, thetav
      INTEGER :: ij, l

      DO ij = 1, ngrid
         z = 0.
         zcoords = 1
     CALL baroclinic_wave_test(deep, moist, perturbation, scale_factor, lon(ij), lat(ij), p, z, zcoords, ulon(ij, 1), ulat(ij, 1), &
                                   temp, thetav, phis(ij), ps(ij), rho, q(ij, 1, 1))

         zcoords = 0
         DO l = ll_begin, ll_end
            p = 0.5*(ap(l) + ap(l + 1) + (bp(l) + bp(l + 1))*ps(ij))
     CALL baroclinic_wave_test(deep, moist, perturbation, scale_factor, lon(ij), lat(ij), p, z, zcoords, ulon(ij, l), ulat(ij, l), &
                                      temp, thetav, phis(ij), ps(ij), rho, q(ij, l, 1))

            vol(ij, l) = 1./rho ! initialize specific volume
            q(ij, l, 2) = 0.
            q(ij, l, 3) = 0.
            CALL initial_value_Terminator(lat(ij), lon(ij), q(ij, l, 4), q(ij, l, 5))
         END DO
      END DO

   END SUBROUTINE compute_etat0

END MODULE etat0_dcmip2016_baroclinic_wave_mod
