diff --git a/prep_flexpart/main.f90 b/prep_flexpart/main.f90
index 475d20ace5a31e592c58f81f485185b8709a9d7d..a975fb6910413d7a5dc02aa5c2993f00d07b9d4d 100644
--- a/prep_flexpart/main.f90
+++ b/prep_flexpart/main.f90
@@ -127,7 +127,11 @@ program main
call prep_outgrid(settings, jd, nr)
! write AGECLASS
- call prep_ageclass(settings, jd, nr)
+ if ( settings%lFPversion.eq.0 ) then
+ call prep_ageclass_FP10(settings, jd, nr)
+ else
+ call prep_ageclass(settings, jd, nr)
+ endif
! read and process observations
if ( settings%obsformat.eq.'obspack' ) then
diff --git a/prep_flexpart/makefile b/prep_flexpart/makefile
index cabcc3ca5da1f447b65f3b655f9b17ae145d10fb..48216e5ca708606f422da1bdb663edcaefc23a3b 100644
--- a/prep_flexpart/makefile
+++ b/prep_flexpart/makefile
@@ -23,6 +23,7 @@ SRCS = mod_var.f90 \
prep_command.f90 \
prep_outgrid.f90 \
prep_ageclass.f90 \
+ prep_ageclass_FP10.f90 \
prep_releases.f90 \
prep_releases_reg.f90 \
process_obs.f90 \
diff --git a/prep_flexpart/mod_settings.f90 b/prep_flexpart/mod_settings.f90
index 5b740c57ed8dbfbe40f1ae5d69f87f17b51134ba..924313e172b5917adc15a4b46b9c0afd67669a89 100644
--- a/prep_flexpart/mod_settings.f90
+++ b/prep_flexpart/mod_settings.f90
@@ -49,6 +49,7 @@ module mod_settings
character(len=max_name_len) :: obsformat
character(len=max_name_len) :: suffix
+ integer :: lFPversion
integer :: lselect
integer :: lrelease
real :: relfreq
@@ -350,6 +351,9 @@ module mod_settings
if ( match ) settings%windfield = cc
! General settings
+ identifier = "lFPversion:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%lselect = int(cn)
identifier = "lselect:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lselect = int(cn)
diff --git a/prep_flexpart/prep_ageclass.f90 b/prep_flexpart/prep_ageclass.f90
index 6bd7f22285a43324a1bba4a532b5a81e2a5865ae..3813158000cc0fd1e309ce7563b0abaa5a344d69 100644
--- a/prep_flexpart/prep_ageclass.f90
+++ b/prep_flexpart/prep_ageclass.f90
@@ -52,9 +52,14 @@ subroutine prep_ageclass(settings, jd, nr)
integer :: ierr
integer :: nageclass, lage
+
+ ! FP 11
namelist /ageclass/ &
- nageclass, &
- lage
+ lage
+
+ namelist /nage/ &
+ nageclass
+
! preset namelist variables
nageclass = 1
@@ -67,8 +72,9 @@ subroutine prep_ageclass(settings, jd, nr)
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
if( settings%lnamelist.eq.1 ) then
- ! use namelist file format
+ ! use namelist file format, FP11
write(100,nml=ageclass)
+ write(100,nml=nage)
else
! use old file format
write(100,fmt='(A48)') '************************************************'
diff --git a/prep_regions/makefile b/prep_regions/makefile
index eec17b9bd3c83cf85956beb6b473cc29cb8b553d..d9e772a5c3f397aa318cb1f2e44e33d6356873e1 100644
--- a/prep_regions/makefile
+++ b/prep_regions/makefile
@@ -1,37 +1,44 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff -llapack
-#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form
-FFLAGS = -O0 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
- -fbacktrace
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -lm -DUSE_NCF -lnetcdff
+
+FFLAGS = -I$(INCPATH1) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) -DUSE_NCF -lnetcdff -fbacktrace -Warray-bounds -fcheck=all # -march=native
+
+LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS)
MAIN = prep_regions
-SRCS = mod_var.f90 \
- mod_settings.f90 \
- mod_strings.f90 \
- mod_dates.f90 \
- mod_ncdf.f90 \
- mod_flexpart.f90 \
- sort.f90 \
- read_reclist.f90 \
- initialize.f90 \
- read_lsm.f90 \
- read_obs.f90 \
- get_surfinf.f90 \
- map_grid.f90 \
- prep_regions.f90
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_settings.o \
+ mod_strings.o \
+ mod_dates.o \
+ mod_ncdf.o \
+ mod_flexpart.o
+
+OBJECTS = sort.o \
+ read_reclist.o \
+ initialize.o \
+ read_lsm.o \
+ read_obs.o \
+ get_surfinf.o \
+ map_grid.o \
+ prep_regions.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
+ rm -f *.o *.mod
+
+.SUFFIXES = $(SUFFIXES) .f90
diff --git a/source/error_cov.f90 b/source/error_cov.f90
index 62c201430f931ae279eee94466aab4e0be522072..2ab3dc6b59642318dd3ad438c79e20268fc84f59 100644
--- a/source/error_cov.f90
+++ b/source/error_cov.f90
@@ -74,6 +74,7 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
real(kind=8), dimension(:,:), allocatable :: evecs
integer :: ierr, info, neval, ind
real, parameter :: smallnum=1.e-15
+ integer :: timei, timef, countrate
! calculate error covariance
! --------------------------
@@ -106,6 +107,13 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
! calculate eigendecomposition of cov
! -----------------------------------
+ ! call to dsyevr to find only ILth to IUth eigenvalues and vectors
+ ! thus reducing time needed for eigen decomposition for large B
+ ! IL = index of smallest eigenvalue to be found
+ ! IU = index of largest eigenvalue to be found
+ ! where 1 <= IL <= IU <=N
+ ! IU = N and IL = max(1, N - f*N) where f is fraction of eigenvalues to keep
+
! B = U*D*U^(T)
lwork = 26*ndvar
liwork = 10*ndvar
@@ -121,24 +129,73 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
allocate( zwork(ndvar,ndvar), stat=ierr )
if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ zwork = cov
+
+ ! arguments to dsyevr
+ ! jobsz = N : N only eigenvalues, V both
+ ! range = A : A all eigenvalues, V values in interval (VL,VU], I from ILth to IUth values
+ ! uplo = L : lower triangle of matrix A stored
+ ! n = ndvar : order of matrix A
+ ! matrix = zwork : matrix A
+ ! lda = ndvar : leading dimension of A
+ ! vl = used if range is V : lower bound of interval to search for eigenvalues
+ ! vu = used if range is V : upper bound of interval to search for eigenvalues
+ ! il = used if range is I : index of smallest eigenvalue to be found
+ ! iu = used if range is I : index of largest eigenvalue to be found
+ ! absol = 0 : absolute tolerance for eigenvalues
+ ! m = neval : number of eigenvalues found
+ ! w = evals : eigenvalues in ascending order
+ ! z = evecs : found eigenvectors as columns
+ ! ldz = ndvar : leading dimension of z
+ ! isuppz = isuppz : matrix of indices to non-zero values in z
+ ! work = temporary working array
+ ! lwork = dimension of work
+ ! iwork = temporary working array
+ ! liwork = dimension of iwork
+ ! info : 0 success, < 0 if info = -i ith argument had an illegal value, > 0 internal error
+
+ ! first calculate only eigenvalues
+ call system_clock(timei, countrate)
+ call dsyevr('N','A','L',ndvar,zwork,ndvar,1.,1.,1,1,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
+ if ( info.ne.0 ) then
+ write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
+ stop
+ endif
+ call system_clock(timef, countrate)
+ print*, 'error_cov: time to calculate eigenvalues (s) = ',(timef-timei)/countrate
+ ! now calculate eigenvalues in range IU = N and IL = max(1, N - f*N)
+ ! where f is fraction of eigenvalues to keep
+ ! and their corresponding eigenvectors
+ ind = minloc(evals, dim=1, mask=evals.ge.(config%trunc*maxval(evals)))
+ print*, 'error_cov: ind, evals(ind) = ',ind, evals(ind)
+ call system_clock(timei, countrate)
+ ! because on exit the lower triangle (if UPLO='L') or the upper triangle
+ ! (if UPLO='U') of A, including the diagonal, is destroyed
+ ! need to reassign matrix A
zwork = cov
- call dsyevr('V','A','L',ndvar,zwork,ndvar,1.,1.,1,1,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
+ evals(:) = 0.
+ evecs(:,:) = 0.
+ call dsyevr('V','I','L',ndvar,zwork,ndvar,1.,1.,ind,ndvar,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
if ( info.ne.0 ) then
write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
stop
endif
+ call system_clock(timef, countrate)
+ print*, 'error_cov: time to calculate eigenvalues and eigenvectors (s) = ',(timef-timei)/countrate
+ print*, 'error_cov: neval = ',neval
! truncate eigenvalues
! eigenvalues are sorted in ascending order
- ind = minloc(evals, dim=1, mask=evals.ge.(config%trunc*maxval(evals)))
- write(logid,*) 'Truncating eigenvalues at: ',config%trunc*maxval(evals)
- neig = ndvar - ind + 1
+ neig = neval
write(logid,*) 'Number of retained eigenvalues: ',neig
call alloc_covar(covar)
- covar%evals = evals(ind:ndvar)
- covar%evecs = evecs(:,ind:ndvar)
+ ! eigenvalues in ascending order
+ ! covar%evals = evals(ind:ndvar)
+ covar%evals = evals(1:neig)
+ ! calculated eigenvectors in columns 1:neig
+ covar%evecs = evecs(:,1:neig)
deallocate(work)
deallocate(zwork)
diff --git a/source/makefile b/source/makefile
index 9e9e29fe8e762d223dcc280ed33190f8910c0150..4f2c631185349124d609c532fc833bb0b473e394 100644
--- a/source/makefile
+++ b/source/makefile
@@ -1,63 +1,68 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff -llapack
-#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form -ffree-line-length-0
-FFLAGS = -O0 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \
- -fbacktrace -ffree-line-length-0 -fcheck=all -Wall
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -llapack -lm -DUSE_NCF -lnetcdff
+
+FFLAGS = -I$(INCPATH1) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) -DUSE_NCF -lnetcdff -fbacktrace -Warray-bounds -fcheck=all -ffree-form -ffree-line-length-0 # -march=native
+
+LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS)
MAIN = main
-SRCS = mod_var.f90 \
- mod_settings.f90 \
- mod_dates.f90 \
- mod_strings.f90 \
- mod_tools.f90 \
- mod_fluxes.f90 \
- mod_obs.f90 \
- mod_states.f90 \
- mod_covar.f90 \
- mod_save.f90 \
- mod_flexpart.f90 \
- mod_transform.f90 \
- mod_perturb.f90 \
- read_reclist.f90 \
- read_ncdf.f90 \
- read_obs.f90 \
- read_bg.f90 \
- initialize.f90 \
- correl.f90 \
- error_cov.f90 \
- read_lsm.f90 \
- read_orog.f90 \
- mod_regions.f90 \
- interpz.f90 \
- convertgrid.f90 \
- get_initconc.f90 \
- init_cini.f90 \
- init_cini_month.f90 \
- init_co2.f90 \
- init_ghg.f90 \
- calc_conc.f90 \
- simulate.f90 \
- average_fp.f90 \
- congrad.f90 \
- mod_analytic.f90 \
- mod_m1qn3.f90 \
- m1qn3_interface.f90 \
- retrieval.f90 \
- main.f90
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_settings.o \
+ mod_dates.o \
+ mod_strings.o \
+ mod_tools.o \
+ mod_fluxes.o \
+ mod_obs.o \
+ mod_states.o \
+ mod_covar.o \
+ mod_save.o \
+ mod_flexpart.o \
+ mod_transform.o \
+ mod_regions.o \
+ mod_analytic.o \
+ mod_m1qn3.o \
+ mod_perturb.o
+
+OBJECTS = read_reclist.o \
+ read_ncdf.o \
+ read_obs.o \
+ read_bg.o \
+ initialize.o \
+ correl.o \
+ error_cov.o \
+ read_lsm.o \
+ read_orog.o \
+ interpz.o \
+ convertgrid.o \
+ get_initconc.o \
+ init_cini.o \
+ init_cini_month.o \
+ init_co2.o \
+ init_ghg.o \
+ calc_conc.o \
+ simulate.o \
+ average_fp.o \
+ congrad.o \
+ m1qn3_interface.o \
+ retrieval.o \
+ main.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
-
+ rm -f *.o *.mod
+.SUFFIXES = $(SUFFIXES) .f90
diff --git a/source/read_obs.f90 b/source/read_obs.f90
index a2377047c6733c91442c50dca8b2eb94b9484792..d987322ac3b4ebb4029352e7f9ba6d6e4f760232 100755
--- a/source/read_obs.f90
+++ b/source/read_obs.f90
@@ -258,9 +258,9 @@ subroutine read_obs(config, files)
hl = hl - 24
endif
if ( alt.le.1000. ) then
- if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
+! if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
else
- if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
+! if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
endif
if ( conc.le.-999. ) cycle read_loop
measerr = max(config%measerr, err)
diff --git a/source/simulate.f90 b/source/simulate.f90
index 6f7759cf538cc23d252a35db8fbe7228b5648aec..999c53244f8aca7b04b6e5318d482c88c5a837f5 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -173,12 +173,14 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
ftime(n) = juldate(int(dates(n)/1d6),int(dates(n)-floor(dates(n)/1d6)*1d6))
end do
factor(:,:,:) = 1.
- filefactor = trim(path_flexrec)//'factor_drygrid'
- inquire ( file=trim(filefactor),exist=lexist )
- if ( lexist ) then
- write(logid,*) 'Reading correction factor :'//trim(filefactor)
- call read_factor(filefactor, nread, nxgrid, nygrid, nxshift, factor)
- endif
+ if ( config%spec.ne.'bca' ) then
+ filefactor = trim(path_flexrec)//'factor_drygrid'
+ inquire ( file=trim(filefactor),exist=lexist )
+ if ( lexist ) then
+ write(logid,*) 'Reading correction factor :'//trim(filefactor)
+ call read_factor(filefactor, nread, nxgrid, nygrid, nxshift, factor)
+ endif
+ endif
endif
! read flexpart footprints
@@ -260,11 +262,13 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( allocated(factor_nest) ) deallocate(factor_nest)
allocate ( factor_nest(nnxgrid,nnygrid,nread) )
factor(:,:,:) = 1.
- filefactor = trim(path_flexrec)//'factor_drygrid_nest'
- inquire ( file=trim(filefactor),exist=lexist )
- if ( lexist ) then
- write(logid,*) 'Reading correction factor for nest :'//trim(filefactor)
- call read_factor(filefactor, nread, nnxgrid, nnygrid, nnxshift, factor_nest)
+ if ( config%spec.ne.'bca' ) then
+ filefactor = trim(path_flexrec)//'factor_drygrid_nest'
+ inquire ( file=trim(filefactor),exist=lexist )
+ if ( lexist ) then
+ write(logid,*) 'Reading correction factor for nest :'//trim(filefactor)
+ call read_factor(filefactor, nread, nnxgrid, nnygrid, nnxshift, factor_nest)
+ endif
endif
endif
if ( lsatellite ) then