Index: trunk/LMDZ.GENERIC/libf/phygeneric/aerosol_optical_properties_averaging.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/aerosol_optical_properties_averaging.F	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/aerosol_optical_properties_averaging.F	(revision 4081)
@@ -4,4 +4,5 @@
      & epir,omegir,gir,qref,omegaref)
 
+      USE rad_blackbody_mod, ONLY: rad_blackbody_planck_law_wavelength
 
       IMPLICIT NONE
Index: trunk/LMDZ.GENERIC/libf/phygeneric/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/physiq_mod.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/physiq_mod.F90	(revision 4081)
@@ -25,4 +25,5 @@
       use rad_correlatedk_ini_aerosol_mod, only: rad_correlatedk_ini_aerosol
       use rad_correlatedk_init_stellar_mod, only: rad_correlatedk_init_stellar
+      use rad_correlatedk_init_thermal_mod, only: rad_correlatedk_init_thermal
       use aerosol_radius, only: aerosol_radius_h2o_liquid_ice_mixture, aerosol_radius_co2
       use aerosol_global_variables , only: aerosol_init, iaero_co2, iaero_h2o
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_blackbody.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_blackbody.F	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_blackbody.F	(revision 4081)
@@ -1,43 +1,59 @@
+      module rad_blackbody_mod
+      
+      implicit none
+      
+      contains
+      
       subroutine rad_blackbody_planck_law_wavelength(blalong,blat,blae)
 
-      implicit double precision (a-h,o-z)
+      implicit none
+      
+      double precision, intent(in) :: blalong
+      double precision, intent(in) :: blat
+      double precision, intent(out) :: blae
+      
+!      double precision :: sigma,pi,c0,h,cbol,rind,c,c1,c2
 
       ! physical constants
-      sigma=5.670374D-8
-      pi=datan(1.d0)*4.d0
-      c0=2.997925d+08
-      h=6.62607d-34
-      cbol=1.380649d-23
-      rind=1.d0
-      c=c0/rind
-      c1=h*(c**2)
-      c2=h*c/cbol
+      double precision, parameter :: sigma=5.670374D-8
+      double precision, parameter :: pi=atan(1.d0)*4.d0
+      double precision, parameter :: c0=2.997925d+08
+      double precision, parameter :: h=6.62607d-34
+      double precision, parameter :: cbol=1.380649d-23
+      double precision, parameter :: rind=1.d0
+      double precision, parameter :: c=c0/rind
+      double precision, parameter :: c1=h*(c**2)
+      double precision, parameter :: c2=h*c/cbol
 
 
-      blae=2.d0*pi*c1/blalong**5/(dexp(c2/blalong/blat)-1.d0)
+      blae=2.d0*pi*c1/blalong**5/(exp(c2/blalong/blat)-1.d0)
 
-
-      return
-      end
+      end subroutine rad_blackbody_planck_law_wavelength
 
       subroutine rad_blackbody_planck_law_wavenumber(blalong,blat,blae)
 
-      implicit double precision (a-h,o-z)
+      implicit none
+
+      double precision, intent(in) :: blalong
+      double precision, intent(in) :: blat
+      double precision, intent(out) :: blae
+      
+!      double precision :: sigma,pi,c0,h,cbol,rind,c,c1,c2
 
       ! physical constants
-      sigma=5.670374D-8
-      pi=datan(1.d0)*4.d0
-      c0=2.997925d+08
-      h=6.62607d-34
-      cbol=1.380649d-23
-      rind=1.d0
-      c=c0/rind
-      c1=h*(c**2)
-      c2=h*c/cbol
+      double precision, parameter :: sigma=5.670374D-8
+      double precision, parameter :: pi=atan(1.d0)*4.d0
+      double precision, parameter :: c0=2.997925d+08
+      double precision, parameter :: h=6.62607d-34
+      double precision, parameter :: cbol=1.380649d-23
+      double precision, parameter :: rind=1.d0
+      double precision, parameter :: c=c0/rind
+      double precision, parameter :: c1=h*(c**2)
+      double precision, parameter :: c2=h*c/cbol
 
 
-      blae=2.d0*pi*c1*blalong**3/(dexp(c2*blalong/blat)-1.d0)
+      blae=2.d0*pi*c1*blalong**3/(exp(c2*blalong/blat)-1.d0)
 
-
-      return
-      end
+      end subroutine rad_blackbody_planck_law_wavenumber
+      
+      end module rad_blackbody_mod
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_stellar.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_stellar.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_stellar.F90	(revision 4081)
@@ -34,4 +34,5 @@
       use callkeys_mod, only: Fat1AU
       use rad_correlatedk_stellar_spectrum_mod, only: rad_correlatedk_stellar_spectrum 
+      use mod_phys_lmdz_para, only : is_master, bcast
 
       implicit none
@@ -74,41 +75,43 @@
       endif
 	
-!$OMP MASTER        
-      nb=0
-      ierr=0
-      ! check that the file contains the right number of bands 
-      open(131,file=file_path,form='formatted')
-      read(131,*,iostat=ierr) file_entries
-      do while (ierr==0)
+      if (is_master) then ! only the master needs to read in BWNV(:) from file 
+       nb=0
+       ierr=0
+       ! check that the file contains the right number of bands 
+       open(131,file=file_path,form='formatted')
+       read(131,*,iostat=ierr) file_entries
+       do while (ierr==0)
         read(131,*,iostat=ierr) dummy
         if (ierr==0) nb=nb+1
-      enddo
-      close(131)
+       enddo
+       close(131)
 
-      write(*,*) 'rad_correlatedk_init_stellar: L_NSPECTV = ',L_NSPECTV, 'in the model '
-      write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
-      if(nb.ne.L_NSPECTV) then
+       write(*,*) 'rad_correlatedk_init_stellar: L_NSPECTV = ',L_NSPECTV, 'in the model '
+       write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
+       if(nb.ne.L_NSPECTV) then
          write(*,*) 'MISMATCH !! I stop here'
          call abort_physic("rad_correlatedk_init_stellar","The number of entries in narrowbands_VI.in does not match with L_NSPECTV",1)
-      endif
+       endif
 
-      ! load and display the data
-      open(111,file=file_path,form='formatted')
-      read(111,*) 
-       do M=1,L_NSPECTV-1
+       ! load and display the data
+       open(111,file=file_path,form='formatted')
+       read(111,*) 
+        do M=1,L_NSPECTV-1
          read(111,*) BWNV(M)
-      end do
-      read(111,*) lastband
-      close(111)
-      BWNV(L_NSPECTV)  =lastband(1)
-      BWNV(L_NSPECTV+1)=lastband(2)
-!$OMP END MASTER
-!$OMP BARRIER
+       end do
+       read(111,*) lastband
+       close(111)
+       BWNV(L_NSPECTV)  =lastband(1)
+       BWNV(L_NSPECTV+1)=lastband(2)
 
-      print*,'rad_correlatedk_init_stellar: VI band limits:'
-      do M=1,L_NSPECTV+1
+       print*,'rad_correlatedk_init_stellar: VI band limits:'
+       do M=1,L_NSPECTV+1
          print*,m,'-->',BWNV(M),' cm^-1'
-      end do
-      print*,' '
+       end do
+       print*,' '
+      endif ! of if (is_master)
+
+      ! Broadcast BWNV to all cores
+      call bcast(BWNV)
 
 !     Set up mean wavenumbers and wavenumber deltas.  Units of 
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_thermal.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_thermal.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_init_thermal.F90	(revision 4081)
@@ -1,2 +1,8 @@
+module rad_correlatedk_init_thermal_mod
+      
+implicit none
+      
+contains
+      
       subroutine rad_correlatedk_init_thermal
 
@@ -26,4 +32,5 @@
       use datafile_mod, only: datadir
       use comcstfi_mod, only: pi
+      use mod_phys_lmdz_para, only : is_master, bcast
 
       implicit none
@@ -90,42 +97,44 @@
       endif
     
-!$OMP MASTER    
-      nb=0
-      ierr=0
-      ! check that the file contains the right number of bands 
-      open(131,file=file_path,form='formatted')
-      read(131,*,iostat=ierr) file_entries
-      do while (ierr==0)
+      if (is_master) then ! only the master needs to read in BWNI(:) from file  
+       nb=0
+       ierr=0
+       ! check that the file contains the right number of bands 
+       open(131,file=file_path,form='formatted')
+       read(131,*,iostat=ierr) file_entries
+       do while (ierr==0)
         read(131,*,iostat=ierr) dummy
 !        write(*,*) 'rad_correlatedk_init_thermal: file_entries:',dummy,'ierr=',ierr
         if (ierr==0) nb=nb+1
-      enddo
-      close(131)
-
-      write(*,*) 'rad_correlatedk_init_thermal: L_NSPECTI = ',L_NSPECTI, 'in the model '
-      write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
-      if(nb.ne.L_NSPECTI) then
+       enddo
+       close(131)
+
+       write(*,*) 'rad_correlatedk_init_thermal: L_NSPECTI = ',L_NSPECTI, 'in the model '
+       write(*,*) '        there are   ',nb, 'entries in ',TRIM(file_path)
+       if(nb.ne.L_NSPECTI) then
          write(*,*) 'MISMATCH !! I stop here'
          call abort_physic("rad_correlatedk_init_thermal","The number of entries in narrowbands_IR.in does not match with L_NSPECTI",1)
-      endif
-
-      ! load and display the data
-      open(111,file=file_path,form='formatted')
-      read(111,*) 
-      do M=1,L_NSPECTI-1
+       endif
+
+       ! load and display the data
+       open(111,file=file_path,form='formatted')
+       read(111,*) 
+       do M=1,L_NSPECTI-1
          read(111,*) BWNI(M)
-      end do
-      read(111,*) lastband
-      close(111)
-      BWNI(L_NSPECTI)  =lastband(1)
-      BWNI(L_NSPECTI+1)=lastband(2)
-!$OMP END MASTER
-!$OMP BARRIER
-
-      print*,''
-      print*,'rad_correlatedk_init_thermal: IR band limits:'
-      do M=1,L_NSPECTI+1
+       end do
+       read(111,*) lastband
+       close(111)
+       BWNI(L_NSPECTI)  =lastband(1)
+       BWNI(L_NSPECTI+1)=lastband(2)
+
+       print*,''
+       print*,'rad_correlatedk_init_thermal: IR band limits:'
+       do M=1,L_NSPECTI+1
          print*,m,'-->',BWNI(M),' cm^-1'
-      end do
+       end do
+      endif ! of if (is_master)
+      
+      ! Broadcast BWNI to all cores
+      call bcast(BWNI)
 
 !     Set up mean wavenumbers and wavenumber deltas.  Units of 
@@ -151,5 +160,5 @@
       print*,'T = ',dble(NTstart)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
 
-      IF(.NOT.ALLOCATED(planckir)) ALLOCATE(planckir(L_NSPECTI,NTstop-NTstart+1))
+      ALLOCATE(planckir(L_NSPECTI,NTstop-NTstart+1))
 
       do NW=1,L_NSPECTI
@@ -223,4 +232,5 @@
       endif
 
-      return
     end subroutine rad_correlatedk_init_thermal
+    
+end module rad_correlatedk_init_thermal_mod
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_online_recombination.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_online_recombination.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_online_recombination.F90	(revision 4081)
@@ -38,7 +38,10 @@
     ! Following arrays are allocated in rad_correlatedk_recombination_setup (excepted otf2tot_idx, in rad_correlatedk_recombination_init) and deallocated in callcork lastcall
     REAL*8, save,  DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_recomb, gasv_recomb
+  !$OMP THREADPRIVATE(gasi_recomb,gasv_recomb)
     REAL*8, save,  DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi_otf, gasv_otf
+  !$OMP THREADPRIVATE(gasi_otf,gasv_otf)
     REAL*8, save,  DIMENSION(:),         ALLOCATABLE :: w_cum
     REAL*8,  save, DIMENSION(:),         ALLOCATABLE :: wtwospec
+  !$OMP THREADPRIVATE(w_cum,wtwospec)
   
     INTEGER, save, DIMENSION(:),         ALLOCATABLE :: otf2tot_idx 
@@ -47,5 +50,5 @@
   
     INTEGER, save, DIMENSION(:),         ALLOCATABLE :: permut_idx
-  !$OMP THREADPRIVATE(gasi_recomb,gasv_recomb,w_cum,otf2tot_idx,rcb2tot_idx,otf2rcb_idx,permut_idx)
+  !$OMP THREADPRIVATE(otf2tot_idx,rcb2tot_idx,otf2rcb_idx,permut_idx)
   
   CONTAINS
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_rayleigh_scattering_opacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_rayleigh_scattering_opacity.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_rayleigh_scattering_opacity.F90	(revision 4081)
@@ -38,4 +38,5 @@
       use comcstfi_mod, only: g, pi
       use callkeys_mod, only: strictboundrayleigh
+      use rad_blackbody_mod, only: rad_blackbody_planck_law_wavelength
 
       implicit none
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_stellar_spectrum.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_stellar_spectrum.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_stellar_spectrum.F90	(revision 4081)
@@ -35,8 +35,10 @@
       use callkeys_mod, only: stelbbody,stelTbb
       use ioipsl_getin_p_mod, only: getin_p
+      use rad_blackbody_mod, only: rad_blackbody_planck_law_wavelength
+      use mod_phys_lmdz_para, only : is_master, bcast
 
       implicit none
 
-      real*8 STELLAR(L_NSPECTV)
+      real*8,intent(out) :: STELLAR(L_NSPECTV)
 
       integer Nfine
@@ -44,6 +46,6 @@
       integer ifine,band
 
-      real,allocatable,save :: lam(:),stel_f(:) ! read by master thread
-                                                ! but used by all threads
+      real,allocatable :: lam(:),stel_f(:) ! read by master 
+                                                
       real lamm,lamp
       real dl
@@ -61,5 +63,6 @@
       STELLAR(:)=0.0
 
-      print*,'enter ave_stellspec'
+      if (is_master) print*,'entering rad_correlatedk_stellar_spectrum'
+
       if(stelbbody)then
          tstellar=stelTbb
@@ -89,6 +92,8 @@
          end if
 	 
-         write(*,*) "Input stellar temperature is:"
-         write(*,*) "tstellar = ",tstellar
+         if (is_master) then
+           write(*,*) "Input stellar temperature is:"
+           write(*,*) "tstellar = ",tstellar
+         endif
 
          ! load high resolution stellar data
@@ -97,14 +102,15 @@
          call getin_p("stelspec_file",stelspec_file) ! default path
 	 
-         write(*,*) "Input stellar spectrum file is:"
-         write(*,*) "stelspec_file = ",trim(stelspec_file)
-         write(*,*) 'Please use ',1,' and only ',1,' header line in ',trim(stelspec_file)
+         if (is_master) then ! only the master needs to do this
+          write(*,*) "Input stellar spectrum file is:"
+          write(*,*) "stelspec_file = ",trim(stelspec_file)
+          write(*,*) 'Please use ',1,' and only ',1,' header line in ',trim(stelspec_file)
 
-         ! Check the target file is there
-         file_path = trim(datadir)//'/stellar_spectra/'//stelspec_file
-         print*, 'stellar flux : ', file_path
-         inquire(FILE=file_path,EXIST=file_exists)         
+          ! Check the target file is there
+          file_path = trim(datadir)//'/stellar_spectra/'//stelspec_file
+          print*, 'stellar flux : ', file_path
+          inquire(FILE=file_path,EXIST=file_exists)         
    
-         if (.not.file_exists) THEN
+          if (.not.file_exists) THEN
 	   write(*,*)'Beware that startype is now deprecated, you should use '
 	   write(*,*)'stelspec_file and tstellar to define the input stellar spectrum.'
@@ -122,29 +128,37 @@
 	   write(*,*)'https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/stellar_spectra/'
            call abort_physic("rad_correlatedk_stellar_spectrum", "Unable to read stellar flux file", 1)
-         end if
+          end if
 
-!$OMP MASTER
-         ! Open the file
-         OPEN(UNIT=110,FILE=file_path,STATUS='old',iostat=ios)
-         ! Get number of line in the file
-         READ(110,*) ! skip first line header just in case
-         Nfine = 0
-         do
+          ! Open the file
+          OPEN(UNIT=110,FILE=file_path,STATUS='old',iostat=ios)
+          ! Get number of line in the file
+          READ(110,*) ! skip first line header just in case
+          Nfine = 0
+          do
            read(110,*,iostat=ios)
            if (ios<0) exit
            Nfine = Nfine + 1
-         end do
-         rewind(110) ! Rewind file after counting lines
-         READ(110,*) ! skip first line header just in case
+          end do
+          rewind(110) ! Rewind file after counting lines
+          READ(110,*) ! skip first line header just in case
 
-	 allocate(lam(Nfine),stel_f(Nfine))
+          ! allocate arrays
+	  allocate(lam(Nfine),stel_f(Nfine))
 
-         do ifine=1,Nfine
+          ! read arrays
+          do ifine=1,Nfine
            read(110,*) lam(ifine), stel_f(ifine) ! lam [um] stel_f [per unit of wavelength] (integrated and normalized by Fat1AU)
-         enddo
+          enddo
 
-!$OMP END MASTER
-!$OMP BARRIER
+         endif ! of if(is_master)
+
+         ! brodcast Nfine, allocate arrays for all except master
+         call bcast(Nfine)
+	 if (.not.is_master) allocate(lam(Nfine),stel_f(Nfine))
+         ! broadcast arrays
+         call bcast(lam)
+         call bcast(stel_f)
 	 
+
          ! sum data by band
          band=1
@@ -166,10 +180,9 @@
 	 
          STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
-!$OMP BARRIER
-!$OMP MASTER
-	 deallocate(lam)
+
+	 ! cleanup: deallocate arrays
+         deallocate(lam)
 	 deallocate(stel_f)
-!$OMP END MASTER
-!$OMP BARRIER         
+
       endif !stelbbody
 
Index: trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90	(revision 4079)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/radcommon_h.F90	(revision 4081)
@@ -63,10 +63,13 @@
 !
 
-      REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in rad_correlatedk_init_thermal
-      REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in rad_correlatedk_init_stellar
+      REAL*8 BWNI(L_NSPECTI+1) !BWNI read by master in rad_correlatedk_init_thermal
+!$OMP THREADPRIVATE(BWNI)
+      REAL*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) 
+!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI)
+      REAL*8 BWNV(L_NSPECTV+1) !BWNV read by master in rad_correlatedk_init_stellar
+!$OMP THREADPRIVATE(BWNV)
+      REAL*8 WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
       REAL*8 STELLARF(L_NSPECTV)
-!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
-	!$OMP WNOV,DWNV,WAVEV,&
-	!$OMP STELLARF)
+!$OMP THREADPRIVATE(WNOV,DWNV,WAVEV,STELLARF)
 
       REAL*8 blami(L_NSPECTI+1)