MODULE etat0_dcmip3_mod

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

! Questions
! Replace ps0 by preff ??

  USE genmod
  USE dcmip_initial_conditions_test_1_2_3

  PRIVATE

  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
      IF (.NOT. assigned_domain(ind)) CYCLE
      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),PARAMETER :: u0=20.         ! Maximum amplitude of the zonal wind (m.s-1)
  REAL(rstd),PARAMETER :: N=0.01         ! Brunt-Vaisala 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

  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) :: T(iim*jjm,llm)
  REAL(rstd) :: p(iim*jjm,llm+1)
  REAL(rstd) :: theta(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
  REAL(rstd) :: C0, C1, GG
  REAL(rstd) :: p0psk, pspsk,r,zz, thetab, thetap
  REAL(rstd) :: dummy, pp
  LOGICAL    :: use_dcmip_routine
  
  Rd=cpp*kappa
  
  GG=(g/N)**2/cpp
  C0=0.25*u0*(u0+2.*Omega*radius)
  
  q(:,:,:)=0
  
!  use_dcmip_routine=.TRUE.
  use_dcmip_routine=.FALSE.
  dummy=0.

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

        IF(use_dcmip_routine) THEN
           CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
        ELSE
           C1=C0*(cos(2*lat_i(ij))-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_i(ij))+cos(latc)*cos(lat_i(ij))*cos(lon_i(ij)-lonc))
           s(ij)= d**2/(d**2+r**2)
        END IF
     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
           pp=0.5*(p(ij,l+1)+p(ij,l))  ! full-layer pressure
           IF(use_dcmip_routine) THEN
              CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, &
                   dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy)
           ELSE
              pspsk=(pp/ps(ij))**kappa
              p0psk=(Peq/ps(ij))**kappa
              thetab = Ts(ij)*p0psk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background pot. temp.
              zz     = -g/N**2*log( Ts(ij)/GG * (pspsk -1)+1)         ! altitude
              thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij)                ! perturbation pot. temp.
              theta(ij,l) = thetab + thetap
              T(ij,l) = theta(ij,l)* ((pp/peq)**kappa)
              ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background temp.
           END IF
        ENDDO
     ENDDO
  ENDDO
  
  IF(use_dcmip_routine) THEN
     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
  ELSE
     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
  END IF
  
  pp=peq
  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(use_dcmip_routine) THEN
              CALL test3_gravity_wave(lon_e(ij+u_right),lat_e(ij+u_right), &
                   pp,0.,0, ulon(ij+u_right,l),ulat(ij+u_right,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
              CALL test3_gravity_wave(lon_e(ij+u_lup),lat_e(ij+u_lup), &
                   pp,0.,0, ulon(ij+u_lup,l),ulat(ij+u_lup,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
              CALL test3_gravity_wave(lon_e(ij+u_ldown),lat_e(ij+u_ldown), &
                   pp,0.,0, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l),&
                   dummy,dummy,dummy,dummy,dummy,dummy)
           ELSE
              ulon(ij+u_right,l) = u0*cos(lat_e(ij+u_right))
              ulat(ij+u_right,l) = 0
              ulon(ij+u_lup,l)   = u0*cos(lat_e(ij+u_lup))
              ulat(ij+u_lup,l)   = 0
              ulon(ij+u_ldown,l) = u0*cos(lat_e(ij+u_ldown))
              ulat(ij+u_ldown,l) = 0
           END IF
        ENDDO
     ENDDO
  ENDDO
  
  CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)    
  
END SUBROUTINE compute_etat0_DCMIP3


END MODULE etat0_DCMIP3_mod