MODULE etat0_jablonowsky06_mod
  USE icosa
  PRIVATE
  REAL(rstd),PARAMETER :: eta0=0.252
  REAL(rstd),PARAMETER :: etat=0.2
  REAL(rstd),PARAMETER :: ps0=1e5
  REAL(rstd),PARAMETER :: u0=35
  REAL(rstd),PARAMETER :: T0=288
  REAL(rstd),PARAMETER :: DeltaT=4.8e5
  REAL(rstd),PARAMETER :: Rd=287
  REAL(rstd),PARAMETER :: Gamma=0.005
  REAL(rstd),PARAMETER :: up0=1

  PUBLIC compute_etat0

CONTAINS
  
  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, temp, ulon, ulat)
    USE disvert_mod
    IMPLICIT NONE  
    INTEGER, INTENT(IN) :: ngrid
    REAL(rstd),INTENT(IN) :: lat(ngrid)
    REAL(rstd),INTENT(IN) :: lon(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)
    
    INTEGER :: l,ij
    REAL(rstd) :: eta(llm)
    REAL(rstd) :: etav(llm)
    REAL(rstd) :: etas, etavs, Tave, phis_ave, r2
    
    DO l=1,llm
       eta(l)  = 0.5 *( ap(l)/ps0+bp(l) + ap(l+1)/ps0+bp(l+1) )
       etav(l) = (eta(l)-eta0)*Pi/2
    ENDDO
    etas  = ap(1)+bp(1)
    etavs = (etas-eta0)*Pi/2
        
    phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
    DO ij=1,ngrid
       ps(ij)=ps0
       phis(ij) = phis_ave + u0*cos(etavs)**1.5*( (-2*sin(lat(ij))**6 * (cos(lat(ij))**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
            +(8./5*cos(lat(ij))**3 * (sin(lat(ij))**2 + 2./3) - Pi/4)*radius*Omega )
    ENDDO

    DO l=1,llm
       Tave=T0*eta(l)**(Rd*Gamma/g)
       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
       DO ij=1,ngrid
          r2 = arc(Pi/9.,2.*Pi/9., lon(ij),lat(ij))**2
          temp(ij,l) = Tave + 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5  &
               * ( (-2*sin(lat(ij))**6*(cos(lat(ij))**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5 &
               + (8./5*cos(lat(ij))**3*(sin(lat(ij))**2+2./3)-Pi/4)*radius*Omega)
          ulon(ij,l) = u0*cos(etav(l))**1.5*sin(2*lat(ij))**2 + up0*exp(-r2/0.01)
          ulat(ij,l) = 0
       ENDDO
    ENDDO
    
    
  END SUBROUTINE compute_etat0
  
END MODULE etat0_jablonowsky06_mod