MODULE etat0_dcmip2_mod
! test cases DCMIP 2012, category 2 : Orographic gravity waves
   USE icosa
   IMPLICIT NONE
   PRIVATE

   INTEGER, SAVE :: testcase
   INTEGER, PARAMETER :: mountain = 0, schaer_noshear = 1, schaer_shear = 2
   REAL(rstd):: Teq_test2 = 300.d0, h0_test2 = 250.d0

   PUBLIC getin_etat0, compute_etat0, Teq_test2, h0_test2

CONTAINS

   SUBROUTINE getin_etat0
      CHARACTER(len=255) :: etat0_type
      etat0_type = 'jablonowsky06'
      CALL getin("etat0", etat0_type)
      CALL getin("dcmip2_Teq", Teq_test2)
      CALL getin("dcmip2_h0", h0_test2)
      SELECT CASE (TRIM(etat0_type))
      CASE ('dcmip2_mountain')
         testcase = mountain
      CASE ('dcmip2_schaer_noshear')
         testcase = schaer_noshear
      CASE ('dcmip2_schaer_shear')
         testcase = schaer_shear
      CASE DEFAULT
         PRINT *, 'This should not happen : etat0_type =', TRIM(etat0_type), ' in etat0_dcmip2.f90/etat0'
         STOP
      END SELECT
      PRINT *, 'Orographic gravity-wave test case :', TRIM(etat0_type)
   END SUBROUTINE getin_etat0

   SUBROUTINE compute_etat0(ngrid, lon, lat, phis, ps, Temp, ulon, ulat)
      USE disvert_mod
      USE omp_para
      INTEGER, INTENT(IN)    :: ngrid
      REAL(rstd), INTENT(IN)  :: lon(ngrid)
      REAL(rstd), INTENT(IN)  :: lat(ngrid)
      REAL(rstd), INTENT(OUT) :: ps(ngrid)
      REAL(rstd), INTENT(OUT) :: phis(ngrid)
      REAL(rstd), INTENT(OUT) :: Temp(ngrid, llm)
      REAL(rstd), INTENT(OUT) :: ulon(ngrid, llm)
      REAL(rstd), INTENT(OUT) :: ulat(ngrid, llm)
      INTEGER :: l, ij
      REAL(rstd) :: hyam, hybm

      ! Hexagons : ps,phis,temp
      DO l = ll_begin, ll_end
         ! The surface pressure is not set yet so we provide the hybrid coefficients
         hyam = .5*(ap(l) + ap(l + 1))/preff
         hybm = .5*(bp(l) + bp(l + 1))
         DO ij = 1, ngrid
            CALL comp_all(hyam, hybm, lon(ij), lat(ij), &
                          ps(ij), phis(ij), Temp(ij, l), ulon(ij, l), ulat(ij, l))
         END DO
      END DO
   END SUBROUTINE compute_etat0

   SUBROUTINE comp_all(hyam, hybm, lon, lat, psj, phisj, tempj, ulonj, ulatj)
      USE dcmip_initial_conditions_test_1_2_3
      REAL(rstd), INTENT(IN) :: hyam, hybm, lon, lat
      REAL(rstd), INTENT(OUT) :: psj, phisj, tempj, ulonj, ulatj
      REAL :: dummy
      dummy = 0.
      SELECT CASE (testcase)
      CASE (mountain)
         CALL test2_steady_state_mountain(lon, lat, dummy, dummy, 0, .TRUE., hyam, hybm, &
                                          ulonj, ulatj, dummy, tempj, phisj, psj, dummy, dummy, h0_test2)
      CASE (schaer_noshear)
         CALL test2_schaer_mountain(lon, lat, dummy, dummy, 0, .TRUE., hyam, hybm, 0, &
                                    ulonj, ulatj, dummy, tempj, phisj, psj, dummy, dummy, h0_test2, Teq_test2)
      CASE (schaer_shear)
         CALL test2_schaer_mountain(lon, lat, dummy, dummy, 0, .TRUE., hyam, hybm, 1, &
                                    ulonj, ulatj, dummy, tempj, phisj, psj, dummy, dummy, h0_test2, Teq_test2)
      END SELECT
   END SUBROUTINE comp_all

END MODULE etat0_dcmip2_mod
