Index: codes/icosagcm/devel/src/base/earth_const.f90
===================================================================
--- codes/icosagcm/devel/src/base/earth_const.f90	(revision 700)
+++ codes/icosagcm/devel/src/base/earth_const.f90	(revision 714)
@@ -13,5 +13,5 @@
   REAL(rstd),SAVE :: Treff=273.
   REAL(rstd),SAVE :: preff=101325.
-  REAL(rstd),SAVE :: pa=50000.
+  REAL(rstd),SAVE :: pa=50000. ! default value set to preff/2 by disvert_std
   REAL(rstd),SAVE :: scale_height=8000. ! atmospheric scale height (m)
   REAL(rstd),SAVE :: scale_factor=1.
Index: codes/icosagcm/devel/src/base/profiling.f90
===================================================================
--- codes/icosagcm/devel/src/base/profiling.f90	(revision 714)
+++ codes/icosagcm/devel/src/base/profiling.f90	(revision 714)
@@ -0,0 +1,90 @@
+MODULE profiling_mod
+  IMPLICIT NONE
+  SAVE
+  PRIVATE
+
+  INTEGER, PARAMETER :: max_id=20, max_depth=10
+  INTEGER :: nb_id=-10, depth=-10, current_id(max_depth)
+  CHARACTER(10), DIMENSION(max_id) :: name
+  REAL :: chrono(max_depth), elapsed(max_id)
+
+  PUBLIC :: init_profiling, register_id, enter_profile, exit_profile, print_profile
+
+CONTAINS
+
+  SUBROUTINE init_profiling
+    !$OMP MASTER
+    nb_id=0
+    depth=0
+    elapsed(:)=0
+    !$OMP END MASTER
+  END SUBROUTINE init_profiling
+
+  SUBROUTINE register_id(thename, id)
+    CHARACTER(*), INTENT(IN) :: thename
+    INTEGER, INTENT(OUT) :: id
+    !$OMP MASTER
+    nb_id = nb_id+1
+    id = nb_id
+    name(id)=thename
+    !$OMP END MASTER
+  END SUBROUTINE register_id
+
+  FUNCTION get_elapsed(start)
+    INTEGER(kind=8) :: count, count_rate
+    REAL :: start,get_elapsed
+    CALL SYSTEM_CLOCK(count,count_rate)
+    get_elapsed = (1.*count)/(1.*count_rate) - start
+    IF(get_elapsed<0.) get_elapsed=0.
+  END FUNCTION get_elapsed
+
+  SUBROUTINE enter_profile(id)
+    INTEGER, INTENT(IN) :: id
+    REAL :: my_chrono
+    !$OMP MASTER
+    depth = depth+1
+    chrono(depth) = get_elapsed(0.)
+    current_id(depth) = id
+    !$OMP END MASTER
+  END SUBROUTINE enter_profile
+   
+  SUBROUTINE exit_profile(id)
+    INTEGER, INTENT(IN) :: id
+    INTEGER :: parent_id
+    REAL :: my_elapsed
+    !$OMP MASTER
+    IF(depth<=0) THEN
+       PRINT *, 'exit_profile called at depth=0 !!'
+       STOP
+    END IF
+    IF(id /= current_id(depth)) THEN
+       PRINT *,' exit_profile : at depth ', depth, ' exiting ', TRIM(name(id)), ' after entering ', TRIM(name(current_id(depth)))
+       STOP
+    END IF
+    my_elapsed = get_elapsed(chrono(depth))
+    ! add elapsed to current profile
+    elapsed(id) = elapsed(id) + my_elapsed
+    depth = depth-1
+    ! and substract from parent profile
+    IF(depth>0) THEN
+       parent_id = current_id(depth)
+       elapsed(parent_id) = elapsed(parent_id) - my_elapsed
+    END IF
+    !$OMP END MASTER
+  END SUBROUTINE exit_profile
+
+  SUBROUTINE print_profile
+    INTEGER :: i
+    REAL :: total
+    !$OMP MASTER
+    PRINT *, '---------------------- Profiling -----------------------'
+    total = SUM(elapsed(1:nb_id))
+    PRINT *, 'Total (s) : ', total
+    DO i=1,nb_id
+       PRINT *, name(i), INT(elapsed(i)), INT(1000.*elapsed(i)/total)
+    END DO
+    PRINT *, '---------------------- Profiling -----------------------'
+    !$OMP END MASTER
+  END SUBROUTINE print_profile
+    
+END MODULE profiling_mod
Index: codes/icosagcm/devel/src/diagnostics/check_conserve.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 700)
+++ codes/icosagcm/devel/src/diagnostics/check_conserve.f90	(revision 714)
@@ -14,5 +14,7 @@
   
   REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, &
+       AAM_velp, AAM_velp_old, AAM_velm, AAM_velm_old, &
        AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master
+  REAL(rstd),SAVE :: AAM_vel_plus_source(3), AAM_vel_minus_source(3)
   REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
   !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
@@ -66,5 +68,5 @@
     INTEGER :: ind,ierr
     REAL(rstd) :: mtot, angtot, rmsdpdt
-    REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot
+    REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot
 
     CALL transfert_request(f_ue,req_e1_vect)
@@ -84,5 +86,5 @@
     CALL check_PV(ztot) 
     CALL exner(f_ps,f_p,f_pks,f_pk)
-    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot) 
+    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot) 
 
     IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master
@@ -90,4 +92,6 @@
        AAM_mass = ang_mass
        AAM_vel = ang_vel
+       AAM_velp = ang_velp
+       AAM_velm = ang_velm
        angtot  = ang_mass + ang_vel
        IF ( it == itau0  ) THEN 
@@ -99,6 +103,10 @@
           AAM_mass_old=AAM_mass
           AAM_vel_old=AAM_vel
+          AAM_velp_old=AAM_velp
+          AAM_velm_old=AAM_velm
           AAM_mass_source=0.
           AAM_vel_source=0.
+          AAM_vel_plus_source=0.
+          AAM_vel_minus_source=0.
        END IF
        rmsvtot=SQRT(rmsvtot/mtot) 
@@ -127,13 +135,17 @@
           WRITE(*,'(A,6E12.4)') 'AAM_vel  : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, &
                AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys)
-          WRITE(*,'(A,6E12.4)') 'AAM_tot  : time,old,new,dissip,dyn,phys', dt*it, &
-               AAM_vel_old+AAM_mass_old, AAM_vel+AAM_mass, &
-               AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), &
-               AAM_vel_source(AAM_dyn)   +AAM_mass_source(AAM_dyn), &
-               AAM_vel_source(AAM_phys)  +AAM_mass_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dt*it, AAM_velp_old, AAM_velp, &
+               AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dt*it, AAM_velm_old, AAM_velm, &
+               AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn)
           AAM_mass_old=AAM_mass
           AAM_vel_old=AAM_vel
+          AAM_velp_old=AAM_velp
+          AAM_velm_old=AAM_velm
           AAM_mass_source=0.
           AAM_vel_source=0.
+          AAM_vel_minus_source=0.
+          AAM_vel_plus_source=0.
        END IF
     END IF
@@ -160,5 +172,5 @@
     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
     INTEGER::ind,ierr
-    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel
+    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm
     
     IF(check_type == check_detailed) THEN
@@ -175,5 +187,5 @@
        END DO
        CALL exner(f_ps,f_p,f_pks,f_pk)
-       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv) 
+       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv) 
        
        IF (is_master) THEN
@@ -186,6 +198,15 @@
           AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass
           AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel
+          AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp
+          AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm
+          !if ((ang_vel - AAM_vel) .gt. 0.) then
+          !  AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel
+          !else
+          !  AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel
+          !endif
           AAM_mass = ang_mass
           AAM_vel = ang_vel
+          AAM_velp = ang_velp
+          AAM_velm = ang_velm
        END IF
     END IF
@@ -252,5 +273,5 @@
 
   SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
-                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)  
+                      stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot)  
   USE pression_mod
   USE vorticity_mod
@@ -258,5 +279,5 @@
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
     TYPE(t_field), POINTER :: f_phis(:)
-    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot
+    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot
     REAL(rstd), POINTER :: ue(:,:)
     REAL(rstd), POINTER :: p(:,:) 
@@ -265,8 +286,9 @@
     REAL(rstd), POINTER :: phis(:) 
     REAL(rstd), POINTER :: rhodz(:,:) 
-    REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv
+    REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
     INTEGER :: ind
 
     e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
+    AAM_velp = 0 ; AAM_velm = 0
 
     DO ind=1,ndomain 
@@ -281,5 +303,5 @@
       phis=f_phis(ind) 
       CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
-           e, s, AAM_mass, AAM_vel, rmsv)
+           e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
     END DO 
 
@@ -289,7 +311,9 @@
     CALL global_sum(AAM_mass, AAM_mass_tot)
     CALL global_sum(AAM_vel, AAM_vel_tot)
+    CALL global_sum(AAM_velp, AAM_velp_tot)
+    CALL global_sum(AAM_velm, AAM_velm_tot)
   END SUBROUTINE check_energy
 
-  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)
+  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
     USE disvert_mod
     USE wind_mod
@@ -302,5 +326,5 @@
     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
+    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
     
     REAL(rstd) :: ucenter(iim*jjm,llm,3)
@@ -335,4 +359,9 @@
                 AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat)
                 AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                if (ulon(ij,l) > 0.) then
+                  AAM_velp  = AAM_velp  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                else
+                  AAM_velm  = AAM_velm  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                endif
              END IF
           END DO
Index: codes/icosagcm/devel/src/diagnostics/observable.f90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 700)
+++ codes/icosagcm/devel/src/diagnostics/observable.f90	(revision 714)
@@ -15,5 +15,6 @@
   TYPE(t_field),POINTER, SAVE :: f_theta(:)
 
-  PUBLIC init_observable, write_output_fields_basic, f_theta, f_buf_i
+  PUBLIC init_observable, write_output_fields_basic, & 
+       f_theta, f_buf_i, f_buf_ulon, f_buf_ulat
 
 CONTAINS
Index: codes/icosagcm/devel/src/diagnostics/wind.F90
===================================================================
--- codes/icosagcm/devel/src/diagnostics/wind.F90	(revision 700)
+++ codes/icosagcm/devel/src/diagnostics/wind.F90	(revision 714)
@@ -12,10 +12,12 @@
 CONTAINS
 
-  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
+  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat, scale_)
     TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
     TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons    
-    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    REAL(rstd),POINTER     :: u(:,:),  ulon(:,:), ulat(:,:)
+    REAL(rstd), OPTIONAL   :: scale_ 
+    REAL(rstd)             :: scale
     INTEGER :: ind
-
+    scale = MERGE(scale_, 1., PRESENT(scale_))
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -25,5 +27,5 @@
        ulon=f_ulon(ind)
        ulat=f_ulat(ind)
-       CALL compute_un2ulonlat(u,ulon, ulat)
+       CALL compute_un2ulonlat(u,ulon, ulat, scale)
     END DO
 
@@ -49,10 +51,11 @@
   END SUBROUTINE ulonlat2un
  
-  SUBROUTINE compute_wind_centered(ue,ucenter)
+  SUBROUTINE compute_wind_centered(ue,ucenter,scale_)
   REAL(rstd) :: ue(3*iim*jjm,llm)
   REAL(rstd) :: ucenter(iim*jjm,llm,3)
+  REAL(rstd), INTENT(IN), OPTIONAL :: scale_
   INTEGER :: ij,l
-  REAL(rstd), PARAMETER :: scale=1.
-  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz
+  REAL(rstd) :: scale,fac, ue_le, cx,cy,cz, ux,uy,uz
+  scale = MERGE(scale_, 1., PRESENT(scale_))
 #include "../kernels_hex/wind_centered.k90"
  END SUBROUTINE compute_wind_centered
@@ -328,12 +331,11 @@
 
 
- SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
+ SUBROUTINE compute_un2ulonlat(un, ulon, ulat, scale)
   REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
   REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
-
   REAL(rstd)             :: uc(iim*jjm,llm,3)
-    
-  CALL compute_wind_centered(un,uc)  
+  REAL(rstd)             :: scale
+  CALL compute_wind_centered(un,uc,scale)
   CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
  
Index: codes/icosagcm/devel/src/dissip/dissip_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/dissip/dissip_gcm.f90	(revision 700)
+++ codes/icosagcm/devel/src/dissip/dissip_gcm.f90	(revision 714)
@@ -35,6 +35,7 @@
 !$OMP THREADPRIVATE(rayleigh_friction_type)
 !$OMP THREADPRIVATE(rayleigh_shear)
-  REAL, SAVE    :: rayleigh_tau
+  REAL, SAVE    :: rayleigh_tau, rayleigh_limlat
 !$OMP THREADPRIVATE(rayleigh_tau)
+!$OMP THREADPRIVATE(rayleigh_limlat)
   
   REAL,SAVE    :: dtdissip
@@ -103,4 +104,7 @@
       rayleigh_friction_type=2
       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+   CASE('giant_liu_schneider') 
+      rayleigh_friction_type=99
+      IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010'
    CASE DEFAULT
       IF (is_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), &
@@ -117,4 +121,9 @@
          STOP
       END IF
+      IF(rayleigh_friction_type == 99) THEN
+         rayleigh_limlat=0.
+         CALL getin("rayleigh_limlat",rayleigh_limlat)
+         rayleigh_limlat = rayleigh_limlat*3.14159/180.
+      ENDIF
    END IF
 
@@ -444,13 +453,18 @@
 
     IF(mintau>0) THEN
-       itau_dissip=INT(mintau/dt)
-       dtdissip=itau_dissip*dt
+       IF (itau_dissip==0) THEN
+         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
+         itau_dissip=INT(mintau/dt)
+       ENDIF
     ELSE
-       IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
+       IF (is_master) PRINT *,"init_dissip: No dissipation time set, setting itau_dissip to 1000000000"
        itau_dissip=100000000
     END IF
     itau_dissip=MAX(1,itau_dissip)
-    IF (is_master) PRINT *,"rayleigh_tau",rayleigh_tau, "mintau ",mintau, &
-         "itau_dissip",itau_dissip," dtdissip ",dtdissip
+    dtdissip=itau_dissip*dt
+    IF (is_master) THEN
+      PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
+      PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
+    ENDIF
 
   END SUBROUTINE init_dissip 
@@ -479,5 +493,5 @@
 
     INTEGER :: ind, shear
-    INTEGER :: l,ij
+    INTEGER :: l,ij,nn
 
 !$OMP BARRIER
@@ -515,4 +529,5 @@
 
       IF(rayleigh_friction_type>0) THEN
+       IF(rayleigh_friction_type<99) THEN
          phi=f_geopot(ind)
          ue=f_ue(ind)
@@ -524,4 +539,22 @@
             ENDDO
          END DO
+       ELSE
+         ue=f_ue(ind)
+            DO ij=ij_begin,ij_end
+              nn = ij+u_right
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                !print*, "latitude", lat_e(nn)*180./3.14159
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+              nn = ij+u_lup
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+              nn = ij+u_ldown
+              IF (ABS(lat_e(nn)) .gt. rayleigh_limlat) THEN
+                due(nn,ll_begin:ll_begin+1) = due(nn,ll_begin:ll_begin+1) - (ue(nn,ll_begin:ll_begin+1)/rayleigh_tau)
+              ENDIF
+            ENDDO
+       ENDIF
       END IF
    END DO
@@ -529,7 +562,9 @@
    CALL trace_end("dissip")
 
+   CALL write_dissip_tendencies
 !$OMP BARRIER
 
     CONTAINS
+
       SUBROUTINE relax(shift_t, shift_u)
         USE dcmip_initial_conditions_test_1_2_3
@@ -559,4 +594,20 @@
        END SUBROUTINE relax
       
+       SUBROUTINE write_dissip_tendencies
+         USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
+         USE wind_mod
+         USE output_field_mod
+
+         CALL transfert_request(f_due_diss1,req_e1_vect)
+         CALL un2ulonlat(f_due_diss1, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
+         CALL output_field("dulon_diss1",f_buf_ulon)
+         CALL output_field("dulat_diss1",f_buf_ulat)
+!
+         CALL transfert_request(f_due_diss2,req_e1_vect)
+         CALL un2ulonlat(f_due_diss2, f_buf_ulon, f_buf_ulat, (1./(tau_graddiv(1))))
+         CALL output_field("dulon_diss2",f_buf_ulon)
+         CALL output_field("dulat_diss2",f_buf_ulat)
+       END SUBROUTINE write_dissip_tendencies
+
   END SUBROUTINE dissip
 
Index: codes/icosagcm/devel/src/icosa_init.f90
===================================================================
--- codes/icosagcm/devel/src/icosa_init.f90	(revision 700)
+++ codes/icosagcm/devel/src/icosa_init.f90	(revision 714)
@@ -21,6 +21,8 @@
   USE etat0_mod
   USE diagflux_mod
+  USE profiling_mod
   IMPLICIT NONE
   
+    CALL init_profiling
     CALL init_mpipara
     CALL trace_off
@@ -56,4 +58,5 @@
 
     CALL init_diagflux
+    CALL zero_du_phys
     CALL timeloop
     CALL switch_omp_no_distrib_level
Index: codes/icosagcm/devel/src/output/output_field.f90
===================================================================
--- codes/icosagcm/devel/src/output/output_field.f90	(revision 700)
+++ codes/icosagcm/devel/src/output/output_field.f90	(revision 714)
@@ -1,9 +1,16 @@
 MODULE output_field_mod
-USE genmod
+  USE genmod
+  USE xios_mod
+  USE profiling_mod
+  IMPLICIT NONE
+  SAVE
   PRIVATE
+
   LOGICAL,SAVE :: xios_output 
 !$OMP THREADPRIVATE(xios_output)     
   LOGICAL,SAVE :: enable_io    
 !$OMP THREADPRIVATE(enable_io)     
+
+  INTEGER :: id_output
 
   PUBLIC enable_io,xios_output,output_field_init,output_field,output_field_finalize
@@ -13,7 +20,7 @@
   SUBROUTINE output_field_init
   USE getin_mod
-  USE xios_mod
-  USE write_field_mod
   IMPLICIT NONE
+
+    CALL register_id('output',id_output)
 
     enable_io=.TRUE.
@@ -30,31 +37,29 @@
       CALL xios_init_write_field
     ENDIF
-
-    
   END SUBROUTINE output_field_init
 
   SUBROUTINE output_field(name_in,field)
     USE field_mod
-    USE xios_mod
     USE write_field_mod
     IMPLICIT NONE  
-     CHARACTER(LEN=*),INTENT(IN) :: name_in
-      TYPE(t_field),POINTER :: field(:)
-
-      IF (xios_output) THEN
-        CALL xios_write_field(name_in,field)
-      ELSE
+    CHARACTER(LEN=*),INTENT(IN) :: name_in
+    TYPE(t_field),POINTER :: field(:)
+    
+    CALL enter_profile(id_output)
+    IF (xios_output) THEN
+       CALL xios_write_field(name_in,field)
+    ELSE
        CALL writeField(name_in,field)
-      ENDIF
+    ENDIF
+    CALL exit_profile(id_output)
 
   END SUBROUTINE output_field
 
   SUBROUTINE output_field_finalize
-  USE ioipsl
-  USE xios_mod
-  IMPLICIT NONE
+    USE ioipsl
+    IMPLICIT NONE
     
     IF (xios_output) THEN
-      CALL xios_write_field_finalize
+       CALL xios_write_field_finalize
     ENDIF
     
Index: codes/icosagcm/devel/src/parallel/mpipara.F90
===================================================================
--- codes/icosagcm/devel/src/parallel/mpipara.F90	(revision 700)
+++ codes/icosagcm/devel/src/parallel/mpipara.F90	(revision 714)
@@ -12,5 +12,6 @@
   INTEGER,SAVE :: mpi_master
   
-  
+  INTEGER,SAVE :: id_mpi ! id for profiling
+
   INTERFACE allocate_mpi_buffer
     MODULE PROCEDURE allocate_mpi_buffer_r2, allocate_mpi_buffer_r3,allocate_mpi_buffer_r4
Index: codes/icosagcm/devel/src/parallel/transfert_mpi.f90
===================================================================
--- codes/icosagcm/devel/src/parallel/transfert_mpi.f90	(revision 700)
+++ codes/icosagcm/devel/src/parallel/transfert_mpi.f90	(revision 714)
@@ -2,4 +2,5 @@
 USE genmod
 USE field_mod
+IMPLICIT NONE
   
   TYPE array
@@ -85,10 +86,9 @@
   END INTERFACE
 
-
-  
 CONTAINS
        
       
   SUBROUTINE init_transfert
+  USE profiling_mod
   USE domain_mod
   USE dimensions
@@ -100,5 +100,7 @@
   INTEGER :: ind,i,j
   LOGICAL ::ok
- 
+
+    CALL register_id('MPI', id_mpi)
+
     CALL create_request(field_t,req_i1)
 
@@ -1094,4 +1096,5 @@
 
   SUBROUTINE send_message_mpi(field,message)
+  USE profiling_mod
   USE field_mod
   USE domain_mod
@@ -1124,4 +1127,6 @@
 
 !    CALL trace_start("send_message_mpi")
+
+    CALL enter_profile(id_mpi)
 
 !$OMP BARRIER
@@ -1482,4 +1487,6 @@
 !$OMP BARRIER
 !    CALL trace_end("send_message_mpi")
+
+    CALL exit_profile(id_mpi)
     
   END SUBROUTINE send_message_mpi
@@ -1499,4 +1506,5 @@
    
   SUBROUTINE wait_message_mpi(message)
+  USE profiling_mod
   USE field_mod
   USE domain_mod
@@ -1527,4 +1535,6 @@
     message%open=.FALSE.
     IF (.NOT. message%pending) RETURN
+
+    CALL enter_profile(id_mpi)
 
 !    CALL trace_start("wait_message_mpi")
@@ -1667,4 +1677,6 @@
 !$OMP BARRIER
     
+    CALL exit_profile(id_mpi)
+
   END SUBROUTINE wait_message_mpi
 
Index: codes/icosagcm/devel/src/physics/physics.f90
===================================================================
--- codes/icosagcm/devel/src/physics/physics.f90	(revision 700)
+++ codes/icosagcm/devel/src/physics/physics.f90	(revision 714)
@@ -1,6 +1,8 @@
 MODULE physics_mod
-
+  USE icosa
   USE field_mod
-
+  USE physics_interface_mod
+  USE omp_para
+  IMPLICIT NONE
   PRIVATE
 
@@ -14,9 +16,10 @@
   TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
   TYPE(t_field),POINTER,SAVE :: f_temp(:)
+  TYPE(t_field),POINTER,SAVE :: f_du_phys(:)
 
   CHARACTER(LEN=255),SAVE :: physics_type
 !$OMP THREADPRIVATE(physics_type)
 
-  PUBLIC :: physics, init_physics
+  PUBLIC :: physics, init_physics, zero_du_phys
 
 CONTAINS
@@ -25,6 +28,4 @@
     USE mpipara
     USE etat0_mod
-    USE icosa
-    USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
     USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
@@ -32,5 +33,4 @@
     USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
     USE physics_external_mod, ONLY : init_physics_external=>init_physics
-    IMPLICIT NONE
 
     physics_inout%dt_phys = dt*itau_physics
@@ -84,10 +84,37 @@
     END SELECT
 
+    CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys')
+
     IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
   END SUBROUTINE init_physics
 
+  SUBROUTINE zero_du_phys()
+    REAL(rstd), DIMENSION(:,:), POINTER :: du
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       du=f_du_phys(ind)
+       du(:,ll_begin:ll_end) = 0.
+    END DO
+  END SUBROUTINE zero_du_phys
+
+  SUBROUTINE add_du_phys(coef, f_u)
+    REAL(rstd), INTENT(IN) :: coef  ! -1 before physics, +1 after physics
+    TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics
+    REAL(rstd), DIMENSION(:,:), POINTER :: u, du
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       du=f_du_phys(ind)
+       u=f_u(ind)
+       du(:,ll_begin:ll_end) = du(:,ll_begin:ll_end) + coef*u(:,ll_begin:ll_end)
+    END DO
+  END SUBROUTINE add_du_phys
+
   SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
-    USE icosa
-    USE physics_interface_mod
     USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
     USE physics_external_mod, ONLY : physics_external => physics
@@ -96,5 +123,4 @@
     USE etat0_heldsz_mod
     USE etat0_venus_mod, ONLY : phys_venus => physics
-    IMPLICIT NONE
     INTEGER, INTENT(IN)   :: it
     TYPE(t_field),POINTER :: f_phis(:)
@@ -109,9 +135,11 @@
     TYPE(t_physics_inout) :: args
 
-    IF(MOD(it,itau_physics)==0) THEN
-    
+    IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN
+
+       ! as a result of the the two calls to add_du_phys,
+       ! du_phys increases by u(after physics) - u (before physics)
+       CALL add_du_phys(-1., f_ue)
+
        SELECT CASE(phys_type)
-       CASE (phys_none)
-          ! No physics, do nothing
        CASE(phys_HS94)
           CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
@@ -129,7 +157,11 @@
        CALL transfert_request(f_ue,req_e0_vect)
        CALL transfert_request(f_q,req_i0)
+
+       CALL add_du_phys(1., f_ue)
     END IF
 
     IF (mod(it,itau_out)==0 ) THEN
+       CALL write_physics_tendencies
+       CALL zero_du_phys
        SELECT CASE(phys_type)
        CASE (phys_DCMIP)
@@ -142,14 +174,20 @@
   END SUBROUTINE physics
 
+  SUBROUTINE write_physics_tendencies
+    USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
+    USE wind_mod
+    USE output_field_mod
+    CALL transfert_request(f_du_phys,req_e1_vect)
+    CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out)))
+    CALL output_field("dulon_phys",f_buf_ulon)
+    CALL output_field("dulat_phys",f_buf_ulat)
+  END SUBROUTINE write_physics_tendencies
+    
   SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-    USE icosa
-    USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
     USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
     USE theta2theta_rhodz_mod
     USE mpipara
-    USE omp_para
     USE checksum_mod
-    IMPLICIT NONE
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_ps(:)
@@ -229,12 +267,8 @@
 
   SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
-    USE icosa
     USE wind_mod
     USE pression_mod
     USE theta2theta_rhodz_mod
-    USE physics_interface_mod
     USE exner_mod
-    USE omp_para
-    IMPLICIT NONE
     TYPE(t_pack_info) :: info
     REAL(rstd) :: phis(iim*jjm)
@@ -272,9 +306,5 @@
 
   SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
-    USE icosa
-    USE physics_interface_mod
     USE theta2theta_rhodz_mod
-    USE omp_para
-    IMPLICIT NONE
     TYPE(t_pack_info) :: info
     REAL(rstd) :: ps(iim*jjm)
@@ -303,9 +333,5 @@
 
   SUBROUTINE compute_update_velocity(dulon, dulat, ue)
-    USE icosa
-    USE physics_interface_mod
     USE wind_mod
-    USE omp_para
-    IMPLICIT NONE
     REAL(rstd) :: dulon(iim*jjm,llm)
     REAL(rstd) :: dulat(iim*jjm,llm)
Index: codes/icosagcm/devel/src/time/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 700)
+++ codes/icosagcm/devel/src/time/timeloop_gcm.f90	(revision 714)
@@ -1,3 +1,4 @@
 MODULE timeloop_gcm_mod
+  USE profiling_mod
   USE icosa
   USE disvert_mod
@@ -13,5 +14,5 @@
   TYPE(t_message),SAVE :: req_ps0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
   LOGICAL, SAVE :: positive_theta
-
+  INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
   PUBLIC :: init_timeloop, timeloop
 
@@ -32,8 +33,18 @@
     CHARACTER(len=255) :: def
 
+    CALL register_id('timeloop',id_timeloop)
+    CALL register_id('dyn',id_dyn)
+    CALL register_id('dissip',id_dissip)
+    CALL register_id('phys',id_phys)
+    CALL register_id('adv',id_adv)
+    CALL register_id('diags',id_diags)
+
     CALL init_caldyn
     
 !    IF (xios_output) itau_out=1
     IF (.NOT. enable_io) itau_out=HUGE(itau_out)
+
+    itau_prof=100
+    CALL getin('itau_profiling',itau_prof)
 
     positive_theta=.FALSE.
@@ -231,11 +242,11 @@
 
     !$OMP MASTER
-    CALL SYSTEM_CLOCK(start_clock)
-    CALL SYSTEM_CLOCK(count_rate=rate_clock)
+    CALL SYSTEM_CLOCK(start_clock, rate_clock)
     !$OMP END MASTER   
 
     DO it=itau0+1,itau0+itaumax
-
        IF (is_master) CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)
+
+       CALL enter_profile(id_timeloop)
 
        IF (xios_output) THEN
@@ -267,4 +278,5 @@
        CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
 
+       CALL enter_profile(id_dyn)
        SELECT CASE(scheme_family)
        CASE(explicit)
@@ -273,5 +285,7 @@
           CALL HEVI_scheme(it, fluxt_zero)
        END SELECT
-
+       CALL exit_profile(id_dyn)
+
+       CALL enter_profile(id_dissip)
        IF (MOD(it,itau_dissip)==0) THEN
           
@@ -288,6 +302,8 @@
           ENDIF
           
+          CALL enter_profile(id_diags)
           CALL check_conserve_detailed(it, AAM_dyn, &
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
        
           CALL dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_u, f_dtheta_rhodz,f_du)
@@ -299,8 +315,12 @@
           END IF
 
+          CALL enter_profile(id_diags)
           CALL check_conserve_detailed(it, AAM_dissip, &
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL exit_profile(id_diags)
        END IF
+       CALL exit_profile(id_dissip)
        
+       CALL enter_profile(id_adv)
        IF(MOD(it,itau_adv)==0) THEN
           CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step
@@ -322,9 +342,13 @@
           IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
        END IF
+       CALL exit_profile(id_adv)
        
+       CALL enter_profile(id_diags)
        IF (MOD(it,itau_physics)==0) THEN
           CALL check_conserve_detailed(it, AAM_dyn, &
             f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
+          CALL enter_profile(id_phys)
           CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)        
+          CALL exit_profile(id_phys)
           CALL check_conserve_detailed(it, AAM_phys, &
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
@@ -345,5 +369,7 @@
           CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
        ENDIF
-
+       CALL exit_profile(id_diags)
+
+       CALL exit_profile(id_timeloop)
     END DO
     
@@ -382,4 +408,6 @@
             '  -- Completion in (min) : ', INT((total-elapsed)/60.)
     END IF
+    IF(MOD(it,itau_prof)==0) CALL print_profile
+
   END SUBROUTINE print_iteration
 
