MODULE etat0_williamson_mod
  USE icosa
  PRIVATE
  REAL(rstd), PARAMETER :: h0=8.E3
  REAL(rstd), PARAMETER :: R0=4
  REAL(rstd), PARAMETER :: K0=7.848E-6
  REAL(rstd), PARAMETER :: omega0=K0
  
  PUBLIC getin_etat0, compute_etat0
  
CONTAINS

  SUBROUTINE getin_etat0
    USE mpipara, ONLY : is_mpi_root
    USE disvert_mod, ONLY : caldyn_eta, eta_lag
    IF(caldyn_eta /= eta_lag) THEN
       IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with caldyn_eta=eta_lag'
       STOP
    END IF

    IF(llm>1) THEN
       IF(is_mpi_root) PRINT *, 'etat0_type=williamson91.6 (Williamson,1991) must be used with llm=1 but llm =',llm
       STOP
    END IF
  END SUBROUTINE getin_etat0

  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,mass,rhodz,ulon,ulat)
    IMPLICIT NONE  
    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) :: mass(ngrid)
    REAL(rstd),INTENT(OUT) :: rhodz(ngrid)
    REAL(rstd),INTENT(OUT) :: ulon(ngrid)
    REAL(rstd),INTENT(OUT) :: ulat(ngrid)

    REAL(rstd) :: coslat,sinlat, A,B,C    
    INTEGER :: ij
    
    DO ij=1,ngrid
       coslat=cos(lat(ij))
       sinlat=sin(lat(ij))
       A = 0.5*omega0*(2*omega+omega0)*coslat**2   &
            + 0.25*K0**2*coslat**(2*R0)*((R0+1)*coslat**2+(2*R0**2-R0-2)-2*R0**2/coslat**2)
       B = 2*(omega+omega0)*K0/((R0+1)*(R0+2))*coslat**R0*((R0**2+2*R0+2)-(R0+1)**2*coslat**2)
       C = 0.25*K0**2*coslat**(2*R0)*((R0+1)*coslat**2-(R0+2))
       
       mass(ij) = (g*h0+radius**2*A+radius**2*B*cos(R0*lon(ij))+radius**2*C*cos(2*R0*lon(ij)))/g
       ulon(ij) = radius*omega0*coslat+radius*K0*coslat**(R0-1)*(R0*sinlat**2-coslat**2)*cos(R0*lon(ij))
       ulat(ij) = -radius*K0*R0*coslat**(R0-1)*sinlat*sin(R0*lon(ij))
    ENDDO

    phis(:)=0.
    rhodz(:)=mass(:)

  END SUBROUTINE compute_etat0
  
END MODULE etat0_williamson_mod