Index: codes/icosagcm/devel/src/initial/etat0_collocated.f90
===================================================================
--- codes/icosagcm/devel/src/initial/etat0_collocated.f90	(revision 848)
+++ codes/icosagcm/devel/src/initial/etat0_collocated.f90	(revision 856)
@@ -1,4 +1,5 @@
 MODULE etat0_collocated_mod
   USE icosa
+  USE grid_param
   USE omp_para, ONLY : is_omp_level_master
   USE caldyn_vars_mod, ONLY : hydrostatic
@@ -6,7 +7,6 @@
   PRIVATE
 
-  LOGICAL :: autoinit_mass, autoinit_NH
   CHARACTER(len=255),SAVE :: etat0_type
-!$OMP THREADPRIVATE(autoinit_mass, autoinit_NH, etat0_type)
+!$OMP THREADPRIVATE(etat0_type)
 
     PUBLIC :: etat0_type, etat0_collocated
@@ -36,5 +36,5 @@
     TYPE(t_field),POINTER :: f_q(:)
   
-    TYPE(t_field),POINTER,SAVE :: f_temp(:)
+    TYPE(t_field),POINTER,SAVE :: f_temp(:) ! SAVE => all threads see the same f_temp
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: mass(:,:)
@@ -64,12 +64,13 @@
       q=f_q(ind)
 
-      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
-         CALL compute_etat0_collocated_hex(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q)
-      ELSE
-         CALL compute_etat0_collocated_hex(ps,mass, phis, temp, u, geopot, W, q)
-      ENDIF
-
-      IF( TRIM(etat0_type)/='williamson91.6' ) CALL compute_temperature2entropy(ps,temp,q,theta_rhodz, 1)
-    
+      SELECT CASE(grid_type)
+      CASE(grid_ico)
+         CALL compute_etat0_collocated_hex(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q)
+      CASE(grid_unst)
+         CALL compute_etat0_collocated_unst(phis, ps, mass, theta_rhodz(:,:,1), u, geopot, W, q)
+      CASE DEFAULT
+         STOP 'Unexpected value of grid_type encountered in etat0_collocated.'
+      END SELECT
+
     ENDDO
     
@@ -78,62 +79,9 @@
   END SUBROUTINE etat0_collocated
 
-  SUBROUTINE compute_temperature2entropy(ps,temp,q,theta_rhodz,offset)
-    USE icosa
-    USE pression_mod
-    USE exner_mod
-    USE omp_para
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: q(iim*jjm,llm,nqtot)
-    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-    INTEGER,INTENT(IN) :: offset
-
-    REAL(rstd) :: p(iim*jjm,llm+1)
-    REAL(rstd) :: cppd,Rd, mass, p_ij, q_ij,r_ij, chi,nu, entropy, theta
-    INTEGER :: i,j,ij,l
-
-    cppd=cpp
-    Rd=kappa*cppd
-
-    CALL compute_pression(ps,p,offset)
-    ! flush p
-    DO    l    = ll_begin, ll_end
-       DO j=jj_begin-offset,jj_end+offset
-          DO i=ii_begin-offset,ii_end+offset
-             ij=(j-1)*iim+i
-             mass = (p(ij,l)-p(ij,l+1))/g ! dry+moist mass
-             p_ij = .5*(p(ij,l)+p(ij,l+1))  ! pressure at full level
-             SELECT CASE(caldyn_thermo)
-             CASE(thermo_theta)
-                theta = temp(ij,l)*(p_ij/preff)**(-kappa) 
-                theta_rhodz(ij,l) = mass * theta
-             CASE(thermo_entropy)
-                nu = log(p_ij/preff)
-                chi = log(temp(ij,l)/Treff)
-                entropy = cppd*chi-Rd*nu
-                theta_rhodz(ij,l) = mass * entropy
-!             CASE(thermo_moist)
-!                q_ij=q(ij,l,1)
-!                r_ij=1.-q_ij
-!                mass=mass*(1-q_ij) ! dry mass
-!                nu = log(p_ij/preff)
-!                chi = log(temp(ij,l)/Treff)
-!                entropy = r_ij*(cppd*chi-Rd*nu) + q_ij*(cppv*chi-Rv*nu)
-!                theta_rhodz(ij,l) = mass * entropy                
-                CASE DEFAULT
-                   STOP
-             END SELECT
-          ENDDO
-       ENDDO
-    ENDDO
-  END SUBROUTINE compute_temperature2entropy
-
-  SUBROUTINE compute_etat0_collocated_hex(ps,mass,phis,temp_i,u, geopot,W, q)
-    USE wind_mod
-    USE disvert_mod
+  SUBROUTINE compute_etat0_collocated_hex(phis,ps,mass,theta_rhodz,u, geopot,W, q)
     REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: W(iim*jjm,llm+1)
@@ -141,49 +89,47 @@
     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
 
-    REAL(rstd) :: ulon_i(iim*jjm,llm)
-    REAL(rstd) :: ulat_i(iim*jjm,llm)
-
     REAL(rstd) :: ps_e(3*iim*jjm)
+    REAL(rstd) :: phis_e(3*iim*jjm)
     REAL(rstd) :: mass_e(3*iim*jjm,llm)
-    REAL(rstd) :: phis_e(3*iim*jjm)
-    REAL(rstd) :: temp_e(3*iim*jjm,llm)
     REAL(rstd) :: geopot_e(3*iim*jjm,llm+1)
-    REAL(rstd) :: ulon_e(3*iim*jjm,llm)
-    REAL(rstd) :: ulat_e(3*iim*jjm,llm)
     REAL(rstd) :: q_e(3*iim*jjm,llm,nqtot)
 
-    INTEGER :: l,i,j,ij
-    REAL :: p_ik, v_ik, mass_ik
-
-    ! For NH geopotential and vertical momentum must be initialized.
-    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
-    ! to hydrostatic geopotential and zero
-    autoinit_mass = .TRUE.
-    autoinit_NH = .NOT. hydrostatic
     w(:,:) = 0
-
-    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
-    IF(autoinit_NH) THEN
-       geopot(:,1) = phis(:) ! surface geopotential
-       DO l = 1, llm
-          DO ij=1,iim*jjm
-             ! hybrid pressure coordinate
-             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
-             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
-             ! v=R.T/p, R=kappa*cpp
-             v_ik = kappa*cpp*temp_i(ij,l)/p_ik
-             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
-          END DO
-       END DO
-    END IF
-
-    CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
+    CALL compute_etat0_collocated(iim*jjm  , lon_i, lat_i, phis,   ps,   mass,   theta_rhodz, geopot,   q)
+    CALL compute_etat0_collocated(3*iim*jjm, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u)
 
   END SUBROUTINE compute_etat0_collocated_hex
 
-  SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, temp, ulon, ulat, geopot, q)
+  SUBROUTINE compute_etat0_collocated_unst(phis,ps,mass,theta_rhodz,u, geopot,W, q)
+    REAL(rstd),INTENT(INOUT) :: ps(primal_num)
+    REAL(rstd),INTENT(INOUT) :: mass(llm,primal_num)
+    REAL(rstd),INTENT(OUT) :: theta_rhodz(llm,primal_num)
+    REAL(rstd),INTENT(OUT) :: phis(primal_num)
+    REAL(rstd),INTENT(OUT) :: u(llm, edge_num)
+    REAL(rstd),INTENT(OUT) :: W(llm+1, primal_num)
+    REAL(rstd),INTENT(OUT) :: geopot(llm+1, primal_num)
+    REAL(rstd),INTENT(OUT) :: q(llm, primal_num, nqtot)
+    ! The 2D/3D arrays below have the shape expected by compute_etat0_collocated
+    REAL(rstd) :: ps_e(edge_num)
+    REAL(rstd) :: phis_e(edge_num)
+    REAL(rstd) :: u_e(edge_num, llm)
+    REAL(rstd) :: mass_i(primal_num, llm), theta_rhodz_i(primal_num, llm), mass_e(edge_num, llm)
+    REAL(rstd) :: geopot_i(edge_num, llm+1), geopot_e(edge_num, llm+1)
+    REAL(rstd) :: q_i(primal_num, llm, nqtot), q_e(edge_num, llm, nqtot)
+    INTEGER :: iq
+
+    w(:,:) = 0
+    CALL compute_etat0_collocated(primal_num  , lon_i, lat_i, phis,   ps,   mass_i,   theta_rhodz_i, geopot_i,   q_i)
+    CALL compute_etat0_collocated(edge_num, lon_e, lat_e, phis_e, ps_e, mass_e, mass_e, geopot_e, q_e, ep_e, u_e)
+    mass = TRANSPOSE(mass_i)
+    theta_rhodz = TRANSPOSE(theta_rhodz_i)
+    geopot = TRANSPOSE(geopot_i)
+    u = TRANSPOSE(u_e)
+    DO iq=1,nqtot
+       q(:,:,iq) = TRANSPOSE(q_i(:,:,iq))
+    END DO
+  END SUBROUTINE compute_etat0_collocated_unst
+  
+  SUBROUTINE compute_etat0_collocated(ngrid, lon, lat, phis, ps, mass, theta_rhodz, geopot, q, ep, uperp)
     USE etat0_isothermal_mod, ONLY : compute_isothermal => compute_etat0
     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
@@ -199,14 +145,29 @@
     USE etat0_dcmip2016_cyclone_mod, ONLY : compute_dcmip2016_cyclone => compute_etat0
     USE etat0_dcmip2016_supercell_mod, ONLY : compute_dcmip2016_supercell => compute_etat0
+    USE disvert_mod, ONLY : ptop, ap, bp, mass_ak, mass_bk, mass_dak, mass_dbk
     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)   :: theta_rhodz(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)
+    REAL(rstd),INTENT(OUT), OPTIONAL   :: ep(ngrid,3), uperp(ngrid,llm)
+
+    REAL(rstd)   :: ulon(ngrid,llm)
+    REAL(rstd)   :: ulat(ngrid,llm)
+    REAL(rstd)   :: temp(ngrid,llm)
+
+    LOGICAL :: autoinit_mass, autoinit_NH
+    INTEGER :: ij,l
+    REAL(rstd) :: clat, slat, clon, slon, u3d(3), p_ik, mass_ik, v_ik
+    REAL(rstd) :: q_ik, r_ik, chi, entropy, theta, log_p_preff
+
+    ! For NH, geopotential and vertical momentum must be initialized.
+    ! Unless autoinit_NH is set to .FALSE. , they will be initialized
+    ! to hydrostatic geopotential and zero
+    autoinit_mass = .TRUE.
+    autoinit_NH = .NOT. hydrostatic
 
     SELECT CASE (TRIM(etat0_type))
@@ -232,5 +193,5 @@
 !       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))
+       CALL compute_w91_6(ngrid, lon, lat, phis, mass(:,1), theta_rhodz(:,1), ulon(:,1), ulat(:,1))
        autoinit_mass = .FALSE. ! do not overwrite mass
     CASE('dcmip2016_baroclinic_wave')
@@ -242,4 +203,63 @@
     END SELECT
 
+    IF(PRESENT(uperp)) THEN
+       DO l = 1, llm
+          DO ij=1,ngrid
+             clat = COS(lat(ij))
+             slat = SIN(lat(ij))
+             clon = COS(lon(ij))
+             slon = SIN(lon(ij))
+             u3d(1) = -clon*slat*ulat(ij,l) - slon*ulon(ij,l)
+             u3d(2) = -slon*slat*ulat(ij,l) + clon*ulon(ij,l)
+             u3d(3) =       clat*ulat(ij,l)    
+             uperp(ij,l) = SUM(u3d*ep(ij,:))
+          END DO
+       END DO
+    END IF
+
+    IF(autoinit_mass) THEN
+       DO l = 1, llm
+          DO ij=1,ngrid
+             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
+             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
+             mass(ij,l) = mass_ik
+             SELECT CASE(caldyn_thermo)
+             CASE(thermo_theta)
+                theta = temp(ij,l)*(p_ik/preff)**(-kappa)
+                theta_rhodz(ij,l) = mass_ik * theta
+             CASE(thermo_entropy)
+                log_p_preff = log(p_ik/preff)
+                chi = log(temp(ij,l)/Treff)
+                entropy = cpp*chi-Rd*log_p_preff
+                theta_rhodz(ij,l) = mass_ik * entropy
+             CASE(thermo_moist)
+                q_ik=q(ij,l,1)
+                r_ik=1.-q_ik
+                mass_ik = mass_ik*(1.-q_ik) ! dry mass
+                log_p_preff = log(p_ik/preff)
+                chi = log(temp(ij,l)/Treff)
+                entropy = r_ik*(cpp*chi-Rd*nu) + q_ik*(cppv*chi-Rv*log_p_preff)
+                theta_rhodz(ij,l) = mass_ik * entropy
+             CASE DEFAULT
+                STOP
+             END SELECT
+          END DO
+       END DO
+    END IF
+
+    IF(autoinit_NH) THEN
+       geopot(:,1) = phis(:) ! surface geopotential
+       DO l = 1, llm
+          DO ij=1,ngrid
+             ! hybrid pressure coordinate
+             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij)
+             mass_ik = (mass_dak(l) + mass_dbk(l)*ps(ij))/g
+             ! v=R.T/p, R=kappa*cpp
+             v_ik = kappa*cpp*temp(ij,l)/p_ik
+             geopot(ij,l+1) = geopot(ij,l) + mass_ik*v_ik*g
+          END DO
+       END DO
+    END IF
+
   END SUBROUTINE compute_etat0_collocated
 
