Commit d9f0585f authored by Sabine's avatar Sabine
Browse files

Merge branch 'dev' of git.nilu.no:flexpart/flexpart into dev

parents d404d981 c8fc7249
......@@ -66,9 +66,10 @@ program flexpart
end do
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
flexversion='Version '//trim(flexversion_major)//'.0beta (2015-05-01)'
flexversion='Version '//trim(flexversion_major)//'.1beta (2016-11-02)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -383,6 +384,17 @@ program flexpart
end do
end do
! Inform whether output kernel is used or not
!*********************************************
if (lroot) then
if (lnokernel) then
write(*,*) "Concentrations are calculated without using kernel"
else
write(*,*) "Concentrations are calculated using kernel"
end if
end if
! Calculate particle trajectories
!********************************
......@@ -401,10 +413,15 @@ program flexpart
call timemanager
! NIK 16.02.2005
do i=1,nspec
write(*,*) '**********************************************'
write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
write(*,*) 'Total number of occurences of in-cloud scavenging', tot_inc_count
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
end do
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
&XPART MODEL RUN!'
......
......@@ -54,15 +54,17 @@ program flexpart
character(len=256) :: inline_options !pathfile, flexversion, arg2
! Initialize mpi
!*********************
! Initialize mpi
!*********************
call mpif_init
if (mp_measure_time) call mpif_mtime('flexpart',0)
! Initialize arrays in com_mod
!*****************************
call com_mod_allocate_part(maxpart_mpi)
! Initialize arrays in com_mod
!*****************************
if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
! Generate a large number of random numbers
!******************************************
......@@ -78,7 +80,7 @@ program flexpart
! FLEXPART version string
flexversion_major = '10' ! Major version number, also used for species file names
! flexversion='Ver. 10 Beta MPI (2015-05-01)'
flexversion='Ver. '//trim(flexversion_major)//' Beta MPI (2015-05-01)'
flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
verbosity=0
! Read the pathnames where input/output files are stored
......@@ -305,9 +307,11 @@ program flexpart
print*,'Initialize all particles to non-existent'
endif
if (.not.(lmpreader.and.lmp_use_reader)) then
do j=1, size(itra1) ! maxpart_mpi
itra1(j)=-999999999
end do
end if
! For continuation of previous run, read in particle positions
!*************************************************************
......@@ -317,7 +321,7 @@ program flexpart
print*,'call readpartpositions'
endif
! readwind process skips this step
if (lmp_use_reader.and..not.lmpreader) call readpartpositions
if (.not.(lmpreader.and.lmp_use_reader)) call readpartpositions
else
if (verbosity.gt.0 .and. lroot) then
print*,'numpart=0, numparticlecount=0'
......@@ -424,6 +428,16 @@ program flexpart
end do
end do
! Inform whether output kernel is used or not
!*********************************************
if (lroot) then
if (lnokernel) then
write(*,*) "Concentrations are calculated without using kernel"
else
write(*,*) "Concentrations are calculated using kernel"
end if
end if
! Calculate particle trajectories
!********************************
......@@ -447,24 +461,29 @@ program flexpart
! NIK 16.02.2005
if (lroot) then
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
else
if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
call MPI_Reduce(tot_blc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(tot_inc_count, 0, 1, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
end if
end if
if (lroot) then
do i=1,nspec
write(*,*) '**********************************************'
write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
write(*,*) 'Total number of occurences of in-cloud scavenging', tot_inc_count
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
write(*,*) '**********************************************'
end do
write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
&XPART MODEL RUN!'
......
This diff is collapsed.
......@@ -131,15 +131,16 @@ module com_mod
! gdomainfill .T., if domain-filling is global, .F. if not
!ZHG SEP 2015 wheather or not to read clouds from GRIB
logical :: readclouds
logical :: readclouds=.false.
!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
logical :: sumclouds
logical :: sumclouds=.false.
logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
!NIK 16.02.2015
integer(selected_int_kind(16)) :: tot_blc_count=0, tot_inc_count=0
integer(selected_int_kind(16)), dimension(maxspec) :: tot_blc_count=0, &
&tot_inc_count=0
!*********************************************************************
......@@ -575,7 +576,8 @@ module com_mod
integer :: numxgridn,numygridn
real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
!real outheight(maxzgrid),outheighthalf(maxzgrid)
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
& OHREA,ASSSPEC
! numxgrid,numygrid number of grid points in x,y-direction
! numxgridn,numygridn number of grid points in x,y-direction for nested output grid
......@@ -592,6 +594,7 @@ module com_mod
! DRYDEP .true., if dry deposition is switched on
! DRYDEPSPEC .true., if dry deposition is switched on for that species
! WETDEP .true., if wet deposition is switched on
! WETDEPSPEC .true., if wet deposition is switched on for that species
! OHREA .true., if OH reaction is switched on
! ASSSPEC .true., if there are two species asscoiated
! (i.e. transfer of mass between these two occurs
......@@ -742,7 +745,7 @@ module com_mod
logical, parameter :: nmlout=.true.
! These variables are used to avoid having separate versions of
! files in cases where differences with MPI version is minor (eso)
! files in cases where differences with MPI version are minor (eso)
!*****************************************************************
integer :: mpi_mode=0 ! .gt. 0 if running MPI version
logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
......
......@@ -125,7 +125,7 @@ subroutine conccalc(itime,weight)
! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
!*****************************************************************************
do ind=indz,indzp
rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,memind(2)) &
+p2*rho(ixp,jy ,ind,2) &
+p3*rho(ix ,jyp,ind,2) &
+p4*rho(ixp,jyp,ind,2)
......@@ -180,7 +180,7 @@ subroutine conccalc(itime,weight)
! If a particle is close to the domain boundary, do not use the kernel either.
!*****************************************************************************
if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgrid-1)-0.5).or. &
(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
......
......@@ -193,7 +193,7 @@ subroutine conccalc(itime,weight)
! If a particle is close to the domain boundary, do not use the kernel either.
!*****************************************************************************
if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgrid-1)-0.5).or. &
(yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
......
......@@ -625,24 +625,24 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Reinitialization of grid
!*************************
! do ks=1,nspec
! do kp=1,maxpointspec_act
! do i=1,numreceptor
! creceptor(i,ks)=0.
! end do
! do jy=0,numygrid-1
! do ix=0,numxgrid-1
! do l=1,nclassunc
! do nage=1,nageclass
! do kz=1,numzgrid
! gridunc(ix,jy,kz,ks,kp,l,nage)=0.
! end do
! end do
! end do
! end do
! end do
! end do
! end do
! do ks=1,nspec
! do kp=1,maxpointspec_act
! do i=1,numreceptor
! creceptor(i,ks)=0.
! end do
! do jy=0,numygrid-1
! do ix=0,numxgrid-1
! do l=1,nclassunc
! do nage=1,nageclass
! do kz=1,numzgrid
! gridunc(ix,jy,kz,ks,kp,l,nage)=0.
! end do
! end do
! end do
! end do
! end do
! end do
! end do
creceptor(:,:)=0.
gridunc(:,:,:,:,:,:,:)=0.
......
......@@ -104,9 +104,6 @@ subroutine concoutput_nest(itime,outnum)
! Measure execution time
if (mp_measure_time) call mpif_mtime('iotime',0)
! call cpu_time(mp_root_time_beg)
! mp_root_wtime_beg = mpi_wtime()
! end if
if (verbosity.eq.1) then
print*,'inside concoutput_surf '
......
......@@ -21,41 +21,41 @@
subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
drygridtotalunc)
! i i o o
! o
!*****************************************************************************
! *
! Output of the concentration grid and the receptor concentrations. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output *
! in sparse matrix format; additional output of *
! uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in file par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
!*****************************************************************************
! *
! Variables: *
! outnum number of samples *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
! *
!*****************************************************************************
! i i o o
! o
!*****************************************************************************
! *
! Output of the concentration grid and the receptor concentrations. *
! *
! Author: A. Stohl *
! *
! 24 May 1995 *
! *
! 13 April 1999, Major update: if output size is smaller, dump output *
! in sparse matrix format; additional output of *
! uncertainty *
! *
! 05 April 2000, Major update: output of age classes; output for backward*
! runs is time spent in grid cell times total mass of *
! species. *
! *
! 17 February 2002, Appropriate dimensions for backward and forward runs *
! are now specified in file par_mod *
! *
! June 2006, write grid in sparse matrix with a single write command *
! in order to save disk space *
! *
! 2008 new sparse matrix format *
! *
!*****************************************************************************
! *
! Variables: *
! outnum number of samples *
! ncells number of cells with non-zero concentrations *
! sparse .true. if in sparse matrix format, else .false. *
! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
! *
!*****************************************************************************
use unc_mod
use point_mod
......@@ -72,24 +72,24 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
real :: sp_fact
real :: outnum,densityoutrecept(maxreceptor),xl,yl
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
! + maxageclass)
!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
! + maxageclass)
!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
! + maxpointspec_act,maxageclass),
! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass),
! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
!real sparse_dump_r(numxgrid*numygrid*numzgrid)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
! + maxageclass)
!real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
! + maxageclass)
!real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
! + maxpointspec_act,maxageclass),
! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass),
! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
! + maxpointspec_act,maxageclass)
!real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
!real sparse_dump_r(numxgrid*numygrid*numzgrid)
!integer sparse_dump_i(numxgrid*numygrid*numzgrid)
!real sparse_dump_u(numxgrid*numygrid*numzgrid)
real(dep_prec) :: auxgrid(nclassunc)
real(sp) :: gridtotal,gridsigmatotal,gridtotalunc
real(dep_prec) :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
......@@ -108,23 +108,23 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Determine current calendar date, needed for the file name
!**********************************************************
! Determine current calendar date, needed for the file name
!**********************************************************
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,ihmmss)
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
!write(unitdates,'(a)') adate//atime
!write(unitdates,'(a)') adate//atime
open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
write(unitdates,'(a)') adate//atime
close(unitdates)
! For forward simulations, output fields have dimension MAXSPEC,
! for backward simulations, output fields have dimension MAXPOINT.
! Thus, make loops either about nspec, or about numpoint
!*****************************************************************
! For forward simulations, output fields have dimension MAXSPEC,
! for backward simulations, output fields have dimension MAXPOINT.
! Thus, make loops either about nspec, or about numpoint
!*****************************************************************
if (ldirect.eq.1) then
......@@ -148,12 +148,12 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
!*******************************************************************
! Compute air density: sufficiently accurate to take it
! from coarse grid at some time
! Determine center altitude of output layer, and interpolate density
! data to that altitude
!*******************************************************************
!*******************************************************************
! Compute air density: sufficiently accurate to take it
! from coarse grid at some time
! Determine center altitude of output layer, and interpolate density
! data to that altitude
!*******************************************************************
do kz=1,numzgrid
if (kz.eq.1) then
......@@ -192,7 +192,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
! Output is different for forward and backward simulations
! Output is different for forward and backward simulations
do kz=1,numzgrid
do jy=0,numygrid-1
do ix=0,numxgrid-1
......@@ -205,10 +205,10 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
end do
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
!*********************************************************************
! Determine the standard deviation of the mean concentration or mixing
! ratio (uncertainty of the output) and the dry and wet deposition
!*********************************************************************
if (verbosity.eq.1) then
print*,'concoutput_surf 3 (sd)'
......@@ -252,18 +252,18 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
do jy=0,numygrid-1
do ix=0,numxgrid-1
! WET DEPOSITION
! WET DEPOSITION
if ((WETDEP).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,wetgrid(ix,jy), &
wetgridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
! Multiply by number of classes to get total concentration
wetgrid(ix,jy)=wetgrid(ix,jy) &
*nclassunc
wetgridtotal=wetgridtotal+wetgrid(ix,jy)
! Calculate standard deviation of the mean
! Calculate standard deviation of the mean
wetgridsigma(ix,jy)= &
wetgridsigma(ix,jy)* &
sqrt(real(nclassunc))
......@@ -271,18 +271,18 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
wetgridsigma(ix,jy)
endif
! DRY DEPOSITION
! DRY DEPOSITION
if ((DRYDEP).and.(ldirect.gt.0)) then
do l=1,nclassunc
auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
end do
call mean(auxgrid,drygrid(ix,jy), &
drygridsigma(ix,jy),nclassunc)
! Multiply by number of classes to get total concentration
! Multiply by number of classes to get total concentration
drygrid(ix,jy)=drygrid(ix,jy)* &
nclassunc
drygridtotal=drygridtotal+drygrid(ix,jy)
! Calculate standard deviation of the mean
! Calculate standard deviation of the mean
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
......@@ -290,18 +290,18 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
drygridsigma(ix,jy)
endif
! CONCENTRATION OR MIXING RATIO
! CONCENTRATION OR MIXING RATIO
do kz=1,numzgrid
do l=1,nclassunc
auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
end do
call mean(auxgrid,grid(ix,jy,kz), &
gridsigma(ix,jy,kz),nclassunc)
! Multiply by number of classes to get total concentration
! Multiply by number of classes to get total concentration
grid(ix,jy,kz)= &
grid(ix,jy,kz)*nclassunc
gridtotal=gridtotal+grid(ix,jy,kz)
! Calculate standard deviation of the mean
! Calculate standard deviation of the mean
gridsigma(ix,jy,kz)= &
gridsigma(ix,jy,kz)* &
sqrt(real(nclassunc))
......@@ -312,13 +312,13 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
end do
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
!*******************************************************************
! Generate output: may be in concentration (ng/m3) or in mixing
! ratio (ppt) or both
! Output the position and the values alternated multiplied by
! 1 or -1, first line is number of values, number of positions
! For backward simulations, the unit is seconds, stored in grid_time
!*******************************************************************
if (verbosity.eq.1) then
print*,'concoutput_surf 4 (output)'
......@@ -326,8 +326,8 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Concentration output
!*********************
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
......@@ -337,7 +337,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Wet deposition
! Wet deposition
sp_count_i=0
sp_count_r=0
sp_fact=-1.
......@@ -345,7 +345,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
if ((ldirect.eq.1).and.(WETDEP)) then
do jy=0,numygrid-1
do ix=0,numxgrid-1
! concentraion greater zero
! concentraion greater zero
if (wetgrid(ix,jy).gt.smallnum) then
if (sp_zer.eqv..true.) then ! first non zero value
sp_count_i=sp_count_i+1
......@@ -378,7 +378,7 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
CALL SYSTEM_CLOCK(count_clock)
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Dry deposition
! Dry deposition
sp_count_i=0
sp_count_r=0
sp_fact=-1.
......@@ -420,9 +420,9 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
endif
! Concentrations
! Concentrations
! surf_only write only 1st layer
! surf_only write only 1st layer
sp_count_i=0
sp_count_r=0
......@@ -444,8 +444,8 @@ subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact* &
grid(ix,jy,kz)* &
factor3d(ix,jy,kz)/tot_mu(ks,kp)
! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)