Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 343)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 344)
@@ -1,4 +1,5 @@
 MODULE etat0_mod
   USE icosa
+  IMPLICIT NONE 	
   PRIVATE
 
@@ -16,4 +17,6 @@
     USE disvert_mod
     ! New interface
+    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0
+    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0
     USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
     USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
@@ -21,6 +24,4 @@
     ! Old interface
     USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
-    USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0
-    USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0
     USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
@@ -38,5 +39,5 @@
     
     REAL(rstd),POINTER :: ps(:), mass(:,:)
-    LOGICAL :: init_mass
+    LOGICAL :: init_mass, collocated
     INTEGER :: ind,i,j,ij,l
 
@@ -52,22 +53,27 @@
     CALL getin("etat0",etat0_type)
     
+    !------------------- New interface ---------------------
+    collocated=.TRUE.
     SELECT CASE (TRIM(etat0_type))
-       !------------------- New interface ---------------------
     CASE ('isothermal')
        CALL getin_etat0_isothermal
-       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('temperature_profile')
        CALL getin_etat0_temperature
-       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('jablonowsky06')
-       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip5')
+    CASE ('dcmip1')
+        CALL getin_etat0_dcmip1
+    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
+       CALL getin_etat0_dcmip2
+    CASE ('dcmip5')
         CALL getin_etat0_dcmip5
-        CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
     CASE ('williamson91.6')
        init_mass=.FALSE.
        CALL getin_etat0_williamson
-       CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
-       !------------------- Old interface --------------------
+    CASE DEFAULT
+       collocated=.FALSE.
+    END SELECT
+
+    !------------------- Old interface --------------------
+    SELECT CASE (TRIM(etat0_type))
     CASE ('start_file')
        CALL etat0_start_file(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
@@ -80,8 +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 ('dcmip1')
-       CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
-       CALL etat0_dcmip2(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)
@@ -95,7 +97,11 @@
         CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE DEFAULT
-       PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
-            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
-       STOP
+      IF(collocated) THEN
+         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+      ELSE
+         PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
+              '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
+         STOP
+      END IF
     END SELECT
 
@@ -163,4 +169,6 @@
     USE wind_mod
     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
+    USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
+    USE etat0_dcmip2_mod, ONLY : compute_dcmip2 => compute_etat0
     USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
     USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
@@ -197,4 +205,10 @@
        CALL compute_jablonowsky06(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i)
        CALL compute_jablonowsky06(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e)
+    CASE('dcmip1')
+       CALL compute_dcmip1(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
+       CALL compute_dcmip1(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
+    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
+       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('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_dcmip1.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 343)
+++ codes/icosagcm/trunk/src/etat0_dcmip1.f90	(revision 344)
@@ -1,5 +1,12 @@
 MODULE etat0_dcmip1_mod
   USE icosa
+  IMPLICIT NONE 	
   PRIVATE
+
+  INTEGER, PARAMETER :: const=1, cos_bell=2, slotted_cyl=3, &
+       dbl_cos_bell_q1=4, dbl_cos_bell_q2=5, complement=6, hadley=7, &
+       dcmip11=-1
+  INTEGER, SAVE :: shape
+!$OMP THREADPRIVATE(shape)
 
   REAL(rstd), SAVE  :: h0=1.
@@ -36,103 +43,76 @@
 !$OMP THREADPRIVATE(zc)
 
-  PUBLIC etat0
+  PUBLIC getin_etat0, compute_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 :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
+
+  SUBROUTINE getin_etat0
     CHARACTER(len=255) :: dcmip1_adv_shape
-    INTEGER :: ind
-    
     R0=radius*0.5
     rt=radius*0.5
     dcmip1_adv_shape='cos_bell'
     CALL getin('dcmip1_shape',dcmip1_adv_shape)
-    
-    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)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      q=f_q(ind)
-
-
-      SELECT CASE(TRIM(dcmip1_adv_shape))
-      CASE('const')
-         CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('cos_bell')
-         CALL compute_etat0_ncar(2,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('slotted_cyl')
-         CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('dbl_cos_bell_q1')
-         CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('dbl_cos_bell_q2')
-         CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('complement')
-         CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('hadley')  ! hadley like meridional circulation 
-         CALL compute_etat0_ncar(7,ps, phis, theta_rhodz, u, q(:,:,1))
-      CASE('dcmip11')
-         IF(nqtot==5) THEN
-            CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1))
-            CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,2))
-            CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,3))
-            CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,4))
-            CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,5))
-         ELSE
-            PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .'
-            PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11'
-            STOP
-         END IF
-      CASE DEFAULT
-         PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   &
-              '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
-              '<dbl_cos_bell_q2>, <complement>, <hadley>'
-         STOP
-      END SELECT
-
-    ENDDO
-
-  END SUBROUTINE etat0
-  
-  SUBROUTINE compute_etat0_ncar(icase, ps, phis, theta_rhodz, u, q)
+    SELECT CASE(TRIM(dcmip1_adv_shape))
+    CASE('const')
+       shape=const
+    CASE('cos_bell')
+       shape=cos_bell
+    CASE('slotted_cyl')
+       shape=slotted_cyl
+    CASE('dbl_cos_bell_q1')
+       shape=dbl_cos_bell_q1
+    CASE('dbl_cos_bell_q2')
+       shape=dbl_cos_bell_q2
+    CASE('complement')
+       shape=complement
+    CASE('hadley')  ! hadley like meridional circulation 
+       shape=hadley
+    CASE('dcmip11')
+       IF(nqtot<5) THEN
+          PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .'
+          PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11'
+          STOP
+       END IF
+       shape=dcmip11
+    CASE DEFAULT
+       PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   &
+            '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
+            '<dbl_cos_bell_q2>, <complement>, <hadley>'
+       STOP
+    END SELECT
+  END SUBROUTINE getin_etat0
+
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
+    INTEGER, INTENT(IN) :: ngrid
+    REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid),phis(ngrid)
+    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm),ulon(ngrid,llm),ulat(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
+    ps = ncar_p0
+    phis = 0. ; temp = 0. ; 
+    ulon = 0. ; ulat=0.
+    SELECT CASE(shape)
+    CASE(dcmip11)
+       CALL compute_etat0_ncar(4,ngrid,lon,lat,q(:,:,1))
+       CALL compute_etat0_ncar(5,ngrid,lon,lat,q(:,:,2))
+       CALL compute_etat0_ncar(3,ngrid,lon,lat,q(:,:,3))
+       CALL compute_etat0_ncar(6,ngrid,lon,lat,q(:,:,4))
+       CALL compute_etat0_ncar(1,ngrid,lon,lat,q(:,:,5))
+    CASE DEFAULT
+       CALL compute_etat0_ncar(shape,ngrid,lon,lat,q(:,:,1))
+    END SELECT
+  END SUBROUTINE compute_etat0
+
+  SUBROUTINE compute_etat0_ncar(icase,ngrid,lon,lat, q)
   USE icosa
   USE disvert_mod
-  USE pression_mod
-  USE exner_mod
-  USE geopotential_mod
-  USE theta2theta_rhodz_mod
   IMPLICIT NONE  
-  INTEGER, INTENT(in) :: icase
-  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm)  
-
-  REAL(rstd) :: qxt1(iim*jjm,llm)
-  REAL(rstd) :: lon, lat
-  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
-  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
-  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
-  REAL(rstd) :: X2(3),X1(3)
-  INTEGER :: i,j,n,l
-
-  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
+  INTEGER, INTENT(IN) :: icase, ngrid
+  REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid)  
+  REAL(rstd),INTENT(OUT) :: q(ngrid,llm)  
+  REAL(rstd) :: zr(llm+1), zrl(llm), qxt1(ngrid,llm)
+  REAL(rstd) :: pr
+  !  REAL(rstd) :: lon, lat
+  INTEGER :: n,l
   
   DO l=1, llm+1
@@ -149,10 +129,8 @@
      q=1
   CASE(2)
-     !--------------------------------------------- SINGLE COSINE BELL
      CALL cosine_bell_1(q)      
   CASE(3)
      CALL slotted_cylinders(q) 	
   CASE(4)
-     PRINT *, 'Double cosine bell'
      CALL cosine_bell_2(q)      
   CASE(5)
@@ -176,109 +154,94 @@
   END SELECT
       
-CONTAINS
-
+  CONTAINS
 !======================================================================
 
     SUBROUTINE cosine_bell_1(hx)
     IMPLICIT NONE 
-    REAL(rstd) :: hx(iim*jjm,llm)
+    REAL(rstd) :: hx(ngrid,llm)
+    REAL(rstd) :: rr1,rr2   
+    INTEGER :: n,l
+    DO l=1,llm 
+       DO n=1,ngrid
+          CALL dist_lonlat(lon0,lat0,lon(n),lat(n),rr1)  ! GC distance from center 
+          rr1 = radius*rr1
+          IF ( rr1 .LT. R0 ) then 
+             hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
+          ELSE
+             hx(n,l)=0.0
+          END IF
+       END DO
+    END DO
+  END SUBROUTINE cosine_bell_1
+
+!==============================================================================
+  
+  SUBROUTINE cosine_bell_2(hx) 
+    REAL(rstd) :: hx(ngrid,llm)
+    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1
+    INTEGER :: n,l
+    DO l=1,llm
+       DO n=1,ngrid
+          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center 
+          rr1 = radius*rr1
+          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center 
+          rr2 = radius*rr2
+          dd1t1 = rr1/rt
+          dd1t1 = dd1t1*dd1t1
+          dd1t2 = (zrl(l) - zc)/zt
+          dd1t2 = dd1t2*dd1t2 
+          dd1 = dd1t1 + dd1t2 
+          dd1 = Min(1.0,dd1)  
+          dd2t1 = rr2/rt
+          dd2t1 = dd2t1*dd2t1
+          dd2 = dd2t1 + dd1t2 
+          dd2 = Min(1.0,dd2) 
+          hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
+       END DO
+    END DO  
+  END SUBROUTINE cosine_bell_2
+  
+  !=============================================================================
+
+  SUBROUTINE slotted_cylinders(hx) 
+    REAL(rstd) :: hx(ngrid,llm)
+    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1
+    INTEGER :: n,l
+    DO l=1,llm
+       DO n=1,ngrid
+          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center 
+          rr1 = radius*rr1
+          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center 
+          rr2 = radius*rr2
+          dd1t1 = rr1/rt
+          dd1t1 = dd1t1*dd1t1
+          dd1t2 = (zrl(l) - zc)/zt
+          dd1t2 = dd1t2*dd1t2 
+          dd1 = dd1t1 + dd1t2  
+          dd2t1 = rr2/rt
+          dd2t1 = dd2t1*dd2t1
+          dd2 = dd2t1 + dd1t2
+          IF ( dd1 .LT. 0.5 ) Then 
+             hx(n,l) = 1.0
+          ELSEIF ( dd2 .LT. 0.5 ) Then 
+             hx(n,l) = 1.0
+          ELSE
+             hx(n,l) = 0.1 
+          END IF
+          IF ( zrl(l) .GT. zc ) Then 
+             IF ( ABS(latc1 - lat_i(n)) .LT. 0.125 ) Then 
+                hx(n,l)= 0.1
+             ENDIF
+          ENDIF
+       END DO
+    END DO
+  END SUBROUTINE slotted_cylinders
     
-      DO l=1,llm 
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL dist_lonlat(lon0,lat0,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            IF ( rr1 .LT. R0 ) then 
-              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
-	    ELSE
-              hx(n,l)=0.0
-	    END IF
-	  END DO
-        END DO
-      END DO
-   
-    END SUBROUTINE cosine_bell_1
-
 !==============================================================================
 
-    SUBROUTINE   cosine_bell_2(hx) 
-    IMPLICIT NONE 	
-    REAL(rstd) :: hx(iim*jjm,llm)
-   
-      DO l=1,llm
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center 
-            rr2 = radius*rr2
-            dd1t1 = rr1/rt
-            dd1t1 = dd1t1*dd1t1
-            dd1t2 = (zrl(l) - zc)/zt
-            dd1t2 = dd1t2*dd1t2 
-            dd1 = dd1t1 + dd1t2 
-            dd1 = Min(1.0,dd1)  
-            dd2t1 = rr2/rt
-            dd2t1 = dd2t1*dd2t1
-            dd2 = dd2t1 + dd1t2 
-            dd2 = Min(1.0,dd2) 
-            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
-          END DO
-        END DO
-      END DO
-
-    END SUBROUTINE cosine_bell_2
-
-!=============================================================================
-
-    SUBROUTINE slotted_cylinders(hx) 
-    IMPLICIT NONE
-    REAL(rstd) :: hx(iim*jjm,llm)
-    
-      DO l=1,llm
-        DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-            n=(j-1)*iim+i
-            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center 
-            rr1 = radius*rr1
-            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center 
-            rr2 = radius*rr2
-            dd1t1 = rr1/rt
-            dd1t1 = dd1t1*dd1t1
-            dd1t2 = (zrl(l) - zc)/zt
-            dd1t2 = dd1t2*dd1t2 
-            dd1 = dd1t1 + dd1t2  
-            dd2t1 = rr2/rt
-            dd2t1 = dd2t1*dd2t1
-            dd2 = dd2t1 + dd1t2
-            
-            IF ( dd1 .LT. 0.5 ) Then 
-              hx(n,l) = 1.0
-            ELSEIF ( dd2 .LT. 0.5 ) Then 
-              hx(n,l) = 1.0
-            ELSE
-	      hx(n,l) = 0.1 
-            END IF  
-    
-            IF ( zrl(l) .GT. zc ) Then 
-              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then 
-                hx(n,l)= 0.1
- 	      ENDIF 
-            ENDIF  
-  
-          ENDDO
-       END DO
-      END DO
-
-    END SUBROUTINE slotted_cylinders 
-
-!==============================================================================
-
     SUBROUTINE hadleyq(hx) 
-      IMPLICIT NONE 
-      REAL(rstd)::hx(iim*jjm,llm) 
+      REAL(rstd)::hx(ngrid,llm) 
       REAL(rstd),PARAMETER:: zz1=2000.,zz2=5000.,zz0=0.5*(zz1+zz2)
+      INTEGER :: n,l
       
       DO l=1,llm
Index: codes/icosagcm/trunk/src/etat0_dcmip2.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 343)
+++ codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 344)
@@ -3,42 +3,25 @@
 ! test cases DCMIP 2012, category 2 : Orographic gravity waves
 
-  USE genmod
-
+  USE icosa
   PRIVATE
 
-  PUBLIC  etat0
+  INTEGER, SAVE :: testcase
+  INTEGER, PARAMETER :: mountain=0, schaer_noshear=1, schaer_shear=2
+
+  PUBLIC getin_etat0, 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, icase
+  SUBROUTINE getin_etat0
     CHARACTER(len=255) :: etat0_type
     etat0_type='jablonowsky06'
     CALL getin("etat0",etat0_type)
-
     SELECT CASE (TRIM(etat0_type))
     CASE('dcmip2_mountain')
-       icase=0
+       testcase = mountain
     CASE('dcmip2_schaer_noshear')
-       icase=1
+       testcase = schaer_noshear
     CASE('dcmip2_schaer_shear')
-       icase=2
+       testcase = schaer_shear
     CASE DEFAULT
        PRINT *, 'This should not happen : etat0_type =', TRIM(etat0_type), ' in etat0_dcmip2.f90/etat0'
@@ -46,44 +29,20 @@
     END SELECT
     PRINT *, 'Orographic gravity-wave test case :', TRIM(etat0_type)
+  END SUBROUTINE getin_etat0
 
-    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis, ps, Temp, ulon, ulat)
+    USE disvert_mod
+    IMPLICIT NONE
+    INTEGER, INTENT(IN)    :: ngrid
+    REAL(rstd),INTENT(IN)  :: lon(ngrid)
+    REAL(rstd),INTENT(IN)  :: lat(ngrid)
+    REAL(rstd),INTENT(OUT) :: ps(ngrid)
+    REAL(rstd),INTENT(OUT) :: phis(ngrid)
+    REAL(rstd),INTENT(OUT) :: Temp(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm)
+    REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm)
+    INTEGER :: l,ij
+    REAL(rstd) :: hyam, hybm
 
-    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_DCMIP2(icase,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_DCMIP2(icase, ps, phis, u, temp,q)
-    USE icosa
-    USE disvert_mod, ONLY : ap,bp
-    USE pression_mod
-    USE theta2theta_rhodz_mod
-    USE wind_mod
-    IMPLICIT NONE
-    INTEGER, INTENT(IN) :: icase
-    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)
-    REAL(rstd) :: ulon(3*iim*jjm,llm), ulat(3*iim*jjm,llm)
-    
-    INTEGER :: i,j,l,ij
-    REAL(rstd) :: hyam, hybm, pp, dummy1, dummy2, dummy3  
-    
     ! Hexagons : ps,phis,temp
     DO l=1,llm
@@ -91,56 +50,29 @@
        hyam = .5*(ap(l)+ap(l+1))/preff
        hybm = .5*(bp(l)+bp(l+1))
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             CALL comp_all(lon_i(ij),lat_i(ij), ps(ij),phis(ij),temp(ij,l), dummy1,dummy2)
-          END DO
+       DO ij=1,ngrid
+          CALL comp_all(hyam,hybm,lon(ij),lat(ij), &
+               ps(ij),phis(ij), Temp(ij,l),ulon(ij,l), ulat(ij,l) )
        END DO
     END DO
-    
-    ! Edges : ulon,ulat
-    DO l=1,llm
-       ! The surface pressure is not set yet so we provide the hybrid coefficients
-       hyam = .5*(ap(l)+ap(l+1))/preff
-       hybm = .5*(bp(l)+bp(l+1))
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             CALL comp_all(lon_e(ij+u_right), lat_e(ij+u_right), &
-                  dummy1,dummy2,dummy3, ulon(ij+u_right,l),ulat(ij+u_right,l))
-             CALL comp_all(lon_e(ij+u_lup), lat_e(ij+u_lup), &
-                  dummy1,dummy2,dummy3, ulon(ij+u_lup,l),ulat(ij+u_lup,l))
-             CALL comp_all(lon_e(ij+u_ldown), lat_e(ij+u_ldown), &
-                  dummy1,dummy2,dummy3, ulon(ij+u_ldown,l),ulat(ij+u_ldown,l))
-          END DO
-       END DO
-    END DO
-    CALL compute_wind_perp_from_lonlat_compound(ulon,ulat,u)
+  END SUBROUTINE compute_etat0
 
-    q=1.
-
-    CONTAINS
-      
-      SUBROUTINE comp_all(lon,lat, psj,phisj,tempj, ulonj,ulatj)
-        USE dcmip_initial_conditions_test_1_2_3
-        REAL(rstd), INTENT(IN) :: lon, lat
-        REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj
-        REAL :: dummy
-        dummy=0.
-        SELECT CASE (icase)
-        CASE(0)
-           CALL test2_steady_state_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm, &
-                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
-        CASE(1)
-           CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,0,&
-                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
-        CASE(2)
-           CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,1, &
-                ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
-        END SELECT
-      END SUBROUTINE comp_all
-
-  END SUBROUTINE compute_etat0_DCMIP2
-
+  SUBROUTINE comp_all(hyam,hybm,lon,lat, psj,phisj,tempj, ulonj,ulatj)
+    USE dcmip_initial_conditions_test_1_2_3
+    REAL(rstd), INTENT(IN) :: hyam, hybm, lon, lat
+    REAL(rstd), INTENT(OUT) :: psj,phisj,tempj,ulonj,ulatj
+    REAL :: dummy
+    dummy=0.
+    SELECT CASE (testcase)
+    CASE(mountain)
+       CALL test2_steady_state_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm, &
+            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
+    CASE(schaer_noshear)
+       CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,0,&
+            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
+    CASE(schaer_shear)
+       CALL test2_schaer_mountain(lon,lat,dummy,dummy,0,.TRUE.,hyam,hybm,1, &
+            ulonj,ulatj,dummy,tempj,phisj,psj,dummy,dummy)
+    END SELECT
+  END SUBROUTINE comp_all
 
 END MODULE etat0_dcmip2_mod
