Index: /codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 125)
+++ /codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 126)
@@ -20,10 +20,10 @@
     CHARACTER(len=255) :: def
   
-    def='direct'
+    def='enstrophy'
     CALL getin('caldyn_conserv',def)
     SELECT CASE(TRIM(def))
     CASE('energy')
        caldyn_conserv=1
-    CASE('direct')
+    CASE('enstrophy')
        caldyn_conserv=2
     CASE DEFAULT
@@ -79,118 +79,81 @@
   END SUBROUTINE allocate_caldyn
    
-  SUBROUTINE check_mass_conservation(f_ps,f_dps)
-  USE icosa
-  IMPLICIT NONE
+  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
+    USE icosa
+    USE vorticity_mod
+    USE kinetic_mod
+    USE theta2theta_rhodz_mod
+    IMPLICIT NONE
+    INTEGER,INTENT(IN)    :: it
+    TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
     TYPE(t_field),POINTER :: f_dps(:)
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: dps(:)
-    REAL(rstd) :: mass_tot,dmass_tot
-    INTEGER :: ind,i,j,ij
-    
-    mass_tot=0
-    dmass_tot=0
-    
-    CALL transfert_request(f_dps,req_i1)
-    CALL transfert_request(f_ps,req_i1)
-
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-
-      ps=f_ps(ind)
-      dps=f_dps(ind)
-
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
-            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
-          ENDIF
-        ENDDO
-      ENDDO
-   
-    ENDDO
-    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        
-
-  END SUBROUTINE check_mass_conservation
-  
-  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q, f_dps, f_dtheta_rhodz, f_du)
-  USE icosa
-  USE vorticity_mod
-  USE kinetic_mod
-  USE theta2theta_rhodz_mod
-  IMPLICIT NONE
-  INTEGER,INTENT(IN)    :: it
-  TYPE(t_field),POINTER :: f_phis(:)
-  TYPE(t_field),POINTER :: f_ps(:)
-  TYPE(t_field),POINTER :: f_theta_rhodz(:)
-  TYPE(t_field),POINTER :: f_u(:)
-  TYPE(t_field),POINTER :: f_q(:)
-  TYPE(t_field),POINTER :: f_dps(:)
-  TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
-  TYPE(t_field),POINTER :: f_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
-
-  
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+    TYPE(t_field),POINTER :: f_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
+    
     CALL transfert_request(f_phis,req_i1) 
     CALL transfert_request(f_ps,req_i1) 
     CALL transfert_request(f_theta_rhodz,req_i1) 
     CALL transfert_request(f_u,req_e1)
-    
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(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)
-      du=f_du(ind)
-      out_u=f_out_u(ind) 
-!$OMP PARALLEL DEFAULT(SHARED)
-      CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
-!$OMP END PARALLEL
-    ENDDO
-
-    IF (mod(it,itau_out)==0 ) THEN
-       PRINT *,'CALL write_output_fields'
-       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)
-    END IF
-
-!    CALL check_mass_conservation(f_ps,f_dps)
-
-  END SUBROUTINE caldyn
-  
-  
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+
+    SELECT CASE(caldyn_conserv)
+    CASE(1) ! energy-conserving
+    CASE(2) ! enstrophy-conserving
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(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)
+          du=f_du(ind)
+          out_u=f_out_u(ind) 
+          !$OMP PARALLEL DEFAULT(SHARED)
+!          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          CALL compute_pvort(ps, u, p,rhodz,qu)
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, dps, dtheta_rhodz, du)
+          !$OMP END PARALLEL
+       ENDDO
+       
+       IF (mod(it,itau_out)==0 ) THEN
+          PRINT *,'CALL write_output_fields'
+          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)
+       END IF
+    
+    CASE DEFAULT
+       STOP
+    END SELECT
+
+    !    CALL check_mass_conservation(f_ps,f_dps)
+    
+END SUBROUTINE caldyn
+  
+  
+SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
   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(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)
     
     INTEGER :: i,j,ij,l
@@ -200,17 +163,12 @@
 
     REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
-    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
-    REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
-! 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 :: Fe(:,:)   ! mass flux   
-    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
-    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence
-    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity      
-    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity  
-    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term
+    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
+    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential
+    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux
+    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
+    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity      
+    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
+    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
     
     LOGICAL,SAVE :: first=.TRUE.
@@ -220,52 +178,29 @@
 !$OMP MASTER  
 !    IF (first) THEN
-      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(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
-!      first=.FALSE.
-!    ENDIF
-!$OMP END MASTER
-!$OMP BARRIER        
- !    du(:,:)=0
-!    theta=1e10
-!    p=1e10
-!    pk=1e10
-!    pks=1e10
-!    alpha=1e10
-!    beta=1e10
-!    phi=1e10
-!    mass=1e10
-!    rhodz=1e10
-!    Fe=1e10
-!    Ftheta=1e10
-!    convm=1e10
-!    w=1e10
-!    qv=1e10
-!    berni=1e10
-    
- !!! Compute pressure
-
-!    PRINT *, 'Computing pressure'
+    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(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
+    
+!!! Compute pressure
     DO    l    = 1, llm+1
-!$OMP DO
-      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 DO
+       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
+    
 !!! Compute mass
-   ! PRINT *, 'Computing mass, theta'
    DO l = 1, llm
 !$OMP DO
@@ -274,5 +209,4 @@
             ij=(j-1)*iim+i
             rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
-            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
          ENDDO
       ENDDO
@@ -309,16 +243,97 @@
       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
+      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
+
+!!$OMP BARRIER
+!!$OMP MASTER  
+    DEALLOCATE(theta)  ! potential temperature
+    DEALLOCATE(pk)   ! Exner function
+    DEALLOCATE(pks)
+    DEALLOCATE(alpha)
+    DEALLOCATE(beta)
+    DEALLOCATE(phi)    ! geopotential
+    DEALLOCATE(Fe)   ! mass flux   
+    DEALLOCATE(Ftheta) ! theta flux   
+    DEALLOCATE(convm)    ! mass flux convergence
+    DEALLOCATE(w)       ! vertical velocity      
+    DEALLOCATE(qv)       ! potential velocity  
+    DEALLOCATE(berni)   ! bernouilli term
+!!$OMP END MASTER
+!!$OMP BARRIER                                                      
+  END SUBROUTINE compute_pvort
+  
 
 !----------- needed next : rhodz, p --------------!
+
+  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)  :: 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)
+    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
+    
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: ww,uu 
+    REAL(rstd) :: delta
+    REAL(rstd) :: etav,hv, du2
+
+    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
+    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
+    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
+    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)              ! geopotential
+    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:), Ftheta(:,:)  ! mass flux, theta flux
+    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence
+    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)      ! vertical velocity      
+    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
+    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
+    
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+        
+!$OMP BARRIER      
+!$OMP MASTER  
+!    IF (first) THEN
+    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(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
+
+!!! Compute theta
+   DO l = 1, llm
+!$OMP DO
+      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 = 1, llm
@@ -577,4 +592,44 @@
 !!$OMP BARRIER                                                      
   END SUBROUTINE compute_caldyn
+
+!-------------------------------- Diagnostics ----------------------------
+
+  SUBROUTINE check_mass_conservation(f_ps,f_dps)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: dps(:)
+    REAL(rstd) :: mass_tot,dmass_tot
+    INTEGER :: ind,i,j,ij
+    
+    mass_tot=0
+    dmass_tot=0
+    
+    CALL transfert_request(f_dps,req_i1)
+    CALL transfert_request(f_ps,req_i1)
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+
+      ps=f_ps(ind)
+      dps=f_dps(ind)
+
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
+            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
+          ENDIF
+        ENDDO
+      ENDDO
+   
+    ENDDO
+    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        
+
+  END SUBROUTINE check_mass_conservation  
   
   SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
