Index: codes/icosagcm/devel/src/initial/etat0_collocated.f90
===================================================================
--- codes/icosagcm/devel/src/initial/etat0_collocated.f90	(revision 847)
+++ codes/icosagcm/devel/src/initial/etat0_collocated.f90	(revision 848)
@@ -6,6 +6,7 @@
   PRIVATE
 
-    CHARACTER(len=255),SAVE :: etat0_type
-!$OMP THREADPRIVATE(etat0_type)
+  LOGICAL :: autoinit_mass, autoinit_NH
+  CHARACTER(len=255),SAVE :: etat0_type
+!$OMP THREADPRIVATE(autoinit_mass, autoinit_NH, etat0_type)
 
     PUBLIC :: etat0_type, etat0_collocated
@@ -64,7 +65,7 @@
 
       IF( TRIM(etat0_type)=='williamson91.6' ) THEN
-         CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q)
+         CALL compute_etat0_collocated_hex(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q)
       ELSE
-         CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
+         CALL compute_etat0_collocated_hex(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
 
@@ -128,20 +129,7 @@
   END SUBROUTINE compute_temperature2entropy
 
-  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
+  SUBROUTINE compute_etat0_collocated_hex(ps,mass,phis,temp_i,u, geopot,W, q)
     USE wind_mod
     USE disvert_mod
-    USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0
-    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_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
-    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
-    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
-    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0  
-    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
-    USE etat0_temperature_mod, ONLY: compute_temperature => compute_etat0
-    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
-    USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
-    USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
@@ -167,5 +155,4 @@
     INTEGER :: l,i,j,ij
     REAL :: p_ik, v_ik, mass_ik
-    LOGICAL :: autoinit_mass, autoinit_NH
 
     ! For NH geopotential and vertical momentum must be initialized.
@@ -176,48 +163,6 @@
     w(:,:) = 0
 
-    SELECT CASE (TRIM(etat0_type))
-    CASE ('isothermal')
-       CALL compute_isothermal(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_isothermal(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    CASE ('temperature_profile')
-       CALL compute_temperature(iim*jjm, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_temperature(3*iim*jjm, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    CASE('jablonowsky06')
-       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('dcmip3')
-       CALL compute_dcmip3(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
-       CALL compute_dcmip3(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
-       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
-    CASE('dcmip4')
-       CALL compute_dcmip4(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_dcmip4(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)
-       CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    CASE('bubble')
-       CALL compute_bubble(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, geopot, q)
-       CALL compute_bubble(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
-!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
-    CASE('williamson91.6')
-       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
-       CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1))
-       autoinit_mass = .FALSE. ! do not overwrite mass
-    CASE('dcmip2016_baroclinic_wave')
-       CALL compute_dcmip2016_baroclinic_wave(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_dcmip2016_baroclinic_wave(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    CASE('dcmip2016_cyclone')
-       CALL compute_dcmip2016_cyclone(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_dcmip2016_cyclone(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    CASE('dcmip2016_supercell')
-       CALL compute_dcmip2016_supercell(iim*jjm,lon_i,lat_i, phis, ps, temp_i, ulon_i, ulat_i, q)
-       CALL compute_dcmip2016_supercell(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
-    END SELECT
+    CALL compute_etat0_collocated(iim*jjm  , lon_i, lat_i, phis,   ps,   mass,   temp_i, ulon_i, ulat_i, geopot,   q)
+    CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, temp_e, ulon_e, ulat_e, geopot_e, q_e)
 
     IF(autoinit_mass) CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
@@ -238,4 +183,63 @@
     CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
 
+  END SUBROUTINE compute_etat0_collocated_hex
+
+  SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, temp, ulon, ulat, geopot, q)
+    USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0
+    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_dcmip3_mod, ONLY : compute_dcmip3 => compute_etat0
+    USE etat0_dcmip4_mod, ONLY : compute_dcmip4 => compute_etat0
+    USE etat0_dcmip5_mod, ONLY : compute_dcmip5 => compute_etat0
+    USE etat0_bubble_mod, ONLY : compute_bubble => compute_etat0  
+    USE etat0_williamson_mod, ONLY : compute_w91_6 => compute_etat0
+    USE etat0_temperature_mod, ONLY: compute_temperature => compute_etat0
+    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : compute_dcmip2016_baroclinic_wave => compute_etat0
+    USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
+    USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
+    INTEGER :: ngrid
+    REAL(rstd),INTENT(IN)    :: lon(ngrid), lat(ngrid)
+    REAL(rstd),INTENT(INOUT) :: ps(ngrid)
+    REAL(rstd),INTENT(INOUT) :: mass(ngrid,llm)
+    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)
+    REAL(rstd),INTENT(OUT)   :: geopot(ngrid,llm+1)
+    REAL(rstd),INTENT(OUT)   :: q(ngrid,llm,nqtot)
+
+    SELECT CASE (TRIM(etat0_type))
+    CASE ('isothermal')
+       CALL compute_isothermal(ngrid, phis, ps, temp, ulon, ulat, q)
+    CASE ('temperature_profile')
+       CALL compute_temperature(ngrid, phis, ps, temp, ulon, ulat, q)
+    CASE('jablonowsky06')
+       CALL compute_jablonowsky06(ngrid, lon, lat, phis, ps, temp, ulon, ulat)
+    CASE('dcmip1')
+       CALL compute_dcmip1(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
+       CALL compute_dcmip2(ngrid, lon, lat, phis, ps, temp, ulon, ulat)
+    CASE('dcmip3')
+       CALL compute_dcmip3(ngrid, lon, lat, phis, ps, temp, ulon, ulat, geopot, q)
+       autoinit_NH = .FALSE. ! compute_dcmip3 initializes geopot
+    CASE('dcmip4')
+       CALL compute_dcmip4(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    CASE('dcmip5')
+       CALL compute_dcmip5(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    CASE('bubble')
+       CALL compute_bubble(ngrid, lon, lat, phis, ps, temp, ulon, ulat, geopot, q)
+!       autoinit_NH = .FALSE. ! compute_bubble initializes geopot
+    CASE('williamson91.6')
+       CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), temp(:,1), ulon(:,1), ulat(:,1))
+       autoinit_mass = .FALSE. ! do not overwrite mass
+    CASE('dcmip2016_baroclinic_wave')
+       CALL compute_dcmip2016_baroclinic_wave(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    CASE('dcmip2016_cyclone')
+       CALL compute_dcmip2016_cyclone(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    CASE('dcmip2016_supercell')
+       CALL compute_dcmip2016_supercell(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
+    END SELECT
+
   END SUBROUTINE compute_etat0_collocated
 
