...
 
Commits (10)
......@@ -13,6 +13,15 @@ path_flexpart: /mypath/
# Path where to write ncdf files
path_output: /mypath/
# Receptor list file
file_recept: /myfile.txt
# Run dates: format yyyymmdd
# datei = start date
# datef = end date
datei: 20140701
datef: 20140831
# Nested output (true/false)
lnested: .true.
......
......@@ -23,7 +23,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine grid_convert(files, config)
subroutine grid_convert(files, config, nr, amonth)
use mod_var
use mod_settings
......@@ -34,28 +34,40 @@ subroutine grid_convert(files, config)
type (files_t), intent(in out) :: files
type (config_t), intent(in) :: config
integer, intent(in) :: nr
character(max_path_len) :: filename, filedates
character(max_path_len) :: filename, filedates, path_flexpart, path_output
character(len=max_path_len), dimension(:), allocatable :: filelist
character(len=8) :: adate
character(len=6) :: atime
character(len=6) :: atime, amonth
logical :: lexist
integer :: nfiles, i, ierr, n
integer :: jjjjmmdd, hhmiss
real, dimension(nxgrid,nygrid,ngrid) :: grid
real(kind=8), dimension(ngrid) :: gtime, time
real, dimension(:,:,:), allocatable :: grid
real(kind=8), dimension(:), allocatable :: gtime, time
real(kind=8) :: jdate
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely
integer :: numx, numy
! list all grid files in directory
path_flexpart = trim(files%path_flexpart)//trim(recname(nr))//'/'//trim(amonth)//'/'
path_output = trim(files%path_output)//trim(recname(nr))//'/'
print*, path_flexpart
print*, path_output
call system('mkdir '//trim(path_output))
if ( config%lnested ) then
call system('ls '//trim(files%path_flexpart)//' | grep grid_time_nest | wc -l > '//trim(files%path_output)//'files.txt')
call system('ls '//trim(files%path_flexpart)//' | grep grid_time_nest >> '//trim(files%path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' | grep grid_time_nest_ | wc -l > '//trim(path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' | grep grid_time_nest_ >> '//trim(path_output)//'files.txt')
else
call system('ls '//trim(files%path_flexpart)//' -I "*nest*" | grep grid_time | wc -l > '//trim(files%path_output)//'files.txt')
call system('ls '//trim(files%path_flexpart)//' -I "*nest*" | grep grid_time >> '//trim(files%path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' -I "*nest*" | grep grid_time_ | wc -l > '//trim(path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' -I "*nest*" | grep grid_time_ >> '//trim(path_output)//'files.txt')
endif
! read file list
open(100,file=trim(files%path_output)//'files.txt',status="old",action="read",iostat=ierr)
open(100,file=trim(path_output)//'files.txt',status="old",action="read",iostat=ierr)
if ( ierr.ne.0 ) then
print*, 'ERROR: problem read files.txt'
stop
......@@ -68,8 +80,32 @@ subroutine grid_convert(files, config)
end do
close(100)
! read header for this month to get ngrid
if ( config%lnested ) then
filename = trim(path_flexpart)//'header_nest'
else
filename = trim(path_flexpart)//'header'
endif
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then
print*, 'ERROR: cannot find flexpart header'
stop
endif
print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
if ( config%linversionout ) then
ngrid = trajdays*ndgrid+2
else
ngrid = numpoint
endif
print*, 'ngrid = ',ngrid
allocate( grid(nxgrid,nygrid,ngrid) )
allocate( gtime(ngrid) )
allocate( time(ngrid) )
! loop over files
filedates = trim(files%path_flexpart)//'dates'
filedates = trim(path_flexpart)//'dates'
do i = 1, nfiles
! read binary grid files
if ( config%lnested ) then
......@@ -79,25 +115,25 @@ subroutine grid_convert(files, config)
adate = filelist(i)(11:18)
atime = filelist(i)(19:24)
endif
print*, 'adate = ',adate
print*, 'atime = ',atime
! print*, 'adate = ',adate
! print*, 'atime = ',atime
read(adate,*) jjjjmmdd
read(atime,*) hhmiss
print*, 'jjjjmmdd = ',jjjjmmdd
print*, 'hhmiss = ',hhmiss
! print*, 'jjjjmmdd = ',jjjjmmdd
! print*, 'hhmiss = ',hhmiss
if ( config%linversionout ) then
jdate = juldate(jjjjmmdd, hhmiss)
else
jdate = juldate(datei,timei)
endif
print*, 'jdate = ',jdate
filename = trim(files%path_flexpart)//trim(filelist(i))
! print*, 'jdate = ',jdate
filename = trim(path_flexpart)//trim(filelist(i))
if ( config%linversionout ) then
call read_grid(config, filename, filedates, jdate, nxgrid, nygrid, grid, gtime)
else
call read_grid_std(config, filename, jdate, nxgrid, nygrid, grid, gtime)
endif
print*, 'gtime = ',gtime
! print*, 'gtime = ',gtime
! write to ncdf
do n = 1, ngrid
call caldate(gtime(n), jjjjmmdd, hhmiss)
......@@ -105,13 +141,16 @@ subroutine grid_convert(files, config)
end do
print*, 'time = ',time
if ( config%lnested ) then
filename = trim(files%path_output)//'grid_time_nest_'//adate//atime//'_001.nc'
filename = trim(path_output)//'grid_time_nest_'//adate//atime//'_001.nc'
else
filename = trim(files%path_output)//'grid_time_'//adate//atime//'_001.nc'
filename = trim(path_output)//'grid_time_'//adate//atime//'_001.nc'
endif
call write_ncdf(filename, time, grid)
end do
deallocate( grid )
deallocate( gtime )
deallocate( time )
end subroutine grid_convert
......@@ -27,6 +27,7 @@ program grid_to_ncdf
use mod_flexpart
use mod_var
use mod_dates
use mod_settings
implicit none
......@@ -36,6 +37,10 @@ program grid_to_ncdf
type ( config_t ) :: config
character(max_path_len) :: settings_files
character(len=6) :: adate
integer :: nr
integer :: yyyymmdd, hhmiss, eomday
real(kind=8) :: juldatei, juldatef, jd
call getarg(1,settings_files)
if (settings_files == '') then
......@@ -47,9 +52,20 @@ program grid_to_ncdf
! initialize parameters
call initialize(files, config)
juldatei = juldate(config%datei, 0)
juldatef = juldate(config%datef, 0)
! convert grid files
call grid_convert(files, config)
do nr = 1, nrec
jd = juldatei
do while ( jd.le.juldatef )
call caldate(jd, yyyymmdd, hhmiss)
write(adate,'(I6.6)') yyyymmdd/100
call grid_convert(files, config, nr, adate)
eomday = calceomday(yyyymmdd/100)
jd = jd + dble(eomday)
end do
end do
end program
......
......@@ -36,6 +36,8 @@ subroutine initialize(files, config)
character(max_path_len) :: filename
logical :: lexist
integer :: yyyymm
character(len=6) :: adate
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
......@@ -43,19 +45,30 @@ subroutine initialize(files, config)
integer :: numx, numy
integer :: i, ierr
! read list of receptors
! ----------------------
filename = trim(files%file_recept)
call read_reclist(filename)
! initialize flexpart variables
! -----------------------------
if ( .not.config%lnested ) filename = trim(files%path_flexpart)//'header'
if ( config%lnested ) filename = trim(files%path_flexpart)//'header_nest'
print*, filename
inquire(file=trim(filename),exist=lexist)
yyyymm=config%datei/100
write(adate,'(i6)') yyyymm
lexist = .false.
i = 1
do while ( (.not.lexist).and.(i.le.nrec) )
if ( .not.config%lnested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
if ( config%lnested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
inquire(file=trim(filename),exist=lexist)
i = i + 1
end do
if ( .not.lexist ) then
print*, 'ERROR initialize: cannot find flexpart header'
write(*,*) 'ERROR initialize: cannot find flexpart header'
stop
endif
print*, 'Reading flexpart header: '//trim(filename)
write(*,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
nxgrid = numx
......
......@@ -18,6 +18,7 @@ SRCS = mod_var.f90 \
initialize.f90 \
grid_convert.f90 \
write_ncdf.f90 \
read_reclist.f90 \
grid_to_ncdf.f90
OBJECTS = $(SRCS:.f90=.o)
......
......@@ -36,6 +36,7 @@ module mod_settings
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) :: file_recept ! receptor list file
end type files_t
......@@ -43,6 +44,8 @@ module mod_settings
type :: config_t
integer :: datei ! start date (yyyymmdd)
integer :: datef ! end date (yyyymmdd)
logical :: lnested ! use nested flexpart output (true or false)
logical :: linversionout ! inversion formatted flexpart output (true or false)
......@@ -274,8 +277,17 @@ module mod_settings
identifier = "path_output:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%path_output = cc
identifier = "file_recept:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%file_recept = cc
! read config settings
identifier = "datei:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%datei = cn
identifier = "datef:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%datef = cn
identifier = "lnested:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lnested = cl
......
......@@ -35,6 +35,7 @@ module mod_var
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 :: recname_len=7 ! length of receptor names
! flexpart variables
! ------------------
......@@ -58,5 +59,7 @@ module mod_var
real(kind=8), dimension(:), allocatable :: reltime ! release times
integer :: datei ! start date of run
integer :: timei ! start time of run
integer :: nrec ! number of receptors
character(recname_len), dimension(:), allocatable :: recname ! names of receptors
end module mod_var
!---------------------------------------------------------------------------------------
! GRID_TO_NCDF: read_reclist
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> read_reclist
!!
!---------------------------------------------------------------------------------------
subroutine read_reclist(filename)
use mod_var
implicit none
character(len=max_path_len), intent(in) :: filename
character(len=200) :: line
integer :: ierr
integer :: cnt
! count number of receptors
open(100,file=trim(filename),action='read',status='old',iostat=ierr)
if(ierr.gt.0) then
write(*,*) 'ERROR: cannot open: '//trim(filename)
stop
endif
write(*,*) 'Reading receptors file: '//trim(filename)
cnt = 0
do while (ierr.eq.0)
read(100,*,iostat=ierr,end=10) line
cnt = cnt + 1
end do
10 continue
close(100)
nrec = cnt
write(*,*) 'Number of receptors: ',nrec
! read receptors
allocate ( recname(nrec) )
open(100,file=trim(filename),action='read',status='old',iostat=ierr)
write(*,*) 'Receptors: '
do cnt = 1, nrec
read(100,*) recname(cnt)
write(*,*) recname(cnt)
end do
close(100)
end subroutine read_reclist
......@@ -20,7 +20,7 @@ path_ohfield: /mypath_to_ohfields/
path_obs: /mypath/TEST_INPUT/OBS/CO2/
# path where to write observation output
path_obsout: /mypath/TEST_OUTPUT/OBS/CO2/
# observation file format (one of obspack, wdcgg, noaa, basic)
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
obsformat: obspack
# observation file suffix
suffix: .txt
......
......@@ -20,7 +20,7 @@ path_ohfield: /mypath_to_ohfields/
path_obs: /mypath/TEST_INPUT/OBS/CO2/
# path where to write observation output
path_obsout: /mypath/TEST_OUTPUT/OBS/CO2/
# observation file format (one of obspack, wdcgg, noaa, basic)
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
obsformat: obspack
# observation file suffix
suffix: .txt
......
......@@ -20,7 +20,7 @@ path_ohfield: /mypath_to_ohfields/
path_obs: /mypath/TEST_INPUT/OBS/GHG/
# path where to write observation output
path_obsout: /mypath/TEST_OUTPUT/OBS/GHG/
# observation file format (one of obspack, wdcgg, noaa, basic)
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
obsformat: wdcgg
# observation file suffix
suffix: .dat
......
......@@ -20,7 +20,7 @@ path_ohfield: /mypath_to_ohfields/
path_obs: /mypath/TEST_INPUT/OBS/GHG/
# path where to write observation output
path_obsout: /mypath/TEST_OUTPUT/OBS/GHG/
# observation file format (one of obspack, wdcgg, noaa, basic)
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
obsformat: noaa
# observation file suffix
suffix: event.txt
......
#!/bin/bash
#---------------------------------------------------
file='./SETTINGS_eSTICC_Aki_data.txt'
#file='./SETTINGS_eSTICC_NILU_data.txt'
#file='./SETTINGS_eSTICC_Aki_data.txt'
file='./SETTINGS_eSTICC_NILU_data.txt'
./prep_flexpart ${file}
......@@ -298,6 +298,8 @@ program main
call read_obs_esticc(settings, jd, nr, nobs, size(obs_in,1),obs_in, &
obs_name_short, obs_date, obs_lat, obs_lon, obs_alt, obs, obs_mdm, firstNetCDF)
endif
elseif ( settings%obsformat.eq.'icos' ) then
call read_icos(settings, jd, nr, nobs, obs)
elseif ( settings%obsformat.eq.'basic' ) then
call read_basic(settings, jd, nr, nobs, obs)
else
......
......@@ -34,6 +34,7 @@ SRCS = mod_var.f90 \
read_obs_esticc.f90 \
read_nilu.f90 \
read_basic.f90 \
read_icos.f90 \
main.f90
......
!---------------------------------------------------------------------------------------
! PREP_FLEXPART: read_icos
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> read_icos
!!
!! Purpose: Reads observations from ICOS files and averages and/or filters them.
!! Writes the averaged/filtered observations to file in the format
!! required by FLEXINVERT.
!!
!! Interface:
!!
!! Inputs
!! settings - settings data structure
!! jd - julian day for start of month
!! nr - index to the receptors list
!!
!! Outputs
!! nobs - number of observations
!! obs - observations data structure
!!
!! Externals
!! caldate
!! alloc_obs
!! parse_string
!! sort
!!
!---------------------------------------------------------------------------------------
subroutine read_icos(settings, jd, nr, nobs, obs)
use mod_settings
use mod_var
use mod_dates
use mod_obs
use mod_strings
use mod_tools
implicit none
type (settings_t), intent(in) :: settings
real(kind=8), intent(in) :: jd
integer, intent(in) :: nr
integer, intent(out) :: nobs
type (obs_t), intent(out) :: obs
character(len=max_path_len) :: file_out
character(len=max_path_len) :: file_obs
character(len=200) :: line
character(len=200), dimension(20) :: args
character(len=1) :: flag
real, dimension(16) :: temp
logical :: lexist
integer :: ierr
integer :: cnt
integer :: narg
integer :: i, l, m, n, p, q, r
integer :: jjjjmmdd, hhmiss, yyyy, mm, dd, hh, mi
integer :: eomday, hloc
real :: conc, err, freq, lat, lon, alt, hgt
real(kind=8) :: jdate
real(kind=8), dimension(maxobs) :: jdobs
real, dimension(maxobs) :: latobs, lonobs, altobs, concobs, errobs
call caldate(jd, jjjjmmdd, hhmiss)
eomday = calceomday(jjjjmmdd/100)
freq = 1.
! find file matching receptor name
do n = 1, nfiles
do i = 1, len_trim(filelist(n))-2
if ( filelist(n)(i:(i+2)).eq.trim(recname(nr)) ) go to 10
end do
end do
10 continue
if ( n.gt.nfiles ) then
print*, 'WARNING: file not found for receptor '//trim(recname(nr))
go to 20
endif
file_obs = filelist(n)
! open obs file
open(100,file=trim(settings%path_obs)//trim(file_obs),action='read',status='old',iostat=ierr)
if ( ierr.ne.0 ) then
print*, 'ERROR: cannot open: '//trim(settings%path_obs)//trim(file_obs)
stop
endif
print*, 'Reading file: '//trim(settings%path_obs)//trim(file_obs)
! read header
l = len_trim('# LATITUDE:')
m = len_trim('# LONGITUDE:')
p = len_trim('# SAMPLING HEIGHTS:')
r = len_trim('# ALTITUDE:')
q = len_trim('#Site')
do while (ierr.eq.0)
read (100, fmt='(A)', iostat=ierr) line
if ( line(1:q) == '#Site' ) ierr = 1
print*, 'line(1:l): ',line(1:l)
if ( line(1:l) == '# LATITUDE:' ) then
call parse_string(line, " ", args(:), narg)
read(args(3),*) lat
endif
if ( line(1:m) == '# LONGITUDE:' ) then
call parse_string(line, " ", args(:), narg)
read(args(3),*) lon
endif
! height metres above ground
if ( line(1:p) == '# SAMPLING HEIGHTS:' ) then
call parse_string(line, " ", args(:), narg)
read(args(4),*) hgt
endif
! height metres above sea level
if ( line(1:r) == '# ALTITUDE:' ) then
call parse_string(line, " ", args(:), narg)
read(args(3),*) alt
endif
end do
if ( settings%zref.eq.1 ) then
! metres above ground
alt = hgt
else if ( settings%zref.eq.2 ) then
! metres above sea level
alt = alt + hgt
endif
print*, 'lat, lon, alt: ',lat, lon, alt
! read data
cnt = 0
read_loop: do
read(100,fmt='(A)',iostat=ierr) line
if(ierr.gt.0) exit read_loop
call parse_string(line, ";", args(:), narg)
do n = 2, 11
read(args(n),*) temp(n)
end do
print*, temp
yyyy = int(temp(3))
mm = int(temp(4))
dd = int(temp(5))
hh = int(temp(6))
mi = int(temp(7))
conc = temp(9)
err = temp(10)
flag = args(12)
print*, flag
! calculate julian date
jjjjmmdd = yyyy*10000+mm*100+dd
hhmiss = hh*10000+mi*100
jdate = juldate(jjjjmmdd, hhmiss)
if ( flag.eq.'N' ) cycle read_loop
if ( conc.le.-999. ) cycle read_loop
if ( (jdate.lt.jd).or.(jdate.ge.min(jreldatef,(jd+real(eomday,kind=8)))) ) cycle read_loop
cnt = cnt + 1
jdobs(cnt) = jdate
concobs(cnt) = conc
errobs(cnt) = err
if ( errobs(cnt).le.-9.99 ) errobs(cnt) = 0.
latobs(cnt) = lat
lonobs(cnt) = lon
altobs(cnt) = alt
end do read_loop
nobs = cnt
print*, 'lat, lon, alt = ',latobs(1),lonobs(1),altobs(1)
! close input file
close(100)
print*, 'Number of raw observations: ',nobs
! no observations go to end
if ( nobs.eq.0 ) then
print*, 'WARNING: no data for '//trim(recname(nr))
go to 20
endif
! allocate obs
call alloc_obs(nobs, obs)
! average and/or select observations
call process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, concobs, errobs, obs)
! no observations after selection
if ( nobs.eq.0 ) go to 20
! write formatted obs to file
file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
inquire(file=trim(file_out),exist=lexist)
if( lexist ) then
! append to existing
print*, 'WARNING: appending to file '//trim(file_out)
open(100,file=trim(file_out),status='old',action='write',access='append',iostat=ierr)
do i = 1, nobs
call caldate(obs%jdi(i), jjjjmmdd, hhmiss)
write(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)') &
trim(recname(nr)), jjjjmmdd, hhmiss, obs%jdi(i), obs%con(i), obs%err(i), obs%num(i)
end do
close(100)
else
! create new file
open(100,file=trim(file_out),status='new',action='write',iostat=ierr)
write(100,fmt='(A)') 'rec yyyymmdd hhmmss juldate conc err num'
do i = 1, nobs
call caldate(obs%jdi(i), jjjjmmdd, hhmiss)
write(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)') &
trim(recname(nr)), jjjjmmdd, hhmiss, obs%jdi(i), obs%con(i), obs%err(i), obs%num(i)
end do
close(100)
endif
! write detailed station list file
file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
inquire(file=trim(file_out),exist=lexist)
if( lexist ) then
! append to existing
open(100,file=trim(file_out),status='old',action='write',access='append',iostat=ierr)
write(100,fmt='(A3,1X,F6.2,1X,F6.2,1X,F6.1)') trim(recname(nr)), obs%lat(1), obs%lon(1), obs%alt(1)
close(100)
else
! create new file
open(100,file=trim(file_out),status='new',action='write',iostat=ierr)
write(100,fmt='(A3,1X,F6.2,1X,F6.2,1X,F6.1)') trim(recname(nr)), obs%lat(1), obs%lon(1), obs%alt(1)
close(100)
endif
20 continue
end subroutine read_icos
......@@ -24,15 +24,23 @@ datei: 20120101
datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: congrad
method: m1qn3
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false) (ghg species only)
# Include ocean boxes in inversion (true or false)
# currently only for ghg species
inc_ocean: .false.
# Optimize initial mixing ratios (true or false)
opt_cini: .true.
# Use spatial correlation in error covariance matrix (true or false)
# if use regions based ecosystems then should be false
spa_corr: .true.
......@@ -42,10 +50,12 @@ spa_corr: .true.
# 1 = use best guess (file must be specified in SETTINGS_files)
prior_bg: 0
# Restart a crashed run from last iteration (congrad only)
# Restart a crashed run
# for congrad/m1qn3 will pick-up from last iteration
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
restart: 0
restart: 1
# Verbose output
# only use for debugging small runs
......@@ -79,7 +89,7 @@ xres: 1.0
yres: 1.0
# Spatial aggregation of grid (true or false)
regions: .true.
regions: .false.
# State vector time resolution:
# time resolution at which NEE fluxes are optimized
......@@ -121,6 +131,5 @@ sigma_ocean: 1000.
sigmatime: 7.
# Total error for domain (Tg/y)
# if globerr <= 0 prior error covariance matrix not scaled
globerr: 300.
......@@ -24,15 +24,23 @@ datei: 20120101
datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: congrad
method: analytic
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false) (ghg species only)
# Include ocean boxes in inversion (true or false)
# currently only for ghg species
inc_ocean: .false.
# Optimize initial mixing ratios (true or false)
opt_cini: .true.
# Use spatial correlation in error covariance matrix (true or false)
# if use regions based ecosystems then should be false
spa_corr: .true.
......@@ -42,14 +50,16 @@ spa_corr: .true.
# 1 = use best guess (file must be specified in SETTINGS_files)
prior_bg: 0
# Restart a crashed run from last iteration (congrad only)
# Restart a crashed run
# for congrad/m1qn3 will pick-up from last iteration
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
restart: 0
restart: 1
# Verbose output
# only use for debugging small runs
verbose: .true.
verbose: .false.
# Species ("co2" or "ghg")
spec: co2
......@@ -121,6 +131,6 @@ sigma_ocean: 1000.
sigmatime: 7.
# Total error for domain (Tg/y)
# if globerr <= 0 prior error covariance matrix not scaled
# globerr <= 0: prior error covariance matrix not scaled
globerr: 300.
......@@ -24,14 +24,21 @@ datei: 20120101
datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: congrad
method: analytic
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false) (ghg species only)
inc_ocean: .false.
# Include ocean boxes in the optimization (true or false)
inc_ocean: .true.
# Optimize initial mixing ratios (true or false)
opt_cini: .true.
# Use spatial correlation in error covariance matrix (true or false)
# if use regions based ecosystems then should be false
......@@ -42,7 +49,9 @@ spa_corr: .true.
# 1 = use prior best guess (file must be specified in SETTINGS_files)
prior_bg: 0
# Restart a crashed run from last iteration (congrad only)
# Restart a crashed run
# for congrad/m1qn3 will pick-up from last iteration
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
restart: 0
......@@ -54,8 +63,8 @@ verbose: .true.
# Species ("co2" or "ghg")
spec: ghg
# Molar mass (in flux files, e.g. CO2-C=12, CH4=16)
molarmass: 16.
# Molar mass in flux files (e.g. CO2-C=12, CH4=16)
molarmass: 16.
# Coefficient to convert from grid units of ppt to ppb
coeff: 1.e-3
......@@ -113,6 +122,5 @@ sigma_ocean: 1000.
sigmatime: 30.
# Total error for inversion domain (Tg/y)
# if globerr <= 0 prior error covariance matrix not scaled
globerr: 10.
......@@ -26,12 +26,19 @@ datef: 20120131
# Inversion method ('analytic', 'congrad' or 'm1qn3')
method: congrad
# Number of iterations
# only used if method is 'congrad' or 'm1qn3'
maxiter: 2
# Optimize flux offsets (true or false)
# optimizes the offsets rather than the fluxes themselves (ghg species only)
offsets: .false.
# Include ocean boxes in inversion (true or false) (ghg species only)
inc_ocean: .false.
# Include ocean boxes in the optimization (true or false)
inc_ocean: .true.
# Optimize initial mixing ratios (true or false)
opt_cini: .false.
# Use spatial correlation in error covariance matrix (true or false)
# if use regions based ecosystems then should be false
......@@ -42,10 +49,12 @@ spa_corr: .true.
# 1 = use prior best guess (file must be specified in SETTINGS_files)
prior_bg: 0
# Restart a crashed run from last iteration (congrad only)
# Restart a crashed run
# for congrad/m1qn3 will pick-up from last iteration
# for analytic will only use pre-calculated covariance matrix and boundary conditions
# 0 = new run
# 1 = restart crashed run
restart: 0
restart: 1
# Verbose output
# only use for debugging small runs
......@@ -61,7 +70,7 @@ molarmass: 16.
coeff: 1.e-3
# Use nested flexpart output (true or false)
nested: .nested.
nested: .true.
# Inversion domain:
# if using nested output it must match the inversion domain
......@@ -112,6 +121,7 @@ sigma_ocean: 1000.
# Temporal correlation length: unit (days)
sigmatime: 30.
# Total error for inversion domain (Tg/y), make sure this value is realistic if you change domain size or spezies!
globerr: 10.
# Total error for inversion domain (Tg/y)
# globerr <= 0: prior error covariance matrix not scaled
globerr: -1.
......@@ -125,10 +125,10 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
x0(:,1) = dble(states%px0)
! transform prior state vector chi space
z0 = phi2chi(states, covar, x0)
z0 = phi2chi(config, states, covar, x0)
! transform state vector chi space
zbg = phi2chi(states, covar, xbg)
zbg = phi2chi(config, states, covar, xbg)
! initial lanczos vector is initial gradient in chi space
zgrad0(:,1) = grad(:,1)
......@@ -170,7 +170,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
close(100)
else
! transform zw to physical space
zx = chi2phi(states, covar, zw)
zx = chi2phi(config, states, covar, zw)
states%px = real(zx(:,1),kind=4)
cost_o = 0.
grad_o(:) = 0.
......@@ -178,7 +178,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
! calculate cost
! Jp = (p - p0)^TB^(-1)(p - p0)
delx = states%px - states%px0
tmp = prod_icov(covar, delx)
tmp = prod_icov(config, covar, delx)
cost_p = dot_product( delx, real(tmp,kind=4) )
! J = 0.5*(Jp + Jo)
cost = 0.5*(cost_p + cost_o)
......@@ -188,7 +188,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
! transform gradient to chi space
x_ad(:,1) = dble(grad_o)
z_ad = phi2chi_ad(covar, x_ad )
z_ad = phi2chi_ad(config, covar, x_ad )
! new gradient J' = J'p + J'o
grad(:,1) = zw(:,1) + z_ad(:,1)
......@@ -317,7 +317,7 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
grad(:,1) = grad(:,1) - matmul(zcglwk(:,1:iter), zqg0(1:iter,1))
! transform state vector to physical space
xa = chi2phi(states, covar, za)
xa = chi2phi(config, states, covar, za)
states%px = real(xa(:,1),kind=4)
if (iter.eq.1) go to 10
......@@ -402,6 +402,9 @@ subroutine congrad(iter, grad, files, config, fluxes, obs, states, covar)
end do
end do
end do
if ( config%opt_cini ) then
tmp(npvar+1:nvar) = tmp(npvar+1:nvar) + config%cinierr * pevecs(jk,npvar+1:nvar)
endif
end do
! calculate last part
do i = 1, nvar
......
......@@ -50,14 +50,14 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
implicit none
logical, intent(in out) :: midpoints
logical, intent(in) :: midpoints
integer, intent(in) :: nx_in
integer, intent(in) :: ny_in
integer, intent(in) :: nx_out
integer, intent(in) :: ny_out
real, dimension(nx_in,ny_in,1), intent(in out) :: work_in
real, dimension(nx_in), intent(in out) :: lon_in
real, dimension(ny_in), intent(in out) :: lat_in
real, dimension(nx_in,ny_in,1), intent(in) :: work_in
real, dimension(nx_in), intent(in) :: lon_in
real, dimension(ny_in), intent(in) :: lat_in
real, dimension(nx_out), intent(in out) :: lon_out
real, dimension(ny_out), intent(in out) :: lat_out
real, dimension(nx_out,ny_out,1), intent(in out) :: work_out
......@@ -79,11 +79,12 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
real :: lly_out
real(kind=8), dimension(nx_in) :: lono
real(kind=8), dimension(ny_in) :: lato
real, dimension(nx_in) :: wlon_in
real, dimension(ny_in) :: slat_in
real, dimension(nx_in) :: lon, wlon_in
real, dimension(ny_in) :: lat, slat_in
real, dimension(nx_out) :: wlon_out
real, dimension(ny_out) :: slat_out
real, dimension(nx_in,ny_in) :: area_xy
real, dimension(nx_in,ny_in,1) :: work
integer, dimension(:), allocatable :: ind
ntime_in = 1
......@@ -92,6 +93,9 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
dy_out = dy
llx_out = llx
lly_out = lly
lon = lon_in
lat = lat_in
work = work_in
tot_in = 0.
tot_out = 0.
......@@ -114,66 +118,65 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
endif
! check lon and lat dimensions
lato = dble(lat_in)
lato = dble(lat)
allocate( ind(ny_in) )
call sort(ny_in, lato, ind)
deallocate(ind)
! indices to lat
do jy = 1, ny_in
indy(jy) = minloc(sngl(lato),dim=1,mask=sngl(lato).eq.lat_in(jy))
indy(jy) = minloc(sngl(lato),dim=1,mask=sngl(lato).eq.lat(jy))
end do
lat_in = sngl(lato)
where (lon_in.gt.180.) lon_in = lon_in - 360.
lono = dble(lon_in)
lat = sngl(lato)
where (lon.gt.180.) lon = lon - 360.
lono = dble(lon)
allocate( ind(nx_in) )
call sort(nx_in, lono, ind)
deallocate(ind)
! indices to lon
do ix = 1, nx_in
indx(ix) = minloc(sngl(lono),dim=1,mask=sngl(lono).eq.lon_in(ix))
indx(ix) = minloc(sngl(lono),dim=1,mask=sngl(lono).eq.lon(ix))
end do
lon_in = sngl(lono)
lon = sngl(lono)
! reorder flx
work_in(:,:,:) = work_in(indx,indy,:)
work(:,:,:) = work(indx,indy,:)
! input resolution
do ix = 2, nx_in
dx_in(ix) = abs(lon_in(ix)-lon_in(ix-1))
dx_in(ix) = abs(lon(ix)-lon(ix-1))
end do
dx_in(1) = dx_in(2)
do jy = 2, ny_in
dy_in(jy) = abs(lat_in(jy)-lat_in(jy-1))
dy_in(jy) = abs(lat(jy)-lat(jy-1))
end do
dy_in(1) = dy_in(2)
if(midpoints) then
! western edges
do ix = 1, nx_in
wlon_in(ix) = lon_in(ix)-0.5*dx_in(ix)
wlon_in(ix) = lon(ix)-0.5*dx_in(ix)
end do
! southern edges
do jy = 1, ny_in
slat_in(jy) = lat_in(jy)-0.5*dy_in(jy)
slat_in(jy) = lat(jy)-0.5*dy_in(jy)
end do
else
wlon_in = lon_in
wlon_in = lon
! longitude midpoints
do ix = 1, nx_in
lon_in(ix) = lon_in(ix)+0.5*dx_in(ix)
lon(ix) = lon(ix)+0.5*dx_in(ix)
end do
slat_in = lat_in
slat_in = lat
! latitude midpoints
do jy = 1, ny_in
lat_in(jy)=lat_in(jy)+0.5*dy_in(jy)
lat(jy)=lat(jy)+0.5*dy_in(jy)
end do
midpoints = .true.
endif
if( (nx_in.eq.nx_out) .and. (ny_in.eq.ny_out) ) then
write(logid,*) 'WARNING convertgrid: input and output grids are the same'
wlon_out = wlon_in
slat_out = slat_in
work_out = work_in
work_out = work
go to 20
endif
......@@ -203,11 +206,7 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
yolap2 = min(slat_in(jj)+dy_in(jj),slat_out(jy)+dy_out)
! calculate overlapping area
aolap = grid_area((yolap1+yolap2)/2., abs(yolap2-yolap1), abs(xolap2-xolap1))
mass = mass + aolap*work_in(ii,jj,nt)
!cgz debug
! if (jy.gt.ny_out-4) then
! write(logid, *) jy, ix, jj, ii, mass, work_in(ii,jj,nt), xolap1 , xolap2, yolap1 , yolap2
! endif
mass = mass + aolap*work(ii,jj,nt)
aolaptot = aolaptot + aolap
end do
end do
......@@ -227,7 +226,7 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
yolap2 = min(slat_in(jj)+dy_in(jj),slat_out(jy)+dy_out)
! calculate overlapping area
aolap = grid_area((yolap1+yolap2)/2., abs(yolap2-yolap1), abs(xolap2-xolap1))
mass = mass + aolap*work_in(ix2,jj,nt)
mass = mass + aolap*work(ix2,jj,nt)
aolaptot = aolaptot + aolap
end do
if(aolaptot.ne.0) then
......@@ -246,7 +245,7 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
xolap2 = min(wlon_in(ii)+dx_in(ii),wlon_out(ix)+dx_out)
! calculate overlapping area
aolap = grid_area((yolap1+yolap2)/2., abs(yolap2-yolap1), abs(xolap2-xolap1))
mass = mass + aolap*work_in(ii,jy2,nt)
mass = mass + aolap*work(ii,jy2,nt)
aolaptot = aolaptot + aolap
end do
if(aolaptot.ne.0) then
......@@ -273,14 +272,14 @@ subroutine convertgrid(midpoints, nx_in, ny_in, lon_in, lat_in, work_in, &
area_xy(:,:) = 0.
do jy = jy1, jy2
do ix = ix1, ix2
area = grid_area(lat_in(jy), dy_in(jy), dx_in(ix))
area = grid_area(lat(jy), dy_in(jy), dx_in(ix))
area_xy(ix,jy) = area
end do
end do
! calculate total in
do nt = 1, ntime_in
tot_in = tot_in + sum( sum(work_in(:,:,nt)*area_xy(:,:),dim=1) )
tot_in = tot_in + sum( sum(work(:,:,nt)*area_xy(:,:),dim=1) )
end do
write(logid,fmt='(1x,''Total mass in:'',1x,E12.4,1x,''Tg/y'')') tot_in*24.*365./float(ntime_in)/1.e9
......
......@@ -26,18 +26,20 @@
!! Interface:
!!
!! Inputs
!! files - file data structure
!! config - configuration settings data structure
!! corr - correlation matrix (zero matrix on input)
!!
!---------------------------------------------------------------------------------------
subroutine correl(config, corr)
subroutine correl(files, config, corr)
use mod_var
use mod_settings
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
real, dimension(nbox,nbox), intent(in out) :: corr
......@@ -47,6 +49,8 @@ subroutine correl(config, corr)
real :: pih, cor, dd, ddx, num
real :: lat1, lat2, lon1, lon2
real :: sigma1, sigma2
character(len=max_path_len) :: rowfmt
integer :: ierr
pih = pi/180.
......@@ -114,6 +118,14 @@ subroutine correl(config, corr)
end do ! ind2
end do ! ind1
! save correlations to file
open(100,file=trim(files%path_output)//'correl.txt',action='write',status='replace',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
print*, rowfmt
do ix = 1, nbox
write(100,fmt=rowfmt) corr(ix,:)
end do
close(100)
end subroutine correl
......@@ -241,10 +241,18 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
end do
end do
do i = 1, nvar
! uncertainty for flux variables
do i = 1, npvar
states%pxerr0(i) = real(sqrt(var(i)),kind=4)
end do
! uncertainty for scalars of initial mixing ratio
if ( config%opt_cini ) then
do i = 1, ntcini*ncini
states%pxerr0(npvar+i) = config%cinierr
end do
endif
! save output
! -----------
......
......@@ -56,11 +56,11 @@ subroutine init_cini(files, obs)
logical :: lexist
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
integer :: ix, jy, kz, i
integer :: ix, jy, kz, i, n
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs) :: cini
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
......@@ -73,14 +73,18 @@ subroutine init_cini(files, obs)
concini(:,:,:) = 0.
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:) = 0.
cini(:,:) = 0.
! sort obs chronologically for efficiency
jdates = obs%jdate(:)
call sort(nobs,jdates,ind)
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: jdates = ',jdates
! print*, 'init_cini: obs%jdate = ',obs%jdate
!print*, 'before sort, init_cini: ind = ',ind
!print*, 'before sort, init_cini: jdates = ',jdates
!print*, 'before sort, init_cini: obs%jdate = ',obs%jdate
call sort(nobs,jdates,ind)
!print*, 'after sort, init_cini: ind = ',ind
!print*, 'after sort,init_cini: jdates = ',jdates
!print*, 'after sort, init_cini: obs%jdate = ',obs%jdate
call flush()
do i = 1, nobs
......@@ -90,8 +94,8 @@ subroutine init_cini(files, obs)
write(areldate,'(I8.8)') jjjjmmdd
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'
! print*, 'init_cini: filename = ',filename
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
! read flexpart initial concentrations weighting
inquire(file=trim(filename),exist=lexist)
......@@ -225,6 +229,12 @@ subroutine init_cini(files, obs)
! 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
! print*, 'init_cini: gbl_lat(jy), n = ',gbl_lat(jy), n
do ix = 1, nxgrid
! interpolate between concentrations
if( (jdates(i) - trajdays).lt.jdmid ) then
......@@ -234,7 +244,7 @@ subroutine init_cini(files, obs)
conc = concini(ix,jy,kz) + (proconcini(ix,jy,kz) - concini(ix,jy,kz)) * &
real((jdates(i) - trajdays) - jdmid, kind=4)/30.
endif
cini(i) = cini(i) + gridini(ix,jy,kz)*conc
cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*conc
end do
end do
!cgz: check cini
......@@ -249,7 +259,7 @@ subroutine init_cini(files, obs)
lastjjjjmm = jjjjmm
write(logid,*) 'CGZ: finished obs: nr=', i, ', obs=', obs%conc(i), ', cini=', obs%cini(i)
write(logid,*) 'CGZ: finished obs: nr=', i, ', obs=', obs%conc(ind(i)), ', cini=', cini(i)
end do
......@@ -257,7 +267,7 @@ subroutine init_cini(files, obs)
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: cini = ',cini
do i = 1, nobs
obs%cini(ind(i)) = cini(i)
obs%cini(ind(i),:) = cini(i,:)
end do
! print*, 'init_cini: obs%cini = ',obs%cini
......
......@@ -60,6 +60,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
use mod_states
use mod_covar
use mod_perturb
use mod_save, only : save_for_restart
implicit none
......@@ -92,12 +93,17 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
integer :: ix, jy
integer :: n, nb, nd, nt, i
integer :: ierr
real(kind=8), dimension(:), allocatable :: jdate
character(len=3), dimension(:), allocatable :: recs
real :: conc
character(len=20) :: dump
! global prior fluxes
! -------------------
! global fluxes are used for calculating contribution to mixing ratios
! from outside the inversion domain and have monthly resolution
! from outside the inversion domain and for the contribution from
! ocean and biomass burning fluxes inside the domain
call alloc_fluxes(config, fluxes)
fluxes%flxnee(:,:,:) = 0.
......@@ -154,6 +160,9 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
nread = eomday*24
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
print*, 'filename_ff = ',files%filename_ff
print*, 'adate = ',adate
print*, 'amonth = ',amonth
filename = str_replace(files%filename_ff, 'YYYY', adate)
filename = str_replace(filename, 'MM', amonth)
filename = trim(files%path_prior)//trim(filename)
......@@ -546,19 +555,34 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! ------------------
! number of state variables
nvar = ndt*ntstate*nbox
npvar = ndt*ntstate*nbox
ndvar = ndt*nbox
if ( config%opt_cini ) then
nvar = npvar + ntcini*ncini
else
nvar = npvar
endif
write(logid,*) 'Total number of state variables: ',nvar
write(logid,*) 'Number of flux state variables: ',npvar
write(logid,*) 'Number of state variables per flux time step: ',ndvar
write(logid,*) 'Total number of state variables: ',nvar
call alloc_states(states)
allocate( xerr(ndvar) )
! assign prior state vector
states%px0(:) = 0.
if ( config%opt_cini ) states%px0(npvar+1:nvar) = 1.
! assign time stamps to flux time steps
states%xtime(:) = statime
! assign time stamps to mixing ratio scalar time steps
do n = 1, ntcini
states%cinitime(n) = juldatei + dble((n-1)*cinires)
end do
print*, 'init_ghg: cinitime = ',states%cinitime
! assign best guess of state vector
if ( (config%prior_bg.eq.1).and.(config%method.eq.'congrad') ) then
! from previous inversion