Index: codes/icosagcm/trunk/src/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn.f90	(revision 156)
+++ codes/icosagcm/trunk/src/caldyn.f90	(revision 157)
@@ -1,8 +1,9 @@
 MODULE caldyn_mod
   USE icosa
+  USE caldyn_gcm_mod, ONLY : caldyn_BC
   PRIVATE
   CHARACTER(LEN=255),SAVE :: caldyn_type
   
-  PUBLIC init_caldyn, caldyn
+  PUBLIC init_caldyn, caldyn, caldyn_BC
   
 CONTAINS
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 156)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 157)
@@ -5,5 +5,5 @@
 
   INTEGER, PARAMETER :: energy=1, enstrophy=2
-  TYPE(t_field),POINTER :: f_out_u(:), f_p(:), f_rhodz(:), f_qu(:)
+  TYPE(t_field),POINTER :: f_out_u(:), f_rhodz(:), f_qu(:)
   REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:)
 
@@ -14,6 +14,4 @@
   TYPE(t_field),POINTER :: f_theta(:)
   TYPE(t_field),POINTER :: f_pk(:)
-!  TYPE(t_field),POINTER :: f_pks(:)
-  TYPE(t_field),POINTER :: f_phi(:)
   TYPE(t_field),POINTER :: f_geopot(:)
   TYPE(t_field),POINTER :: f_divm(:)
@@ -24,5 +22,5 @@
   TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu
 
-  PUBLIC init_caldyn, caldyn, write_output_fields,req_ps
+  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields,req_ps
 
 CONTAINS
@@ -48,28 +46,4 @@
     IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
 
-    def='direct'
-    CALL getin('caldyn_exner',def)
-    SELECT CASE(TRIM(def))
-    CASE('lmdz')
-       caldyn_exner=lmdz
-    CASE('direct')
-       caldyn_exner=direct
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>'
-       STOP
-    END SELECT
-
-    def='direct'
-    CALL getin('caldyn_hydrostat',def)
-    SELECT CASE(TRIM(def))
-    CASE('lmdz')
-       caldyn_hydrostat=lmdz
-    CASE('direct')
-       caldyn_hydrostat=direct
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>'
-       STOP
-    END SELECT
-    
     CALL allocate_caldyn
   
@@ -81,5 +55,4 @@
 
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
-    CALL allocate_field(f_p,field_t,type_real,llm+1)
     CALL allocate_field(f_rhodz,field_t,type_real,llm) 
     CALL allocate_field(f_qu,field_u,type_real,llm) 
@@ -93,13 +66,56 @@
     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_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_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_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) 
 
   END SUBROUTINE allocate_caldyn
+
+  SUBROUTINE caldyn_BC(f_phis, f_wflux)
+    USE icosa
+    USE mpipara
+    USE omp_para
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_wflux(:)
+    REAL(rstd),POINTER  :: phis(:)
+    REAL(rstd),POINTER  :: wflux(:,:)
+    REAL(rstd),POINTER  :: geopot(:,:)
+    REAL(rstd),POINTER  :: wwuu(:,:)
+
+    INTEGER :: ind,i,j,ij,l
+
+    IF (omp_first) THEN
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          geopot=f_geopot(ind)
+          phis=f_phis(ind)
+          wflux=f_wflux(ind)
+          wwuu=f_wwuu(ind)
+          
+          DO j=jj_begin-1,jj_end+1
+             DO i=ii_begin-1,ii_end+1
+                ij=(j-1)*iim+i
+                ! lower BCs : geopot=phis, wflux=0, wwuu=0
+                geopot(ij,1) = phis(ij)
+                wflux(ij,1) = 0.
+                wwuu(ij+u_right,1)=0   
+                wwuu(ij+u_lup,1)=0   
+                wwuu(ij+u_ldown,1)=0
+                ! top BCs : wflux=0, wwuu=0
+                wflux(ij,llm+1)  = 0.
+                wwuu(ij+u_right,llm+1)=0   
+                wwuu(ij+u_lup,llm+1)=0   
+                wwuu(ij+u_ldown,llm+1)=0
+             ENDDO
+          ENDDO
+       END DO
+    ENDIF
+
+    !$OMP BARRIER
+  END SUBROUTINE caldyn_BC
    
   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
@@ -127,19 +143,16 @@
     REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:), qu(:,:)
 
 ! temporary shared variable
     REAL(rstd),POINTER  :: theta(:,:)  
-    REAL(rstd),POINTER  :: pk(:,:) !, pks(:)
-    REAL(rstd),POINTER  :: phi(:,:)    
+    REAL(rstd),POINTER  :: pk(:,:)
     REAL(rstd),POINTER  :: geopot(:,:)
     REAL(rstd),POINTER  :: divm(:,:) 
     REAL(rstd),POINTER  :: wwuu(:,:)
-    
-    
-    INTEGER :: ind,ij
+        
+    INTEGER :: ind
     LOGICAL,SAVE :: first=.TRUE.
 !$OMP THREADPRIVATE(first)
-
     
     IF (first) THEN
@@ -163,11 +176,10 @@
           CALL swap_geometry(ind)
           ps=f_ps(ind)
+          u=f_u(ind)
+          theta_rhodz = f_theta_rhodz(ind)
           rhodz=f_rhodz(ind)
-          p=f_p(ind)
+          theta = f_theta(ind)
           qu=f_qu(ind)
-          u=f_u(ind)
-
-          CALL compute_pvort(ps, u, p,rhodz,qu)
-
+          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
        ENDDO
 
@@ -177,28 +189,23 @@
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
+          ps=f_ps(ind)
+          u=f_u(ind)
+          theta_rhodz=f_theta_rhodz(ind)
+          rhodz=f_rhodz(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)
           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)
           wflux=f_wflux(ind)
-          ps=f_ps(ind)
+          wwuu=f_wwuu(ind)
           dps=f_dps(ind)
-          theta_rhodz=f_theta_rhodz(ind)
-          dtheta_rhodz=f_dtheta_rhodz(ind)
-          rhodz=f_rhodz(ind)
-          p=f_p(ind)
-          qu=f_qu(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          out_u=f_out_u(ind)
-          theta = f_theta(ind)
-          pk = f_pk(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, geopot, divm, wwuu)
-!                              theta,pk, pks, phi, geopot, divm, wwuu)
+          CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)
        ENDDO       
        
@@ -207,24 +214,24 @@
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
+          ps=f_ps(ind)
+          u=f_u(ind)
+          theta_rhodz=f_theta_rhodz(ind)
+          rhodz=f_rhodz(ind)
+          theta = f_theta(ind)
+          qu=f_qu(ind)
+          CALL compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
           phis=f_phis(ind)
-          phi=f_phi(ind)
+          pk = f_pk(ind)
           geopot = f_geopot(ind)  
-          ps=f_ps(ind)
+          CALL compute_geopot(ps,phis,rhodz,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)
+          wflux=f_wflux(ind)
+          wwuu=f_wwuu(ind)
           dps=f_dps(ind)
-          hflux=f_hflux(ind)
-          wflux=f_wflux(ind)
-          theta_rhodz=f_theta_rhodz(ind)
-          dtheta_rhodz=f_dtheta_rhodz(ind)
-          rhodz=f_rhodz(ind)
-          p=f_p(ind)
-          qu=f_qu(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          out_u=f_out_u(ind) 
-          wwuu=f_wwuu(ind)
-          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, phi, divm, wwuu)
-!                              theta,pk, pks, phi, geopot, divm, wwuu)
+          CALL compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps, dtheta_rhodz, du)
        ENDDO
        
@@ -254,5 +261,5 @@
 END SUBROUTINE caldyn
     
-SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
+SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu)
   USE icosa
   USE disvert_mod
@@ -263,6 +270,7 @@
   REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
   REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
+  REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
   
@@ -274,20 +282,8 @@
   CALL trace_start("compute_pvort")  
 
-  CALL wait_message(req_ps)   
-!!! Compute pressure
-   DO    l    = ll_begin, ll_endp1
-    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
-           p(ij,l) = ap(l) + bp(l) * ps(ij)
-        ENDDO
-     ENDDO
-  ENDDO
-
-!$OMP BARRIER  
-
-!!! Compute mass
+  CALL wait_message(req_ps)  
+  CALL wait_message(req_theta_rhodz) 
+
+!!! Compute mass & theta
    DO l = ll_begin,ll_end
      CALL test_message(req_u) 
@@ -295,5 +291,6 @@
         DO i=ii_begin-1,ii_end+1
            ij=(j-1)*iim+i
-           rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
+           rhodz(ij,l) = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
+           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
         ENDDO
      ENDDO
@@ -304,5 +301,4 @@
 !!! Compute shallow-water potential vorticity
   DO l = ll_begin,ll_end
-    CALL test_message(req_theta_rhodz) 
 
     DO j=jj_begin-1,jj_end+1
@@ -347,9 +343,6 @@
                                                        
   END SUBROUTINE compute_pvort
-
-  
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
-                            theta,pk, geopot, divm, wwuu) 
-!                            theta,pk, pks, phi, geopot, divm, wwuu) 
+  
+  SUBROUTINE compute_geopot(ps,phis,rhodz,theta,pk,geopot) 
   USE icosa
   USE disvert_mod
@@ -359,44 +352,63 @@
   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
+    
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: p_ik, exner_ik
+
+    CALL trace_start("compute_geopot")
+
+!!! Compute exner function and geopotential
+    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
+
+  CALL trace_end("compute_geopot")
+
+  END SUBROUTINE compute_geopot
+
+  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,divm, dtheta_rhodz, du)
+  USE icosa
+  USE disvert_mod
+  USE exner_mod
+  USE trace
+  USE omp_para
+  IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
-
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
+    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm) ! Exner function
+    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
+
+    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) :: dtheta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
-
-! 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) :: 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),INTENT(OUT) :: du(iim*3*jjm,llm)
+
     REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
     
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: ww,uu, p_ik, exner_ik
-
-    CALL trace_start("compute_caldyn")
+    REAL(rstd) :: ww,uu
+
+    CALL trace_start("compute_caldyn_horiz")
 
     CALL wait_message(req_theta_rhodz) 
-
-!!! Compute theta
-   DO l = ll_begin, ll_end
-      IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-      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
 
   DO l = ll_begin, ll_end
@@ -438,59 +450,9 @@
       ENDDO
     ENDDO
-  ENDDO
-
-!$OMP BARRIER   
-!!! cumulate mass flux divergence from top to bottom
-  DO  l = llm-1, 1, -1
-    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-!$OMP DO SCHEDULE(STATIC) 
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
-      ENDDO
-    ENDDO
-  ENDDO
-
-! IMPLICIT FLUSH on divm
-!!!!!!!!!!!!!!!!!!!!!!!!!
-  
-!!! Compute vertical mass flux
-!  DO l = 2,llm
-  DO l=ll_beginp1,ll_end
-    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ! 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 )
-      ENDDO
-    ENDDO
-  ENDDO
-
-! compute dps, set vertical mass flux at the surface to 0
-  IF (omp_first) THEN
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        wflux(ij,1)  = 0.
-        ! dps/dt = -int(div flux)dz
-        dps(ij)=-divm(ij,1) * g 
-      ENDDO
-    ENDDO
-  ENDIF
-
-
-  IF (omp_last) THEN
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        wflux(ij,llm+1)  = 0.
-      ENDDO
-    ENDDO
-  ENDIF
+
+ END DO
+
 !!! Compute potential vorticity (Coriolis) contribution to du
-
+ 
     SELECT CASE(caldyn_conserv)
     CASE(energy) ! energy-conserving TRiSK
@@ -596,91 +558,4 @@
     END SELECT
   
-!!!! 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
-             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
-
 !!! Compute bernouilli term = Kinetic Energy + geopotential
    DO l=ll_begin,ll_end
@@ -702,5 +577,5 @@
 
  
-!!! gradients of Bernoulli and Exner functions
+!!! Add gradients of Bernoulli and Exner functions to du
 
    DO l=ll_begin,ll_end
@@ -729,5 +604,72 @@
   ENDDO
 
-   
+  CALL trace_end("compute_caldyn_horiz")
+
+END SUBROUTINE compute_caldyn_horiz
+
+SUBROUTINE compute_caldyn_vert(u,theta,rhodz,divm, wflux,wwuu, dps,dtheta_rhodz,du)
+  USE icosa
+  USE disvert_mod
+  USE exner_mod
+  USE trace
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
+    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(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
+    REAL(rstd),INTENT(OUT) :: wwuu(iim*3*jjm,llm+1)
+    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+
+! temporary variable    
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: p_ik, exner_ik
+
+  CALL trace_start("compute_caldyn_vert")
+
+!$OMP BARRIER   
+!!! cumulate mass flux divergence from top to bottom
+  DO  l = llm-1, 1, -1
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+!$OMP DO SCHEDULE(STATIC) 
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
+      ENDDO
+    ENDDO
+  ENDDO
+
+! IMPLICIT FLUSH on divm
+!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! compute dps
+  IF (omp_first) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        ! dps/dt = -int(div flux)dz
+        dps(ij)=-divm(ij,1) * g 
+      ENDDO
+    ENDDO
+  ENDIF
+  
+!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
+  DO l=ll_beginp1,ll_end
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        ! 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 )
+      ENDDO
+    ENDDO
+  ENDDO
+
   DO l=ll_begin,ll_endm1
     DO j=jj_begin,jj_end
@@ -748,26 +690,5 @@
   ENDDO
   
-  IF (omp_first) THEN
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        wwuu(ij+u_right,1)=0   
-        wwuu(ij+u_lup,1)=0   
-        wwuu(ij+u_ldown,1)=0
-      ENDDO
-    ENDDO
-  ENDIF   
-
-  IF (omp_last) THEN
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        wwuu(ij+u_right,llm+1)=0   
-        wwuu(ij+u_lup,llm+1)=0   
-        wwuu(ij+u_ldown,llm+1)=0
-      ENDDO
-    ENDDO
-  ENDIF   
-  
+! Compute vertical transport
   DO l=ll_beginp1,ll_end
     DO j=jj_begin,jj_end
@@ -784,4 +705,5 @@
  !$OMP BARRIER
 
+! Add vertical transport to du
   DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
@@ -795,7 +717,7 @@
   ENDDO      
   
-  CALL trace_end("compute_caldyn")
-
-  END SUBROUTINE compute_caldyn
+  CALL trace_end("compute_caldyn_vert")
+
+  END SUBROUTINE compute_caldyn_vert
 
 !-------------------------------- Diagnostics ----------------------------
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 156)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 157)
@@ -199,4 +199,6 @@
     LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
     LOGICAL, PARAMETER :: check=.FALSE.
+
+    CALL caldyn_BC(f_phis, f_wflux)
 
 !$OMP BARRIER
