Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 344)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 345)
@@ -16,5 +16,5 @@
     USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
-    ! New interface
+    ! Generic interface
     USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0
     USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0
@@ -22,7 +22,6 @@
     USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
     USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
-    ! Old interface
+    ! Ad hoc interfaces
     USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
-    USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
     USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
@@ -53,5 +52,5 @@
     CALL getin("etat0",etat0_type)
     
-    !------------------- New interface ---------------------
+    !------------------- Generic interface ---------------------
     collocated=.TRUE.
     SELECT CASE (TRIM(etat0_type))
@@ -65,4 +64,5 @@
     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
        CALL getin_etat0_dcmip2
+    CASE ('dcmip3')
     CASE ('dcmip5')
         CALL getin_etat0_dcmip5
@@ -74,5 +74,5 @@
     END SELECT
 
-    !------------------- Old interface --------------------
+    !------------------- Ad hoc interfaces --------------------
     SELECT CASE (TRIM(etat0_type))
     CASE ('start_file')
@@ -86,6 +86,4 @@
        CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
        PRINT *, "Venus (Lebonnois et al., 2012) test case"
-    CASE ('dcmip3')
-       CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('dcmip4')
         IF(nqtot<2) THEN
@@ -171,4 +169,5 @@
     USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
     USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
+    USE etat0_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
     USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
     USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
@@ -211,4 +210,7 @@
        CALL compute_dcmip2(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
        CALL compute_dcmip2(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)      
+    CASE('dcmip3')
+       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
+       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
     CASE('dcmip5')
        CALL compute_dcmip5(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 344)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 345)
@@ -1,197 +1,40 @@
 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
-
+  ! test cases DCMIP 2012, category 3 : Non-hydrostatic gravity waves
+  IMPLICIT NONE
   PRIVATE
-
-  PUBLIC  etat0
-
+  PUBLIC :: compute_etat0
+  
 CONTAINS
-
-
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  USE theta2theta_rhodz_mod
-  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(:)
-    TYPE(t_field),POINTER,SAVE :: f_temp(:)
-    
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: Temp(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-
-    INTEGER :: ind
-
-    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 
-
-    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)
-      q=f_q(ind)
-      temp=f_temp(ind)
-      CALL compute_etat0_DCMIP3(ps,phis,u,Temp,q)
-    ENDDO
-
-    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
-    CALL deallocate_field(f_temp)
-            
-  END SUBROUTINE etat0
   
-
-  SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, temp,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) :: Temp(iim*jjm,llm)
-  REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
-  
-  REAL(rstd) :: Ts(iim*jjm)
-  REAL(rstd) :: s(iim*jjm)
-  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
-  
-!$OMP BARRIER
-  CALL compute_pression(ps,p,0)
-!$OMP BARRIER
-  
-  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,Temp(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
-              Temp(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
-
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
+    USE genmod
+    USE dcmip_initial_conditions_test_1_2_3
+    USE disvert_mod
+    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) :: q(ngrid,llm,nqtot)
+    REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
+    REAL(rstd) :: dummy, pp
+    INTEGER :: l,ij
+    pp=peq
+    DO ij=1,ngrid
+       CALL test3_gravity_wave(lon(ij),lat(ij),pp,dummy,0, &
+            dummy,dummy,dummy,dummy,phis(ij),ps(ij),dummy,dummy)
+    END DO
+    DO l=1,llm
+       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(lon(ij),lat(ij),pp,dummy,0, &
+               ulon(ij,l),ulat(ij,l),dummy,Temp(ij,l),dummy,dummy,dummy,dummy)
+       END DO
+    END DO    
+    q=0.
+  END SUBROUTINE compute_etat0
 
 END MODULE etat0_DCMIP3_mod
