MODULE etat0_dcmip5_mod
   USE icosa, ONLY: rstd, pi, nqtot, getin
   IMPLICIT NONE
   PRIVATE

   REAL(rstd), PARAMETER :: zt = 15000
   REAL(rstd), PARAMETER :: q0 = 0.021
   REAL(rstd), PARAMETER :: qt = 1e-11
   REAL(rstd), PARAMETER :: T0 = 302.15
   REAL(rstd), PARAMETER :: Ts = 302.15
   REAL(rstd), PARAMETER :: zq1 = 3000
   REAL(rstd), PARAMETER :: zq2 = 8000
   REAL(rstd), PARAMETER :: Gamma = 0.007
   REAL(rstd), PARAMETER :: pb = 101500
   REAL(rstd), PARAMETER :: Deltap = 1115
   REAL(rstd), PARAMETER :: rp = 282000
   REAL(rstd), PARAMETER :: zp = 7000
   REAL(rstd), PARAMETER :: epsilon = 1e-25
   REAL(rstd), PARAMETER :: Rd = 287

   REAL(rstd), PARAMETER :: lonc = Pi, latc = Pi/18, &
                            Tv0 = T0*(1 + 0.608*q0), Tvt = Tv0 - Gamma*zt

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

   PUBLIC getin_etat0, compute_etat0

CONTAINS

   SUBROUTINE getin_etat0
      USE mpipara, ONLY: is_mpi_root
      ! @getin_doc dcmip5_testcase realm=initial type=integer default=1
      ! Active when `etat0=="dcmip5"`. Selects a variant of DCMIP 2012 test case 5
      ! (tropical cyclone). Valid values :
      !  * 1 : simple physics
      !  * 2 : full physics (not implemented)
      ! See DCMIP-TestCaseDocument_v1.6_23Jul2012.pdf pp 48-53.
      dcmip5_testcase = 1
      CALL getin("dcmip5_testcase", dcmip5_testcase)
      IF (nqtot < 1) THEN
         IF (is_mpi_root) THEN
            PRINT *, "nqtot must be at least 1 for test case DCMIP5"
         END IF
         STOP
      END IF
   END SUBROUTINE getin_etat0

   SUBROUTINE compute_etat0(ngrid, lon, lat, phis, ps, Temp, ulon, ulat, q)
      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)
      REAL(rstd), INTENT(OUT) :: q(ngrid, llm, nqtot)

      INTEGER :: l, ij
      REAL(rstd) :: r, d, d1, d2, zz, Tave, Tp, vt, aa, bb

      phis(:) = 0.

      DO ij = 1, ngrid
         r = radius*arc(lonc, latc, lon(ij), lat(ij))
         ps(ij) = pb - DeltaP*exp(-(r/rp)**1.5)
      END DO

      DO l = ll_begin, ll_end
         aa = .5*(ap(l) + ap(l + 1))
         bb = .5*(bp(l) + bp(l + 1))
         DO ij = 1, ngrid
            r = radius*arc(lonc, latc, lon(ij), lat(ij))
            CALL compute_z(aa + bb*ps(ij), ps(ij), r, zz)
            IF (zz <= zt) THEN
               q(ij, l, 1) = q0*exp(-zz/zq1)*exp(-(zz/zq2)**2)
               Tave = Tv0 - Gamma*zz
               Tp = (Tv0 - Gamma*zz)*((1 + 2*Rd*(Tv0 - Gamma*zz)*zz/ &
                                       (g*zp**2*(1 - pb/Deltap*exp((r/rp)**1.5 + (zz/zp)**2))))**(-1) - 1)
               vt = -omega*sin(latc)*r + sqrt((omega*sin(latc)*r)**2 - 1.5*(r/rp)**1.5*(Tv0 - Gamma*zz)*Rd &
                                              /(1 + 2*Rd*(Tv0 - Gamma*zz)*zz/(g*zp**2) &
                                                - pb/Deltap*exp((r/rp)**1.5 + (zz/zp)**2)))
            ELSE
               q(ij, l, 1) = qt
               Tave = Tvt
               Tp = 0
               vt = 0
            END IF
            Temp(ij, l) = Tave + Tp

            d1 = sin(latc)*cos(lat(ij)) - cos(latc)*sin(lat(ij))*cos(lon(ij) - lonc)
            d2 = cos(latc)*sin(lon(ij) - lonc)
            d = MAX(1e-25, sqrt(d1**2 + d2**2))
            ulon(ij, l) = vt*d1/d
            ulat(ij, l) = vt*d2/d
         END DO
      END DO
   END SUBROUTINE compute_etat0

   SUBROUTINE compute_z(pmodel, ps, r, z)
      USE icosa
      IMPLICIT NONE
      REAL(rstd), PARAMETER   :: epsilon0 = 2e-13
      REAL(rstd), INTENT(IN)  :: pmodel
      REAL(rstd), INTENT(IN)  :: ps
      REAL(rstd), INTENT(IN)  :: r
      REAL(rstd), INTENT(OUT) :: z

      REAL(rstd) :: pt, pave, pp, dpdz
      REAL(rstd) :: znp1
      REAL(rstd) :: epsilon
      REAL(rstd) :: p
      INTEGER  :: n

      pt = pb*(Tvt/Tv0)**(g/(Rd*Gamma))

      IF (pmodel > pt) THEN
         z = Tv0/Gamma*(1 - (pmodel/ps)**(Rd*Gamma/g))
      ELSE
         z = zt + Rd*Tvt/g*log(Pt/pmodel)
      END IF

      IF (z > zt .OR. r > 1000000) return

      epsilon = 1
      n = 0

      DO WHILE (epsilon > epsilon0 .AND. n < 20)

         IF (z < zt) THEN
            pave = pb*((Tv0 - Gamma*z)/Tv0)**(g/(Rd*Gamma))
            pp = -Deltap*exp(-(r/rp)**1.5 - (z/zp)**2)*((Tv0 - Gamma*z)/Tv0)**(g/(Rd*Gamma))
         ELSE
            pave = pt*exp(g*(zt - z)/(Rd*Tvt))
            pp = 0
         END IF
         p = pave + pp

         dpdz = 2*Deltap*z/zp**2*exp(-(r/rp)**1.5)*exp(-(z/zp)**2)*((Tv0 - Gamma*z)/Tv0)**(g/(Rd*Gamma)) &
                - g/(Rd*Tv0)*(preff - Deltap*exp(-(r/rp)**1.5)*exp(-(z/zp)**2))*((Tv0 - Gamma*z)/Tv0)**(g/(Rd*Gamma) - 1)

         znp1 = z - (pmodel - p)/(-dpdz)

         epsilon = ABS((znp1 - z)/znp1)
         z = znp1
         n = n + 1
      END DO

   END SUBROUTINE compute_z

END MODULE etat0_dcmip5_mod
