Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 325)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 326)
@@ -7,8 +7,10 @@
     TYPE(t_field),POINTER,SAVE :: f_vort(:)
     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
-   
-  PUBLIC init_check_conserve, check_conserve 
+
+    INTEGER, PARAMETER :: check_basic=1, check_detailed=2
+    INTEGER, SAVE :: check_type
+  PUBLIC init_check_conserve, check_conserve, check_conserve_detailed
     REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
-!$OMP THREADPRIVATE(mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
+!$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
      
 CONTAINS 
@@ -17,6 +19,9 @@
 
   SUBROUTINE init_check_conserve
-  USE icosa
-  IMPLICIT NONE
+    USE icosa
+    USE getin_mod
+    USE omp_para, ONLY : is_master
+    IMPLICIT NONE
+    CHARACTER(LEN=255) :: check_type_str
     CALL allocate_field(f_pk,field_t,type_real,llm)
     CALL allocate_field(f_p,field_t,type_real,llm+1)
@@ -24,5 +29,20 @@
     CALL allocate_field(f_rhodz,field_t,type_real,llm)
     CALL allocate_field(f_vort,field_z,type_real,llm)
+    
+    check_type_str='basic'
+    CALL getin("check_conservation",check_type_str)
+    SELECT CASE(TRIM(check_type_str))
+    CASE('basic')
+       check_type=check_basic
+    CASE('detailed')
+       check_type=check_detailed
+    CASE DEFAULT
+       IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), &
+            '> for variable check_conservation : options are <basic>, <detailed>'
+       STOP
+    END SELECT
   END SUBROUTINE init_check_conserve
+
+!--------------------------------- Basic check --------------------------------
 
   SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
@@ -40,10 +60,10 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_phis(:)
-    INTEGER::it
+    INTEGER, INTENT(IN) :: it
 
     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
     INTEGER::ind,ierr
-    REAL(rstd) :: mtot, rmsdpdt
-    REAL(rstd) :: etot, stot, angtot, rmsvtot, ztot
+    REAL(rstd) :: mtot, angtot, rmsdpdt
+    REAL(rstd) :: etot, stot, ang_mass, ang_vit, rmsvtot, ztot
 
     CALL transfert_request(f_ue,req_e1_vect)
@@ -63,5 +83,6 @@
     CALL check_PV(ztot) 
     CALL exner(f_ps,f_p,f_pks,f_pk)
-    CALL check_EN(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 
+    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vit, rmsvtot) 
+    angtot  = ang_mass + ang_vit 
 
     IF (is_master) THEN 
@@ -70,5 +91,5 @@
           mtot0 = mtot
           etot0 = etot 
-          angtot0  = angtot 
+          angtot0  = angtot
           stot0 = stot 
        END IF
@@ -92,6 +113,81 @@
     END IF
   END SUBROUTINE check_conserve
-  
+
+!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
+
+  SUBROUTINE check_conserve_detailed(tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+    USE icosa 
+    USE pression_mod
+    USE vorticity_mod
+    USE caldyn_gcm_mod
+    USE exner_mod 
+    USE mpipara, ONLY : is_mpi_root, comm_icosa
+    IMPLICIT NONE
+    CHARACTER(LEN=*), INTENT(IN) :: tag
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    INTEGER, INTENT(IN) ::it
+
+    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
+    INTEGER::ind,ierr
+    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, AAM_mass, AAM_vel
+    
+    IF(check_type == check_detailed) THEN
+       etot=0.0; stot=0.0; rmsv=0.0 
+       AAM_vel=0. ; AAM_vel=0.
+       ztot = 0.0 
+
+       CALL transfert_request(f_ue,req_e1_vect)
+       CALL pression(f_ps,f_p)
+       DO ind=1,ndomain 
+          IF (.NOT. assigned_domain(ind)) CYCLE
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          p=f_p(ind) 
+          rhodz=f_rhodz(ind)
+          CALL compute_rhodz(p,rhodz) 
+       END DO
+       CALL exner(f_ps,f_p,f_pks,f_pk)
+       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, AAM_mass, AAM_vel, rmsv) 
+       
+       IF (is_mpi_root) THEN 
+          !$OMP MASTER       
+          IF ( it == itau0  ) Then 
+             ztot0 = ztot
+             mtot0 = mtot
+             etot0 = etot 
+             angtot0  = AAM_mass + AAM_vel
+             stot0 = stot 
+          END IF
+          
+          PRINT *, TRIM(tag), AAM_mass, AAM_vel
+          !$OMP END MASTER 
+       END IF
+    END IF
+  END SUBROUTINE check_conserve_detailed
+ 
 !---------------------------------------------------------------------
+
+  SUBROUTINE global_sum(loc_sum, glob_sum)
+    USE mpi_mod
+    USE mpipara
+    USE transfert_omp_mod
+    IMPLICIT NONE
+    REAL(rstd), INTENT(IN) :: loc_sum
+    REAL(rstd), INTENT(OUT) :: glob_sum
+    REAL(rstd) :: sum
+
+    CALL reduce_sum_omp(loc_sum, sum)
+    IF (using_mpi) THEN
+!$OMP MASTER    
+      CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+!$OMP END MASTER
+    ELSE
+      glob_sum=loc_sum
+    ENDIF 
+  END SUBROUTINE global_sum
 
   SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
@@ -128,33 +224,21 @@
       ENDDO
     ENDDO
-
-    IF (using_mpi) THEN
-        CALL reduce_sum_omp(mloc, mloc_mpi)
-        CALL reduce_sum_omp(rmsloc, rmsloc_mpi)
-!$OMP MASTER    
-      CALL MPI_REDUCE(mloc_mpi, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-      CALL MPI_REDUCE(rmsloc_mpi, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-!$OMP END MASTER
-    ELSE
-      mtot=mloc
-      rmsdpdt=rmsloc
-    ENDIF
-    
+    
+    CALL global_sum(mloc, mtot)
+    CALL global_sum(rmsloc, rmsdpdt)
   END SUBROUTINE check_mass_conserve
 
 !---------------------------------------------------------------------
 
-  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot)  
+  SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
+                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)  
   USE icosa
   USE pression_mod
   USE vorticity_mod
-  USE mpi_mod
-  USE mpipara
-  USE transfert_omp_mod
   IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ue(:)
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
     TYPE(t_field), POINTER :: f_phis(:)
-    REAL(rstd),INTENT(OUT) :: etot, stot, angtot, rmsvtot
+    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot
     REAL(rstd), POINTER :: ue(:,:)
     REAL(rstd), POINTER :: p(:,:) 
@@ -163,9 +247,8 @@
     REAL(rstd), POINTER :: phis(:) 
     REAL(rstd), POINTER :: rhodz(:,:) 
-    REAL(rstd) :: e, s, ang, rmsv
-    REAL(rstd) :: e_mpi, s_mpi, ang_mpi, rmsv_mpi
+    REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv
     INTEGER :: ind
 
-    e = 0 ; s = 0 ; ang = 0 ; rmsv = 0
+    e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
 
     DO ind=1,ndomain 
@@ -179,33 +262,21 @@
       pk=f_pk(ind) 
       phis=f_phis(ind) 
-      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
+      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz,pk,phis, &
+           e, s, AAM_mass, AAM_vel, rmsv)
     END DO 
 
-    IF (using_mpi) THEN
-      CALL reduce_sum_omp(e, e_mpi)
-      CALL reduce_sum_omp(s, s_mpi)
-      CALL reduce_sum_omp(ang, ang_mpi)
-      CALL reduce_sum_omp(rmsv, rmsv_mpi)
-!$OMP MASTER    
-      CALL MPI_REDUCE(e_mpi, etot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-      CALL MPI_REDUCE(s_mpi, stot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-      CALL MPI_REDUCE(ang_mpi, angtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-      CALL MPI_REDUCE(rmsv_mpi, rmsvtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
-!$OMP END MASTER
-    ELSE
-      etot=e
-      stot=s
-      angtot=ang
-      rmsvtot=rmsv
-    ENDIF 
-
-  END SUBROUTINE check_en 
-
-  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
-  USE icosa
-  USE disvert_mod
-  USE wind_mod
-  USE omp_para
-  IMPLICIT NONE
+    CALL global_sum(e, etot)
+    CALL global_sum(s, stot)
+    CALL global_sum(rmsv, rmsvtot)
+    CALL global_sum(AAM_mass, AAM_mass_tot)
+    CALL global_sum(AAM_vel, AAM_vel_tot)
+  END SUBROUTINE check_energy
+
+  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)
+    USE icosa
+    USE disvert_mod
+    USE wind_mod
+    USE omp_para
+    IMPLICIT NONE
     INTEGER,INTENT(IN)::ind 
     REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
@@ -215,81 +286,43 @@
     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd),INTENT(INOUT) :: e
-    REAL(rstd),INTENT(INOUT) :: s
-    REAL(rstd),INTENT(INOUT) :: ang
-    REAL(rstd),INTENT(INOUT) :: rmsv
-    
-    REAL(rstd):: theta(iim*jjm,llm)
-    REAL(rstd)::KE(iim*jjm,llm) 
+    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
+    
     REAL(rstd) :: ucenter(iim*jjm,3,llm)
     REAL(rstd) :: ulon(iim*jjm,llm)
     REAL(rstd) :: ulat(iim*jjm,llm)
-
-    REAL(rstd)::masse(iim*jjm,llm)  
-    REAL(rstd)::rad,radomeg,lat,lon
-    REAL(rstd) ::xyz(3)
-    INTEGER :: i,j,ij,l,ij2
-
-  
-    DO l = ll_begin, ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-	  ij=(j-1)*iim+i
-	  masse(ij,l) = rhodz(ij,l)*Ai(ij) 
-          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-          IF (domain(ind)%own(i,j)) THEN
-	    s = s + Ai(ij)*theta_rhodz(ij,l)
-          END IF 
-        END DO 
-      END DO
-    END DO 
-
-    DO l = ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-	    KE(ij,l)=   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 +           &
-                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
-                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
-                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
-                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
-          END IF 
-        END DO
-      END DO 
-    END DO  
-
-    DO l = ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-	  ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-	    e = e + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
-          END IF 
-        END DO 
-      END DO 
-    END DO 
-
-!------------------------------ ANGULAR VELOCITY 
+    
+    REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat
+    INTEGER :: i,j,ij,l
+    
     CALL compute_wind_centered(u,ucenter)
     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
     rad = radius 
     radomeg = rad*omega
-
+    
     DO l = ll_begin,ll_end
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-	  ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-            lat=lat_i(ij)  
-            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
-          END IF 
-        END DO
-      END DO
-    END DO    
-          
-  END SUBROUTINE compute_en
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             IF (domain(ind)%own(i,j)) THEN
+                lat = lat_i(ij)  
+                mass_ij = rhodz(ij,l)*Ai(ij) 
+                theta_ij = theta_rhodz(ij,l)/rhodz(ij,l)
+                KE_ij = 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 +           &
+                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
+                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
+                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
+                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+                rmsv = rmsv + 2*mass_ij*KE_ij
+                s    = s + Ai(ij)*theta_rhodz(ij,l)
+                e    = e + mass_ij*(phis(ij)+theta_ij*pk(ij,l)+KE_ij)
+                AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat)
+                AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l)
+             END IF
+          END DO
+       END DO
+    END DO
+    
+END SUBROUTINE compute_energy
 
 !---------------------------------------------------------------------
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 325)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 326)
@@ -211,5 +211,8 @@
   
   DO it=itau0+1,itau0+itaumax
-    
+     
+    CALL check_conserve_detailed('detailed_budget 0', &
+         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
+
     IF (xios_output) CALL xios_update_calendar(it)
     IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN
@@ -256,4 +259,7 @@
        END SELECT
     END DO
+
+    CALL check_conserve_detailed('detailed_budget 1', &
+         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
 
     IF (MOD(it,itau_dissip)==0) THEN
@@ -281,4 +287,7 @@
     END IF
 
+    CALL check_conserve_detailed('detailed_budget 2', &
+         f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)
+
     IF(MOD(it,itau_adv)==0) THEN
 
