Index: codes/icosagcm/devel/src/dissip/dissip_gcm.f90
===================================================================
--- codes/icosagcm/devel/src/dissip/dissip_gcm.f90	(revision 844)
+++ codes/icosagcm/devel/src/dissip/dissip_gcm.f90	(revision 845)
@@ -1,5 +1,7 @@
 MODULE dissip_gcm_mod
   USE icosa
-
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
   PRIVATE
 
@@ -41,11 +43,10 @@
   REAL,SAVE    :: dtdissip
 !$OMP THREADPRIVATE(dtdissip)
- 
+  
   PUBLIC init_dissip, dissip
+
 CONTAINS
-
+  
   SUBROUTINE allocate_dissip
-  USE icosa
-  IMPLICIT NONE  
     CALL allocate_field(f_due_diss1,field_u,type_real,llm)
     CALL allocate_field(f_due_diss2,field_u,type_real,llm)
@@ -58,376 +59,180 @@
  
   SUBROUTINE init_dissip
-  USE icosa
-  USE disvert_mod
-  USE mpi_mod
-  USE mpipara
-  USE transfert_mod
-  USE time_mod
-  USE transfert_omp_mod
-  USE omp_para
-  IMPLICIT NONE
-  
-  TYPE(t_field),POINTER,SAVE  :: f_u(:)
-  TYPE(t_field),POINTER,SAVE  :: f_du(:)
-  REAL(rstd),POINTER     :: u(:)
-  REAL(rstd),POINTER     :: du(:)
-  TYPE(t_field),POINTER,SAVE  :: f_theta(:)
-  TYPE(t_field),POINTER ,SAVE :: f_dtheta(:)
-  REAL(rstd),POINTER    :: theta(:)
-  REAL(rstd),POINTER    :: dtheta(:)
-  REAL(rstd)            :: dumax,dumax1
-  REAL(rstd)            :: dthetamax,dthetamax1
-  REAL(rstd)            :: r
-  REAL(rstd)            :: tau
-  REAL(rstd)            :: zz, zvert, fact
-  INTEGER               :: l
-  CHARACTER(len=255)    :: rayleigh_friction_key
-  REAL(rstd)            :: mintau
-  INTEGER               :: seed_size
-  INTEGER,ALLOCATABLE   :: seed(:)
-  
-            
-  INTEGER :: i,j,ij,ind,it,iter,M
-
-   rayleigh_friction_key='none'
-   CALL getin("rayleigh_friction_type",rayleigh_friction_key)
-   SELECT CASE(TRIM(rayleigh_friction_key))
-   CASE('none')
-      rayleigh_friction_type=0
-      IF (is_master) PRINT *, 'No Rayleigh friction'
-   CASE('dcmip2_schaer_noshear')
-      rayleigh_friction_type=1
-      rayleigh_shear=0
-      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
-   CASE('dcmip2_schaer_shear')
-      rayleigh_shear=1
-      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), &
-        ' in dissip_gcm.f90/init_dissip'
-      STOP
-   END SELECT
-
-   IF(rayleigh_friction_type>0) THEN
-      rayleigh_tau=0.
-      CALL getin("rayleigh_friction_tau",rayleigh_tau)
-      rayleigh_tau = rayleigh_tau / scale_factor
-      IF(rayleigh_tau<=0) THEN
-         IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
-         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
-
+    REAL(rstd)            :: tau  
+    CHARACTER(len=255)    :: rayleigh_friction_key
+    rayleigh_friction_key='none'
+    CALL getin("rayleigh_friction_type",rayleigh_friction_key)
+    SELECT CASE(TRIM(rayleigh_friction_key))
+    CASE('none')
+       rayleigh_friction_type=0
+       IF (is_master) PRINT *, 'No Rayleigh friction'
+    CASE('dcmip2_schaer_noshear')
+       rayleigh_friction_type=1
+       rayleigh_shear=0
+       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
+    CASE('dcmip2_schaer_shear')
+       rayleigh_shear=1
+       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), &
+            ' in dissip_gcm.f90/init_dissip'
+       STOP
+    END SELECT
+    
+    IF(rayleigh_friction_type>0) THEN
+       rayleigh_tau=0.
+       CALL getin("rayleigh_friction_tau",rayleigh_tau)
+       rayleigh_tau = rayleigh_tau / scale_factor
+       IF(rayleigh_tau<=0) THEN
+          IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
+          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
+    
     CALL allocate_dissip
+
+    CALL init_message(f_due_diss1,req_e1_vect,req_due)
+    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
+
+    tau_graddiv(:)=5000
+    CALL getin("tau_graddiv",tau)
+    tau_graddiv(:)=tau/scale_factor
+
+    CALL getin("nitergdiv",nitergdiv)
+    
+    tau_gradrot(:)=5000
+    CALL getin("tau_gradrot",tau)
+    tau_gradrot(:)=tau/scale_factor
+
+    CALL getin("nitergrot",nitergrot)
+    
+
+    tau_divgrad(:)=5000
+    CALL getin("tau_divgrad",tau)
+    tau_divgrad(:)=tau/scale_factor
+
+    CALL getin("niterdivgrad",niterdivgrad)
+
+    CALL dissip_constants
+    CALL dissip_profile
+  END SUBROUTINE init_dissip
+
+  SUBROUTINE dissip_constants
+    ! The SAVE attribute is required here, otherwise allocate_field will not work with OpenMP
+    TYPE(t_field),POINTER, SAVE  :: f_u(:)
+    TYPE(t_field),POINTER, SAVE  :: f_du(:)
+    TYPE(t_field),POINTER, SAVE  :: f_theta(:)
+    TYPE(t_field),POINTER, SAVE  :: f_dtheta(:)
+    REAL(rstd),POINTER :: u(:)
+    REAL(rstd),POINTER :: du(:)
+    REAL(rstd),POINTER :: theta(:)
+    REAL(rstd),POINTER :: dtheta(:)
+    REAL(rstd) :: dumax, dthetamax
+    INTEGER :: it, iter, ind
     CALL allocate_field(f_u,field_u,type_real)
     CALL allocate_field(f_du,field_u,type_real)
     CALL allocate_field(f_theta,field_t,type_real)
     CALL allocate_field(f_dtheta,field_t,type_real)
-    
-    CALL init_message(f_due_diss1,req_e1_vect,req_due)
-    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
-
-    tau_graddiv(:)=5000
-    CALL getin("tau_graddiv",tau)
-    tau_graddiv(:)=tau/scale_factor
-
-    CALL getin("nitergdiv",nitergdiv)
-    
-    tau_gradrot(:)=5000
-    CALL getin("tau_gradrot",tau)
-    tau_gradrot(:)=tau/scale_factor
-
-    CALL getin("nitergrot",nitergrot)
-    
-
-    tau_divgrad(:)=5000
-    CALL getin("tau_divgrad",tau)
-    tau_divgrad(:)=tau/scale_factor
-
-    CALL getin("niterdivgrad",niterdivgrad)
-
-
-    cgraddiv=-1
-    cdivgrad=-1
-    cgradrot=-1
-
-!$OMP BARRIER
-!$OMP MASTER
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      u=f_u(ind)
-
-! set random seed to get reproductibility when using a different number of process
-      CALL RANDOM_SEED(size=M)
-      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
-
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i   
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_right)=r-0.5
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_lup)=r-0.5
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_ldown)=r-0.5
-        ENDDO
-      ENDDO        
-    ENDDO
-!$OMP END MASTER  
-!$OMP BARRIER
-
+
+    PRINT *, '=========', SHAPE(f_u)
+
+    !--------------------------- Compute cgraddiv ----------------------------
+    cgraddiv=-1.
+    CALL random_vel(f_u)    
     DO it=1,20
-      
-      dumax=0
-      DO iter=1,nitergdiv
-        CALL transfert_request(f_u,req_e1_vect)
-        DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          CALL compute_gradiv(u,du,1,1)
-          u=du
-        ENDDO
-      ENDDO
-
-      CALL transfert_request(f_du,req_e1_vect)
-      
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        u=f_u(ind)
-        du=f_du(ind)
-         
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i   
-            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
-            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup)))
-            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
+       DO iter=1,nitergdiv
+          CALL transfert_request(f_u,req_e1_vect)
+          DO ind=1,ndomain
+             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+             CALL swap_dimensions(ind)
+             CALL swap_geometry(ind)
+             u=f_u(ind)
+             du=f_du(ind)
+             CALL compute_gradiv(u,du,1,1)
           ENDDO
-        ENDDO
-      ENDDO
-
-      IF (using_mpi) THEN 
-        CALL reduce_max_omp(dumax,dumax1)
-!$OMP MASTER        
-        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
-!$OMP END MASTER     
-        CALL bcast_omp(dumax)  
-      ELSE
-        CALL allreduce_max_omp(dumax,dumax1)
-        dumax=dumax1
-      ENDIF  
-			
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        u=f_u(ind)
-        du=f_du(ind)
-        u=du/dumax
-      ENDDO
-      IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
-
-    ENDDO 
-    IF (is_master) PRINT *,"gradiv : dumax",dumax
-    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
-                              'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
-    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
-                               3./340.*dumax**(-.5/nitergdiv)
+       ENDDO
+       CALL transfert_request(f_du,req_e1_vect)       
+       dumax = max_vel(f_du)       
+       CALL rescale_field(1./dumax, f_du, f_u)
+       IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
+    ENDDO
 
     cgraddiv=dumax**(-1./nitergdiv)
-    IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
-
-!$OMP BARRIER
-!$OMP MASTER
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      u=f_u(ind)
-
-! set random seed to get reproductibility when using a different number of process
-      CALL RANDOM_SEED(size=M)
-      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
- 
-       DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i   
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_right)=r-0.5
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_lup)=r-0.5
-          CALL RANDOM_NUMBER(r)
-          u(ij+u_ldown)=r-0.5
-        ENDDO
-      ENDDO        
-    ENDDO
-!$OMP END MASTER  
-!$OMP BARRIER
-
-
+    IF (is_master) THEN
+       PRINT *,"gradiv : dumax",dumax
+       PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
+         'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
+       PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
+         3./340.*dumax**(-.5/nitergdiv)
+       PRINT *,"cgraddiv : ",cgraddiv
+    END IF
+    
+    !----------------- Compute cgradrot --------------------
+    cgradrot=-1.
+    CALL random_vel(f_u)
     DO it=1,20
- 
-      dumax=0
-      DO iter=1,nitergrot
-        CALL transfert_request(f_u,req_e1_vect)
-        DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          CALL compute_gradrot(u,du,1,1)
-          u=du
-        ENDDO
-      ENDDO
-
-      CALL transfert_request(f_du,req_e1_vect)
-      
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        u=f_u(ind)
-        du=f_du(ind)
-        
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i   
-            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
-            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
-            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
+       DO iter=1,nitergrot
+          CALL transfert_request(f_u,req_e1_vect)
+          DO ind=1,ndomain
+             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+             CALL swap_dimensions(ind)
+             CALL swap_geometry(ind)
+             u=f_u(ind)
+             du=f_du(ind)
+             CALL compute_gradrot(u,du,1,1)
+             u=du
           ENDDO
-        ENDDO
-      ENDDO
-
-      IF (using_mpi) THEN 
-        CALL reduce_max_omp(dumax,dumax1)
-!$OMP MASTER        
-        CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
-!$OMP END MASTER     
-        CALL bcast_omp(dumax)  
-      ELSE
-        CALL allreduce_max_omp(dumax,dumax1)
-        dumax=dumax1
-      ENDIF  
-
-      
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        u=f_u(ind)
-        du=f_du(ind)
-        u=du/dumax
-      ENDDO
-      
-      IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
-
-    ENDDO 
-    IF (is_master) PRINT *,"gradrot : dumax",dumax
-  
+       ENDDO
+       CALL transfert_request(f_du,req_e1_vect)
+       dumax = max_vel(f_du)
+       CALL rescale_field(1./dumax, f_du, f_u)
+       IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
+    ENDDO
+
     cgradrot=dumax**(-1./nitergrot) 
+    IF (is_master) PRINT *,"gradrot : dumax",dumax  
     IF (is_master) PRINT *,"cgradrot : ",cgradrot
    
-
-
-!$OMP BARRIER
-!$OMP MASTER
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      theta=f_theta(ind)
-
-! set random seed to get reproductibility when using a different number of process
-      CALL RANDOM_SEED(size=M)
-      CALL RANDOM_SEED(put=(/(i,i=1,M)/))
-
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i   
-          CALL RANDOM_NUMBER(r)
-          theta(ij)=r-0.5
-        ENDDO
-      ENDDO  
-    ENDDO
-!$OMP END MASTER
-!$OMP BARRIER
-
-    DO it=1,20
-
-      dthetamax=0
-      DO iter=1,niterdivgrad
-        CALL transfert_request(f_theta,req_i1)
-        DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          theta=f_theta(ind)
-          dtheta=f_dtheta(ind)
-          CALL compute_divgrad(theta,dtheta,1,1)
-          theta=dtheta
-        ENDDO
-      ENDDO
-
-      CALL transfert_request(f_dtheta,req_i1)
-      
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        theta=f_theta(ind)
-        dtheta=f_dtheta(ind)
-        
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i   
-            dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
+    !----------------- Compute cdivgrad --------------------
+
+    cdivgrad=-1.
+    CALL random_scalar(f_theta)    
+    DO it=1,20       
+       DO iter=1,niterdivgrad
+          CALL transfert_request(f_theta,req_i1)
+          DO ind=1,ndomain
+             IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+             CALL swap_dimensions(ind)
+             CALL swap_geometry(ind)
+             theta=f_theta(ind)
+             dtheta=f_dtheta(ind)
+             CALL compute_divgrad(theta,dtheta,1,1)
           ENDDO
-        ENDDO
-      ENDDO
-
-      IF (using_mpi) THEN 
-        CALL reduce_max_omp(dthetamax ,dthetamax1)
-!$OMP MASTER        
-        CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
-!$OMP END MASTER     
-        CALL bcast_omp(dthetamax)  
-      ELSE
-        CALL allreduce_max_omp(dthetamax,dthetamax1)
-        dumax=dumax1
-      ENDIF  
-      
-      IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
-
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        theta=f_theta(ind)
-        dtheta=f_dtheta(ind)
-        theta=dtheta/dthetamax
-      ENDDO
-    ENDDO 
-
+       ENDDO       
+       CALL transfert_request(f_dtheta,req_i1)
+       dthetamax = max_scalar(f_dtheta)      
+       IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
+       CALL rescale_field(1./dthetamax, f_dtheta, f_theta)
+    END DO
+    
+    cdivgrad=dthetamax**(-1./niterdivgrad) 
     IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
-  
-    cdivgrad=dthetamax**(-1./niterdivgrad) 
     IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
 
-     
+  END SUBROUTINE dissip_constants
+
+  SUBROUTINE dissip_profile
+    USE disvert_mod
+    INTEGER    :: l
+    REAL(rstd) :: mintau, zz, zvert, fact
     fact=2
     DO l=1,llm
@@ -441,19 +246,19 @@
        tau_gradrot(l) = tau_gradrot(l)/zvert
        tau_divgrad(l) = tau_divgrad(l)/zvert
-    ENDDO
-
+    END DO
+    
     mintau=tau_graddiv(1)
     DO l=1,llm
-      mintau=MIN(mintau,tau_graddiv(l))
-      mintau=MIN(mintau,tau_gradrot(l))
-      mintau=MIN(mintau,tau_divgrad(l))
-    ENDDO
-
+       mintau=MIN(mintau,tau_graddiv(l))
+       mintau=MIN(mintau,tau_gradrot(l))
+       mintau=MIN(mintau,tau_divgrad(l))
+    END DO
+    
     IF(rayleigh_friction_type>0) mintau=MIN(mintau,rayleigh_tau)
-
+    
     IF(mintau>0) THEN
        IF (itau_dissip==0) THEN
-         IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
-         itau_dissip=INT(mintau/dt)
+          IF (is_master) PRINT *,"init_dissip: Automatic computation of itau_dissip..."
+          itau_dissip=INT(mintau/dt)
        ENDIF
     ELSE
@@ -464,21 +269,11 @@
     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
+       PRINT *,"init_dissip: rayleigh_tau",rayleigh_tau, "mintau ",mintau
+       PRINT *,"init_dissip: itau_dissip",itau_dissip," dtdissip ",dtdissip
     ENDIF
-
-  END SUBROUTINE init_dissip 
- 
+    
+  END SUBROUTINE dissip_profile
   
   SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
-  USE icosa
-  USE theta2theta_rhodz_mod
-  USE pression_mod
-  USE exner_mod
-  USE geopotential_mod
-  USE trace
-  USE time_mod
-  USE omp_para
-  IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:), f_dtheta_rhodz(:)
@@ -492,6 +287,5 @@
     REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
 
-    INTEGER :: ind, shear
-    INTEGER :: l,ij,nn
+    INTEGER :: ind, l,ij,nn
 
 !$OMP BARRIER
@@ -500,5 +294,4 @@
     CALL gradiv(f_ue,f_due_diss1)
     CALL gradrot(f_ue,f_due_diss2)
-
     CALL divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz_diss)
     
@@ -516,115 +309,103 @@
 !DIR$ SIMD
         DO ij=ij_begin,ij_end
-
             due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip 
             due(ij+u_lup,l)   = -0.5*( due_diss1(ij+u_lup,l)  /tau_graddiv(l) + due_diss2(ij+u_lup,l)  /tau_gradrot(l))*itau_dissip
             due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip
-
             dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
         ENDDO
       ENDDO
-
-!      dtheta_rhodz=0
-!      due=0
-
+      
       IF(rayleigh_friction_type>0) THEN
-       IF(rayleigh_friction_type<99) THEN
-         phi=f_geopot(ind)
-         ue=f_ue(ind)
-         DO l=ll_begin,ll_end
+         IF(rayleigh_friction_type<99) THEN
+            phi=f_geopot(ind)
+            ue=f_ue(ind)
+            DO l=ll_begin,ll_end
+               DO ij=ij_begin,ij_end
+                  CALL relax(t_right, u_right)
+                  CALL relax(t_lup,   u_lup)
+                  CALL relax(t_ldown, u_ldown)
+               ENDDO
+            END DO
+         ELSE
+            ue=f_ue(ind)
             DO ij=ij_begin,ij_end
-               CALL relax(t_right, u_right)
-               CALL relax(t_lup,   u_lup)
-               CALL relax(t_ldown, u_ldown)
+               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
-         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
+         ENDIF
       END IF
    END DO
-
+   
    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
-        REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
-             p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
-             fz, u3d(3), uref
-        REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
-        LOGICAL :: hybrid_eta
-        INTEGER :: shift_u, shift_t, zcoords, nn
-        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
-        IF(z>zh) THEN  ! relax only in the sponge layer z>zh
-
-           nn = ij+shift_u
-           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
-           CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
-                hyam,hybm,shear,ulon,ulat,w,t,phis,ps,rho,q)
-!           u3d = ulon*elon_e(nn,:) + ulat*elat_e(nn,:)
-           u3d = ulon*elon_e(nn,:) ! ulat=0
-           uref = sum(u3d*ep_e(nn,:))
-
-           fz = sin((pi/2)*(z-zh)/(ztop-zh))
-           fz = fz*fz/rayleigh_tau
-           due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
-!           fz = 1./rayleigh_tau
-!           due(nn,l) = due(nn,l) - fz*ue(nn,l)
-         END IF
-       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
-
-  SUBROUTINE gradiv(f_ue,f_due)
-  USE icosa
-  USE trace
-  USE omp_para
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ue(:)
-    TYPE(t_field),POINTER :: f_due(:)
-    REAL(rstd),POINTER    :: ue(:,:)
-    REAL(rstd),POINTER    :: due(:,:)
-    INTEGER :: ind
-    INTEGER :: it,l,ij
-       
-    CALL trace_start("gradiv")
-
-    DO ind=1,ndomain
+ CONTAINS
+   
+   SUBROUTINE relax(shift_t, shift_u)
+     USE dcmip_initial_conditions_test_1_2_3
+     REAL(rstd) :: z, ulon,ulat, lon,lat, & ! input to test2_schaer_mountain
+          p,hyam,hybm,w,t,phis,ps,rho,q, &   ! unused input/output to test2_schaer_mountain
+          fz, u3d(3), uref
+     REAL(rstd), PARAMETER :: zh=2e4,ztop=3e4  ! DCMIP values
+     LOGICAL :: hybrid_eta
+     INTEGER :: shift_u, shift_t, zcoords, nn
+     z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
+     IF(z>zh) THEN  ! relax only in the sponge layer z>zh
+        nn = ij+shift_u
+        zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
+        CALL test2_schaer_mountain(lon_e(nn),lat_e(nn),p,z,zcoords,hybrid_eta, &
+             hyam,hybm,rayleigh_shear,ulon,ulat,w,t,phis,ps,rho,q)
+        u3d = ulon*elon_e(nn,:) ! ulat=0
+        uref = sum(u3d*ep_e(nn,:))
+        
+        fz = sin((pi/2)*(z-zh)/(ztop-zh))
+        fz = fz*fz/rayleigh_tau
+        due(nn,l) = due(nn,l) - itau_dissip*fz*(ue(nn,l)-uref)
+     END IF
+   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
+
+ 
+ SUBROUTINE gradiv(f_ue,f_due)
+   TYPE(t_field), POINTER :: f_ue(:)
+   TYPE(t_field), POINTER :: f_due(:)
+   REAL(rstd),POINTER :: ue(:,:)
+   REAL(rstd),POINTER :: due(:,:)
+   INTEGER :: ind
+   INTEGER :: it,l,ij
+   
+   CALL trace_start("gradiv")
+   
+   DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
@@ -633,44 +414,35 @@
       due=f_due(ind) 
       DO  l = ll_begin, ll_end
-!DIR$ SIMD
-        DO ij=ij_begin,ij_end
-             due(ij+u_right,l)=ue(ij+u_right,l)
-             due(ij+u_lup,l)=ue(ij+u_lup,l)
-             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
-        ENDDO
-      ENDDO
-    ENDDO
-
-    DO it=1,nitergdiv
-        
+         !DIR$ SIMD
+         DO ij=ij_begin,ij_end
+            due(ij+u_right,l)=ue(ij+u_right,l)
+            due(ij+u_lup,l)=ue(ij+u_lup,l)
+            due(ij+u_ldown,l)=ue(ij+u_ldown,l)
+         ENDDO
+      ENDDO
+   ENDDO
+   
+   DO it=1,nitergdiv     
       CALL send_message(f_due,req_due)
       CALL wait_message(req_due)
-        
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        due=f_due(ind) 
-        CALL compute_gradiv(due,due,ll_begin,ll_end)
-      ENDDO
-    ENDDO
-
+         IF (.NOT. assigned_domain(ind)) CYCLE
+         CALL swap_dimensions(ind)
+         CALL swap_geometry(ind)
+         due=f_due(ind) 
+         CALL compute_gradiv(due,due,ll_begin,ll_end)
+      ENDDO
+   ENDDO
+   
    CALL trace_end("gradiv")
-
-  END SUBROUTINE gradiv
-  
-
-  SUBROUTINE gradrot(f_ue,f_due)
-  USE icosa
-  USE trace
-  USE omp_para
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ue(:)
-    TYPE(t_field),POINTER :: f_due(:)
-    REAL(rstd),POINTER    :: ue(:,:)
-    REAL(rstd),POINTER    :: due(:,:)
-    INTEGER :: ind
-    INTEGER :: it,l,ij
-       
+ END SUBROUTINE gradiv
+ 
+  
+ SUBROUTINE gradrot(f_ue,f_due)
+   TYPE(t_field), POINTER :: f_ue(:)
+    TYPE(t_field), POINTER :: f_due(:)
+    REAL(rstd),POINTER :: ue(:,:)
+    REAL(rstd),POINTER :: due(:,:)
+    INTEGER :: ind, it,l,ij
     CALL trace_start("gradrot")
 
@@ -692,33 +464,24 @@
 
     DO it=1,nitergrot
-        
-      CALL send_message(f_due,req_due)
-      CALL wait_message(req_due)
-        
-      DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-        due=f_due(ind) 
-        CALL compute_gradrot(due,due,ll_begin,ll_end)
-      ENDDO
-
-    ENDDO
-
+       CALL send_message(f_due,req_due)
+       CALL wait_message(req_due)
+       DO ind=1,ndomain
+          IF (.NOT. assigned_domain(ind)) CYCLE
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          due=f_due(ind) 
+          CALL compute_gradrot(due,due,ll_begin,ll_end)
+       ENDDO
+    ENDDO
+    
     CALL trace_end("gradrot")
-
   END SUBROUTINE gradrot
   
   SUBROUTINE divgrad(f_theta,f_dtheta)
-  USE icosa
-  USE trace
-  USE omp_para
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_theta(:)
-    TYPE(t_field),POINTER :: f_dtheta(:)
-    REAL(rstd),POINTER    :: theta(:,:)
-    REAL(rstd),POINTER    :: dtheta(:,:)
-    INTEGER :: ind
-    INTEGER :: it
+    TYPE(t_field), POINTER :: f_theta(:)
+    TYPE(t_field), POINTER :: f_dtheta(:)
+    REAL(rstd),POINTER :: theta(:,:)
+    REAL(rstd),POINTER :: dtheta(:,:)
+    INTEGER :: ind, it
 
     CALL trace_start("divgrad")
@@ -733,8 +496,6 @@
     ENDDO
 
-    DO it=1,niterdivgrad
-        
-      CALL transfert_request(f_dtheta,req_i1)
-        
+    DO it=1,niterdivgrad        
+      CALL transfert_request(f_dtheta,req_i1)        
       DO ind=1,ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
@@ -744,19 +505,13 @@
         CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
       ENDDO
-
     ENDDO
 
     CALL trace_end("divgrad")
-
   END SUBROUTINE divgrad
    
   SUBROUTINE divgrad_theta_rhodz(f_mass,f_theta_rhodz,f_dtheta_rhodz)
-  USE icosa
-  USE trace
-  USE omp_para
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_mass(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+    TYPE(t_field), POINTER :: f_mass(:)
+    TYPE(t_field), POINTER :: f_theta_rhodz(:)
+    TYPE(t_field), POINTER :: f_dtheta_rhodz(:)
     
     REAL(rstd),POINTER :: mass(:,:)
@@ -813,12 +568,8 @@
     ENDDO
 
-
     CALL trace_end("divgrad")
-
   END SUBROUTINE divgrad_theta_rhodz  
   
   SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -843,12 +594,8 @@
     DO l=llb,lle
 !DIR$ SIMD
-      DO ij=ij_begin,ij_end
- 
+      DO ij=ij_begin,ij_end 
           gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )       
-
-          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))       
-    
+          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))    
           gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
-
       ENDDO
     ENDDO
@@ -867,6 +614,4 @@
   
   SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
-  USE icosa
-  IMPLICIT NONE
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -874,26 +619,17 @@
     REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
     REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
-
-    INTEGER :: ij,l
-
-       
+    INTEGER :: ij,l       
     DO l=llb,lle
 !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
- 
-          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )        
-
-          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))       
-    
+          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )
+          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))
           grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) )
-
-      ENDDO
-    ENDDO
-    
+      ENDDO
+    ENDDO    
     
     DO l=llb,lle
 !DIR$ SIMD
       DO ij=ij_begin,ij_end
-
           divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
                              ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &  
@@ -912,9 +648,6 @@
 
   END SUBROUTINE compute_divgrad
-
     
   SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -943,11 +676,7 @@
 !DIR$ SIMD
       DO ij=ij_begin,ij_end
-  
-          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 
-         
-          gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 
-        
-          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
-        
+          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))         
+          gradrot_e(ij+u_lup,l)  =1/le(ij+u_lup)  *ne(ij,lup)  *(rot_v(ij+z_up,l)   -rot_v(ij+z_lup,l))
+          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))        
       ENDDO
     ENDDO
@@ -955,10 +684,8 @@
     DO l=llb,lle
 !DIR$ SIMD
-      DO ij=ij_begin,ij_end
-    
+      DO ij=ij_begin,ij_end    
           gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
           gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
           gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
-        
       ENDDO
     ENDDO  
@@ -966,5 +693,142 @@
   END SUBROUTINE compute_gradrot
 
+!----------------------- Utility routines ------------------
+
+  SUBROUTINE global_max(dumax)
+    USE mpi_mod
+    USE mpipara
+    REAL(rstd) :: dumax, dumax1
+    IF (using_mpi) THEN 
+       CALL reduce_max_omp(dumax,dumax1)
+       !$OMP MASTER        
+       CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
+       !$OMP END MASTER     
+       CALL bcast_omp(dumax)  
+    ELSE
+       CALL allreduce_max_omp(dumax,dumax1)
+       dumax=dumax1
+    ENDIF
+  END SUBROUTINE global_max    
+
+  FUNCTION max_vel(f_du)
+    TYPE(t_field)       :: f_du(:)
+    REAL(rstd),POINTER  :: du(:)
+    INTEGER :: ind, i,j,ij
+    REAL(rstd) :: max_vel, dumax
+
+    dumax=0.
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       du=f_du(ind)
+       
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i   
+             if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
+             if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
+             if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
+          ENDDO
+       ENDDO
+    ENDDO
+    CALL global_max(dumax)
+    max_vel=dumax
+  END FUNCTION max_vel
+
+  FUNCTION max_scalar(f_dtheta)
+    TYPE(t_field)       :: f_dtheta(:)
+    REAL(rstd),POINTER  :: dtheta(:)
+    INTEGER :: ind, i,j,ij
+    REAL(rstd) :: max_scalar, dthetamax
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+       CALL swap_dimensions(ind)
+       dtheta=f_dtheta(ind)
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
+          ENDDO
+       ENDDO
+    ENDDO
+    CALL global_max(dthetamax)
+    max_scalar=dthetamax
+  END FUNCTION max_scalar
+
+  SUBROUTINE random_vel(f_u)
+    TYPE(t_field)      :: f_u(:)
+    REAL(rstd),POINTER :: u(:)
+    INTEGER :: M, ind, i,j,ij
+    REAL(rstd) :: r
+
+    !$OMP BARRIER
+    !$OMP MASTER
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       
+       ! set random seed to get reproductibility when using a different number of process
+       CALL RANDOM_SEED(size=M)
+       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
+       
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i   
+             CALL RANDOM_NUMBER(r)
+             u(ij+u_right)=r-0.5
+             CALL RANDOM_NUMBER(r)
+             u(ij+u_lup)=r-0.5
+             CALL RANDOM_NUMBER(r)
+             u(ij+u_ldown)=r-0.5
+          ENDDO
+       ENDDO
+    ENDDO
+    !$OMP END MASTER  
+    !$OMP BARRIER
+    
+  END SUBROUTINE random_vel
+  
+  SUBROUTINE random_scalar(f_theta)
+    TYPE(t_field)      :: f_theta(:)
+    REAL(rstd),POINTER :: theta(:)
+    INTEGER :: M, ind, i,j,ij
+    REAL(rstd) :: r
+    !$OMP BARRIER
+    !$OMP MASTER
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       theta=f_theta(ind)
+       ! set random seed to get reproductibility when using a different number of process
+       CALL RANDOM_SEED(size=M)
+       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
+       
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i   
+             CALL RANDOM_NUMBER(r)
+             theta(ij)=r-0.5
+          ENDDO
+       ENDDO
+    ENDDO
+    !$OMP END MASTER
+    !$OMP BARRIER
+  END SUBROUTINE random_scalar
+  
+  SUBROUTINE rescale_field(scal, f_du, f_u)
+    TYPE(t_field)      :: f_du(:), f_u(:)
+    REAL(rstd),POINTER :: du(:), u(:)
+    REAL(rstd) :: scal
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+       CALL swap_dimensions(ind)
+       u=f_u(ind)
+       du=f_du(ind)
+       u=scal*du
+    ENDDO
+  END SUBROUTINE rescale_field
 
 END MODULE dissip_gcm_mod
-       
