source: codes/icosagcm/trunk/src/initial/etat0_dcmip3.f90 @ 548

MODULE etat0_dcmip3_mod
  ! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: compute_etat0
  
CONTAINS
  
  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q)
    USE dcmip_initial_conditions_test_1_2_3
    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) :: ulon(ngrid,llm)
    REAL(rstd), INTENT(OUT) :: ulat(ngrid,llm)
    REAL(rstd), INTENT(OUT) :: temp(ngrid,llm)
    REAL(rstd), INTENT(OUT) :: geopot(ngrid,llm+1)
    REAL(rstd), INTENT(OUT) :: q(ngrid,llm,nqtot)
    REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
    REAL(rstd) :: dummy, pp, zz
    INTEGER :: l,ij
    pp=peq
    DO ij=1,ngrid
       CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, &
            dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
    END DO
    DO l=ll_begin,ll_endp1
       DO ij=1,ngrid
          pp = ap(l) + bp(l)*ps(ij) ! half-layer pressure
          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,zz,0, &
               dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy)
          geopot(ij,l) = g*zz ! initialize geopotential for NH
       END DO
    END DO
    DO l=ll_begin,ll_end
       DO ij=1,ngrid
          pp = .5*(ap(l)+ap(l+1)) + .5*(bp(l)+bp(l+1))*ps(ij) ! full-layer pressure
          CALL test3_gravity_wave(scale_factor, lon(ij),lat(ij),pp,dummy,0, &
               ulon(ij,l),ulat(ij,l),dummy,Temp(ij,l),dummy,dummy,dummy,dummy)
       END DO
       q(:,l,:)=0.
    END DO
    
  END SUBROUTINE compute_etat0

END MODULE etat0_DCMIP3_mod