MODULE etat0_dcmip2016_supercell_mod
  USE icosa
  IMPLICIT NONE
  PRIVATE
  
  PUBLIC getin_etat0, compute_etat0

CONTAINS

  SUBROUTINE getin_etat0
    USE mpipara, ONLY : is_mpi_root
    USE tracer_icosa_mod
    USE dcmip2016_supercell_mod, ONLY : supercell_init, supercell_test
    IF(nqtot<5) THEN
       IF (is_mpi_root)  THEN
          PRINT *, "nqtot must be at least 5 for test case dcmip2016_supercell"
       END IF
       STOP
    END IF

    CALL supercell_init
    
  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_supercell_mod, ONLY : supercell_init, supercell_test, supercell_z
    USE terminator, ONLY: initial_value_Terminator
    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
    INTEGER,PARAMETER :: perturbation = 1
    
    REAL(rstd) :: p,z, thetav,rho
 
    phis(:)=0
    DO ij=1,ngrid
       z=0.
       CALL supercell_z(lon(ij), lat(ij), z, ps(ij), thetav, rho, q(ij,1,1), perturbation)
       DO l=ll_begin,ll_end
         p=0.5*(ap(l)+ap(l+1) + (bp(l)+bp(l+1)) * ps(ij))
         CALL supercell_test(lon(ij),lat(ij),p,z,0,ulon(ij,l),ulat(ij,l),temp(ij,l),thetav, ps(ij),rho, q(ij,l,1),perturbation) 
        
        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
    END DO
 
    
  END SUBROUTINE compute_etat0
  
END MODULE etat0_dcmip2016_supercell_mod