MODULE etat0_dcmip5_mod
  USE icosa
  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) :: lonc
  REAL(rstd) :: latc
  TYPE(t_field),POINTER :: f_out_i(:)
  REAL(rstd),POINTER :: out_i(:,:)

  PUBLIC  etat0
CONTAINS
  
   
    
  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
  
    REAL(rstd),POINTER :: ps(:)
    REAL(rstd),POINTER :: phis(:)
    REAL(rstd),POINTER :: theta_rhodz(:,:)
    REAL(rstd),POINTER :: u(:,:)
    REAL(rstd),POINTER :: q(:,:,:)
    INTEGER :: ind
 
    CALL allocate_field(f_out_i,field_t,type_real,llm)
 
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      phis=f_phis(ind)
      theta_rhodz=f_theta_rhodz(ind)
      u=f_u(ind)
      q=f_q(ind)
      out_i=f_out_i(ind)
      CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q)
    ENDDO
    CALL writefield("out_i",f_out_i)
    CALL deallocate_field(f_out_i)
  END SUBROUTINE etat0
  
  SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q)
  USE icosa
  USE pression_mod
  USE wind_mod
  USE theta2theta_rhodz_mod
  USE exner_mod
  IMPLICIT NONE  
  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
  REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm,llm)
  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
  
  INTEGER :: i,j,l,ij
  INTEGER :: testcase
  REAL(rstd) :: Tv0, Tvt 
  REAL(rstd) :: r 
  REAL(rstd) :: z(iim*jjm,llm),zz 
  REAL(rstd) :: p(iim*jjm,llm+1) 
  REAL(rstd) :: Tave,Tp 
  REAL(rstd) :: T(iim*jjm,llm) 
  REAL(rstd) :: vt 
  REAL(rstd) :: d1,d2,d 
  REAL(rstd) :: ulon(3*iim*jjm,llm) 
  REAL(rstd) :: ulat(3*iim*jjm,llm)
  REAL(rstd) :: lon,lat
  REAL(rstd) :: pks(iim*jjm) 
  REAL(rstd) :: pk(iim*jjm,llm) 


    latc=Pi/18
    lonc=Pi 

    Tv0=T0/(1+0.608*q0)
    Tvt=Tv0-Gamma*zt
  
    testcase=1
    CALL getin("dcmip5_testcase",testcase)
    
    phis(:)=0.
   
    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
        r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
        ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
      ENDDO
    ENDDO
    
    CALL compute_pression(ps,p,0)

    DO l=1,llm
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l))
        ENDDO
      ENDDO
    ENDDO
    out_i=z
    
    DO l=1,llm
      DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
          IF (z(ij,l)<=zt) THEN
            q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2)
          ELSE
            q(ij,l,1)=qt
          ENDIF
        ENDDO
      ENDDO
    ENDDO

    
     DO l=1,llm
      DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          IF (z(ij,l)<=zt) THEN
            Tave=Tv0-Gamma*z(ij,l)
            Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) /                           &
                                            (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) 
          ELSE
            Tave=Tvt
            Tp=0
          ENDIF
          T(ij,l)=Tave+Tp 
          out_i(ij,l)=Tave+Tp
        ENDDO
      ENDDO
    ENDDO

    DO l=1,llm
      DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i

          CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
          zz=0.5*(z(ij,l)+z(ij+t_right,l))
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          IF (zz<=zt) THEN
            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
            vt = 0
          ENDIF
          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
          d2=cos(latc)*sin(lon-lonc)
          d=MAX(1e-25,sqrt(d1**2+d2**2))
          ulon(ij+u_right,l)=vt*d1/d
          ulat(ij+u_right,l)=vt*d2/d


          
          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
          zz=0.5*(z(ij,l)+z(ij+t_lup,l))
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          IF (zz<=zt) THEN
            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
            vt = 0
          ENDIF
          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
          d2=cos(latc)*sin(lon-lonc)
          d=MAX(1e-25,sqrt(d1**2+d2**2))
          ulon(ij+u_lup,l)=vt*d1/d
          ulat(ij+u_lup,l)=vt*d2/d
          
          
          
          
          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
          zz=0.5*(z(ij,l)+z(ij+t_ldown,l))
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          IF (zz<=zt) THEN
            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
            vt = 0
          ENDIF
          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
          d2=cos(latc)*sin(lon-lonc)
          d=MAX(1e-25,sqrt(d1**2+d2**2))
          ulon(ij+u_ldown,l)=vt*d1/d
          ulat(ij+u_ldown,l)=vt*d2/d
          


          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
          zz=z(ij,l)
          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
          IF (zz<=zt) THEN
            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
            vt = 0
          ENDIF
          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
          d2=cos(latc)*sin(lon-lonc)
          d=MAX(1e-25,sqrt(d1**2+d2**2))
!          ulon(ij+u_ldown,l)=vt*d1/d
!          ulat(ij+u_ldown,l)=vt*d2/d
!          out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2)
        ENDDO
      ENDDO
    ENDDO
    
    CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue)
    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)

  END SUBROUTINE compute_etat0_dcmip5
  
  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) :: Tv0,Tvt,pt
     REAL(rstd) :: pave, pp, dpdz
     REAL(rstd) :: znp1
     REAL(rstd) :: epsilon
     REAL(rstd) :: p
     INTEGER  :: n
          
     Tv0=T0/(1+0.608*q0)
     Tvt=Tv0-Gamma*zt
     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
     
!     PRINT*,pmodel,pt,z,zt
     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 )
!       PRINT *,"epsilon",epsilon,r,z,znp1
!       PRINT *,"----",n,pmodel,ps,pave,dpdz
       z=znp1
       n=n+1
     ENDDO  
   
   END SUBROUTINE compute_z

END MODULE etat0_dcmip5_mod