Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn_gcm.f90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn_gcm.f90	(revision 298)
@@ -33,4 +33,5 @@
     USE exner_mod
     USE mpipara
+    USE omp_para, ONLY: omp_master
     IMPLICIT NONE
     CHARACTER(len=255) :: def
@@ -50,5 +51,5 @@
        STOP
     END SELECT
-    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
+    IF (is_mpi_root .AND. omp_master) PRINT *, 'caldyn_conserv=',def
 
     CALL allocate_caldyn
@@ -294,5 +295,5 @@
     IF (write_out) THEN
 
-       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
+!       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
 
 ! ---> for openMP test to fix later
@@ -892,4 +893,5 @@
   USE icosa
   USE mpipara
+  USE omp_para, ONLY: omp_master
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -924,5 +926,5 @@
    
     ENDDO
-    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        
+    IF (is_mpi_root .AND. omp_master) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        
 
   END SUBROUTINE check_mass_conservation  
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/check_conserve.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/check_conserve.f90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/check_conserve.f90	(revision 298)
@@ -88,5 +88,5 @@
        WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
        WRITE(134,*)"=================================================="
-       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+!       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
        
 4000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  &
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/dissip_gcm.f90	(revision 298)
@@ -68,4 +68,5 @@
   USE mpi_mod
   USE mpipara
+  USE omp_para, ONLY: omp_master
   USE transfert_mod
   USE time_mod
@@ -100,15 +101,15 @@
    CASE('none')
       rayleigh_friction_type=0
-      IF (is_mpi_root) PRINT *, 'No Rayleigh friction'
+      IF (is_mpi_root .AND. omp_master) PRINT *, 'No Rayleigh friction'
    CASE('dcmip2_schaer_noshear')
       rayleigh_friction_type=1
       rayleigh_shear=0
-      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
+      IF (is_mpi_root .AND. omp_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
    CASE('dcmip2_schaer_shear')
       rayleigh_shear=1
       rayleigh_friction_type=2
-      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+      IF (is_mpi_root .AND. omp_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
    CASE DEFAULT
-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
+      IF (is_mpi_root .AND. omp_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
       STOP
    END SELECT
@@ -119,5 +120,5 @@
       rayleigh_tau = rayleigh_tau / scale_factor
       IF(rayleigh_tau<=0) THEN
-         IF (is_mpi_root) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
+         IF (is_mpi_root .AND. omp_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
          STOP
       END IF
@@ -232,15 +233,15 @@
         u=du/dumax
       ENDDO
-      IF (is_mpi_root) PRINT *,"gradiv : it :",it ,": dumax",dumax
+      IF (is_mpi_root .AND. omp_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
 
     ENDDO 
-    IF (is_mpi_root) PRINT *,"gradiv : dumax",dumax
-    IF (is_mpi_root) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
+    IF (is_mpi_root .AND. omp_master) PRINT *,"gradiv : dumax",dumax
+    IF (is_mpi_root .AND. omp_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_mpi_root)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', &
+    IF (is_mpi_root .AND. omp_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', &
                                2.8/340.*dumax**(-.5/nitergdiv)
 
     cgraddiv=dumax**(-1./nitergdiv)
-    IF (is_mpi_root) PRINT *,"cgraddiv : ",cgraddiv
+    IF (is_mpi_root .AND. omp_master) PRINT *,"cgraddiv : ",cgraddiv
 
     DO ind=1,ndomain
@@ -320,11 +321,11 @@
       ENDDO
       
-      IF (is_mpi_root) PRINT *,"gradrot : it :",it ,": dumax",dumax
+      IF (is_mpi_root .AND. omp_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
 
     ENDDO 
-    IF (is_mpi_root) PRINT *,"gradrot : dumax",dumax
+    IF (is_mpi_root .AND. omp_master) PRINT *,"gradrot : dumax",dumax
   
     cgradrot=dumax**(-1./nitergrot) 
-    IF (is_mpi_root) PRINT *,"cgradrot : ",cgradrot
+    IF (is_mpi_root .AND. omp_master) PRINT *,"cgradrot : ",cgradrot
    
 
@@ -389,5 +390,5 @@
       ENDIF  
       
-      IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
+      IF (is_mpi_root .AND. omp_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
 
       DO ind=1,ndomain
@@ -401,8 +402,8 @@
     ENDDO 
 
-    IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
+    IF (is_mpi_root .AND. omp_master) PRINT *,"divgrad : divgrad",dthetamax
   
     cdivgrad=dthetamax**(-1./niterdivgrad) 
-    IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
+    IF (is_mpi_root .AND. omp_master) PRINT *,"cdivgrad : ",cdivgrad
 
      
@@ -431,9 +432,9 @@
        dtdissip=itau_dissip*dt
     ELSE
-       IF (is_mpi_root) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
+       IF (is_mpi_root .AND. omp_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
        itau_dissip=100000000
     END IF
     itau_dissip=MAX(1,itau_dissip)
-    IF (is_mpi_root) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
+    IF (is_mpi_root .AND. omp_master) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
 
   END SUBROUTINE init_dissip 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/sponge.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/sponge.f90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/sponge.f90	(revision 298)
@@ -29,16 +29,16 @@
     tau_sponge = 1.e-4
     CALL getin("tau_sponge",tau_sponge)
-    PRINT*,'tau_sponge = ',tau_sponge
+    IF (is_mpi_master .AND. omp_master) PRINT*,'tau_sponge = ',tau_sponge
 
     iflag_sponge = 0
     CALL getin("iflag_sponge",iflag_sponge)
-    PRINT*,'iflag_sponge = ',iflag_sponge
+    IF (is_mpi_master .AND. omp_master) PRINT*,'iflag_sponge = ',iflag_sponge
     
     mode_sponge = 1
     CALL getin("mode_sponge",mode_sponge)
-    PRINT*,'mode_sponge = ',mode_sponge
+    IF (is_mpi_master .AND. omp_master) PRINT*,'mode_sponge = ',mode_sponge
 
     IF (iflag_sponge == 0) THEN
-      PRINT*,'init_sponge: no sponge'
+      IF (is_mpi_master .AND. omp_master) PRINT*,'init_sponge: no sponge'
       RETURN
     ENDIF
Index: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/time.f90	(revision 298)
@@ -51,4 +51,5 @@
   USE getin_mod
   USE mpipara
+  USE omp_para, ONLY: omp_master
   IMPLICIT NONE
   REAL(rstd) :: run_length
@@ -74,5 +75,5 @@
     write_period=write_period/scale_factor
     itau_out=FLOOR(.5+write_period/dt)
-    IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
+    IF (is_mpi_root .AND. omp_master) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
     ENDIF 
 
@@ -86,5 +87,5 @@
     CALL getin('itau_physics',itau_physics)
 
-    IF (is_mpi_root)  THEN
+    IF (is_mpi_root .AND. omp_master)  THEN
        PRINT *, 'itaumax=',itaumax
        PRINT *, 'itau_adv=',itau_adv, 'itau_dissip=',itau_dissip, 'itau_physics=',itau_physics
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F	(revision 298)
@@ -12,4 +12,5 @@
       use ioipsl_getincom
 #endif
+      use mod_phys_lmdz_para, only : is_master
 
       IMPLICIT NONE
@@ -59,5 +60,5 @@
       hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
       CALL getin('hybrid',hybrid)
-      write(lunout,*) trim(modname),': hybrid=',hybrid
+      if (is_master) write(lunout,*) trim(modname),': hybrid=',hybrid
 
 ! Ouverture possible de fichiers typiquement E.T.
@@ -82,7 +83,7 @@
 c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
 
-         write(lunout,*)'*****************************'
-         write(lunout,*)'WARNING reading esasig.def'
-         write(lunout,*)'*****************************'
+         if (is_master) write(lunout,*)'*****************************'
+         if (is_master) write(lunout,*)'WARNING reading esasig.def'
+         if (is_master) write(lunout,*)'*****************************'
          READ(99,*) scaleheight
          READ(99,*) dz0
@@ -131,7 +132,7 @@
 
       ELSE IF(ierr4.eq.0) then
-         write(lunout,*)'****************************'
-         write(lunout,*)'Reading z2sig.def'
-         write(lunout,*)'****************************'
+         if (is_master) write(lunout,*)'****************************'
+         if (is_master) write(lunout,*)'Reading z2sig.def'
+         if (is_master) write(lunout,*)'****************************'
 
          READ(99,*) scaleheight
@@ -174,7 +175,8 @@
 
       if (hybrid) then  ! use hybrid coordinates
-         write(lunout,*) "*********************************"
-         write(lunout,*) "Using hybrid vertical coordinates"
-         write(lunout,*) 
+         if (is_master) write(lunout,*) "***************************"
+         if (is_master) write(lunout,*) "Using hybrid vertical",
+     &		" coordinates"
+         if (is_master) write(lunout,*) 
 c        Coordonnees hybrides avec mod
          DO l = 1, llm
@@ -187,7 +189,8 @@
          ap(llmp1) = 0.
       else ! use sigma coordinates
-         write(lunout,*) "********************************"
-         write(lunout,*) "Using sigma vertical coordinates"
-         write(lunout,*) 
+         if (is_master) write(lunout,*) "***************************"
+         if (is_master) write(lunout,*) "Using sigma vertical",
+     &		" coordinates"
+         if (is_master) write(lunout,*) 
 c        Pour ne pas passer en coordonnees hybrides
          DO l = 1, llm
@@ -200,8 +203,8 @@
       bp(llmp1) =   0.
 
-      write(lunout,*) trim(modname),': BP '
-      write(lunout,*)  bp
-      write(lunout,*) trim(modname),': AP '
-      write(lunout,*)  ap
+      if (is_master) write(lunout,*) trim(modname),': BP '
+      if (is_master) write(lunout,*)  bp
+      if (is_master) write(lunout,*) trim(modname),': AP '
+      if (is_master) write(lunout,*)  ap
 
 c     Calcul au milieu des couches :
@@ -226,8 +229,8 @@
       end if
 
-      write(lunout,*) trim(modname),': BPs '
-      write(lunout,*)  bps
-      write(lunout,*) trim(modname),': APs'
-      write(lunout,*)  aps
+      if (is_master) write(lunout,*) trim(modname),': BPs '
+      if (is_master) write(lunout,*)  bps
+      if (is_master) write(lunout,*) trim(modname),': APs'
+      if (is_master) write(lunout,*)  aps
 
       DO l = 1, llm
@@ -236,9 +239,9 @@
       ENDDO
 
-      write(lunout,*)trim(modname),' : PRESNIVS' 
-      write(lunout,*)presnivs 
-      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
-     &                'height of ',scaleheight,' km)' 
-      write(lunout,*)pseudoalt
+      if (is_master) write(lunout,*)trim(modname),' : PRESNIVS' 
+      if (is_master) write(lunout,*)presnivs 
+      if (is_master) write(lunout,*)'Pseudo altitude of Presnivs : ',
+     &		'(for a scale height of ',scaleheight,' km)' 
+      if (is_master) write(lunout,*)pseudoalt
 
 c     --------------------------------------------------
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/inigeom.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/inigeom.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d_common/inigeom.F	(revision 298)
@@ -16,4 +16,5 @@
 c
 c
+      use mod_phys_lmdz_para, only : is_master
       IMPLICIT NONE
 c
@@ -160,5 +161,5 @@
 c
 c
-      WRITE(6,3) 
+      if (is_master) WRITE(6,3) 
  3    FORMAT( // 10x,' ....  INIGEOM  date du 01/06/98   ..... ',
      * //5x,'   Calcul des elongations cu et cv  comme sommes des 4 ' /
@@ -183,10 +184,10 @@
       ENDIF
 
-      WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
-     *  nitergdiv,nitergrot,niterh
+      if (is_master) WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,
+     *  gamdi_h,coefdis,nitergdiv,nitergrot,niterh
 c
       pi    = 2.* ASIN(1.)
 c
-      WRITE(6,990) 
+      if (is_master) WRITE(6,990) 
 
 c     ----------------------------------------------------------------
@@ -197,5 +198,6 @@
        IF( ysinus )  THEN
 c
-        WRITE(6,*) ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
+        if (is_master) WRITE(6,*) ' ***  Inigeom ,  Y = Sinus ',
+     *		'( Latitude ) *** '
 c
 c   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
@@ -207,5 +209,6 @@
        ELSE
 c
-        WRITE(6,*) '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
+        if (is_master) WRITE(6,*) '*** Inigeom ,  Y = Latitude  ,',
+     *		' der. sinusoid . ***'
 
 c  .... utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ...
@@ -262,5 +265,6 @@
 c   .....................................................................
 
-      WRITE(6,*)'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
+      if (is_master) WRITE(6,*)'*** Inigeom , Y = Latitude  ,',
+     *		' der.tg. hyperbolique ***'
  
        CALL fxyhyper( clat, grossismy, dzoomy, tauy    , 
@@ -656,4 +660,5 @@
 c-----------------------------------------------------------------------
 c
+       if (is_master) then
        WRITE(6,*) '   ***  Coordonnees de la grille  *** '
        WRITE(6,995)
@@ -661,20 +666,23 @@
        WRITE(6,*) '   LONGITUDES  aux pts.   V  ( degres )  '
        WRITE(6,995)
+       end if
         DO i=1,iip1
          rlonvv(i) = rlonv(i)*180./pi
         ENDDO
-       WRITE(6,400) rlonvv
-c
-       WRITE(6,995)
-       WRITE(6,*) '   LATITUDES   aux pts.   V  ( degres )  '
-       WRITE(6,995)
+       if (is_master) WRITE(6,400) rlonvv
+c
+       if (is_master) WRITE(6,995)
+       if (is_master) WRITE(6,*) '   LATITUDES   aux pts.   V  ',
+     *		'( degres )  '
+       if (is_master) WRITE(6,995)
         DO i=1,jjm
          rlatuu(i)=rlatv(i)*180./pi
         ENDDO
-       WRITE(6,400) (rlatuu(i),i=1,jjm)
+       if (is_master) WRITE(6,400) (rlatuu(i),i=1,jjm)
 c
         DO i=1,iip1
           rlonvv(i)=rlonu(i)*180./pi
         ENDDO
+       if (is_master) then
        WRITE(6,995)
        WRITE(6,*) '   LONGITUDES  aux pts.   U  ( degres )  '
@@ -685,9 +693,11 @@
        WRITE(6,*) '   LATITUDES   aux pts.   U  ( degres )  '
        WRITE(6,995)
+       end if
+       if (is_master) WRITE(6,995)
         DO i=1,jjp1
          rlatuu(i)=rlatu(i)*180./pi
         ENDDO
-       WRITE(6,400) (rlatuu(i),i=1,jjp1)
-       WRITE(6,995)
+       if (is_master) WRITE(6,400) (rlatuu(i),i=1,jjp1)
+       if (is_master) WRITE(6,995)
 c
 444    format(f10.3,f6.0)
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 298)
@@ -19,4 +19,5 @@
       use assert_m, only: assert
       use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
+      use mod_phys_lmdz_para, only : is_master
       IMPLICIT NONE
 !-----------------------------------------------------------------------
@@ -952,4 +953,5 @@
 
 
+      if (is_master) then
       write(lunout,*)' #########################################'
       write(lunout,*)' Configuration des parametres lus via run.def '
@@ -1007,4 +1009,5 @@
        write(lunout,*)' moyzon_ch = ', moyzon_ch
       endif
+      end if !of if (is_master)
 
       RETURN
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F	(revision 298)
@@ -11,4 +11,5 @@
       use netcdf95, only: NF95_PUT_VAR
       use control_mod, only : planet_type
+      use mod_phys_lmdz_para, only : is_master
 
       IMPLICIT NONE
@@ -81,9 +82,11 @@
       !!!     .... while keeping everything OK for LMDZ EARTH
       if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then
-          write(lunout,*) trim(modname),' : Planeto-like start file'
+          if (is_master) write(lunout,*) trim(modname),
+     &		' : Planeto-like start file'
           start_file_type="planeto"
           idecal = 4
       else
-          write(lunout,*) trim(modname),' : Earth-like start file'
+          if (is_master) write(lunout,*) trim(modname),
+     &		' : Earth-like start file'
           idecal = 5
       endif
@@ -573,6 +576,6 @@
       ierr = NF_CLOSE(nid) ! fermer le fichier
 
-      PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
-      PRINT*,'rad,omeg,g,cpp,kappa',
+      if (is_master) PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
+      if (is_master) PRINT*,'rad,omeg,g,cpp,kappa',
      ,        rad,omeg,g,cpp,kappa
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeropacity.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeropacity.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeropacity.F90	(revision 298)
@@ -6,4 +6,5 @@
        USE comgeomfi_h
        USE tracer_h, only: noms,rho_co2,rho_ice
+       use mod_phys_lmdz_para, only : is_master
                   
        implicit none
@@ -86,41 +87,41 @@
       IF (firstcall) THEN
 
-        write(*,*) "Tracers found in aeropacity:"
+        if (is_master) write(*,*) "Tracers found in aeropacity:"
         do iq=1,nq
           tracername=noms(iq)
           if (tracername.eq."co2_ice") then
             i_co2ice=iq
-          write(*,*) "i_co2ice=",i_co2ice
+          if (is_master) write(*,*) "i_co2ice=",i_co2ice
 
           endif
           if (tracername.eq."h2o_ice") then
             i_h2oice=iq
-            write(*,*) "i_h2oice=",i_h2oice
+            if (is_master) write(*,*) "i_h2oice=",i_h2oice
           endif
         enddo
 
         if (noaero) then
-          print*, "No active aerosols found in aeropacity"
+          if (is_master) print*, "No active aerosols found in aeropacity"
         else
-          print*, "If you would like to use aerosols, make sure any old"
-          print*, "start files are updated in newstart using the option"
-          print*, "q=0"
-          write(*,*) "Active aerosols found in aeropacity:"
+          if (is_master) print*, "If you would like to use aerosols, make sure any old"
+          if (is_master) print*, "start files are updated in newstart using the option"
+          if (is_master) print*, "q=0"
+          if (is_master) write(*,*) "Active aerosols found in aeropacity:"
         endif
 
         if ((iaero_co2.ne.0).and.(.not.noaero)) then
-          print*, 'iaero_co2=  ',iaero_co2
+          if (is_master) print*, 'iaero_co2=  ',iaero_co2
         endif
         if (iaero_h2o.ne.0) then
-          print*,'iaero_h2o=  ',iaero_h2o    
+          if (is_master) print*,'iaero_h2o=  ',iaero_h2o    
         endif
         if (iaero_dust.ne.0) then
-          print*,'iaero_dust= ',iaero_dust
+          if (is_master) print*,'iaero_dust= ',iaero_dust
         endif
         if (iaero_h2so4.ne.0) then
-          print*,'iaero_h2so4= ',iaero_h2so4
+          if (is_master) print*,'iaero_h2so4= ',iaero_h2so4
         endif
         if (iaero_back2lay.ne.0) then
-          print*,'iaero_back2lay= ',iaero_back2lay
+          if (is_master) print*,'iaero_back2lay= ',iaero_back2lay
         endif
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_cpp_mugaz.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_cpp_mugaz.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_cpp_mugaz.F90	(revision 298)
@@ -18,4 +18,5 @@
 
       use gases_h
+      use mod_phys_lmdz_para, only : is_master
       implicit none
 
@@ -86,11 +87,13 @@
       cpp_c = 1000.0*cpp_c
 
+      if (is_master) then
       print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
       print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
       print*,'Predefined Cp in physics is ',cpp,'J kg^-1 K^-1'
       print*,'Predefined Mg in physics is ',mugaz,'amu'
+      end if
 
       if (check_cpp_match) then
-         print*,'REQUEST TO CHECK cpp_match :'
+         if (is_master) print*,'REQUEST TO CHECK cpp_match :'
          if((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
               (abs(1.-mugaz/mugaz_c).gt.1.e-6)) then
@@ -101,14 +104,14 @@
             stop
          else
-            print*,'--> OK. Settings match composition.'
+            if (is_master) print*,'--> OK. Settings match composition.'
          endif
       endif
 
       if (.not.force_cpp) then
-          print*,'*** Setting cpp & mugaz to computations in calc_cpp_mugaz.'
+          if (is_master) print*,'*** Setting cpp & mugaz to computations in calc_cpp_mugaz.'
           mugaz = mugaz_c
           cpp = cpp_c
       else
-          print*,'*** Setting cpp & mugaz to predefined values.'
+          if (is_master) print*,'*** Setting cpp & mugaz to predefined values.'
       endif
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_rayleigh.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_rayleigh.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/calc_rayleigh.F90	(revision 298)
@@ -28,4 +28,5 @@
       use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
       use gases_h
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -55,8 +56,8 @@
       do igas=1,ngasmx
          if(igas.eq.vgas)then
-            print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
+            if (is_master) print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
 	 endif
 	 if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
-            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
+            if (is_master) print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
             'as its mixing ratio is less than 0.05.' 
             ! ignore variable gas in Rayleigh calculation
@@ -84,5 +85,5 @@
 
             if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then
-               print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
+               if (is_master) print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
                typeknown=.true.
             endif
@@ -92,5 +93,5 @@
 
       if(.not.typeknown)then
-         print*,'Rayleigh scattering is for a mixed gas atmosphere.'
+         if (is_master) print*,'Rayleigh scattering is for a mixed gas atmosphere.'
          typeknown=.true.
       endif
@@ -159,7 +160,7 @@
       IF (L_NSPECTV > 6) THEN
         icantbewrong = L_NSPECTV-6
-        print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
-        print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
-        print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
+        if (is_master) print*,'At 1 atm and lambda = ',WAVEV(icantbewrong),' um'
+        if (is_master) print*,'tau_R = ',TAURAY(icantbewrong)*1013.25
+        if (is_master) print*,'sig_R = ',TAURAY(icantbewrong)*g*mugaz*1.67e-27*100, &
                'cm^2 molecule^-1'
       ENDIF
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callcorrk.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callcorrk.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/callcorrk.F90	(revision 298)
@@ -18,4 +18,5 @@
       use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay
       USE tracer_h
+      USE mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -218,7 +219,7 @@
 !--------------------------------------------------
 !     set up correlated k
-         print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
+         if (is_master) print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
          call getin_p("corrkdir",corrkdir)
-         print*, "corrkdir = ",corrkdir
+         if (is_master) print*, "corrkdir = ",corrkdir
          write( tmp1, '(i3)' ) L_NSPECTI
          write( tmp2, '(i3)' ) L_NSPECTV
@@ -272,20 +273,20 @@
          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
 	    call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
-            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            if (is_master) print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            if (is_master) print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
 	 end if
          if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
 	    call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
                              reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            if (is_master) print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
+            if (is_master) print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
          endif
          if(iaer.eq.iaero_dust)then
 	    call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
-            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
+            if (is_master) print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          if(iaer.eq.iaero_h2so4)then
 	    call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
-            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
+            if (is_master) print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
           if(iaer.eq.iaero_back2lay)then
@@ -548,7 +549,7 @@
             muvarrad(2*nlayer+1)=muvar(ig,1)
 
-            print*,'Recalculating qvar with VARIABLE epsi for kastprof'
-            print*,'Assumes that the variable gas is H2O!!!'
-            print*,'Assumes that there is only one tracer'
+            if (is_master) print*,'Recalculating qvar with VARIABLE epsi for kastprof'
+            if (is_master) print*,'Assumes that the variable gas is H2O!!!'
+            if (is_master) print*,'Assumes that there is only one tracer'
             !i_var=igcm_h2o_vap
             i_var=1
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniaerosol.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniaerosol.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniaerosol.F	(revision 298)
@@ -4,4 +4,5 @@
       use radinc_h, only: naerkind
       use aerosol_mod
+      use mod_phys_lmdz_para, only : is_master
 
       IMPLICIT NONE
@@ -26,5 +27,5 @@
          iaero_co2=ia
       endif
-      write(*,*) '--- CO2 aerosol = ', iaero_co2
+      if (is_master) write(*,*) '--- CO2 aerosol = ', iaero_co2
  
       if (aeroh2o) then
@@ -32,5 +33,5 @@
          iaero_h2o=ia
       endif
-      write(*,*) '--- H2O aerosol = ', iaero_h2o
+      if (is_master) write(*,*) '--- H2O aerosol = ', iaero_h2o
 
       if (dusttau.gt.0) then
@@ -38,5 +39,5 @@
          iaero_dust=ia
       endif
-      write(*,*) '--- Dust aerosol = ', iaero_dust
+      if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust
 
       if (aeroh2so4) then
@@ -44,5 +45,5 @@
          iaero_h2so4=ia
       endif
-      write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
+      if (is_master) write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
       
       if (aeroback2lay) then
@@ -50,7 +51,8 @@
          iaero_back2lay=ia
       endif
-      write(*,*) '--- Two-layer aerosol = ', iaero_back2lay
+      if (is_master) write(*,*) '--- Two-layer aerosol = ',
+     &		 iaero_back2lay
 
-      write(*,*) '=== Number of aerosols= ', ia
+      if (is_master) write(*,*) '=== Number of aerosols= ', ia
       
 ! For the zero aerosol case, we currently make a dummy co2 aerosol which is zero everywhere.
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F	(revision 298)
@@ -52,4 +52,5 @@
 !      USE ioipsl_getincom 
       USE ioipsl_getincom_p
+      use mod_phys_lmdz_para, only : is_master
       IMPLICIT NONE
 !#include "dimensions.h"
@@ -89,4 +90,5 @@
       rcp=r/cpp
       ! Ehouarn debug:
+      if (is_master) then
       write(*,*) "inifis: rad=",rad
       write(*,*) "        daysec=",daysec
@@ -96,4 +98,5 @@
       write(*,*) "        r=",r
       write(*,*) "        rcp=",rcp
+      end if
 
       avocado = 6.02214179e23   ! added by RW
@@ -121,4 +124,5 @@
 !!!      IF(ierr.EQ.0) THEN
       IF(iscallphys) THEN
+         if (is_master) then
          PRINT*
          PRINT*
@@ -126,231 +130,254 @@
          PRINT*,' inifis: Parametres pour la physique (callphys.def)'
          PRINT*,'--------------------------------------------'
-
-         write(*,*) "Directory where external input files are:"
+         end if
+
+         if (is_master) write(*,*) "Directory where external input",
+     &		" files are:"
+
          ! default 'datadir' is set in "datadir_mod"
          call getin_p("datadir",datadir) ! default path
-         write(*,*) " datadir = ",trim(datadir)
-
-         write(*,*) "Run with or without tracer transport ?"
+         if (is_master) write(*,*) " datadir = ",trim(datadir)
+
+         if (is_master) write(*,*) "Run with or without tracer",
+     &		" transport ?"
          tracer=.false. ! default value
          call getin_p("tracer",tracer)
-         write(*,*) " tracer = ",tracer
-
-         write(*,*) "Run with or without atm mass update ",
-     &      " due to tracer evaporation/condensation?"
+         if (is_master) write(*,*) " tracer = ",tracer
+
+         if (is_master) write(*,*) "Run with or without atm mass",
+     &		" update due to tracer evaporation/condensation?"
          mass_redistrib=.false. ! default value
          call getin_p("mass_redistrib",mass_redistrib)
-         write(*,*) " mass_redistrib = ",mass_redistrib
-
-         write(*,*) "Diurnal cycle ?"
-         write(*,*) "(if diurnal=false, diurnal averaged solar heating)"
+         if (is_master) write(*,*) " mass_redistrib = ",mass_redistrib
+
+         if (is_master) write(*,*) "Diurnal cycle ?"
+         if (is_master) write(*,*) "(if diurnal=false, diurnal",
+     &		" averaged solar heating)"
          diurnal=.true. ! default value
          call getin_p("diurnal",diurnal)
-         write(*,*) " diurnal = ",diurnal
-
-         write(*,*) "Seasonal cycle ?"
-         write(*,*) "(if season=false, Ls stays constant, to value ",
-     &   "set in 'start'"
+         if (is_master) write(*,*) " diurnal = ",diurnal
+
+         if (is_master) write(*,*) "Seasonal cycle ?"
+         if (is_master) write(*,*) "(if season=false, Ls stays",
+     &		" constant, to value set in 'start'"
          season=.true. ! default value
          call getin_p("season",season)
-         write(*,*) " season = ",season
-
-         write(*,*) "Tidally resonant rotation ?"
+         if (is_master) write(*,*) " season = ",season
+
+         if (is_master) write(*,*) "Tidally resonant rotation ?"
          tlocked=.false. ! default value
          call getin_p("tlocked",tlocked)
-         write(*,*) "tlocked = ",tlocked
-
-         write(*,*) "Saturn ring shadowing ?"
+         if (is_master) write(*,*) "tlocked = ",tlocked
+
+         if (is_master) write(*,*) "Saturn ring shadowing ?"
          rings_shadow = .false.
          call getin_p("rings_shadow", rings_shadow)
-         write(*,*) "rings_shadow = ", rings_shadow
+         if (is_master) write(*,*) "rings_shadow = ", rings_shadow
          
-         write(*,*) "Compute latitude-dependent gravity field?"
+         if (is_master) write(*,*) "Compute latitude-dependent",
+     &		" gravity field?"
          oblate = .false.
          call getin_p("oblate", oblate)
-         write(*,*) "oblate = ", oblate
-
-         write(*,*) "Flattening of the planet (a-b)/a "
+         if (is_master) write(*,*) "oblate = ", oblate
+
+         if (is_master) write(*,*) "Flattening of the planet (a-b)/a "
          flatten = 0.0
          call getin_p("flatten", flatten)
-         write(*,*) "flatten = ", flatten
+         if (is_master) write(*,*) "flatten = ", flatten
          
 
-         write(*,*) "Needed if oblate=.true.: J2"
+         if (is_master) write(*,*) "Needed if oblate=.true.: J2"
          J2 = 0.0
          call getin_p("J2", J2)
-         write(*,*) "J2 = ", J2
+         if (is_master) write(*,*) "J2 = ", J2
          
-         write(*,*) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
+         if (is_master) write(*,*) "Needed if oblate=.true.: Planet",
+     &		" mass (*1e24 kg)"
          MassPlanet = 0.0
          call getin_p("MassPlanet", MassPlanet)
-         write(*,*) "MassPlanet = ", MassPlanet         
-
-         write(*,*) "Needed if oblate=.true.: Planet mean radius (m)"
+         if (is_master) write(*,*) "MassPlanet = ", MassPlanet         
+
+         if (is_master) write(*,*) "Needed if oblate=.true.: Planet",
+     &		" mean radius (m)"
          Rmean = 0.0
          call getin_p("Rmean", Rmean)
-         write(*,*) "Rmean = ", Rmean
+         if (is_master) write(*,*) "Rmean = ", Rmean
          
 ! Test of incompatibility:
 ! if tlocked, then diurnal should be false
          if (tlocked.and.diurnal) then
-           print*,'If diurnal=true, we should turn off tlocked.'
+           if (is_master) print*,'If diurnal=true, we should turn off',
+     &		' tlocked.'
            stop
          endif
 
-         write(*,*) "Tidal resonance ratio ?"
+         if (is_master) write(*,*) "Tidal resonance ratio ?"
          nres=0          ! default value
          call getin_p("nres",nres)
-         write(*,*) "nres = ",nres
-
-         write(*,*) "Write some extra output to the screen ?"
+         if (is_master) write(*,*) "nres = ",nres
+
+         if (is_master) write(*,*) "Write some extra output to the",
+     &		" screen ?"
          lwrite=.false. ! default value
          call getin_p("lwrite",lwrite)
-         write(*,*) " lwrite = ",lwrite
-
-         write(*,*) "Save statistics in file stats.nc ?"
+         if (is_master) write(*,*) " lwrite = ",lwrite
+
+         if (is_master) write(*,*) "Save statistics in file stats.nc ?"
          callstats=.true. ! default value
          call getin_p("callstats",callstats)
-         write(*,*) " callstats = ",callstats
-
-         write(*,*) "Test energy conservation of model physics ?"
+         if (is_master) write(*,*) " callstats = ",callstats
+
+         if (is_master) write(*,*) "Test energy conservation of model",
+     &		" physics ?"
          enertest=.false. ! default value
          call getin_p("enertest",enertest)
-         write(*,*) " enertest = ",enertest
-
-         write(*,*) "Check to see if cpp values used match gases.def ?"
+         if (is_master) write(*,*) " enertest = ",enertest
+
+         if (is_master) write(*,*) "Check to see if cpp values used",
+     &		" match gases.def ?"
          check_cpp_match=.true. ! default value
          call getin_p("check_cpp_match",check_cpp_match)
-         write(*,*) " check_cpp_match = ",check_cpp_match
-
-         write(*,*) "call radiative transfer ?"
+         if (is_master) write(*,*) " check_cpp_match = ",check_cpp_match
+
+         if (is_master) write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value
          call getin_p("callrad",callrad)
-         write(*,*) " callrad = ",callrad
-
-         write(*,*) "call correlated-k radiative transfer ?"
+         if (is_master) write(*,*) " callrad = ",callrad
+
+         if (is_master) write(*,*) "call correlated-k radiative",
+     &		" transfer ?"
          corrk=.true. ! default value
          call getin_p("corrk",corrk)
-         write(*,*) " corrk = ",corrk
-
-         write(*,*) "prohibit calculations outside corrk T grid?"
+         if (is_master) write(*,*) " corrk = ",corrk
+
+         if (is_master) write(*,*) "prohibit calculations outside",
+     &		" corrk T grid?"
          strictboundcorrk=.true. ! default value
          call getin_p("strictboundcorrk",strictboundcorrk)
-         write(*,*) "strictboundcorrk = ",strictboundcorrk
-
-         write(*,*) "call gaseous absorption in the visible bands?",
-     &              "(matters only if callrad=T)"
+         if (is_master) write(*,*) "strictboundcorrk = ",
+     &		strictboundcorrk
+
+         if (is_master) write(*,*) "call gaseous absorption in the",
+     &		" visible bands? (matters only if callrad=T)"
          callgasvis=.false. ! default value
          call getin_p("callgasvis",callgasvis)
-         write(*,*) " callgasvis = ",callgasvis
+         if (is_master) write(*,*) " callgasvis = ",callgasvis
         
-         write(*,*) "call continuum opacities in radiative transfer ?",
-     &              "(matters only if callrad=T)"
+         if (is_master) write(*,*) "call continuum opacities in",
+     &		" radiative transfer ? (matters only if callrad=T)"
          continuum=.true. ! default value
          call getin_p("continuum",continuum)
-         write(*,*) " continuum = ",continuum
-
-         write(*,*) "use analytic function for H2O continuum ?"
+         if (is_master) write(*,*) " continuum = ",continuum
+
+         if (is_master) write(*,*) "use analytic function for H2O",
+     &		" continuum ?"
          H2Ocont_simple=.false. ! default value
          call getin_p("H2Ocont_simple",H2Ocont_simple)
-         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
+         if (is_master) write(*,*) " H2Ocont_simple = ",H2Ocont_simple
  
-         write(*,*) "call turbulent vertical diffusion ?"
+         if (is_master) write(*,*) "call turbulent vertical",
+     &		" diffusion ?"
          calldifv=.true. ! default value
          call getin_p("calldifv",calldifv)
-         write(*,*) " calldifv = ",calldifv
-
-         write(*,*) "use turbdiff instead of vdifc ?"
+         if (is_master) write(*,*) " calldifv = ",calldifv
+
+         if (is_master) write(*,*) "use turbdiff instead of vdifc ?"
          UseTurbDiff=.true. ! default value
          call getin_p("UseTurbDiff",UseTurbDiff)
-         write(*,*) " UseTurbDiff = ",UseTurbDiff
-
-         write(*,*) "call convective adjustment ?"
+         if (is_master) write(*,*) " UseTurbDiff = ",UseTurbDiff
+
+         if (is_master) write(*,*) "call convective adjustment ?"
          calladj=.true. ! default value
          call getin_p("calladj",calladj)
-         write(*,*) " calladj = ",calladj
-
-         write(*,*) "call CO2 condensation ?"
+         if (is_master) write(*,*) " calladj = ",calladj
+
+         if (is_master) write(*,*) "call CO2 condensation ?"
          co2cond=.false. ! default value
          call getin_p("co2cond",co2cond)
-         write(*,*) " co2cond = ",co2cond
+         if (is_master) write(*,*) " co2cond = ",co2cond
 ! Test of incompatibility
          if (co2cond.and.(.not.tracer)) then
-            print*,'We need a CO2 ice tracer to condense CO2'
+            if (is_master) print*,'We need a CO2 ice tracer to',
+     &		' condense CO2'
             call abort
          endif 
  
-         write(*,*) "CO2 supersaturation level ?"
+         if (is_master) write(*,*) "CO2 supersaturation level ?"
          co2supsat=1.0 ! default value
          call getin_p("co2supsat",co2supsat)
-         write(*,*) " co2supsat = ",co2supsat
-
-         write(*,*) "Radiative timescale for Newtonian cooling ?"
+         if (is_master) write(*,*) " co2supsat = ",co2supsat
+
+         if (is_master) write(*,*) "Radiative timescale for Newtonian",
+     &		" cooling ?"
          tau_relax=30. ! default value
          call getin_p("tau_relax",tau_relax)
-         write(*,*) " tau_relax = ",tau_relax
+         if (is_master) write(*,*) " tau_relax = ",tau_relax
          tau_relax=tau_relax*24*3600 ! convert Earth days --> seconds
 
-         write(*,*)"call thermal conduction in the soil ?"
+         if (is_master) write(*,*)"call thermal conduction in",
+     &		" the soil ?"
          callsoil=.true. ! default value
          call getin_p("callsoil",callsoil)
-         write(*,*) " callsoil = ",callsoil
+         if (is_master) write(*,*) " callsoil = ",callsoil
          
-         write(*,*)"Rad transfer is computed every iradia",
-     &             " physical timestep"
+         if (is_master) write(*,*)"Rad transfer is computed",
+     &		" every iradia physical timestep"
          iradia=1 ! default value
          call getin_p("iradia",iradia)
-         write(*,*)" iradia = ",iradia
+         if (is_master) write(*,*)" iradia = ",iradia
        
-         write(*,*)"Rayleigh scattering ?"
+         if (is_master) write(*,*)"Rayleigh scattering ?"
          rayleigh=.false.
          call getin_p("rayleigh",rayleigh)
-         write(*,*)" rayleigh = ",rayleigh
-
-         write(*,*) "Use blackbody for stellar spectrum ?"
+         if (is_master) write(*,*)" rayleigh = ",rayleigh
+
+         if (is_master) write(*,*) "Use blackbody for stellar",
+     & 		" spectrum ?"
          stelbbody=.false. ! default value
          call getin_p("stelbbody",stelbbody)
-         write(*,*) " stelbbody = ",stelbbody
-
-         write(*,*) "Stellar blackbody temperature ?"
+         if (is_master) write(*,*) " stelbbody = ",stelbbody
+
+         if (is_master) write(*,*) "Stellar blackbody temperature ?"
          stelTbb=5800.0 ! default value
          call getin_p("stelTbb",stelTbb)
-         write(*,*) " stelTbb = ",stelTbb
-
-         write(*,*)"Output mean OLR in 1D?"
+         if (is_master) write(*,*) " stelTbb = ",stelTbb
+
+         if (is_master) write(*,*)"Output mean OLR in 1D?"
          meanOLR=.false.
          call getin_p("meanOLR",meanOLR)
-         write(*,*)" meanOLR = ",meanOLR
-
-         write(*,*)"Output spectral OLR in 3D?"
+         if (is_master) write(*,*)" meanOLR = ",meanOLR
+
+         if (is_master) write(*,*)"Output spectral OLR in 3D?"
          specOLR=.false.
          call getin_p("specOLR",specOLR)
-         write(*,*)" specOLR = ",specOLR
-
-         write(*,*)"Operate in kastprof mode?"
+         if (is_master) write(*,*)" specOLR = ",specOLR
+
+         if (is_master) write(*,*)"Operate in kastprof mode?"
          kastprof=.false.
          call getin_p("kastprof",kastprof)
-         write(*,*)" kastprof = ",kastprof
-
-         write(*,*)"Uniform absorption in radiative transfer?"
+         if (is_master) write(*,*)" kastprof = ",kastprof
+
+         if (is_master) write(*,*)"Uniform absorption in",
+     &		" radiative transfer?"
          graybody=.false.
          call getin_p("graybody",graybody)
-         write(*,*)" graybody = ",graybody
+         if (is_master) write(*,*)" graybody = ",graybody
 
 ! Slab Ocean 
-         write(*,*) "Use slab-ocean ?"
+         if (is_master) write(*,*) "Use slab-ocean ?"
          ok_slab_ocean=.false.         ! default value
          call getin_p("ok_slab_ocean",ok_slab_ocean)
-         write(*,*) "ok_slab_ocean = ",ok_slab_ocean
-
-         write(*,*) "Use slab-sea-ice ?"
+         if (is_master) write(*,*) "ok_slab_ocean = ",ok_slab_ocean
+
+         if (is_master) write(*,*) "Use slab-sea-ice ?"
          ok_slab_sic=.true.         ! default value
          call getin_p("ok_slab_sic",ok_slab_sic)
-         write(*,*) "ok_slab_sic = ",ok_slab_sic
-
-         write(*,*) "Use heat transport for the ocean ?"
+         if (is_master) write(*,*) "ok_slab_sic = ",ok_slab_sic
+
+         if (is_master) write(*,*) "Use heat transport for the ocean ?"
          ok_slab_heat_transp=.true.   ! default value
          call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
-         write(*,*) "ok_slab_heat_transp = ",ok_slab_heat_transp
+         if (is_master) write(*,*) "ok_slab_heat_transp = ",
+     &		ok_slab_heat_transp
 
 
@@ -359,94 +386,99 @@
 ! if kastprof used, we must be in 1D
          if (kastprof.and.(ngrid.gt.1)) then
-           print*,'kastprof can only be used in 1D!'
+           if (is_master) print*,'kastprof can only be used in 1D!'
            call abort
          endif
 
-         write(*,*)"Stratospheric temperature for kastprof mode?"
+         if (is_master) write(*,*)"Stratospheric temperature",
+     &		" for kastprof mode?"
          Tstrat=167.0
          call getin_p("Tstrat",Tstrat)
-         write(*,*)" Tstrat = ",Tstrat
-
-         write(*,*)"Remove lower boundary?"
+         if (is_master) write(*,*)" Tstrat = ",Tstrat
+
+         if (is_master) write(*,*)"Remove lower boundary?"
          nosurf=.false.
          call getin_p("nosurf",nosurf)
-         write(*,*)" nosurf = ",nosurf
+         if (is_master) write(*,*)" nosurf = ",nosurf
 
 ! Tests of incompatibility:
          if (nosurf.and.callsoil) then
-           print*,'nosurf not compatible with soil scheme!'
-           print*,'... got to make a choice!'
+           if (is_master) print*,'nosurf not compatible with soil',
+     &		' scheme! ... got to make a choice!'
            call abort
          endif
 
-         write(*,*)"Add an internal heat flux?",
+         if (is_master) write(*,*)"Add an internal heat flux?",
      .             "... matters only if callsoil=F"
          intheat=0.
          call getin_p("intheat",intheat)
-         write(*,*)" intheat = ",intheat
-
-         write(*,*)"Use Newtonian cooling for radiative transfer?"
+         if (is_master) write(*,*)" intheat = ",intheat
+
+         if (is_master) write(*,*)"Use Newtonian cooling for",
+     &		" radiative transfer?"
          newtonian=.false.
          call getin_p("newtonian",newtonian)
-         write(*,*)" newtonian = ",newtonian
+         if (is_master) write(*,*)" newtonian = ",newtonian
 
 ! Tests of incompatibility:
          if (newtonian.and.corrk) then
-           print*,'newtonian not compatible with correlated-k!'
+           if (is_master) print*,'newtonian not compatible with ',
+     &		'correlated-k!'
            call abort
          endif
          if (newtonian.and.calladj) then
-           print*,'newtonian not compatible with adjustment!'
+           if (is_master) print*,'newtonian not compatible with ',
+     &		'adjustment!'
            call abort
          endif
          if (newtonian.and.calldifv) then
-           print*,'newtonian not compatible with a boundary layer!'
+           if (is_master) print*,'newtonian not compatible with a ',
+     &		'boundary layer!'
            call abort
          endif
 
-         write(*,*)"Test physics timescale in 1D?"
+         if (is_master) write(*,*)"Test physics timescale in 1D?"
          testradtimes=.false.
          call getin_p("testradtimes",testradtimes)
-         write(*,*)" testradtimes = ",testradtimes
+         if (is_master) write(*,*)" testradtimes = ",testradtimes
 
 ! Test of incompatibility:
 ! if testradtimes used, we must be in 1D
          if (testradtimes.and.(ngrid.gt.1)) then
-           print*,'testradtimes can only be used in 1D!'
+           if (is_master) print*,'testradtimes can only be used in 1D!'
            call abort
          endif
 
-         write(*,*)"Default planetary temperature?"
+         if (is_master) write(*,*)"Default planetary temperature?"
          tplanet=215.0
          call getin_p("tplanet",tplanet)
-         write(*,*)" tplanet = ",tplanet
-
-         write(*,*)"Which star?"
+         if (is_master) write(*,*)" tplanet = ",tplanet
+
+         if (is_master) write(*,*)"Which star?"
          startype=1 ! default value = Sol
          call getin_p("startype",startype)
-         write(*,*)" startype = ",startype
-
-         write(*,*)"Value of stellar flux at 1 AU?"
+         if (is_master) write(*,*)" startype = ",startype
+
+         if (is_master) write(*,*)"Value of stellar flux at 1 AU?"
          Fat1AU=1356.0 ! default value = Sol today
          call getin_p("Fat1AU",Fat1AU)
-         write(*,*)" Fat1AU = ",Fat1AU
+         if (is_master) write(*,*)" Fat1AU = ",Fat1AU
 
 
 ! TRACERS:
 
-         write(*,*)"Varying H2O cloud fraction?"
+         if (is_master) write(*,*)"Varying H2O cloud fraction?"
          CLFvarying=.false.     ! default value
          call getin_p("CLFvarying",CLFvarying)
-         write(*,*)" CLFvarying = ",CLFvarying
-
-         write(*,*)"Value of fixed H2O cloud fraction?"
+         if (is_master) write(*,*)" CLFvarying = ",CLFvarying
+
+         if (is_master) write(*,*)"Value of fixed H2O cloud fraction?"
          CLFfixval=1.0                ! default value
          call getin_p("CLFfixval",CLFfixval)
-         write(*,*)" CLFfixval = ",CLFfixval
-
-         write(*,*)"fixed radii for Cloud particles?"
+         if (is_master) write(*,*)" CLFfixval = ",CLFfixval
+
+         if (is_master) write(*,*)"fixed radii for Cloud particles?"
          radfixed=.false. ! default value
          call getin_p("radfixed",radfixed)
-         write(*,*)" radfixed = ",radfixed
+         if (is_master) write(*,*)" radfixed = ",radfixed
 
          if(kastprof)then
@@ -454,8 +486,8 @@
          endif  
 
-	 write(*,*)"Number mixing ratio of CO2 ice particles:"
+	 if (is_master) write(*,*)"Number mixing ratio of CO2 ice particles:"
          Nmix_co2=1.e6 ! default value
          call getin_p("Nmix_co2",Nmix_co2)
-         write(*,*)" Nmix_co2 = ",Nmix_co2
+         if (is_master) write(*,*)" Nmix_co2 = ",Nmix_co2
 
 !         write(*,*)"Number of radiatively active aerosols:"
@@ -464,107 +496,120 @@
 !         write(*,*)" naerkind = ",naerkind
 
-         write(*,*)"Opacity of dust (if used):"
+         if (is_master) write(*,*)"Opacity of dust (if used):"
          dusttau=0. ! default value
          call getin_p("dusttau",dusttau)
-         write(*,*)" dusttau = ",dusttau
-
-         write(*,*)"Radiatively active CO2 aerosols?"
+         if (is_master) write(*,*)" dusttau = ",dusttau
+
+         if (is_master) write(*,*)"Radiatively active CO2 aerosols?"
          aeroco2=.false.     ! default value
          call getin_p("aeroco2",aeroco2)
-         write(*,*)" aeroco2 = ",aeroco2
-
-         write(*,*)"Fixed CO2 aerosol distribution?"
+         if (is_master) write(*,*)" aeroco2 = ",aeroco2
+
+         if (is_master) write(*,*)"Fixed CO2 aerosol distribution?"
          aerofixco2=.false.     ! default value
          call getin_p("aerofixco2",aerofixco2)
-         write(*,*)" aerofixco2 = ",aerofixco2
-
-         write(*,*)"Radiatively active water ice?"
+         if (is_master) write(*,*)" aerofixco2 = ",aerofixco2
+
+         if (is_master) write(*,*)"Radiatively active water ice?"
          aeroh2o=.false.     ! default value
          call getin_p("aeroh2o",aeroh2o)
-         write(*,*)" aeroh2o = ",aeroh2o
-
-         write(*,*)"Fixed H2O aerosol distribution?"
+         if (is_master) write(*,*)" aeroh2o = ",aeroh2o
+
+         if (is_master) write(*,*)"Fixed H2O aerosol distribution?"
          aerofixh2o=.false.     ! default value
          call getin_p("aerofixh2o",aerofixh2o)
-         write(*,*)" aerofixh2o = ",aerofixh2o
-
-         write(*,*)"Radiatively active sulfuric acid aersols?"
+         if (is_master) write(*,*)" aerofixh2o = ",aerofixh2o
+
+         if (is_master) write(*,*)"Radiatively active sulfuric",
+     &		" acid aersols?"
          aeroh2so4=.false.     ! default value
          call getin_p("aeroh2so4",aeroh2so4)
-         write(*,*)" aeroh2so4 = ",aeroh2so4
+         if (is_master) write(*,*)" aeroh2so4 = ",aeroh2so4
 	 
 !=================================
 
-	 write(*,*)"Radiatively active two-layer aersols?"
+	 if (is_master) write(*,*)"Radiatively active two-layer aersols?"
          aeroback2lay=.false.     ! default value
          call getin_p("aeroback2lay",aeroback2lay)
-         write(*,*)" aeroback2lay = ",aeroback2lay
-
-         write(*,*)"TWOLAY AEROSOL: total optical depth ",
+         if (is_master) write(*,*)" aeroback2lay = ",aeroback2lay
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: total optical depth ",
      &              "in the tropospheric layer (visible)"
          obs_tau_col_tropo=8.D0
          call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
-         write(*,*)" obs_tau_col_tropo = ",obs_tau_col_tropo
-
-         write(*,*)"TWOLAY AEROSOL: total optical depth ",
+         if (is_master) write(*,*)" obs_tau_col_tropo = ",
+     &		obs_tau_col_tropo
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: total optical depth ",
      &              "in the stratospheric layer (visible)"
          obs_tau_col_strato=0.08D0
          call getin_p("obs_tau_col_strato",obs_tau_col_strato)
-         write(*,*)" obs_tau_col_strato = ",obs_tau_col_strato
-
-         write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
+         if (is_master) write(*,*)" obs_tau_col_strato = ",
+     &		obs_tau_col_strato
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: pres_bottom_tropo?",
+     &		" in pa"
          pres_bottom_tropo=66000.0
          call getin_p("pres_bottom_tropo",pres_bottom_tropo)
-         write(*,*)" pres_bottom_tropo = ",pres_bottom_tropo
-
-         write(*,*)"TWOLAY AEROSOL: pres_top_tropo? in pa"
+         if (is_master) write(*,*)" pres_bottom_tropo = ",
+     &		pres_bottom_tropo
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: pres_top_tropo?",
+     &		" in pa"
          pres_top_tropo=18000.0
          call getin_p("pres_top_tropo",pres_top_tropo)
-         write(*,*)" pres_top_tropo = ",pres_top_tropo
-
-         write(*,*)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
+         if (is_master) write(*,*)" pres_top_tropo = ",pres_top_tropo
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: pres_bottom_strato?",
+     &		" in pa"
          pres_bottom_strato=2000.0
          call getin_p("pres_bottom_strato",pres_bottom_strato)
-         write(*,*)" pres_bottom_strato = ",pres_bottom_strato
-
-         write(*,*)"TWOLAY AEROSOL: pres_top_strato? in pa"
+         if (is_master) write(*,*)" pres_bottom_strato = ",
+     &		pres_bottom_strato
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: pres_top_strato?",
+     &		" in pa"
          pres_top_strato=100.0
          call getin_p("pres_top_strato",pres_top_strato)
-         write(*,*)" pres_top_strato = ",pres_top_strato
-
-         write(*,*)"TWOLAY AEROSOL: particle size in the ",
-     &              "tropospheric layer, in meters"
+         if (is_master) write(*,*)" pres_top_strato = ",pres_top_strato
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: particle size",
+     &		" in the tropospheric layer, in meters"
          size_tropo=2.e-6
          call getin_p("size_tropo",size_tropo)
-         write(*,*)" size_tropo = ",size_tropo
-
-         write(*,*)"TWOLAY AEROSOL: particle size in the ",
-     &              "stratospheric layer, in meters"
+         if (is_master) write(*,*)" size_tropo = ",size_tropo
+
+         if (is_master) write(*,*)"TWOLAY AEROSOL: particle size",
+     &		 "in the stratospheric layer, in meters"
          size_strato=1.e-7
          call getin_p("size_strato",size_strato)
-         write(*,*)" size_strato = ",size_strato
+         if (is_master) write(*,*)" size_strato = ",size_strato
 
 !=================================
 
-         write(*,*)"Cloud pressure level (with kastprof only):"
+         if (is_master) write(*,*)"Cloud pressure level (with ",
+     &		"kastprof only):"
          cloudlvl=0. ! default value
          call getin_p("cloudlvl",cloudlvl)
-         write(*,*)" cloudlvl = ",cloudlvl
-
-         write(*,*)"Is the variable gas species radiatively active?"
+         if (is_master) write(*,*)" cloudlvl = ",cloudlvl
+
+         if (is_master) write(*,*)"Is the variable gas species",
+     &		" radiatively active?"
          Tstrat=167.0
          varactive=.false.
          call getin_p("varactive",varactive)
-         write(*,*)" varactive = ",varactive
-
-         write(*,*)"Is the variable gas species distribution set?"
+         if (is_master) write(*,*)" varactive = ",varactive
+
+         if (is_master) write(*,*)"Is the variable gas species",
+     &		" distribution set?"
          varfixed=.false.
          call getin_p("varfixed",varfixed)
-         write(*,*)" varfixed = ",varfixed
-
-         write(*,*)"What is the saturation % of the variable species?"
+         if (is_master) write(*,*)" varfixed = ",varfixed
+
+         if (is_master) write(*,*)"What is the saturation % of",
+     &		" the variable species?"
          satval=0.8
          call getin_p("satval",satval)
-         write(*,*)" satval = ",satval
+         if (is_master) write(*,*)" satval = ",satval
 
 
@@ -572,100 +617,106 @@
 ! if varactive, then varfixed should be false
          if (varactive.and.varfixed) then
-           print*,'if varactive, varfixed must be OFF!'
+           if (is_master) print*,'if varactive, varfixed must be OFF!'
            stop
          endif
 
-         write(*,*) "Gravitationnal sedimentation ?"
+         if (is_master) write(*,*) "Gravitationnal sedimentation ?"
          sedimentation=.false. ! default value
          call getin_p("sedimentation",sedimentation)
-         write(*,*) " sedimentation = ",sedimentation
-
-         write(*,*) "Compute water cycle ?"
+         if (is_master) write(*,*) " sedimentation = ",sedimentation
+
+         if (is_master) write(*,*) "Compute water cycle ?"
          water=.false. ! default value
          call getin_p("water",water)
-         write(*,*) " water = ",water
+         if (is_master) write(*,*) " water = ",water
          
 ! Test of incompatibility:
 ! if water is true, there should be at least a tracer
          if (water.and.(.not.tracer)) then
-           print*,'if water is ON, tracer must be ON too!'
+           if (is_master) print*,'if water is ON, tracer must be ',
+     &		'ON too!'
            stop
          endif
 
-         write(*,*) "Include water condensation ?"
+         if (is_master) write(*,*) "Include water condensation ?"
          watercond=.false. ! default value
          call getin_p("watercond",watercond)
-         write(*,*) " watercond = ",watercond
+         if (is_master) write(*,*) " watercond = ",watercond
 
 ! Test of incompatibility:
 ! if watercond is used, then water should be used too
          if (watercond.and.(.not.water)) then
-           print*,'if watercond is used, water should be used too'
+           if (is_master) print*,'if watercond is used, water should ',
+     &		'be used too'
            stop
          endif
 
-         write(*,*) "Include water precipitation ?"
+         if (is_master) write(*,*) "Include water precipitation ?"
          waterrain=.false. ! default value
          call getin_p("waterrain",waterrain)
-         write(*,*) " waterrain = ",waterrain
-
-         write(*,*) "Include surface hydrology ?"
+         if (is_master) write(*,*) " waterrain = ",waterrain
+
+         if (is_master) write(*,*) "Include surface hydrology ?"
          hydrology=.false. ! default value
          call getin_p("hydrology",hydrology)
-         write(*,*) " hydrology = ",hydrology
-
-         write(*,*) "Evolve surface water sources ?"
+         if (is_master) write(*,*) " hydrology = ",hydrology
+
+         if (is_master) write(*,*) "Evolve surface water sources ?"
          sourceevol=.false. ! default value
          call getin_p("sourceevol",sourceevol)
-         write(*,*) " sourceevol = ",sourceevol
-
-         write(*,*) "Ice evolution timestep ?"
+         if (is_master) write(*,*) " sourceevol = ",sourceevol
+
+         if (is_master) write(*,*) "Ice evolution timestep ?"
          icetstep=100.0 ! default value
          call getin_p("icetstep",icetstep)
-         write(*,*) " icetstep = ",icetstep
-
-         write(*,*) "Snow albedo ?"
+         if (is_master) write(*,*) " icetstep = ",icetstep
+
+         if (is_master) write(*,*) "Snow albedo ?"
          albedosnow=0.5         ! default value
          call getin_p("albedosnow",albedosnow)
-         write(*,*) " albedosnow = ",albedosnow
-
-         write(*,*) "Maximum ice thickness ?"
+         if (is_master) write(*,*) " albedosnow = ",albedosnow
+
+         if (is_master) write(*,*) "Maximum ice thickness ?"
          maxicethick=2.0         ! default value
          call getin_p("maxicethick",maxicethick)
-         write(*,*) " maxicethick = ",maxicethick
-
-         write(*,*) "Freezing point of seawater ?"
+         if (is_master) write(*,*) " maxicethick = ",maxicethick
+
+         if (is_master) write(*,*) "Freezing point of seawater ?"
          Tsaldiff=-1.8          ! default value
          call getin_p("Tsaldiff",Tsaldiff)
-         write(*,*) " Tsaldiff = ",Tsaldiff
-
-         write(*,*) "Does user want to force cpp and mugaz?"
+         if (is_master) write(*,*) " Tsaldiff = ",Tsaldiff
+
+         if (is_master) write(*,*) "Does user want to force",
+     &		" cpp and mugaz?"
          force_cpp=.false. ! default value
          call getin_p("force_cpp",force_cpp)
-         write(*,*) " force_cpp = ",force_cpp
+         if (is_master) write(*,*) " force_cpp = ",force_cpp
 
          if (force_cpp) then
            mugaz = -99999.
-           PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
+           if (is_master) PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
            call getin_p("mugaz",mugaz)
            IF (mugaz.eq.-99999.) THEN
-               PRINT *, "mugaz must be set if force_cpp = T"
+               if (is_master) PRINT *, "mugaz must be set if ',
+     &		'force_cpp = T"
                STOP
            ELSE
-               write(*,*) "inifis: mugaz=",mugaz
+               if (is_master) write(*,*) "inifis: mugaz=",mugaz
            ENDIF
            !Ehouarn: once mugaz has been set, r, the specific
            ! gas constant must be computed
            r=8.3144622/(mugaz*1.e-3)
-           write(*,*) "inifis: r=",r
+           if (is_master) write(*,*) "inifis: r=",r
            
            cpp = -99999.
-           PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
+           if (is_master) PRINT *,'SPECIFIC HEAT CAPACITY in ',
+     &		'J K-1 kg-1 ?'
            call getin_p("cpp",cpp)
            IF (cpp.eq.-99999.) THEN
-               PRINT *, "cpp must be set if force_cpp = T"
+               if (is_master) PRINT *, "cpp must be set if ',
+     &		'force_cpp = T"
                STOP
            ELSE
-               write(*,*) "inifis: cpp=",cpp
+               if (is_master) write(*,*) "inifis: cpp=",cpp
            ENDIF
 !         else
@@ -675,10 +726,11 @@
          call calc_cpp_mugaz
 
-         PRINT*,'--------------------------------------------'
-         PRINT*
-         PRINT*
+         if (is_master) PRINT*,'------------------------------------'
+         if (is_master) PRINT*
+         if (is_master) PRINT*
       ELSE
-         write(*,*)
-         write(*,*) 'Cannot read file callphys.def. Is it here ?'
+         if (is_master) write(*,*)
+         if (is_master) write(*,*) 'Cannot read file callphys.def.",
+     &		" Is it here ?'
          stop
       ENDIF
@@ -687,4 +739,5 @@
 8001  FORMAT(t5,a12,i8)
 
+      if (is_master) then
       PRINT*
       PRINT*,'inifis: daysec',daysec
@@ -694,4 +747,5 @@
       PRINT*,'        or each ',iradia*dtphys,' seconds'
       PRINT*
+      end if
 
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniorbit.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniorbit.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/iniorbit.F	(revision 298)
@@ -4,4 +4,5 @@
       USE planete_mod, only: apoastr, periastr, year_day, obliquit,
      &                       peri_day, e_elips, p_elips, timeperi
+      USE mod_phys_lmdz_para, only : is_master
       IMPLICIT NONE
 
@@ -67,4 +68,5 @@
       !!!! We hope that all cases are above 25 Mkm [OK with Gliese 581d]
       IF ( apoastr .gt. 25.) THEN
+        IF (is_master) THEN
         PRINT*,'!!!!! WARNING !!!!!'
         PRINT*,'!!!!! YOU ARE ABOUT TO WITNESS A DIRT HACK !!!!!'
@@ -73,4 +75,5 @@
         PRINT*,'So I am assuming units are in Mkm here'
         PRINT*,'and I am performing a conversion towards AU.'
+        END IF
         periastr = periastr / 149.598 ! Mkm to AU
         apoastr = apoastr / 149.598 ! Mkm to AU
@@ -80,7 +83,7 @@
 
  
-      PRINT*,'Periastron in AU  ',periastr
-      PRINT*,'Apoastron in AU  ',apoastr 
-      PRINT*,'Obliquity in degrees  :',obliquit
+      IF (is_master) PRINT*,'Periastron in AU  ',periastr
+      IF (is_master) PRINT*,'Apoastron in AU  ',apoastr 
+      IF (is_master) PRINT*,'Obliquity in degrees  :',obliquit
 
 
@@ -88,6 +91,6 @@
       p_elips=0.5*(periastr+apoastr)*(1-e_elips*e_elips)
 
-      print*,'e_elips',e_elips
-      print*,'p_elips',p_elips
+      IF (is_master) print*,'e_elips',e_elips
+      IF (is_master) print*,'p_elips',p_elips
 
 c-----------------------------------------------------------------------
@@ -100,5 +103,5 @@
       zanom=2.*pi*(zz-nint(zz))
       zxref=abs(zanom)
-      PRINT*,'zanom  ',zanom
+      IF (is_master) PRINT*,'zanom  ',zanom
 
 c  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
@@ -114,10 +117,11 @@
       zx0=zx0+zdx
       if(zanom.lt.0.) zx0=-zx0
-      PRINT*,'zx0   ',zx0
+      IF (is_master) PRINT*,'zx0   ',zx0
 
 c zteta est la longitude solaire
 
       timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
-      PRINT*,'Solar longitude of periastron timeperi = ',timeperi
+      IF (is_master) PRINT*,'Solar longitude of periastron',
+     &		' timeperi = ',timeperi
 
       RETURN
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2H2.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2H2.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2H2.F90	(revision 298)
@@ -15,4 +15,5 @@
 
       use datafile_mod, only: datadir
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -66,6 +67,6 @@
 
       if(firstcall)then ! called by sugas_corrk only
-         print*,'----------------------------------------------------'
-         print*,'Initialising H2-H2 continuum from HITRAN database...'
+         if (is_master) print*,'----------------------------------------------------'
+         if (is_master) print*,'Initialising H2-H2 continuum from HITRAN database...'
 
 !     1.1 Open the ASCII files
@@ -106,7 +107,7 @@
 !$OMP BARRIER
 
-         print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
-         print*,'   temperature ',temp,' K'
-         print*,'   pressure ',pres,' Pa'
+         if (is_master) print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
+         if (is_master) print*,'   temperature ',temp,' K'
+         if (is_master) print*,'   pressure ',pres,' Pa'
 
       endif
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2He.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2He.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/interpolateH2He.F90	(revision 298)
@@ -15,4 +15,5 @@
 
       use datafile_mod, only: datadir
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -68,6 +69,6 @@
 
       if(firstcall)then ! called by sugas_corrk only
-         print*,'----------------------------------------------------'
-         print*,'Initialising H2-He continuum from HITRAN database...'
+         if (is_master) print*,'----------------------------------------------------'
+         if (is_master) print*,'Initialising H2-He continuum from HITRAN database...'
 
 !     1.1 Open the ASCII files
@@ -108,8 +109,10 @@
 !$OMP BARRIER
 
+         if (is_master) then 
          print*,'interpolateH2He: At wavenumber ',wn,' cm^-1'
          print*,'   temperature                 ',temp,' K'
          print*,'   H2 partial pressure         ',presH2,' Pa'
          print*,'   and He partial pressure     ',presHe,' Pa'
+         end if
 
       endif
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/phyetat0_academic.F90	(revision 298)
@@ -11,4 +11,5 @@
                      inquire_dimension, inquire_dimension_length
   use slab_ice_h, only: noceanmx
+  use mod_phys_lmdz_para, only : is_master
 
   implicit none
@@ -102,11 +103,11 @@
 ! Ehouarn, if file not found, then call tabfi with nid_start==0
 if (.not.found_file) then
-  write(*,*) 'phyetat0_academic: call tabfi with nid_start=0'
+  if (is_master) write(*,*) 'phyetat0_academic: call tabfi with nid_start=0'
   call tabfi (ngrid,0,Lmodif,tab0,day_ini,lmax,p_rad, &
                    p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
 else
   ! possibility to modify tab_cntrl in tabfi
-  write(*,*)
-  write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
+  if (is_master) write(*,*)
+  if (is_master) write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
   call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
                    p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
@@ -176,8 +177,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <phisfi>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <phisfi>"
   phisfi(:)=0
 else
-  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
+  if (is_master) write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
              minval(phisfi), maxval(phisfi)
 endif
@@ -190,10 +191,10 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <albedodat>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <albedodat>"
   do ig=1,ngrid
     albedodat(ig)=0.
   enddo
 else
-  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
+  if (is_master) write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
              minval(albedodat), maxval(albedodat)
 endif
@@ -208,5 +209,5 @@
     zmea(:)=0.
 else
-  write(*,*) "phyetat0: <ZMEA> range:", &
+  if (is_master) write(*,*) "phyetat0: <ZMEA> range:", &
              minval(zmea), maxval(zmea)
 endif
@@ -219,8 +220,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSTD>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <ZSTD>"
    zstd(:)=0.
 else
-  write(*,*) "phyetat0: <ZSTD> range:", &
+  if (is_master) write(*,*) "phyetat0: <ZSTD> range:", &
              minval(zstd), maxval(zstd)
 endif
@@ -233,8 +234,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSIG>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <ZSIG>"
   zsig(:)=0.
 else
-  write(*,*) "phyetat0: <ZSIG> range:", &
+  if (is_master) write(*,*) "phyetat0: <ZSIG> range:", &
              minval(zsig), maxval(zsig)
 endif
@@ -247,8 +248,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZGAM>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <ZGAM>"
   zgam(:)=0.
 else
-  write(*,*) "phyetat0: <ZGAM> range:", &
+  if (is_master) write(*,*) "phyetat0: <ZGAM> range:", &
              minval(zgam), maxval(zgam)
 endif
@@ -261,8 +262,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZTHE>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <ZTHE>"
   zthe(:)=0.
 else
-  write(*,*) "phyetat0: <ZTHE> range:", &
+  if (is_master) write(*,*) "phyetat0: <ZTHE> range:", &
              minval(zthe), maxval(zthe)
 endif
@@ -278,5 +279,5 @@
   !tsurf(:)=175.0
 else
-  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
+  if (is_master) write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
              minval(tsurf), maxval(tsurf)
 endif
@@ -289,8 +290,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <emis>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <emis>"
   emis(:)=0.5
 else
-  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
+  if (is_master) write(*,*) "phyetat0: Surface emissivity <emis> range:", &
              minval(emis), maxval(emis)
 endif
@@ -303,8 +304,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <cloudfrac>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <cloudfrac>"
   cloudfrac(:,:)=0.
 else
-  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
+  if (is_master) write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
              minval(cloudfrac), maxval(cloudfrac)
 endif
@@ -317,8 +318,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <totcloudfrac>"
   totcloudfrac(:)=0.5
 else
-  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
+  if (is_master) write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
              minval(totcloudfrac), maxval(totcloudfrac)
 endif
@@ -331,5 +332,5 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <hice>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <hice>"
 !  call abort
       do ig=1,ngrid
@@ -337,5 +338,5 @@
       enddo
 else
-  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
+  if (is_master) write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
              minval(hice), maxval(hice)
 endif
@@ -349,5 +350,5 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <rnat>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <rnat>"
       do ig=1,ngrid
         rnat(ig)=1.
@@ -362,5 +363,5 @@
       enddo
 
-  write(*,*) "phyetat0: Nature of surface <rnat> range:", &
+  if (is_master) write(*,*) "phyetat0: Nature of surface <rnat> range:", &
              minval(rnat), maxval(rnat)
 endif
@@ -372,10 +373,10 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <pctsrf_sic>"
       do ig=1,ngrid
       pctsrf_sic(ig)=0.
       enddo
 else
-  write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
+  if (is_master) write(*,*) "phyetat0: Pourcentage of sea ice cover <pctsrf_sic> range:", &
              minval(pctsrf_sic), maxval(pctsrf_sic)
 endif
@@ -387,5 +388,5 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tslab>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <tslab>"
       do ig=1,ngrid
       do iq=1,noceanmx
@@ -394,5 +395,5 @@
       enddo
 else
-  write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
+  if (is_master) write(*,*) "phyetat0: Slab ocean temperature <tslab> range:", &
              minval(tslab), maxval(tslab)
 endif
@@ -404,10 +405,10 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tsea_ice>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <tsea_ice>"
       do ig=1,ngrid
       tsea_ice(ig)=273.15-1.8
       enddo
 else
-  write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
+  if (is_master) write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
              minval(tsea_ice), maxval(tsea_ice)
 endif
@@ -419,10 +420,10 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <sea_ice>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <sea_ice>"
       do ig=1,ngrid
       tsea_ice(ig)=0.
       enddo
 else
-  write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
+  if (is_master) write(*,*) "phyetat0: Oceanic ice quantity <sea_ice> range:", &
              minval(sea_ice), maxval(sea_ice)
 endif
@@ -438,8 +439,8 @@
 endif
 if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <q2>"
+  if (is_master) write(*,*) "phyetat0: Failed loading <q2>"
   q2(:,:)=0.001
 else
-  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
+  if (is_master) write(*,*) "phyetat0: PBL wind variance <q2> range:", &
              minval(q2), maxval(q2)
 endif
@@ -453,5 +454,5 @@
       ! "h2o_ice" should be loaded instead
       txt="h2o_ice"
-      write(*,*) 'phyetat0: loading surface tracer', &
+      if (is_master) write(*,*) 'phyetat0: loading surface tracer', &
                            ' h2o_ice instead of h2o_vap'
     endif
@@ -462,8 +463,8 @@
     endif
     if (.not.found) then
-      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
-      write(*,*) "         ",trim(txt)," is set to zero"
+      if (is_master) write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+      if (is_master) write(*,*) "         ",trim(txt)," is set to zero"
     else
-      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
+      if (is_master) write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
                  minval(qsurf(:,iq)), maxval(qsurf(:,iq))
     endif
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/physiq.F90	(revision 298)
@@ -569,5 +569,5 @@
 !ym               rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
 
-         write(*,*) "physiq: firstcall, call phyetat0_academic"
+         if (is_master) write(*,*) "physiq: firstcall, call phyetat0_academic"
          call phyetat0_academic(ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,   &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
@@ -576,18 +576,18 @@
 
 !mi initialising tsurf with pt
-         write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
+         if (is_master) write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
          tsurf(:)=pt(:,1)
 
          if (pday.ne.day_ini) then
-           write(*,*) "ERROR in physiq.F90:"
-           write(*,*) "bad synchronization between physics and dynamics"
-           write(*,*) "dynamics day: ",pday
-           write(*,*) "physics day:  ",day_ini
-           write(*,*) "taking dynamics day for physics:  ",day_ini
+           if (is_master) write(*,*) "ERROR in physiq.F90:"
+           if (is_master) write(*,*) "bad synchronization between physics and dynamics"
+           if (is_master) write(*,*) "dynamics day: ",pday
+           if (is_master) write(*,*) "physics day:  ",day_ini
+           if (is_master) write(*,*) "taking dynamics day for physics:  ",day_ini
            day_ini=pday
 !ym           stop
          endif
 
-         write (*,*) 'In physiq day_ini =', day_ini
+         if (is_master) write (*,*) 'In physiq day_ini =', day_ini
 
 !        Initialize albedo and orbital calculation
@@ -598,5 +598,5 @@
          albedo(:)=albedo0(:)
 
-         if(tlocked)then
+         if(tlocked .and. is_master)then
             print*,'Planet is tidally locked at resonance n=',nres
             print*,'Make sure you have the right rotation rate!!!'
@@ -652,8 +652,8 @@
 
          else
-            print*,'WARNING! Thermal conduction in the soil turned off'
+            if (is_master) print*,'WARNING! Thermal conduction in the soil turned off'
             capcal(:)=1.e6
             fluxgrd(:)=intheat
-            print*,'Flux from ground = ',intheat,' W m^-2'
+            if (is_master) print*,'Flux from ground = ',intheat,' W m^-2'
          endif
          icount=1
@@ -667,6 +667,6 @@
                      status="old",iostat=ierr)
             if (ierr.ne.0) then
-              write(*,*) "physiq: Error! No num_run file!"
-              write(*,*) " (which is needed for sourceevol option)"
+              if (is_master) write(*,*) "physiq: Error! No num_run file!"
+              if (is_master) write(*,*) " (which is needed for sourceevol option)"
               stop
             endif
@@ -678,5 +678,5 @@
             if(num_run.ne.0.and.mod(num_run,2).eq.0)then
             !if(num_run.ne.0.and.mod(num_run,3).eq.0)then
-               print*,'Updating ice at end of this year!'
+               if (is_master) print*,'Updating ice at end of this year!'
                ice_update=.true.
                ice_min(:)=1.e4
@@ -696,6 +696,6 @@
            enddo
 
-           print*,'WARNING! Surface type currently decided by surface inertia'
-           print*,'This should be improved e.g. in newstart.F'
+           if (is_master) print*,'WARNING! Surface type currently decided by surface inertia'
+           if (is_master) print*,'This should be improved e.g. in newstart.F'
          endif!(.not.ok_slab_ocean)
 
@@ -711,5 +711,5 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if(nearco2cond) then
-            write(*,*)' WARNING! Starting at Tcond+1K'
+            if (is_master) write(*,*)' WARNING! Starting at Tcond+1K'
             do l=1, nlayer
                do ig=1,ngrid
@@ -892,5 +892,5 @@
 
            if(rings_shadow) then 
-                write(*,*) 'Rings shadow activated'
+                if (is_master) write(*,*) 'Rings shadow activated'
                 if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation 
                     pas = 1./nb_hours
@@ -1920,5 +1920,5 @@
 
 
-           print*,'--> Ls =',zls*180./pi
+           if (is_master) print*,'--> Ls =',zls*180./pi
 !        -------------------------------------------------------------------
 !        Writing NetCDF file  "RESTARTFI" at the end of the run
@@ -2359,5 +2359,5 @@
       endif
 
-      write(*,*) "physiq: done, zday=",zday
+      if (is_master) write(*,*) "physiq: done, zday=",zday
       return
     end subroutine physiq
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radii_mod.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radii_mod.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radii_mod.F90	(revision 298)
@@ -37,4 +37,5 @@
       use radinc_h, only: naerkind
       use aerosol_mod
+      use mod_phys_lmdz_para, only : is_master
 !      USE tracer_h
       Implicit none
@@ -54,5 +55,5 @@
       integer :: iaer   
       
-      print*,'enter su_aer_radii'
+      if (is_master) print*,'enter su_aer_radii'
           do iaer=1,naerkind
 !     these values will change once the microphysics gets to work
@@ -100,28 +101,28 @@
          if (radfixed) then
 
-	    write(*,*)"radius of H2O water particles:"
+	    if (is_master) write(*,*)"radius of H2O water particles:"
             rad_h2o=13. ! default value
             call getin_p("rad_h2o",rad_h2o)
-            write(*,*)" rad_h2o = ",rad_h2o
-
-	    write(*,*)"radius of H2O ice particles:"
+            if (is_master) write(*,*)" rad_h2o = ",rad_h2o
+
+	    if (is_master) write(*,*)"radius of H2O ice particles:"
             rad_h2o_ice=35. ! default value
             call getin_p("rad_h2o_ice",rad_h2o_ice)
-            write(*,*)" rad_h2o_ice = ",rad_h2o_ice
+            if (is_master) write(*,*)" rad_h2o_ice = ",rad_h2o_ice
 
 	 else
 
-            write(*,*)"Number mixing ratio of H2O water particles:"
+            if (is_master) write(*,*)"Number mixing ratio of H2O water particles:"
             Nmix_h2o=1.e6 ! default value
             call getin_p("Nmix_h2o",Nmix_h2o)
-            write(*,*)" Nmix_h2o = ",Nmix_h2o
-
-            write(*,*)"Number mixing ratio of H2O ice particles:"
+            if (is_master) write(*,*)" Nmix_h2o = ",Nmix_h2o
+
+            if (is_master) write(*,*)"Number mixing ratio of H2O ice particles:"
             Nmix_h2o_ice=Nmix_h2o ! default value
             call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
-            write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
+            if (is_master) write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
 	 endif
 
-      print*,'exit su_aer_radii'
+      if (is_master) print*,'exit su_aer_radii'
 
    end subroutine su_aer_radii
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspi.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspi.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspi.F90	(revision 298)
@@ -25,4 +25,5 @@
       use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma
       use datafile_mod, only: datadir
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -105,6 +106,6 @@
       close(131)
 
-      write(*,*) 'setspi: L_NSPECTI = ',L_NSPECTI, 'in the model '
-      write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
+      if (is_master) write(*,*) 'setspi: L_NSPECTI = ',L_NSPECTI, 'in the model '
+      if (is_master) write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
       if(nb.ne.L_NSPECTI) then
          write(*,*) 'MISMATCH !! I stop here'
@@ -125,4 +126,5 @@
 !$OMP BARRIER
 
+      if (is_master) then
       print*,''
       print*,'setspi: IR band limits:'
@@ -130,4 +132,5 @@
          print*,m,'-->',BWNI(M),' cm^-1'
       end do
+      end if
 
 !     Set up mean wavenumbers and wavenumber deltas.  Units of 
@@ -149,7 +152,9 @@
 !     original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
 
+      if (is_master) then
       print*,''
       print*,'setspi: Current Planck integration range:'
       print*,'T = ',dble(NTstar)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
+      end if
 
       do NW=1,L_NSPECTI
@@ -173,5 +178,5 @@
       ! force planck=sigma*eps*T^4 for each temperature in array
       if(forceEC)then
-         print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
+         if (is_master) print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
          do nt=NTstar,NTstop
             plancksum=0.0D0
@@ -198,7 +203,7 @@
             plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
          end do
-         print*,'setspi: At lower limit:'
-         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
-         print*,'actual sig*T^4   = ',sigma*(dble(nt)/NTfac)**4,' W m^-2'
+         if (is_master) print*,'setspi: At lower limit:'
+         if (is_master) print*,'in model sig*T^4 = ',plancksum,' W m^-2'
+         if (is_master) print*,'actual sig*T^4   = ',sigma*(dble(nt)/NTfac)**4,' W m^-2'
          
          ! check energy conservation at upper temperature boundary
@@ -208,8 +213,8 @@
             plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
          end do
-         print*,'setspi: At upper limit:'
-         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
-         print*,'actual sig*T^4   = ',sigma*(dble(nt)/NTfac)**4,' W m^-2'
-         print*,''
+         if (is_master) print*,'setspi: At upper limit:'
+         if (is_master) print*,'in model sig*T^4 = ',plancksum,' W m^-2'
+         if (is_master) print*,'actual sig*T^4   = ',sigma*(dble(nt)/NTfac)**4,' W m^-2'
+         if (is_master) print*,''
       endif
 
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspv.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspv.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/setspv.F90	(revision 298)
@@ -27,4 +27,5 @@
                              STELLARF,TAURAY
       use datafile_mod, only: datadir
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -82,6 +83,6 @@
       close(131)
 
-      write(*,*) 'setspv: L_NSPECTV = ',L_NSPECTV, 'in the model '
-      write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
+      if (is_master) write(*,*) 'setspv: L_NSPECTV = ',L_NSPECTV, 'in the model '
+      if (is_master) write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
       if(nb.ne.L_NSPECTV) then
          write(*,*) 'MISMATCH !! I stop here'
@@ -102,4 +103,5 @@
 !$OMP BARRIER
 
+      if (is_master) then
       print*,'setspv: VI band limits:'
       do M=1,L_NSPECTV+1
@@ -107,4 +109,5 @@
       end do
       print*,' '
+      end if
 
 !     Set up mean wavenumbers and wavenumber deltas.  Units of 
@@ -123,5 +126,5 @@
 !     Set up stellar spectrum
 
-      write(*,*)'setspv: Interpolating stellar spectrum from the hires data...'
+      if (is_master) write(*,*)'setspv: Interpolating stellar spectrum from the hires data...'
       call ave_stelspec(STELLAR)
 
@@ -132,6 +135,6 @@
          sum         = sum+STELLARF(N)
       end do
-      write(6,'("setspv: Stellar flux at 1 AU = ",f7.2," W m-2")') sum
-      print*,' '
+      if (is_master) write(6,'("setspv: Stellar flux at 1 AU = ",f7.2," W m-2")') sum
+      if (is_master) print*,' '
 
 
@@ -144,5 +147,5 @@
          call calc_rayleigh
       else
-         print*,'setspv: No Rayleigh scattering, check for NaN in output!'
+         if (is_master) print*,'setspv: No Rayleigh scattering, check for NaN in output!'
          do N=1,L_NSPECTV
             TAURAY(N) = 1E-16
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/suaer_corrk.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/suaer_corrk.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/suaer_corrk.F90	(revision 298)
@@ -11,4 +11,5 @@
       use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
       use aerosol_mod
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -170,5 +171,5 @@
 
        if (iaer.eq.iaero_back2lay) then
-         print*, 'naerkind= back2lay', iaer
+         if (is_master) print*, 'naerkind= back2lay', iaer
 
 !     visible
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/sugas_corrk.F90
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/sugas_corrk.F90	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/sugas_corrk.F90	(revision 298)
@@ -30,4 +30,5 @@
 !      use ioipsl_getincom
       use ioipsl_getincom_p 
+      use mod_phys_lmdz_para, only : is_master
       implicit none
 
@@ -105,5 +106,5 @@
       do igas=1,ngas
          read(111,*) gastype(igas)
-         print*,'Gas ',igas,' is ',gastype(igas)
+         if (is_master) print*,'Gas ',igas,' is ',gastype(igas)
       enddo
 
@@ -135,4 +136,5 @@
 
       ! display the values
+      if (is_master) then 
       print*,'Variable gas volume mixing ratios:'
       do n=1,L_REFVAR
@@ -141,4 +143,5 @@
       end do
       print*,''
+      end if
 
 !=======================================================================
@@ -172,4 +175,5 @@
  
       ! display the values
+      if (is_master) then
       print*,'Correlated-k g-space grid:'
       do n=1,L_NGAUSS
@@ -177,4 +181,5 @@
       end do
       print*,''
+      end if
 
 !=======================================================================
@@ -210,4 +215,5 @@
 
       ! display the values
+      if (is_master) then
       print*,'Correlated-k pressure grid (mBar):'
       do n=1,L_NPREF
@@ -215,4 +221,5 @@
       end do
       print*,''
+      end if
 
       ! save the min / max matrix values
@@ -255,8 +262,10 @@
 
       ! display the values
+      if (is_master) then
       print*,'Correlated-k temperature grid:'
       do n=1,L_NTREF
          print*,n,'.',tgasref(n),' K'
       end do
+      end if
 
       ! save the min / max matrix values
@@ -275,7 +284,7 @@
 
       ! display the values
-      print*,''
-      print*,'Correlated-k matrix size:' 
-      print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']' 
+      if (is_master) print*,''
+      if (is_master) print*,'Correlated-k matrix size:' 
+      if (is_master) print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']' 
 
 !=======================================================================
@@ -285,5 +294,5 @@
 !        wavelength used to separate IR from VI in graybody. We will need that anyway
          IR_VI_wnlimit=3000.
-         write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
+         if (is_master) write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
 	 
 	 nVI_limit=0
@@ -304,18 +313,18 @@
       if (graybody) then
 !        constant absorption coefficient in visible
-         write(*,*)"graybody: constant absorption coefficient in visible:"
+         if (is_master) write(*,*)"graybody: constant absorption coefficient in visible:"
          kappa_VI=-100000.
          call getin_p("kappa_VI",kappa_VI)
-         write(*,*)" kappa_VI = ",kappa_VI
+         if (is_master) write(*,*)" kappa_VI = ",kappa_VI
 	 kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule         
       
 !        constant absorption coefficient in IR
-         write(*,*)"graybody: constant absorption coefficient in InfraRed:"
+         if (is_master) write(*,*)"graybody: constant absorption coefficient in InfraRed:"
          kappa_IR=-100000.
          call getin_p("kappa_IR",kappa_IR)
-         write(*,*)" kappa_IR = ",kappa_IR	 
+         if (is_master) write(*,*)" kappa_IR = ",kappa_IR	 
          kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27	 ! conversion from m^2/kg to cm^2/molecule 
 
-         write(*,*)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
+         if (is_master) write(*,*)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
 	       
       Else
@@ -330,6 +339,6 @@
          if ((corrkdir(1:4).eq.'null'))then   !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
-            print*,'using no corrk data'
-            print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
+            if (is_master) print*,'using no corrk data'
+            if (is_master) print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
 	 else
             file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
@@ -363,5 +372,5 @@
 	 end if
       else
-         print*,'Visible corrk gaseous absorption is set to zero.'
+         if (is_master) print*,'Visible corrk gaseous absorption is set to zero.'
          gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
       endif
@@ -371,5 +380,5 @@
       ! INFRA-RED
       if ((corrkdir(1:4).eq.'null'))then       !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
-         print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
+         if (is_master) print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
 !$OMP MASTER         
          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
@@ -662,4 +671,5 @@
       endif
 
+      if (is_master) then
       print*,'----------------------------------------------------'
       print*,'And that`s all we have. It`s possible that other'
@@ -667,4 +677,5 @@
       print*,'don`t yet have data for it...'
       print*,''
+      end if
 
 !     Deallocate local arrays
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/surfini.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/surfini.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/surfini.F	(revision 298)
@@ -5,4 +5,5 @@
       use comgeomfi_h, only: lati
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
+      use mod_phys_lmdz_para, only : is_master
 
       IMPLICIT NONE
@@ -36,6 +37,8 @@
       call planetwide_minval(albedodat,min_albedo)
       call planetwide_maxval(albedodat,max_albedo)
+      if (is_master) then
       write(*,*) 'surfini: minimum corrected albedo',min_albedo
       write(*,*) 'surfini: maximum corrected albedo',max_albedo
+      end if
 
       if (igcm_co2_ice.ne.0) then
@@ -52,12 +55,16 @@
         ENDDO ! of DO ig=1,ngrid     
       else
+        if (is_master) then
         write(*,*) "surfini: No CO2 ice tracer on surface  ..."
         write(*,*) "         and therefore no albedo change."
+        end if
       endif
 
       call planetwide_minval(psolaralb,min_albedo)
       call planetwide_maxval(psolaralb,max_albedo)
+      if (is_master) then
       write(*,*) 'surfini: minimum corrected albedo',min_albedo
       write(*,*) 'surfini: maximum corrected albedo',max_albedo
+      end if
 
       END
Index: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F
===================================================================
--- codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F	(revision 297)
+++ codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/tabfi.F	(revision 298)
@@ -53,4 +53,5 @@
       use planete_mod, only: year_day, periastr, apoastr, peri_day,
      &                       obliquit, z0, lmixmin, emin_turb
+      use mod_phys_lmdz_para, only : is_master
 
       implicit none
@@ -82,5 +83,6 @@
       LOGICAL :: found
       
-      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
+      if (is_master) write(*,*)"tabfi: nid=",nid," tab0=",tab0,
+     &		" Lmodif=",Lmodif
 
       IF (nid.eq.0) then
@@ -90,5 +92,5 @@
         ! Ehouarn: Default Saturn values:
         tab_cntrl(:)=0
-        write(*,*) "Using default Saturn values..."
+        if (is_master) write(*,*) "Using default Saturn values..."
         ! these should be read in a def file I guess...
         lmax=0 ! not used anyways
@@ -164,5 +166,5 @@
          call abort
        else
-         write(*,*)'tabfi: tab_cntrl',tab_cntrl
+         if (is_master) write(*,*)'tabfi: tab_cntrl',tab_cntrl
        endif
 c
@@ -227,4 +229,5 @@
    5  FORMAT(a20,f12.2,f12.2)
  
+      if (is_master) then
       write(*,*) '*****************************************************'
       write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
@@ -265,4 +268,5 @@
       write(*,*)
       write(*,*) 'Lmodif in tabfi!!!!!!!',Lmodif
+      end if !of if (is_master)
 
 c-----------------------------------------------------------------------
