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
