Index: codes/icosagcm/trunk/src/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn.f90	(revision 158)
+++ codes/icosagcm/trunk/src/caldyn.f90	(revision 159)
@@ -31,5 +31,5 @@
   END SUBROUTINE init_caldyn
   
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+  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)
   USE icosa
@@ -40,4 +40,5 @@
   TYPE(t_field),POINTER :: f_phis(:)
   TYPE(t_field),POINTER :: f_ps(:)
+  TYPE(t_field),POINTER :: f_mass(:)
   TYPE(t_field),POINTER :: f_theta_rhodz(:)
   TYPE(t_field),POINTER :: f_u(:)
@@ -51,5 +52,5 @@
     SELECT CASE (TRIM(caldyn_type))
       CASE('gcm')
-        CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+        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)
       CASE('adv')
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 158)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 159)
@@ -5,6 +5,6 @@
 
   INTEGER, PARAMETER :: energy=1, enstrophy=2
-  TYPE(t_field),POINTER :: f_out_u(:), f_rhodz(:), f_qu(:)
-  REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:)
+  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:)
+  REAL(rstd),POINTER :: out_u(:,:), p(:,:), qu(:,:)
 
   TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
@@ -18,9 +18,10 @@
   TYPE(t_field),POINTER :: f_wwuu(:)
    
-  INTEGER :: caldyn_hydrostat, caldyn_conserv
+  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
+  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps, &
+       caldyn_eta, eta_mass, eta_lag
 
 CONTAINS
@@ -41,5 +42,6 @@
        caldyn_conserv=enstrophy
     CASE DEFAULT
-       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>'
+       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
+            TRIM(def),'> options are <energy>, <enstrophy>'
        STOP
     END SELECT
@@ -55,20 +57,19 @@
 
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
-    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
     CALL allocate_field(f_qu,field_u,type_real,llm) 
   
-    CALL allocate_field(f_buf_i,field_t,type_real,llm)
-    CALL allocate_field(f_buf_p,field_t,type_real,llm+1) 
-    CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
+    CALL allocate_field(f_buf_i,   field_t,type_real,llm)
+    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
+    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
     CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
     CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
-    CALL allocate_field(f_buf_v,field_z,type_real,llm)
-    CALL allocate_field(f_buf_s,field_t,type_real)
-
-    CALL allocate_field(f_theta,field_t,type_real,llm)     ! potential temperature
-    CALL allocate_field(f_pk,field_t,type_real,llm)
-    CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential
-    CALL allocate_field(f_divm,field_t,type_real,llm)      ! mass flux divergence
-    CALL allocate_field(f_wwuu,field_u,type_real,llm+1) 
+    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
+    CALL allocate_field(f_buf_s,   field_t,type_real)
+
+    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
+    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')
 
   END SUBROUTINE allocate_caldyn
@@ -119,5 +120,5 @@
   END SUBROUTINE caldyn_BC
    
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
+  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)
     USE icosa
@@ -132,4 +133,5 @@
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
@@ -140,8 +142,8 @@
     TYPE(t_field),POINTER :: f_du(:)
     
-    REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: ps(:), dps(:)
+    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: rhodz(:,:), qu(:,:)
+    REAL(rstd),POINTER :: qu(:,:)
 
 ! temporary shared variable
@@ -178,8 +180,8 @@
           u=f_u(ind)
           theta_rhodz = f_theta_rhodz(ind)
-          rhodz=f_rhodz(ind)
+          mass=f_mass(ind)
           theta = f_theta(ind)
           qu=f_qu(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
+          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu)
        ENDDO
 
@@ -192,20 +194,19 @@
           u=f_u(ind)
           theta_rhodz=f_theta_rhodz(ind)
-          rhodz=f_rhodz(ind)
+          mass=f_mass(ind)
           theta = f_theta(ind)
           qu=f_qu(ind)
-          phis=f_phis(ind)
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot)
+          CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
           divm = f_divm(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)
+          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,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)
+          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du)
        ENDDO       
        
@@ -217,21 +218,20 @@
           u=f_u(ind)
           theta_rhodz=f_theta_rhodz(ind)
-          rhodz=f_rhodz(ind)
+          mass=f_mass(ind)
           theta = f_theta(ind)
           qu=f_qu(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
-          phis=f_phis(ind)
+          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu)
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          CALL compute_geopot(ps,phis,rhodz,theta, pk,geopot)
+          CALL compute_geopot(ps,mass,theta, pk,geopot)
           hflux=f_hflux(ind)
           divm = f_divm(ind)
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm,dtheta_rhodz,du)
+          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,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)
+          CALL compute_caldyn_vert(u,theta,mass,divm, wflux,wwuu, dps, dtheta_rhodz, du)
        ENDDO
        
@@ -248,6 +248,8 @@
 !       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
 !            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+
        CALL writefield("ps",f_ps)
        CALL writefield("dps",f_dps)
+       CALL writefield("mass",f_mass)
        CALL writefield("vort",f_qu)
        CALL writefield("theta",f_theta)
@@ -344,5 +346,5 @@
   END SUBROUTINE compute_pvort
   
-  SUBROUTINE compute_geopot(ps,phis,rhodz,theta,pk,geopot) 
+  SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot) 
   USE icosa
   USE disvert_mod
@@ -351,10 +353,9 @@
   USE omp_para
   IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    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(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential
+    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
+    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(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
     
     INTEGER :: i,j,ij,l
@@ -363,19 +364,69 @@
     CALL trace_start("compute_geopot")
 
+    IF(caldyn_eta==eta_mass) THEN
+    
 !!! Compute exner function and geopotential
-    DO l = 1,llm
-!$OMP DO SCHEDULE(STATIC) 
+       DO l = 1,llm
+          !$OMP DO SCHEDULE(STATIC) 
+          DO j=jj_begin-1,jj_end+1
+             DO i=ii_begin-1,ii_end+1
+                ij=(j-1)*iim+i
+                p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
+                !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
+                exner_ik = cpp * (p_ik/preff) ** kappa
+                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 
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ELSE
+
+       ! We are using a Lagrangian vertical coordinate
+       ! Pressure must be 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
+
+       ! uppermost layer
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
              ij=(j-1)*iim+i
-             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
-             exner_ik = cpp * (p_ik/preff) ** kappa
-             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 
+             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
+          END DO
+       END DO
+       ! other layers
+       DO l = llm-1, 1, -1
+          !$OMP DO SCHEDULE(STATIC) 
+          DO j=jj_begin-1,jj_end+1
+             DO i=ii_begin-1,ii_end+1
+                ij=(j-1)*iim+i
+                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
+             END DO
+          END DO
+       END DO
+       ! surface pressure (for diagnostics)
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
+          END DO
+       END DO
+
+       DO l = 1,llm
+          !$OMP DO SCHEDULE(STATIC) 
+          DO j=jj_begin-1,jj_end+1
+             DO i=ii_begin-1,ii_end+1
+                ij=(j-1)*iim+i
+                p_ik = pk(ij,l)
+                exner_ik = cpp * (p_ik/preff) ** kappa
+                ! 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 
+                pk(ij,l) = exner_ik
+             ENDDO
           ENDDO
        ENDDO
-    ENDDO
+
+    END IF
 
   CALL trace_end("compute_geopot")
@@ -798,6 +849,6 @@
     
     ! Temperature
-    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
-    
+!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
+     
     CALL getin('physics',physics_type)
     IF (TRIM(physics_type)=='dcmip') THEN
@@ -828,11 +879,10 @@
     ENDIF
     
-    ! geopotential
+    ! geopotential ! FIXME
     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
     CALL writefield("p",f_buf_p)
-    CALL writefield("phi",f_buf_i)
+    CALL writefield("phi",f_geopot)   ! geopotential
     CALL writefield("theta",f_buf1_i) ! potential temperature
-    CALL writefield("pk",f_buf2_i) ! Exner pressure
-  
+    CALL writefield("pk",f_buf2_i)    ! Exner pressure
   
   END SUBROUTINE write_output_fields
Index: codes/icosagcm/trunk/src/disvert.f90
===================================================================
--- codes/icosagcm/trunk/src/disvert.f90	(revision 158)
+++ codes/icosagcm/trunk/src/disvert.f90	(revision 159)
@@ -5,5 +5,9 @@
   REAL(rstd), SAVE, POINTER :: presnivs(:)
   REAL(rstd), SAVE, POINTER :: mass_al(:), mass_bl(:), mass_ak(:), mass_bk(:), mass_dak(:), mass_dbk(:)
+
   REAL(rstd) :: ptop ! pressure at top of atmosphere l=llm+1
+  ! vertical coordinate : Lagrangian or mass-based
+  INTEGER, SAVE :: caldyn_eta
+  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2
 
 CONTAINS
@@ -91,4 +95,42 @@
   END SUBROUTINE init_disvert  
   
+  SUBROUTINE compute_rhodz(comp, ps, rhodz)
+    USE icosa
+    USE omp_para
+    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
+    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
+    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd) :: m, err
+    INTEGER :: l,i,j,ij,dd
+    err=0.
+
+    dd=0
+
+    DO l = ll_begin, ll_end
+       DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
+             ij=(j-1)*iim+i
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
+             IF(comp) THEN
+                rhodz(ij,l) = m
+             ELSE
+                err = MAX(err,abs(m-rhodz(ij,l)))
+             END IF
+          ENDDO
+       ENDDO
+    ENDDO
+
+    IF(.NOT. comp) THEN
+       IF(err>1e-10) THEN
+          PRINT *, 'Discrepancy between ps and rhodz detected', err
+          STOP
+       ELSE
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
+       END IF
+    END IF
+
+  END SUBROUTINE compute_rhodz
+
+ 
   SUBROUTINE write_apbp
   USE icosa
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 158)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 159)
@@ -4,6 +4,7 @@
 CONTAINS
   
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
     USE icosa
+    USE disvert_mod
     USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0
     USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
@@ -19,9 +20,21 @@
     IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_mass(:)
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
-    
+    REAL(rstd),POINTER :: ps(:), mass(:,:)
+    LOGICAL :: init_mass
+    INTEGER :: ind,i,j,ij,l
+
+    ! most etat0 routines set ps and not mass
+    ! in that case and if caldyn_eta == eta_lag
+    ! 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.
+    ! otherwise mass will be overwritten 
+    init_mass = (caldyn_eta == eta_lag)
+
     etat0_type='jablonowsky06'
     CALL getin("etat0",etat0_type)
@@ -56,6 +69,17 @@
        STOP
     END SELECT
-    
+
+    IF(init_mass) THEN ! initialize mass distribution using ps
+       !$OMP BARRIER
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          mass=f_mass(ind); ps=f_ps(ind)
+          CALL compute_rhodz(.TRUE., ps, mass)
+       END DO
+    END IF
+
   END SUBROUTINE etat0
          
 END MODULE etat0_mod
+
Index: codes/icosagcm/trunk/src/field.f90
===================================================================
--- codes/icosagcm/trunk/src/field.f90	(revision 158)
+++ codes/icosagcm/trunk/src/field.f90	(revision 159)
@@ -63,5 +63,5 @@
          field(ind)%name = name
       ELSE
-         field(ind)%name = '(unkown)'
+         field(ind)%name = '(undefined)'
       END IF
 
@@ -232,5 +232,5 @@
     
     IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN
-       PRINT *, 'get_val_r2d : bad pointer assignation with ' // TRIM(field%name) 
+       PRINT *, 'get_val_r2d : bad pointer assignment with ' // TRIM(field%name) 
        STOP
     END IF
@@ -244,6 +244,7 @@
     
     IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN
-       PRINT *, 'get_val_r3d : bad pointer assignation with ' // TRIM(field%name) 
+       PRINT *, 'get_val_r3d : bad pointer assignment with ' // TRIM(field%name) 
        STOP
+!       CALL ABORT
     END IF
     field_pt=>field%rval3d
@@ -256,5 +257,5 @@
     
     IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN
-       PRINT *, 'get_val_r4d : bad pointer assignation with ' // TRIM(field%name)
+       PRINT *, 'get_val_r4d : bad pointer assignment with ' // TRIM(field%name)
        STOP
     END IF
@@ -268,5 +269,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignation'        
+    IF (field%ndim/=2 .OR. field%data_type/=type_integer) STOP 'get_val_i2d : bad pointer assignment'        
     field_pt=>field%ival2d
   END SUBROUTINE  getval_i2d
@@ -277,5 +278,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignation'        
+    IF (field%ndim/=3 .OR. field%data_type/=type_integer) STOP 'get_val_i3d : bad pointer assignment'        
     field_pt=>field%ival3d
   END SUBROUTINE  getval_i3d
@@ -286,5 +287,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignation'        
+    IF (field%ndim/=4 .OR. field%data_type/=type_integer) STOP 'get_val_i4d : bad pointer assignment'        
     field_pt=>field%ival4d
   END SUBROUTINE  getval_i4d
@@ -295,5 +296,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignation'        
+    IF (field%ndim/=2 .OR. field%data_type/=type_logical) STOP 'get_val_l2d : bad pointer assignment'        
     field_pt=>field%lval2d
   END SUBROUTINE  getval_l2d
@@ -304,5 +305,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignation'        
+    IF (field%ndim/=3 .OR. field%data_type/=type_logical) STOP 'get_val_l3d : bad pointer assignment'        
     field_pt=>field%lval3d
   END SUBROUTINE  getval_l3d
@@ -313,5 +314,5 @@
     TYPE(t_field),INTENT(IN) :: field
     
-    IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignation'        
+    IF (field%ndim/=4 .OR. field%data_type/=type_logical) STOP 'get_val_l4d : bad pointer assignment'        
     field_pt=>field%lval4d
   END SUBROUTINE  getval_l4d    
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 158)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 159)
@@ -11,13 +11,9 @@
   TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0
 
-  TYPE(t_field),POINTER :: f_phis(:)
   TYPE(t_field),POINTER :: f_q(:)
-  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
-  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
-  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
-  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
-  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
-  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
-  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_rhodz(:), f_mass(:), f_massm1(:), f_massm2(:), f_dmass(:)
+  TYPE(t_field),POINTER :: f_phis(:), f_ps(:),f_psm1(:), f_psm2(:), f_dps(:)
+  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:), f_du(:)
+  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:), f_dtheta_rhodz(:)
   TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)  
 
@@ -34,4 +30,5 @@
   USE caldyn_mod
   USE etat0_mod
+  USE disvert_mod
   USE guided_mod
   USE advect_tracer_mod
@@ -45,27 +42,5 @@
   IMPLICIT NONE
 
-    CHARACTER(len=255) :: scheme_name
-
-    CALL allocate_field(f_dps,field_t,type_real)
-    CALL allocate_field(f_du,field_u,type_real,llm)
-    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
-! Model state at current time step (RK/MLF/Euler)
-    CALL allocate_field(f_phis,field_t,type_real)
-    CALL allocate_field(f_ps,field_t,type_real)
-    CALL allocate_field(f_u,field_u,type_real,llm)
-    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
-! Model state at previous time step (RK/MLF)
-    CALL allocate_field(f_psm1,field_t,type_real)
-    CALL allocate_field(f_um1,field_u,type_real,llm)
-    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
-! Tracers
-    CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
-    CALL allocate_field(f_rhodz,field_t,type_real,llm)
-! Mass fluxes
-    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
-    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
-    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
-    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
-
+    CHARACTER(len=255) :: def
 
 !----------------------------------------------------
@@ -122,11 +97,60 @@
 
 
-
-
-
-    scheme_name='runge_kutta'
-    CALL getin('scheme',scheme_name)
-
-    SELECT CASE (TRIM(scheme_name))
+    def='eta_mass'
+    CALL getin('caldyn_eta',def)
+    SELECT CASE(TRIM(def))
+    CASE('eta_mass')
+       caldyn_eta=eta_mass
+    CASE('eta_lag')
+       caldyn_eta=eta_lag
+    CASE DEFAULT
+       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_eta : <', TRIM(def),'> options are <eta_mass>, <eta_lag>'
+       STOP
+    END SELECT
+    IF (is_mpi_root) PRINT *, 'caldyn_eta=',TRIM(def), caldyn_eta, eta_lag, eta_mass
+
+! Time-independant orography
+    CALL allocate_field(f_phis,field_t,type_real,name='phis')
+! Trends
+    CALL allocate_field(f_du,field_u,type_real,llm,name='du')
+    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz')
+! Model state at current time step (RK/MLF/Euler)
+    CALL allocate_field(f_ps,field_t,type_real, name='ps')
+    CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
+    CALL allocate_field(f_u,field_u,type_real,llm,name='u')
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
+! Model state at previous time step (RK/MLF)
+    CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
+    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')
+! Tracers
+    CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
+! Mass fluxes
+    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
+    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
+
+    IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
+       CALL allocate_field(f_dps,field_t,type_real,name='dps')
+       CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
+       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
+    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)
+       ! the following are unused but must point to something
+       f_wfluxt => f_wflux
+       f_dps  => f_phis
+       f_psm1 => f_phis
+    END IF
+
+
+    def='runge_kutta'
+    CALL getin('scheme',def)
+
+    SELECT CASE (TRIM(def))
       CASE('euler')
         scheme=euler
@@ -141,9 +165,18 @@
         nb_stage=matsuno_period+1
      ! Model state 2 time steps ago (MLF)
-        CALL allocate_field(f_psm2,field_t,type_real)
         CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
         CALL allocate_field(f_um2,field_u,type_real,llm)
+        IF(caldyn_eta == eta_mass) THEN ! eta = mass coordinate (default)
+           CALL allocate_field(f_psm2,field_t,type_real)
+           ! the following are unused but must point to something
+           f_massm2 => f_mass
+        ELSE ! eta = Lagrangian vertical coordinate
+           CALL allocate_field(f_massm2,field_t,type_real,llm)
+           ! the following are unused but must point to something
+           f_psm2 => f_phis
+        END IF
+
     CASE default
-       PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
+       PRINT*,'Bad selector for variable scheme : <', TRIM(def),             &
             ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
        STOP
@@ -157,5 +190,5 @@
     CALL init_check_conserve
     CALL init_physics
-    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+    CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
 
     CALL transfert_request(f_phis,req_i0) 
@@ -173,5 +206,7 @@
   USE icosa
   USE dissip_gcm_mod
+  USE disvert_mod
   USE caldyn_mod
+  USE caldyn_gcm_mod, ONLY : req_ps
   USE etat0_mod
   USE guided_mod
@@ -184,21 +219,17 @@
   USE check_conserve_mod
   IMPLICIT NONE  
-    REAL(rstd),POINTER :: phis(:)
     REAL(rstd),POINTER :: q(:,:,:)
-    REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
-    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
-    REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
-    REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
-    REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
-    REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+    REAL(rstd),POINTER :: phis(:), ps(:) ,psm1(:), psm2(:), dps(:)
+    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:), du(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), massm1(:,:), massm2(:,:), dmass(:,:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:), dtheta_rhodz(:,:)
     REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
 
     INTEGER :: ind
     INTEGER :: it,i,j,n,  stage
-    CHARACTER(len=255) :: scheme_name
     LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
     LOGICAL, PARAMETER :: check=.FALSE.
 
-    CALL caldyn_BC(f_phis, f_wflux)
+    CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
 
 !$OMP BARRIER
@@ -234,5 +265,5 @@
     DO stage=1,nb_stage
        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
-            f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
+            f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q, &
             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
        SELECT CASE (scheme)
@@ -333,5 +364,4 @@
 
     SUBROUTINE RK_scheme(stage)
-      USE caldyn_gcm_mod
       IMPLICIT NONE
       INTEGER :: ind, stage
@@ -344,38 +374,37 @@
       tau = dt*coef(stage)
 
-      DO ind=1,ndomain
-         CALL swap_dimensions(ind)
-         CALL swap_geometry(ind)
-         ps=f_ps(ind)   
-         psm1=f_psm1(ind) 
-         dps=f_dps(ind) 
-         
-         IF (stage==1) THEN ! first stage : save model state in XXm1
-           IF (omp_first) THEN
-             DO j=jj_begin,jj_end
-               DO i=ii_begin,ii_end
-                 ij=(j-1)*iim+i
-                 psm1(ij)=ps(ij)
+      ! if mass coordinate, deal with ps first on one core
+      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) 
+               
+               IF (stage==1) THEN ! first stage : save model state in XXm1
+                  DO j=jj_begin,jj_end
+                     DO i=ii_begin,ii_end
+                        ij=(j-1)*iim+i
+                        psm1(ij)=ps(ij)
+                     ENDDO
+                  ENDDO
+               ENDIF
+            
+               ! updates are of the form x1 := x0 + tau*f(x1)
+               DO j=jj_begin,jj_end
+                  DO i=ii_begin,ii_end
+                     ij=(j-1)*iim+i
+                     ps(ij)=psm1(ij)+tau*dps(ij)
+                  ENDDO
                ENDDO
-             ENDDO
-           ENDIF
-         END IF
-         
-         ! updates are of the form x1 := x0 + tau*f(x1)
-           IF (omp_first) THEN
-             DO j=jj_begin,jj_end
-               DO i=ii_begin,ii_end
-                 ij=(j-1)*iim+i
-                 ps(ij)=psm1(ij)+tau*dps(ij)
-               ENDDO
-             ENDDO
-           ENDIF
-       
-       ENDDO
-
-       CALL send_message(f_ps,req_ps)      
-
-
-      
+            ENDDO
+         ENDIF
+   
+         CALL send_message(f_ps,req_ps)
+      END IF
+
+      ! now deal with other prognostic variables
       DO ind=1,ndomain
          CALL swap_dimensions(ind)
@@ -386,6 +415,6 @@
          
          IF (stage==1) THEN ! first stage : save model state in XXm1
-              
-         DO l=ll_begin,ll_end
+            
+           DO l=ll_begin,ll_end
              DO j=jj_begin,jj_end
                DO i=ii_begin,ii_end
@@ -398,4 +427,15 @@
              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
@@ -413,5 +453,15 @@
            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)
@@ -476,45 +526,9 @@
   END SUBROUTINE timeloop    
 
-  SUBROUTINE compute_rhodz(comp, ps, rhodz)
-    USE icosa
-    USE disvert_mod
-    USE omp_para
-    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
-    REAL(rstd), INTENT(IN) :: ps(iim*jjm)
-    REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
-    REAL(rstd) :: m, err
-    INTEGER :: l,i,j,ij,dd
-    err=0.
-
-    dd=0
-
-    DO l = ll_begin, ll_end
-       DO j=jj_begin-dd,jj_end+dd
-          DO i=ii_begin-dd,ii_end+dd
-             ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-             IF(comp) THEN
-                rhodz(ij,l) = m
-             ELSE
-                err = MAX(err,abs(m-rhodz(ij,l)))
-             END IF
-          ENDDO
-       ENDDO
-    ENDDO
-
-    IF(.NOT. comp) THEN
-       IF(err>1e-10) THEN
-          PRINT *, 'Discrepancy between ps and rhodz detected', err
-          STOP
-       ELSE
-!          PRINT *, 'No discrepancy between ps and rhodz detected'
-       END IF
-    END IF
-
-  END SUBROUTINE compute_rhodz
-
-  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
+ SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
     USE icosa
     USE omp_para
+    USE disvert_mod
+    IMPLICIT NONE
     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
     REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
@@ -538,12 +552,14 @@
        ENDDO
 
-       DO l=ll_begin,ll_endp1
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             wfluxt(ij,l) = tau*wflux(ij,l)
-           ENDDO
-         ENDDO
-       ENDDO
+       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
+          DO l=ll_begin,ll_endp1
+             DO j=jj_begin,jj_end
+                DO i=ii_begin,ii_end
+                   ij=(j-1)*iim+i
+                   wfluxt(ij,l) = tau*wflux(ij,l)
+                ENDDO
+             ENDDO
+          ENDDO
+       END IF
     ELSE
 
@@ -559,14 +575,16 @@
        ENDDO
 
-       DO l=ll_begin,ll_endp1
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
-           ENDDO
-         ENDDO
-       ENDDO
-
-    END IF
+       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
+          DO l=ll_begin,ll_endp1
+             DO j=jj_begin,jj_end
+                DO i=ii_begin,ii_end
+                   ij=(j-1)*iim+i
+                   wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
+                ENDDO
+             ENDDO
+          ENDDO
+       END IF
+
+   END IF
 
   END SUBROUTINE accumulate_fluxes
