Commit 08a38b55 authored by Espen Sollum's avatar Espen Sollum
Browse files

Added (hardcoded) option to output number of particles per grid cell.

parent bb579a9d
......@@ -192,17 +192,24 @@ subroutine conccalc(itime,weight)
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
end do
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
if (lparticlecountoutput) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
end do
else
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
end if
endif
endif
......@@ -335,14 +342,21 @@ subroutine conccalc(itime,weight)
xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
if (lparticlecountoutput) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
end do
else
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
endif
endif
endif
else ! attribution via uniform kernel
ddx=xl-real(ix) ! distance to left cell border
......
......@@ -204,11 +204,26 @@ subroutine conccalc(itime,weight)
(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. &
(jy.le.numygrid-1)) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
end do
else
if (lparticlecountoutput) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
end do
else
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
end if
endif
endif
else ! attribution via uniform kernel
......@@ -238,20 +253,36 @@ subroutine conccalc(itime,weight)
if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
w=wx*wy
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
w=wx*(1.-wy)
do ks=1,nspec
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
endif
......@@ -259,20 +290,36 @@ subroutine conccalc(itime,weight)
if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
w=(1.-wx)*(1.-wy)
do ks=1,nspec
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
w=(1.-wx)*wy
do ks=1,nspec
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
endif
endif
......@@ -303,11 +350,26 @@ subroutine conccalc(itime,weight)
(yl.gt.real(numygridn-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
(jy.le.numygridn-1)) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
end do
else
if (lparticlecountoutput) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
end do
else
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
end do
endif
endif
endif
else ! attribution via uniform kernel
......@@ -337,20 +399,36 @@ subroutine conccalc(itime,weight)
if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
if ((jy.ge.0).and.(jy.le.numygridn-1)) then
w=wx*wy
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
w=wx*(1.-wy)
do ks=1,nspec
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
endif
......@@ -358,20 +436,36 @@ subroutine conccalc(itime,weight)
if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
w=(1.-wx)*(1.-wy)
do ks=1,nspec
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
if ((jy.ge.0).and.(jy.le.numygridn-1)) then
w=(1.-wx)*wy
do ks=1,nspec
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
end do
else
do ks=1,nspec
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
end do
end do
endif
endif
endif
endif
......
......@@ -445,10 +445,18 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact=sp_fact*(-1.)
endif
sp_count_r=sp_count_r+1
sparse_dump_r(sp_count_r)= &
sp_fact* &
grid(ix,jy,kz)* &
factor3d(ix,jy,kz)/tot_mu(ks,kp)
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)=
......
......@@ -59,8 +59,13 @@ module par_mod
logical, parameter :: lnokernel=.false.
!*********************************************************************
! Set to T to change output units to number of particles per grid cell
!*********************************************************************
logical, parameter :: lparticlecountoutput=.false.
!***********************************************************
! Number of directories/files used for FLEXPART input/output
! number of directories/files used for FLEXPART input/output
!***********************************************************
integer,parameter :: numpath=4
......@@ -140,21 +145,21 @@ module par_mod
! Maximum dimensions of the input mother grids
!*********************************************
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new
integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
! integer,parameter :: nxmax=181,nymax=91,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
! ECMWF
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=359 ! 1.0 degree 92 level
integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 1.0 degree 138 level
! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=359 ! 0.5 degree 138 level
! integer,parameter :: nxmax=181,nymax=91,nuvzmax=92,nwzmax=92,nzmax=92,nxshift=0 ! CERA 2.0 degree 92 level
! INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data
! GFS
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138,nxshift=0
! integer,parameter :: nxshift=359 ! for ECMWF
integer,parameter :: nxshift=0 ! for GFS
!*********************************************
! Maximum dimensions of the nested input grids
!*********************************************
integer,parameter :: maxnests=1,nxmaxn=451,nymaxn=226
integer,parameter :: maxnests=0,nxmaxn=451,nymaxn=226
! nxmax,nymax maximum dimension of wind fields in x and y
! direction, respectively
......@@ -212,7 +217,7 @@ module par_mod
! Maximum number of particles, species, and similar
!**************************************************
integer,parameter :: maxpart=10000000
integer,parameter :: maxpart=3000000
integer,parameter :: maxspec=4
real,parameter :: minmass=0.0001
......
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