MODULE etat0_bubble_mod
  USE icosa
  IMPLICIT NONE
  PRIVATE

  REAL(rstd), PARAMETER :: latc=0., lonc=0., &
       width        = 0.3, & ! width of mountain (radian)
       width_bubble = 0.5, & ! width of bubble (radian), curvature effects O(width^2)
       zc           = 350., & ! height at bubble center
       zw           = 250.    ! bubble vertical scale => radius=4000m for circular bubble 
  REAL(rstd), SAVE :: T0 ! adiabatic atmosphere at theta=T0 (default 300K)
  !$OMP THREADPRIVATE(T0)
  REAL(rstd), SAVE :: ztop ! mountain top (default 0.)
  !$OMP THREADPRIVATE(ztop)
  REAL(rstd), SAVE :: delta_T ! bubble extra temperature (default 5K)
  !$OMP THREADPRIVATE(delta_T)
  
  PUBLIC getin_etat0, compute_etat0

CONTAINS

  SUBROUTINE getin_etat0
    T0=300.
    CALL getin("bubble_T0",T0)
    ztop=0.
    CALL getin("bubble_mountain_top",ztop)
    delta_T=5.
    CALL getin("bubble_delta_T",delta_T)
  END SUBROUTINE getin_etat0

  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q)
    USE icosa
    USE disvert_mod
    USE omp_para
    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) :: ps(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) :: geopot(ngrid,llm+1)
    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
    
    INTEGER :: l,ij
    REAL(rstd) :: sinlat, coslat, K, Ts, p_ij, T_ij, z_ij
    ! adiabatic profile
    ! p=preff, T=T0 at Phi=0
    ! => theta = T0, cp.T + Phi = cst = cp.T0
    DO ij=1,ngrid
       sinlat=SIN(lat(ij))
       coslat=COS(lat(ij))
       K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
       K=acos(K)/width
       phis(ij) = g*ztop*exp(-(K**2))
       Ts = T0-phis(ij)/cpp ! adiabatic 
       ps(ij) = preff * (Ts/T0)**(1/kappa) ! theta=T0=T.(p/preff)^-kappa
    ENDDO
    
    DO l=1,llm+1
       DO ij=1,ngrid
          p_ij = ap(l)+bp(l)*ps(ij)
          T_ij = T0*(p_ij/preff)**kappa ! adiabatic 
          geopot(ij,l) = cpp*(T0-T_ij)
       END DO
    END DO
    
    DO l=1, llm
       DO ij=1,ngrid
          sinlat=SIN(lat(ij))
          coslat=COS(lat(ij))
          K=sin(latc)*sinlat+cos(latc)*coslat
          K=acos(K)/width_bubble
          p_ij = 0.5 *( ap(l)+bp(l)*ps(ij) + ap(l+1)+bp(l+1)*ps(ij) )
          T_ij = T0*(p_ij/preff)**kappa
          z_ij = (T0-T_ij)*cpp/g
          z_ij = (z_ij - zc)/zw
          K=MIN(1.,sqrt(K**2+z_ij**2))
          K=.5*(1.+COS(pi*K))
          Temp(ij,l) = T_ij + delta_T*K
          ulat(ij,l)=0.
          ulon(ij,l)=0.
          q(ij,l,:)=0.
       END DO
    END DO
    
  END SUBROUTINE compute_etat0
  
END MODULE etat0_bubble_mod