Commit 20963b15 authored by Espen Sollum's avatar Espen Sollum
Browse files

Backwards deposition for the MPI version

parent 3f149ccf
......@@ -66,7 +66,8 @@ program flexpart
character(len=256) :: inline_options !pathfile, flexversion, arg2
integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
integer :: detectformat
integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
& tot_i=0
! Initialize mpi
......@@ -205,11 +206,11 @@ program flexpart
metdata_format = detectformat()
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
print *,'ECMWF metdata detected'
if (lroot) print *,'ECMWF metdata detected'
elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
print *,'NCEP metdata detected'
if (lroot) print *,'NCEP metdata detected'
else
print *,'Unknown metdata format'
if (lroot) print *,'Unknown metdata format'
stop
endif
......@@ -482,28 +483,24 @@ program flexpart
! NIK 16.02.2005
if (lroot) then
call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
call MPI_Reduce(tot_inc_count, tot_i, 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, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
& mp_comm_used, mp_ierr)
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(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& tot_blc_count(i)
& tot_b(i)
! & tot_blc_count(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& tot_inc_count(i)
& tot_i(i)
! & tot_inc_count(i)
write(*,*) '**********************************************'
end do
......
......@@ -139,6 +139,8 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
write(unitdates,'(a)') adate//atime
close(unitdates)
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
IF (init) THEN
file_stat='OLD'
init=.false.
......@@ -312,7 +314,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
drygridsigma(ix,jy)= &
drygridsigma(ix,jy)* &
sqrt(real(nclassunc))
125 drygridsigmatotal=drygridsigmatotal+ &
drygridsigmatotal=drygridsigmatotal+ &
drygridsigma(ix,jy)
endif
......
......@@ -103,7 +103,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
real,parameter :: weightair=28.97
logical :: sp_zer
LOGICAL,save :: init=.true.
logical,save :: init=.true.
character :: adate*8,atime*6
character(len=3) :: anspec
integer :: mind
......@@ -127,7 +127,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
! If 'dates' file exists, make a backup
! If 'dates' file exists in output directory, make a backup
inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
if (ldates_file.and.init) then
open(unit=unitdates, file=path(2)(1:length(2))//'dates',form='formatted', &
......@@ -256,23 +256,34 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
do ks=1,nspec
write(anspec,'(i3.3)') ks
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
if (DRYBKDEP.or.WETBKDEP) then !scavdep output
if (DRYBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
atime//'_'//anspec,form='unformatted')
if (WETBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgrid) itime
else
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
write(unitoutgridppt) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime
endif
endif ! if deposition output
do kp=1,maxpointspec_act
do nage=1,nageclass
......@@ -352,7 +363,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Concentration output
!*********************
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
! Wet deposition
sp_count_i=0
......@@ -447,10 +458,18 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
if (lparticlecountoutput) then
sparse_dump_r(sp_count_r)= &
sp_fact* &
grid(ix,jy,kz)
else
sparse_dump_r(sp_count_r)= &
sp_fact* &
grid(ix,jy,kz)* &
factor3d(ix,jy,kz)/tot_mu(ks,kp)
end if
! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
! sparse_dump_u(sp_count_r)=
......@@ -636,24 +655,8 @@ 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
creceptor(:,:)=0.
gridunc(:,:,:,:,:,:,:)=0.
if (mp_measure_time) call mpif_mtime('rootonly',1)
......
......@@ -232,6 +232,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end do
40 indmin=indj
if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
call writeprecip(itime,memind(1))
endif
else
! No wind fields, which can be used, are currently in memory
......@@ -286,7 +290,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end do
60 indmin=indj
mind3=memstat
! if (WETBKDEP.and.lroot) then
if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
call writeprecip(itime,memind(1))
endif
endif
......
......@@ -331,16 +331,20 @@ subroutine readcommand
ind_rel = 0
case (3) ! 3 .. wet deposition in outputfield
ind_rel = 3
write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE #### '
write(*,*) ' #### Releaseheight is forced to 0 - 20km #### '
write(*,*) ' #### Release is performed above ground lev #### '
if (lroot) then
write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE #### '
write(*,*) ' #### Releaseheight is forced to 0 - 20km #### '
write(*,*) ' #### Release is performed above ground lev #### '
end if
WETBKDEP=.true.
allocate(xscav_frac1(maxpart,maxspec))
case (4) ! 4 .. dry deposition in outputfield
ind_rel = 4
write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### '
write(*,*) ' #### Releaseheight is forced to 0 - 2*href #### '
write(*,*) ' #### Release is performed above ground lev #### '
if (lroot) then
write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### '
write(*,*) ' #### Releaseheight is forced to 0 - 2*href #### '
write(*,*) ' #### Release is performed above ground lev #### '
end if
DRYBKDEP=.true.
allocate(xscav_frac1(maxpart,maxspec))
end select
......
......@@ -102,12 +102,12 @@ subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
! coordinate parameters
integer :: isec1(56),isec2(22+nxmax+nymax)
real(kind=4) :: zsec4(jpunp)
real(kind=4) :: xaux,yaux
real(kind=8) :: xauxin,yauxin
real,parameter :: eps=1.e-4
real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
real(sp) :: zsec4(jpunp)
real(sp) :: xaux,yaux
real(dp) :: xauxin,yauxin
real(sp),parameter :: eps=1.e-4
real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
logical :: hflswitch,strswitch !,readcloud
......@@ -122,7 +122,7 @@ subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
strswitch=.false.
!ZHG test the grib fields that have lcwc without using them
! readcloud=.false.
!hg end
levdiff2=nlev_ec-nwz+1
iumax=0
iwmax=0
......
......@@ -215,7 +215,12 @@ subroutine releaseparticles(itime)
do k=1,nspec
xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
*timecorrect(k)/average_timecorrect
! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
! if ( henry(k).gt.0 .or. &
! crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
! ccn_aero(k).gt.0. .or. in_aero(k).gt.0. ) then
xscav_frac1(ipart,k)=-1.
endif
! Assign certain properties to particle
!**************************************
end do
......@@ -365,7 +370,7 @@ subroutine releaseparticles(itime)
!Af ind_rel is defined in readcommand.f
if (ind_rel .eq. 1) then
if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
! Interpolate the air density
!****************************
......
......@@ -114,16 +114,16 @@ subroutine timemanager(metdata_format)
! integer :: ksp
integer :: ip
integer :: loutnext,loutstart,loutend
integer :: ix,jy,ldeltat,itage,nage
integer :: ix,jy,ldeltat,itage,nage,idummy
integer :: i_nan=0,ii_nan,total_nan_intl=0 !added by mc to check instability in CBL scheme
integer :: numpart_tot_mpi ! for summing particles on all processes
real :: outnum,weight,prob(maxspec)
real :: decfact
real :: outnum,weight,prob(maxspec), prob_rec(maxspec), decfact,wetscav
real(sp) :: gridtotalunc
real(dep_prec) :: drygridtotalunc=0_dep_prec,wetgridtotalunc=0_dep_prec,&
& drydeposit(maxspec)=0_dep_prec
real :: xold,yold,zold,xmassfract
real :: grfraction(3)
real, parameter :: e_inv = 1.0/exp(1.0)
! Measure time spent in timemanager
......@@ -158,6 +158,8 @@ subroutine timemanager(metdata_format)
! print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
if (.not.lusekerneloutput) write(*,*) 'Not using the kernel'
if (turboff) write(*,*) 'Turbulence switched off'
do itime=0,ideltas,lsynctime
......@@ -707,6 +709,42 @@ subroutine timemanager(metdata_format)
yold=ytra1(j)
zold=ztra1(j)
! RECEPTOR: dry/wet depovel
!****************************
! Before the particle is moved
! the calculation of the scavenged mass shall only be done once after release
! xscav_frac1 was initialised with a negative value
if (DRYBKDEP) then
do ks=1,nspec
if ((xscav_frac1(j,ks).lt.0)) then
call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
if (DRYDEPSPEC(ks)) then ! dry deposition
xscav_frac1(j,ks)=prob_rec(ks)
else
xmass1(j,ks)=0.
xscav_frac1(j,ks)=0.
endif
endif
enddo
endif
if (WETBKDEP) then
do ks=1,nspec
if ((xscav_frac1(j,ks).lt.0)) then
call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
if (wetscav.gt.0) then
xscav_frac1(j,ks)=wetscav* &
(zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
else
xmass1(j,ks)=0.
xscav_frac1(j,ks)=0.
endif
endif
enddo
endif
! Integrate Lagevin equation for lsynctime seconds
!*************************************************
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment