diff --git a/FLEXINVERT_userguide.docx b/FLEXINVERT_userguide.docx
deleted file mode 100644
index 991d85d2f0e2cdf90908a67a845822098a36e875..0000000000000000000000000000000000000000
Binary files a/FLEXINVERT_userguide.docx and /dev/null differ
diff --git a/grid_to_ncdf/SETTINGS b/grid_to_ncdf/SETTINGS
index c4f81b08ad1776246500a72eb8455271f2e9b606..0020b7786a6d42fa85c5a893138ed2c5f337d460 100644
--- a/grid_to_ncdf/SETTINGS
+++ b/grid_to_ncdf/SETTINGS
@@ -19,8 +19,8 @@ file_recept: /myfile.txt
# Run dates: format yyyymmdd
# datei = start date
# datef = end date
-datei: 20120101
-datef: 20120131
+datei: 20140701
+datef: 20140831
# Nested output (true/false)
lnested: .true.
diff --git a/prep_flexpart/process_obs.f90 b/prep_flexpart/process_obs.f90
index 3f2812b5a00908c0b8960e3c5cceb6fae78c72dd..81fe389e87280d098f08da64b040381b5e4dd4dd 100644
--- a/prep_flexpart/process_obs.f90
+++ b/prep_flexpart/process_obs.f90
@@ -94,7 +94,7 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
lonobs(1:nobs) = lonobs_out
altobs(1:nobs) = altobs_out
errobs(1:nobs) = errobs_out
- print*, 'process_obs: ind = ',ind
+! print*, 'process_obs: ind = ',ind
! if average and/or select obs
if( settings%laverage.eq.1.or.settings%lselect.eq.1 ) then
@@ -106,8 +106,8 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
eomday = calceomday(jjjjmmdd/100)
jdsom = juldate((jjjjmmdd/100)*100+1, 0)
jdeom = jdsom + real(eomday,kind=8)
- print*, 'jdsom: ',jdsom
- print*, 'jdeom: ',jdeom
+! print*, 'jdsom: ',jdsom
+! print*, 'jdeom: ',jdeom
! loop over all observations in current month
do while( jdate.le.jdobs(nobs) )
jdt = dnint(jdate*1.e6)/1.e6
@@ -141,23 +141,23 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
jdateend = min(jdateend,jdeom)
! check that not beyond end of run time
jdateend = min(jdateend,jdatef)
- print*, 'jdate: ',jdate
- print*, 'jdatestart: ',jdatestart
- print*, 'jdateend: ',jdateend
+! print*, 'jdate: ',jdate
+! print*, 'jdatestart: ',jdatestart
+! print*, 'jdateend: ',jdateend
cnt = 0
buffer(:,:) = 0.
! loop over observations in averaging interval (fill buffer)
do while( (jdt.ge.jdatestart).and.(jdt.lt.jdateend) )
call caldate(jdt, jjjjmmdd, hhmiss)
- print*, 'jdt:',jdt
- print*, 'jjjjmmdd, hhmiss:',jjjjmmdd,hhmiss
+! print*, 'jdt:',jdt
+! print*, 'jjjjmmdd, hhmiss:',jjjjmmdd,hhmiss
if( settings%lselect.eq.1 ) then
! select obs
hloc = hhmiss/10000 + int(lonobs(n)*24./360.)
- print*, 'hh:',hhmiss/10000
+! print*, 'hh:',hhmiss/10000
if( hloc.gt.24) hloc = hloc - 24
if( hloc.lt.0) hloc = hloc + 24
- print*, 'hloc:',hloc
+! print*, 'hloc:',hloc
if( altobs(n).lt.1000. ) then
! select daytime (12-16h)
if( hloc.ge.12.and.hloc.le.16) then
@@ -194,17 +194,17 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
nprev = n
n = n + 1
if ( n.gt.nobs ) go to 10
- print*, 'jdt: ',jdt
+! print*, 'jdt: ',jdt
! jdt = jdobs(n)
jdt = dnint(jdobs(n)*1.e6)/1.e6
- print*, 'next jdt: ',jdt
+! print*, 'next jdt: ',jdt
end do
if( nprev.eq.0 ) nprev = 1
10 continue
- print*, 'cnt:',cnt
- print*, 'buffer:'
+! print*, 'cnt:',cnt
+! print*, 'buffer:'
do c = 1, cnt
- print*, buffer(c,1), buffer(c,2), buffer(c,3), buffer(c,4), buffer(c,5), buffer(c,6)
+! print*, buffer(c,1), buffer(c,2), buffer(c,3), buffer(c,4), buffer(c,5), buffer(c,6)
end do
if( settings%laverage.eq.1 ) then
! average obs in buffer
diff --git a/prep_flexpart/read_basic.f90 b/prep_flexpart/read_basic.f90
index fa4629062d3d6e222315eb846f3361a6ce82ad3e..efbf336a19ccc08603981bd43bd29b2f4de5b35e 100644
--- a/prep_flexpart/read_basic.f90
+++ b/prep_flexpart/read_basic.f90
@@ -89,7 +89,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
! find file matching receptor name
do n = 1, nfiles
do i = 1, len_trim(filelist(n))-2
- if ( filelist(n)(i:(i+2)).eq.to_lower(trim(recname(nr))) ) go to 10
+ if ( to_lower(filelist(n)(i:(i+2))).eq.to_lower(trim(recname(nr))) ) go to 10
end do
end do
10 continue
@@ -151,8 +151,8 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
jdobs(cnt) = jdate
concobs(cnt) = conc
errobs(cnt) = err
- lonobs(cnt) = temp(9)
- latobs(cnt) = temp(10)
+ latobs(cnt) = temp(9)
+ lonobs(cnt) = temp(10)
altobs(cnt) = temp(11)
if ( errobs(cnt).le.-9.99 ) errobs(cnt) = 0.
end do read_loop
diff --git a/prep_flexpart/run_flexpart.sh b/prep_flexpart/run_flexpart.sh
index a64469aaf2deff00c76435b7af8374a15cfce283..ad6f3906ec09930d0659afc97526b779c1901b46 100755
--- a/prep_flexpart/run_flexpart.sh
+++ b/prep_flexpart/run_flexpart.sh
@@ -14,11 +14,9 @@
# To run tests reference to current directory
DIR=$PWD
-LEN=`expr length $DIR - 11`
+LEN=`expr length $DIR - 13`
ROOTDIR=${DIR:0:$LEN}
-echo ${ROOTDIR}
-
#---------------------------------------------------
# User settings
@@ -28,9 +26,9 @@ TIMELIM=12:00:00
PARTITION="nilu"
EXENAME=FP_ecmwf_gfortran
METEO=ecmwf
-DIRFLEX=${ROOTDIR}TEST_CASE_CH4/FLEXPART/
-DIROPTIONS=${ROOTDIR}TEST_CASE_CH4/FP_FILES/
-STNLIST=(PAL)
+DIRFLEX=/home/rthompson/REPOS/GITHUB/FLEXPART/
+DIROPTIONS=${ROOTDIR}TEST_INPUT/FLEXPART/GHG/NEST/
+STNLIST=(SSL)
MONLIST=(01)
YEAR=2012
diff --git a/prep_flexpart/sbatch_prep_flexpart.sh b/prep_flexpart/sbatch_prep_flexpart.sh
index 4373726b13c67a61611388bbcfeb8de9c3b4390c..403d6405511a1e92dd803df82d782bdb59ba1c34 100755
--- a/prep_flexpart/sbatch_prep_flexpart.sh
+++ b/prep_flexpart/sbatch_prep_flexpart.sh
@@ -4,8 +4,6 @@ partition=debug
settings_files='./SETTINGS_ghg'
#---------------------------------------------------
-cp ../../prep_flexpart/prep_flexpart .
-
cat <<EOF > run_job.sh
#!/bin/bash
./prep_flexpart ${settings_files}
diff --git a/prep_fluxes/FLUXES_ghg b/prep_fluxes/FLUXES_ghg
index 879ec1205f8bcfac43dd4a6086a195e7bab0ae24..f85dbb92658ddb787124d453c3a402988d25a60c 100644
--- a/prep_fluxes/FLUXES_ghg
+++ b/prep_fluxes/FLUXES_ghg
@@ -15,7 +15,7 @@
# ntime = number of time steps
# timeref = ref date if timestamp is julian days (yyyymmdd) otherwise 0
# timestamp = sec or month (if julian days and timeref is not 0 this is ignored)
-file_in: ../TEST_CASE_CH4/INPUT/FLUX/CH4_TOTAL_2012_05x05.nc
+file_in: /mypath/TEST_INPUT/FLUXES/GHG/CH4_TOTAL_2012_05x05.nc
varname_in: emisch4
lonname_in: longitude
latname_in: latitude
@@ -38,17 +38,17 @@ timestamp:
# ntime_out = number of time steps
# llx_out = left longitude of domain
# lly_out = lower latitude of domain
-file_out: ../TEST_CASE_CH4/INPUT/FLUX/CH4_TOTAL_2012_20x20.nc
+file_out: /mypath/TEST_OUTPUT/FLUXES/GHG/CH4_TOTAL_2012_10x10.nc
varname_out: emisch4
varunit: kgC/m2/h
lonname_out: longitude
latname_out: latitude
timename_out: time
year: 2012
-nx_out: 180
-ny_out: 90
-dx_out: 2.0
-dy_out: 2.0
+nx_out: 360
+ny_out: 180
+dx_out: 1.0
+dy_out: 1.0
ntime_out: 12
llx_out: -180
lly_out: -90
diff --git a/prep_fluxes/adjust_time.f90 b/prep_fluxes/adjust_time.f90
index 72a675009e858faebd729db8fed9b74208862d09..8b2f670c03beef544592879b9969aaea057a54c2 100644
--- a/prep_fluxes/adjust_time.f90
+++ b/prep_fluxes/adjust_time.f90
@@ -53,7 +53,7 @@ subroutine adjust_time(settings, time, time_out)
! input timestamp is julian days
do i = 1, ntime_out
time_out(i) = time(i) + juldate(int(settings%timeref),0) - &
- juldate(19000101,0) + 1.
+ juldate(19000101,0)
end do
! adjust time stamp to given year
call caldate(time_out(1) + juldate(19000101,0), yyyymmdd, hhmiss)
diff --git a/prep_fluxes/convertgrid.f90 b/prep_fluxes/convertgrid.f90
index 05d0dd96e3593ec05e97c63324279e3498d945dc..f46c18d750b4b22b7629527b5d336a3f59f63d81 100644
--- a/prep_fluxes/convertgrid.f90
+++ b/prep_fluxes/convertgrid.f90
@@ -111,9 +111,12 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
end do
dy_in(1) = dy_in(2)
- ! round dx and dy to 7dp
- dx_in = anint(dx_in*1.e7)/1.e7
- dy_in = anint(dy_in*1.e7)/1.e7
+ ! round dx and dy to 4dp
+ dx_in = anint(dx_in*1.e4)/1.e4
+ dy_in = anint(dy_in*1.e4)/1.e4
+ ! round lon_in and lat_in to 4dp
+ lon_in = anint(lon_in*1.e4)/1.e4
+ lat_in = anint(lat_in*1.e4)/1.e4
! if not midpoints check if west and south boundaries given
if(.not.midpoints) then
@@ -159,6 +162,10 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
end do
endif
+ ! round wlon_in and slat_in to 4dp
+ wlon_in = anint(wlon_in*1.e4)/1.e4
+ slat_in = anint(slat_in*1.e4)/1.e4
+
if( (nx_in.eq.nx_out) .and. (ny_in.eq.ny_out) ) then
write(*,*) 'Comment: input and output grids are the same'
wlon_out = wlon_in
@@ -268,7 +275,7 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
! calculate total in
do nt = 1, ntime_out
- tot_in = tot_in + sum( sum(flx_tmp(:,:,nt)*area_xy(:,:),dim=1) )
+ tot_in = tot_in + sum( flx_tmp(:,:,nt)*area_xy(:,:) )
end do
print*, tot_in*24.*365./float(ntime_out)/1.e9
diff --git a/prep_fluxes/rev_test.out b/prep_fluxes/rev_test.out
deleted file mode 100644
index de4df400863aca9f0ce1bb5c23033c99bb3ea0ad..0000000000000000000000000000000000000000
--- a/prep_fluxes/rev_test.out
+++ /dev/null
@@ -1,41 +0,0 @@
- /home/rthompson/REPOS/GITHUB/FLEXINVERTplus/TEST_OUTPUT/FLUXES/CH4_TOTAL_2012_10x10.nc
- emisch4
- longitude
- latitude
- time
- T
- 360
- 180
- 12
- 19000101
- /home/rthompson/REPOS/GITHUB/FLEXINVERTplus/TEST_OUTPUT/FLUXES/CH4_TOTAL_2012_05x05_test.nc
- emisch4
- longitude
- latitude
- time
- 2012
- 720
- 360
- 0.500000000
- 0.500000000
- 12
- -180.000000
- -90.0000000
- nx_in, ny_in: 360 180
- ntime: 12
- nx_out, ny_out: 720 360
- dx_out, dy_out: 0.500000000 0.500000000
- llx_out, lly_out: -180.000000 -90.0000000
- indextime = 1
- Reading file: /home/rthompson/REPOS/GITHUB/FLEXINVERTplus/TEST_OUTPUT/FLUXES/CH4_TOTAL_2012_10x10.nc
- time = 40907.0000 40938.0000 40966.0000 40997.0000 41027.0000 41058.0000 41088.0000 41119.0000 41150.0000 41180.0000 41211.0000 41241.0000
- lat_in: -89.5000000 -88.5000000 -87.5000000 -86.5000000 -85.5000000 -84.5000000 -83.5000000 -82.5000000 -81.5000000 -80.5000000 -79.5000000 -78.5000000 -77.5000000 -76.5000000 -75.5000000 -74.5000000 -73.5000000 -72.5000000 -71.5000000 -70.5000000 -69.5000000 -68.5000000 -67.5000000 -66.5000000 -65.5000000 -64.5000000 -63.5000000 -62.5000000 -61.5000000 -60.5000000 -59.5000000 -58.5000000 -57.5000000 -56.5000000 -55.5000000 -54.5000000 -53.5000000 -52.5000000 -51.5000000 -50.5000000 -49.5000000 -48.5000000 -47.5000000 -46.5000000 -45.5000000 -44.5000000 -43.5000000 -42.5000000 -41.5000000 -40.5000000 -39.5000000 -38.5000000 -37.5000000 -36.5000000 -35.5000000 -34.5000000 -33.5000000 -32.5000000 -31.5000000 -30.5000000 -29.5000000 -28.5000000 -27.5000000 -26.5000000 -25.5000000 -24.5000000 -23.5000000 -22.5000000 -21.5000000 -20.5000000 -19.5000000 -18.5000000 -17.5000000 -16.5000000 -15.5000000 -14.5000000 -13.5000000 -12.5000000 -11.5000000 -10.5000000 -9.50000000 -8.50000000 -7.50000000 -6.50000000 -5.50000000 -4.50000000 -3.50000000 -2.50000000 -1.50000000 -0.500000000 0.500000000 1.50000000 2.50000000 3.50000000 4.50000000 5.50000000 6.50000000 7.50000000 8.50000000 9.50000000 10.5000000 11.5000000 12.5000000 13.5000000 14.5000000 15.5000000 16.5000000 17.5000000 18.5000000 19.5000000 20.5000000 21.5000000 22.5000000 23.5000000 24.5000000 25.5000000 26.5000000 27.5000000 28.5000000 29.5000000 30.5000000 31.5000000 32.5000000 33.5000000 34.5000000 35.5000000 36.5000000 37.5000000 38.5000000 39.5000000 40.5000000 41.5000000 42.5000000 43.5000000 44.5000000 45.5000000 46.5000000 47.5000000 48.5000000 49.5000000 50.5000000 51.5000000 52.5000000 53.5000000 54.5000000 55.5000000 56.5000000 57.5000000 58.5000000 59.5000000 60.5000000 61.5000000 62.5000000 63.5000000 64.5000000 65.5000000 66.5000000 67.5000000 68.5000000 69.5000000 70.5000000 71.5000000 72.5000000 73.5000000 74.5000000 75.5000000 76.5000000 77.5000000 78.5000000 79.5000000 80.5000000 81.5000000 82.5000000 83.5000000 84.5000000 85.5000000 86.5000000 87.5000000 88.5000000 89.5000000
- lon_in: -179.500000 -178.500000 -177.500000 -176.500000 -175.500000 -174.500000 -173.500000 -172.500000 -171.500000 -170.500000 -169.500000 -168.500000 -167.500000 -166.500000 -165.500000 -164.500000 -163.500000 -162.500000 -161.500000 -160.500000 -159.500000 -158.500000 -157.500000 -156.500000 -155.500000 -154.500000 -153.500000 -152.500000 -151.500000 -150.500000 -149.500000 -148.500000 -147.500000 -146.500000 -145.500000 -144.500000 -143.500000 -142.500000 -141.500000 -140.500000 -139.500000 -138.500000 -137.500000 -136.500000 -135.500000 -134.500000 -133.500000 -132.500000 -131.500000 -130.500000 -129.500000 -128.500000 -127.500000 -126.500000 -125.500000 -124.500000 -123.500000 -122.500000 -121.500000 -120.500000 -119.500000 -118.500000 -117.500000 -116.500000 -115.500000 -114.500000 -113.500000 -112.500000 -111.500000 -110.500000 -109.500000 -108.500000 -107.500000 -106.500000 -105.500000 -104.500000 -103.500000 -102.500000 -101.500000 -100.500000 -99.5000000 -98.5000000 -97.5000000 -96.5000000 -95.5000000 -94.5000000 -93.5000000 -92.5000000 -91.5000000 -90.5000000 -89.5000000 -88.5000000 -87.5000000 -86.5000000 -85.5000000 -84.5000000 -83.5000000 -82.5000000 -81.5000000 -80.5000000 -79.5000000 -78.5000000 -77.5000000 -76.5000000 -75.5000000 -74.5000000 -73.5000000 -72.5000000 -71.5000000 -70.5000000 -69.5000000 -68.5000000 -67.5000000 -66.5000000 -65.5000000 -64.5000000 -63.5000000 -62.5000000 -61.5000000 -60.5000000 -59.5000000 -58.5000000 -57.5000000 -56.5000000 -55.5000000 -54.5000000 -53.5000000 -52.5000000 -51.5000000 -50.5000000 -49.5000000 -48.5000000 -47.5000000 -46.5000000 -45.5000000 -44.5000000 -43.5000000 -42.5000000 -41.5000000 -40.5000000 -39.5000000 -38.5000000 -37.5000000 -36.5000000 -35.5000000 -34.5000000 -33.5000000 -32.5000000 -31.5000000 -30.5000000 -29.5000000 -28.5000000 -27.5000000 -26.5000000 -25.5000000 -24.5000000 -23.5000000 -22.5000000 -21.5000000 -20.5000000 -19.5000000 -18.5000000 -17.5000000 -16.5000000 -15.5000000 -14.5000000 -13.5000000 -12.5000000 -11.5000000 -10.5000000 -9.50000000 -8.50000000 -7.50000000 -6.50000000 -5.50000000 -4.50000000 -3.50000000 -2.50000000 -1.50000000 -0.500000000 0.500000000 1.50000000 2.50000000 3.50000000 4.50000000 5.50000000 6.50000000 7.50000000 8.50000000 9.50000000 10.5000000 11.5000000 12.5000000 13.5000000 14.5000000 15.5000000 16.5000000 17.5000000 18.5000000 19.5000000 20.5000000 21.5000000 22.5000000 23.5000000 24.5000000 25.5000000 26.5000000 27.5000000 28.5000000 29.5000000 30.5000000 31.5000000 32.5000000 33.5000000 34.5000000 35.5000000 36.5000000 37.5000000 38.5000000 39.5000000 40.5000000 41.5000000 42.5000000 43.5000000 44.5000000 45.5000000 46.5000000 47.5000000 48.5000000 49.5000000 50.5000000 51.5000000 52.5000000 53.5000000 54.5000000 55.5000000 56.5000000 57.5000000 58.5000000 59.5000000 60.5000000 61.5000000 62.5000000 63.5000000 64.5000000 65.5000000 66.5000000 67.5000000 68.5000000 69.5000000 70.5000000 71.5000000 72.5000000 73.5000000 74.5000000 75.5000000 76.5000000 77.5000000 78.5000000 79.5000000 80.5000000 81.5000000 82.5000000 83.5000000 84.5000000 85.5000000 86.5000000 87.5000000 88.5000000 89.5000000 90.5000000 91.5000000 92.5000000 93.5000000 94.5000000 95.5000000 96.5000000 97.5000000 98.5000000 99.5000000 100.500000 101.500000 102.500000 103.500000 104.500000 105.500000 106.500000 107.500000 108.500000 109.500000 110.500000 111.500000 112.500000 113.500000 114.500000 115.500000 116.500000 117.500000 118.500000 119.500000 120.500000 121.500000 122.500000 123.500000 124.500000 125.500000 126.500000 127.500000 128.500000 129.500000 130.500000 131.500000 132.500000 133.500000 134.500000 135.500000 136.500000 137.500000 138.500000 139.500000 140.500000 141.500000 142.500000 143.500000 144.500000 145.500000 146.500000 147.500000 148.500000 149.500000 150.500000 151.500000 152.500000 153.500000 154.500000 155.500000 156.500000 157.500000 158.500000 159.500000 160.500000 161.500000 162.500000 163.500000 164.500000 165.500000 166.500000 167.500000 168.500000 169.500000 170.500000 171.500000 172.500000 173.500000 174.500000 175.500000 176.500000 177.500000 178.500000 179.500000
- ix1, ix2, jy1, jy2: 1 360 1 180
- 555.843140
- 558.101868
- Total mass in: 555.8431 Tg/y
- Total mass out: 558.1019 Tg/y
- time: 40908.0000 40939.0000 40967.0000 40998.0000 41028.0000 41059.0000 41089.0000 41120.0000 41151.0000 41181.0000 41212.0000 41242.0000
- ntime_out: 12
- Writing file: /home/rthompson/REPOS/GITHUB/FLEXINVERTplus/TEST_OUTPUT/FLUXES/CH4_TOTAL_2012_05x05_test.nc
diff --git a/prep_regions/SETTINGS_regions b/prep_regions/SETTINGS_regions
index b24f8b4fba412f8f0e7251146ec09149a059982e..d3b696f02419c3a426b99107e8ff2956bd5a8fcd 100644
--- a/prep_regions/SETTINGS_regions
+++ b/prep_regions/SETTINGS_regions
@@ -7,9 +7,6 @@
#
# ==================================================
-# Define map globally with "zoom" (true) or regionally (false)
-usezoom: .false.
-
# Use prior fluxes
# weight surface influnce by annual fluxes (true/false)
useflux: .true.
@@ -33,3 +30,10 @@ nmap: 3
# 3) 1 x 0.5 = 0.5 degree
mapdec: 1,2,2
+# Optional: snow-ice mask file
+# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
+file_ice:
+varname_ice:
+lonname_ice:
+latname_ice:
+
diff --git a/prep_regions/SETTINGS_zoom_regions b/prep_regions/SETTINGS_zoom_regions
index 772d3880a81e9c32cbbbb4fb0a9b2cd7126841af..14a38992dc550d5ad7dccea7d5556c7950227ea0 100644
--- a/prep_regions/SETTINGS_zoom_regions
+++ b/prep_regions/SETTINGS_zoom_regions
@@ -43,3 +43,10 @@ zlly: 38
zurx: 28
zury: 70
+# Optional: snow-ice mask file
+# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
+file_ice:
+varname_ice:
+lonname_ice:
+latname_ice:
+
diff --git a/prep_regions/map_grid.f90 b/prep_regions/map_grid.f90
index 75e371a8102561da64089a25f336c674bb4a3f87..a50b887248ee133c08b1e9cb61bdefcad3db9d2b 100644
--- a/prep_regions/map_grid.f90
+++ b/prep_regions/map_grid.f90
@@ -78,7 +78,7 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
integer :: numx, numy, num
integer :: nres, ndec, nfx, nfy, nfbox, nsort
integer, dimension(2*maxbox) :: map_box
- integer :: nbox_invert, nb, nbox_land, nbox_water
+ integer :: nbox_invert, nb, nbprev, nbox_land, nbox_water
real, dimension(2*maxbox) :: x_map, y_map, sens_map, water_box
real, dimension(nxregrid,nyregrid) :: sens_lump, contrib_lump, water
real, dimension(nxregrid*nyregrid) :: sort_contrib, sort_contrib_inv
@@ -258,6 +258,29 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
print*, 'nbox_invert, nb = ',nbox_invert, nb
+ ! if a box contains land and water split it
+! nbprev = 0
+! do nb = 1, nbox_invert
+! if ( water_box(nb).le.0.99 ) then
+! ! land - check if any water in this box
+! do jy = 1, nyregrid
+! do ix = 1, nxregrid
+! if ( (nbox_xy(ix,jy).eq.nb).and.lsm(ix,jy).lt.0.01 ) then
+! ! make a new box for water
+! if ( nb.ne.nbprev ) then
+! nbox_invert = nbox_invert + 1
+! water_box(nbox_invert) = 1.
+! endif
+! nbox_xy(ix,jy) = nbox_invert
+! nbprev = nb
+! endif
+! end do
+! end do
+! endif
+! end do
+
+ print*, 'after splitting land/water: nbox_invert, nb = ',nbox_invert, nb
+
! renumber boxes for land (1:nbox_land) and water (-1:-nbox_water)
nbox_water = 0
nbox_land = 0
diff --git a/prep_regions/mod_settings.f90 b/prep_regions/mod_settings.f90
index 44563c5c49cafd9e7ed4f861c114c8e875295933..f928a2996f09f0f0e5e923d9a76ce7e58d038132 100644
--- a/prep_regions/mod_settings.f90
+++ b/prep_regions/mod_settings.f90
@@ -90,7 +90,7 @@ module mod_settings
real :: n_edge_lat ! lat of northern edge of inversion grid
real :: xres ! longitudinal grid resolution
real :: yres ! latitudinal grid resolution
- integer :: nt_flx ! number of flux time steps
+ integer :: nstep_flx ! number of flux time steps
end type config_t
@@ -289,9 +289,6 @@ module mod_settings
logical :: match, cl
- ! default values
- files%file_ice = ""
-
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
if(ierr.gt.0) then
@@ -363,20 +360,6 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%latname_lsm = cc
- ! read snow-ice mask settings
- identifier = "file_ice:"
- call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) files%file_ice = cc
- identifier = "varname_ice:"
- call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) files%varname_ice = cc
- identifier = "lonname_ice:"
- call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) files%lonname_ice = cc
- identifier = "latname_ice:"
- call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) files%latname_ice = cc
-
! read flux prior settings
identifier = "filename_flx:"
call read_content (line, identifier, cc, cn, cl, match)
@@ -523,9 +506,9 @@ module mod_settings
identifier = "yres:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%yres = real(cn)
- identifier = "nt_flx:"
+ identifier = "nstep_flx:"
call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) config%nt_flx = int(cn)
+ if ( match ) config%nstep_flx = int(cn)
endif
@@ -537,9 +520,10 @@ module mod_settings
! read region settings
! --------------------------------------------------
- subroutine read_region_settings(filename)
+ subroutine read_region_settings(filename, files)
- character(len=200), intent(in) :: filename
+ type(files_t), intent(in out) :: files
+ character(len=200), intent(in) :: filename
integer :: ierr, i, n
character(len=200) :: line, identifier, cc
@@ -554,6 +538,9 @@ module mod_settings
stop
endif
+ ! default values
+ files%file_ice = ""
+
! loop over lines
read_loop: do
read (100, fmt='(A)', iostat=ierr) line
@@ -610,6 +597,20 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) zury = int(cn)
+ ! read snow-ice mask settings
+ identifier = "file_ice:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%file_ice = cc
+ identifier = "varname_ice:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%varname_ice = cc
+ identifier = "lonname_ice:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%lonname_ice = cc
+ identifier = "latname_ice:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%latname_ice = cc
+
endif
end do read_loop
diff --git a/prep_regions/mod_var.f90 b/prep_regions/mod_var.f90
index 5b63e81b18cd4b7882b41786240c2ba924ce4fdf..14a5a58f85fb7eb6fa9cf18f3be5fe8a7df2318c 100644
--- a/prep_regions/mod_var.f90
+++ b/prep_regions/mod_var.f90
@@ -35,7 +35,7 @@ module mod_var
integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
- integer, parameter :: maxobs=10000 ! max number of observations
+ integer, parameter :: maxobs=100000 ! max number of observations
! flexpart variables
! ------------------
@@ -101,7 +101,10 @@ module mod_var
integer :: nrec ! number of receptor sites or campaigns
integer :: nobs ! total number of observations
character(recname_len), dimension(:), allocatable :: recname ! names of receptors
-
+ real, dimension(:), allocatable :: reclat ! latitude of receptor
+ real, dimension(:), allocatable :: reclon ! longitude of receptor
+ real, dimension(:), allocatable :: recalt ! altitude of receptor
+
end module mod_var
diff --git a/prep_regions/prep_regions.f90 b/prep_regions/prep_regions.f90
index 29f5efe3795e07402500353de3a5e045b12a4e8a..0c25d3709d1c9ae15ccc509d90fe432875d4bece 100644
--- a/prep_regions/prep_regions.f90
+++ b/prep_regions/prep_regions.f90
@@ -78,7 +78,7 @@ program prep_regions
call read_config_settings(settings_config, config)
! read region settings
- call read_region_settings(settings_regions)
+ call read_region_settings(settings_regions, files)
! initialize parameters
call initialize(files, config)
diff --git a/prep_regions/read_lsm.f90 b/prep_regions/read_lsm.f90
index 8b74d050718b29877e6b02e16b5d5a546efe0a8a..c5577ba31d53029d1b308c0ea7bc32fe6e1b0045 100644
--- a/prep_regions/read_lsm.f90
+++ b/prep_regions/read_lsm.f90
@@ -175,7 +175,7 @@ subroutine read_lsm(files)
! integrate ice into lsm
ice_av = sum(ice,dim=3)/12.
- where( ice_av(:,:).gt.0.5 ) lsm = 0.
+ where( ice_av(:,:).gt.0.99 ) lsm = 0.
deallocate(lon)
deallocate(lat)
diff --git a/prep_regions/read_obs.f90 b/prep_regions/read_obs.f90
index c083efdcf336abb0ddf454c18cf3dc9afbb446bd..33a87bbcd29754d49b6178669d084451291cca6b 100644
--- a/prep_regions/read_obs.f90
+++ b/prep_regions/read_obs.f90
@@ -53,10 +53,10 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
integer :: ierr
integer :: cnt
integer :: narg, n, reclen
- integer :: nf, nfiles
- integer :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss
+ integer :: nf, nfiles, nr
+ integer :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss, hl
real(kind=8) :: jdate, avetime
- real :: conc, err
+ real :: conc, err, lon, alt
integer :: num
! list observation files
@@ -108,6 +108,16 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
endif
write(*,*) 'Reading file: '//trim(files%path_obs)//trim(filelist(nf))
+ ! station coordinates for data selection
+ do nr = 1, nrec
+ if ( to_upper(recs).eq.recname(nr) ) then
+ lon = reclon(nr)
+ alt = recalt(nr)
+ go to 20
+ endif
+ end do
+20 continue
+
! read header
read (100, fmt='(A)', iostat=ierr) header
print*, 'read_obs: header = ',header
@@ -122,19 +132,34 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
endif
- ! if obs file doesn't contain avetime (backwards compatability)
+ ! if obs file doesn't contain avetime (backwards compatability)
! avetime = 0d0
! if ( reclen.eq.4 ) then
! read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
! recs, yyyymmdd, hhmiss, jdate, conc, err, num
! else
-! read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
+! read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
! recs, yyyymmdd, hhmiss, jdate, conc, err, num
! endif
if ( ierr.gt.0 ) exit read_loop
if ( jdate.ge.(juldatef+1d0) ) exit read_loop
if ( jdate.lt.juldatei ) cycle read_loop
if ( conc.le.-999. ) cycle read_loop
+ ! select day/night for low/high alt sites
+ hh = hhmiss/10000
+ hl = hh + int(lon*24./360.)
+ if ( hl.lt.0 ) then
+ hl = hl + 24
+ else if ( hl.ge.24 ) then
+ hl = hl - 24
+ endif
+ if ( alt.le.1000. ) then
+ if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
+ else
+ if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
+ endif
+ if ( conc.le.-999. ) cycle read_loop
+
cnt = cnt + 1
obstimes(cnt) = jdate
avetimes(cnt) = avetime
diff --git a/prep_regions/read_reclist.f90 b/prep_regions/read_reclist.f90
index 977900bed4c03f500823640ef679d064d8ee5b1c..82bb1f4e7fcd86a09036135dcf2bd823a4a6c77d 100644
--- a/prep_regions/read_reclist.f90
+++ b/prep_regions/read_reclist.f90
@@ -33,7 +33,6 @@ subroutine read_reclist(filename)
character(len=200) :: line
integer :: ierr
integer :: cnt
- real, dimension(:), allocatable :: reclat, reclon, recalt
! count number of receptors
diff --git a/prep_syndata/README_prepsyndata.txt b/prep_syndata/README_prepsyndata.txt
index 862fc1d44d5457d115daa6a90e09fd677803f790..a51cd40bebd4d7c8fcfc965127f2d34e36cb5e40 100644
--- a/prep_syndata/README_prepsyndata.txt
+++ b/prep_syndata/README_prepsyndata.txt
@@ -1,6 +1,6 @@
================================================================
- FLEXINVERT-Plus SYNTHETIC DATA PREPARATION
+ FLEXINVERT SYNTHETIC DATA PREPARATION
================================================================
@@ -9,7 +9,7 @@ Description:
Notes:
1) Output from a forward run of FLEXINVERT (run_mode = 0)
- are needed before running prep_syndata
+ are needed before running synthetic data test.
2) For CO2 only NEE fluxes are perturbed.
3) Fluxes are returned on the same domain as the flexpart
grid_time files (if nested is true this is the nested
@@ -17,8 +17,7 @@ Notes:
greater domain than the flexpart files.
4) If the inversion domain (covered by state vector) is
smaller than the flexpart domain, only the fluxes inside
- the inversion domain are perturbed (fluxes outside this
- domain are equal to the prior).
+ the inversion domain are perturbed.
Synthetic Data:
1) Prior fluxes: prepared by randomly perturbing prior fluxes
diff --git a/prep_syndata/calc_error_cov.f90 b/prep_syndata/calc_error_cov.f90
index d091d2792caea26d62bb703e95dbf2b6ecd53e52..61119866cb34150f6298406ed7c11204d413e78e 100644
--- a/prep_syndata/calc_error_cov.f90
+++ b/prep_syndata/calc_error_cov.f90
@@ -187,8 +187,8 @@ subroutine calc_error_cov(config, files, state)
! truncate eigenvalues
! eigenvalues are sorted in ascending order
- ind = minloc(evaltmp, dim=1, mask=evaltmp.ge.(trunc*maxval(evaltmp)))
- write(logid,*) 'Truncating eigenvalues at: ',trunc*maxval(evaltmp)
+ ind = minloc(evaltmp, dim=1, mask=evaltmp.ge.(config%trunc*maxval(evaltmp)))
+ write(logid,*) 'Truncating eigenvalues at: ',config%trunc*maxval(evaltmp)
neig = ndvar - ind + 1
write(logid,*) 'Number of retained eigenvalues: ',neig
diff --git a/prep_syndata/get_flux.f90 b/prep_syndata/get_flux.f90
index bb88d6c89682131d48bdc57e98e661ab514a7dbe..8554d8a11201e6b94e708634090f73c428ee0464 100644
--- a/prep_syndata/get_flux.f90
+++ b/prep_syndata/get_flux.f90
@@ -137,7 +137,7 @@ subroutine get_flux(config, files, flx, state)
deallocate( tmp )
! flux time stamp is start of time step
flx_time(n) = jd
- jd = jd + flxtres
+ jd = jd + dble(nt_flx)/24d0
end do
else if ( config%spec.eq.'co2' ) then
diff --git a/prep_syndata/get_obs.f90 b/prep_syndata/get_obs.f90
index 3cd38b922e89314a6a17f883377ca84e1193aa0d..e248e4c4582ae1adbdbd62e51d3de49ea072665a 100644
--- a/prep_syndata/get_obs.f90
+++ b/prep_syndata/get_obs.f90
@@ -92,7 +92,9 @@ subroutine get_obs(config, files, obs, obserr)
! check obserr
do i = 1, nobs
- if( obserr(i).eq.9.99 .or. obserr(i).eq.0. ) obserr(i) = config%measerr
+! if( obserr(i).eq.9.99 .or. obserr(i).eq.0. ) obserr(i) = config%measerr
+ ! use only measerr
+ obserr(i) = config%measerr
end do
if ( config%spec.eq.'co2' ) then
diff --git a/prep_syndata/initialize.f90 b/prep_syndata/initialize.f90
index 18ab77fc5eef53bcf165875d44ecae9b97dff4ac..81f4c267bc9935426b921715ca16bfe837603ade 100644
--- a/prep_syndata/initialize.f90
+++ b/prep_syndata/initialize.f90
@@ -39,7 +39,7 @@ subroutine initialize(files, config)
character(max_path_len) :: filename
character(max_path_len) :: rowfmt
logical :: lexist
- integer :: yyyymm
+ integer :: yyyymm, yyyymmdd, hhmiss
character(6) :: adate
character(4) :: ayear
integer :: numpoint
@@ -114,37 +114,54 @@ subroutine initialize(files, config)
! number days in inversion interval
nday = int(juldatef - juldatei + 1.)
write(logid,*) 'initialize: nday = ',nday
-
+
+ ! timestep of fluxes (hours)
+ nt_nee = config%nstep_nee
+ nt_flx = config%nstep_flx
+
! number of NEE fields per day
if ( config%spec.eq.'co2' ) then
- nd_nee = config%num_nee_day
+ if ( nt_nee.eq.0 ) then
+ write(logid,*) 'ERROR: nt_nee cannot be zero'
+ stop
+ endif
+ nd_nee = 24/nt_nee
endif
if( config%spec.eq.'co2' ) then
! statres is time interval over which state variables are averaged
! statres_hr is the sub-daily time interval for state variables
! ndt is number of time steps per 24h in prior NEE
- ! flxtres is temporal resolution of low frequency fluxes (biomass burning and ocean)
! ntstate is number of averaged time invervals for state variables
+ ! nflx is number of prior flux fields in inversion interval
statres = config%statres
statres_hr = config%statres_hr
ndt = int(24./statres_hr)
- flxtres = aint(365./real(config%nt_flx))
- ntstate = nint(real(nday)/statres)
- nt_flx = nday*nd_nee
- nflx = nt_flx
+ ntstate = int(real(nday)/statres)
+! nflx = nday*nd_nee
else
! statres is time resolution of state variables
- ! flxtres is temporal resolution of fluxes (same as state variables for non-CO2 tracers)
- ! nflx is the number of input flux timesteps falling within the inversion interval
! ntstate is also number of state time invervals
+ ! nflx is number of prior flux fields in inversion interval
statres = config%statres
- flxtres = aint(365./real(config%nt_flx))
- ntstate = nint(real(nday)/statres)
+ ntstate = int(real(nday)/statres)
ndt = 1
- nt_flx = config%nt_flx
- nflx = max(1, nday/int(365./real(config%nt_flx)))
+! nflx = nday*24/nt_flx
endif
+
+ ! adjust nday to match number of state time steps
+ nday = ntstate*int(statres)
+ write(logid,*) 'initialize: adjusted nday to match number of state time steps = ',nday
+ juldatef = juldatei + dble(nday) - 1d0
+ call caldate(juldatef, yyyymmdd, hhmiss)
+ write(logid,*) 'initialize: adjusted end date to match number of state time steps = ',yyyymmdd
+
+ if ( config%spec.eq.'co2' ) then
+ nflx = nday*nd_nee
+ else
+ nflx = nday*24/nt_flx
+ endif
+
if ( ntstate.lt.1 ) then
write(logid,*) 'ERROR initialize: number flux time steps < 1 -> increase time interval'
stop
@@ -159,6 +176,8 @@ subroutine initialize(files, config)
print*, 'ntstate: ',ntstate
print*, 'nt_flx: ',nt_flx
print*, 'nflx: ',nflx
+ print*, 'nt_nee: ',nt_nee
+ print*, 'nd_nee: ',nd_nee
print*, 'nday: ',nday
print*, 'nyr: ',nyr
diff --git a/prep_syndata/makefile b/prep_syndata/makefile
index 827e75ab9bcd6cad951ef676bce2e0b158d2b8a2..1fb0149330a17e737df4f3e6a43cb749e13b622a 100644
--- a/prep_syndata/makefile
+++ b/prep_syndata/makefile
@@ -6,7 +6,7 @@ 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
+ -fbacktrace -ffree-line-length-0 -fcheck=all -Wall
LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
MAIN = prep_syndata
diff --git a/prep_syndata/mod_save.f90 b/prep_syndata/mod_save.f90
index d63a2f20da88537d030c5330e4831802faead529..a06c67832721028e568c2054131f6cbd8c36c511 100644
--- a/prep_syndata/mod_save.f90
+++ b/prep_syndata/mod_save.f90
@@ -60,7 +60,7 @@ module mod_save
type (config_t), intent (in) :: config
type (files_t), intent (in) :: files
- real, dimension(nxgrid,nygrid,nflx), intent (in out) :: flx
+ real, dimension(nxgrid,nygrid,nflx), intent (in out) :: flx
real, dimension(nvar), intent (in out) :: state
character(len=max_name_len) :: filename, lonname, latname, timename, varname
@@ -164,9 +164,6 @@ module mod_save
! interpolate between previous and current timesteps
if ( config%lognormal ) then
! lognormal
-! flx_out(ix,jy,n) = state((nt-2)*nbox+nbox_xy(ix,jy)) + &
-! (state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
-! (flx_time(n) - stat_time(nt))/statres
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * ( state((nt-2)*nbox+nbox_xy(ix,jy)) + &
(state((nt-1)*nbox+nbox_xy(ix,jy))-state((nt-2)*nbox+nbox_xy(ix,jy))) * &
(flx_time(n) - stat_time(nt))/statres )
@@ -181,9 +178,6 @@ module mod_save
! interpolate between current and next timesteps
if ( config%lognormal ) then
! lognormal
-! flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + &
-! (state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
-! (stat_time(nt) - flx_time(n))/statres
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * ( state((nt-1)*nbox+nbox_xy(ix,jy)) + &
(state(nt*nbox+nbox_xy(ix,jy))-state((nt-1)*nbox+nbox_xy(ix,jy))) * &
(stat_time(nt) - flx_time(n))/statres )
@@ -197,7 +191,6 @@ module mod_save
else
! take current timestep
if ( config%lognormal ) then
-! flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy))
flx_out(ix,jy,n) = flx(ix1+ix-1,jy1+jy-1,n) * state((nt-1)*nbox+nbox_xy(ix,jy))
else
flx_out(ix,jy,n) = state((nt-1)*nbox+nbox_xy(ix,jy)) + flx(ix1+ix-1,jy1+jy-1,n)
@@ -215,9 +208,6 @@ module mod_save
! replace fluxes in inversion domain with perturbed fluxes
flx(ix1:ix2,jy1:jy2,:) = flx_out(:,:,:)
-!*** FOR SF6
- where ( flx(:,:,:).lt.0. ) flx(:,:,:) = 0.
-
! undo numerical scaling
flx = flx/numscale
! convert to kg/m2/h
@@ -226,7 +216,7 @@ module mod_save
print*, 'mod_save: min/max flx = ',minval(flx), maxval(flx)
! adjust time to days since 1-Jan-1900
- flx_time = flx_time - juldate(19000101, 0) + 1d0
+ flx_time = flx_time - juldate(19000101, 0) !+ 1d0
! write flux to file
! ------------------
@@ -279,7 +269,7 @@ module mod_save
call check( nf90_def_dim(ncid, latname, nygrid, lat_dimid) )
call check( nf90_def_var(ncid, latname, nf90_real, lat_dimid, lat_varid) )
call check( nf90_put_att(ncid, lat_varid, units, dimunit) )
- call check( nf90_def_dim(ncid, timename, nt_flx, time_dimid) )
+ call check( nf90_def_dim(ncid, timename, nflx, time_dimid) )
call check( nf90_def_var(ncid, timename, nf90_real, time_dimid, time_varid) )
call check( nf90_put_att(ncid, time_varid, units, timeunit) )
dimids = (/ lon_dimid, lat_dimid, time_dimid /)
@@ -322,7 +312,7 @@ module mod_save
character(len=recname_len), dimension(nobs) :: recs
integer, dimension(nobs) :: jjjjmmdd, hhmiss
real(kind=8), dimension(nobs) :: jdate
- real :: conc, cini, bkg, nee, fff, ocn, bbg, ghg, prior, post, delta, err
+ real :: conc, cini, cinipos, bkg, nee, fff, ocn, ghg, prior, post, delta, err
integer :: i, ierr
logical :: lexist
@@ -341,19 +331,14 @@ module mod_save
if ( config%spec.eq.'co2' ) then
write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',11,'(F10.4,1X))'
do i = 1, nobs
- read(100,fmt=rowfmt) recs(i), jjjjmmdd(i), hhmiss(i), jdate(i), conc, cini, bkg, &
- nee, fff, ocn, bbg, prior, post, delta, err
- end do
+ read(100,fmt=rowfmt) recs(i), jjjjmmdd(i), hhmiss(i), jdate(i), conc, cini, cinipos, bkg, &
+ nee, fff, ocn, prior, post, delta, err
+ end do
else
- write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',8,'(F10.4,1X))'
+ write(rowfmt,'(A,I6,A)') '(A3,1X,I8,1X,I6,1X,F14.6,1X,',9,'(F10.4,1X))'
do i = 1, nobs
- if ( config%offsets ) then
- read(100,fmt=rowfmt) recs(i), jjjjmmdd(i), hhmiss(i), jdate(i), conc, cini, bkg, &
- ghg, prior, post, delta, err
- else
- read(100,fmt=rowfmt) recs(i), jjjjmmdd(i), hhmiss(i), jdate(i), conc, cini, bkg, &
- ocn, prior, post, delta, err
- endif
+ read(100,fmt=rowfmt) recs(i), jjjjmmdd(i), hhmiss(i), jdate(i), conc, cini, cinipos, bkg, &
+ ghg, prior, post, delta, err
end do
endif
close(100)
diff --git a/prep_syndata/mod_settings.f90 b/prep_syndata/mod_settings.f90
index 479269b226872038988bbf5659450ddbc4570a4d..15d09a1c47c455ce726b001c8137a370982825bc 100644
--- a/prep_syndata/mod_settings.f90
+++ b/prep_syndata/mod_settings.f90
@@ -91,12 +91,13 @@ module mod_settings
real :: n_edge_lat ! lat of northern edge of inversion grid
real :: xres ! longitudinal grid resolution
real :: yres ! latitudinal grid resolution
- integer :: nt_flx ! number of flux time steps
- integer :: num_nee_day ! number of NEE time steps per day
+ integer :: nstep_flx ! number of flux time steps
+ integer :: nstep_nee ! number of NEE time steps per day
real :: statres ! number of flux time steps
real :: statres_hr ! sub-daily temporal resolution of state vector (e.g. 6, 12 or 24h)
real :: measerr ! measurement error
logical :: lognormal ! use lognormal distribution in state space (true or false)
+ real :: trunc ! truncation factor for eigenvalues of B
end type config_t
@@ -461,6 +462,7 @@ module mod_settings
! logical defaults
config%lognormal = .false.
+ config%trunc = 1.e-3
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
@@ -512,6 +514,9 @@ module mod_settings
identifier = "lognormal:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lognormal = cl
+ identifier = "trunc:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%trunc = cn
! read inversion domain settings
identifier = "w_edge_lon:"
@@ -532,12 +537,12 @@ module mod_settings
identifier = "yres:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%yres = real(cn)
- identifier = "nt_flx:"
+ identifier = "nstep_flx:"
call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) config%nt_flx = int(cn)
- identifier = "num_nee_day:"
+ if ( match ) config%nstep_flx = int(cn)
+ identifier = "nstep_nee:"
call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) config%num_nee_day = int(cn)
+ if ( match ) config%nstep_nee = int(cn)
identifier = "statres:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%statres = real(cn)
diff --git a/prep_syndata/mod_var.f90 b/prep_syndata/mod_var.f90
index 8a5f3d0cf56d8772d4c36c961f6fd48a0050c3bd..5d6e10b18c598cd2eccb2b76d9bf7b7c42878c4a 100644
--- a/prep_syndata/mod_var.f90
+++ b/prep_syndata/mod_var.f90
@@ -35,8 +35,7 @@ module mod_var
integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
- integer, parameter :: maxobs=10000 ! max number of observations
- real, parameter :: trunc=1.e-9 ! truncation threshold for eigenvalues of error covariance matrix
+ integer, parameter :: maxobs=50000 ! max number of observations
! flexpart variables
! ------------------
@@ -66,8 +65,9 @@ module mod_var
integer :: nday ! number of days in inversion interval
integer :: nyr ! number of years in inversion interval
integer :: nd_nee ! number of NEE values per day (in prior NEE file)
- integer :: nt_flx ! number of flux values per year in prior flux files
- integer :: nflx ! number of input flux time steps falling within inversion interval
+ integer :: nt_nee ! timestep of NEE fluxes (hours)
+ integer :: nt_flx ! timestep of other prior fluxes (hours)
+ integer :: nflx ! number of prior flux fields in inversion interval
integer :: ntstate ! number of time steps in state vector
real :: statres ! time resolution of the state vector (for CO2 1 day) (days)
real :: statres_hr ! sub-daily time resolution of the state vector (hours)
diff --git a/prep_syndata/prep_syndata.f90 b/prep_syndata/prep_syndata.f90
index 3e30d7d125183a3d5396645d45fa923eb21d187a..673a13dec0268cf3adbf5b42a68437f54f9b26db 100644
--- a/prep_syndata/prep_syndata.f90
+++ b/prep_syndata/prep_syndata.f90
@@ -79,7 +79,7 @@ program prep_syndata
call get_error_cov(files, config)
! get prior fluxes
- allocate( flx(nxgrid,nygrid,nflx))
+ allocate( flx(nxgrid,nygrid,nflx) )
allocate( state(nvar) )
flx(:,:,:) = 0.
state(:) = 0.
diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config
index 2b571e311fdae09d5f6a01a21721da13ff7b3ace..011642e1d65954c9f1372960119b3ab469094301 100644
--- a/settings/SETTINGS_co2_config
+++ b/settings/SETTINGS_co2_config
@@ -23,6 +23,13 @@ seed: 100
datei: 20120101
datef: 20120131
+# Use lognormal distribution (true or false) (only for non-CO2 species)
+# if true need to use m1qn3 method
+lognormal: .false.
+
+# Truncation of eigenvalues of B (cut at trunc x largest eigenvalue)
+trunc: 1.e-6
+
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: m1qn3
@@ -31,10 +38,10 @@ average_fp: .false.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
-maxiter: 2
+maxiter: 20
-# Include ocean boxes in inversion (true or false)
-# currently only for ghg species
+# Optimize ocean boxes (true or false)
+# currently only for ghg species (CO2 ocean fluxes not optimized)
inc_ocean: .false.
# Optimize initial mixing ratios (true or false)
@@ -54,11 +61,11 @@ prior_bg: 0
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
-restart: 1
+restart: 0
# Verbose output
# only use for debugging small runs
-verbose: .true.
+verbose: .false.
# Species ("co2" or "ghg")
spec: co2
@@ -92,43 +99,63 @@ regions: .false.
# State vector time resolution:
# time resolution at which NEE fluxes are optimized
-# statres determines averaging interval over 1 or more days (given in days)
-# statres_hr determines time resolution within 1 day (given in hours)
-# example: statres = 5, statres_hr = 6
-# gives 4 time intervals per day that are averaged over 5 days
-statres: 5
+# statres = averaging interval over 1 or more days (given in days)
+# statres_hr = time resolution within 1 day (given in hours)
+# e.g. statres = 10, statres_hr = 6 give 4 time intervals per day that are averaged over 10 days
+statres: 10
statres_hr: 12
-# non-NEE prior fluxes:
-# number of flux fields per year in input
-nt_flx: 12
-
# NEE prior fluxes:
-# num_nee_day = number of NEE values per 24h
-num_nee_day: 8
+# nstep_nee = time step of NEE fluxes (integer hours)
+nstep_nee: 3
+nstep_nee_reg: 3
+
+# Fossil fuel prior fluxes:
+# nstep_ff = time step of fossil fuel emissions (integer hours, for monthly data use 720)
+# coef_ff = coefficient to convert from input flux units to kg/m2/h
+nstep_ff: 720
+nstep_ff_reg: 720
+coeff_ff: 1.
+coeff_ff_reg: 1.
+
+# other prior fluxes:
+# nstep_ocn = time step of other fluxes (integer hours, for monthly data use 720)
+nstep_ocn: 720
# Measurement error: unit same as obs
# used if error in obs input <= zero
measerr: 0.5
-# Scalar of initial conc error: fraction
-cinierr: 0.01
+# Initial mixing ratio error: fraction
+cinierr: 0.005
# Prior flux error: fraction
flxerr: 0.25
# Fossil fuel error: fraction
+# only used to calculate uncertainty projected into observation space
ffferr: 0.05
# Spatial correlation length for land: unit (km)
-sigma_land: 20.
+sigma_land: 50.
# Spatial correlation length for ocean: unit (km)
-sigma_ocean: 1000.
+sigma_ocean: 2000.
# Temporal correlation length: unit (days)
-sigmatime: 7.
+sigmatime: 30.
# Total error for domain (Tg/y)
+# globerr <= 0: prior error covariance matrix not scaled
globerr: 300.
+# Settings for optimization of initial mixing ratios
+# comma separated list of northern edges of latitude bands
+cini_lat: -30.,0.,30.,90.
+# comma separated list of upper level of vertical bands (upper most level > outheight(nzgrid))
+#cini_alt: 2000., 10000., 55000.
+cini_alt: 55000.
+# time resolution for initial mixing ratio scalars (days)
+cinires: 30.
+
+
diff --git a/settings/SETTINGS_co2_files b/settings/SETTINGS_co2_files
index 520ac09fe5bff31cff1d164076f49d96263899c7..4af6a1f2485fe67ca43b71828122f33773561dcb 100644
--- a/settings/SETTINGS_co2_files
+++ b/settings/SETTINGS_co2_files
@@ -33,6 +33,7 @@ path_obs: /mypath/TEST_OUTPUT/OBS/CO2/
path_prior: /mypath/TEST_OUTPUT/FLUXES/CO2/
path_output: /mypath/TEST_OUTPUT/RESULTS/CO2/NO_NEST/
path_flexpart: /mypath/TEST_OUTPUT/FLEXOUT/CO2/NO_NEST/
+path_flexncdf: /mypath/TEST_OUTPUT/RESULTS/CO2/NO_NEST/
# Prior best guess file (if prior_bg = 1)
file_bg:
@@ -50,6 +51,13 @@ varname_regs: regions
lonname_regs: longitude
latname_regs: latitude
+# Orography file:
+# needs to correspond to global flexpart domain and resolution
+file_orog: /mypath/TEST_INPUT/elev.1-deg.nc
+varname_orog: data
+lonname_orog: lon
+latname_orog: lat
+
# Land-sea mask file:
# needs to correspond to flexpart resolution (for nested output needs to correspond to nested resolution)
file_lsm: /mypath/TEST_INPUT/lsm_1x1.nc
@@ -66,7 +74,7 @@ lonname_nee: longitude
latname_nee: latitude
timename_nee: time
-# Global prior fossil fuel flux file: flux kg/m2/h
+# Global prior fossil fuel (and bio-fuel/mass burning) flux file: flux kg/m2/h
# needs to correspond to flexpart domain and resolution
# generic name with year and month YYYYMM
filename_ff: CO2_FF_YYYYMM_10x10.nc
@@ -84,15 +92,6 @@ lonname_ocn: longitude
latname_ocn: latitude
timename_ocn: time
-# Global prior biomass burning flux file: flux kg/m2/h
-# needs to correspond to flexpart domain and resolution
-# generic name with year YYYY
-filename_bbg: CO2_BBG_YYYY_10x10.nc
-varname_bbg: emisco2
-lonname_bbg: longitude
-latname_bbg: latitude
-timename_bbg: time
-
# Regional prior NEE flux file: flux kg/m2/h
# only specify if using nested output (needs to correspond to nested resolution)
# generic name with year YYYY
@@ -102,7 +101,7 @@ lonnest_nee:
latnest_nee:
timenest_nee:
-# Regional prior fossil fuel flux file: flux kg/m2/h
+# Regional prior fossil fuel (and bio-fuel/mass burning) flux file: flux kg/m2/h
# only specify if using nested output (needs to correspond to nested resolution)
# generic name with year and month YYYYMM
filenest_ff:
@@ -120,21 +119,14 @@ lonnest_ocn:
latnest_ocn:
timenest_ocn:
-# Regional prior biomass burning flux file: flux kg/m2/h
-# only specify if using nested output (needs to correspond to nested resolution)
-# generic name with year YYYY
-filenest_bbg:
-varnest_bbg:
-lonnest_bbg:
-latnest_bbg:
-timenest_bbg:
-
-# Initial concentrations: conc ppm
+# Initial concentrations: conc (default is ppm)
# file_initconc given as generic name with year and month YYYYMM
# also need months preceding and proceeding inversion interval
path_initconc: /mypath/TEST_INPUT/INIT_CONC/CO2/
file_initconc: co2_noaa_YYYYMM.nc
varname_init: CO2
+# type of product (1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O and CO2, 4 = FP-CTM (daily), 5 = NOAA interpolated (monthly), 6 = FP-CTM (monthly))
+init_type: 5
# Receptor list file:
file_recept: /mypath/TEST_INPUT/reclist_co2.txt
diff --git a/settings/SETTINGS_co2_nest_config b/settings/SETTINGS_co2_nest_config
index f5237730ea0e6e5fd7cdb9714ff9e9961710985f..54ca08a3b3d15be5999451898ecacb522c0abca9 100644
--- a/settings/SETTINGS_co2_nest_config
+++ b/settings/SETTINGS_co2_nest_config
@@ -23,18 +23,25 @@ seed: 100
datei: 20120101
datef: 20120131
+# Use lognormal distribution (true or false) (only for non-CO2 species)
+# if true need to use m1qn3 method
+lognormal: .false.
+
+# Truncation of eigenvalues of B (cut at trunc x largest eigenvalue)
+trunc: 1.e-6
+
# Inversion method ('analytic', 'congrad' or 'm1qn3')
-method: analytic
+method: m1qn3
# Average/select flexpart releases (true or false)
average_fp: .false.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
-maxiter: 2
+maxiter: 20
-# Include ocean boxes in inversion (true or false)
-# currently only for ghg species
+# Optimize ocean boxes (true or false)
+# currently only for ghg species (CO2 ocean fluxes not optimized)
inc_ocean: .false.
# Optimize initial mixing ratios (true or false)
@@ -54,7 +61,7 @@ prior_bg: 0
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
-restart: 1
+restart: 0
# Verbose output
# only use for debugging small runs
@@ -94,42 +101,61 @@ regions: .true.
# time resolution at which NEE fluxes are optimized
# statres determines averaging interval over 1 or more days (given in days)
# statres_hr determines time resolution within 1 day (given in hours)
-# example: statres = 5, statres_hr = 6
-# gives 4 time intervals per day that are averaged over 5 days
-statres: 5
+# e.g. statres = 10, statres_hr = 6 give 4 time intervals per day that are averaged over 10 days
+statres: 10
statres_hr: 12
-# non-NEE prior fluxes:
-# number of flux fields per year in input
-nt_flx: 12
-
# NEE prior fluxes:
-# num_nee_day = number of NEE values per 24h
-num_nee_day: 8
+# nstep_nee = time step of NEE fluxes (integer hours)
+nstep_nee: 3
+nstep_nee_reg: 3
+
+# Fossil fuel prior fluxes:
+# nstep_ff = time step of fossil fuel emissions (integer hours, for monthly data use 720)
+# coef_ff = coefficient to convert from input flux units to kg/m2/h
+nstep_ff: 720
+nstep_ff_reg: 720
+coeff_ff: 1.
+coeff_ff_reg: 1.
+
+# Ocean prior fluxes:
+# nstep_ocn = time step of other fluxes (integer hours, for monthly data use 720)
+nstep_ocn: 720
# Measurement error: unit same as obs
# used if error in obs input <= zero
measerr: 0.5
-# Scalar of initial conc error: fraction
-cinierr: 0.01
+# Initial mixing ratio error: fraction
+cinierr: 0.005
# Prior flux error: fraction
flxerr: 0.25
# Fossil fuel error: fraction
+# only used to calculate uncertainty projected into observation space
ffferr: 0.05
# Spatial correlation length for land: unit (km)
-sigma_land: 20.
+sigma_land: 50.
# Spatial correlation length for ocean: unit (km)
-sigma_ocean: 1000.
+sigma_ocean: 2000.
# Temporal correlation length: unit (days)
-sigmatime: 7.
+sigmatime: 30.
# Total error for domain (Tg/y)
# globerr <= 0: prior error covariance matrix not scaled
globerr: 300.
+# Settings for optimization of initial mixing ratios
+# comma separated list of northern edges of latitude bands
+cini_lat: -30.,0.,30.,90.
+# comma separated list of upper level of vertical bands (upper most level > outheight(nzgrid))
+#cini_alt: 2000., 10000., 55000.
+cini_alt: 55000.
+# time resolution for initial mixing ratio scalars (days)
+cinires: 30.
+
+
diff --git a/settings/SETTINGS_co2_nest_files b/settings/SETTINGS_co2_nest_files
index c4a5342494a1ed7b55bec00f8326b570d3d7bbbc..2420049355f3d3bad311ff9d4958d788323c9b54 100644
--- a/settings/SETTINGS_co2_nest_files
+++ b/settings/SETTINGS_co2_nest_files
@@ -50,6 +50,13 @@ varname_regs: regions
lonname_regs: longitude
latname_regs: latitude
+# Orography file:
+# needs to correspond to global flexpart domain and resolution
+file_orog: /mypath/TEST_INPUT/elev.2-deg.nc
+varname_orog: data
+lonname_orog: lon
+latname_orog: lat
+
# Land-sea mask file:
# needs to correspond to flexpart resolution (for nested output needs to correspond to nested resolution)
file_lsm: /mypath/TEST_INPUT/lsm_0.5x0.5.nc
@@ -66,7 +73,7 @@ lonname_nee: longitude
latname_nee: latitude
timename_nee: time
-# Global prior fossil fuel flux file: flux kg/m2/h
+# Global prior fossil fuel (and bio-fuel/mass burning) flux file: flux kg/m2/h
# needs to correspond to flexpart domain and resolution
# generic name with year and month YYYYMM
filename_ff: CO2_FF_YYYYMM_20x20.nc
@@ -84,15 +91,6 @@ lonname_ocn: longitude
latname_ocn: latitude
timename_ocn: time
-# Global prior biomass burning flux file: flux kg/m2/h
-# needs to correspond to flexpart domain and resolution
-# generic name with year YYYY
-filename_bbg: CO2_BBG_YYYY_20x20.nc
-varname_bbg: emisco2
-lonname_bbg: longitude
-latname_bbg: latitude
-timename_bbg: time
-
# Regional prior NEE flux file: flux kg/m2/h
# only specify if using nested output (needs to correspond to nested resolution)
# generic name with year YYYY
@@ -102,7 +100,7 @@ lonnest_nee: longitude
latnest_nee: latitude
timenest_nee: time
-# Regional prior fossil fuel flux file: flux kg/m2/h
+# Regional prior fossil fuel (and bio-fuel/mass burning) flux file: flux kg/m2/h
# only specify if using nested output (needs to correspond to nested resolution)
# generic name with year and month YYYYMM
filenest_ff: CO2_FF_YYYYMM_05x05.nc
@@ -120,21 +118,14 @@ lonnest_ocn: longitude
latnest_ocn: latitude
timenest_ocn: time
-# Regional prior biomass burning flux file: flux kg/m2/h
-# only specify if using nested output (needs to correspond to nested resolution)
-# generic name with year YYYY
-filenest_bbg: CO2_BBG_YYYY_05x05.nc
-varnest_bbg: emisco2
-lonnest_bbg: longitude
-latnest_bbg: latitude
-timenest_bbg: time
-
-# Initial concentrations: conc ppm
+# Initial concentrations: conc (default is ppm)
# file_initconc given as generic name with year and month YYYYMM
# also need months preceding and proceeding inversion interval
path_initconc: /mypath/TEST_INPUT/INIT_CONC/CO2/
file_initconc: co2_noaa_YYYYMM.nc
varname_init: CO2
+# type of product (1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O and CO2, 4 = FP-CTM (daily), 5 = NOAA interpolated (monthly), 6 = FP-CTM (monthly))
+init_type: 5
# Receptor list file:
file_recept: /mypath/TEST_INPUT/reclist_co2.txt
diff --git a/settings/SETTINGS_ghg_config b/settings/SETTINGS_ghg_config
index 05349ab2d0ce59f7e3aba29bc87d80e9df1e5c9d..f684ff0963ee8d6cc01e0e245739cd3289f0d0fa 100644
--- a/settings/SETTINGS_ghg_config
+++ b/settings/SETTINGS_ghg_config
@@ -23,21 +23,24 @@ seed: 100
datei: 20120101
datef: 20120131
-# Inversion method ('analytic', 'congrad' or 'm1qn3')
-method: analytic
-
# Use lognormal distribution (true or false) (only for non-CO2 species)
# if true need to use m1qn3 method
-lognormal: .false.
+lognormal: .false.
+
+# Truncation of eigenvalues of B (cut at trunc x largest eigenvalue)
+trunc: 1.e-4
+
+# Inversion method ('analytic', 'congrad' or 'm1qn3')
+method: analytic
# Average/select flexpart releases (true or false)
-average_fp: .false.
+average_fp: .true.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
-# Include ocean boxes in the optimization (true or false)
+# Optimize ocean boxes (true or false)
inc_ocean: .true.
# Optimize initial mixing ratios (true or false)
@@ -61,7 +64,7 @@ restart: 0
# Verbose output
# only use for debugging small runs
-verbose: .false.
+verbose: .true.
# Species ("co2" or "ghg")
spec: ghg
@@ -98,16 +101,15 @@ regions: .true.
# must be an integer multiple of 1 day
statres: 10.
-# non-NEE prior fluxes:
-# number of flux fields per year in input
-nt_flx: 12
+# prior fluxes:
+# nstep_flx = time step of other fluxes (integer hours, for monthly data use 720)
+nstep_flx: 720
# Measurement error: unit same as obs
-# used if error in obs input <= zero
measerr: 5.0
-# Scalar of initial conc error: fraction
-cinierr: 0.01
+# Initial mixing ratio error: fraction
+cinierr: 0.005
# Prior flux error: fraction
flxerr: 0.5
@@ -125,5 +127,16 @@ sigma_ocean: 1000.
sigmatime: 30.
# Total error for inversion domain (Tg/y)
+# globerr <= 0: prior error covariance matrix not scaled
globerr: 10.
+# Settings for optimization of initial mixing ratios
+# comma separated list of northern edges of latitude bands
+cini_lat: -30.,0.,30.,90.
+# comma separated list of upper level of vertical bands (upper most level > outheight(nzgrid))
+#cini_alt: 2000., 10000., 55000.
+cini_alt: 55000.
+# time resolution for initial mixing ratio scalars (days)
+cinires: 30.
+
+
diff --git a/settings/SETTINGS_ghg_files b/settings/SETTINGS_ghg_files
index d143d969cf7516afcac3e3bcab78ee107469f182..663abb952c1580610a0bd1798b365a81bf08f7ea 100644
--- a/settings/SETTINGS_ghg_files
+++ b/settings/SETTINGS_ghg_files
@@ -33,6 +33,7 @@ path_obs: /mypath/TEST_OUTPUT/OBS/GHG/
path_prior: /mypath/TEST_OUTPUT/FLUXES/GHG/
path_output: /mypath/TEST_OUTPUT/RESULTS/GHG/NO_NEST/congrad/
path_flexpart: /mypath/TEST_OUTPUT/FLEXOUT/GHG/NO_NEST/
+path_flexncdf: /mypath/TEST_OUTPUT/RESULTS/GHG/NO_NEST/congrad/
# Prior best guess file (if prior_bg = 1 in SETTINGS_config)
file_bg:
@@ -50,6 +51,13 @@ varname_regs: regions
lonname_regs: longitude
latname_regs: latitude
+# Orography file:
+# needs to correspond to global flexpart domain and resolution
+file_orog: /xnilu_wrk/users/rlt/LANDCOVER/elev.1-deg.nc
+varname_orog: data
+lonname_orog: lon
+latname_orog: lat
+
# Land-sea mask file:
# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
file_lsm: /mypath/TEST_INPUT/lsm_1x1.nc
@@ -57,14 +65,6 @@ varname_lsm: lsm
lonname_lsm: lon
latname_lsm: lat
-# Snow-ice mask file:
-# optional: only used to define aggregated grid in "prep_regions"
-# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
-#file_ice:
-#varname_ice:
-#lonname_ice:
-#latname_ice:
-
# Global prior flux file: flux kg/m2/h
# needs to correspond to flexpart domain and resolution
# generic name with year YYYY
@@ -88,6 +88,8 @@ timenest_flx:
path_initconc: /mypath/TEST_INPUT/INIT_CONC/GHG/
file_initconc: ch4_noaa_YYYYMM.nc
varname_init: CH4
+# type of product (1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O, 4 = FP-CTM (daily), 5 = NOAA interpolated (monthly), 6 = FP-CTM (monthly))
+init_type: 5
# Receptor list file:
file_recept: /mypath/TEST_INPUT/reclist_ghg.txt
diff --git a/settings/SETTINGS_ghg_nest_config b/settings/SETTINGS_ghg_nest_config
index edbef48446ae44dce8578c7a04cd54d4157c4dde..ca532816da452cd072591b09bab18c5959a7b17d 100644
--- a/settings/SETTINGS_ghg_nest_config
+++ b/settings/SETTINGS_ghg_nest_config
@@ -23,21 +23,24 @@ seed: 100
datei: 20120101
datef: 20120131
-# Inversion method ('analytic', 'congrad' or 'm1qn3')
-method: congrad
-
# Use lognormal distribution (true or false) (only for non-CO2 species)
# if true need to use m1qn3 method
-lognormal: .false.
+lognormal: .false.
+
+# Truncation of eigenvalues of B (cut at trunc x largest eigenvalue)
+trunc: 1.e-4
+
+# Inversion method ('analytic', 'congrad' or 'm1qn3')
+method: congrad
# Average/select flexpart releases (true or false)
-average_fp: .false.
+average_fp: .true.
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
-maxiter: 2
+maxiter: 10
-# Include ocean boxes in the optimization (true or false)
+# Optimize ocean boxes (true or false)
inc_ocean: .true.
# Optimize initial mixing ratios (true or false)
@@ -57,11 +60,11 @@ prior_bg: 0
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
-restart: 1
+restart: 0
# Verbose output
# only use for debugging small runs
-verbose: .true.
+verbose: .false.
# Species ("co2" or "ghg")
spec: ghg
@@ -98,16 +101,16 @@ regions: .true.
# must be an integer multiple of 1 day
statres: 10.
-# non-NEE prior fluxes:
-# number of flux fields per year in input
-nt_flx: 12
+# prior fluxes:
+# nstep_flx = time step of other fluxes (integer hours, for monthly data use 720)
+nstep_flx: 720
# Measurement error: unit same as obs
# used if error in obs input <= zero
measerr: 5.0
-# Scalar of initial conc error: fraction
-cinierr: 0.01
+# Initial mixing ratio error: fraction
+cinierr: 0.005
# Prior flux error: fraction
flxerr: 0.5
@@ -128,3 +131,11 @@ sigmatime: 30.
# globerr <= 0: prior error covariance matrix not scaled
globerr: 10.
+# Settings for optimization of initial mixing ratios
+# comma separated list of northern edges of latitude bands
+cini_lat: -30.,0.,30.,90.
+# comma separated list of upper level of vertical bands (upper most level > outheight(nzgrid))
+#cini_alt: 2000., 10000., 55000.
+cini_alt: 55000.
+# time resolution for initial mixing ratio scalars (days)
+cinires: 30.
diff --git a/settings/SETTINGS_ghg_nest_files b/settings/SETTINGS_ghg_nest_files
index 92e0700e8cf5c426d6596bceae4c330009851f1a..9d1574b00c5fc46ddf60ba36d6948a3f06ffc677 100644
--- a/settings/SETTINGS_ghg_nest_files
+++ b/settings/SETTINGS_ghg_nest_files
@@ -33,6 +33,7 @@ path_obs: /mypath/TEST_OUTPUT/OBS/GHG/
path_prior: /mypath/TEST_OUTPUT/FLUXES/GHG/
path_output: /mypath/TEST_OUTPUT/RESULTS/GHG/NEST/
path_flexpart: /mypath/TEST_OUTPUT/FLEXOUT/GHG/NEST/
+path_flexncdf: /mypath/TEST_OUTPUT/RESULTS/GHG/NEST/
# Prior best guess file (if prior_bg = 1 in SETTINGS_config)
file_bg:
@@ -50,6 +51,13 @@ varname_regs: regions
lonname_regs: longitude
latname_regs: latitude
+# Orography file:
+# needs to correspond to global flexpart domain and resolution
+file_orog: /mypath/TEST_OUTPUT/elev.2-deg.nc
+varname_orog: data
+lonname_orog: lon
+latname_orog: lat
+
# Land-sea mask file:
# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
file_lsm: /mypath/TEST_INPUT/lsm_0.5x0.5.nc
@@ -57,14 +65,6 @@ varname_lsm: lsm
lonname_lsm: lon
latname_lsm: lat
-# Snow-ice mask file:
-# optional: if specified it is only used to define aggregated grid in "prep_regions"
-# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
-#file_ice:
-#varname_lsm:
-#lonname_lsm:
-#latname_lsm:
-
# Global prior flux file: flux kg/m2/h
# needs to correspond to flexpart domain and resolution
# generic name with year YYYY
@@ -88,6 +88,8 @@ timenest_flx: time
path_initconc: /mypath/TEST_INPUT/INIT_CONC/GHG/
file_initconc: ch4_noaa_YYYYMM.nc
varname_init: CH4
+# type of product (1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O, 4 = FP-CTM (daily), 5 = NOAA interpolated (monthly), 6 = FP-CTM (monthly))
+init_type: 5
# Receptor list file:
file_recept: /mypath/TEST_INPUT/reclist_ghg.txt
diff --git a/source/average_fp.f90 b/source/average_fp.f90
index 152a236140068c5756bdf04844c35e30e9dab0e7..fc854fafd0e3a89ec321ae9c045ac9e952952f6b 100755
--- a/source/average_fp.f90
+++ b/source/average_fp.f90
@@ -105,6 +105,13 @@ subroutine average_fp(files, config, obs)
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(i))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'header'
+ ! check if netcdf file for this observation already exists
+ write(areldate,'(I8.8)') jjjjmmdd
+ write(areltime,'(I6.6)') hhmiss
+ inquire ( file=trim(files%path_flexncdf)//trim(obs%recs(i))//'/'//trim(adate)//'/'//&
+ 'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc', exist=lexist )
+ if ( lexist ) go to 10
+
! Get header and release information related to current observation location
! (if not read and same for previous observation)
if ( i.eq.1 .or. month.ne.prevmonth ) then
@@ -150,6 +157,7 @@ subroutine average_fp(files, config, obs)
gtime_tmp = gtime ! Save the output date of each field in footprint grid
grid_tmp = grid ! Save the footprint to temporary grid for summation
nr_footprints = 1 ! Start counting number of footprints to be averaged
+ print*, 'average_fp: gtime = ', gtime
! Same for nested grid if applicable
if ( config%nested ) then
@@ -216,6 +224,7 @@ subroutine average_fp(files, config, obs)
grid = grid_tmp/real(nr_footprints)
gridnest = gridnest_tmp/real(nr_footprints)
gtime = gtime_tmp
+ print*, 'average_fp: gtime after average = ',gtime
! Save as netcdf
! output filename same as start date and time of observation
@@ -223,10 +232,10 @@ subroutine average_fp(files, config, obs)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
- ! create subdirectory in path_output for files
+ ! create subdirectory in path_flexncdf for files
subdir_fp = trim(obs%recs(i))//'/'//trim(adate)//'/'
- call system('mkdir -p '//trim(files%path_output)//trim(obs%recs(i)))
- call system('mkdir -p '//trim(files%path_output)//trim(subdir_fp))
+ call system('mkdir -p '//trim(files%path_flexncdf)//trim(obs%recs(i)))
+ call system('mkdir -p '//trim(files%path_flexncdf)//trim(subdir_fp))
filename_out = trim(subdir_fp)//'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc'
call save_average_fp(files, filename_out, nxgrid, nygrid, maxngrid, ngrid, grid, gbl_lon, gbl_lat, gtime)
! If applicable same for nested grid
diff --git a/source/error_cov.f90 b/source/error_cov.f90
index d8a8886340c15db51b254db8677392f704410906..62c201430f931ae279eee94466aab4e0be522072 100644
--- a/source/error_cov.f90
+++ b/source/error_cov.f90
@@ -131,8 +131,8 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
! truncate eigenvalues
! eigenvalues are sorted in ascending order
- ind = minloc(evals, dim=1, mask=evals.ge.(trunc*maxval(evals)))
- write(logid,*) 'Truncating eigenvalues at: ',trunc*maxval(evals)
+ 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
write(logid,*) 'Number of retained eigenvalues: ',neig
@@ -252,7 +252,8 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
if ( config%opt_cini ) then
do i = 1, ntcini*ncini
if ( config%lognormal ) then
- states%pxerr0(npvar+i) = log(max(smallnum,config%cinierr))
+ ! convert to error of z
+ states%pxerr0(npvar+i) = log(1.+config%cinierr)
else
states%pxerr0(npvar+i) = config%cinierr
endif
diff --git a/source/get_initcams.f90 b/source/get_initcams.f90
deleted file mode 100644
index b89a396371a456d102cf43fee0275e7dbdde21ac..0000000000000000000000000000000000000000
--- a/source/get_initcams.f90
+++ /dev/null
@@ -1,260 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: get_initcams
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> get_initcams
-!!
-!! Purpose: Reads the initial concentration files and interpolates these to
-!! the flexpart grid
-!!
-!! Interface:
-!!
-!! Inputs
-!! filename - name of 4D concentration file
-!! varname - variable name
-!! concini - concentrations on flexpart grid (zero array on input)
-!!
-!! Externals
-!! convertgrid
-!! interpz
-!! check
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine get_initcams(filename, varname, ntime, concini)
-
- use mod_var
- use mod_tools, only : sort
- use mod_save, only : check
- use netcdf
-
- implicit none
-
- character(len=max_path_len), intent(in) :: filename
- character(len=max_name_len), intent(in) :: varname
- integer, intent(in) :: ntime
- real, dimension(nxgrid,nygrid,nzgrid,ntime), intent(in out) :: concini
-
- character(len=*), parameter :: lonname='longitude'
- character(len=*), parameter :: latname='latitude'
- character(len=*), parameter :: timename='time'
- character(len=*), parameter :: prsname='hlevel' ! hlevel for CH4, level for N2O
- logical, parameter :: midpoints=.true. ! true for CH4, false for N2O
- integer :: ncid, xid, yid, zid, tid, apid, bpid, varid
- integer :: n, kz, status
- integer :: nx, ny, nz, nt, nlev, nzmax
- integer :: nx_out, ny_out
- logical :: flag_hyb
- real, dimension(:), allocatable :: lon, lat, lonin, latin
- real(kind=8), dimension(:), allocatable :: lono, lato
- real, dimension(:), allocatable :: ap, bp, pres, mpres, alt
- integer, dimension(:), allocatable :: indx, indy
- real, dimension(nxgrid) :: lon_out
- real, dimension(nygrid) :: lat_out
- real, dimension(:,:,:), allocatable :: work1, work2, work3, work4
- real, dimension(:,:,:,:), allocatable :: conc, concin
- real :: sclfact, offset
-
- ! initialize
- ! ----------
-
- ! don't pass globals to subroutines
- lon_out = gbl_lon
- lat_out = gbl_lat
- nx_out = nxgrid
- ny_out = nygrid
-
- ! read netcdf file
- ! ----------------
-
- write(logid,*) 'Reading file: '//trim(filename)
-
- ! read dimensions
- call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
- call check( nf90_inq_dimid(ncid,lonname,xid) )
- call check( nf90_inq_dimid(ncid,latname,yid) )
- call check( nf90_inq_dimid(ncid,prsname,zid) )
- call check( nf90_inq_dimid(ncid,timename,tid) )
- call check( nf90_inquire_dimension(ncid,xid,len=nx) )
- call check( nf90_inquire_dimension(ncid,yid,len=ny) )
- call check( nf90_inquire_dimension(ncid,zid,len=nlev) )
- call check( nf90_inquire_dimension(ncid,tid,len=nt) )
-
- if ( ntime.ne.nt ) write(logid,*) 'ERROR get_initcams: inconsistent time dimension'
-
- allocate( lon(nx) )
- allocate( lat(ny) )
-
- ! read variables
- call check( nf90_inq_varid(ncid,lonname,xid) )
- call check( nf90_get_var(ncid,xid,lon) )
- call check( nf90_inq_varid(ncid,latname,yid) )
- call check( nf90_get_var(ncid,yid,lat) )
- ! check if contains hybrid pressure coordinates
- status = nf90_inq_varid(ncid,'hyai',apid)
- print*, 'get_initcams: status = ',status
- if (status.eq.-49) then
- ! no hybrid pressure coordinates
- flag_hyb = .false.
- nz = nlev
- allocate( mpres(nz) )
- call check( nf90_inq_varid(ncid,prsname,zid) )
- call check( nf90_get_var(ncid,zid,mpres) )
- else
- ! contains hybrid pressure coordinates
- flag_hyb = .true.
- nz = nlev - 1
- allocate( ap(nlev) )
- allocate( bp(nlev) )
- allocate( pres(nlev) )
- allocate( mpres(nz) )
- call check( nf90_inq_varid(ncid,'hyai',apid) )
- call check( nf90_get_var(ncid,apid,ap) )
- call check( nf90_inq_varid(ncid,'hybi',bpid) )
- call check( nf90_get_var(ncid,bpid,bp) )
- endif
- allocate( alt(nz) )
- allocate( conc(nx,ny,nz,nt) )
- call check( nf90_inq_varid(ncid,trim(varname),varid) )
- call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
- ! check if file contains scalar and offset
- status = nf90_inquire_attribute(ncid,varid,'scale_factor')
- if ( status.eq.nf90_noerr ) then
- status = nf90_inquire_attribute(ncid,varid,'add_offset')
- if ( status.eq.nf90_noerr ) then
- call check( nf90_get_att(ncid, varid, 'scale_factor', sclfact) )
- call check( nf90_get_att(ncid, varid, 'add_offset', offset) )
- conc = conc*sclfact + offset
- endif
- endif
- call check( nf90_close(ncid) )
-
- ! adjust lat and lon of grid
- allocate( indy(ny) )
- allocate( indx(nx) )
- allocate( lato(ny) )
- allocate( lono(nx) )
- lato = dble(lat)
- lono = dble(lon)
- ! sort from south to north and west to east
- call sort(ny, lato, indy)
- call sort(nx, lono, indx)
- conc = conc(indx,indy,:,:)
- print*, 'get_initcams: indy = ',indy
- print*, 'get_initcams: indx = ',indx
- ! cut extra longitude grid cells (assumes grid is -180 to 180 and not 0 to 360)
- if ( lono(nx).ge.180. ) then
- allocate( concin((nx-1),ny,nz,nt) )
- allocate( lonin((nx-1)) )
- concin = conc(1:(nx-1),:,:,:)
- lonin = sngl(lono(1:(nx-1)))
- else
- allocate( concin(nx,ny,nz,nt) )
- allocate( lonin(nx) )
- concin = conc(:,:,:,:)
- lonin = sngl(lono(:))
- endif
- allocate( latin(ny) )
- if ( lato(ny).eq.90. ) then
- do n = 1, ny
- latin(n) = -90. + real(n - 1)*(180./real(ny))
- end do
- else
- latin = sngl(lato)
- endif
- print*, 'get_initcams: latin = ',latin
- print*, 'get_initcams: lonin = ',lonin
-
- ! calculate mid layer pressure
- ! (for simplicity currently use standard surface pressure
- ! should be good enough for background calculation,
- ! future update would be to read surf pres from file and do
- ! vertical interpolation to flexpart levels grid by grid)
- if ( flag_hyb ) then
- do n = 1, nlev
- pres(n) = ap(n) + bp(n)*psurf
- end do
- ! avoid pressure of zero
- pres(nlev) = 2.
- do n = 1, nz
- mpres(n) = 0.5*(pres(n) + pres(n+1))
- end do
- endif
-
- do n = 1, nz
- alt(n)=7500*log(psurf/mpres(n))
- end do
- if ( flag_hyb ) print*, 'get_initcams: pres = ',pres
- print*, 'get_initcams: mpres = ',mpres
- print*, 'get_initcams: alt = ',alt
-
- ! interpolate to flexpart grid
- ! ----------------------------
-
- if ( lono(nx).ge.180. ) then
- allocate( work1((nx-1),ny,1) )
- nx = nx - 1
- else
- allocate( work1(nx,ny,1) )
- endif
-
- ! determine upper alt level of conc field for vertical interpolation
- nzmax = minloc(alt, dim=1, mask=abs(alt-outheight(nzgrid)).eq.minval(abs(alt-outheight(nzgrid))))
- nzmax = min(nz, nzmax+1)
- print*, 'get_initcams: nzmax = ',nzmax
-
- allocate( work2(nxgrid,nygrid,1) )
- allocate( work3(nxgrid,nygrid,nzmax) )
- allocate( work4(nxgrid,nygrid,nzgrid) )
-
- do n = 1, nt
- ! interpolate horizontally
- do kz = 1, nzmax
- work1(:,:,1) = concin(:,:,kz,n)
- ! gbl_lon and gbl_lat are lon and lat of global flexpart domain
- call convertgrid(midpoints, nx, ny, lonin, latin, work1, nx_out, ny_out, lon_out, lat_out, work2)
- work3(:,:,kz) = work2(:,:,1)
- end do
- ! interpolate vertically
- call interpz(nzmax, alt(1:nzmax), work3, work4)
- concini(:,:,:,n) = work4
- end do
-
- deallocate( lon )
- deallocate( lat )
- deallocate( lono )
- deallocate( lato )
- deallocate( lonin )
- deallocate( latin )
- deallocate( indx )
- deallocate( indy )
- if ( allocated(ap) ) deallocate( ap )
- if ( allocated(bp) ) deallocate( bp )
- if ( allocated(pres) ) deallocate( pres )
- deallocate( mpres )
- deallocate( alt )
- deallocate( conc )
- deallocate( concin )
- deallocate( work1 )
- deallocate( work2 )
- deallocate( work3 )
- deallocate( work4 )
-
-end subroutine get_initcams
-
-
diff --git a/source/get_initconc.f90 b/source/get_initconc.f90
index d06ac3eb2167c20e572fe11f734017d742edf4c0..35dc1ab7605c05923718468d2f300aab8f910e2a 100644
--- a/source/get_initconc.f90
+++ b/source/get_initconc.f90
@@ -87,7 +87,7 @@ subroutine get_initconc(files, config, filename, ntime, concini)
ny_out = nygrid
! which initial concentration field product
- ! 1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O, 4 = FP-CTM, 5 = NOAA interpolated, 6 = FP-CTM (monthly)
+ ! 1 = CAMS CH4, 2 = EGG4, 3 = CAMS N2O and CO2, 4 = FP-CTM, 5 = NOAA interpolated, 6 = FP-CTM (monthly)
iproduct = int(files%init_type)
varname = files%varname_init
@@ -105,7 +105,7 @@ subroutine get_initconc(files, config, filename, ntime, concini)
timename = 'time'
prsname = 'level'
case (3)
- ! CAMS N2O
+ ! CAMS N2O and CO2
lonname = 'longitude'
latname = 'latitude'
timename = 'time'
@@ -206,11 +206,25 @@ subroutine get_initconc(files, config, filename, ntime, concini)
! convert to unit of Pa
mpres = mpres*100.
case (3)
- ! CAMS N2O
- ! vertical coordinate is pressure
- allocate( mpres(nz) )
- call check( nf90_inq_varid(ncid,trim(prsname),varid) )
- call check( nf90_get_var(ncid,varid,mpres) )
+ ! CAMS N2O and CO2
+ if ( config%spec.eq.'n2o' ) then
+ ! vertical coordinate is pressure
+ allocate( mpres(nz) )
+ call check( nf90_inq_varid(ncid,trim(prsname),varid) )
+ call check( nf90_get_var(ncid,varid,mpres) )
+ else if ( config%spec.eq.'co2' ) then
+ ! hybrid pressure coordinates
+ allocate( ap(nz+1) )
+ allocate( bp(nz+1) )
+ call check( nf90_inq_varid(ncid,'ap',varid) )
+ call check( nf90_get_var(ncid,varid,ap) )
+ call check( nf90_inq_varid(ncid,'bp',varid) )
+ call check( nf90_get_var(ncid,varid,bp) )
+ ! read surface pressure
+ allocate( surfp(nx,ny,nt) )
+ call check( nf90_inq_varid(ncid,'Psurf',varid) )
+ call check( nf90_get_var(ncid,varid,surfp) )
+ endif
case (4)
! FP-CTM
allocate( mpres(nz) )
@@ -260,8 +274,14 @@ subroutine get_initconc(files, config, filename, ntime, concini)
deallocate( pw )
case (2)
! EGG4
- ! adjust units from mass mixing ratio to mole fraction
+ ! adjust units from mass mixing ratio to ppb
conc = conc*1.e9*mmair/config%molarmass
+ case (3)
+ ! CAMS N2O and CO2
+ if ( config%spec.eq.'co2' ) then
+ ! adjust units from mole fraction to ppm
+ conc = conc*1.e6
+ endif
end select
! calculate geopotential height
@@ -292,11 +312,27 @@ subroutine get_initconc(files, config, filename, ntime, concini)
end do
zkind = 0
case (3)
- ! CAMS N2O
- do n = 1, nz
- ! at present temp not included so use estimate
- geoh(:,:,n) = (gasc*ts/gc/mmol)*log(psurf/mpres(n))
- end do
+ ! CAMS N2O and CO2
+ if ( config%spec.eq.'n2o' ) then
+ do n = 1, nz
+ ! at present temp not included so use estimate
+ geoh(:,:,n) = (gasc*ts/gc/mmol)*log(psurf/mpres(n))
+ end do
+ else if ( config%spec.eq.'co2' ) then
+ ! calculate midlayer pressure
+ allocate( pres(nx,ny,nz) )
+ allocate( rpres(nx,ny) )
+ do n = 1, nz
+ pres(:,:,n) = 0.5*(ap(n)+ap(n+1)) + 0.5*(bp(n)+bp(n+1))*sum(surfp(:,:,:),dim=3)/real(nt)
+ end do
+ do n = 1, nz
+ ! at present temp not included so use estimate
+ rpres = sum(surfp(:,:,:),dim=3)/real(nt)/pres(:,:,n)
+ geoh(:,:,n) = (gasc*ts/gc/mmol)*log(rpres)
+ end do
+ endif
+ deallocate( rpres )
+ deallocate( surfp )
zkind = 0
case (4)
! FP-CTM
@@ -339,7 +375,7 @@ subroutine get_initconc(files, config, filename, ntime, concini)
ny = ny - 1
midpoints = .true.
case (3)
- ! CAMS N2O
+ ! CAMS N2O and CO2
! file contains extra longitude
nx = nx - 1
! latitude given from 90 to -90
diff --git a/source/get_initfpctm.f90 b/source/get_initfpctm.f90
deleted file mode 100644
index 3bdc173861db1c93874f4fb7faf4fcce34a110b8..0000000000000000000000000000000000000000
--- a/source/get_initfpctm.f90
+++ /dev/null
@@ -1,185 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: get_initfpctm
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> get_initfpctm
-!!
-!! Purpose: Reads the initial concentration files and interpolates these to
-!! the flexpart grid
-!!
-!! Interface:
-!!
-!! Inputs
-!! filename - name of 4D concentration file
-!! varname - variable name
-!! concini - concentrations on flexpart grid (zero array on input)
-!!
-!! Externals
-!! convertgrid
-!! interpz
-!! check
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine get_initfpctm(filename, varname, concini)
-
- use mod_var
- use mod_save, only : check
- use netcdf
-
- implicit none
-
- character(len=max_path_len), intent(in) :: filename
- character(len=max_name_len), intent(in) :: varname
- real, dimension(nxgrid,nygrid,nzgrid), intent(in out) :: concini
-
- character(len=*), parameter :: lonname='longitude'
- character(len=*), parameter :: latname='latitude'
- character(len=*), parameter :: prsname='pressure'
- character(len=*), parameter :: levname='lev'
- character(len=*), parameter :: timename='time'
- integer :: ncid, xid, yid, zid, tid, varid
- integer :: kz
- integer :: nx, ny, nz, nt
- integer :: nx_out, ny_out
- logical :: midpoints
- real, dimension(:), allocatable :: lon, lat, prs, geoh
- real, dimension(nxgrid) :: lon_out
- real, dimension(nygrid) :: lat_out
- real, dimension(:,:,:), allocatable :: work, workout, concout
- real, dimension(:,:,:,:), allocatable :: conc
- integer :: status
- logical :: preslevel
-
- ! initialize
- ! ----------
-
- ! don't pass globals to subroutines
- lon_out = gbl_lon
- lat_out = gbl_lat
- nx_out = nxgrid
- ny_out = nygrid
-
- ! read netcdf file
- ! ----------------
-
- write(logid,*) 'Reading file: '//trim(filename)
- preslevel = .true.
- midpoints = .false.
-
- ! read dimensions
- call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
- call check( nf90_inq_dimid(ncid,lonname,xid) )
- call check( nf90_inq_dimid(ncid,latname,yid) )
-
- ! the data may be on pressure or z-levels, check this (info lost when using "check")
- status = nf90_inq_dimid(ncid,prsname,zid)
- if (status.eq.-46) then
- ! error code -46 is given when the dimension name does not exist, try if level rather than pressure is available.
- write(logid,*) 'WARNING get_initfpctm: pressure not available, test if initial concentration on z-level.'
- status = nf90_inq_dimid(ncid,levname,zid)
- call check( nf90_inq_dimid(ncid,levname,zid) )
- preslevel = .false.
- midpoints = .true.
- else
- call check( nf90_inq_dimid(ncid,prsname,zid) )
- endif
-
- ! time does not have to be available, check this here (info lost when using "check")
- status=nf90_inq_dimid(ncid,timename,tid)
- if (status.eq.-46) then
- ! no time variable, assume it's a mean value
- nt=1
- else
- ! write error if it was not -46 and get size of time dimension
- call check( nf90_inq_dimid(ncid,timename,tid) )
- call check( nf90_inquire_dimension(ncid,tid,len=nt) )
- endif
-
- call check( nf90_inquire_dimension(ncid,xid,len=nx) )
- call check( nf90_inquire_dimension(ncid,yid,len=ny) )
- call check( nf90_inquire_dimension(ncid,zid,len=nz) )
-
- ! mean so time dimension should be 1
- if( nt.ne.1 ) then
- write(logid,*) 'ERROR get_initfpctm: time dimension should have length of 1 in '//trim(filename)
- stop
- endif
-
- allocate( lon(nx) )
- allocate( lat(ny) )
- allocate( prs(nz) )
- allocate( conc(nx,ny,nz,nt) )
- allocate( work(nx,ny,1) )
- allocate( workout(nxgrid,nygrid,1) )
- allocate( concout(nxgrid,nygrid,nz) )
- allocate( geoh(nz) )
-
- ! read variables
- call check( nf90_inq_varid(ncid,lonname,xid) )
- call check( nf90_get_var(ncid,xid,lon) )
- call check( nf90_inq_varid(ncid,latname,yid) )
- call check( nf90_get_var(ncid,yid,lat) )
- if(preslevel)then
- call check( nf90_inq_varid(ncid,prsname,zid) )
- call check( nf90_get_var(ncid,zid,prs) )
- else
- call check( nf90_inq_varid(ncid,levname,zid) )
- ! save height as geopotential height
- call check( nf90_get_var(ncid,zid,geoh) )
- endif
-
- call check( nf90_inq_varid(ncid,trim(varname),varid) )
- call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
-
- call check( nf90_close(ncid) )
-
- ! calculate geopotential height for each layer
- if ( preslevel ) then
- do kz = 1, nz
- geoh(kz) = (gasc*ts/gc/mmol)*log(psurf/prs(kz))
- enddo
- endif
-
-
- ! interpolate to flexpart grid
- ! ----------------------------
-
- ! interpolate horizontally
- do kz = 1, nz
- work(:,:,1) = conc(:,:,kz,1)
- ! lon_out and lat_out are lon and lat of global flexpart domain
- call convertgrid(midpoints, nx, ny, lon, lat, work, nx_out, ny_out, lon_out, lat_out, workout)
- concout(:,:,kz) = workout(:,:,1)
- end do
-
- ! interpolate vertically
- call interpz(nz, geoh, concout, concini)
-
- deallocate( lon )
- deallocate( lat )
- deallocate( prs )
- deallocate( conc )
- deallocate( work )
- deallocate( workout )
- deallocate( concout )
- deallocate( geoh )
-
-end subroutine get_initfpctm
-
-
diff --git a/source/get_initfpctm_hour.f90 b/source/get_initfpctm_hour.f90
deleted file mode 100644
index de4498e67b477d4f3ba7eb180f434ea255edb1b5..0000000000000000000000000000000000000000
--- a/source/get_initfpctm_hour.f90
+++ /dev/null
@@ -1,142 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: get_initfpctm_hour
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> get_initfpctm_hour
-!!
-!! Purpose: Reads the initial concentration files and interpolates these to
-!! the flexpart grid
-!!
-!! Interface:
-!!
-!! Inputs
-!! filename - name of 4D concentration file
-!! varname - variable name
-!! concini - concentrations on flexpart grid (zero array on input)
-!!
-!! Externals
-!! convertgrid
-!! interpz
-!! check
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine get_initfpctm_hour(filename, varname, concini)
-
- use mod_var
- use mod_save, only : check
- use netcdf
-
- implicit none
-
- character(len=max_path_len), intent(in) :: filename
- character(len=max_name_len), intent(in) :: varname
- real, dimension(nxgrid,nygrid,nzgrid,nztime), intent(in out) :: concini
-
- character(len=*), parameter :: lonname='longitude'
- character(len=*), parameter :: latname='latitude'
- character(len=*), parameter :: altname='lev'
- character(len=*), parameter :: timename='Date'
- integer :: ncid, xid, yid, zid, tid, varid
- integer :: kz, it
- integer :: nx, ny, nz, nt
- integer :: nx_out, ny_out
- logical :: midpoints
- real, dimension(:), allocatable :: lon, lat, alt, time
- real, dimension(nxgrid) :: lon_out
- real, dimension(nygrid) :: lat_out
- real, dimension(:,:,:), allocatable :: work, workout, workout2, concout
- real, dimension(:,:,:,:), allocatable :: conc
-
- ! initialize
- ! ----------
-
- ! don't pass globals to subroutines
- lon_out = gbl_lon
- lat_out = gbl_lat
- nx_out = nxgrid
- ny_out = nygrid
-
- ! read netcdf file
- ! ----------------
-
- write(logid,*) 'Reading file: '//trim(filename)
-
- ! read dimensions
- call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
- call check( nf90_inq_dimid(ncid,lonname,xid) )
- call check( nf90_inq_dimid(ncid,latname,yid) )
- call check( nf90_inq_dimid(ncid,altname,zid) )
- call check( nf90_inq_dimid(ncid,timename,tid) )
- call check( nf90_inquire_dimension(ncid,xid,len=nx) )
- call check( nf90_inquire_dimension(ncid,yid,len=ny) )
- call check( nf90_inquire_dimension(ncid,zid,len=nz) )
- call check( nf90_inquire_dimension(ncid,tid,len=nt) )
-
- allocate( lon(nx) )
- allocate( lat(ny) )
- allocate( alt(nz) )
- allocate( time(nt) )
- allocate( conc(nx,ny,nz,nt) )
- allocate( work(nx,ny,1) )
- allocate( workout(nxgrid,nygrid,1) )
- allocate( workout2(nxgrid,nygrid,nz) )
- allocate( concout(nxgrid,nygrid,nzgrid) )
-
- ! read variables
- call check( nf90_inq_varid(ncid,lonname,xid) )
- call check( nf90_get_var(ncid,xid,lon) )
- call check( nf90_inq_varid(ncid,latname,yid) )
- call check( nf90_get_var(ncid,yid,lat) )
- call check( nf90_inq_varid(ncid,altname,zid) )
- call check( nf90_get_var(ncid,zid,alt) )
- call check( nf90_inq_varid(ncid,timename,tid) )
- call check( nf90_get_var(ncid,tid,time) )
- call check( nf90_inq_varid(ncid,trim(varname),varid) )
- call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
- call check( nf90_close(ncid) )
-
- ! interpolate to flexpart grid
- ! ----------------------------
-
- midpoints = .true.
- do it = 1, nt
- ! interpolate horizontally
- do kz = 1, nz
- work(:,:,1) = conc(:,:,kz,it)
- ! lon_out and lat_out are lon and lat of global flexpart domain
- call convertgrid(midpoints, nx, ny, lon, lat, work, nx_out, ny_out, lon_out, lat_out, workout)
- workout2(:,:,kz) = workout(:,:,1)
- end do
- ! interpolate vertically
- call interpz(nz, alt, workout2, concout)
- concini(:,:,:,it) = concout
- end do
-
- deallocate( lon )
- deallocate( lat )
- deallocate( alt )
- deallocate( conc )
- deallocate( work )
- deallocate( workout )
- deallocate( workout2 )
- deallocate( concout )
-
-end subroutine get_initfpctm_hour
-
-
diff --git a/source/init_cams.f90 b/source/init_cams.f90
deleted file mode 100644
index 4d6f33a952a8363ee79ce4b54136f904c0d08b52..0000000000000000000000000000000000000000
--- a/source/init_cams.f90
+++ /dev/null
@@ -1,265 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: init_cams
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> init_cams
-!!
-!! Purpose: Calculates the initial concentration for each observation using the
-!! 4D initial concentration fields and the flexpart grid_initial files.
-!! Note: adapted from init_cini for CAMS CO2 inversion data.
-!!
-!! Interface:
-!!
-!! Inputs
-!! files - file data structure
-!! config - configuration settings data structure
-!! obs - observation data structure
-!!
-!! Externals
-!! caldate
-!! read_init
-!! get_initconc
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine init_cams(files, config, obs)
-
- use mod_flexpart
- use mod_settings
- use mod_strings
- use mod_var
- use mod_dates
- use mod_obs
- use mod_tools
-
- implicit none
-
- type (files_t), intent (in) :: files
- type (config_t), intent (in) :: config
- type (obs_t), intent (in out) :: obs
-
- character(max_path_len) :: path_flexrec, filename
- character(max_name_len) :: varname
- character(len=8) :: adate, areldate, areltime
- character(len=4) :: ayear
- character(len=2) :: amonth
- logical :: lexist
-
- integer, parameter :: hfreq = 24 ! CAMS N2O is 3-hourly but CH4 is 24-hourly
- integer :: jjjjmm, jjjjmmdd, hhmiss, eomday
- integer :: lasttime, startmonth
- integer :: month, prevmonth
- integer :: ix, jy, kz, i, n, nt, nhours, ntime, ierr
- real(kind=8) :: jdi, hh
- real(kind=8), dimension(nobs) :: jdates
- real, dimension(nobs,ncini) :: cini
- integer, dimension(nobs) :: ind
- real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
- real, dimension(:,:,:,:), allocatable :: concini
-
- real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
- integer :: numpoint, release_nr, nr_init
- integer :: ibdate, ibtime
- real(kind=8), dimension(maxpoint,2) :: releases
-
-
- ! loop over observations
- ! ----------------------
-
- cini(:,:) = 0.
- lasttime = 0
- prevmonth = 0
-
- ! sort obs chronologically for efficiency
- jdates = obs%jdate(:)
- call sort(nobs,jdates,ind)
-
- do i = 1, nobs
-
- ! reset values
- gridini(:,:,:) = 0.
- nr_init = 0
-
- ! check month of observation
- call caldate(jdates(i), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- month = jjjjmmdd/100
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
-
- print*, 'init_cams: path_flexrec = ',path_flexrec
- print*, 'init_cams: average_fp = ',config%average_fp
-
- if ( config%average_fp ) then
-
- ! match and, if needed, average grid_initial files to observation timestamp
-
- ! get header and release information related to current observation location
- ! (if not already read and same for previous observation)
- filename = trim(path_flexrec)//'header'
- if ( i.eq.1 .or. month.ne.prevmonth ) then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- else
- if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- endif
- endif
- ! find the first release number corresponding to the current observation
- release_nr=1
- do while ( int(releases(release_nr,1)*1.e4,kind=8).lt.int(jdates(i)*1.e4,kind=8) .and. release_nr.lt.numpoint )
- release_nr = release_nr+1
- end do
- ! define first filename
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- ! round seconds
- hhmiss = (hhmiss/100)*100
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- ! read first flexpart grid_init files
- inquire ( file=trim(filename), exist=lexist )
- if ( lexist ) then
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- nr_init = 1
- ! current release_nr was already read, start with next
- release_nr = release_nr + 1
- do while ( int(releases(release_nr,2)*1.e4,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1.e4,kind=8) .and. release_nr.lt.numpoint )
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- write(areldate,'(I8.8)') jjjjmmdd
- hhmiss=int(hhmiss/100)*100
- write(areltime,'(I6.6)') hhmiss
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- ! sum the grid with previous initial grids
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- release_nr = release_nr+1
- nr_init = nr_init + 1
- end do ! finished summing grid_initial
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- gridini = gridini/real(nr_init)
- else
- ! filename not found
- write(logid,*) 'WARNING: cannot find '//trim(filename)
- go to 10
- endif
-
- else
-
- ! timestamps of grid_init files match observations
-
- ! define file name
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- print*, 'init_cams: filename = ',filename
- ! read flexpart grid_init file
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- call read_init(filename, gridini)
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
-
- endif
-
- ! calculate termination time of trajectory
- jdi = jdates(i) - trajdays
- call caldate(jdi, jjjjmmdd, hhmiss)
- jjjjmm = jjjjmmdd/100
- ! round time to nearest 3-hours
-! hh = (jdi - floor(jdi))*24d0
-! hh = floor(hh/3d0)*3d0
-
- ! read concentration field only if not already in memory
- if ( jjjjmm.ne.lasttime ) then
- ! allocate concini knowing data freq (from hfreq)
- if (allocated( concini ) ) deallocate( concini )
- eomday = calceomday(jjjjmm)
- ntime = eomday*24/hfreq
- allocate( concini(nxgrid,nygrid,nzgrid,ntime), stat = ierr )
- if ( ierr.ne.0 ) then
- write(logid,*) 'ERROR init_cams: not enough memory'
- stop
- endif
- concini(:,:,:,:) = 0.
- ! read new concentration field
- write(ayear,'(I4.4)') jjjjmm/100
- write(amonth,'(I2.2)') jjjjmm-(jjjjmm/100)*100
- filename = str_replace(files%file_initconc, 'YYYY', ayear)
- filename = str_replace(filename, 'MM', amonth)
- filename = trim(files%path_initconc)//trim(filename)
- print*, 'init_cams: filename = ',filename
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initcams(filename, varname, ntime, concini)
- else
- write(logid,*) 'ERROR init_cams: cannot find file ',filename
- stop
- endif
- endif
-
- ! calculate contribution of initial concentration to receptor
- ! index time dimension of concini
- ! for CAMS N2O 3-hourly data, for CAMS CH4 daily data
- startmonth = (jjjjmmdd/100)*100 + 1
- nhours = int((jdi - juldate(startmonth, 0))*24d0)
- nt = nhours/hfreq + 1
- print*, 'init_cams: startmonth = ',startmonth
- print*, 'init_cams: nhours = ',nhours
- print*, 'init_cams: nt = ',nt
- do kz = 1, nzgrid
- do jy = 1, nygrid
- n = 0
- do while(.true.)
- n = n + 1
- if ( gbl_lat(jy).lt.cini_lat(n) ) exit
- end do
- do ix = 1, nxgrid
- cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*concini(ix,jy,kz,nt)
- end do
- end do
- end do
-
-10 continue
-
- lasttime = jjjjmm
-
- end do
-
- ! reorder cini for observations
- do i = 1, nobs
- obs%cini(ind(i),:) = cini(i,:)
- end do
-
- if (allocated( concini )) deallocate( concini )
-
-end subroutine init_cams
-
-
diff --git a/source/init_co2.f90 b/source/init_co2.f90
index 987645ba932d94d4fa5356d25305b572c2454fff..2c3032f486932aac201f291f2362d219a083fcff 100644
--- a/source/init_co2.f90
+++ b/source/init_co2.f90
@@ -64,7 +64,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
implicit none
- type (files_t), intent (in) :: files
+ type (files_t), intent (in out) :: files
type (config_t), intent (in) :: config
type (fluxes_t), intent (in out) :: fluxes
type (obs_t), intent (in out) :: obs
@@ -79,22 +79,20 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
real, dimension(nbox,ntstate,ndt) :: err_box
real, dimension(nbox,ntstate,nd_nee) :: nee_box
real, dimension(nbox,nbox) :: corr
- real(kind=8), dimension(ntstate) :: statime
- real(kind=8), dimension(:), allocatable :: timeff, timefp
+ real(kind=8), dimension(:), allocatable :: timefp
real, dimension(nxregrid,nyregrid,nd_nee) :: reg_nee
- real(kind=8), dimension(int(statres)*nd_nee) :: timenee
- real, dimension(nxregrid,nyregrid,int(statres)*nd_nee) :: reg_flx
- real, dimension(:,:,:), allocatable :: flx
+ real(kind=8), dimension(:), allocatable :: timenee, timeff, timeocn
+ real, dimension(:,:,:), allocatable :: flx, flxnee, flxff, flxocn
real, dimension(:), allocatable :: xerr
integer, dimension(ndt) :: cnt
integer :: nread, num
- integer :: eomday
+ integer :: eomday, eomday_save
integer :: jjjjmmdd, hhmiss, hh
real :: area
real(kind=8) :: jd, jdi
integer :: ix, jy, ix0, jy0
- integer :: n, nb, nd, nt, i
- integer :: ierr
+ integer :: n, nb, nd, nt, nts, i
+ integer :: ierr, lenstr
real(kind=8), dimension(:), allocatable :: jdate
character(len=3), dimension(:), allocatable :: recs
real :: conc
@@ -106,7 +104,13 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! -------------------
! global fluxes are used for calculating contribution to mixing ratios
- ! from outside the inversion domain i.e. background
+ ! from outside the inversion domain i.e. background and are averaged to
+ ! the temporal resolution of the state vector
+
+ ! note: for global inversions the fossil fuel and NEE fluxes
+ ! averaged to the state vector temporal resolution are not needed
+ ! (no background contribution from fluxes) but are calculated any how to
+ ! avoid additional exceptions in the code
call alloc_fluxes(config, fluxes)
fluxes%flxnee(:,:,:) = 0.
@@ -119,25 +123,24 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
n = 0
jd = juldatei
-! do while ( (jd.le.juldatef).and.(n.lt.ntflx) )
do while ( (jd.le.juldatef).and.(n.lt.ntstate) )
n = n + 1
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
write(amonth,'(I2.2)') (jjjjmmdd-(jjjjmmdd/10000)*10000)/100
-! eomday = calceomday(jjjjmmdd/100)
!--- NEE fluxes
! number of fields to read for each flux time step
-! nread = eomday*nd_nee
- nread = int(statres)*nd_nee
- print*, 'init_co2: n = ',n
+ nread = int(statres)*24/nt_nee
print*, 'init_co2: nee nread = ',nread
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
filename = str_replace(files%filename_nee, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -154,7 +157,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nxgrid, nygrid, &
jd, nread, num, &
flx)
- print*, 'init_co2: num = ',num
+ print*, 'init_co2: nee num = ',num
! convert to kg/m2/s
flx = flx/3600.
! apply numerical scaling
@@ -166,13 +169,15 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
!--- fossil fuel fluxes
! number of fields to read for each flux time step
- ! nread = eomday*24
- nread = int(statres)*24
- print*, 'init_co2: ff nread = ',nread
+ nread = int(statres*24./real(nt_ff))
+ if ( nread.eq.0 ) nread = 1
+ print*, 'init_co2: ff nread = ',nread
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
filename = str_replace(files%filename_ff, 'YYYY', adate)
- filename = str_replace(filename, 'MM', amonth)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -180,8 +185,12 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
stop
endif
write(logid,*) 'Reading fossil fuel :'//trim(filename)
- ! timestamp is hours since 1-Jan each year
- jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
+ if ( nt_ff.lt.24 ) then
+ ! assume timestamp is hours since 1-Jan each year
+ jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
+ else
+ jdi = jd
+ endif
call read_ncdf(filename, &
files%varname_ff, &
files%lonname_ff, &
@@ -191,9 +200,9 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nxgrid, nygrid, &
jdi, nread, num, &
flx)
- print*, 'init_co2: num = ',num
- ! convert from micromol/m2/s to kg/m2/s
- flx = flx*12./1.e9
+ print*, 'init_co2: ff num = ',num
+ ! convert from input flux units to kg/m2/s
+ flx = flx*config%coeff_ff/3600.
! apply numerical scaling
flx = flx*numscale
! average field to flux time step
@@ -203,12 +212,15 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
!--- ocean fluxes
! number of fields to read for each flux time step
-! nread = eomday/int(365/config%nt_flx)
- nread = int(statres)/int(365/config%nt_flx)
+ nread = int(statres*24./real(nt_ocn))
if ( nread.eq.0 ) nread = 1
+ print*, 'init_co2: ocn nread = ',nread
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
filename = str_replace(files%filename_ocn, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -225,6 +237,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nxgrid, nygrid, &
jd, nread, num, &
flx)
+ print*, 'init_co2: ocn num = ',num
! convert to kg/m2/s
flx = flx/3600.
! apply numerical scaling
@@ -236,34 +249,71 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! flux time stamp is start of time step
fluxes%time(n) = jd
-! jd = jd + dble(eomday)
jd = jd + dble(statres)
end do
+ print*, 'init_co2: fluxes%time = ',fluxes%time
! regional prior fluxes
! ---------------------
+ ! these fluxes are used to calculate the contribution from fluxes inside the domain
+ ! fossil fuel and NEE fluxes are averaged/interpolated to temporal resolution of the footprints
+ ! ocean fluxes are average/interpolated to temporal resolution of state vector
+
+ ! note: if nested is false use same fluxes as for global but
+ ! average/interpolate to temporal resolution of the footprints
+
+ if ( .not.config%nested ) then
+ lenstr=len_trim(files%filenest_nee)
+ if ( lenstr.eq.0 ) then
+ files%filenest_nee = files%filename_nee
+ files%varnest_nee = files%varname_nee
+ files%lonnest_nee = files%lonname_nee
+ files%latnest_nee = files%latname_nee
+ files%timenest_nee = files%timename_nee
+ endif
+ lenstr=len_trim(files%filenest_ff)
+ if (lenstr.eq.0) then
+ files%filenest_ff = files%filename_ff
+ files%varnest_ff = files%varname_ff
+ files%lonnest_ff = files%lonname_ff
+ files%latnest_ff = files%latname_ff
+ files%timenest_ff = files%timename_ff
+ endif
+ endif
+
!--- fossil fuel
+ print*, 'nday = ',nday
+ print*, 'nt_ff_reg = ',nt_ff_reg
fluxes%flxff_nest(:,:,:) = 0.
- n = 0
+ print*, 'init_co2: ffreg timesteps input to read = ',(nday*24)/nt_ff_reg
+ allocate( flxff(nxregrid,nyregrid,(nday*24)/nt_ff_reg), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
+ allocate( timeff((nday*24)/nt_ff_reg) )
jd = juldatei
+ n = 0
+ ! loop over months reading all input fluxes
do while ( jd.le.juldatef )
- n = n + 1
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
write(amonth,'(I2.2)') (jjjjmmdd-(jjjjmmdd/10000)*10000)/100
-! eomday = calceomday(jjjjmmdd/100)
- ! number of fields to read for each flux time step
-! nread = eomday*24
- nread = int(statres)*24
+ eomday = calceomday(jjjjmmdd/100)
+ eomday_save = eomday
+ ! no leap years in input data
+ if ( eomday.eq.29 ) eomday = 28
+ print*, 'init_co2: ffreg eomday = ',eomday
+ ! number of fields to read this month
+ nread = eomday*24/nt_ff_reg
print*, 'init_co2: ffreg nread = ',nread
+ if ( nread.eq.0 ) nread = 1
allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
- allocate( timeff(nread) )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
filename = str_replace(files%filenest_ff, 'YYYY', adate)
- filename = str_replace(filename, 'MM', amonth)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -271,8 +321,12 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
stop
endif
write(logid,*) 'Reading fossil fuel :'//trim(filename)
- ! timestamp is hours since 1-Jan each year
- jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
+ if ( nt_ff_reg.lt.24 ) then
+ ! assume timestamp is hours since 1-Jan each year
+ jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
+ else
+ jdi = jd
+ endif
call read_ncdf(filename, &
files%varnest_ff, &
files%lonnest_ff, &
@@ -282,87 +336,90 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nxregrid, nyregrid, &
jdi, nread, num, &
flx)
- print*, 'init_co2: num = ',num
- ! convert from micromol/m2/s to kg/m2/s
- flx = flx*12./1.e9
- ! apply numerical scaling
- flx = flx*numscale
- ! fossil fuel timestamp julian days
-! do i = 1, 24*eomday
- do i = 1, 24*int(statres)
- timeff(i) = jd + dble(i-1)/24d0
+ print*, 'init_co2: ffreg num = ',num
+ ! timestamp julian days
+ do i = 1, nread
+ if ( (n+i).gt.int(real(nday*24)/real(nt_ff_reg)) ) exit
+ timeff(n+i) = jd + dble((i-1)*nt_ff_reg)/24d0
+ flxff(:,:,n+i) = flx(:,:,i)
end do
- print*, 'init_co2: timeff = ',timeff
- ! average/interpolate to footprint timestamp
- if ( ndgrid.eq.24 ) then
- ! no averaging/interpolation needed
-! fluxes%flxff_nest(:,:,(n-1)*ndgrid*eomday+1:n*ndgrid*eomday) = flx
-! fluxes%timefp((n-1)*eomday*ndgrid+1:n*eomday*ndgrid) = timeff
- fluxes%flxff_nest(:,:,(n-1)*ndgrid*int(statres)+1:n*ndgrid*int(statres)) = flx
- fluxes%timefp((n-1)*int(statres)*ndgrid+1:n*int(statres)*ndgrid) = timeff
- else
- ! do averaging/interpolation
-! allocate( datain(nxregrid,eomday*24) )
-! allocate( dataout(nxregrid,eomday*ndgrid) )
-! allocate( timefp(eomday*ndgrid) )
- allocate( datain(nxregrid,int(statres)*24) )
- allocate( dataout(nxregrid,int(statres)*ndgrid) )
- allocate( timefp(int(statres)*ndgrid) )
-! do i = 1, eomday*ndgrid
- do i = 1, int(statres)*ndgrid
- timefp(i) = jd + dble(i-1)/dble(ndgrid)
- enddo
- if ( ndgrid.lt.24 ) then
- ! average flux
- do jy = 1, nyregrid
- datain = flx(:,jy,:)
-! call average(nxregrid, eomday*24, timeff, datain, eomday*ndgrid, timefp, dataout)
-! fluxes%flxff_nest(:,jy,(n-1)*ndgrid*eomday+1:n*ndgrid*eomday) = dataout(:,:)
- call average(nxregrid, int(statres)*24, timeff, datain, int(statres)*ndgrid, timefp, dataout)
- fluxes%flxff_nest(:,jy,(n-1)*ndgrid*int(statres)+1:n*ndgrid*int(statres)) = dataout(:,:)
- end do
- else if ( ndgrid.gt.24 ) then
- ! interpolate flux
- do jy = 1, nyregrid
- datain = flx(:,jy,:)
-! call interp(nxregrid, eomday*24, timeff, datain, eomday*ndgrid, timefp, dataout)
-! fluxes%flxff_nest(:,jy,(n-1)*ndgrid*eomday+1:n*ndgrid*eomday) = dataout(:,:)
- call interp(nxregrid, int(statres)*24, timeff, datain, int(statres)*ndgrid, timefp, dataout)
- fluxes%flxff_nest(:,jy,(n-1)*ndgrid*int(statres)+1:n*ndgrid*int(statres)) = dataout(:,:)
- end do
- endif
-! fluxes%timefp((n-1)*eomday*ndgrid+1:n*eomday*ndgrid) = timefp
- fluxes%timefp((n-1)*int(statres)*ndgrid+1:n*int(statres)*ndgrid) = timefp
- deallocate( datain )
- deallocate( dataout )
- deallocate( timefp )
- endif
deallocate( flx )
- deallocate( timeff )
-! jd = jd + dble(eomday)
- jd = jd + dble(statres)
+ n = n + nread
+ jd = jd + dble(eomday_save)
end do
+ n = min(n, int(real(nday*24)/real(nt_ff_reg)))
+ print*, 'init_co2: ffreg n = ',n
+ print*, 'init_co2: timeff = ',timeff
+
+ ! convert from input flux units to kg/m2/s
+ flxff = flxff*config%coeff_ff_reg/3600.
+ ! apply numerical scaling
+ flxff = flxff*numscale
+
+ ! average/interpolate to footprint timestamp
+ allocate( datain(nxregrid,n) )
+ allocate( dataout(nxregrid,nday*ndgrid) )
+ allocate( timefp(nday*ndgrid) )
+ do i = 1, nday*ndgrid
+ timefp(i) = juldatei + dble(i-1)/dble(ndgrid)
+ enddo
+ print*, 'init_co2: timefp = ',timefp
+ if ( ndgrid.lt.(24/nt_ff_reg) ) then
+ ! average flux
+ do jy = 1, nyregrid
+ datain = flxff(:,jy,:)
+ call average(nxregrid, n, timeff, datain, nday*ndgrid, timefp, dataout)
+ fluxes%flxff_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( ndgrid.gt.(24/nt_ff_reg) ) then
+ ! interpolate flux
+ do jy = 1, nyregrid
+ datain = flxff(:,jy,:)
+ call interp(nxregrid, n, timeff, datain, nday*ndgrid, timefp, dataout)
+ fluxes%flxff_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( ndgrid.eq.(24/nt_ff_reg) ) then
+ fluxes%flxff_nest(:,:,:) = flxff
+ endif
+ fluxes%timefp(:) = timefp
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( timefp )
+ deallocate( flxff )
+ deallocate( timeff )
!--- ocean fluxes
+ ! ocean contribution from ocean fluxes inside the inversion domain
+ ! is calculated at same resolution as state vector and the background
+ ! contribution, hence, only need to do following if nested = true
+
if ( config%nested ) then
+ write(logid,*) 'Reading regional ocean fluxes'
! nested
fluxes%flxocn_nest(:,:,:) = 0.
+ print*, 'init_co2: ocn reg timesteps input to read = ',(nday*24)/nt_ocn
+ allocate( flxocn(nxregrid,nyregrid,(nday*24)/nt_ocn), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
+ allocate( timeocn((nday*24)/nt_ocn) )
n = 0
jd = juldatei
-! do while ( (jd.le.juldatef).and.(n.lt.ntflx) )
- do while ( (jd.le.juldatef).and.(n.lt.ntstate) )
- n = n + 1
+ do while ( jd.le.juldatef )
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
-! eomday = calceomday(jjjjmmdd/100)
- ! number of fields to read for each flux time step
-! nread = eomday/int(365/config%nt_flx)
- nread = int(statres)/int(365/config%nt_flx)
+ eomday = calceomday(jjjjmmdd/100)
+ eomday_save = eomday
+ ! no leap years in input data
+ if ( eomday.eq.29 ) eomday = 28
+ ! number of fields to read this month
+ nread = eomday*24/nt_ocn
if ( nread.eq.0 ) nread = 1
allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
filename = str_replace(files%filenest_ocn, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -379,41 +436,83 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nxregrid, nyregrid, &
jd, nread, num, &
flx)
- ! convert to kg/m2/s
- flx = flx/3600.
- ! apply numerical scaling
- flx = flx*numscale
- ! average field to flux time step
- fluxes%flxocn_nest(:,:,n) = sum(flx,dim=3)/real(nread)
+ ! timestamp in julian days
+ do i = 1, nread
+ if ( (n+i).gt.(nday*24)/nt_ocn ) exit
+ timeocn(n+i) = jd + dble((i-1)*nt_ocn)/24d0
+ flxocn(:,:,n+i) = flx(:,:,i)
+ end do
deallocate( flx )
-! jd = jd + dble(eomday)
- jd = jd + dble(statres)
+ n = n + nread
+ jd = jd + dble(eomday_save)
end do
+ n = min(n, (nday*24)/nt_ocn)
+ print*, 'init_co2: ocnreg n = ',n
+ print*, 'init_co2: timeocn = ',timeocn
+
+ ! convert to kg/m2/s
+ flxocn = flxocn/3600.
+ ! apply numerical scaling
+ flxocn = flxocn*numscale
+
+ ! average/interpolate to state vector timestamp
+ allocate( datain(nxregrid,n) )
+ allocate( dataout(nxregrid,ntstate) )
+ if ( statres.gt.real(nt_ocn/24) ) then
+ ! average flux
+ do jy = 1, nyregrid
+ datain = flxocn(:,jy,:)
+ call average(nxregrid, n, timeocn, datain, ntstate, fluxes%time, dataout)
+ fluxes%flxocn_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.lt.real(nt_ocn/24) ) then
+ ! interpolate flux
+ do jy = 1, nyregrid
+ datain = flxocn(:,jy,:)
+ call interp(nxregrid, n, timeocn, datain, ntstate, fluxes%time, dataout)
+ fluxes%flxocn_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.eq.real(nt_ocn/24) ) then
+ fluxes%flxocn_nest(:,:,:) = flxocn(:,:,:)
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( flxocn )
+ deallocate( timeocn )
+
endif
!--- NEE fluxes
- ! read hourly NEE and calculate uncertainty for boxes
+ ! read regional NEE and average/interpolate to footprint timestamp
write(logid,*) 'Reading hourly NEE'
+
fluxes%flxnee_nest(:,:,:) = 0.
- err_box(:,:,:) = 0.
- nee_box(:,:,:) = 0.
+ allocate( flxnee(nxregrid,nyregrid,nday*nd_nee), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
+ allocate( timenee(nday*nd_nee) )
+ n = 0
jd = juldatei
- nd = 1
-
- do while ( (jd.le.(juldatef-dble(statres)+1d0)).and.(nd.le.ntstate) )
-
- print*, 'init_co2: hourly nee nd = ',nd
- reg_flx(:,:,:) = 0.
+ do while ( jd.le.juldatef )
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,'(I4.4)') jjjjmmdd/10000
-
- ! read NEE fluxes for current day (starting from 00:00)
+ eomday = calceomday(jjjjmmdd/100)
+ eomday_save = eomday
+ print*, 'init_co2: eomday = ',eomday
+ ! no leap years in input data
+ if ( eomday.eq.29 ) eomday = 28
+ ! number of fields to read this month
+ nread = eomday*nd_nee
+ if ( nread.eq.0 ) nread = 1
+ print*, 'init_co2: neereg nread = ',nread
+ allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
filename = str_replace(files%filenest_nee, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
- print*, filename
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
write(logid,*) 'ERROR: cannot find '//trim(filename)
@@ -426,81 +525,71 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
files%timenest_nee,&
rllx, rlly, &
nxregrid, nyregrid, &
- jd, int(statres)*nd_nee, num, &
- reg_flx)
- ! check if everything read
- print*, 'init_co2: num = ',num
- if ( num.ne.int(statres)*nd_nee ) then
- jdi = jd + num/nd_nee
- nread = int(statres)*nd_nee - num
- allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
- if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
- ! read remaining days from next file
- call caldate(jdi, jjjjmmdd, hhmiss)
- write(adate,'(I4.4)') jjjjmmdd/10000
- filename = str_replace(files%filenest_nee, 'YYYY', adate)
- filename = trim(files%path_prior)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if (.not.lexist) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
- stop
- endif
- call read_ncdf(filename, &
- files%varnest_nee, &
- files%lonnest_nee, &
- files%latnest_nee, &
- files%timenest_nee,&
- rllx, rlly, &
- nxregrid, nyregrid, &
- jdi, nread, num, &
- flx)
- reg_flx(:,:,int(statres)*nd_nee-nread+1:int(statres)*nd_nee) = flx(:,:,:)
- deallocate( flx )
- endif
- ! convert to kg/m2/s
- reg_flx = reg_flx/3600.
- ! apply numerical scaling
- reg_flx = reg_flx*numscale
+ jd, nread, num, &
+ flx)
+ do i = 1, nread
+ if ( (n+i).gt.(nday*nd_nee) ) exit
+ timenee(n+i) = jd + dble((i-1)*nt_nee)/24d0
+ flxnee(:,:,n+i) = flx(:,:,i)
+ end do
+ deallocate( flx )
+ n = n + nread
+ print*, 'init_co2: n = ',n
+ jd = jd + dble(eomday_save)
+ end do
+ n = min(n, (nday*nd_nee))
+ print*, 'init_co2: neereg n = ',n
+ print*, 'init_co2: timenee = ',timenee
+ print*, 'init_co2: ndgrid, nd_nee = ',ndgrid, nd_nee
+
+ ! convert to kg/m2/s
+ flxnee = flxnee/3600.
+ ! apply numerical scaling
+ flxnee = flxnee*numscale
- do i = 1, int(statres)*nd_nee
- timenee(i) = jd + dble(i-1)/dble(nd_nee)
+ ! average/interpolate to footprint time
+ allocate( datain(nxregrid,n) )
+ allocate( dataout(nxregrid,nday*ndgrid) )
+ if ( ndgrid.lt.nd_nee ) then
+ ! average flux
+ print*, 'init_co2: averaging nee to footprint time resolution'
+ do jy = 1, nyregrid
+ datain = flxnee(:,jy,:)
+ call average(nxregrid, n, timenee, datain, nday*ndgrid, fluxes%timefp, dataout)
+ fluxes%flxnee_nest(:,jy,:) = dataout(:,:)
end do
- print*, 'init_co2: timenee = ',timenee
+ else if ( ndgrid.gt.nd_nee ) then
+ ! interpolate flux
+ print*, 'init_co2: interpolating nee to footprint time resolution'
+ do jy = 1, nyregrid
+ datain = flxnee(:,jy,:)
+ call interp(nxregrid, n, timenee, datain, nday*ndgrid, fluxes%timefp, dataout)
+ fluxes%flxnee_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( ndgrid.eq.nd_nee ) then
+ print*, 'init_co2: fluxes have already footprint time resolution'
+ fluxes%flxnee_nest(:,:,:) = flxnee(:,:,:)
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( timenee )
- ! average/interpolate to footprint time
- if ( ndgrid.eq.nd_nee ) then
- ! no averaging/interpolation needed
- fluxes%flxnee_nest(:,:,(nd-1)*int(statres)*ndgrid+1:nd*int(statres)*ndgrid) = reg_flx(:,:,:)
- else
- ! do averaging/interpolation
- allocate( datain(nxregrid,int(statres)*nd_nee) )
- allocate( dataout(nxregrid,int(statres)*ndgrid) )
- allocate( timefp(int(statres)*ndgrid) )
- timefp(:) = fluxes%timefp((nd-1)*int(statres)*ndgrid+1:nd*int(statres)*ndgrid)
- if ( ndgrid.lt.nd_nee ) then
- ! average flux
- do jy = 1, nyregrid
- datain = reg_flx(:,jy,:)
- call average(nxregrid, int(statres)*nd_nee, timenee, datain, int(statres)*ndgrid, timefp, dataout)
- fluxes%flxnee_nest(:,jy,(nd-1)*int(statres)*ndgrid+1:nd*int(statres)*ndgrid) = dataout(:,:)
- end do
- else if ( ndgrid.gt.nd_nee ) then
- ! interpolate flux
- do jy = 1, nyregrid
- datain = reg_flx(:,jy,:)
- call interp(nxregrid, int(statres)*nd_nee, timenee, datain, int(statres)*ndgrid, timefp, dataout)
- fluxes%flxnee_nest(:,jy,(nd-1)*int(statres)*ndgrid+1:nd*int(statres)*ndgrid) = dataout(:,:)
- end do
- endif
- deallocate( datain )
- deallocate( dataout )
- deallocate( timefp )
- endif
+ ! uncertainty for NEE
+ ! -------------------
+
+ ! calculate state vector uncertainty for NEE
+ ! on aggregated grid and for statres and statres_hr timestamps
+
+ err_box(:,:,:) = 0.
+ nee_box(:,:,:) = 0.
+
+ do nts = 1, ntstate
! average hourly NEE over statres dimension
reg_nee(:,:,:) = 0.
do i = 1, int(statres)
- reg_nee(:,:,:) = reg_nee(:,:,:) + reg_flx(:,:,(i-1)*nd_nee+1:i*nd_nee)
+ n = (nts-1)*int(statres) + i
+ reg_nee(:,:,:) = reg_nee(:,:,:) + flxnee(:,:,(n-1)*nd_nee+1:n*nd_nee)
end do
reg_nee(:,:,:) = reg_nee(:,:,:)/statres
@@ -510,11 +599,11 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).gt.0 ) then
- nee_box(nbox_xy(ix,jy),nd,n) = nee_box(nbox_xy(ix,jy),nd,n) + reg_nee(ix,jy,n) * area
+ nee_box(nbox_xy(ix,jy),nts,n) = nee_box(nbox_xy(ix,jy),nts,n) + reg_nee(ix,jy,n) * area
endif
end do
end do
- nee_box(:,nd,n) = nee_box(:,nd,n)/area_box
+ nee_box(:,nts,n) = nee_box(:,nts,n)/area_box
end do
! calculate prior NEE uncertainty for each box
@@ -530,19 +619,15 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
do nt = 1, ndt
if ( (hh.ge.(nt-1)*24/ndt).and.(hh.lt.nt*24/ndt) ) then
cnt(nt) = cnt(nt) + 1
- err_box(nb,nd,nt) = err_box(nb,nd,nt) + nee_box(nb,nd,n)*config%flxerr
+ err_box(nb,nts,nt) = err_box(nb,nts,nt) + nee_box(nb,nts,n)*config%flxerr
endif
end do
end do
do nt = 1, ndt
- err_box(nb,nd,nt) = err_box(nb,nd,nt)/cnt(nt)
+ err_box(nb,nts,nt) = err_box(nb,nts,nt)/cnt(nt)
end do
end do
- statime(nd) = jd
- nd = nd + 1
- jd = jd + real(statres,kind=8)
-
end do
! verbose output for testing
@@ -559,13 +644,13 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
end do
end do
close(100)
- ! NEE for region before interpolation (last statres timestep only)
- filename = trim(files%path_output)//'reg_flx.txt'
+ ! NEE for region before interpolation (first statres timestep only)
+ filename = trim(files%path_output)//'flxnee.txt'
write(rowfmt,'(A,I6,A)') '(',nxregrid,'(E11.4,1X))'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
do nd = 1, int(statres)*nd_nee
do n = 1, nyregrid
- write(100,fmt=rowfmt) reg_flx(:,n,nd)
+ write(100,fmt=rowfmt) flxnee(:,n,nd)
end do
end do
close(100)
@@ -573,22 +658,24 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
filename = trim(files%path_output)//'nee_box.txt'
write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
- do nd = 1, ntstate
+ do nts = 1, ntstate
do n = 1, nd_nee
- write(100,fmt=rowfmt) nee_box(:,nd,n)
+ write(100,fmt=rowfmt) nee_box(:,nts,n)
end do
end do
close(100)
! NEE uncertainty for boxes
write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
open(100,file=trim(files%path_output)//'err_box.txt',status='replace',action='write',iostat=ierr)
- do nd = 1, ntstate
+ do nts = 1, ntstate
do n = 1, ndt
- write(100,fmt=rowfmt) err_box(:,nd,n)
+ write(100,fmt=rowfmt) err_box(:,nts,n)
end do
end do
close(100)
endif
+
+ deallocate( flxnee )
! prior state vector
! ------------------
@@ -614,7 +701,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
if ( config%opt_cini ) states%px0(npvar+1:nvar) = 1.
! assign time stamps to flux time steps
- states%xtime(:) = statime
+ states%xtime(:) = fluxes%time(:)
print*, 'init_co2: states%xtime = ',states%xtime
! assign time stamps to mixing ratio scalar time steps
diff --git a/source/init_fpctm.f90 b/source/init_fpctm.f90
deleted file mode 100644
index 3d5e376e07d0e1ae6ee5fe7091e39915f2b27682..0000000000000000000000000000000000000000
--- a/source/init_fpctm.f90
+++ /dev/null
@@ -1,320 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: init_fpctm
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> init_fpctm
-!!
-!! Purpose: Calculates the initial concentration for each observation using the
-!! 4D initial concentration fields and the flexpart grid_initial files.
-!!
-!! Interface:
-!!
-!! Inputs
-!! files - file data structure
-!! config - configuration settings data structure
-!! obs - observation data structure
-!!
-!! Externals
-!! caldate
-!! read_init
-!! get_initfpctm
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine init_fpctm(files, config, obs)
-
- use mod_flexpart
- use mod_settings
- use mod_strings
- use mod_var
- use mod_dates
- use mod_obs
- use mod_tools
-
- implicit none
-
- type (files_t), intent (in) :: files
- type (config_t), intent (in) :: config
- type (obs_t), intent (in out) :: obs
-
- character(max_path_len) :: path_flexrec, pathname, filename
- character(max_name_len) :: varname
- character(len=8) :: adate, ayear, areldate, areltime
- logical :: lexist
- integer :: lastjjjjmm, prejjjjmm, projjjjmm
- integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
- integer :: ix, jy, kz, i, n
- integer :: month, prevmonth
- real :: conc
- real(kind=8) :: jdmid, jdpre, jdlast
- real(kind=8), dimension(nobs) :: jdates
- real, dimension(nobs,ncini) :: cini
- integer, dimension(nobs) :: ind
- real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
- real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
-
- real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
- integer :: numpoint, release_nr, nr_init
- integer :: ibdate, ibtime
- real(kind=8), dimension(maxpoint,2) :: releases
-
-
- ! loop over observations
- ! ----------------------
-
- lastjjjjmm = 0
- concini(:,:,:) = 0.
- preconcini(:,:,:) = 0.
- proconcini(:,:,:) = 0.
- cini(:,:) = 0.
- prevmonth = 0
-
- ! sort obs chronologically for efficiency
- jdates = obs%jdate(:)
- call sort(nobs,jdates,ind)
-
- do i = 1, nobs
-
- ! reset values
- gridini(:,:,:) = 0.
- nr_init = 0
-
- ! check month of observation
- call caldate(jdates(i), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- month = jjjjmmdd/100
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
- print*, 'init_fpctm: pathflexrec = ',path_flexrec
-
- if ( config%average_fp ) then
-
- ! match and, if needed, average grid_initial files to observation timestamp
-
- ! get header and release information related to current observation location
- ! (if not already read and same for previous observation)
- filename = trim(path_flexrec)//'header'
- if ( i.eq.1 .or. month.ne.prevmonth ) then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- else
- if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- endif
- endif
- ! find the first release number corresponding to the current observation
- release_nr=1
- do while ( int(releases(release_nr,1)*1e4,kind=8).lt.int(jdates(i)*1e4,kind=8) .and. release_nr.lt.numpoint )
- release_nr = release_nr+1
- end do
- ! define first filename
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- ! round seconds
- hhmiss = (hhmiss/100)*100
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- print*, 'init_fpctm: filename = ',filename
- ! read first flexpart grid_init files
- inquire ( file=trim(filename), exist=lexist )
- if ( lexist ) then
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- nr_init = 1
- ! current release_nr was already read, start with next
- release_nr = release_nr + 1
- do while ( int(releases(release_nr,2)*1e4,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1e4,kind=8) .and. release_nr.lt.numpoint )
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- write(areldate,'(I8.8)') jjjjmmdd
- hhmiss=int(hhmiss/100)*100
- write(areltime,'(I6.6)') hhmiss
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- ! sum the grid with previous initial grids
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- release_nr = release_nr+1
- nr_init = nr_init + 1
- end do ! finished summing grid_initial
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- gridini = gridini/real(nr_init)
- else
- ! filename not found
- write(logid,*) 'WARNING: cannot find '//trim(filename)
- go to 10
- endif
-
- else
-
- ! timestamps of grid_init files match observations
-
- ! define file name
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- print*, 'init_fpctm: filename = ',filename
- ! read flexpart grid_init file
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- call read_init(filename, gridini)
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
-
- endif
-
- ! calculate termination time of trajectory
- call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
- jjjjmm = jjjjmmdd/100
-
- ! preceding and proceding months
- mm = jjjjmm - (jjjjmm/100)*100
- if( mm.gt.1 ) then
- prejjjjmm = jjjjmm - 1
- else
- prejjjjmm = (jjjjmm/100 - 1)*100 + 12
- endif
- if( mm.lt.12 ) then
- projjjjmm = jjjjmm + 1
- else
- projjjjmm = (jjjjmm/100 + 1)*100 + 1
- endif
-
- jdmid = juldate(jjjjmm*100+14, 0)
- jdpre = juldate(prejjjjmm*100+14, 0)
- ! prevent 0 date in first time step (need to check this)
- if( lastjjjjmm.ne.0 ) then
- jdlast = juldate(lastjjjjmm*100+14, 0)
- else
- jdlast = 0
- endif
- eomday = calceomday(prejjjjmm)
-
- ! get new concentrations if not same month as previous observation
- if( jjjjmm.ne.lastjjjjmm ) then
- if( (lastjjjjmm.ne.0).and.(jdmid.eq.(jdlast+dble(eomday))) ) then
- ! this is proceeding month
- preconcini = concini
- concini = proconcini
- ! get proceeding concentration field
- write(adate,'(I6.6)') projjjjmm
- write(ayear,'(I4.4)') projjjjmm/100
- pathname = str_replace(files%path_initconc, 'YYYY', ayear)
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(pathname)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initfpctm(filename, varname, proconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
- proconcini = concini
- endif
- else
- ! get current concentration field
- write(adate,'(I6.6)') jjjjmm
- write(ayear,'(I4.4)') jjjjmm/100
- pathname = str_replace(files%path_initconc, 'YYYY', ayear)
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(pathname)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initfpctm(filename, varname, concini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
- ! get previous concentration field
- write(adate,'(I6.6)') prejjjjmm
- write(ayear,'(I4.4)') prejjjjmm/100
- pathname = str_replace(files%path_initconc, 'YYYY', ayear)
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(pathname)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initfpctm(filename, varname, preconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
- ! get next concentration field
- write(adate,'(I6.6)') projjjjmm
- write(ayear,'(I4.4)') projjjjmm/100
- pathname = str_replace(files%path_initconc, 'YYYY', ayear)
- filename = str_replace(files%file_initconc, 'YYYYMM', adate)
- filename = trim(pathname)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- varname = files%varname_init
- call get_initfpctm(filename, varname, proconcini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
- proconcini=concini
- endif
- endif
- endif
-
- ! calculate contribution of initial concentration to receptor
- do kz = 1, nzgrid
- do jy = 1, nygrid
- n = 0
- do while(.true.)
- n = n + 1
- if ( gbl_lat(jy).lt.cini_lat(n) ) exit
- end do
- do ix = 1, nxgrid
- ! interpolate between concentrations
- if( (jdates(i) - trajdays).lt.jdmid ) then
- conc = preconcini(ix,jy,kz) + (concini(ix,jy,kz) - preconcini(ix,jy,kz)) * &
- real((jdates(i) - trajdays) - jdpre, kind=4)/30.
- else
- conc = concini(ix,jy,kz) + (proconcini(ix,jy,kz) - concini(ix,jy,kz)) * &
- real((jdates(i) - trajdays) - jdmid, kind=4)/30.
- endif
- cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*conc
- end do
- end do
- end do
-
-10 continue
-
- lastjjjjmm = jjjjmm
- prevmonth = month
-
- end do
-
- ! reorder cini for observations
- do i = 1, nobs
- obs%cini(ind(i),:) = cini(i,:)
- end do
-
-end subroutine init_fpctm
-
-
diff --git a/source/init_fpctm_hour.f90 b/source/init_fpctm_hour.f90
deleted file mode 100644
index 108cefc0a8c1ddd4778293e7de796b7b1f666b6a..0000000000000000000000000000000000000000
--- a/source/init_fpctm_hour.f90
+++ /dev/null
@@ -1,259 +0,0 @@
-!---------------------------------------------------------------------------------------
-! FLEXINVERT: init_fpctm_hour
-!---------------------------------------------------------------------------------------
-! FLEXINVERT is free software: you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! FLEXINVERT is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License
-! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
-!
-! Copyright 2017, Rona Thompson
-!---------------------------------------------------------------------------------------
-!
-!> init_fpctm_hour
-!!
-!! Purpose: Calculates the initial concentration for each observation using the
-!! 4D initial concentration fields and the flexpart grid_initial files.
-!! Note: adapted from init_cini for FP-CTM hourly output
-!!
-!! Interface:
-!!
-!! Inputs
-!! files - file data structure
-!! config - configuration settings data structure
-!! obs - observation data structure
-!!
-!! Externals
-!! caldate
-!! read_init
-!! get_initfpctm_hour
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine init_fpctm_hour(files, config, obs)
-
- use mod_flexpart
- use mod_settings
- use mod_strings
- use mod_var
- use mod_dates
- use mod_obs
- use mod_tools
- use netcdf
- use mod_save, only : check
-
- implicit none
-
- type (files_t), intent (in) :: files
- type (config_t), intent (in) :: config
- type (obs_t), intent (in out) :: obs
-
- character(max_path_len) :: path_flexrec, pathname, filename
- character(max_name_len) :: varname
- character(len=8) :: adate, ayear, areldate, areltime, amonth, atime
- logical :: lexist
- integer :: jjjjmm, jjjjmmdd, hhmiss, lasttime, currenttime
- integer :: ix, jy, kz, i, n, ncid, tid, it
- integer :: month, prevmonth
- real(kind=8) :: jdi, hh
- real(kind=8), dimension(nobs) :: jdates
- real(kind=8), dimension(:), allocatable:: time
- real, dimension(nobs,ncini) :: cini
- integer, dimension(nobs) :: ind
- real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
- real, dimension(:,:,:,:), allocatable :: concini
-
- real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
- integer :: numpoint, release_nr, nr_init
- integer :: ibdate, ibtime
- real(kind=8), dimension(maxpoint,2) :: releases
-
-
- ! loop over observations
- ! ----------------------
-
- cini(:,:) = 0.
- lasttime = 0
- prevmonth = 0
-
- ! sort obs chronologically for efficiency
- jdates = obs%jdate(:)
- call sort(nobs,jdates,ind)
-
- do i = 1, nobs
-
- ! reset values
- gridini(:,:,:) = 0.
- nr_init = 0
-
- ! check month of observation
- call caldate(jdates(i), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- month = jjjjmmdd/100
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
- print*, 'init_fpctm: pathflexrec = ',path_flexrec
-
- if ( config%average_fp ) then
-
- ! match and, if needed, average grid_initial files to observation timestamp
-
- ! get header and release information related to current observation location
- ! (if not already read and same for previous observation)
- filename = trim(path_flexrec)//'header'
- if ( i.eq.1 .or. month.ne.prevmonth ) then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- else
- if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
- write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
- call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
- endif
- endif
- ! find the first release number corresponding to the current observation
- release_nr=1
- do while ( int(releases(release_nr,1)*1e4,kind=8).lt.int(jdates(i)*1e4,kind=8) .and. release_nr.lt.numpoint )
- release_nr = release_nr+1
- end do
- ! define first filename
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- ! round seconds
- hhmiss = (hhmiss/100)*100
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- print*, 'init_fpctm: filename = ',filename
- ! read first flexpart grid_init files
- inquire ( file=trim(filename), exist=lexist )
- if ( lexist ) then
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- nr_init = 1
- ! current release_nr was already read, start with next
- release_nr = release_nr + 1
- do while ( int(releases(release_nr,2)*1e4,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1e4,kind=8) .and. release_nr.lt.numpoint )
- call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
- write(adate,'(I6.6)') jjjjmmdd/100
- write(areldate,'(I8.8)') jjjjmmdd
- hhmiss=int(hhmiss/100)*100
- write(areltime,'(I6.6)') hhmiss
- path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- ! sum the grid with previous initial grids
- call read_init(filename, gridini_tmp)
- gridini = gridini + gridini_tmp
- release_nr = release_nr+1
- nr_init = nr_init + 1
- end do ! finished summing grid_initial
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- gridini = gridini/real(nr_init)
- else
- ! filename not found
- write(logid,*) 'WARNING: cannot find '//trim(filename)
- go to 10
- endif
-
- else
-
- ! timestamps of grid_init files match observations
-
- ! define file name
- write(areldate,'(I8.8)') jjjjmmdd
- write(areltime,'(I6.6)') hhmiss
- filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
- print*, 'init_fpctm: filename = ',filename
- ! read flexpart grid_init file
- inquire(file=trim(filename),exist=lexist)
- if(lexist) then
- call read_init(filename, gridini)
- ! convert from ppt to fraction of one
- gridini = gridini*1.e-12
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
-
- endif
-
- ! calculate termination time of trajectory
- jdi = jdates(i) - trajdays
- call caldate(jdi, jjjjmmdd, hhmiss)
- jjjjmm = jjjjmmdd/100
- ! round time to nearest 3-hours
- hh = (jdi - floor(jdi))*24d0
- hh = floor(hh/3d0)*3d0
-
- ! read concentration field only if not already in memory
- currenttime = jjjjmm*100
- if ( currenttime.ne.lasttime ) then
- ! read new concentration field
- write(adate,'(I8.8)') jjjjmmdd
- write(ayear,'(I4.4)') jjjjmm/100
- write(amonth,'(I6.6)') jjjjmm
- write(atime,'(I2.2)') int(hh)
- pathname = str_replace(files%path_initconc,'YYYY', trim(ayear))
- filename = str_replace(files%file_initconc, 'YYYYMM', trim(amonth))
- filename = trim(pathname)//trim(filename)
- inquire(file=trim(filename),exist=lexist)
- print*, filename
- if(lexist) then
- ! read time dimension from file (nztime is global)
- call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
- call check( nf90_inq_dimid(ncid,'Date',tid) )
- call check( nf90_inquire_dimension(ncid,tid,len=nztime) )
- if( allocated(time) ) deallocate(time)
- allocate( time(nztime) )
- call check( nf90_inq_varid(ncid,'Date',tid) )
- call check( nf90_get_var(ncid,tid,time) )
- call check( nf90_close(ncid) )
- ! read conc and interpolate to flexpart grid
- if( allocated(concini) ) deallocate(concini)
- allocate( concini(nxgrid,nygrid,nzgrid,nztime) )
- varname = files%varname_init
- call get_initfpctm_hour(filename, varname, concini)
- else
- write(logid,*) 'WARNING: cannot find ',filename
- go to 10
- endif
- endif
-
- ! calculate contribution of initial concentration to receptor
-! it = minloc(abs(time-jdates(i)),dim=1,mask=abs(time-jdates(i)).eq.minval(abs(time-jdates(i))))
- it = minloc(abs(time-jdi),dim=1,mask=abs(time-jdi).eq.minval(abs(time-jdi)))
- do kz = 1, nzgrid
- do jy = 1, nygrid
- n = 0
- do while(.true.)
- n = n + 1
- if ( gbl_lat(jy).lt.cini_lat(n) ) exit
- end do
- do ix = 1, nxgrid
- cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*concini(ix,jy,kz,it)
- end do
- end do
- end do
-
-10 continue
-
- lasttime = currenttime
-
- end do
-
- ! reorder cini for observations
- do i = 1, nobs
- obs%cini(ind(i),:) = cini(i,:)
- end do
-
-end subroutine init_fpctm_hour
-
-
diff --git a/source/init_ghg.f90 b/source/init_ghg.f90
index c652bbb3470a00ac9e7fcad7a353a0065542dac8..3ad231446b01b28839ef7b7682c1daba10d7a8fc 100644
--- a/source/init_ghg.f90
+++ b/source/init_ghg.f90
@@ -74,8 +74,10 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
character(max_path_len) :: filename
character(max_path_len) :: rowfmt
character(len=4) :: adate
+ character(len=2) :: amonth
logical :: lexist
- real, dimension(:,:,:), allocatable :: flx
+ real, dimension(:,:,:), allocatable :: flx, flxghg
+ real, dimension(:,:), allocatable :: datain, dataout
real, dimension(nbox,ntstate) :: flx_box
real, dimension(nbox,ntstate) :: err_box
real, dimension(nbox) :: xerr
@@ -84,9 +86,9 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
real :: area
real(kind=8) :: jd
integer :: ix, jy, ix0, jy0, ix1, ix2, jy1, jy2
- integer :: n, nb, nread, num
+ integer :: i, n, nb, nread, num, eomday
integer :: ierr
- real(kind=8), dimension(:), allocatable :: jdate
+ real(kind=8), dimension(:), allocatable :: jdate, timeghg
character(len=3), dimension(:), allocatable:: recs
real :: conc
character(len=20) :: dump
@@ -95,8 +97,8 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! global prior fluxes
! -------------------
- ! If nested run global fluxes are only used to calculate the contribution to mixing ratios
- ! from fluxes outside the inversion domain. Otherwise, global fluxes are used to
+ ! if nested run global fluxes are only used to calculate the contribution to mixing ratios
+ ! from fluxes outside the inversion domain, otherwise, global fluxes are used to
! calculate the contribution to mixing ratios and as prior fluxes inside the domain.
call alloc_fluxes(config, fluxes)
@@ -105,20 +107,26 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! loop over prior flux files
+ allocate( flxghg(nxgrid,nygrid,(nday*24/nt_flx)), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_ghg: not enough memory'
+ allocate( timeghg(nday*24/nt_flx) )
+
n = 0
jd = juldatei
-
- do while ( (jd.le.juldatef).and.(n.lt.ntstate) )
-
- n = n + 1
+ do while ( jd.le.juldatef )
call caldate(jd, yyyymmdd, hhmiss)
write(adate,'(I4.4)') yyyymmdd/10000
-
- nread = nint(statres/(real(365)/real(config%nt_flx)))
+ write(amonth,'(I2.2)') (yyyymmdd-(yyyymmdd/10000)*10000)/100
+ eomday = calceomday(yyyymmdd/100)
+ nread = eomday*24/nt_flx
if ( nread.eq.0 ) nread = 1
+ print*, 'init_ghg: nread = ',nread
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_ghg: not enough memory'
filename = str_replace(files%filename_flx, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
write(logid,*) 'Reading prior fluxes :'//trim(filename)
inquire(file=trim(filename),exist=lexist)
@@ -135,51 +143,90 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
nxgrid, nygrid, &
jd, nread, num, &
flx)
- ! convert to kg/m2/s
- flx = flx/3600.
- ! apply numerical scaling
- flx = flx*numscale
- fluxes%flx(:,:,n) = sum(flx,dim=3)/real(nread)
+ do i = 1, nread
+ if ( (n+i).gt.(nday*24/nt_flx) ) exit
+ flxghg(:,:,n+i) = flx(:,:,i)
+ timeghg(n+i) = jd + dble((i-1)*nt_flx)/24d0
+ end do
deallocate( flx )
-
- ! flux time stamp is start of time step
- fluxes%time(n) = jd
-
- ! ocean fluxes for lognormal case when inc_ocean is false
- ! cover only inversion domain
- if ( config%lognormal.and..not.config%inc_ocean.and..not.config%nested ) then
- ! use land-sea mask to set all land fluxes to zero
- ix1 = int((rllx - llx + dx)/dx)
- ix2 = int((rurx - llx)/dx)
- jy1 = int((rlly - lly + dy)/dy)
- jy2 = int((rury - lly)/dy)
- fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx(ix1:ix2,jy1:jy2,n)
- endif
-
- jd = jd + dble(statres)
-
+ n = n + nread
+ jd = jd + dble(eomday)
end do
-
+ n = min(n, nday*24/nt_flx)
+ print*, 'init_ghg: n = ',n
+ print*, 'init_ghg: timeghg = ',timeghg
+
+ ! convert to kg/m2/s
+ flxghg = flxghg/3600.
+ ! apply numerical scaling
+ flxghg = flxghg*numscale
+
+ ! average/interpolate to state vector timestamp
+ allocate( datain(nxgrid,n) )
+ allocate( dataout(nxgrid,ntstate) )
+ do i = 1, ntstate
+ fluxes%time(i) = juldatei + dble(i-1)*dble(statres)
+ enddo
+ if ( statres.gt.real(nt_flx/24) ) then
+ ! average flux
+ do jy = 1, nygrid
+ datain = flxghg(:,jy,:)
+ call average(nxgrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.gt.real(nt_flx/24) ) then
+ ! interpolate flux
+ do jy = 1, nygrid
+ datain = flxghg(:,jy,:)
+ call interp(nxgrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.eq.real(nt_flx/24) ) then
+ fluxes%flx(:,:,:) = flxghg
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( flxghg )
+ deallocate( timeghg )
+
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ ! cover only inversion domain
+ if ( config%lognormal.and..not.config%inc_ocean.and..not.config%nested ) then
+ ! use land-sea mask to set all land fluxes to zero
+ ix1 = int((rllx - llx + dx)/dx)
+ ix2 = int((rurx - llx)/dx)
+ jy1 = int((rlly - lly + dy)/dy)
+ jy2 = int((rury - lly)/dy)
+ do n = 1, ntstate
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx(ix1:ix2,jy1:jy2,n)
+ end do
+ endif
! regional prior fluxes
! ---------------------
if ( config%nested ) then
+ allocate( flxghg(nxregrid,nyregrid,(nday*24/nt_flx)), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_ghg: not enough memory'
+ allocate( timeghg(nday*24/nt_flx) )
+
n = 0
jd = juldatei
-
- do while ( (jd.le.juldatef).and.(n.lt.ntstate) )
-
- n = n + 1
+ do while ( jd.le.juldatef )
call caldate(jd, yyyymmdd, hhmiss)
write(adate,'(I4.4)') yyyymmdd/10000
-
- nread = int(statres/(real(365)/real(config%nt_flx)))
+ write(amonth,'(I2.2)') (yyyymmdd-(yyyymmdd/10000)*10000)/100
+ eomday = calceomday(yyyymmdd/100)
+ nread = eomday*24/nt_flx
if ( nread.eq.0 ) nread = 1
+ print*, 'init_ghg: reg nread = ',nread
allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_ghg: not enough memory'
filename = str_replace(files%filenest_flx, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
filename = trim(files%path_prior)//trim(filename)
write(logid,*) 'Reading prior fluxes :'//trim(filename)
inquire(file=trim(filename),exist=lexist)
@@ -196,28 +243,64 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
nxregrid, nyregrid, &
jd, nread, num, &
flx)
- ! convert to kg/m2/s
- flx = flx/3600.
- ! apply numerical scaling
- flx = flx*numscale
- fluxes%flx_nest(:,:,n) = sum(flx,dim=3)/real(nread)
+ do i = 1, nread
+ if ( (n+i).gt.(nday*24)/nt_flx ) exit
+ timeghg(n+i) = jd + dble((i-1)*nt_flx)/24d0
+ flxghg(:,:,n+i) = flx(:,:,i)
+ end do
deallocate( flx )
+ n = n + nread
+ jd = jd + dble(eomday)
+ end do
+ n = min(n, (nday*24)/nt_flx)
+ print*, 'init_ghg: reg n = ',n
+ print*, 'init_ghg: reg timeghg = ',timeghg
- ! ocean fluxes for lognormal case when inc_ocean is false
- if ( config%lognormal.and..not.config%inc_ocean ) then
- ! use land-sea mask to set all land fluxes to zero
- fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx_nest(:,:,n)
- endif
+ ! convert to kg/m2/s
+ flxghg = flxghg/3600.
+ ! apply numerical scaling
+ flxghg = flxghg*numscale
- jd = jd + dble(statres)
+ ! average/interpolate to state vector timestamp
+ allocate( datain(nxregrid,n) )
+ allocate( dataout(nxregrid,ntstate) )
+ if ( statres.gt.real(nt_flx/24) ) then
+ ! average flux
+ do jy = 1, nyregrid
+ datain = flxghg(:,jy,:)
+ call average(nxregrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.gt.real(nt_flx/24) ) then
+ ! interpolate flux
+ do jy = 1, nyregrid
+ datain = flxghg(:,jy,:)
+ call interp(nxregrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.eq.real(nt_flx/24) ) then
+ fluxes%flx_nest(:,:,:) = flxghg
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( flxghg )
+ deallocate( timeghg )
- end do
-
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ if ( config%lognormal.and..not.config%inc_ocean ) then
+ ! use land-sea mask to set all land fluxes to zero
+ do n = 1, ntstate
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx_nest(:,:,n)
+ end do
+ endif
+
+ ! nest -> use nested fluxes for state vector
allocate( flx(nxregrid,nyregrid,ntstate) )
flx(:,:,:) = fluxes%flx_nest
else
+ ! no nest -> use global fluxes for state vector
allocate( flx(nxgrid,nygrid,ntstate) )
flx(:,:,:) = fluxes%flx
@@ -334,7 +417,8 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
if ( config%lognormal ) then
! lognormal
! uncertainty of a scalar of fluxes -> use constant uncertainty
- xerr(:) = config%flxerr
+ ! transform to error of z
+ xerr(:) = log(1.+config%flxerr)
else
! normal
! average error over time steps
@@ -342,7 +426,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
endif
! assign best guess of state vector
- if ( (config%prior_bg.eq.1).and.(config%method.eq.'congrad') ) then
+ if ( (config%prior_bg.eq.1).and.((config%method.eq.'congrad').or.(config%method.eq.'m1qn3')) ) then
! from previous inversion
call read_bg(files, config, states)
else
diff --git a/source/initialize.f90 b/source/initialize.f90
index 39e640d5e5ee0692f2d0b0e49e6f372e7f1fe3b4..64e4ef441588e39ab2a97a2ee39e5375da189a7d 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -202,8 +202,22 @@ subroutine initialize(files, config)
! number days in inversion interval
nday = int(juldatef - juldatei + 1.)
- ! number of NEE fields per day
- nd_nee = config%num_nee_day
+ ! timestep of fluxes (hours)
+ nt_nee = config%nstep_nee
+ nt_nee_reg = config%nstep_nee_reg
+ nt_ff = config%nstep_ff
+ nt_ff_reg = config%nstep_ff_reg
+ nt_ocn = config%nstep_ocn
+ nt_flx = config%nstep_flx
+
+ ! number of NEE values per day
+ if ( config%spec.eq.'co2' ) then
+ if ( nt_nee_reg.eq.0 ) then
+ write(logid,*) 'ERROR: nt_nee cannot be zero'
+ stop
+ endif
+ nd_nee = 24/nt_nee_reg
+ endif
! number of initial mixing ratio time steps to optimize
ntcini = max(1, nday/cinires)
@@ -217,13 +231,13 @@ subroutine initialize(files, config)
statres = config%statres
statres_hr = config%statres_hr
ndt = int(24./statres_hr)
- ntstate = int(floor(real(nday)/statres))
+ ntstate = int(real(nday)/statres)
else
! statres is time resolution of state variables
! ntflx is number of state time invervals
! ntstate is also number of state time invervals
statres = config%statres
- ntstate = nint(real(nday)/statres)
+ ntstate = int(real(nday)/statres)
ndt = 1
endif
if ( ntstate.lt.1 ) then
@@ -322,8 +336,11 @@ subroutine initialize(files, config)
write(logid,*) 'ntstate: ',ntstate
if ( trim(config%spec).eq.'co2' ) then
write(logid,*) 'nd_nee: ',nd_nee
+ write(logid,*) 'nt_nee: ',nt_nee
+ write(logid,*) 'nt_ff: ',nt_ff
write(logid,*) 'ndt: ',ndt
endif
+ write(logid,*) 'nt_flx: ',nt_flx
write(logid,*) 'Regional domain:'
write(logid,*) 'w_edge_lon: ',rllx
write(logid,*) 's_edge_lat: ',rlly
diff --git a/source/m1qn3_interface.f90 b/source/m1qn3_interface.f90
index 58f6ea95487581ea7935f9dc8b6e6127127154c5..884d45885b280db8e61edf5528606856956e2994 100644
--- a/source/m1qn3_interface.f90
+++ b/source/m1qn3_interface.f90
@@ -67,7 +67,7 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! local variables for interface with simulate.f90
real :: cost_o, cost_p
- real(kind=8), dimension(nvar,1) :: phi, chi_ad, px, ones
+ real(kind=8), dimension(nvar,1) :: phi, chi, chi_ad, px, ones
real, dimension(nvar) :: grad_o
character(len=2) :: aiter
@@ -92,7 +92,9 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
ones(:,1) = 1d0
n = nvar
! initial guess at vector that minimizes the cost in chi space
- x(:) = 0d0
+ phi(:,1) = real(states%px,kind=8)
+ chi = phi2chi(config, states, covar, phi)
+ x(:) = chi(:,1)
! initial gradient in chi space
g(:) = grad(:,1)
! initial cost
@@ -189,8 +191,12 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! calculate cost
if ( config%lognormal ) then
! lognormal
+ ! for mode
! 2Jp = 2*sum(z) + x^Tx (where x = state vector in chi space and z is state vector in physical space)
- cost_p = 2.*sum(states%px) + real(dot_product(x,x),kind=4)
+! cost_p = 2.*sum(states%px) + real(dot_product(x,x),kind=4)
+ ! for median
+ ! 2Jp = x^Tx
+ cost_p = real(dot_product(x,x),kind=4)
else
! normal
! 2Jp = x^Tx (where x = state vector in chi space)
@@ -217,12 +223,16 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
! J' = J'p + J'o
if ( config%lognormal ) then
! lognormal
- ! J'p = B^(1/2)ones + chi
+ ! for mode (where x = state vector chi space)
+ ! J' = B^(1/2)ones + chi + J'o
! misuse phi2chi_ad to calculate B^(1/2)ones
- phi = phi2chi_ad(config, covar, ones )
- g(:) = phi(:,1) + x(:) + chi_ad(:,1)
+! phi = phi2chi_ad(config, covar, ones )
+! g(:) = phi(:,1) + x(:) + chi_ad(:,1)
+ ! for median
+ ! J' = chi + J'o
+ g(:) = x(:) + chi_ad(:,1)
else
- ! normal
+ ! normal (where x = state vector chi space)
g(:) = x(:) + chi_ad(:,1)
endif
diff --git a/source/mod_perturb.f90 b/source/mod_perturb.f90
index f454b2747557acd85cf71920efd2249b35df384d..d443d497abe3c1e5beb3e083ce72d3d1a46a3618 100644
--- a/source/mod_perturb.f90
+++ b/source/mod_perturb.f90
@@ -98,12 +98,10 @@ module mod_perturb
! add perturbation to state vector
if ( config%lognormal ) then
! lognormal
- states%px0(:) = exp(real(perturb(:),kind=4))
+ states%px(:) = real(perturb(:),kind=4)
else
! normal
states%px0(:) = states%px0(:) + real(perturb(:),kind=4)
- endif
- if ( config%prior_bg.ne.1 ) then
states%px(:) = states%px0(:)
endif
diff --git a/source/mod_regions.f90 b/source/mod_regions.f90
index 353f03bbdc1c9f6bd7315f8714115cb947582285..4c4dcd7b7d412fae4c9a928797fe9039252881e6 100644
--- a/source/mod_regions.f90
+++ b/source/mod_regions.f90
@@ -181,7 +181,7 @@ module mod_regions
open(100,file=trim(files%path_output)//'hloc_box.txt',action='write',status='replace',iostat=ierr)
do n = 1, nbox
- write(100,fmt='(F4.1)') hloc(n)
+ write(100,fmt='(F5.1)') hloc(n)
end do
close(100)
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index f2298c1031b0857964a82cbce11945aafe1fd1e9..94a61e6ae41d82dc8ae01c12216433dac871a3a9 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -503,7 +503,7 @@ module mod_save
! regrid without adding prior fluxes
if ( config%lognormal ) then
! lognormal distribution
- ! errors are for z = ln(x/xb) so errors in x = xb*(exp(err) - 1)
+ ! zerr = ln(1 + xerr) so xerr = exp(zerr) - 1 which is fraction of the flux
if ( config%nested ) then
epri(ix,jy,n) = fluxes%flx_nest(ix,jy,n)*abs(exp(states%pxerr0((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
epos(ix,jy,n) = fluxes%flx_nest(ix,jy,n)*abs(exp(states%pxerr((n-1)*nbox+nbox_xy(ix,jy))) - 1.)/numscale
@@ -756,7 +756,7 @@ module mod_save
call check( nf90_def_dim(ncid, timename, nday, time_dimid) )
call check( nf90_def_var(ncid, timename, nf90_real, time_dimid, time_varid) )
call check( nf90_put_att(ncid, time_varid, units, timeunit) )
- call check( nf90_def_dim(ncid, hourname, nd_nee, hour_dimid) )
+ call check( nf90_def_dim(ncid, hourname, ndgrid, hour_dimid) )
call check( nf90_def_var(ncid, hourname, nf90_real, hour_dimid, hour_varid) )
call check( nf90_put_att(ncid, hour_varid, units, hourunit) )
@@ -805,6 +805,7 @@ module mod_save
integer :: n, ierr
! write to file
+ ! for lognormal these are the log transformed scalars
open(100,file=trim(files%path_output)//'scalars_cini.txt',status='replace',action='write',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(',4,'(E11.4,1X))'
write(100,*) 'prior, pri_err, post, pos_err'
@@ -884,7 +885,7 @@ module mod_save
integer, dimension(3) :: dimids
integer :: varid, lon_varid, lat_varid, time_varid
- call check( nf90_create(trim(files%path_output)//filename, nf90_clobber, ncid) )
+ call check( nf90_create(trim(files%path_flexncdf)//filename, nf90_clobber, ncid) )
call check( nf90_def_dim(ncid, lonname, nx_in, lon_dimid) )
call check( nf90_def_var(ncid, lonname, nf90_real, lon_dimid, lon_varid) )
@@ -893,7 +894,7 @@ module mod_save
call check( nf90_def_var(ncid, latname, nf90_real, lat_dimid, lat_varid) )
call check( nf90_put_att(ncid, lat_varid, units, dimunit) )
call check( nf90_def_dim(ncid, timename, ngrid, time_dimid) )
- call check( nf90_def_var(ncid, timename, nf90_real, time_dimid, time_varid) )
+ call check( nf90_def_var(ncid, timename, nf90_double, time_dimid, time_varid) )
dimids = (/ lon_dimid, lat_dimid, time_dimid /)
@@ -904,6 +905,7 @@ module mod_save
call check( nf90_put_var(ncid, lon_varid, lon_in) )
call check( nf90_put_var(ncid, lat_varid, lat_in) )
+ print*, 'save_average_fp: gtime = ',gtime
call check( nf90_put_var(ncid, time_varid, gtime(1:ngrid)) )
call check( nf90_put_var(ncid, varid, grid(:,:,1:ngrid)) )
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index 32a996cc593a4a5e074d5082269e241a006b2790..2df124883200f574b16a69cfa2579d7de3927e92 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -41,6 +41,7 @@ module mod_settings
character(len=max_path_len) :: path_prior ! path to prior fluxes
character(len=max_path_len) :: path_output ! path to output
character(len=max_path_len) :: path_flexpart ! path to flexpart output
+ character(len=max_path_len) :: path_flexncdf ! path to flexpart netcdf output (from previous run)
character(len=max_path_len) :: file_bg ! prior best guess file
character(len=max_name_len) :: suffix ! observation file suffix
character(len=max_name_len) :: file_log ! log file name
@@ -116,6 +117,7 @@ module mod_settings
character(len=3) :: spec ! species ('co2' or 'ghg')
character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad' or 'm1qn3')
logical :: lognormal ! use log-normal distribution in state space
+ real :: trunc ! truncation factor for eigenvalues of B
integer :: maxiter ! max number of iterations for congrad and m1qn3
logical :: inc_ocean ! include ocean boxes in inversion (true or false)
logical :: opt_cini ! optimize background mixing ratios (true or false)
@@ -133,10 +135,16 @@ module mod_settings
real :: xres ! longitudinal grid resolution
real :: yres ! latitudinal grid resolution
logical :: regions ! spatial aggregation of grid (true or false)
- integer :: num_nee_day ! number of NEE values per 24h
+ integer :: nstep_nee ! timestep of global NEE fluxes in hours
+ integer :: nstep_nee_reg ! timestep of regional NEE fluxes in hours
+ integer :: nstep_ff ! timestep of global fossil fuel emissions in hours
+ integer :: nstep_ff_reg ! timestep of regional fossil fuel emissions in hours
+ integer :: nstep_ocn ! timestep of ocean prior fluxes in hours
+ integer :: nstep_flx ! timestep of GHG prior fluxes in hours
+ real :: coeff_ff ! coefficient to convert from input global flux units to kg/m2/h
+ real :: coeff_ff_reg ! coefficient to convert from input regional flux units to kg/m2/h
real :: statres ! temporal resolution of state vector (days)
real :: statres_hr ! sub-daily temporal resolution of state vector (e.g. 6, 12 or 24h)
- integer :: nt_flx ! number of prior flux fields per year
real :: measerr ! measurement error (ppm or ppb)
real :: cinierr ! initial concentration error (ppm or ppb)
real :: flxerr ! flux error (fraction)
@@ -376,6 +384,9 @@ module mod_settings
identifier = "path_flexpart:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%path_flexpart = cc
+ identifier = "path_flexncdf:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) files%path_flexncdf = cc
identifier = "path_output:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%path_output = cc
@@ -590,24 +601,6 @@ module mod_settings
end do read_loop
- ! default settings
- lenstr=len_trim(files%filenest_nee)
- if (lenstr == 0) then
- files%filenest_nee = files%filename_nee
- files%varnest_nee = files%varname_nee
- files%lonnest_nee = files%lonname_nee
- files%latnest_nee = files%latname_nee
- files%timenest_nee = files%timename_nee
- endif
- lenstr=len_trim(files%filenest_ff)
- if (lenstr == 0) then
- files%filenest_ff = files%filename_ff
- files%varnest_ff = files%varname_ff
- files%lonnest_ff = files%lonname_ff
- files%latnest_ff = files%latname_ff
- files%timenest_ff = files%timename_ff
- endif
-
end subroutine read_file_settings
! --------------------------------------------------
@@ -632,6 +625,7 @@ module mod_settings
config%spa_corr = .false.
config%verbose = .false.
config%lognormal = .false.
+ config%trunc = 1.e-3
! open file
open (100, file = trim (filename), status = 'old', iostat=ierr)
@@ -684,6 +678,9 @@ module mod_settings
identifier = "lognormal:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lognormal = cl
+ identifier = "trunc:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%trunc = cn
identifier = "maxiter:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%maxiter = int(cn)
@@ -749,21 +746,41 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%regions = cl
- ! read NEE time settings
- identifier = "num_nee_day:"
+ ! read prior flux time settings
+ identifier = "nstep_nee:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_nee = int(cn)
+ identifier = "nstep_nee_reg:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_nee_reg = int(cn)
+ identifier = "nstep_ff:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_ff = int(cn)
+ identifier = "nstep_ff_reg:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_ff_reg = int(cn)
+ identifier = "nstep_ocn:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_ocn = int(cn)
+ identifier = "nstep_flx:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_flx = int(cn)
+
+ ! coefficient to convert from fossil fuel input units to kg/m2/h
+ identifier = "coeff_ff:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%coeff_ff = cn
+ identifier = "coeff_ff_reg:"
call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) config%num_nee_day = int(cn)
+ if ( match ) config%coeff_ff_reg = cn
- ! read flux resolution settings
+ ! read state vector resolution settings
identifier = "statres:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%statres = real(cn,kind=4)
identifier = "statres_hr:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%statres_hr = real(cn,kind=4)
- identifier = "nt_flx:"
- call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) config%nt_flx = real(cn)
! read prior error settings
identifier = "measerr:"
diff --git a/source/mod_transform.f90 b/source/mod_transform.f90
index 98289585c57f4da2c755649d6291b7c89f38d204..0fa1038825ba41bdf62ceb9f16645974a26b947e 100644
--- a/source/mod_transform.f90
+++ b/source/mod_transform.f90
@@ -75,7 +75,11 @@ module mod_transform
end do
end do
if ( config%opt_cini ) then
- chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * config%cinierr
+ if ( config%lognormal ) then
+ chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * log(1.+config%cinierr)
+ else
+ chi2phi(npvar+1:nvar,1) = chi(npvar+1:nvar,1) * config%cinierr
+ endif
endif
chi2phi(:,1) = chi2phi(:,1) + dble(states%px0)
@@ -120,7 +124,11 @@ module mod_transform
end do
end do
if ( config%opt_cini ) then
- phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./config%cinierr
+ if ( config%lognormal ) then
+ phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./log(1.+config%cinierr)
+ else
+ phi2chi(npvar+1:nvar,1) = (phi(npvar+1:nvar,1) - states%px0(npvar+1:nvar)) * 1./config%cinierr
+ endif
endif
end function phi2chi
@@ -164,7 +172,11 @@ module mod_transform
end do
end do
if ( config%opt_cini ) then
- phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * config%cinierr
+ if ( config%lognormal ) then
+ phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * log(1.+config%cinierr)
+ else
+ phi2chi_ad(npvar+1:nvar,1) = phi_ad(npvar+1:nvar,1) * config%cinierr
+ endif
endif
end function phi2chi_ad
@@ -207,7 +219,11 @@ module mod_transform
end do
! B^(-1)*x for initial mixing ratio scalars
if ( config%opt_cini ) then
- prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./config%cinierr**2
+ if ( config%lognormal ) then
+ prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./(log(1.+config%cinierr))**2
+ else
+ prod_icov(npvar+1:nvar) = x(npvar+1:nvar)* 1./config%cinierr**2
+ endif
endif
end function prod_icov
diff --git a/source/mod_var.f90 b/source/mod_var.f90
index 62773d88c7585fd17953a046575e301b48e332c9..77144ea3b866cda08d5b996f3b9246487fc7dd59 100644
--- a/source/mod_var.f90
+++ b/source/mod_var.f90
@@ -36,7 +36,6 @@ module mod_var
integer, parameter :: maxpoint=10000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
- real, parameter :: trunc=1.e-9 ! truncation threshold for eigenvalues of error covariance matrix
real, parameter :: reqrd = 1.e-5 ! required reduction in gradient norm for congrad
! flexpart variables
@@ -94,12 +93,16 @@ module mod_var
real(kind=8) :: juldatei ! julian start date
real(kind=8) :: juldatef ! julian end date
integer :: nday ! number of days in inversion interval
- integer :: nd_nee ! number of NEE values per day (in prior NEE file)
+ integer :: nd_nee ! number of NEE values per day global fluxes
+ integer :: nt_nee ! timestep of global NEE fluxes (hours)
+ integer :: nt_nee_reg ! timestep of regional NEE fluxes (hours)
+ integer :: nt_ff ! timestep of global fossil fuel emissions (hours)
+ integer :: nt_ff_reg ! timestep of regional fossil fuel emissions (hours)
+ integer :: nt_ocn ! timestep of ocean prior fluxes (hours)
+ integer :: nt_flx ! timestep of other prior fluxes (hours)
integer :: ntstate ! number of time steps in state vector
-! integer :: ntflx ! number of time steps in non-NEE fluxes
real :: statres ! time resolution of the state vector (for CO2 1 day) (days)
real :: statres_hr ! subdaily time resolution of the state vector (for CO2 e.g. 6, 12 or 24 h)
-! real :: flxtres ! time resolution of fluxes (for GHG same as statres) (days)
! grid variables
! --------------
diff --git a/source/read_bg.f90 b/source/read_bg.f90
index 8cd803d617401a56fdc20e7ddd803d1fdb8cf781..205896b7cf7643e552c85511f816ca20aee85098 100644
--- a/source/read_bg.f90
+++ b/source/read_bg.f90
@@ -48,6 +48,7 @@ subroutine read_bg(files, config, states)
real :: area
real, dimension(:,:,:,:), allocatable :: xpos
real, dimension(:,:,:), allocatable :: xpos_box
+ integer, dimension(:), allocatable :: cnt
integer :: nlon, nlat, nt, nh
integer :: ix, jy, n
@@ -100,7 +101,13 @@ subroutine read_bg(files, config, states)
endif
! numerical scaling
- xpos = xpos*numscale
+ if ( config%lognormal ) then
+ ! xpos is x/xb
+ xpos = log(xpos)
+ else
+ ! xpos is p - p0
+ xpos = xpos*numscale
+ endif
! aggregate to boxes
if ( config%spec.eq.'co2' ) then
@@ -111,19 +118,30 @@ subroutine read_bg(files, config, states)
xpos_box(:,:,:) = 0.
ndt = 1
endif
+ allocate( cnt(nbox) )
+ cnt(:) = 0
do n = 1, ntstate
do jy = 1, nyregrid
area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).gt.0 ) then
do nt = 1, ndt
- xpos_box(nbox_xy(ix,jy),n,nt) = xpos_box(nbox_xy(ix,jy),n,nt) + xpos(ix,jy,n,nt) * area
+ if ( config%lognormal ) then
+ xpos_box(nbox_xy(ix,jy),n,nt) = xpos_box(nbox_xy(ix,jy),n,nt) + xpos(ix,jy,n,nt)
+ cnt(nbox_xy(ix,jy)) = cnt(nbox_xy(ix,jy)) + 1
+ else
+ xpos_box(nbox_xy(ix,jy),n,nt) = xpos_box(nbox_xy(ix,jy),n,nt) + xpos(ix,jy,n,nt) * area
+ endif
end do
endif
end do
end do
do nt = 1, ndt
- xpos_box(:,n,nt) = xpos_box(:,n,nt)/area_box
+ if ( config%lognormal ) then
+ where( cnt.gt.0 ) xpos_box(:,n,nt) = xpos_box(:,n,nt)/real(cnt(:))
+ else
+ xpos_box(:,n,nt) = xpos_box(:,n,nt)/area_box
+ endif
end do
end do
@@ -140,6 +158,7 @@ subroutine read_bg(files, config, states)
deallocate( xpos )
deallocate( xpos_box )
+ deallocate( cnt )
end subroutine read_bg
diff --git a/source/read_obs.f90 b/source/read_obs.f90
index e1b276071f5af64189af62408f0f946520acf863..28c5cf01e8168b4e625e71b187ec02b69fece594 100755
--- a/source/read_obs.f90
+++ b/source/read_obs.f90
@@ -144,9 +144,9 @@ subroutine read_obs(config, files)
hl = hl - 24
endif
if ( alt.le.1000. ) then
-! if ( (hl.lt.11).or.(hl.gt.17) ) cycle read_loop
+! if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
else
-! if ( (hl.gt.4).and.(hl.lt.22) ) 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/retrieval.f90 b/source/retrieval.f90
index 93be9898226634d6442a5fe7b5c499ee2b6fda83..596245fb48eaadbe653d7d818ae62dbccf530343 100644
--- a/source/retrieval.f90
+++ b/source/retrieval.f90
@@ -100,9 +100,13 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
! gradient state space
if ( config%lognormal ) then
- ! lognormal (note p is log transformed variable)
- ! J'p = ones + B^(-1)(p - p0)
- grad_p = ones(:,1) + prod_icov(config, covar, states%px - states%px0)
+ ! lognormal (note p is the log transformed variable z)
+ ! for mode
+ ! J'p = ones + B^(-1)p
+! grad_p = ones(:,1) + prod_icov(config, covar, states%px)
+ ! for median
+ ! J'p = B'(-1)p
+ grad_p = prod_icov(config, covar, states%px)
else
! normal
! J'p = B^(-1)(p - p0)
@@ -111,10 +115,14 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
! cost state space
if ( config%lognormal ) then
- ! lognormal (note p is log transformed variable)
- ! 2Jp = 2*sum(p) + (p - p0)^TB^(-1)(p - p0)
- icovp = prod_icov(config, covar, states%px - states%px0)
- cost_p = 2.*sum(states%px) + dot_product( states%px - states%px0, icovp )
+ ! lognormal (note p is log transformed variable z)
+ ! for mode
+ ! 2Jp = 2*sum(p) + p^TB^(-1)p
+! icovp = prod_icov(config, covar, states%px)
+! cost_p = 2.*sum(states%px) + dot_product( states%px, icovp )
+ ! for median
+ ! 2Jp = p^TB^(-1)p
+ cost_p = dot_product( states%px, grad_p )
else
! normal
! 2Jp = (p - p0)^TB^(-1)(p - p0)
@@ -187,10 +195,15 @@ subroutine retrieval(files, config, fluxes, obs, states, covar)
! cost state space
if ( config%lognormal ) then
- ! lognormal (note p is log transformed variable)
- ! 2Jp = 2*sum(p) + (p - p0)^TB^(-1)(p - p0)
- icovp = prod_icov(config, covar, states%px - states%px0)
- cost_p = 2.*sum(states%px) + dot_product( states%px - states%px0, icovp )
+ ! lognormal (note p is log transformed variable z)
+ ! for mode
+ ! 2Jp = 2*sum(p) + p^TB^(-1)p
+! icovp = prod_icov(config, covar, states%px )
+! cost_p = 2.*sum(states%px) + dot_product( states%px, icovp )
+ ! for median
+ ! 2Jp = p^TB^(-1)p
+ grad_p = prod_icov(config, covar, states%px)
+ cost_p = dot_product( states%px, real(grad_p,kind=4) )
else
! normal
! 2Jp = (p - p0)^TB^(-1)(p - p0)
diff --git a/source/sbatch_flexinvert.sh b/source/sbatch_flexinvert.sh
index 3b1f2266be192e40de4db88445b8a8d481043bf3..c29c2f760a96a7a4cd0cc11bb7168140bd5d80ba 100755
--- a/source/sbatch_flexinvert.sh
+++ b/source/sbatch_flexinvert.sh
@@ -1,9 +1,9 @@
#!/bin/bash
#---------------------------------------------------
partition=nilu
-jobname=test
-settings_files='/mypath/settings/SETTINGS_ghg_files'
-settings_config='/mypath/settings/SETTINGS_ghg_config'
+jobname=testco2
+settings_files='/home/rthompson/REPOS/GITHUB/FLEXINVERTplus_clean/settings/SETTINGS_co2_files'
+settings_config='/home/rthompson/REPOS/GITHUB/FLEXINVERTplus_clean/settings/SETTINGS_co2_config'
#---------------------------------------------------
cat <<EOF > run_job.sh
diff --git a/source/simulate.f90 b/source/simulate.f90
index d6f3ed2fcb810b6d47bd8b6b9be20169579b299b..d55f4c6d35863e1d886b2a76dee5452946ea2079 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -185,7 +185,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
else
! footprints averaged/selected to match observation timestamps
subdir_fp = trim(obs%recs(i))//'/'//trim(adate)//'/'
- filename = trim(files%path_output)//trim(subdir_fp)//'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc'
+ filename = trim(files%path_flexncdf)//trim(subdir_fp)//'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc'
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filename)
@@ -228,7 +228,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
else
! footprints averaged/selected to match observation timestamps
subdir_fp = trim(obs%recs(i))//'/'//trim(adate)//'/'
- filenest = trim(files%path_output)//trim(subdir_fp)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001.nc'
+ filenest = trim(files%path_flexncdf)//trim(subdir_fp)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001.nc'
inquire ( file=trim(filenest),exist=lexist )
if ( lexist ) then
! read averaged footprints from ncdf files
@@ -277,16 +277,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( .not.allocated(hbox) ) allocate( hbox(nbox*ngrid) )
hbox(:) = 0.
do n = 1, ngrid
-! ! vectorise grid
-! do jy = 1, nyregrid
-! if ( config%nested ) then
-! gridnestvec((jy-1)*nxregrid+1:jy*nxregrid) = gridnest(ixn1:(ixn1+nxregrid-1),jyn1+jy-1,n)
-! else
-! gridnestvec((jy-1)*nxregrid+1:jy*nxregrid) = grid(ix1:(ix1+nxregrid-1),jy1+jy-1,n)
-! endif
-! end do
-! ! map to irregular grid
-! hbox((n-1)*nbox+1:n*nbox) = matmul(gridoper, gridnestvec)
do jy = 1, nyregrid
do ix = 1, nxregrid
if ( nbox_xy(ix,jy).gt.0 ) then
@@ -303,6 +293,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if (i.lt.10.and.iter.eq.0) write(logid,*) 'simulate: time to calc hbox (ms) = ',(time2-time1)*1000/countrate
! transport for grid regional domain
+ ! to calculate contribution of fixed fluxes in domain
if ( .not.allocated(hnest) ) allocate( hnest(nxregrid,nyregrid,ngrid) )
if ( config%nested ) then
hnest(:,:,:) = gridnest(ixn1:ixn2,jyn1:jyn2,1:ngrid)
@@ -350,8 +341,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! ----------
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
- ! select state variables for current footprint times and
- ! average transport operator to state variable temporal resolution
+ ! select state variables for current footprint times
allocate( px(ntstep*ndt*nbox), stat=ierr )
if( ierr.ne.0 ) then
write(logid,*) 'ERROR simulate: not enough memory'
@@ -452,14 +442,13 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! model mixing ratios
! -------------------
- if ( iter.eq.0.or.iter.eq.99 ) then
+ if ( iter.eq.0.or.(config%restart.eq.1.and.iter.eq.(lastiter+1)) ) then
call system_clock(time1, countrate)
! compute fixed mixing ratios
call calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, i, ix1, ix2, jy1, jy2)
call system_clock(time2, countrate)
if (i.eq.1.and.iter.eq.0) write(logid,*) 'simulate: time to calc conc (ms) = ',(time2-time1)*1000/countrate
endif ! iter
- if ( i.lt.20 ) print*, 'simulate: obs%ghg(i) = ',obs%ghg(i)
! contribution from state vector fluxes
if ( config%method.eq.'congrad'.or.config%method.eq.'m1qn3' ) then
@@ -601,7 +590,8 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
filename = trim(files%path_output)//'grad_'//aiter//'.txt'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
do i = 1, nvar
- write(100,fmt='(E11.4)') grad_o(i)
+! write(100,fmt='(E11.4)') grad_o(i)
+ write(100,fmt='(E13.6)') grad_o(i)
end do
close(100)
endif