Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 124)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 125)
@@ -2,6 +2,6 @@
   USE icosa
 
-  TYPE(t_field),POINTER :: f_out_u(:)
-  REAL(rstd),POINTER :: out_u(:,:)
+  TYPE(t_field),POINTER :: f_out_u(:), f_p(:), f_rhodz(:), f_qu(:)
+  REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:)
 
   TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
@@ -10,5 +10,5 @@
   PUBLIC init_caldyn, caldyn, write_output_fields
 
-  INTEGER :: caldyn_hydrostat
+  INTEGER :: caldyn_hydrostat, caldyn_conserv
 
 CONTAINS
@@ -20,4 +20,16 @@
     CHARACTER(len=255) :: def
   
+    def='direct'
+    CALL getin('caldyn_conserv',def)
+    SELECT CASE(TRIM(def))
+    CASE('energy')
+       caldyn_conserv=1
+    CASE('direct')
+       caldyn_conserv=2
+    CASE DEFAULT
+       PRINT*,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>'
+       STOP
+    END SELECT
+
     def='direct'
     CALL getin('caldyn_exner',def)
@@ -53,4 +65,7 @@
 
   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) 
   
   CALL allocate_field(f_buf_i,field_t,type_real,llm)
@@ -118,11 +133,8 @@
   TYPE(t_field),POINTER :: f_du(:)
 
-  REAL(rstd),POINTER :: phis(:)
-  REAL(rstd),POINTER :: ps(:)
-  REAL(rstd),POINTER :: theta_rhodz(:,:)
-  REAL(rstd),POINTER :: u(:,:)
-  REAL(rstd),POINTER :: dps(:)
-  REAL(rstd),POINTER :: dtheta_rhodz(:,:)
-  REAL(rstd),POINTER :: du(:,:)
+  REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
+  REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
+  REAL(rstd),POINTER :: u(:,:), du(:,:)
+  REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
   INTEGER :: ind,ij
 
@@ -137,14 +149,17 @@
       CALL swap_geometry(ind)
       
-      out_u=f_out_u(ind) 
       phis=f_phis(ind)
       ps=f_ps(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)
-      dps=f_dps(ind)
-      dtheta_rhodz=f_dtheta_rhodz(ind)
       du=f_du(ind)
+      out_u=f_out_u(ind) 
 !$OMP PARALLEL DEFAULT(SHARED)
-      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
+      CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
 !$OMP END PARALLEL
     ENDDO
@@ -161,16 +176,21 @@
   
   
-  SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
+  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
   USE icosa
   USE disvert_mod
   USE exner_mod
   IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    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)  :: phis(iim*jjm)
+    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(OUT) :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
+
     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)
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
     
     INTEGER :: i,j,ij,l
@@ -180,13 +200,11 @@
 
     REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
-    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
     REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
     REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
- ! Intermediate variable to compute exner function
+! Intermediate variable to compute exner function
     REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
     REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
- !   
+   
     REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
     REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
     REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
@@ -202,18 +220,16 @@
 !$OMP MASTER  
 !    IF (first) THEN
-      ALLOCATE(theta(iim*jjm,llm))  ! potential temperature
-      ALLOCATE(p(iim*jjm,llm+1))  ! pression
-      ALLOCATE(pk(iim*jjm,llm))   ! Exner function
+      ALLOCATE(theta(iim*jjm,llm))    ! potential temperature
+      ALLOCATE(pk(iim*jjm,llm))       ! Exner function
       ALLOCATE(pks(iim*jjm))
       ALLOCATE(alpha(iim*jjm,llm))
       ALLOCATE(beta(iim*jjm,llm))
-      ALLOCATE(phi(iim*jjm,llm))    ! geopotential
-      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
-      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
+      ALLOCATE(phi(iim*jjm,llm))      ! geopotential
+      ALLOCATE(Fe(3*iim*jjm,llm))     ! mass flux   
       ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
       ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
-      ALLOCATE(w(iim*jjm,llm))       ! vertical velocity      
-      ALLOCATE(qv(2*iim*jjm,llm))       ! potential velocity  
-      ALLOCATE(berni(iim*jjm,llm))   ! bernouilli term
+      ALLOCATE(w(iim*jjm,llm))        ! vertical velocity      
+      ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
+      ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
 !      first=.FALSE.
 !    ENDIF
@@ -292,6 +308,17 @@
         ENDDO
       ENDDO
-    ENDDO
-
+
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
+        qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
+        qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
+     END DO
+  END DO
+  
+ENDDO
+
+!----------- needed next : rhodz, p --------------!
 
   DO l = 1, llm
@@ -378,5 +405,5 @@
         ij=(j-1)*iim+i
 
-       du(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) *      &  
+       du(ij+u_right,l) = qu(ij+u_right,l)/de(ij+u_right) *      &  
                       ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
                         wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
@@ -391,5 +418,5 @@
 
 
-       du(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) *         &  
+       du(ij+u_lup,l) = qu(ij+u_lup,l)/de(ij+u_lup) *         &  
                       ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
                         wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
@@ -404,5 +431,5 @@
 
 
-       du(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) *   &  
+       du(ij+u_ldown,l) = qu(ij+u_ldown,l)/de(ij+u_ldown) *   &  
                       ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
                         wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
@@ -536,5 +563,4 @@
 !!$OMP MASTER  
     DEALLOCATE(theta)  ! potential temperature
-    DEALLOCATE(p)  ! pression
     DEALLOCATE(pk)   ! Exner function
     DEALLOCATE(pks)
@@ -542,6 +568,4 @@
     DEALLOCATE(beta)
     DEALLOCATE(phi)    ! geopotential
-!    DEALLOCATE(mass)   ! mass    
-    DEALLOCATE(rhodz)   ! mass density   
     DEALLOCATE(Fe)   ! mass flux   
     DEALLOCATE(Ftheta) ! theta flux   
