Index: codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90	(revision 683)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_wflux.k90	(revision 686)
@@ -7,6 +7,5 @@
          convm(l,ij) = convm(l,ij) + convm(l+1,ij)
       END DO
-      l=1
-      dmass_col(ij) = convm(l,ij)
+      dmass_col(ij) = convm(1,ij)
       DO l = 2,llm
          ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
Index: codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 683)
+++ codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 686)
@@ -52,10 +52,9 @@
       DO ij=1,primal_num
          ! Compute residual R_il and B_il
-         l=1
          ! bottom interface l=1
-         ml_g2 = gm2*m_il(l,ij)
-         B_il(l,ij) = A_ik(l,ij) + ml_g2 + tau2_g*rho_bot
-         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
-         + tau2_g*( p_ik(l,ij)-pbot+rho_bot*(Phi_il(l,ij)-PHI_BOT(ij)) )
+         ml_g2 = gm2*m_il(1,ij)
+         B_il(1,ij) = A_ik(1,ij) + ml_g2 + tau2_g*rho_bot
+         R_il(1,ij) = ml_g2*( Phi_il(1,ij)-Phi_star_il(1,ij)) &
+         + tau2_g*( p_ik(1,ij)-pbot+rho_bot*(Phi_il(1,ij)-PHI_BOT(ij)) )
          DO l = 2,llm
             ! inner interfaces
@@ -68,18 +67,16 @@
             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
          END DO
-         l=llm+1
          ! top interface l=llm+1
-         ml_g2 = gm2*m_il(l,ij)
-         B_il(l,ij) = A_ik(l-1,ij) + ml_g2
-         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
-         + tau2_g*( ptop-p_ik(l-1,ij) )
+         ml_g2 = gm2*m_il(llm+1,ij)
+         B_il(llm+1,ij) = A_ik(llm+1 -1,ij) + ml_g2
+         R_il(llm+1,ij) = ml_g2*( Phi_il(llm+1,ij)-Phi_star_il(llm+1,ij)) &
+         + tau2_g*( ptop-p_ik(llm+1 -1,ij) )
          !
          ! Forward sweep :
          ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
          ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
-         l=1
-         X_ij = 1./B_il(l,ij)
-         C_ik(l,ij) = -A_ik(l,ij) * X_ij
-         D_il(l,ij) = R_il(l,ij) * X_ij
+         X_ij = 1./B_il(1,ij)
+         C_ik(1,ij) = -A_ik(1,ij) * X_ij
+         D_il(1,ij) = R_il(1,ij) * X_ij
          DO l = 2,llm
             X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
@@ -87,13 +84,12 @@
             D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
          END DO
-         l=llm+1
-         X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
-         D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
+         X_ij = 1./( B_il(llm+1,ij) + A_ik(llm+1 -1,ij)*C_ik(llm+1 -1,ij) )
+         D_il(llm+1,ij) = (R_il(llm+1,ij)+A_ik(llm+1 -1,ij)*D_il(llm+1 -1,ij)) * X_ij
          ! Back substitution :
          ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
          ! + Newton-Raphson update
          ! top interface l=llm+1
-         x_il(l,ij) = D_il(l,ij)
-         Phi_il(l,ij) = Phi_il(l,ij) - x_il(l,ij)
+         x_il(llm+1,ij) = D_il(llm+1,ij)
+         Phi_il(llm+1,ij) = Phi_il(llm+1,ij) - x_il(llm+1,ij)
          DO l = llm,1,-1
             ! Back substitution at lower interfaces
Index: codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90	(revision 683)
+++ codes/icosagcm/devel/src/kernels_unst/compute_geopot.k90	(revision 686)
@@ -6,12 +6,10 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* theta(l,ij,1)*rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* theta(llm,ij,1)*rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( theta(l,ij,1)*rhodz(l,ij) + theta(l+1,ij,1)*rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* theta(l,ij,1)*rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* theta(1,ij,1)*rhodz(1,ij)
          END IF
       END DO
@@ -29,12 +27,10 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
          END IF
       END DO
@@ -54,12 +50,10 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
          END IF
       END DO
@@ -79,12 +73,10 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)*(1.+theta(llm,ij,2))
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij)*(1.+theta(l,ij,2)) + rhodz(l+1,ij)*(1.+theta(l+1,ij,2)) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)*(1.+theta(1,ij,2))
          END IF
       END DO
Index: codes/icosagcm/devel/src/kernels_unst/remap_eta.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/remap_eta.k90	(revision 686)
+++ codes/icosagcm/devel/src/kernels_unst/remap_eta.k90	(revision 686)
@@ -0,0 +1,34 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- remap_eta ----------------------------------
+   ! IN : rhodz, mass_dak, mass_dbk
+   ! TMP : mass_col, cur_lev, new_rhodz_cum
+   ! OUT : rhodz, old_rhodz, eta
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij=1,primal_num
+      rhodz_cum(1,ij)=0.
+      cur_lev(ij)=1
+      eta(1,ij)=1.
+      new_rhodz_cum(ij)=0.
+      DO l = 1,llm
+         rhodz_cum(l+1,ij) = rhodz_cum(l,ij) + rhodz(l,ij)
+      END DO
+      mass_col(ij) = rhodz_cum(llm+1,ij)
+      DO l = 1,llm
+         old_rhodz(l,ij) = rhodz(l,ij)
+         rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij)
+         rhodz_cum_target = new_rhodz_cum(ij) + rhodz(l,ij)
+         DO level = cur_lev(ij), llm
+            rhodz_cum_levp1 = rhodz_cum(level+1,ij)
+            IF(rhodz_cum_target<=rhodz_cum_levp1) EXIT
+         END DO
+         IF(level>llm) level=llm
+         rhodz_cum_lev = rhodz_cum(level,ij)
+         ! now rhodz_cum_lev <= rhodz_cum_target <= rhodz_cum_levp1
+         cur_lev(ij) = level
+         new_rhodz_cum(ij) = rhodz_cum_target
+         eta(l+1,ij) = level + (rhodz_cum_target-rhodz_cum_lev)/(rhodz_cum_levp1-rhodz_cum_lev)
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- remap_eta ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/remap_theta.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/remap_theta.k90	(revision 686)
+++ codes/icosagcm/devel/src/kernels_unst/remap_theta.k90	(revision 686)
@@ -0,0 +1,24 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- remap_theta ----------------------------------
+   ! IN : thetarhodz, eta
+   ! TMP : thetarhodz_cum, new_thetarhodz_cum
+   ! OUT : thetarhodz
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij=1,primal_num
+      thetarhodz_cum(1,ij)=0.
+      DO l = 1,llm
+         thetarhodz_cum(l+1,ij) = thetarhodz_cum(l,ij) + thetarhodz(l,ij)
+      END DO
+      DO l = 1,llm+1
+         X = eta(l,ij)
+         level = MIN(llm,FLOOR(X)) ! eta=llm+1 => level=llm, X=1
+         X = X-level
+         new_thetarhodz_cum(l,ij) = thetarhodz_cum(level,ij)+X*thetarhodz(level,ij)
+      END DO
+      DO l = 1,llm
+         thetarhodz(l,ij) = new_thetarhodz_cum(l+1,ij) - new_thetarhodz_cum(l,ij)
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- remap_theta ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/remap_u.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/remap_u.k90	(revision 686)
+++ codes/icosagcm/devel/src/kernels_unst/remap_u.k90	(revision 686)
@@ -0,0 +1,27 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- remap_u ----------------------------------
+   ! IN : u, old_rhodz, rhodz, eta
+   ! TMP : urhodz_cum, new_urhodz_cum
+   ! OUT : u
+   !$OMP DO SCHEDULE(STATIC)
+   DO edge=1,edge_num
+      ij_left = left(edge)
+      ij_right = right(edge)
+      urhodz_cum(1,edge)=0.
+      DO l = 1,llm
+         urhodz(l,edge) = u(l,edge)*(old_rhodz(l,ij_left)+old_rhodz(l,ij_right))
+         urhodz_cum(l+1,edge) = urhodz_cum(l,edge) + urhodz(l,edge)
+      END DO
+      DO l = 1,llm+1
+         X = .5*(eta(l,ij_left)+eta(l,ij_right))
+         level = MIN(llm,FLOOR(X))
+         X = X-level
+         new_urhodz_cum(l,edge) = urhodz_cum(level,edge)+X*urhodz(level,edge)
+      END DO
+      DO l = 1,llm
+         u(l,edge) = (new_urhodz_cum(l+1,edge) - new_urhodz_cum(l,edge)) / (rhodz(l,ij_left)+rhodz(l,ij_right))
+      END DO
+   END DO
+   !$OMP END DO
+   !---------------------------- remap_u ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/theta.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 683)
+++ codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 686)
@@ -7,5 +7,4 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=0
          mass_col(ij)=0.
          DO l = 1,llm
