From 2d7717b12556da7b12bfb09b176d092e07997f59 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Thu, 25 Apr 2024 21:11:47 +0200
Subject: [PATCH] Changes to make the parallelization of receptoroutput more
efficient and correction to restart a run using partoutput
Former-commit-id: f50831436f74aeed3ab189eccaa6c593927250cf
---
src/FLEXPART.f90 | 5 +-
src/chemistry_mod.f90 | 112 ++++++++++++-------------------
src/com_mod.f90 | 6 +-
src/emissions_mod.f90 | 83 ++++++++++++-----------
src/gitversion.txt | 1 +
src/initdomain_mod.f90 | 27 +++++---
src/netcdf_output_mod.f90 | 14 +++-
src/readoptions_mod.f90 | 4 +-
src/receptor_mod.f90 | 137 ++++++++++++++++++++------------------
src/timemanager_mod.f90 | 24 ++++---
src/totals_mod.f90 | 1 +
11 files changed, 212 insertions(+), 202 deletions(-)
create mode 100644 src/gitversion.txt
diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index e6d865e6..7a20ec48 100644
--- a/src/FLEXPART.f90
+++ b/src/FLEXPART.f90
@@ -44,7 +44,7 @@ program flexpart
real :: s_timemanager
character(len=256) :: &
inline_options ! pathfile, flexversion, arg2
- character(len=256) :: gitversion_tmp="undefined"
+ character(len=256) :: gitversion_tmp="8d9c680 Mon Apr 22 15:37:39 2024 +0200"
! Keeping track of the total running time of FLEXPART, printed out at the end.
!*****************************************************************************
@@ -300,11 +300,10 @@ subroutine read_options_and_initialise_flexpart
numreceptor=0
numsatreceptor=0
nlayermax=1
- ! NOTE: need to change to calc receptors conc without netcdf
#ifdef USE_NCF
call readreceptors_satellite
- call readreceptors
#endif
+ call readreceptors
! Read the landuse inventory
!***************************
diff --git a/src/chemistry_mod.f90 b/src/chemistry_mod.f90
index 19abd5f2..403b9890 100644
--- a/src/chemistry_mod.f90
+++ b/src/chemistry_mod.f90
@@ -373,11 +373,13 @@ module chemistry_mod
integer :: jrx, jry
real(kind=dp) :: jul
!! testing
- character(len=4) :: atime
- character(len=20) :: frmt
- real, dimension(nxjr,nyjr) :: jscalar
+! character(len=4) :: atime
+! character(len=20) :: frmt
+! real, dimension(nxjr,nyjr) :: jscalar, sza_grid
- jscalar(:,:)=1.
+ !! testing
+! jscalar(:,:)=1.
+! sza_grid(:,:)=0.
! current hour of simulation
curhour=itime/3600
@@ -397,7 +399,7 @@ module chemistry_mod
dt2=float(memCLtime(2)-interp_time)
dtt=1./(dt1+dt2)
!! testing
- print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt
+! print*, 'getchemhourly: dt1, dt2, dtt = ',dt1,dt2,dtt
do nr=1,nreagent
write(*,*) 'Interpolating fields for reagent: ',reagents(nr)
clfield_cur(:,:,:,nr)=(dt2*CL_field(:,:,:,nr,1) + dt1*CL_field(:,:,:,nr,2))*dtt
@@ -410,13 +412,15 @@ module chemistry_mod
do jy=1,nyCL
! jrate_average dimensions given as grid midpoints
jry=int((latCL(jy)-(latjr(1)-0.5*dyjr))/dyjr)+1
+ jry=min(jry,nyCL)
!! testing
- if (kz.eq.1.and.jy.lt.10) print*, 'latCL, latjr = ',latCL(jy), latjr(jry)
+! if (kz.eq.1.and.jy.lt.10) print*, 'latCL, latjr = ',latCL(jy), latjr(jry)
do ix=1,nxCL
! jrate_average dimensions given as grid midpoints
jrx=int((lonCL(ix)-(lonjr(1)-0.5*dxjr))/dxjr)+1
+ jrx=min(jrx,nxCL)
!! testing
- if (kz.eq.1.and.jy.lt.10.and.ix.lt.10) print*, 'lonCL, lonjr = ',lonCL(ix), lonjr(jrx)
+! if (kz.eq.1.and.jy.lt.10.and.ix.lt.10) print*, 'lonCL, lonjr = ',lonCL(ix), lonjr(jrx)
! solar zenith angle
sza=zenithangle(latjr(jry),lonjr(jrx),jjjjmmdd,hhmmss)
! calculate J(O1D) (jrate)
@@ -427,7 +431,10 @@ module chemistry_mod
if(jrate_cur.gt.0.) then
clfield_cur(ix,jy,kz,nr)=clfield_cur(ix,jy,kz,nr)*jrate/jrate_cur
!! testing
- if (kz.eq.1) jscalar(ix,jy)=jrate/jrate_cur
+! if (kz.eq.1) then
+! jscalar(ix,jy)=jrate/jrate_cur
+! sza_grid(ix,jy)=sza
+! endif
endif
end do
end do
@@ -439,29 +446,27 @@ module chemistry_mod
endif ! curhour
!! testing
- write(frmt,fmt='(A,I4,A)') '(',nxCL,'(E14.6))'
- write(atime,fmt='(I4.4)') curhour
- open(600,file=path(2)(1:length(2))//'clfield_'//atime//'.txt',action='write',status='replace')
- do kz=1,nzCL
- do jy=1,nyCL
- write(600,fmt=frmt) clfield_cur(:,jy,kz,1)
- end do
- end do
- close(600)
- write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
- open(600,file=path(2)(1:length(2))//'jscalar_'//atime//'.txt',action='write',status='replace')
- do jy=1,nyjr
- write(600,fmt=frmt) jscalar(:,jy)
- end do
- close(600)
- if (itime.eq.0) then
- write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
- open(600,file=path(2)(1:length(2))//'javerage_'//atime//'.txt',action='write',status='replace')
- do jy=1,nyjr
- write(600,fmt=frmt) jrate_average(:,jy,1)
- end do
- endif
- close(600)
+! write(frmt,fmt='(A,I4,A)') '(',nxCL,'(E14.6))'
+! write(atime,fmt='(I4.4)') curhour
+! open(600,file=path(2)(1:length(2))//'clfield_'//atime//'.txt',action='write',status='replace')
+! do kz=1,nzCL
+! do jy=1,nyCL
+! write(600,fmt=frmt) clfield_cur(:,jy,kz,1)
+! end do
+! end do
+! close(600)
+! write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
+! open(600,file=path(2)(1:length(2))//'jscalar_'//atime//'.txt',action='write',status='replace')
+! do jy=1,nyjr
+! write(600,fmt=frmt) jscalar(:,jy)
+! end do
+! close(600)
+! write(frmt,fmt='(A,I4,A)') '(',nxjr,'(E14.6))'
+! open(600,file=path(2)(1:length(2))//'sza_'//atime//'.txt',action='write',status='replace')
+! do jy=1,nyjr
+! write(600,fmt=frmt) sza_grid(:,jy)
+! end do
+! close(600)
end subroutine getchemhourly
@@ -521,15 +526,15 @@ module chemistry_mod
#ifdef _OPENMP
allocate( chem_loss_tmp(nreagent,nspec,numthreads) )
chem_loss_tmp(:,:,:) = 0d0
-#endif
+#endif
+
! Loop over particles
!*****************************************
!$OMP PARALLEL &
!$OMP PRIVATE(ii,jpart,ngrid,j,xtn,ytn,ix,jy, &
!$OMP xlon,ylat,clx,cly,clz,clzm,kz,altCLtop,dz1,dz2,dzz,nr,i, &
-!$OMP cl_cur,indz,temp,ks,clrate,restmass,clreacted,ithread)
-! !$OMP REDUCTION(+:chem_loss)
+!$OMP cl_cur,indz,temp,ks,clrate,restmass,clreacted,ithread)
#ifdef _OPENMP
ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
@@ -576,8 +581,8 @@ module chemistry_mod
endif
! get position in the chem field
! assumes chem field dimensions given as grid midpoints
- clx=min(nxCL, int((xlon-(lonCL(1)-0.5*dxCL))/dxCL)+1)
- cly=min(nyCL, int((ylat-(latCL(1)-0.5*dyCL))/dyCL)+1)
+ clx=min(nxCL,int((xlon-(lonCL(1)-0.5*dxCL))/dxCL)+1)
+ cly=min(nyCL,int((ylat-(latCL(1)-0.5*dyCL))/dyCL)+1)
! get the level of the chem field for the particle
! z is the z-coord of the trajectory above model orography in metres
@@ -655,11 +660,11 @@ module chemistry_mod
!$OMP END DO
!$OMP END PARALLEL
-
+
#ifdef _OPENMP
do ithread=1,numthreads
chem_loss(:,:) = chem_loss(:,:)+chem_loss_tmp(:,:,ithread)
- end do
+ end do
deallocate( chem_loss_tmp )
#endif
@@ -848,37 +853,6 @@ module chemistry_mod
end function zenithangle
- function find_minloc(n,vect,val) result(ind)
-
- !*****************************************************************************
- ! *
- ! This function returns index of a vector corresponding to the closest *
- ! value to variable val *
- ! *
- ! Author: R. Thompson, Mar 2024 *
- ! *
- !*****************************************************************************
-
- implicit none
-
- integer, intent(in) :: n
- real, dimension(n), intent(in) :: vect
- real, intent(in) :: val
- integer :: i, ind
- real :: diff, diff_min
-
- diff_min=huge(1.0)
-
- ind=0
- do i=1,n
- diff=abs(vect(i)-val)
- if ( diff.lt.diff_min ) then
- diff_min=diff
- ind=i
- endif
- end do
-
- end function find_minloc
end module chemistry_mod
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 1a7957d4..3d1eec5e 100644
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -419,7 +419,7 @@ module com_mod
real, allocatable, dimension(:) :: xkreceptor,nnreceptor
character(len=16), allocatable, dimension(:) :: receptorname
integer :: cpointer(maxrecsample)
- integer :: numreceptor
+ integer :: numreceptor, numcurrec
! satellite receptors
real, allocatable, dimension(:) :: xsatellite,ysatellite
@@ -429,7 +429,7 @@ module com_mod
real, allocatable, dimension(:,:,:) :: csatellite, csatuncert
real, allocatable, dimension(:,:) :: xksatellite, nnsatellite
character(len=16), allocatable, dimension(:) :: satellitename
- integer :: numsatreceptor, nlayermax, numsatellite
+ integer :: numsatreceptor, nlayermax, numsatellite, numcursat
integer, allocatable, dimension(:) :: nnsatlayer
integer :: csatpointer(maxrecsample)
@@ -437,7 +437,9 @@ module com_mod
! creceptor concentrations at receptor points
! receptorarea area of 1*1 grid cell at receptor point
! numreceptor number of receptors (non-satellite)
+ ! numcurrec number of receptors in current time interval (updated each time interval)
! numsatreceptor number of satellite receptors (aka. retrievals)
+ ! numcursat number of satellite receptors in current time interval (updated each time interval)
! numsatellite number of satellite instruments
! nlayermax max number of vertical layers in satellite retrievals
! nnsatlayer actual number of vertical layers for each satellite
diff --git a/src/emissions_mod.f90 b/src/emissions_mod.f90
index e0414708..775b894a 100644
--- a/src/emissions_mod.f90
+++ b/src/emissions_mod.f90
@@ -69,13 +69,11 @@ module emissions_mod
character(len=20) :: frmt
real :: em_frac
real, allocatable, dimension(:,:,:) :: em_neg
-#ifdef _OPENMP
real, allocatable, dimension(:,:) :: tot_em_up_tmp
real, allocatable, dimension(:,:,:,:) :: mass_field_tmp, em_neg_tmp
-#endif
real :: f_m
logical :: lexist
- integer, parameter :: unittest=120
+! integer, parameter :: unittest=120
! fraction of stored emissions that is released in a time step
@@ -91,31 +89,28 @@ module emissions_mod
dt2=float(memtime(2)-itime)
dtt=1./(dt1+dt2)
+ tot_em_up(:) = 0.
+
! estimate mass in PBL from particle positions
!**********************************************
+ mass_field(:,:,:) = 0.
allocate( em_neg(nspec-1,nxem,nyem) )
-#ifdef _OPENMP
+#if _OPENMP
allocate( mass_field_tmp(nspec,nxem,nyem,numthreads))
allocate( em_neg_tmp(nspec-1,nxem,nyem,numthreads))
allocate( tot_em_up_tmp(nspec,numthreads) )
-#endif
-
-!$OMP PARALLEL
-!$OMP WORKSHARE
mass_field_tmp(:,:,:,:) = 0.
em_neg_tmp(:,:,:,:) = 0.
tot_em_up_tmp(:,:) = 0.
+#endif
tot_em_up(:) = 0.
mass_field(:,:,:) = 0.
em_neg(:,:,:) = 0.
-!$OMP END WORKSHARE
-!$OMP END PARALLEL
!$OMP PARALLEL &
!$OMP PRIVATE(ii,xlon,ylat,em_ix,em_jy,ix,jy,ixp,jyp,ddx,ddy, &
-!$OMP rddx,rddy,p1,p2,p3,p4,mm,hm,hmixi,ks,ithread) &
-!$OMP SHARED(em_cond,mass_field_tmp)
+!$OMP rddx,rddy,p1,p2,p3,p4,mm,hm,hmixi,ks,ithread)
#ifdef _OPENMP
ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
@@ -124,8 +119,6 @@ module emissions_mod
#endif
!$OMP DO
-
-! REDUCTION(+:mass_field)
do ii=1,count%alive ! loop over all particles
xlon=xlon0+part(ii)%xlon*dx
@@ -135,7 +128,7 @@ module emissions_mod
em_ix=min(nxem, int((xlon-(lonem(1)-0.5*dxem))/dxem)+1)
em_jy=min(nyem, int((ylat-(latem(1)-0.5*dyem))/dyem)+1)
!! testing
- if (ii.lt.20) print*, 'lonem, lon, latem, lat = ',lonem(em_ix),xlon,latem(em_jy),ylat
+! if (ii.lt.20) print*, 'lonem, lon, latem, lat = ',lonem(em_ix),xlon,latem(em_jy),ylat
! interpolate to particle position
ix=int(part(ii)%xlon)
@@ -180,10 +173,13 @@ module emissions_mod
!$OMP END DO
!$OMP END PARALLEL
- ! Manual reduction of mass_field
+#ifdef _OPENMP
+ ! manual reduction of mass_field
do ithread=1,numthreads
mass_field(:,:,:) = mass_field(:,:,:)+mass_field_tmp(:,:,:,ithread)
- end do
+ end do
+#endif
+
f_m = exp(-1.*real(lsynctime)/tau_ipm)
! estimate emissions for each particle
@@ -191,8 +187,7 @@ module emissions_mod
!$OMP PARALLEL &
!$OMP PRIVATE(ii,xlon,ylat,em_ix,em_jy, &
-!$OMP ks,em_cur,tmp,ithread) &
-!$OMP SHARED(em_cond,tot_em_up_tmp,em_neg_tmp)
+!$OMP ks,em_cur,tmp,ithread)
#ifdef _OPENMP
ithread = OMP_GET_THREAD_NUM()+1 ! Starts with 1
@@ -200,8 +195,7 @@ module emissions_mod
ithread = 1
#endif
-!$OMP DO
-! REDUCTION(+:tot_em_up,em_neg)
+!$OMP DO
do ii=1,count%alive ! loop over particles
if (.not.em_cond(ii)) cycle ! skip particles not in PBL
@@ -233,9 +227,13 @@ module emissions_mod
if (-1.*tmp.gt.mass(ii,ks)) then
tmp = tmp + mass(ii,ks)
! subtract mass from atmosphere by setting it to zero below
+#ifdef _OPENMP
tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) - real(mass(ii,ks),kind=dp)
+#else
+ tot_em_up(ks) = tot_em_up(ks) - real(mass(ii,ks),kind=dp)
+#endif
mass(ii,ks) = 0.
- ! add remaining uptake to em_neg
+ ! add remaining uptake to em_neg
#ifdef _OPENMP
em_neg_tmp(ks-1,em_ix,em_jy,ithread) = em_neg_tmp(ks-1,em_ix,em_jy,ithread) + &
tmp/mass(ii,1)*mass_field(1,em_ix,em_jy)
@@ -245,11 +243,19 @@ module emissions_mod
#endif
else
mass(ii,ks)=mass(ii,ks)+tmp
+#ifdef _OPENMP
tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) + real(tmp,kind=dp)
+#else
+ tot_em_up(ks) = tot_em_up(ks) + real(tmp,kind=dp)
+#endif
endif
else
mass(ii,ks)=mass(ii,ks)+tmp
+#ifdef _OPENMP
tot_em_up_tmp(ks,ithread) = tot_em_up_tmp(ks,ithread) + real(tmp,kind=dp)
+#else
+ tot_em_up(ks) = tot_em_up(ks) + real(tmp,kind=dp)
+#endif
endif ! negative emissions
end do ! nspec
@@ -279,14 +285,13 @@ module emissions_mod
!$OMP END DO
!$OMP END PARALLEL
+ ! manual reduction of em_neg and tot_em_up
#ifdef _OPENMP
do ithread=1,numthreads
em_neg(:,:,:) = em_neg(:,:,:)+em_neg_tmp(:,:,:,ithread)
tot_em_up(:) = tot_em_up(:) + tot_em_up_tmp(:,ithread)
- end do
- deallocate( mass_field_tmp,em_neg_tmp,tot_em_up_tmp )
-#else
- tot_em_up(:) = tot_em_up_tmp(:,1)
+ end do
+ deallocate( mass_field_tmp, em_neg_tmp, tot_em_up_tmp )
#endif
! update for negative emissions
@@ -305,18 +310,18 @@ module emissions_mod
!$OMP END PARALLEL
!! test
- inquire(file=path(2)(1:length(2))//'mass_field.txt',exist=lexist)
- write(*,*) 'emissions: lexist = ',lexist
- if (lexist) then
- open(unittest,file=path(2)(1:length(2))//'mass_field.txt',ACCESS='APPEND',STATUS='OLD')
- else
- open(unittest,file=path(2)(1:length(2))//'mass_field.txt',STATUS='NEW')
- endif
- write(frmt, '(A, I4, A)') '(', nxem, 'E12.3)'
- do jy=1,nyem-1
- write(unittest,frmt) (mass_field(2,ix,jy),ix=1,nxem)
- end do
- close(unittest)
+! inquire(file=path(2)(1:length(2))//'mass_field.txt',exist=lexist)
+! write(*,*) 'emissions: lexist = ',lexist
+! if (lexist) then
+! open(unittest,file=path(2)(1:length(2))//'mass_field.txt',ACCESS='APPEND',STATUS='OLD')
+! else
+! open(unittest,file=path(2)(1:length(2))//'mass_field.txt',STATUS='NEW')
+! endif
+! write(frmt, '(A, I4, A)') '(', nxem, 'E12.3)'
+! do jy=1,nyem-1
+! write(unittest,frmt) (mass_field(2,ix,jy),ix=1,nxem)
+! end do
+! close(unittest)
!!
end subroutine emissions
@@ -572,7 +577,7 @@ module emissions_mod
allocate( em_field(nspec-1,nxem,nyem,2) )
em_field(:,:,:,:)=0.
allocate( em_res(nspec-1,nxem,nyem) )
- if (ipin.eq.1) then
+ if (ipin.eq.2) then
! read residual emissions from end of previous run
write(*,*) 'Reading residual emissions from previous run'
call em_res_read
diff --git a/src/gitversion.txt b/src/gitversion.txt
new file mode 100644
index 00000000..bc95db4e
--- /dev/null
+++ b/src/gitversion.txt
@@ -0,0 +1 @@
+8d9c680 Mon Apr 22 15:37:39 2024 +0200
diff --git a/src/initdomain_mod.f90 b/src/initdomain_mod.f90
index e5a8b710..0ea8502d 100644
--- a/src/initdomain_mod.f90
+++ b/src/initdomain_mod.f90
@@ -536,9 +536,18 @@ module initdomain_mod
endif
endif
- ! Exit here if resuming a run from particle dump
- !***********************************************
- if (gdomainfill.and.ipin.ne.0) return
+ ! If resuming a run from particle dump
+ ! calculate total mass each species then exit
+ !********************************************
+ if (gdomainfill.and.ipin.ne.0) then
+ write(*,*) 'Initialising particles from partoutput'
+ tot_mass(:)=0.
+ do ks=2,nspec
+ tot_mass(ks)=sum(mass(1:count%alive,ks))
+ write(*,'(A,E12.4,A)') 'Species '//species(ks)//': ',tot_mass(ks),' (kg)'
+ end do
+ return
+ endif
! Allocate fields used within this subroutine
!*********************************************
@@ -550,7 +559,6 @@ module initdomain_mod
! Do not release particles twice (i.e., not at both in the leftmost and rightmost
! grid cell) for a global domain
!*****************************************************************************
- write(*,*) 'init_domainfill: nx_we(2), nx-2 = ',nx_we(2), nx-2
if (xglobal) nx_we(2)=min(nx_we(2),nx-2)
write(*,*) 'init_domainfill: nx_we, ny_sn = ',nx_we, ny_sn
@@ -741,10 +749,10 @@ module initdomain_mod
! poles), distribute the particles randomly. When only few particles are
! left distribute them randomly
- if ((ncolumn.gt.20).and.(ncolumn-j.gt.20)) then
+ if ((ncolumn.gt.20)) then !.and.(ncolumn-j.gt.20)) then
pnew=pnew-deltacol
- else if (ncolumn.gt.20) then
- pnew=pnew-ran1(idummy,0)*(pnew-pp(nz))
+! else if (ncolumn.gt.20) then
+! pnew=pnew-ran1(idummy,0)*(pnew-pp(nz))
else
pnew=pp(1)-ran1(idummy,0)*(pp(1)-pp(nz))
endif
@@ -871,6 +879,8 @@ module initdomain_mod
! Assumes lon and lat dimensions are midpoints
ixm=int((xl-(lonini(1,indxn)-0.5*dxini(indxn)))/dxini(indxn))+1
jym=int((yl-(latini(1,indxn)-0.5*dyini(indxn)))/dyini(indxn))+1
+ ixm=min(ixm,nxini(indxn))
+ jym=min(jym,nyini(indxn))
!! testing
if (jj.eq.1.and.jy.lt.5.and.ix.lt.5) then
print*, 'init_domainfill: lonini, xl, latini, yl = ',lonini(ixm,indxn),xl,latini(jym,indxn),yl
@@ -970,9 +980,6 @@ module initdomain_mod
! Total mass each species
!************************
- if (.not.allocated(tot_mass)) then
- allocate( tot_mass(nspec) )
- endif
tot_mass(:)=0.
do ks=2,nspec
tot_mass(ks)=sum(mass(1:count%alive,ks))
diff --git a/src/netcdf_output_mod.f90 b/src/netcdf_output_mod.f90
index e3d36b66..d6be602c 100644
--- a/src/netcdf_output_mod.f90
+++ b/src/netcdf_output_mod.f90
@@ -2436,7 +2436,7 @@ subroutine readpartpositions_netcdf(ibtime,ibdate)
endif
! Open partoutput_end.nc file
- call nf90_err(nf90_open(trim('partoutput_end.nc'), mode=NF90_NOWRITE,ncid=ncidend))
+ call nf90_err(nf90_open(path(2)(1:length(2))//trim('partoutput_end.nc'), mode=NF90_NOWRITE,ncid=ncidend))
! Take the positions of the particles at the last timestep in the file
! It needs to be the same as given in the COMMAND file, this is arbitrary
@@ -2467,6 +2467,9 @@ subroutine readpartpositions_netcdf(ibtime,ibdate)
error stop
endif
+ !! testing
+! print*, 'readpartpositions_netcdf: julin, julcommand = ',julin, julcommand
+
! Then the particle dimension
call nf90_err(nf90_inq_dimid(ncid=ncidend,name='particle',dimid=pIDend))
call nf90_err(nf90_inquire_dimension(ncid=ncidend,dimid=pIDend,len=plen))
@@ -2510,7 +2513,9 @@ subroutine readpartpositions_netcdf(ibtime,ibdate)
do i=1,plen
if (part(i)%z.lt.0) then
call terminate_particle(i,0)
- write(*,*) 'Particle ',i,'is not alive in the restart file.'
+ if (mdomainfill.eq.0) then
+ write(*,*) 'Particle ',i,'is not alive in the restart file.'
+ endif
iterminate=iterminate+1
endif
part(i)%nclass=min(int(ran1(idummy,0)*real(nclassunc))+1, &
@@ -2523,6 +2528,11 @@ subroutine readpartpositions_netcdf(ibtime,ibdate)
call nf90_err(nf90_close(ncidend))
+ !! testing
+! print*, 'readpartpositions_netcdf: number alive = ',count%alive
+! print*, 'readpartpositions_netcdf: range(part%z) = ',minval(part(1:count%alive)%z),maxval(part(1:count%alive)%z)
+! print*, 'readpartpositions_netcdf: part(1)%tstart = ',part(1)%tstart
+
end subroutine readpartpositions_netcdf
subroutine readinitconditions_netcdf()
diff --git a/src/readoptions_mod.f90 b/src/readoptions_mod.f90
index 2f2ab8ed..ff0971ef 100644
--- a/src/readoptions_mod.f90
+++ b/src/readoptions_mod.f90
@@ -1969,7 +1969,7 @@ subroutine readreceptors
lon=-999.9
read(unitreceptor,receptors,iostat=ios)
if ((lon.lt.-900).or.(ios.ne.0)) exit
- if ((time.lt.bdate).or.(time.gt.edate)) cycle ! skip receptors not in simulation window
+ if ((time.lt.bdate).or.(time.ge.edate)) cycle ! skip receptors not in simulation window
j=j+1
end do
numreceptor=j
@@ -1992,7 +1992,7 @@ subroutine readreceptors
lon=-999.9
read(unitreceptor,receptors,iostat=ios)
if ((lon.lt.-900).or.(ios.ne.0)) exit ! read error
- if ((time.lt.bdate).or.(time.gt.edate)) cycle ! skip receptors not in simulation window
+ if ((time.lt.bdate).or.(time.ge.edate)) cycle ! skip receptors not in simulation window
j=j+1
receptorname(j)=receptor
xreceptor(j)=(lon-xlon0)/dx ! transform to grid coordinates
diff --git a/src/receptor_mod.f90 b/src/receptor_mod.f90
index dad99ec7..0883b60b 100644
--- a/src/receptor_mod.f90
+++ b/src/receptor_mod.f90
@@ -133,6 +133,9 @@ module receptor_mod
call receptorcalc(itime,weight,lrecoutstart,lrecoutend,recoutnum,recoutnumsat)
endif
+ !! testing
+ print*, 'receptor_mod: itime, lrecoutstart, lrecoutend = ',itime, lrecoutstart, lrecoutend
+
if ((itime.eq.lrecoutend).and.(any(recoutnum.gt.0.).or.any(recoutnumsat.gt.0.))) then
! output receptor concentrations
@@ -174,8 +177,8 @@ module receptor_mod
endif
!! testing
- print*, 'receptor_mod: concvar_id = ',concvar_id
- print*, 'receptor_mod: satvar_id = ',satvar_id
+! print*, 'receptor_mod: concvar_id = ',concvar_id
+! print*, 'receptor_mod: satvar_id = ',satvar_id
! Initialize variables
@@ -188,7 +191,7 @@ module receptor_mod
#endif
!! testing
- print*, 'receptor_mod: nthread = ',nthreads
+! print*, 'receptor_mod: nthread = ',nthreads
allocate(densityoutrecept(nlayermax))
allocate(nnrec(maxrecsample,nlayermax),&
@@ -222,6 +225,8 @@ module receptor_mod
! Loop over general receptors
!****************************
+
+ write(*,fmt='(A,1X,I8,1X,A,1X,I3)') 'Number of receptors output at itime ',itime,'is',numcurrec
!$OMP PARALLEL &
!$OMP PRIVATE(n,k,nn,nr,ix,jy,ixp,jyp,ddx,ddy,rddx,rddy,p1,p2,p3,p4,indz,indzp,il,dz1,dz2,dz, &
@@ -232,29 +237,22 @@ module receptor_mod
#else
thread=1
#endif
- !! testing
- print*, 'receptor_mod: thread = ',thread
nr=0
!$OMP DO
- do n=1,numreceptor
+ do k=1,numcurrec
- if ((treceptor(n).lt.lrecoutstart).or. &
- (treceptor(n).ge.lrecoutend)) cycle ! skip if not in current sampling time interval
+ if (cpointer(k).eq.0.) cycle
- ! find indice to creceptor, xkreceptor, and nnreceptor
- do k=1,maxrecsample
- if (cpointer(k).eq.n) exit
- end do
+ ! number of the receptor
+ n=cpointer(k)
! counter of receptor values this time interval
nr=nr+1
nr_omp(thread)=nr
!! testing
- if (n.lt.40) then
- print*, 'receptor_mod: n, thread, nr, cpointer(k), rpointer = ',n, thread, nr, cpointer(k), rpointer
- endif
+! print*, 'receptor_mod: n, thread, nr, cpointer(k), rpointer = ',n, thread, nr, cpointer(k), rpointer
if (((.not.llcmoutput).and.(iout.eq.2)).or.&
(llcmoutput.and.((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)))) then
@@ -356,7 +354,7 @@ module receptor_mod
timerec_omp(nr,thread)=treceptor(n)
namerec_omp(nr,thread)=receptorname(n)
- end do ! numreceptor
+ end do ! numcurrec
!$OMP END DO
!$OMP END PARALLEL
@@ -364,7 +362,7 @@ module receptor_mod
nr=1
do ithread=1,nthreads
!! testing
- print*, 'receptor_mod: nr_omp(ithread) = ',nr_omp(ithread)
+! print*, 'receptor_mod: nr_omp(ithread) = ',nr_omp(ithread)
if (nr_omp(ithread).eq.0) cycle
crec(:,nr:(nr+nr_omp(ithread)-1),1)=crec_omp(:,1:nr_omp(ithread),1,ithread)
cunc(:,nr:(nr+nr_omp(ithread)-1),1)=cunc_omp(:,1:nr_omp(ithread),1,ithread)
@@ -377,7 +375,7 @@ module receptor_mod
namerec(nr:(nr+nr_omp(ithread)-1))=namerec_omp(1:nr_omp(ithread),ithread)
nr=nr+nr_omp(ithread)
!! testing
- print*, 'receptor_mod: nr = ',nr
+! print*, 'receptor_mod: nr = ',nr
end do
! total number of receptors to write this interval
@@ -397,8 +395,8 @@ module receptor_mod
rpointer=rpointer+nr
!! testing
- print*, 'receptor_mod: nr_omp = ',nr_omp
- print*, 'receptor_mod: rpointer = ',rpointer
+! print*, 'receptor_mod: nr_omp = ',nr_omp
+! print*, 'receptor_mod: rpointer = ',rpointer
! Loop over satellites
!*********************
@@ -413,6 +411,8 @@ module receptor_mod
cunc_omp(:,:,:,:)=0.
nr_omp(:)=0
+ write(*,fmt='(A,1X,I8,1X,A,1X,I3)') 'Number of satellite receptors output at itime ',itime,'is',numcursat
+
!$OMP PARALLEL &
!$OMP PRIVATE(n,nr,nn,nchar,k,ix,jy,ixp,jyp,ddx,ddy,rddx,rddy,p1,p2,p3,p4,kz,numsatlayer,zmid,indz,indzp, &
!$OMP il,dz1,dz2,dz,ind,rho_p,rhoi,densityoutrecept,ks,thread)
@@ -423,27 +423,21 @@ module receptor_mod
thread=1
#endif
- !! testing
- print*, 'receptor_mod: satellites, thread = ',thread
-
nr=0
!$OMP DO
- do n=1,numsatreceptor
+ do k=1,numcursat
- if ((tsatellite(n).lt.lrecoutstart).or. &
- (tsatellite(n).ge.lrecoutend)) cycle ! skip if not in current sampling time interval
+ if (csatpointer(k).eq.0) cycle
- ! find indice to csatellite, xksatellite, and nnsatellite
- do k=1,maxrecsample
- if (csatpointer(k).eq.n) exit
- end do
+ ! number of satellite receptor
+ n=csatpointer(k)
! counter of receptor values this time interval
nr=nr+1
nr_omp(thread)=nr
!! testing
- print*, 'receptor_mod: n, thread, nr, csatpointer(k) = ',n, thread, nr, csatpointer(k)
+! print*, 'receptor_mod: n, thread, nr, csatpointer(k) = ',n, thread, nr, csatpointer(k)
! get actual number vertical layers for this retrieval
do nn=1,numsatellite
@@ -544,7 +538,7 @@ module receptor_mod
timerec_omp(nr,thread)=tsatellite(n)
namerec_omp(nr,thread)=satellitename(n)
- end do ! numsatreceptor
+ end do ! numcursat
!$OMP END DO
!$OMP END PARALLEL
@@ -552,7 +546,7 @@ module receptor_mod
nr=1
do ithread=1,nthreads
!! testing
- print*, 'receptor_mod: satellites, nr_omp(ithread) = ',nr_omp(ithread)
+! print*, 'receptor_mod: satellites, nr_omp(ithread) = ',nr_omp(ithread)
if (nr_omp(ithread).eq.0) cycle
crec(:,nr:nr+nr_omp(ithread)-1,:)=crec_omp(:,1:nr_omp(ithread),:,ithread)
cunc(:,nr:nr+nr_omp(ithread)-1,:)=cunc_omp(:,1:nr_omp(ithread),:,ithread)
@@ -565,7 +559,7 @@ module receptor_mod
namerec(nr:nr+nr_omp(ithread)-1)=namerec_omp(1:nr_omp(ithread),ithread)
nr=nr+nr_omp(ithread)
!! testing
- print*, 'receptor_mod: satellites, nr = ',nr
+! print*, 'receptor_mod: satellites, nr = ',nr
end do
! total number of receptors to write this interval
@@ -585,8 +579,8 @@ module receptor_mod
spointer=spointer+nr
!! testing
- print*, 'receptor_mod: satellites, nr_omp = ',nr_omp
- print*, 'receptor_mod: spointer = ',spointer
+! print*, 'receptor_mod: satellites, nr_omp = ',nr_omp
+! print*, 'receptor_mod: spointer = ',spointer
! close files
if (lnetcdfout.eq.1) then
@@ -600,9 +594,6 @@ module receptor_mod
#endif
endif
- ! Reinitialization
- !*****************
-
deallocate(densityoutrecept)
deallocate(crec,cunc,nnrec,xkrec)
deallocate(crec_omp,cunc_omp,nnrec_omp,xkrec_omp)
@@ -610,17 +601,6 @@ module receptor_mod
deallocate(lonrec_omp,latrec_omp,altrec_omp,timerec_omp,namerec_omp)
deallocate(nr_omp)
- creceptor(:,:)=0.
- crecuncert(:,:)=0.
- nnreceptor(:)=0.
- xkreceptor(:)=0.
- if (numsatreceptor.gt.0) then
- csatellite(:,:,:)=0.
- csatuncert(:,:,:)=0.
- nnsatellite(:,:)=0.
- xksatellite(:,:)=0.
- endif
-
! End of receptor output
!***********************
@@ -629,6 +609,27 @@ module receptor_mod
endif ! output receptor conc
+ endif ! receptor output
+
+ if (itime.eq.lrecoutend) then
+
+ ! Reinitialize
+ !*************
+
+ creceptor(:,:)=0.
+ crecuncert(:,:)=0.
+ nnreceptor(:)=0.
+ xkreceptor(:)=0.
+ if (numsatreceptor.gt.0) then
+ csatellite(:,:,:)=0.
+ csatuncert(:,:,:)=0.
+ nnsatellite(:,:)=0.
+ xksatellite(:,:)=0.
+ endif
+
+ recoutnum(:)=0.
+ recoutnumsat(:,:)=0.
+
! Update output timesteps
!************************
lrecoutnext=lrecoutnext+lrecoutstep
@@ -640,7 +641,7 @@ module receptor_mod
call receptorcalc(itime,weight,lrecoutstart,lrecoutend,recoutnum,recoutnumsat)
endif
- endif ! calculate receptors
+ endif
end subroutine receptoroutput
@@ -711,13 +712,14 @@ module receptor_mod
!********************
! pointer for creceptor, xkreceptor and nnreceptor
+ numcurrec=0
cpointer(:)=0
k=0
do n=1,numreceptor
if ((treceptor(n).lt.lrecoutstart).or. &
- (treceptor(n).gt.lrecoutend)) cycle ! skip if not in current sampling time interval
+ (treceptor(n).ge.lrecoutend)) cycle ! skip if not in current sampling time interval
! update pointer
k=k+1
@@ -727,8 +729,8 @@ module receptor_mod
endif
cpointer(k)=n
- !! test
- print*, 'receptorcalc: n, receptorname(n) = ',n, receptorname(n)
+ !! testing
+! print*, 'receptorcalc: n, receptorname(n) = ',n, receptorname(n)
! Reset concentrations for new receptor
conc(:)=0.
@@ -746,14 +748,12 @@ module receptor_mod
endif
h=hxmax*hymax*hzmax
- !! test
+ !! testing
! print*, 'receptorcalc: rec_ff = ',rec_ff
! Estimate concentration at receptor
!***********************************
-!! TO DO: replace ATOMIC with j variable with dim of thread
-
!$OMP PARALLEL &
!$OMP PRIVATE(i,itage,hz,zd,hx,xd,hy,yd,h,r2,xkern,ks) &
!$OMP SHARED(j,unc,conc) &
@@ -816,7 +816,7 @@ module receptor_mod
do ks=ks_start,nspec
if (llcmoutput) then
- ! special case CTM output use mass ratio species to airtracer
+ ! special case LCM output use mass ratio species to airtracer
! species 1 is always airtracer
conc(ks)=conc(ks) + mass(i,ks)/mass(i,1) * &
weight * xkern
@@ -852,7 +852,10 @@ module receptor_mod
recoutnum(k)=recoutnum(k)+weight
endif
- !! test
+ ! update number receptors this time interval
+ numcurrec=k
+
+ !! testing
! print*, 'receptorcalc: j, conc, xksum = ',j, conc(2), xksum
! print*, 'receptorcalc: n, k, cpointer(k) = ',n, k, cpointer(k)
! print*, 'receptorcalc: nnreceptor(k) = ',nnreceptor(k)
@@ -866,15 +869,16 @@ module receptor_mod
!*********************
! pointer for csatellite, xksatellite and nnsatellite
+ numcursat=0
csatpointer(:)=0
k=0
do n=1,numsatreceptor
if ((tsatellite(n).lt.lrecoutstart).or. &
- (tsatellite(n).gt.lrecoutend)) cycle ! skip if not in current sampling time interval
+ (tsatellite(n).ge.lrecoutend)) cycle ! skip if not in current sampling time interval
- !! test
+ !! testing
! print*, 'receptorcalc: lrecoutstart, lrecoutend, tsatellite(n) = ',lrecoutstart, lrecoutend, tsatellite(n)
! update pointer
@@ -894,7 +898,7 @@ module receptor_mod
endif
end do
- !! test
+ !! testing
! print*, 'receptorcalc: n, satellitename(n) = ',n, satellitename(n)
do kz=1,numsatlayer
@@ -921,7 +925,7 @@ module receptor_mod
h=hxsat*hysat*(0.5*hz)
- !! test
+ !! testing
! print*, 'zmid, hz, eta =',zmid, hz, eta
!$OMP PARALLEL &
@@ -959,7 +963,7 @@ module receptor_mod
do ks=ks_start,nspec
if (llcmoutput) then
- ! special case CTM output use mass ratio species to airtracer
+ ! special case LCM output use mass ratio species to airtracer
! species 1 is always airtracer
conc(ks)=conc(ks) + mass(i,ks)/mass(i,1) * &
weight * xkern
@@ -995,7 +999,7 @@ module receptor_mod
recoutnumsat(kz,k)=recoutnumsat(kz,k)+weight
endif
- !! test
+ !! testing
! print*, 'receptorcalc: for satellite: j, conc, xksum = ',j, conc(2), xksum
! print*, 'receptorcalc: n, k, csatpointer(k) = ',n, k, csatpointer(k)
! print*, 'receptorcalc: nnsatellite(:,k) = ',nnsatellite(:,k)
@@ -1003,6 +1007,9 @@ module receptor_mod
end do ! numsatlayer
+ ! update current number of satellite receptors this time interval
+ numcursat=k
+
end do ! numsatreceptor
end subroutine receptorcalc
diff --git a/src/timemanager_mod.f90 b/src/timemanager_mod.f90
index 0871a411..931477c1 100644
--- a/src/timemanager_mod.f90
+++ b/src/timemanager_mod.f90
@@ -304,6 +304,19 @@ subroutine timemanager
else
call init_domainfill
endif
+ if (ipin.eq.2) then
+ ! Particles initialized from partoutput
+#ifdef ETA
+!$OMP PARALLEL PRIVATE(i,j)
+!$OMP DO
+ do i=1,count%alive
+ j=count%ialive(i)
+ call update_z_to_zeta(itime,j)
+ end do
+!$OMP END DO
+!$OMP END PARALLEL
+#endif
+ endif
else
call boundcond_domainfill(itime,loutend)
endif
@@ -704,17 +717,8 @@ subroutine timemanager
#ifdef USE_NCF
if ((mdomainfill.eq.1).and.(llcmoutput)) then
-
- if (.not.allocated(tot_mass)) then
- allocate( tot_mass(nspec) )
- tot_mass(:)=0.
- endif
-
do ks=1,nspec
- do i=1,count%alive
- j=count%ialive(i)
- tot_mass(ks)=tot_mass(ks)+real(mass(j,ks),kind=dp)
- end do
+ tot_mass(ks)=sum(real(mass(1:count%alive,ks),kind=dp))
end do
call totals_write(itime)
endif
diff --git a/src/totals_mod.f90 b/src/totals_mod.f90
index 912593cd..65491241 100644
--- a/src/totals_mod.f90
+++ b/src/totals_mod.f90
@@ -44,6 +44,7 @@ module totals_mod
if (stat.ne.0) error stop "Could not allocate totals arrays"
chem_loss(:,:)=0.
+ allocate( tot_mass(nspec) )
allocate( tot_em_up(nspec) )
allocate( tot_em_field(nspec) )
allocate( tot_em_res(nspec) )
--
GitLab