MODULE etat0_dcmip3_mod

! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves

! Questions
! Replace ps0 by preff ??

  USE genmod
  PRIVATE
! some of the following should be obtained via getin
!  REAL(rstd),PARAMETER :: ztop=10000.   ! Height position of the model top (m)
!  REAL(rstd),PARAMETER :: ptop=27391.9  ! Pressure at the model top (Pa)
!  REAL(rstd),PARAMETER :: Xscale=1.     ! Small-planet scaling factor (1=Earth)
!  REAL(rstd),PARAMETER :: Omega =0      ! Planetary rotation rate (s-1)
  REAL(rstd),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1)
  REAL(rstd),PARAMETER :: zs=0.          ! Surface elevation (m)
  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaissala frequency (s-1)
  REAL(rstd),PARAMETER :: Teq=300.       ! Surface temperature at the equator (K)
  REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
  REAL(rstd),PARAMETER :: d=5000.        ! Witdth parameter for theta
  REAL(rstd),PARAMETER :: lonc=2*pi/3    ! Longitudinal centerpoint of theta
  REAL(rstd),PARAMETER :: latc=0         ! Longitudinal centerpoint of theta
  REAL(rstd),PARAMETER :: dtheta=1.      ! Maximum amplitude of theta (K)
  REAL(rstd),PARAMETER :: Lz=20000.      ! Vertical wave lenght of the theta perturbation

  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 :: u(:,:)
    REAL(rstd),POINTER :: theta_rhodz(:,:)
    REAL(rstd),POINTER :: q(:,:,:)

    INTEGER :: ind

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)

      ps=f_ps(ind)
      phis=f_phis(ind)
      u=f_u(ind)
      theta_rhodz=f_theta_rhodz(ind)
      q=f_q(ind)
      CALL compute_etat0_DCMIP3(ps,phis,u,theta_rhodz,q)
    ENDDO
            
  END SUBROUTINE etat0
  

  SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, theta_rhodz,q)
  USE icosa
  USE pression_mod
  USE theta2theta_rhodz_mod
  USE wind_mod
  IMPLICIT NONE
  REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
  REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
  REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm)
  REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
  REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
  
  REAL(rstd) :: Ts(iim*jjm)
  REAL(rstd) :: s(iim*jjm)
  REAL(rstd) :: z(iim*jjm,llm)
  REAL(rstd) :: T(iim*jjm,llm)
  REAL(rstd) :: p(iim*jjm,llm+1)
  REAL(rstd) :: thetap(iim*jjm,llm)
  REAL(rstd) :: thetap_rhodz(iim*jjm,llm)
  REAL(rstd) :: ulon(3*iim*jjm,llm)
  REAL(rstd) :: ulat(3*iim*jjm,llm)
  
  
  INTEGER :: i,j,l,ij
  REAL(rstd) :: Rd        ! gas constant of dry air, P=rho.Rd.T
  REAL(rstd) :: alpha, rm, zs
  REAL(rstd) :: lon,lat
  REAL(rstd) :: C0,C1
  REAL(rstd) :: GG
  REAL(rstd) :: pspsk,r
  
  Rd=cpp/kappa
    
!    alpha=g/Rd/lapse_rate              ! g/(Rd*Gamma)

    GG=(g/N)**2/cpp
    C0=0.25*u0*(u0+2.*Omega*radius)

    q(:,:,:)=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)
        C1=C0*(cos(2*lat)-1)
        
        !--- GROUND GEOPOTENTIAL
        phis(ij)=0.

        !--- GROUND TEMPERATURE
        Ts(ij) = GG+(Teq-GG)*EXP(-C1*(N/g)**2)

        !--- GROUND PRESSURE
        Ps(ij) = peq*EXP(C1/GG/Rd)*(Ts(ij)/Teq)**(1/kappa)


        r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
        s(ij)= d**2/(d**2+r**2)
      END DO
    END DO

    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

          pspsk=(0.5*(p(ij,l+1)+p(ij,l))/ps(ij) )**kappa
          T(ij,l)=Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1)
          z(ij,l)=-g/N**2*log( Ts(ij)/GG * (pspsk -1)+1) 
        ENDDO
      ENDDO
    ENDDO
    
    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)

    DO l=1,llm
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          thetap(ij,l)=dtheta *sin(2*Pi*z(ij,l)/Lz) * s(ij)
        ENDDO
      ENDDO
    ENDDO  
    
    CALL compute_theta2theta_rhodz(ps,thetap,thetap_rhodz,0)
    
    
    DO l=1,llm
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          theta_rhodz(ij,l)=theta_rhodz(ij,l)+thetap_rhodz(ij,l)
        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)
          ulon(ij+u_right,l)=u0*cos(lat)
          ulat(ij+u_right,l)=0

          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
          ulon(ij+u_lup,l)=u0*cos(lat)
          ulat(ij+u_lup,l)=0

          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
          ulon(ij+u_ldown,l)=u0*cos(lat)
          ulat(ij+u_ldown,l)=0
        ENDDO
      ENDDO
    ENDDO
    
 
    CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)

    
         
  END SUBROUTINE compute_etat0_DCMIP3


END MODULE etat0_DCMIP3_mod