MODULE etat0_dcmip5_mod
  USE icosa
  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

  PUBLIC getin_etat0, compute_etat0

CONTAINS
  
  SUBROUTINE getin_etat0
    USE mpipara, ONLY : is_mpi_root
    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
    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=1,llm
       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
          ENDIF
          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
    ENDDO
  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)
     ENDIF
     
     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
       ENDIF
       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
     ENDDO  
   
   END SUBROUTINE compute_z

END MODULE etat0_dcmip5_mod