Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 155)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 156)
@@ -14,6 +14,7 @@
   TYPE(t_field),POINTER :: f_theta(:)
   TYPE(t_field),POINTER :: f_pk(:)
-  TYPE(t_field),POINTER :: f_pks(:)
+!  TYPE(t_field),POINTER :: f_pks(:)
   TYPE(t_field),POINTER :: f_phi(:)
+  TYPE(t_field),POINTER :: f_geopot(:)
   TYPE(t_field),POINTER :: f_divm(:)
   TYPE(t_field),POINTER :: f_wwuu(:)
@@ -94,6 +95,7 @@
     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_pks,field_t,type_real) ! Exner function
-    CALL allocate_field(f_phi,field_t,type_real,llm)    ! geopotential
+!    CALL allocate_field(f_pks,field_t,type_real) ! Exner function
+!    CALL allocate_field(f_phi,field_t,type_real,llm)     ! geopotential
+    CALL allocate_field(f_geopot,field_t,type_real,llm+1)  ! geopotential (new)
     CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence
     CALL allocate_field(f_wwuu,field_u,type_real,llm+1) 
@@ -129,6 +131,7 @@
 ! temporary shared variable
     REAL(rstd),POINTER  :: theta(:,:)  
-    REAL(rstd),POINTER  :: pk(:,:), pks(:)
+    REAL(rstd),POINTER  :: pk(:,:) !, pks(:)
     REAL(rstd),POINTER  :: phi(:,:)    
+    REAL(rstd),POINTER  :: geopot(:,:)
     REAL(rstd),POINTER  :: divm(:,:) 
     REAL(rstd),POINTER  :: wwuu(:,:)
@@ -189,11 +192,13 @@
           theta = f_theta(ind)
           pk = f_pk(ind)
-          pks = f_pks(ind)
-          phi = f_phi(ind)  
+!          pks = f_pks(ind)
+!          phi = f_phi(ind)
+          geopot = f_geopot(ind)  
           divm = f_divm(ind)
           wwuu=f_wwuu(ind)
 
           CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, &
-                              theta,pk, pks, phi, divm, wwuu)
+                              theta,pk, geopot, divm, wwuu)
+!                              theta,pk, pks, phi, geopot, divm, wwuu)
        ENDDO       
        
@@ -203,4 +208,6 @@
           CALL swap_geometry(ind)
           phis=f_phis(ind)
+          phi=f_phi(ind)
+          geopot = f_geopot(ind)  
           ps=f_ps(ind)
           dps=f_dps(ind)
@@ -218,5 +225,6 @@
           CALL compute_pvort(ps, u, p,rhodz,qu)
           CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
-                              theta,pk, pks, phi, divm, wwuu)
+                              theta,pk, phi, divm, wwuu)
+!                              theta,pk, pks, phi, geopot, divm, wwuu)
        ENDDO
        
@@ -342,5 +350,6 @@
   
   SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
-                            theta,pk, pks, phi, divm, wwuu) 
+                            theta,pk, geopot, divm, wwuu) 
+!                            theta,pk, pks, phi, geopot, divm, wwuu) 
   USE icosa
   USE disvert_mod
@@ -364,9 +373,9 @@
 ! temporary variable
     REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm), pks(iim*jjm) ! Exner function
-    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) !, pks(iim*jjm) ! Exner function
+!    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential
+    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)    ! geopotential
     REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence
     REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
-
     
     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
@@ -374,5 +383,5 @@
     
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: ww,uu
+    REAL(rstd) :: ww,uu, p_ik, exner_ik
 
     CALL trace_start("compute_caldyn")
@@ -587,65 +596,91 @@
     END SELECT
   
-!!! Compute Exner function
-!    CALL compute_exner(ps,p,pks,pk,1) 
-! replaced in source
-    IF (omp_first) THEN    
+!!!! Compute Exner function
+!!    CALL compute_exner(ps,p,pks,pk,1) 
+!! replaced in source
+!    IF (omp_first) THEN    
+!       DO j=jj_begin-1,jj_end+1
+!          DO i=ii_begin-1,ii_end+1
+!             ij=(j-1)*iim+i
+!             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+!          ENDDO
+!       ENDDO
+!     ENDIF
+!     
+!       ! 3D : pk
+!     DO l = ll_begin, ll_end
+!        DO j=jj_begin-1,jj_end+1
+!           DO i=ii_begin-1,ii_end+1
+!              ij=(j-1)*iim+i
+!              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+!           ENDDO
+!        ENDDO
+!     ENDDO
+!
+!!---> flush pk,pks, theta
+!!$OMP BARRIER
+!
+!!! Compute geopotential (old)
+!
+!  ! for first layer
+!  IF (omp_first) THEN
+!    DO j=jj_begin-1,jj_end+1
+!      DO i=ii_begin-1,ii_end+1
+!        ij=(j-1)*iim+i
+!        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+!      ENDDO
+!    ENDDO
+!  ENDIF
+!!!-> implicit flush on phi(:,1)
+!          
+!  ! for other layers
+!   DO l = ll_beginp1, ll_end
+!     DO j=jj_begin-1,jj_end+1
+!       DO i=ii_begin-1,ii_end+1
+!         ij=(j-1)*iim+i
+!         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+!                          * (  pk(ij,l-1) -  pk(ij,l)    )
+!       ENDDO
+!     ENDDO
+!   ENDDO
+!
+!!$OMP BARRIER
+!   DO l = 2, llm
+!!$OMP DO
+!     DO j=jj_begin-1,jj_end+1
+!       DO i=ii_begin-1,ii_end+1
+!         ij=(j-1)*iim+i
+!         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
+!       ENDDO
+!     ENDDO
+!   ENDDO
+!! --> IMPLICIT FLUSH on phi
+      
+!!! Compute exner function and geopotential
+
+! geopot=phis for first layer
+!$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
+          geopot(ij,1) = phis(ij)
+       ENDDO
+    ENDDO
+! for other layers
+    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
-             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+             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
-     ENDIF
-     
-       ! 3D : pk
-     DO l = ll_begin, ll_end
-        DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-              ij=(j-1)*iim+i
-              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
-           ENDDO
-        ENDDO
-     ENDDO
-
-!---> flush pk,pks, theta
-!$OMP BARRIER
-
-!!! Compute geopotential
-
-  ! for first layer
-  IF (omp_first) THEN
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
-        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
-      ENDDO
-    ENDDO
-  ENDIF
-!!-> implicit flush on phi(:,1)
-          
-  ! for other layers
-   DO l = ll_beginp1, ll_end
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                          * (  pk(ij,l-1) -  pk(ij,l)    )
-       ENDDO
-     ENDDO
-   ENDDO
-
-!$OMP BARRIER
-   DO l = 2, llm
-!$OMP DO
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
-       ENDDO
-     ENDDO
-   ENDDO
-! --> IMPLICIT FLUSH on phi
-      
+    ENDDO
+
 !!! Compute bernouilli term = Kinetic Energy + geopotential
    DO l=ll_begin,ll_end
@@ -654,5 +689,5 @@
         ij=(j-1)*iim+i
 
-        berni(ij,l) =  phi(ij,l) &                                                          
+        berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
                      + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
Index: codes/icosagcm/trunk/src/disvert.f90
===================================================================
--- codes/icosagcm/trunk/src/disvert.f90	(revision 155)
+++ codes/icosagcm/trunk/src/disvert.f90	(revision 156)
@@ -4,4 +4,6 @@
   REAL(rstd), SAVE, POINTER :: bp(:)
   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
 
 CONTAINS
@@ -17,5 +19,6 @@
   IMPLICIT NONE
     CHARACTER(LEN=255) :: disvert_type = 'std'
-    
+    INTEGER :: l
+
     CALL getin("disvert",disvert_type)
     
@@ -61,4 +64,28 @@
         
     END SELECT
+
+    ! Convert from pressure coordinate to mass coordinate
+    ! pk = ap(k) + bp(k)*ps(ij) = ptop + g*( mass_ak(k) + mass_bk(k)*ms(ij) )
+    ptop = ap(llm+1)
+    ALLOCATE(mass_al(llm+1))
+    ALLOCATE(mass_bl(llm+1))
+    ALLOCATE(mass_ak(llm))
+    ALLOCATE(mass_bk(llm))
+    ALLOCATE(mass_dak(llm))
+    ALLOCATE(mass_dbk(llm))
+
+    ! FIXME : leave ps for the moment ; change ps to ms later
+    DO l=1,llm+1
+!       mass_al(l) = (ap(l)-ptop)/g
+!       mass_bl(l) = bp(l)/g
+       mass_al(l) = (ap(l)-ptop)
+       mass_bl(l) = bp(l)
+    END DO
+    DO l=1,llm
+       mass_ak(l) = .5*(mass_al(l)+mass_al(l+1))
+       mass_bk(l) = .5*(mass_bl(l)+mass_bl(l+1))
+       mass_dak(l) = mass_al(l)-mass_al(l+1)  ! >0
+       mass_dbk(l) = mass_bl(l)-mass_bl(l+1)  ! >0
+    END DO
 
   END SUBROUTINE init_disvert  
