Skip to content
Snippets Groups Projects
Commit 03aa365c authored by ronesy's avatar ronesy
Browse files

changes to find median solution for log-normal PDF in state space and changes... 

changes to find median solution for log-normal PDF in state space and changes to reading of prior fluxes
parent c1c9fd9e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 143 additions and 120 deletions
FLEXINVERT_userguide.docx deleted 100644 → 0
+ 0
− 0
View file @ c1c9fd9e
File deleted
grid_to_ncdf/SETTINGS
+ 2
− 2
View file @ 03aa365c
...@@ -19,8 +19,8 @@ file_recept: /myfile.txt ...@@ -19,8 +19,8 @@ file_recept: /myfile.txt
# Run dates: format yyyymmdd # Run dates: format yyyymmdd
# datei = start date # datei = start date
# datef = end date # datef = end date
datei: 20120101 datei: 20140701
datef: 20120131 datef: 20140831
# Nested output (true/false) # Nested output (true/false)
lnested: .true. lnested: .true.
......
prep_flexpart/process_obs.f90
+ 15
− 15
View file @ 03aa365c
...@@ -94,7 +94,7 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc ...@@ -94,7 +94,7 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
lonobs(1:nobs) = lonobs_out lonobs(1:nobs) = lonobs_out
altobs(1:nobs) = altobs_out altobs(1:nobs) = altobs_out
errobs(1:nobs) = errobs_out errobs(1:nobs) = errobs_out
print*, 'process_obs: ind = ',ind ! print*, 'process_obs: ind = ',ind
! if average and/or select obs ! if average and/or select obs
if( settings%laverage.eq.1.or.settings%lselect.eq.1 ) then 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 ...@@ -106,8 +106,8 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
eomday = calceomday(jjjjmmdd/100) eomday = calceomday(jjjjmmdd/100)
jdsom = juldate((jjjjmmdd/100)*100+1, 0) jdsom = juldate((jjjjmmdd/100)*100+1, 0)
jdeom = jdsom + real(eomday,kind=8) jdeom = jdsom + real(eomday,kind=8)
print*, 'jdsom: ',jdsom ! print*, 'jdsom: ',jdsom
print*, 'jdeom: ',jdeom ! print*, 'jdeom: ',jdeom
! loop over all observations in current month ! loop over all observations in current month
do while( jdate.le.jdobs(nobs) ) do while( jdate.le.jdobs(nobs) )
jdt = dnint(jdate*1.e6)/1.e6 jdt = dnint(jdate*1.e6)/1.e6
...@@ -141,23 +141,23 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc ...@@ -141,23 +141,23 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
jdateend = min(jdateend,jdeom) jdateend = min(jdateend,jdeom)
! check that not beyond end of run time ! check that not beyond end of run time
jdateend = min(jdateend,jdatef) jdateend = min(jdateend,jdatef)
print*, 'jdate: ',jdate ! print*, 'jdate: ',jdate
print*, 'jdatestart: ',jdatestart ! print*, 'jdatestart: ',jdatestart
print*, 'jdateend: ',jdateend ! print*, 'jdateend: ',jdateend
cnt = 0 cnt = 0
buffer(:,:) = 0. buffer(:,:) = 0.
! loop over observations in averaging interval (fill buffer) ! loop over observations in averaging interval (fill buffer)
do while( (jdt.ge.jdatestart).and.(jdt.lt.jdateend) ) do while( (jdt.ge.jdatestart).and.(jdt.lt.jdateend) )
call caldate(jdt, jjjjmmdd, hhmiss) call caldate(jdt, jjjjmmdd, hhmiss)
print*, 'jdt:',jdt ! print*, 'jdt:',jdt
print*, 'jjjjmmdd, hhmiss:',jjjjmmdd,hhmiss ! print*, 'jjjjmmdd, hhmiss:',jjjjmmdd,hhmiss
if( settings%lselect.eq.1 ) then if( settings%lselect.eq.1 ) then
! select obs ! select obs
hloc = hhmiss/10000 + int(lonobs(n)*24./360.) 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.gt.24) hloc = hloc - 24
if( hloc.lt.0) hloc = hloc + 24 if( hloc.lt.0) hloc = hloc + 24
print*, 'hloc:',hloc ! print*, 'hloc:',hloc
if( altobs(n).lt.1000. ) then if( altobs(n).lt.1000. ) then
! select daytime (12-16h) ! select daytime (12-16h)
if( hloc.ge.12.and.hloc.le.16) then if( hloc.ge.12.and.hloc.le.16) then
...@@ -194,17 +194,17 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc ...@@ -194,17 +194,17 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
nprev = n nprev = n
n = n + 1 n = n + 1
if ( n.gt.nobs ) go to 10 if ( n.gt.nobs ) go to 10
print*, 'jdt: ',jdt ! print*, 'jdt: ',jdt
! jdt = jdobs(n) ! jdt = jdobs(n)
jdt = dnint(jdobs(n)*1.e6)/1.e6 jdt = dnint(jdobs(n)*1.e6)/1.e6
print*, 'next jdt: ',jdt ! print*, 'next jdt: ',jdt
end do end do
if( nprev.eq.0 ) nprev = 1 if( nprev.eq.0 ) nprev = 1
10 continue 10 continue
print*, 'cnt:',cnt ! print*, 'cnt:',cnt
print*, 'buffer:' ! print*, 'buffer:'
do c = 1, cnt 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 end do
if( settings%laverage.eq.1 ) then if( settings%laverage.eq.1 ) then
! average obs in buffer ! average obs in buffer
......
prep_flexpart/read_basic.f90
+ 3
− 3
View file @ 03aa365c
...@@ -89,7 +89,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -89,7 +89,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
! find file matching receptor name ! find file matching receptor name
do n = 1, nfiles do n = 1, nfiles
do i = 1, len_trim(filelist(n))-2 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
end do end do
10 continue 10 continue
...@@ -151,8 +151,8 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -151,8 +151,8 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
jdobs(cnt) = jdate jdobs(cnt) = jdate
concobs(cnt) = conc concobs(cnt) = conc
errobs(cnt) = err errobs(cnt) = err
lonobs(cnt) = temp(9) latobs(cnt) = temp(9)
latobs(cnt) = temp(10) lonobs(cnt) = temp(10)
altobs(cnt) = temp(11) altobs(cnt) = temp(11)
if ( errobs(cnt).le.-9.99 ) errobs(cnt) = 0. if ( errobs(cnt).le.-9.99 ) errobs(cnt) = 0.
end do read_loop end do read_loop
......
prep_flexpart/run_flexpart.sh
+ 4
− 6
View file @ 03aa365c
...@@ -14,11 +14,9 @@ ...@@ -14,11 +14,9 @@
# To run tests reference to current directory # To run tests reference to current directory
DIR=$PWD DIR=$PWD
LEN=`expr length $DIR - 11` LEN=`expr length $DIR - 13`
ROOTDIR=${DIR:0:$LEN} ROOTDIR=${DIR:0:$LEN}
echo ${ROOTDIR}
#--------------------------------------------------- #---------------------------------------------------
# User settings # User settings
...@@ -28,9 +26,9 @@ TIMELIM=12:00:00 ...@@ -28,9 +26,9 @@ TIMELIM=12:00:00
PARTITION="nilu" PARTITION="nilu"
EXENAME=FP_ecmwf_gfortran EXENAME=FP_ecmwf_gfortran
METEO=ecmwf METEO=ecmwf
DIRFLEX=${ROOTDIR}TEST_CASE_CH4/FLEXPART/ DIRFLEX=/home/rthompson/REPOS/GITHUB/FLEXPART/
DIROPTIONS=${ROOTDIR}TEST_CASE_CH4/FP_FILES/ DIROPTIONS=${ROOTDIR}TEST_INPUT/FLEXPART/GHG/NEST/
STNLIST=(PAL) STNLIST=(SSL)
MONLIST=(01) MONLIST=(01)
YEAR=2012 YEAR=2012
......
prep_flexpart/sbatch_prep_flexpart.sh
+ 0
− 2
View file @ 03aa365c
...@@ -4,8 +4,6 @@ partition=debug ...@@ -4,8 +4,6 @@ partition=debug
settings_files='./SETTINGS_ghg' settings_files='./SETTINGS_ghg'
#--------------------------------------------------- #---------------------------------------------------
cp ../../prep_flexpart/prep_flexpart .
cat <<EOF > run_job.sh cat <<EOF > run_job.sh
#!/bin/bash #!/bin/bash
./prep_flexpart ${settings_files} ./prep_flexpart ${settings_files}
......
prep_fluxes/FLUXES_ghg
+ 6
− 6
View file @ 03aa365c
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
# ntime = number of time steps # ntime = number of time steps
# timeref = ref date if timestamp is julian days (yyyymmdd) otherwise 0 # 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) # 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 varname_in: emisch4
lonname_in: longitude lonname_in: longitude
latname_in: latitude latname_in: latitude
...@@ -38,17 +38,17 @@ timestamp: ...@@ -38,17 +38,17 @@ timestamp:
# ntime_out = number of time steps # ntime_out = number of time steps
# llx_out = left longitude of domain # llx_out = left longitude of domain
# lly_out = lower latitude 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 varname_out: emisch4
varunit: kgC/m2/h varunit: kgC/m2/h
lonname_out: longitude lonname_out: longitude
latname_out: latitude latname_out: latitude
timename_out: time timename_out: time
year: 2012 year: 2012
nx_out: 180 nx_out: 360
ny_out: 90 ny_out: 180
dx_out: 2.0 dx_out: 1.0
dy_out: 2.0 dy_out: 1.0
ntime_out: 12 ntime_out: 12
llx_out: -180 llx_out: -180
lly_out: -90 lly_out: -90
......
prep_fluxes/adjust_time.f90
+ 1
− 1
View file @ 03aa365c
...@@ -53,7 +53,7 @@ subroutine adjust_time(settings, time, time_out) ...@@ -53,7 +53,7 @@ subroutine adjust_time(settings, time, time_out)
! input timestamp is julian days ! input timestamp is julian days
do i = 1, ntime_out do i = 1, ntime_out
time_out(i) = time(i) + juldate(int(settings%timeref),0) - & time_out(i) = time(i) + juldate(int(settings%timeref),0) - &
juldate(19000101,0) + 1. juldate(19000101,0)
end do end do
! adjust time stamp to given year ! adjust time stamp to given year
call caldate(time_out(1) + juldate(19000101,0), yyyymmdd, hhmiss) call caldate(time_out(1) + juldate(19000101,0), yyyymmdd, hhmiss)
......
prep_fluxes/convertgrid.f90
+ 11
− 4
View file @ 03aa365c
...@@ -111,9 +111,12 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, & ...@@ -111,9 +111,12 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
end do end do
dy_in(1) = dy_in(2) dy_in(1) = dy_in(2)
! round dx and dy to 7dp ! round dx and dy to 4dp
dx_in = anint(dx_in*1.e7)/1.e7 dx_in = anint(dx_in*1.e4)/1.e4
dy_in = anint(dy_in*1.e7)/1.e7 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 check if west and south boundaries given
if(.not.midpoints) then if(.not.midpoints) then
...@@ -159,6 +162,10 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, & ...@@ -159,6 +162,10 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
end do end do
endif 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 if( (nx_in.eq.nx_out) .and. (ny_in.eq.ny_out) ) then
write(*,*) 'Comment: input and output grids are the same' write(*,*) 'Comment: input and output grids are the same'
wlon_out = wlon_in wlon_out = wlon_in
...@@ -268,7 +275,7 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, & ...@@ -268,7 +275,7 @@ subroutine convertgrid(midpoints, lon_in, lat_in, flx_tmp, &
! calculate total in ! calculate total in
do nt = 1, ntime_out 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 end do
print*, tot_in*24.*365./float(ntime_out)/1.e9 print*, tot_in*24.*365./float(ntime_out)/1.e9
......
prep_fluxes/rev_test.out deleted 100644 → 0
+ 0
− 41
View file @ c1c9fd9e
/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
prep_regions/SETTINGS_regions
+ 7
− 3
View file @ 03aa365c
...@@ -7,9 +7,6 @@ ...@@ -7,9 +7,6 @@
# #
# ================================================== # ==================================================
# Define map globally with "zoom" (true) or regionally (false)
usezoom: .false.
# Use prior fluxes # Use prior fluxes
# weight surface influnce by annual fluxes (true/false) # weight surface influnce by annual fluxes (true/false)
useflux: .true. useflux: .true.
...@@ -33,3 +30,10 @@ nmap: 3 ...@@ -33,3 +30,10 @@ nmap: 3
# 3) 1 x 0.5 = 0.5 degree # 3) 1 x 0.5 = 0.5 degree
mapdec: 1,2,2 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:
prep_regions/SETTINGS_zoom_regions
+ 7
− 0
View file @ 03aa365c
...@@ -43,3 +43,10 @@ zlly: 38 ...@@ -43,3 +43,10 @@ zlly: 38
zurx: 28 zurx: 28
zury: 70 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:
prep_regions/map_grid.f90
+ 24
− 1
View file @ 03aa365c
...@@ -78,7 +78,7 @@ subroutine map_grid(config, flux, surfinf, nbox_xy) ...@@ -78,7 +78,7 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
integer :: numx, numy, num integer :: numx, numy, num
integer :: nres, ndec, nfx, nfy, nfbox, nsort integer :: nres, ndec, nfx, nfy, nfbox, nsort
integer, dimension(2*maxbox) :: map_box 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(2*maxbox) :: x_map, y_map, sens_map, water_box
real, dimension(nxregrid,nyregrid) :: sens_lump, contrib_lump, water real, dimension(nxregrid,nyregrid) :: sens_lump, contrib_lump, water
real, dimension(nxregrid*nyregrid) :: sort_contrib, sort_contrib_inv real, dimension(nxregrid*nyregrid) :: sort_contrib, sort_contrib_inv
...@@ -258,6 +258,29 @@ subroutine map_grid(config, flux, surfinf, nbox_xy) ...@@ -258,6 +258,29 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
print*, 'nbox_invert, nb = ',nbox_invert, nb 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) ! renumber boxes for land (1:nbox_land) and water (-1:-nbox_water)
nbox_water = 0 nbox_water = 0
nbox_land = 0 nbox_land = 0
......
prep_regions/mod_settings.f90
+ 23
− 22
View file @ 03aa365c
...@@ -90,7 +90,7 @@ module mod_settings ...@@ -90,7 +90,7 @@ module mod_settings
real :: n_edge_lat ! lat of northern edge of inversion grid real :: n_edge_lat ! lat of northern edge of inversion grid
real :: xres ! longitudinal grid resolution real :: xres ! longitudinal grid resolution
real :: yres ! latitudinal 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 end type config_t
...@@ -289,9 +289,6 @@ module mod_settings ...@@ -289,9 +289,6 @@ module mod_settings
logical :: match, cl logical :: match, cl
! default values
files%file_ice = ""
! open file ! open file
open (100, file = trim (filename), status = 'old', iostat=ierr) open (100, file = trim (filename), status = 'old', iostat=ierr)
if(ierr.gt.0) then if(ierr.gt.0) then
...@@ -363,20 +360,6 @@ module mod_settings ...@@ -363,20 +360,6 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%latname_lsm = cc 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 ! read flux prior settings
identifier = "filename_flx:" identifier = "filename_flx:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
...@@ -523,9 +506,9 @@ module mod_settings ...@@ -523,9 +506,9 @@ module mod_settings
identifier = "yres:" identifier = "yres:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%yres = real(cn) if ( match ) config%yres = real(cn)
identifier = "nt_flx:" identifier = "nstep_flx:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%nt_flx = int(cn) if ( match ) config%nstep_flx = int(cn)
endif endif
...@@ -537,9 +520,10 @@ module mod_settings ...@@ -537,9 +520,10 @@ module mod_settings
! read region 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 integer :: ierr, i, n
character(len=200) :: line, identifier, cc character(len=200) :: line, identifier, cc
...@@ -554,6 +538,9 @@ module mod_settings ...@@ -554,6 +538,9 @@ module mod_settings
stop stop
endif endif
! default values
files%file_ice = ""
! loop over lines ! loop over lines
read_loop: do read_loop: do
read (100, fmt='(A)', iostat=ierr) line read (100, fmt='(A)', iostat=ierr) line
...@@ -610,6 +597,20 @@ module mod_settings ...@@ -610,6 +597,20 @@ module mod_settings
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) zury = int(cn) 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 endif
end do read_loop end do read_loop
......
prep_regions/mod_var.f90
+ 5
− 2
View file @ 03aa365c
...@@ -35,7 +35,7 @@ module mod_var ...@@ -35,7 +35,7 @@ module mod_var
integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month 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 :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output 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 ! flexpart variables
! ------------------ ! ------------------
...@@ -101,7 +101,10 @@ module mod_var ...@@ -101,7 +101,10 @@ module mod_var
integer :: nrec ! number of receptor sites or campaigns integer :: nrec ! number of receptor sites or campaigns
integer :: nobs ! total number of observations integer :: nobs ! total number of observations
character(recname_len), dimension(:), allocatable :: recname ! names of receptors 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 end module mod_var
prep_regions/prep_regions.f90
+ 1
− 1
View file @ 03aa365c
...@@ -78,7 +78,7 @@ program prep_regions ...@@ -78,7 +78,7 @@ program prep_regions
call read_config_settings(settings_config, config) call read_config_settings(settings_config, config)
! read region settings ! read region settings
call read_region_settings(settings_regions) call read_region_settings(settings_regions, files)
! initialize parameters ! initialize parameters
call initialize(files, config) call initialize(files, config)
......
prep_regions/read_lsm.f90
+ 1
− 1
View file @ 03aa365c
...@@ -175,7 +175,7 @@ subroutine read_lsm(files) ...@@ -175,7 +175,7 @@ subroutine read_lsm(files)
! integrate ice into lsm ! integrate ice into lsm
ice_av = sum(ice,dim=3)/12. 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(lon)
deallocate(lat) deallocate(lat)
......
prep_regions/read_obs.f90
+ 30
− 5
View file @ 03aa365c
...@@ -53,10 +53,10 @@ subroutine read_obs(files, reclist, obstimes, avetimes) ...@@ -53,10 +53,10 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
integer :: ierr integer :: ierr
integer :: cnt integer :: cnt
integer :: narg, n, reclen integer :: narg, n, reclen
integer :: nf, nfiles integer :: nf, nfiles, nr
integer :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss integer :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss, hl
real(kind=8) :: jdate, avetime real(kind=8) :: jdate, avetime
real :: conc, err real :: conc, err, lon, alt
integer :: num integer :: num
! list observation files ! list observation files
...@@ -108,6 +108,16 @@ subroutine read_obs(files, reclist, obstimes, avetimes) ...@@ -108,6 +108,16 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
endif endif
write(*,*) 'Reading file: '//trim(files%path_obs)//trim(filelist(nf)) 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 header
read (100, fmt='(A)', iostat=ierr) header read (100, fmt='(A)', iostat=ierr) header
print*, 'read_obs: header = ',header print*, 'read_obs: header = ',header
...@@ -122,19 +132,34 @@ subroutine read_obs(files, reclist, obstimes, avetimes) ...@@ -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) & 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 recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
endif endif
! if obs file doesn't contain avetime (backwards compatability) ! if obs file doesn't contain avetime (backwards compatability)
! avetime = 0d0 ! avetime = 0d0
! if ( reclen.eq.4 ) then ! 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) & ! 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 ! recs, yyyymmdd, hhmiss, jdate, conc, err, num
! else ! 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 ! recs, yyyymmdd, hhmiss, jdate, conc, err, num
! endif ! endif
if ( ierr.gt.0 ) exit read_loop if ( ierr.gt.0 ) exit read_loop
if ( jdate.ge.(juldatef+1d0) ) exit read_loop if ( jdate.ge.(juldatef+1d0) ) exit read_loop
if ( jdate.lt.juldatei ) cycle read_loop if ( jdate.lt.juldatei ) cycle read_loop
if ( conc.le.-999. ) 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 cnt = cnt + 1
obstimes(cnt) = jdate obstimes(cnt) = jdate
avetimes(cnt) = avetime avetimes(cnt) = avetime
......
prep_regions/read_reclist.f90
+ 0
− 1
View file @ 03aa365c
...@@ -33,7 +33,6 @@ subroutine read_reclist(filename) ...@@ -33,7 +33,6 @@ subroutine read_reclist(filename)
character(len=200) :: line character(len=200) :: line
integer :: ierr integer :: ierr
integer :: cnt integer :: cnt
real, dimension(:), allocatable :: reclat, reclon, recalt
! count number of receptors ! count number of receptors
......
prep_syndata/README_prepsyndata.txt
+ 3
− 4
View file @ 03aa365c
================================================================ ================================================================
FLEXINVERT-Plus SYNTHETIC DATA PREPARATION FLEXINVERT SYNTHETIC DATA PREPARATION
================================================================ ================================================================
...@@ -9,7 +9,7 @@ Description: ...@@ -9,7 +9,7 @@ Description:
Notes: Notes:
1) Output from a forward run of FLEXINVERT (run_mode = 0) 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. 2) For CO2 only NEE fluxes are perturbed.
3) Fluxes are returned on the same domain as the flexpart 3) Fluxes are returned on the same domain as the flexpart
grid_time files (if nested is true this is the nested grid_time files (if nested is true this is the nested
...@@ -17,8 +17,7 @@ Notes: ...@@ -17,8 +17,7 @@ Notes:
greater domain than the flexpart files. greater domain than the flexpart files.
4) If the inversion domain (covered by state vector) is 4) If the inversion domain (covered by state vector) is
smaller than the flexpart domain, only the fluxes inside smaller than the flexpart domain, only the fluxes inside
the inversion domain are perturbed (fluxes outside this the inversion domain are perturbed.
domain are equal to the prior).
Synthetic Data: Synthetic Data:
1) Prior fluxes: prepared by randomly perturbing prior fluxes 1) Prior fluxes: prepared by randomly perturbing prior fluxes
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment