Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect.f90	(revision 267)
@@ -11,9 +11,9 @@
   !==========================================================================
 
-  SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng)
+  SUBROUTINE init_advect(normal,tangent,sqrt_leng)
     IMPLICIT NONE
     REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
     REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
+    REAL(rstd),INTENT(OUT) :: sqrt_leng(iim*jjm)
 
     INTEGER :: ij
@@ -40,7 +40,7 @@
           tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
 
-          one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
-                                          sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
-                                          sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
+          sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
+                                   sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
+                                   sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
     ENDDO
 
@@ -49,10 +49,10 @@
   !=======================================================================================
 
-  SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
+  SUBROUTINE compute_gradq3d(qi_in,sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
     USE trace
     USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng_in(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: sqrt_leng_in(iim*jjm)
     REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
     REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
@@ -66,5 +66,5 @@
     INTEGER :: ij,k,ind,l
     REAL(rstd)  :: qi(iim*jjm,llm)
-    REAL(rstd)  :: one_over_sqrt_leng(iim*jjm)
+    REAL(rstd)  :: sqrt_leng(iim*jjm)
     REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
     REAL(rstd) :: detx,dety,detz,det
@@ -74,5 +74,5 @@
 
     qi=qi_in
-    one_over_sqrt_leng=one_over_sqrt_leng_in
+    sqrt_leng=sqrt_leng_in
     
     CALL trace_start("compute_gradq3d1")
@@ -252,6 +252,6 @@
              maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 
              maggrd = sqrt(maggrd) 
-             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
-             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
+             maxq_c = qi(ij,l) + maggrd*sqrt_leng(ij)
+             minq_c = qi(ij,l) - maggrd*sqrt_leng(ij)
              maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
                   qi(ij+t_rdown,l),qi(ij+t_ldown,l))
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect_tracer.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/advect_tracer.f90	(revision 267)
@@ -8,5 +8,5 @@
   TYPE(t_field),SAVE,POINTER :: f_gradq3d(:)
   TYPE(t_field),SAVE,POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
-  TYPE(t_field),SAVE,POINTER :: f_one_over_sqrt_leng(:)
+  TYPE(t_field),SAVE,POINTER :: f_sqrt_leng(:)
 
   TYPE(t_message),SAVE :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d
@@ -28,5 +28,5 @@
     REAL(rstd),POINTER :: tangent(:,:)
     REAL(rstd),POINTER :: normal(:,:)
-    REAL(rstd),POINTER :: one_over_sqrt_leng(:)
+    REAL(rstd),POINTER :: sqrt_leng(:)
     INTEGER :: ind
 
@@ -35,5 +35,5 @@
     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
     CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
-    CALL allocate_field(f_one_over_sqrt_leng,field_t,type_real, name='one_over_sqrt_leng')
+    CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng')
     CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
     CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
@@ -47,6 +47,6 @@
        normal=f_normal(ind)
        tangent=f_tangent(ind)
-       one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
-       CALL init_advect(normal,tangent,one_over_sqrt_leng)
+       sqrt_leng=f_sqrt_leng(ind)
+       CALL init_advect(normal,tangent,sqrt_leng)
     END DO
 
@@ -66,5 +66,5 @@
     TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step
 
-    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), one_over_sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
+    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:)
     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
@@ -145,6 +145,6 @@
           q       = f_q(ind)
           gradq3d = f_gradq3d(ind)
-          one_over_sqrt_leng=f_one_over_sqrt_leng(ind)
-          CALL compute_gradq3d(q(:,:,k),one_over_sqrt_leng,gradq3d,xyz_i,xyz_v)
+          sqrt_leng=f_sqrt_leng(ind)
+          CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v)
        END DO
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn.f90	(revision 267)
@@ -1,5 +1,4 @@
 MODULE caldyn_mod
   USE icosa
-  USE caldyn_gcm_mod, ONLY : caldyn_BC
   PRIVATE
   CHARACTER(LEN=255),SAVE :: caldyn_type
@@ -31,4 +30,14 @@
       
   END SUBROUTINE init_caldyn
+
+  SUBROUTINE caldyn_BC(f_phis, f_wflux)
+    USE caldyn_gcm_mod, ONLY : caldyn_gcm_BC=>caldyn_BC
+    IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_phis(:), f_wflux(:)
+    SELECT CASE (TRIM(caldyn_type))
+    CASE('gcm')
+       CALL caldyn_gcm_BC(f_phis, f_wflux)
+    END SELECT
+  END SUBROUTINE caldyn_BC
   
   SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90	(revision 267)
@@ -427,7 +427,12 @@
     ENDDO
        
-    itau_dissip=INT(mintau/dt)
+    IF(mintau>0) THEN
+       itau_dissip=INT(mintau/dt)
+       dtdissip=itau_dissip*dt
+    ELSE
+       IF (is_mpi_root) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
+       itau_dissip=100000000
+    END IF
     itau_dissip=MAX(1,itau_dissip)
-    dtdissip=itau_dissip*dt
     IF (is_mpi_root) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/getin.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/getin.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/getin.f90	(revision 267)
@@ -15,4 +15,5 @@
   USE transfert_omp_mod
   USE omp_para
+  USE mpipara
   IMPLICIT NONE
     CHARACTER(LEN=*) :: name
@@ -21,4 +22,5 @@
 !$OMP MASTER    
     CALL getin_(name,value)
+    IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', TRIM(value)
 !$OMP END MASTER
     IF (omp_in_parallel()) CALL bcast_omp(value)
@@ -29,4 +31,5 @@
   USE transfert_omp_mod
   USE omp_para
+  USE mpipara
   IMPLICIT NONE
     CHARACTER(LEN=*) :: name
@@ -35,4 +38,5 @@
 !$OMP MASTER    
     CALL getin_(name,value)
+    IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', value
 !$OMP END MASTER
     IF (omp_in_parallel()) CALL bcast_omp(value)
@@ -43,4 +47,5 @@
   USE transfert_omp_mod
   USE omp_para
+  USE mpipara
   IMPLICIT NONE
     CHARACTER(LEN=*) :: name
@@ -49,4 +54,5 @@
 !$OMP MASTER    
     CALL getin_(name,value)
+    IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name), ' = ', value
 !$OMP END MASTER
     IF (omp_in_parallel()) CALL bcast_omp(value)
@@ -57,4 +63,5 @@
   USE ioipsl, ONLY : getin_=>getin
   USE omp_para
+  USE mpipara
   USE transfert_omp_mod
   IMPLICIT NONE
@@ -64,4 +71,5 @@
 !$OMP MASTER    
     CALL getin_(name,value)
+    IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name), ' = ', value
 !$OMP END MASTER
     IF (omp_in_parallel()) CALL bcast_omp(value)
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/mpipara.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/mpipara.F90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/mpipara.F90	(revision 267)
@@ -21,9 +21,25 @@
   END INTERFACE free_mpi_buffer
 
+  PRIVATE :: getin
+
 CONTAINS
 
+  SUBROUTINE getin(name,value) ! Copied from getin.f90 to avoid circular dependency
+  USE ioipsl, ONLY : getin_=>getin
+  USE transfert_omp_mod
+  USE omp_para
+  IMPLICIT NONE
+    CHARACTER(LEN=*) :: name
+    CHARACTER(LEN=*) :: value
+
+!$OMP MASTER    
+    CALL getin_(name,value)
+    IF(is_mpi_root) PRINT *,'GETIN ',TRIM(name),' = ', TRIM(value)
+!$OMP END MASTER
+    IF (omp_in_parallel()) CALL bcast_omp(value)
+  END SUBROUTINE getin
+
   SUBROUTINE init_mpipara
   USE mpi_mod
-  USE getin_mod
 #ifdef CPP_USING_XIOS
   USE xios
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/spherical_geom.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/spherical_geom.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/spherical_geom.f90	(revision 267)
@@ -179,5 +179,5 @@
 
 
-  SUBROUTINE circumcenter(A0,B0,C0,Center)
+  SUBROUTINE circumcenter(A0,B0,C0,center)
   USE vector
   IMPLICIT NONE
@@ -185,5 +185,5 @@
     REAL(rstd), INTENT(OUT) :: Center(3)
     
-    REAL(rstd)  :: a(3),b(3),c(3)
+    REAL(rstd)  :: a(3),b(3),c(3), ac(3), ab(3), p1(3), q(3), p2(3)
     
     a=A0/sqrt(sum(A0**2))
@@ -191,7 +191,15 @@
     c=C0/sqrt(sum(C0**2))
     
-    CALL Cross_product2(b-a,c-b,center)
-    center=center/sqrt(sum(center**2))
-    
+    ab=b-a
+    ac=c-a
+    CALL Cross_product2(ab,ac,p1)
+    IF(.FALSE.) THEN ! Direct solution, round-off error
+       center=p1/norm(p1)
+    ELSE ! Two-step solution, stable
+       q = SUM(ac**2)*ab-SUM(ab**2)*ac
+       CALL Cross_product2(p1,q,p2)
+       p2 = a + p2/(2.*SUM(p1**2))
+       center = p2/norm(p2)
+    END IF
   END SUBROUTINE circumcenter
     
@@ -204,22 +212,35 @@
     REAL(rstd), INTENT(OUT) :: Centr(3)
     
-    REAL(rstd) :: p1(3),p2(3),cross(3)
-    REAL(rstd) :: norm_cross
+    REAL(rstd) :: p1(3),p2(3),cross(3), cc(3)
+    REAL(rstd) :: norm_cross, area
     INTEGER :: i,j
-      
-      Centr(:)=0
-      DO i=1,n
-        j=MOD(i,n)+1
-        p1=points(:,i)/norm(points(:,i))
-        p2=points(:,j)/norm(points(:,j))
-        CALL cross_product2(p1,p2,cross)
-        norm_cross=norm(cross)
-        if (norm_cross<1e-10) CYCLE 
-        
-        Centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross
-      ENDDO
-      
-      Centr(:)=centr(:)/norm(centr(:))
-  
+
+    centr(:)=0
+    IF(.FALSE.) THEN 
+       ! Gauss formula (subject to round-off error)
+       DO i=1,n
+          j=MOD(i,n)+1
+          p1=points(:,i)/norm(points(:,i))
+          p2=points(:,j)/norm(points(:,j))
+          CALL cross_product2(p1,p2,cross)
+          norm_cross=norm(cross)
+          if (norm_cross<1e-10) CYCLE 
+          centr(:)=centr(:)+asin(norm_cross)*cross(:)/norm_cross
+       ENDDO       
+    ELSE
+       ! Simple area-weighted average (second-order accurate, stable)
+       cc=SUM(points,2) ! arithmetic average used as center
+       cc=cc/norm(cc)
+       DO i=1,n
+          j=MOD(i,n)+1
+          p1=points(:,i)/norm(points(:,i))
+          p2=points(:,j)/norm(points(:,j))
+          CALL surf_triangle(cc,p1,p2,area)
+          centr(:)=centr(:)+area*(p1+p2+cc)
+       ENDDO
+    END IF
+    
+    centr(:)=centr(:)/norm(centr(:))
+
   END SUBROUTINE compute_centroid
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90	(revision 267)
@@ -68,8 +68,4 @@
     CALL getin('run_length',run_length)
     itaumax=run_length/dt
-    IF (is_mpi_root)  THEN
-       PRINT *, 'itaumax=',itaumax
-       PRINT *, 'itau_adv=',itau_adv, 'itau_dissip=',itau_dissip, 'itau_physics=',itau_physics
-    END IF
     dt=dt/scale_factor
 
@@ -90,4 +86,8 @@
     CALL getin('itau_physics',itau_physics)
 
+    IF (is_mpi_root)  THEN
+       PRINT *, 'itaumax=',itaumax
+       PRINT *, 'itau_adv=',itau_adv, 'itau_dissip=',itau_dissip, 'itau_physics=',itau_physics
+    END IF
     
     CALL create_time_counter_header
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/transfert_omp.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/transfert_omp.f90	(revision 266)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/transfert_omp.f90	(revision 267)
@@ -95,7 +95,8 @@
   IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(INOUT) :: Var
-    
-    CALL bcast_omp_cgen(Var,len(Var),buffer_c)
-    
+    INTEGER :: lenvar
+    lenvar=len(Var)
+    CALL bcast_omp_i(lenvar)
+    CALL bcast_omp_cgen(Var,lenvar,buffer_c)    
   END SUBROUTINE bcast_omp_c
 
