Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 159)
+++ 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
