Commit 8f52cb4e authored by Sabine's avatar Sabine
Browse files

first version of the particledump for backward simulation in netcdf

parent 187be914
......@@ -78,13 +78,14 @@ subroutine boundcond_domainfill(itime,loutend)
!********************************************************************
do i=1,numpart
if (itra1(i).eq.itime) then
if ((ytra1(i).gt.real(ny_sn(2))).or. &
(ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
((xtra1(i).lt.real(nx_we(1))).or. &
(xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
endif
! keep the trajectories alive outside domain
! if (itra1(i).eq.itime) then
! if ((ytra1(i).gt.real(ny_sn(2))).or. &
! (ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
! if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
! ((xtra1(i).lt.real(nx_we(1))).or. &
! (xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
! endif
if (itra1(i).ne.-999999999) numactiveparticles= &
numactiveparticles+1
end do
......
......@@ -20,7 +20,7 @@ module com_mod
use par_mod, only: dp, numpath, maxnests, maxageclass, maxspec, ni, &
numclass, nymax, nxmax, maxcolumn, maxwf, nzmax, nxmaxn, nymaxn, &
maxreceptor, maxpart, maxrand, nwzmax, nuvzmax
maxreceptor, maxpart, maxrand, nwzmax, nuvzmax, maxoutsteps
implicit none
......@@ -641,6 +641,15 @@ module com_mod
integer :: idt(maxpart),itramem(maxpart),itrasplit(maxpart)
integer :: numparticlecount
! for the particledump
integer :: npart_av(maxpart), part_av_itime(maxpart,maxoutsteps), part_av_filled(maxpart)
real :: part_av_cartx(maxpart),part_av_carty(maxpart),&
part_av_cartz(maxpart),part_av_z(maxpart,maxoutsteps),part_av_topo(maxpart,maxoutsteps), &
part_av_x(maxpart,maxoutsteps),part_av_y(maxpart,maxoutsteps)
real :: part_av_pv(maxpart,maxoutsteps),part_av_qv(maxpart,maxoutsteps),part_av_tt(maxpart,maxoutsteps),&
part_av_rho(maxpart,maxoutsteps),part_av_tro(maxpart,maxoutsteps),part_av_hmix(maxpart,maxoutsteps)
real :: part_av_uu(maxpart,maxoutsteps),part_av_vv(maxpart,maxoutsteps)!,part_av_energy(maxpart,maxoutsteps)
real(kind=dp) :: xtra1(maxpart),ytra1(maxpart)
real :: ztra1(maxpart),xmass1(maxpart,maxspec)
......
......@@ -106,7 +106,9 @@ outgrid_init_nest.o \
part0.o \
partdep.o \
partoutput.o \
partoutput_average.o \
partoutput_short.o \
partpos_average.o \
pbl_profile.o \
plumetraj.o \
psih.o \
......
......@@ -66,7 +66,7 @@ module netcdf_output_mod
private
public :: writeheader_netcdf,concoutput_surf_nest_netcdf,concoutput_netcdf,&
&concoutput_nest_netcdf,concoutput_surf_netcdf
&concoutput_nest_netcdf,concoutput_surf_netcdf,writeaveragepart_netcdf
! include 'netcdf.inc'
......@@ -230,6 +230,65 @@ subroutine nf90_err(status)
end subroutine nf90_err
subroutine writeaveragepart_netcdf(adate,atime)
implicit none
character :: adate*8,atime*6
character(len=255) :: fname
integer :: cache_size, ncid, i, partused
integer npointDimID,topoID,zID,xID,yID,pvID,qvID,ttID,uuID,vvID,rhoID,troID,hmixID,itimeID,timestepsDimID
!the other dimesion should be the maximum ageclass/outputtimestep
!the first dimension is maxpart, particles are written at npoint(numpart)
!find out how many particles are filled
partused=0
do i=1,maxpart
if (part_av_filled(i).gt.0) partused=partused+1
enddo
fname=path(2)(1:length(2))//'partposit_average_'//adate//atime//'.nc'
cache_size = 16 * partused
call nf90_err(nf90_create(trim(fname), cmode = nf90_hdf5, ncid = ncid, &
cache_size = cache_size))
call nf90_err(nf90_def_dim(ncid, 'numpart', partused, npointDimID))
call nf90_err(nf90_def_dim(ncid, 'timesteps', 20, timestepsDimID))
call nf90_err(nf90_def_var(ncid, 'topo', nf90_float, (/ npointDimID,timestepsDimID /), topoID))
call nf90_err(nf90_def_var(ncid, 'z', nf90_float, (/ npointDimID,timestepsDimID /), zID))
call nf90_err(nf90_def_var(ncid, 'x', nf90_float, (/ npointDimID,timestepsDimID /), xID))
call nf90_err(nf90_def_var(ncid, 'y', nf90_float, (/ npointDimID,timestepsDimID /), yID))
call nf90_err(nf90_def_var(ncid, 'pv', nf90_float, (/ npointDimID,timestepsDimID /), pvID))
call nf90_err(nf90_def_var(ncid, 'qv', nf90_float, (/ npointDimID,timestepsDimID /), qvID))
call nf90_err(nf90_def_var(ncid, 'tt', nf90_float, (/ npointDimID,timestepsDimID /), ttID))
call nf90_err(nf90_def_var(ncid, 'uu', nf90_float, (/ npointDimID,timestepsDimID /), uuID))
call nf90_err(nf90_def_var(ncid, 'vv', nf90_float, (/ npointDimID,timestepsDimID /), vvID))
call nf90_err(nf90_def_var(ncid, 'rho', nf90_float, (/ npointDimID,timestepsDimID /), rhoID))
call nf90_err(nf90_def_var(ncid, 'tro', nf90_float, (/ npointDimID,timestepsDimID /), troID))
call nf90_err(nf90_def_var(ncid, 'hmix', nf90_float, (/ npointDimID,timestepsDimID /), hmixID))
call nf90_err(nf90_def_var(ncid, 'itime',nf90_float, (/ npointDimID,timestepsDimID /), itimeID))
call nf90_err(nf90_put_var(ncid, topoID, part_av_topo(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, zID, part_av_z(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, xID, part_av_x(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, yID, part_av_y(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, pvID, part_av_pv(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, qvID, part_av_qv(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, ttID, part_av_tt(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, uuID, part_av_uu(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, rhoID, part_av_rho(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, troID, part_av_tro(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, hmixID, part_av_hmix(1:partused,1:maxoutsteps)))
call nf90_err(nf90_put_var(ncid, itimeID, part_av_itime(1:partused,1:maxoutsteps)))
call nf90_err(nf90_close(ncid))
end subroutine writeaveragepart_netcdf
!****************************************************************
! Create netcdf file and write header/metadata information
! lnest = .false. : Create main output file
......
......@@ -229,8 +229,9 @@ module par_mod
! Maximum number of particles, species, and similar
!**************************************************
integer,parameter :: maxpart=20000000
integer,parameter :: maxspec=2
integer,parameter :: maxoutsteps=20
integer,parameter :: maxpart=100000
integer,parameter :: maxspec=1
real,parameter :: minmass=0.0001
......@@ -285,7 +286,7 @@ module par_mod
!************************************
integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93,unitpartout_average=105
integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
integer,parameter :: unitoutgrid=97, unitoutgridppt=99, unitoutinfo=1
......
......@@ -116,8 +116,7 @@ subroutine partoutput_short(itime)
if (xlon.lt.-180.) xlon=xlon+360.
numshortall=numshortall+1
if ((xlon.gt.-140).and.(xlon.lt.60).and.(ylat.gt.10).and. &
(xmass1(i,1).gt.0.)) then
if (xmass1(i,1).gt.0.) then
numshortout=numshortout+1
idump(1,numshortout)=nint(xlon*180.)
idump(2,numshortout)=nint(ylat*360.)
......
......@@ -342,7 +342,7 @@ subroutine readcommand
! For quasilag and domainfill ioutputforechrelease is forbidden
!*****************************************************************************
if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0)) then
if ((ldirect.lt.0).and.(ioutputforeachrelease.eq.0).and.(mdomainfill.eq.0)) then
write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
write(*,*) '#### FOR BACKWARD RUNS, IOUTPUTFOREACHRLEASE ####'
write(*,*) '#### MUST BE SET TO ONE! ####'
......
......@@ -421,6 +421,9 @@ subroutine readreleases
ireleaseend(numpoint)=int((jul2-bdate)*86400.)
endif
endif
!calculate releasetime also for domainfilling
ireleasestart(numpoint)=int((jul1-bdate)*86400.)
ireleaseend(numpoint)=int((jul2-bdate)*86400.)
! Check, whether the total number of particles may exceed totally allowed
......
......@@ -170,8 +170,12 @@ subroutine timemanager
if (mdomainfill.ge.1) then
if (itime.eq.0) then
call init_domainfill
else
call boundcond_domainfill(itime,loutend)
else !only works backward, switched on for partposit!
write(*,*) 'boundconc: ',itime,ireleasestart(1), ireleaseend(1)
if ( (abs(itime).le.abs(ireleasestart(1))).and. &
(abs(itime).gt.abs( ireleaseend(1))) ) then
call boundcond_domainfill(itime,loutend)
endif
endif
else
call releaseparticles(itime)
......@@ -326,7 +330,13 @@ subroutine timemanager
drygridtotalunc
45 format(i9,' SECONDS SIMULATED: ',i8, &
' PARTICLES: Uncertainty: ',3f7.3)
if (ipout.ge.1) call partoutput(itime) ! dump particle positions
! if (ipout.ge.1) call partoutput(itime) ! dump particle positions
! if (ipout.ge.1) call partoutput_short(itime) ! dump particle positions in extremely compressed format
if (ipout.ge.1) then
if (mod(itime,345600).eq.0) call partoutput(itime) ! dump particle positions
call partoutput_average(itime) ! dump particle positions
endif
loutnext=loutnext+loutstep
loutstart=loutnext-loutaver/2
loutend=loutnext+loutaver/2
......@@ -434,6 +444,11 @@ subroutine timemanager
us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
cbt(j))
! Calculate average position for particle dump output
!****************************************************
call partpos_average(itime,j)
! Calculate the gross fluxes across layer interfaces
!***************************************************
......
Supports Markdown
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