diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index ec735eb06f5f63bdc451d7e6fc66d9db4f37d2d7..bab4c82371d6ce5000ced36913ac233f3f04a1df 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="2464aef Fri Jun 21 12:24:26 2024 +0200"
! Keeping track of the total running time of FLEXPART, printed out at the end.
!*****************************************************************************
@@ -180,6 +180,7 @@ subroutine read_options_and_initialise_flexpart
use receptor_netcdf_mod, only: read_satellite_info, receptorout_init
#endif
use receptor_mod, only: alloc_receptor
+ use nudging_mod, only: nudging_init
implicit none
@@ -323,6 +324,13 @@ subroutine read_options_and_initialise_flexpart
!*************************************************************************
call initialise_particles
+ ! Read nudging settings (after ini particles because
+ ! it needs species information)
+ !*************************************************************************
+ if(lnudging) then
+ call nudging_init
+ endif
+
! Initialize variables for totals calculations
!*********************************************
#ifdef USE_NCF
diff --git a/src/chemistry_mod.f90 b/src/chemistry_mod.f90
index 4994cd1ea994f3ef9ad70398851bae3d52f88d29..3156b3c7be3e313579da7e4be4d1b8659b1c21a6 100644
--- a/src/chemistry_mod.f90
+++ b/src/chemistry_mod.f90
@@ -60,6 +60,7 @@ module chemistry_mod
implicit none
character(len=256) :: preag_path
+ character(len=256) :: preag_file
character(len=16) :: preagent, punit
integer :: phourly
integer,parameter :: unitreagents=2, unitreagentsout=3
@@ -68,10 +69,11 @@ module chemistry_mod
! declare namelist
namelist /reagent_params/ &
- preagent, preag_path, phourly, punit
+ preagent, preag_path, preag_file,phourly, punit
preagent="" ! read failure indicator value
preag_path=""
+ preag_file=""
punit=""
phourly=0
reag_hourly(:)=0
@@ -115,6 +117,7 @@ module chemistry_mod
read(unitreagents,reagent_params,iostat=readerror)
reagents(j)=preagent
reag_path(j)=preag_path
+ reag_file(j)=preag_file
reag_hourly(j)=phourly
reag_unit(j)=punit
! namelist output
@@ -208,7 +211,9 @@ module chemistry_mod
! Read new chem field and store in 2nd position
!**********************************************
- write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
+! write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
+ write(CL_name,'(A)') trim(reag_path(nr))//trim(reag_file(nr))//'_'//amonth//'.nc'
+
inquire(file=CL_name,exist=lexist)
if (lexist) then
write(*,*) 'Reading CL field: ',CL_name
@@ -282,7 +287,9 @@ module chemistry_mod
! Read new chem field
!********************
- write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
+ !write(CL_name,'(A)') trim(reag_path(nr))//trim(reagents(nr))//'_'//amonth//'.nc'
+ write(CL_name,'(A)') trim(reag_path(nr))//trim(reag_file(nr))//'_'//amonth//'.nc'
+
inquire(file=CL_name,exist=lexist)
if (lexist) then
write(*,*) 'Reading CL field: ',CL_name
@@ -483,9 +490,9 @@ module chemistry_mod
!$OMP END PARALLEL
jrate_average=jrate_average/24.
end if
- !! testing
-! print*, 'getchemhourly: range(jrate_average) = ',minval(jrate_average),maxval(jrate_average)
+ print*, 'getchemhourly: range(jrate_average) = ',minval(jrate_average),maxval(jrate_average)
endif
+ !! testing
! loop over reagents
do nr=1,nreagent
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 9e87f459e04518081c5c26f381d6151a950e8086..ad7fe4d301d029a0ae3a8c0e146a1847abda736f 100644
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -103,7 +103,7 @@ module com_mod
integer :: surf_only ! deprecated
logical :: turbswitch
integer :: cblflag !added by mc for cbl
- logical :: llcmoutput
+ logical :: llcmoutput, lnudging
! ctl factor, by which time step must be smaller than Lagrangian time scale
! ifine reduction factor for time step used for vertical wind
@@ -229,6 +229,7 @@ module com_mod
integer,allocatable,dimension(:) :: ishape,orient
! chemical reagent variables
character(len=256) :: reag_path(maxreagent)
+ character(len=256) :: reag_file(maxreagent)
character(len=16) :: reagents(maxreagent), reag_unit(maxreagent)
integer :: reag_hourly(maxreagent), nreagent
! reaction rates
diff --git a/src/makefile_gfortran b/src/makefile_gfortran
index 6b0c04d73d4a7ad380e4f7add40fc3f0fe93fb7f..99fadd9a4f635efe140618657122f43bd8ebf0ef 100644
--- a/src/makefile_gfortran
+++ b/src/makefile_gfortran
@@ -103,6 +103,7 @@ LIBS = -leccodes -leccodes_f90 -lm $(NCOPT) $(ETAOPT) #MD -ljasper
FFLAGS = $(INC) -O$(O_LEV) -m64 -cpp -mcmodel=large $(NCOPT) $(ETAOPT) $(FUSER) #-convert little_endian -fpp -WB -check all
+#DBGFLAGS = $(INC) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=large -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) $(NCOPT) $(ETAOPT) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
DBGFLAGS = $(INC) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=large -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) $(NCOPT) $(ETAOPT) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
LDFLAGS = $(FFLAGS) $(LIBPATH) $(LIBRPATH) $(LIBS)
@@ -134,6 +135,7 @@ OBJECTS_NCF = \
netcdf_output_mod.o chemistry_mod.o \
totals_mod.o initdomain_mod.o \
receptor_netcdf_mod.o emissions_mod.o \
+nudging_mod.o
OBJECTS_ETA = coord_ecmwf_mod.o
@@ -293,7 +295,7 @@ receptor_mod.o: par_mod.o com_mod.o point_mod.o particle_mod.o date_mod.o \
endif
ifneq ($(ncf), no)
netcdf_output_mod.o: mean_mod.o outgrid_mod.o readoptions_mod.o drydepo_mod.o
-totals_mod.o: par_mod.o com_mod.o netcdf_output_mod.o
+totals_mod.o: par_mod.o com_mod.o netcdf_output_mod.o nudging_mod.o
chemistry_mod.o: par_mod.o com_mod.o date_mod.o particle_mod.o \
point_mod.o windfields_mod.o totals_mod.o netcdf_output_mod.o
initdomain_mod.o: par_mod.o com_mod.o point_mod.o random_mod.o outgrid_mod.o \
@@ -302,8 +304,11 @@ receptor_netcdf_mod.o: par_mod.o com_mod.o point_mod.o date_mod.o \
windfields_mod.o netcdf_output_mod.o
emissions_mod.o: par_mod.o com_mod.o point_mod.o particle_mod.o date_mod.o \
netcdf_output_mod.o totals_mod.o windfields_mod.o
+nudging_mod.o: par_mod.o com_mod.o point_mod.o readoptions_mod.o
+
endif
+
#8)
conv_mod.o: flux_mod.o class_gribfile_mod.o qvsat_mod.o sort_mod.o
ifeq ($(ncf), no)
@@ -329,7 +334,7 @@ initialise_mod.o: interpol_mod.o xmass_mod.o output_mod.o turbulence_mod.o
ifneq ($(ncf), no)
timemanager_mod.o: advance_mod.o conv_mod.o plume_mod.o getfields_mod.o restart_mod.o \
initialise_mod.o chemistry_mod.o initdomain_mod.o receptor_mod.o \
- emissions_mod.o totals_mod.o
+ emissions_mod.o totals_mod.o nudging_mod.o
else
timemanager_mod.o: advance_mod.o conv_mod.o plume_mod.o getfields_mod.o restart_mod.o \
initialise_mod.o receptor_mod.o
diff --git a/src/readoptions_mod.f90 b/src/readoptions_mod.f90
index c214b2acdfdcb0daecf0403a2fe19a4433d08b45..86ffe17fce20c3d80216613fb9dc0e84b0de3f26 100644
--- a/src/readoptions_mod.f90
+++ b/src/readoptions_mod.f90
@@ -482,7 +482,7 @@ subroutine readcommand
character(len=50) :: line
integer :: ios
- integer :: lturbulence_meso,lcmoutput
+ integer :: lturbulence_meso,lcmoutput, nudging
character(len=50) :: ohfields_path ! deprecated
namelist /command/ &
@@ -529,6 +529,7 @@ subroutine readcommand
logvertinterp, &
ohfields_path, &
lcmoutput, &
+ nudging, &
itsplit ! deprecated: only for IO back compatibility
! Presetting namelist command
@@ -575,13 +576,13 @@ subroutine readcommand
logvertinterp=0
ohfields_path=''
lcmoutput=0
+ nudging=0
itsplit=999999999 ! deprecated: only for IO back compatibility
!Af set release-switch
WETBKDEP=.false.
DRYBKDEP=.false.
-
! Open the command file and read user options
! Namelist input first: try to read as namelist file
!**************************************************************************
@@ -1167,6 +1168,33 @@ subroutine readcommand
write(*,*) 'Switch for LCM output LCMOUTPUT = ',llcmoutput
+ !Switch for nudging
+ !*******************************************************************
+ !CGZ implementation for domain-filling mode and LCM output only
+ if (nudging.eq.1) then
+ if ((mdomainfill.ne.1) )then
+ write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
+ write(*,*) '#### NUDGING ONLY FOR DOMAIN FILLING MODE ####'
+ write(*,*) '#### Set MDOMAINFILL=1. ####'
+ error stop
+ endif
+ endif
+
+ if (nudging.eq.1) then
+ if ((.not.llcmoutput) )then
+ write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
+ write(*,*) '#### NUDGING ONLY WITH LCM OUTPUT ####'
+ error stop
+ endif
+ endif
+
+ if (nudging.eq.0) then
+ lnudging=.false.
+ else
+ lnudging=.true.
+ endif
+ !*******************************************************************
+
! Compute modeling time in seconds and beginning date in Julian date
!*******************************************************************
@@ -2944,7 +2972,6 @@ subroutine readspecies(id_spec,pos_spec)
Fn(pos_spec)=f*f*e*((dquer(pos_spec))**3)/(psa*pia*pla) ! Newton's regime
Fs(pos_spec)=f*e**(1.3)*(dquer(pos_spec)**3/(psa*pia*pla)) ! Stokes' regime
endif
-
! Pre-compute ks and kn values needed for horizontal and average orientation (B&B Figure 12 k_(s,max))
ks1(pos_spec)=(Fs(pos_spec)**(1./3.) + Fs(pos_spec)**(-1./3.))/2.
ks2(pos_spec)=0.5*((Fs(pos_spec)**0.05)+(Fs(pos_spec)**(-0.36)))
diff --git a/src/timemanager_mod.f90 b/src/timemanager_mod.f90
index e7666b688d1c379777df04e925a574a52e7af7ea..de61bae6af248cd40102607065dd0feee5534c4a 100644
--- a/src/timemanager_mod.f90
+++ b/src/timemanager_mod.f90
@@ -111,6 +111,7 @@ subroutine timemanager
use receptor_netcdf_mod, only: readreceptors_satellite, verttransform_satellite
use emissions_mod
use totals_mod
+ use nudging_mod, only: nudging_individual
#endif
use receptor_mod, only: receptoroutput
@@ -526,6 +527,17 @@ subroutine timemanager
!$OMP END DO
!$OMP END PARALLEL
+
+ !Nudge fields to observations
+ !**************************************************
+ if (lnudging .and. itime.gt.2*lsynctime) then
+ call nudging_individual(itime)
+ end if
+ !***************************************
+
+
+
+
#ifdef _OPENMP
call omp_set_num_threads(numthreads_grid)
#endif
@@ -744,7 +756,9 @@ subroutine timemanager
! De-allocate memory and end
!***************************
- call dealloc_all_particles
+ !Switch off as test, function does not work if no particles were terminated during simulation!!!!
+ write(*,*) 'WARNING; switched off dealloc all particles!!!!!!!!!!!!!!!!!!!!!!!!!'
+ !call dealloc_all_particles
call dealloc_windfields
call dealloc_domainfill
call dealloc_drydepo
diff --git a/src/totals_mod.f90 b/src/totals_mod.f90
index bf6af51f64ebcbc7969c6119a3b3b808ef5ac164..69fc521b4a3016f80d4a6580a8e05c5eb175bfef 100644
--- a/src/totals_mod.f90
+++ b/src/totals_mod.f90
@@ -14,12 +14,14 @@ module totals_mod
use par_mod, only: dp
use com_mod
use netcdf_output_mod, only: nf90_err
+ use nudging_mod, only: d_nudge
implicit none
character(len=256) :: fn_totals
integer :: nc_id, specdim_id, reagdim_id, timedim_id, nchardim_id
integer :: time_id, spec_id, cl_id, emis_id, efld_id, eres_id, mass_id
+ integer :: nudge_id
real(kind=dp), dimension(:,:), allocatable :: chem_loss
real(kind=dp), dimension(:), allocatable :: tot_mass
real(kind=dp), dimension(:), allocatable :: tot_em_up
@@ -128,6 +130,14 @@ module totals_mod
call nf90_err( nf90_put_att(nc_id, cl_id, 'units', 'kg s-1') )
call nf90_err( nf90_put_att(nc_id, cl_id, 'long_name', 'Mass loss through chemical reactions') )
+ if (lnudging) then
+ ! nudging
+ call nf90_err( nf90_def_var(nc_id, 'nudging', nf90_double, (/ specdim_id, timedim_id /), nudge_id) )
+ call nf90_err( nf90_put_att(nc_id, cl_id, 'units', 'kg s-1') )
+ call nf90_err( nf90_put_att(nc_id, cl_id, 'long_name', 'Mass flux due to nudging') )
+ end if
+
+
! write global attributes
!************************
@@ -172,6 +182,7 @@ module totals_mod
integer :: var_id
! open file
+ print*, 'cgz db open file totals: ', trim(fn_totals)
call nf90_err( nf90_open(trim(fn_totals), nf90_write, nc_id) )
! get length of time dimension -> increase index by one to write new data
@@ -208,6 +219,13 @@ module totals_mod
(/ 1, 1, tidx /), (/ nreagent, nspec, 1/)) )
endif
+ if (lnudging) then
+ call nf90_err( nf90_inq_varid(nc_id, "nudging", var_id) )
+ call nf90_err( nf90_put_var(nc_id, var_id, d_nudge(1:nspec)/real(lsynctime), &
+ (/ 1, tidx /), (/ nspec, 1/)) )
+ end if
+
+
call nf90_err( nf90_close(nc_id) )
end subroutine totals_write