Index: codes/icosagcm/trunk/src/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn.f90	(revision 161)
+++ codes/icosagcm/trunk/src/caldyn.f90	(revision 162)
@@ -32,5 +32,5 @@
   
   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
-       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
+       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
   USE icosa
   USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn
@@ -47,4 +47,5 @@
   TYPE(t_field),POINTER :: f_wflux(:)
   TYPE(t_field),POINTER :: f_dps(:)
+  TYPE(t_field),POINTER :: f_dmass(:)
   TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
   TYPE(t_field),POINTER :: f_du(:)
@@ -53,5 +54,5 @@
       CASE('gcm')
         CALL caldyn_gcm(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
-             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
+             f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
       CASE('adv')
         CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 161)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 162)
@@ -5,5 +5,5 @@
 
   INTEGER, PARAMETER :: energy=1, enstrophy=2
-  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:)
+  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
   REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:)
 
@@ -15,13 +15,12 @@
   TYPE(t_field),POINTER :: f_pk(:)
   TYPE(t_field),POINTER :: f_geopot(:)
-  TYPE(t_field),POINTER :: f_divm(:)
   TYPE(t_field),POINTER :: f_wwuu(:)
    
   INTEGER :: caldyn_conserv
   
-  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu
-
-  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps, &
-       caldyn_eta, eta_mass, eta_lag
+  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
+
+  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
+       req_ps, req_mass, caldyn_eta, eta_mass, eta_lag
 
 CONTAINS
@@ -58,4 +57,5 @@
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
     CALL allocate_field(f_qu,field_u,type_real,llm) 
+    CALL allocate_field(f_qv,field_z,type_real,llm) 
   
     CALL allocate_field(f_buf_i,   field_t,type_real,llm)
@@ -70,5 +70,4 @@
     CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')  ! geopotential
-    CALL allocate_field(f_divm,  field_t,type_real,llm,  name='divm')    ! mass flux divergence
     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
 
@@ -121,6 +120,7 @@
    
   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
-       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
+       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
     USE icosa
+    USE disvert_mod
     USE vorticity_mod
     USE kinetic_mod
@@ -139,4 +139,5 @@
     TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
     TYPE(t_field),POINTER :: f_dps(:)
+    TYPE(t_field),POINTER :: f_dmass(:)
     TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
     TYPE(t_field),POINTER :: f_du(:)
@@ -146,4 +147,5 @@
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: qu(:,:)
+    REAL(rstd),POINTER :: qv(:,:)
 
 ! temporary shared variable
@@ -151,5 +153,5 @@
     REAL(rstd),POINTER  :: pk(:,:)
     REAL(rstd),POINTER  :: geopot(:,:)
-    REAL(rstd),POINTER  :: divm(:,:) 
+    REAL(rstd),POINTER  :: convm(:,:) 
     REAL(rstd),POINTER  :: wwuu(:,:)
         
@@ -158,11 +160,21 @@
 !$OMP THREADPRIVATE(first)
     
-    IF (first) THEN
+    ! MPI messages need to be sent at first call to caldyn
+    ! This is needed only once : the next ones will be sent by timeloop
+    IF (first) THEN 
       first=.FALSE.
-      CALL init_message(f_ps,req_i1,req_ps)
+      IF(caldyn_eta==eta_mass) THEN
+         CALL init_message(f_ps,req_i1,req_ps)
+      ELSE
+         CALL init_message(f_mass,req_i1,req_mass)
+      END IF
       CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
       CALL init_message(f_u,req_e1_vect,req_u)
       CALL init_message(f_qu,req_e1_scal,req_qu)
-      CALL send_message(f_ps,req_ps) 
+      IF(caldyn_eta==eta_mass) THEN
+         CALL send_message(f_ps,req_ps) 
+      ELSE
+         CALL send_message(f_mass,req_mass) 
+      END IF
     ENDIF
     
@@ -183,5 +195,6 @@
           theta = f_theta(ind)
           qu=f_qu(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu)
+          qv=f_qv(ind)
+          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
        ENDDO
 
@@ -201,12 +214,14 @@
           CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
-          divm = f_divm(ind)
+          convm = f_dmass(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)
-          wflux=f_wflux(ind)
-          wwuu=f_wwuu(ind)
-          dps=f_dps(ind)
-          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du)
+          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
+          IF(caldyn_eta==eta_mass) THEN
+             wflux=f_wflux(ind)
+             wwuu=f_wwuu(ind)
+             dps=f_dps(ind)
+             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+          END IF
        ENDDO       
        
@@ -221,17 +236,19 @@
           theta = f_theta(ind)
           qu=f_qu(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu)
+          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
           CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
-          divm = f_divm(ind)
+          convm = f_dmass(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)
-          wflux=f_wflux(ind)
-          wwuu=f_wwuu(ind)
-          dps=f_dps(ind)
-          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du)
+          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
+          IF(caldyn_eta==eta_mass) THEN
+             wflux=f_wflux(ind)
+             wwuu=f_wwuu(ind)
+             dps=f_dps(ind)
+             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+          END IF
        ENDDO
        
@@ -252,6 +269,9 @@
        CALL writefield("dps",f_dps)
        CALL writefield("mass",f_mass)
+       CALL writefield("dmass",f_dmass)
        CALL writefield("vort",f_qu)
        CALL writefield("theta",f_theta)
+       CALL writefield("exner",f_pk)
+       CALL writefield("pv",f_qv)
 
     END IF
@@ -263,5 +283,5 @@
 END SUBROUTINE caldyn
     
-SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
+SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
   USE icosa
   USE disvert_mod
@@ -273,29 +293,47 @@
   REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
   REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
+  REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
   
   INTEGER :: i,j,ij,l
-  REAL(rstd) :: etav,hv
-  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
+  REAL(rstd) :: etav,hv, m
+!  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
   
 
   CALL trace_start("compute_pvort")  
 
-  CALL wait_message(req_ps)  
+  IF(caldyn_eta==eta_mass) THEN
+     CALL wait_message(req_ps)  
+  ELSE
+     CALL wait_message(req_mass)
+  END IF
   CALL wait_message(req_theta_rhodz) 
 
-!!! Compute mass & theta
-   DO l = ll_begin,ll_end
-     CALL test_message(req_u) 
-     DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-           ij=(j-1)*iim+i
-           rhodz(ij,l) = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
-           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+  IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
+     DO l = ll_begin,ll_end
+        CALL test_message(req_u) 
+        DO j=jj_begin-1,jj_end+1
+           DO i=ii_begin-1,ii_end+1
+              ij=(j-1)*iim+i
+              m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+!              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
+              rhodz(ij,l) = m
+              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+           ENDDO
         ENDDO
      ENDDO
-  ENDDO
+  ELSE ! Compute only theta
+     DO l = ll_begin,ll_end
+        CALL test_message(req_u) 
+        DO j=jj_begin-1,jj_end+1
+           DO i=ii_begin-1,ii_end+1
+              ij=(j-1)*iim+i
+              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+           ENDDO
+        ENDDO
+     ENDDO
+  END IF
 
   CALL wait_message(req_u)   
@@ -377,5 +415,5 @@
                 pk(ij,l) = exner_ik
                 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
              ENDDO
           ENDDO
@@ -413,5 +451,4 @@
           END DO
        END DO
-
        DO l = 1,llm
           !$OMP DO SCHEDULE(STATIC) 
@@ -434,5 +471,5 @@
   END SUBROUTINE compute_geopot
 
-  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm, dtheta_rhodz, du)
+  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
   USE icosa
   USE disvert_mod
@@ -449,5 +486,5 @@
 
     REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
-    REAL(rstd),INTENT(OUT) :: divm(iim*jjm,llm)  ! mass flux divergence
+    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
     REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
@@ -483,6 +520,6 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i  
-        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
+        ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+        convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
                                 ne_rup*hflux(ij+u_rup,l)       +  &  
                                 ne_lup*hflux(ij+u_lup,l)       +  &  
@@ -659,5 +696,5 @@
 END SUBROUTINE compute_caldyn_horiz
 
-SUBROUTINE compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps,dtheta_rhodz,du)
+SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
   USE icosa
   USE disvert_mod
@@ -669,5 +706,5 @@
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(INOUT)  :: divm(iim*jjm,llm)  ! mass flux divergence
+    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
 
     REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
@@ -684,5 +721,5 @@
 
 !$OMP BARRIER   
-!!! cumulate mass flux divergence from top to bottom
+!!! cumulate mass flux convergence from top to bottom
   DO  l = llm-1, 1, -1
     IF (caldyn_conserv==energy) CALL test_message(req_qu) 
@@ -691,10 +728,10 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
+        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
       ENDDO
     ENDDO
   ENDDO
 
-! IMPLICIT FLUSH on divm
+! IMPLICIT FLUSH on convm
 !!!!!!!!!!!!!!!!!!!!!!!!!
 
@@ -705,5 +742,5 @@
         ij=(j-1)*iim+i
         ! dps/dt = -int(div flux)dz
-        dps(ij)=-divm(ij,1) * g 
+        dps(ij) = convm(ij,1) * g 
       ENDDO
     ENDDO
@@ -718,5 +755,5 @@
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
-        wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 )
+        wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
       ENDDO
     ENDDO
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 161)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 162)
@@ -9,5 +9,5 @@
   INTEGER  :: itau_sync=10
 
-  TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0
+  TYPE(t_message) :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0
 
   TYPE(t_field),POINTER :: f_q(:)
@@ -130,4 +130,5 @@
     CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
     CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
+    CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
 
     IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
@@ -136,10 +137,7 @@
        CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')
        ! the following are unused but must point to something
-       f_massm1 => f_mass
-       f_dmass => f_mass
+!       f_massm1 => f_mass
     ELSE ! eta = Lagrangian vertical coordinate
-       CALL allocate_field(f_mass,field_t,type_real,llm)
-       CALL allocate_field(f_massm1,field_t,type_real,llm)
-       CALL allocate_field(f_dmass,field_t,type_real,llm)
+       CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
        ! the following are unused but must point to something
        f_wfluxt => f_wflux
@@ -147,5 +145,4 @@
        f_psm1 => f_phis
     END IF
-
 
     def='runge_kutta'
@@ -197,4 +194,5 @@
 
     CALL init_message(f_ps,req_i0,req_ps0)
+    CALL init_message(f_mass,req_i0,req_mass0)
     CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
     CALL init_message(f_u,req_e0_vect,req_u0)
@@ -208,5 +206,5 @@
   USE disvert_mod
   USE caldyn_mod
-  USE caldyn_gcm_mod, ONLY : req_ps
+  USE caldyn_gcm_mod, ONLY : req_ps, req_mass
   USE etat0_mod
   USE guided_mod
@@ -237,6 +235,10 @@
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
-     rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
+     rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
+     IF(caldyn_eta==eta_mass) THEN
+        CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
+     ELSE
+        rhodz(:,:)=mass(:,:)
+     END IF
   END DO
   fluxt_zero=.TRUE.
@@ -245,8 +247,10 @@
     IF (MOD(it,itau_sync)==0) THEN
       CALL send_message(f_ps,req_ps0)
+      CALL send_message(f_mass,req_mass0)
       CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
       CALL send_message(f_u,req_u0)
       CALL send_message(f_q,req_q0) 
       CALL wait_message(req_ps0)
+      CALL wait_message(req_mass0)
       CALL wait_message(req_theta_rhodz0) 
       CALL wait_message(req_u0)
@@ -266,5 +270,5 @@
        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
             f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
-            f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
+            f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
        SELECT CASE (scheme)
        CASE(euler)
@@ -326,20 +330,32 @@
        CALL swap_dimensions(ind)
        CALL swap_geometry(ind)
-       IF(with_dps) THEN
-         ps=f_ps(ind) ; dps=f_dps(ind) ; 
-
-         IF (omp_first) THEN
-           DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-               ij=(j-1)*iim+i
-               ps(ij)=ps(ij)+dt*dps(ij)
-             ENDDO
-           ENDDO
-         ENDIF
-         
-         hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
-         wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
-         CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
-       END IF
+
+       IF(with_dps) THEN ! update ps/mass
+          IF(caldyn_eta==eta_mass) THEN ! update ps
+             ps=f_ps(ind) ; dps=f_dps(ind) ;              
+             IF (omp_first) THEN
+                DO j=jj_begin,jj_end
+                   DO i=ii_begin,ii_end
+                      ij=(j-1)*iim+i
+                      ps(ij)=ps(ij)+dt*dps(ij)
+                   ENDDO
+                ENDDO
+             ENDIF 
+          ELSE ! update mass
+             mass=f_mass(ind) ; dmass=f_dmass(ind) ;              
+             DO l=1,llm
+                DO j=jj_begin,jj_end
+                   DO i=ii_begin,ii_end
+                      ij=(j-1)*iim+i
+                      mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
+                   ENDDO
+                ENDDO
+             END DO
+          END IF
+
+          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
+          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
+          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind))
+       END IF ! update ps/mass
        
        u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind)
@@ -377,10 +393,9 @@
       IF(caldyn_eta==eta_mass) THEN
          IF (omp_first) THEN
+
             DO ind=1,ndomain
                CALL swap_dimensions(ind)
                CALL swap_geometry(ind)
-               ps=f_ps(ind)   
-               psm1=f_psm1(ind) 
-               dps=f_dps(ind) 
+               ps=f_ps(ind) ; psm1=f_psm1(ind) ; dps=f_dps(ind) 
                
                IF (stage==1) THEN ! first stage : save model state in XXm1
@@ -402,6 +417,35 @@
             ENDDO
          ENDIF
-   
          CALL send_message(f_ps,req_ps)
+      
+      ELSE ! Lagrangian coordinate, deal with mass
+         DO ind=1,ndomain
+            CALL swap_dimensions(ind)
+            CALL swap_geometry(ind)
+            mass=f_mass(ind); dmass=f_dmass(ind); massm1=f_massm1(ind)
+
+            IF (stage==1) THEN ! first stage : save model state in XXm1
+               DO l=ll_begin,ll_end
+                  DO j=jj_begin,jj_end
+                     DO i=ii_begin,ii_end
+                        ij=(j-1)*iim+i
+                        massm1(ij,l)=mass(ij,l)
+                     ENDDO
+                  ENDDO
+               ENDDO
+            END IF
+
+            ! updates are of the form x1 := x0 + tau*f(x1)
+            DO l=ll_begin,ll_end
+               DO j=jj_begin,jj_end
+                  DO i=ii_begin,ii_end
+                     ij=(j-1)*iim+i
+                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
+                  ENDDO
+               ENDDO
+            ENDDO
+         END DO
+         CALL send_message(f_mass,req_mass)
+
       END IF
 
@@ -410,10 +454,10 @@
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
-         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind)
-         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind)
-         dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
+         u=f_u(ind)      ; du=f_du(ind)      ; um1=f_um1(ind) 
+         theta_rhodz=f_theta_rhodz(ind)
+         theta_rhodzm1=f_theta_rhodzm1(ind)
+         dtheta_rhodz=f_dtheta_rhodz(ind)
          
          IF (stage==1) THEN ! first stage : save model state in XXm1
-            
            DO l=ll_begin,ll_end
              DO j=jj_begin,jj_end
@@ -427,19 +471,6 @@
              ENDDO
            ENDDO
-
-           IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable
-              DO l=ll_begin,ll_end
-                 DO j=jj_begin,jj_end
-                    DO i=ii_begin,ii_end
-                       ij=(j-1)*iim+i
-                       massm1(ij,l)=mass(ij,l)
-                    ENDDO
-                 ENDDO
-              ENDDO
-           END IF
-
-         END IF
-         ! updates are of the form x1 := x0 + tau*f(x1)
-         
+         END IF        
+
          DO l=ll_begin,ll_end
            DO j=jj_begin,jj_end
@@ -453,15 +484,5 @@
            ENDDO
          ENDDO
-         IF(caldyn_eta==eta_lag) THEN ! mass = additional prognostic variable
-            DO l=ll_begin,ll_end
-               DO j=jj_begin,jj_end
-                  DO i=ii_begin,ii_end
-                     ij=(j-1)*iim+i
-                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
-                  ENDDO
-               ENDDO
-            ENDDO
-         END IF
-         
+
          IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage
             hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
@@ -562,4 +583,5 @@
           ENDDO
        END IF
+
     ELSE
 
@@ -590,3 +612,34 @@
   END SUBROUTINE accumulate_fluxes
   
+  FUNCTION maxval_i(p)
+    USE icosa
+    IMPLICIT NONE
+    REAL(rstd), DIMENSION(iim*jjm) :: p
+    REAL(rstd) :: maxval_i
+    INTEGER :: j, ij
+    
+    maxval_i=p((jj_begin-1)*iim+ii_begin)
+    
+    DO j=jj_begin-1,jj_end+1
+       ij=(j-1)*iim
+       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
+    END DO
+  END FUNCTION maxval_i
+
+  FUNCTION maxval_ik(p)
+    USE icosa
+    IMPLICIT NONE
+    REAL(rstd) :: p(iim*jjm, llm)
+    REAL(rstd) :: maxval_ik(llm)
+    INTEGER :: l,j, ij
+    
+    DO l=1,llm
+       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
+       DO j=jj_begin-1,jj_end+1
+          ij=(j-1)*iim
+          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
+       END DO
+    END DO
+  END FUNCTION maxval_ik
+
 END MODULE timeloop_gcm_mod
