MODULE etat0_dcmip2016_cyclone_mod
  USE icosa
  IMPLICIT NONE
  PRIVATE


  PUBLIC getin_etat0, compute_etat0

CONTAINS
  
  SUBROUTINE getin_etat0
    USE mpipara, ONLY : is_mpi_root

    IF(nqtot<5) THEN
       IF (is_mpi_root)  THEN
          PRINT *, "nqtot must be at least 5 for test case DCMIP2016 tropical cyclone"
       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
  USE dcmip2016_cyclone_mod, ONLY : tropical_cyclone_test
  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) :: 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 :: zcoords
    REAL :: p,z
    REAL :: rho, thetav
    INTEGER :: l,ij

   
    DO ij=1,ngrid
       z=0.
       zcoords=1                 
       CALL tropical_cyclone_test(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 tropical_cyclone_test(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_cyclone_mod