Index: /codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0.f90	(revision 36)
+++ /codes/icosagcm/trunk/src/etat0.f90	(revision 37)
@@ -9,4 +9,5 @@
   USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
   USE etat0_ncar_mod, ONLY : etat0_ncar=>etat0  
+  USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -28,4 +29,6 @@
       CASE ('ncar')
         CALL etat0_ncar(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+      CASE ('dcmip3')
+        CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
       CASE DEFAULT
         PRINT*, 'Bad selector for varaiable etat0 <',etat0_type,'> options are <jablonowsky06>, <academic>, <ncar> '
Index: /codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- /codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 37)
+++ /codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 37)
@@ -0,0 +1,188 @@
+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 
