MODULE etat0_dcmip2016_baroclinic_wave_mod
  USE icosa
  IMPLICIT NONE
  PRIVATE
  
  REAL(rstd),PARAMETER :: eta0=0.252
  REAL(rstd),PARAMETER :: etat=0.2
  REAL(rstd),PARAMETER :: ps0=1e5
  REAL(rstd),PARAMETER :: u0=35
  REAL(rstd),PARAMETER :: T0=288
  REAL(rstd),PARAMETER :: DeltaT=4.8e5
  REAL(rstd),PARAMETER :: Rd=287
  REAL(rstd),PARAMETER :: Gamma=0.005
  REAL(rstd),PARAMETER :: up0=1
  REAL(rstd),PARAMETER :: lonc=Pi/9, latc=2*Pi/9, latw=2*Pi/9
  REAL(rstd),PARAMETER :: pw=34000
  REAL(rstd),PARAMETER :: q0=0.021
  
  INTEGER,SAVE :: testcase
  !$OMP THREADPRIVATE(testcase)  
  
  PUBLIC getin_etat0, compute_etat0

CONTAINS

  SUBROUTINE getin_etat0
    USE mpipara, ONLY : is_mpi_root
    USE tracer_mod
    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
!    CALL set_advection_scheme(1,advect_none)
!    CALL set_advection_scheme(2,advect_none)
!    CALL set_advection_scheme(3,advect_none)
    
  END SUBROUTINE getin_etat0

  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
    USE icosa
    USE disvert_mod
    USE omp_para
    USE terminator
    USE tracer_mod
    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) :: 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) :: etal, etavl, etas, etavs, sinlat, coslat, &
         Y, Tave, T, phis_ave, vort, r2, utot, &
         dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
    
    etas=ap(1)/preff+bp(1)
    etavs=(etas-eta0)*Pi/2
    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
    
    DO ij=1,ngrid
       sinlat=SIN(lat(ij))
       coslat=COS(lat(ij))
       phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sinlat**6 * (coslat**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
            +(8./5*coslat**3 * (sinlat**2 + 2./3) - Pi/4)*radius*Omega )
       ps(ij)=ps0
    ENDDO
    
    DO l=ll_begin,ll_end
       etal = 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
       etavl=(etal-eta0)*Pi/2
       
       Tave=T0*etal**(Rd*Gamma/g)
       dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* etal**(Rd*Gamma/g-kappa-1)
       IF (etat>etal) THEN
          Tave=Tave+DeltaT*(etat-etal)**5
          dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-etal)**4 * etal**(-kappa)  &
               + kappa * (etat-etal)**5 * etal**(-kappa-1))
       END IF
       
       DO ij=1,ngrid
          sinlat=SIN(lat(ij))
          coslat=COS(lat(ij))
          
          K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
          r=radius*acos(K)
          utot=u0*cos(etavl)**1.5*sin(2*lat(ij))**2 + up0*exp(-(r/(0.1*radius))**2)
          ulon(ij,l) = utot
          ulat(ij,l) = 0.
          Y = ((-2*sinlat**6*(coslat**2+1./3)+10./63)*2*u0*cos(etavl)**1.5     &
               + (8./5*coslat**3*(sinlat**2+2./3)-Pi/4)*radius*Omega)
          T = Tave + 0.75*(etal*Pi*u0/Rd)*sin(etavl)*cos(etavl)**0.5 * Y
          temp(ij,l)=T
          
          q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2)
          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