Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 365)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 366)
@@ -12,5 +12,6 @@
   
   SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
-       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du_slow, f_du_fast)
+       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, &
+       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast) 
     USE icosa
     USE observable_mod
@@ -34,4 +35,5 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_W(:)
     TYPE(t_field),POINTER :: f_geopot(:)
     TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
@@ -41,9 +43,13 @@
     TYPE(t_field) :: f_du_slow(:)
     TYPE(t_field) :: f_du_fast(:)
+    TYPE(t_field) :: f_dW_slow(:)
+    TYPE(t_field) :: f_dW_fast(:)
+    TYPE(t_field) :: f_dPhi_slow(:)
+    TYPE(t_field) :: f_dPhi_fast(:)
     
     REAL(rstd),POINTER :: ps(:), dps(:)
     REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
-    REAL(rstd),POINTER :: du(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: u(:,:), qu(:,:), qv(:,:)
+    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
+    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:)
 
 ! temporary shared variable
@@ -68,4 +74,8 @@
       CALL init_message(f_u,req_e1_vect,req_u)
       CALL init_message(f_qu,req_e1_scal,req_qu)
+      IF(.NOT.hydrostatic) THEN
+         CALL init_message(f_geopot,req_i1,req_geopot)
+         CALL init_message(f_w,req_i1,req_w)
+      END IF
     ENDIF
     
@@ -81,4 +91,11 @@
     CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
     CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
+
+    IF(.NOT.hydrostatic) THEN
+       CALL send_message(f_geopot,req_geopot) ! COM03
+       CALL wait_message(req_geopot) ! COM03
+       CALL send_message(f_w,req_w) ! COM04
+       CALL wait_message(req_w) ! COM04
+    END IF
     
     DO ind=1,ndomain
@@ -87,14 +104,21 @@
        CALL swap_geometry(ind)
        ps=f_ps(ind)
-       u=f_u(ind)
-       du=f_du_fast(ind)
        theta_rhodz=f_theta_rhodz(ind)
        mass=f_mass(ind)
        theta = f_theta(ind)
+       CALL compute_theta(ps,theta_rhodz, mass,theta)
        pk = f_pk(ind)
        geopot = f_geopot(ind)
-       CALL compute_theta(ps,theta_rhodz, mass,theta)
-       CALL compute_geopot(ps,mass,theta, pk,geopot)
-       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot, du) ! computes du_fast and updates u
+       IF(hydrostatic) THEN
+          CALL compute_geopot(ps,mass,theta, pk,geopot)
+       ELSE
+          W = f_W(ind)
+          dW = f_dW_fast(ind)
+          dPhi = f_dPhi_fast(ind)
+          CALL compute_caldyn_solver(tau,mass,theta,pk,geopot,W,dPhi,dW) ! computes d(Phi,W)_fast and updates Phi,W
+       END IF
+       u=f_u(ind)
+       du=f_du_fast(ind)
+       CALL compute_caldyn_fast(tau,u,mass,theta,pk,geopot,du) ! computes du_fast and updates du
     ENDDO
     
Index: codes/icosagcm/trunk/src/caldyn_kernels_base.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 365)
+++ codes/icosagcm/trunk/src/caldyn_kernels_base.f90	(revision 366)
@@ -18,5 +18,5 @@
   !$OMP THREADPRIVATE(caldyn_conserv) 
 
-  TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
+  TYPE(t_message),PUBLIC :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w
 
   PUBLIC :: compute_geopot, compute_caldyn_vert
@@ -51,5 +51,4 @@
 
     IF(caldyn_eta==eta_mass .AND. .NOT. DEC) THEN
-
 !!! Compute exner function and geopotential
        DO l = 1,llm
@@ -64,8 +63,8 @@
        ENDDO
     ELSE 
-       ! We are using a Lagrangian vertical coordinate
-       ! Pressure must be computed first top-down (temporarily stored in pk)
+       ! We are using DEC or a Lagrangian vertical coordinate
+       ! Pressure is computed first top-down (temporarily stored in pk)
        ! Then Exner pressure and geopotential are computed bottom-up
-       ! Notice that the computation below should work also when caldyn_eta=eta_mass          
+       ! Works also when caldyn_eta=eta_mass          
 
        IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
Index: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 365)
+++ codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90	(revision 366)
@@ -6,5 +6,6 @@
   PRIVATE
 
-  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_fast, compute_caldyn_slow
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_slow, &
+       compute_caldyn_solver, compute_caldyn_fast
 
 CONTAINS
@@ -120,5 +121,92 @@
   END SUBROUTINE compute_pvort_only
 
-  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot, du)
+  SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk, geopot,W, dPhi,dW)
+    USE icosa
+    USE disvert_mod
+    USE exner_mod
+    USE trace
+    USE omp_para
+    USE disvert_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
+    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
+
+    INTEGER :: ij,l
+    INTEGER    :: ij_omp_begin_ext, ij_omp_end_ext
+    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij, m_il
+
+    CALL trace_start("compute_caldyn_solver")
+
+    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
+    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
+    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
+
+    IF(tau>0) THEN
+       PRINT *, 'No implicit solver for NH yet'
+       STOP
+    END IF
+
+    ! Compute pressure, stored temporarily in pk
+    ! kappa = R/Cp
+    ! 1-kappa = Cv/Cp
+    ! Cp/Cv = 1/(1-kappa)
+    gamma = 1./(1.-kappa)
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
+          X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
+          ! kappa.theta.rho = p/exner 
+          ! => X = (p/p0)/(exner/Cp)
+          !      = (p/p0)^(1-kappa)
+          pk(ij,l) = preff*(X_ij**gamma)
+       ENDDO
+    ENDDO
+
+    ! Compute tendencies for geopotential and vertical momentum
+    DO l=ll_beginp1,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          m_il = .5*(rhodz(ij,l)+rhodz(ij,l-1))
+          dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il
+          dPhi(ij,l) = g*g*W(ij,l)/m_il
+       ENDDO
+       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
+       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
+    ENDDO
+    ! Lower BC (FIXME : no orography yet !)
+    DO ij=ij_begin,ij_end         
+       dPhi(ij,1)=0
+       W(ij,1)=0
+       dW(ij,1)=0
+       dPhi(ij,llm+1)=0 ! rigid lid
+       W(ij,llm+1)=0
+       dW(ij,llm+1)=0
+    ENDDO
+    ! Upper BC p=ptop
+!    DO ij=ij_omp_begin_ext,ij_omp_end_ext
+!       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
+!       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
+!    ENDDO
+    
+    ! Compute Exner function (needed by compute_caldyn_fast)
+    DO l=ll_begin,ll_end
+       !DIR$ SIMD
+       DO ij=ij_omp_begin_ext,ij_omp_end_ext
+          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
+       ENDDO
+    ENDDO
+    
+    CALL trace_end("compute_caldyn_solver")
+    
+  END SUBROUTINE compute_caldyn_solver
+  
+  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
     USE icosa
     USE disvert_mod
@@ -127,11 +215,11 @@
     USE omp_para
     IMPLICIT NONE
-    REAL(rstd), INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
-    REAL(rstd),INTENT(INOUT)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
-    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
+    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
+    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT)   :: du(iim*3*jjm,llm)
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
 
@@ -140,7 +228,8 @@
 
     CALL trace_start("compute_caldyn_fast")
-    
-    ! Compute bernouilli term
+
+    ! Compute Bernoulli term
     IF(boussinesq) THEN
+
        DO l=ll_begin,ll_end
           !DIR$ SIMD
Index: codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90
===================================================================
--- codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 365)
+++ codes/icosagcm/trunk/src/dcmip_initial_conditions_test_1_2_3_v5.f90	(revision 366)
@@ -1203,14 +1203,14 @@
 !-----------------------------------------------------------------------
 
-	real(rstd), intent(in)  :: lon, &		! Longitude (radians)
-				lat, &		! Latitude (radians)
-				z		! Height (m)
-
-	real(rstd), intent(inout) :: p		! Pressure  (Pa)
-				
+	real(rstd), intent(in)  :: lon, &	! Longitude (radians)
+                                   lat		! Latitude (radians)
+
+	real(rstd), intent(inout) :: p, &	! Pressure  (Pa)
+				     z		! Height (m)
+
 
 	integer,  intent(in) :: zcoords 	! 0 or 1 see below
 
-	real(rstd), intent(out) :: u, & 		! Zonal wind (m s^-1)
+	real(rstd), intent(out) :: u, & 	! Zonal wind (m s^-1)
 				v, &		! Meridional wind (m s^-1)
 				w, &		! Vertical Velocity (m s^-1)
@@ -1230,5 +1230,6 @@
 			    	Om      = 0.d0,			&	! Rotation Rate of Earth
                             	as      = a/X,			&	! New Radius of small Earth     
-			    	u0      = 20.d0,		&	! Reference Velocity 
+!			    	u0      = 20.d0,		&	! Reference Velocity 
+			    	u0      = 0.d0,		        &	! FIXME : no zonal wind for NH tests
                             	Teq     = 300.d0,		&	! Temperature at Equator    
 			    	Peq     = 100000.d0,		&	! Reference PS at Equator
@@ -1297,4 +1298,5 @@
 
 		height = (-g/N2)*log( (TS/bigG)*( (p/ps)**(Rd/cp) - 1.d0  ) + 1.d0 )
+                z = height ! modified : initialize z when pressure-based
 
 	endif
Index: codes/icosagcm/trunk/src/earth_const.f90
===================================================================
--- codes/icosagcm/trunk/src/earth_const.f90	(revision 365)
+++ codes/icosagcm/trunk/src/earth_const.f90	(revision 366)
@@ -16,5 +16,5 @@
   REAL(rstd),SAVE :: mu                 ! molar mass of the atmosphere
 
-  LOGICAL, SAVE :: boussinesq
+  LOGICAL, SAVE :: boussinesq, hydrostatic
 
 CONTAINS
@@ -34,8 +34,14 @@
     CALL getin("scale_height",scale_height)
     
-    mu=kappa/cpp
     boussinesq=.FALSE.
     CALL getin("boussinesq",boussinesq)  
+    hydrostatic=.TRUE.
+    CALL getin("hydrostatic",hydrostatic)  
+    IF(boussinesq .AND. .NOT. hydrostatic) THEN
+       PRINT *, 'boussinesq=.TRUE. and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported'
+       STOP
+    END IF
     
+    mu=kappa/cpp
     radius=radius/scale_factor
     omega=omega*scale_factor
@@ -45,3 +51,2 @@
   
 END MODULE earth_const
-  
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 365)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 366)
@@ -13,5 +13,5 @@
 CONTAINS
   
-  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q)
     USE mpipara, ONLY : is_mpi_root
     USE disvert_mod
@@ -35,8 +35,10 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_w(:)
     TYPE(t_field),POINTER :: f_q(:)
     
     REAL(rstd),POINTER :: ps(:), mass(:,:)
-    LOGICAL :: init_mass, collocated
+    LOGICAL :: autoinit_mass, autoinit_geopot, collocated
     INTEGER :: ind,i,j,ij,l
 
@@ -45,7 +47,7 @@
     ! the initial distribution of mass is taken to be the same
     ! as what the mass coordinate would dictate
-    ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE.
+    ! however if etat0_XXX defines mass then the flag autoinit_mass must be set to .FALSE.
     ! otherwise mass will be overwritten 
-    init_mass = (caldyn_eta == eta_lag)
+    autoinit_mass = (caldyn_eta == eta_lag)
 
     etat0_type='jablonowsky06'
@@ -70,8 +72,9 @@
         CALL getin_etat0_dcmip5
     CASE ('williamson91.6')
-       init_mass=.FALSE.
+       autoinit_mass=.FALSE.
        CALL getin_etat0_williamson
     CASE DEFAULT
        collocated=.FALSE.
+       autoinit_mass = .FALSE.
     END SELECT
 
@@ -90,5 +93,5 @@
    CASE DEFAULT
       IF(collocated) THEN
-         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
       ELSE
          PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
@@ -98,6 +101,6 @@
     END SELECT
 
-    IF(init_mass) THEN ! initialize mass distribution using ps
 !       !$OMP BARRIER
+    IF(autoinit_mass) THEN
        DO ind=1,ndomain
           IF (.NOT. assigned_domain(ind)) CYCLE
@@ -105,11 +108,11 @@
           CALL swap_geometry(ind)
           mass=f_mass(ind); ps=f_ps(ind)
-          CALL compute_rhodz(.TRUE., ps, mass)
+          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
        END DO
     END IF
-
+ 
   END SUBROUTINE etat0
 
-  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
     USE theta2theta_rhodz_mod
     IMPLICIT NONE
@@ -119,4 +122,6 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_geopot(:)
+    TYPE(t_field),POINTER :: f_W(:)
     TYPE(t_field),POINTER :: f_q(:)
   
@@ -128,4 +133,6 @@
     REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: geopot(:,:)
+    REAL(rstd),POINTER :: W(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
     INTEGER :: ind
@@ -143,10 +150,12 @@
       temp=f_temp(ind)
       u=f_u(ind)
+      geopot=f_geopot(ind)
+      w=f_w(ind)
       q=f_q(ind)
 
-      IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
-       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
+      IF( TRIM(etat0_type)=='williamson91.6' ) THEN
+       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, geopot, W, q)
       ELSE
-       CALL compute_etat0_collocated(ps,mass, phis, temp, u, q)
+       CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
     ENDDO
@@ -158,6 +167,7 @@
   END SUBROUTINE etat0_collocated
 
-  SUBROUTINE compute_etat0_collocated(ps,mass, phis, temp_i, u, q)
+  SUBROUTINE compute_etat0_collocated(ps,mass,phis,temp_i,u, geopot,W, q)
     USE wind_mod
+    USE disvert_mod
     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
     USE etat0_dcmip1_mod, ONLY : compute_dcmip1 => compute_etat0
@@ -174,4 +184,6 @@
     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)
+    REAL(rstd),INTENT(OUT) :: geopot(iim*jjm,llm+1)
     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
 
@@ -183,4 +195,5 @@
     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)
@@ -188,4 +201,13 @@
 
     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.
+    ! Unless autoinit_mass is set to .FALSE. , they will be initialized
+    ! to hydrostatic geopotential and zero
+    autoinit_mass = .TRUE.
+    autoinit_NH = .NOT. hydrostatic
+    w(:,:) = 0
 
     !$OMP BARRIER
@@ -208,6 +230,7 @@
        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)
+       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)
@@ -219,5 +242,21 @@
        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
     END SELECT
+
+    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
 
     !$OMP BARRIER
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 365)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 366)
@@ -7,5 +7,5 @@
 CONTAINS
   
-  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
+  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,geopot,q)
     USE genmod
     USE dcmip_initial_conditions_test_1_2_3
@@ -20,7 +20,8 @@
     REAL(rstd), INTENT(OUT) :: ulat(ngrid,llm)
     REAL(rstd), INTENT(OUT) :: temp(ngrid,llm)
+    REAL(rstd), INTENT(OUT) :: geopot(ngrid,llm+1)
     REAL(rstd), INTENT(OUT) :: q(ngrid,llm,nqtot)
     REAL(rstd),PARAMETER :: Peq=1e5        ! Reference surface pressure at the equator (hPa)
-    REAL(rstd) :: dummy, pp
+    REAL(rstd) :: dummy, pp, zz
     INTEGER :: l,ij
     pp=peq
@@ -28,4 +29,13 @@
        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=ll_begin,ll_endp1
+       DO ij=1,ngrid
+          pp = ap(l) + bp(l)*ps(ij) ! half-layer pressure
+          CALL test3_gravity_wave(lon(ij),lat(ij),pp,zz,0, &
+               dummy,dummy,dummy,dummy,dummy,dummy,dummy,dummy)
+          geopot(ij,l) = g*zz ! initialize geopotential for NH
+       END DO
+       q(:,l,:)=0.
     END DO
     DO l=ll_begin,ll_end
Index: codes/icosagcm/trunk/src/euler_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/euler_scheme.f90	(revision 365)
+++ codes/icosagcm/trunk/src/euler_scheme.f90	(revision 366)
@@ -5,12 +5,13 @@
 
   TYPE(t_field),POINTER,SAVE,PUBLIC :: f_geopot(:), f_q(:), &
-       f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:), &
-       f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:), &
-       f_u(:),f_um1(:),f_um2(:), f_du(:), & ! previous state
-       f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:), & ! trends
-       f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:), & ! accumulated mass fluxes
-                                ! slow/fast trend arrays for HEVI time scheme
-       f_dmass_slow(:,:), f_dps_slow(:,:), f_dtheta_rhodz_slow(:,:), &
-       f_du_slow(:,:), f_du_fast(:,:)
+       f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:), & ! current and previous time steps + tendency of mass,
+       f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:), &           ! column-integrated mass
+       f_u(:),f_um1(:),f_um2(:), f_du(:), f_W(:), &                   ! 'horizontal' and vertical (NH) momentum
+       f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:), & ! mass-weighted theta/entropy
+       f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:), &            ! accumulated mass fluxes
+       f_du_slow(:,:), f_du_fast(:,:), &                              ! slow/fast trend arrays 
+       f_dmass_slow(:,:), f_dps_slow(:,:), f_dtheta_rhodz_slow(:,:), &! for HEVI time scheme 
+       f_dPhi_slow(:,:), f_dPhi_fast(:,:), &                          ! geopotential tendencies (NH)
+       f_dW_slow(:,:), f_dW_fast(:,:)                                 ! vertical momentum tendencies (NH)
 
   INTEGER, PARAMETER, PUBLIC :: explicit=1, hevi=2, euler=1, rk4=2, mlf=3, rk25=4, ark23=6, ark33=7
Index: codes/icosagcm/trunk/src/hevi_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 365)
+++ codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 366)
@@ -28,11 +28,13 @@
   SUBROUTINE set_coefs_ark33(dt)
     ! Fully-explicit RK3 scheme disguised as ARK
-    ! To check correctness of caldyn_hevi
     REAL(rstd) :: dt
-    REAL(rstd), PARAMETER, DIMENSION(3,4) :: &
-         ajl=(/ zero, (/.5,0.,0./), (/-1.,2.,0./), (/1./6.,2./3.,1./6./) /)
-    CALL set_coefs_hevi(dt, ajl, ajl)
+    CALL set_coefs_rk(dt, (/ zero, (/.5,0.,0./), (/-1.,2.,0./), (/1./6.,2./3.,1./6./) /) )
   END SUBROUTINE set_coefs_ark33
     
+  SUBROUTINE set_coefs_rk(dt, ajl)
+    REAL(rstd) :: dt, ajl(3,4)
+    CALL set_coefs_hevi(dt,ajl,ajl)
+  END SUBROUTINE set_coefs_rk
+
   SUBROUTINE set_coefs_hevi(dt, ajl_slow, ajl_fast)
     REAL(rstd) :: dt, ajl_slow(3,4), ajl_fast(3,4) ! fast/slow Butcher tableaus
@@ -58,7 +60,9 @@
        CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), &
             f_phis, f_ps,f_mass,f_theta_rhodz,f_u,f_q, &
-            f_geopot, f_hflux, f_wflux, &
+            f_W, f_geopot, f_hflux, f_wflux, &
             f_dps_slow(:,j), f_dmass_slow(:,j), f_dtheta_rhodz_slow(:,j), &
-            f_du_slow(:,j), f_du_fast(:,j) )
+            f_du_slow(:,j), f_du_fast(:,j), &
+            f_dPhi_slow(:,j), f_dPhi_fast(:,j), &
+            f_dW_slow(:,j), f_dW_fast(:,j) )
        ! accumulate mass fluxes for transport scheme
        DO ind=1,ndomain
@@ -79,4 +83,10 @@
           CALL update(bjl(l,j), f_u, f_du_slow(:,l))
           CALL update(cjl(l,j), f_u, f_du_fast(:,l))
+          IF(.NOT. hydrostatic) THEN
+ !            CALL update(bjl(l,j), f_W, f_dW_slow(:,l)) ! slow tendencies of w, Phi not implemented yet
+             CALL update(cjl(l,j), f_W, f_dW_fast(:,l))
+ !            CALL update(bjl(l,j), f_geopot, f_dPhi_slow(:,l))
+             CALL update(cjl(l,j), f_geopot, f_dPhi_fast(:,l))
+          END IF
        END DO
     END DO
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 365)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 366)
@@ -76,5 +76,4 @@
     CALL allocate_field(f_phis,field_t,type_real,name='phis')
     ! Model state at current time step
-    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
     CALL allocate_field(f_ps,field_t,type_real, name='ps')
     CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
@@ -82,4 +81,6 @@
     CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
     CALL allocate_field(f_u,field_u,type_real,llm,name='u')
+    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
+    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W')
     CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
     ! Mass fluxes
@@ -108,4 +109,8 @@
        CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow')
        CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast')
+       CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
+       CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
+       CALL allocate_fields(nb_stage,f_dPhi_slow, field_t,type_real,llm+1,name='dPhi_slow')
+       CALL allocate_fields(nb_stage,f_dPhi_fast, field_t,type_real,llm+1,name='dPhi_fast')
        f_dps => f_dps_slow(:,1)
        f_du => f_du_slow(:,1)
@@ -131,5 +136,5 @@
     CALL init_check_conserve
 
-    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
+    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
 
     CALL transfert_request(f_phis,req_i0) 
