MODULE etat0_dcmip2016_baroclinic_wave_mod
  USE icosa
  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_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
    ENDIF
    
  
  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 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) :: 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 :: deep=0
    INTEGER :: zcoords
    REAL :: p,z
    REAL :: rho, thetav
    INTEGER :: ij,l
    INTEGER :: moist
    
    moist=0
    IF (physics_thermo==thermo_moist .OR. physics_thermo==thermo_fake_moist) moist=1
   
    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(ij,1),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(ij,l),thetav,phis(ij),ps(ij),rho,q(ij,l,1))
                                   
        IF (physics_thermo==thermo_fake_moist) temp(ij,l)=Temp(ij,l)*(1+0.608*q(ij,l,1))
        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
    ENDDO

  END SUBROUTINE compute_etat0

END MODULE etat0_dcmip2016_baroclinic_wave_mod