Commit 859ab8e0 authored by ronesy's avatar ronesy

Bug fix for m1qn3_interface and update to avoid recalculating initial conditions on restart

parent b895297b
......@@ -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,7 +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)
mass = mass + aolap*work(ii,jj,nt)
aolaptot = aolaptot + aolap
end do
end do
......@@ -223,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
......@@ -242,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
......@@ -269,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
......
......@@ -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)
......@@ -634,9 +643,12 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! initial concentrations
! ----------------------
! initial concentrations from termination of trajectories
call init_cini(files, obs)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
call init_cini(files, obs)
endif
! add random perturbation for MC
! ------------------------------
......@@ -648,6 +660,39 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
call perturb_obs(config, obs)
endif
! save datatype obs for restart
! -----------------------------
if ( config%restart.eq.1 ) then
! read datatype obs from file
write(logid,*) 'WARNING init_co2: reading cini from previous run'
open(100,file=trim(files%path_output)//'obsdatatype.txt',status='old',action='read',iostat=ierr)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR init_co2: cannot open file obsdatatype.txt'
stop
endif
allocate( recs(nobs) )
allocate( jdate(nobs) )
write(rowfmt,'(A,I6,A)') '(A3,1X,F14.6,1X,F10.4,1X,',ncini,'(F10.4,1X))'
read(100,*) dump
do n = 1, nobs
read(100,rowfmt,iostat=ierr) recs(n), jdate(n), conc, obs%cini(n,:)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR init_co2: reading obsdatatype.txt'
exit
endif
end do
! check consistency
if ( recs(1).ne.obs%recs(1).or.recs(nobs).ne.obs%recs(nobs).or.jdate(1).ne.obs%jdate(1).or.jdate(nobs).ne.obs%jdate(nobs) ) then
write(logid,*) 'ERROR init_co2: data in obsdatatype.txt inconsistent with run'
stop
endif
deallocate( recs )
deallocate( jdate )
else
! save datatype obs
call save_for_restart(files, obs)
endif
end subroutine init_co2
......
......@@ -60,31 +60,36 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
use mod_states
use mod_covar
use mod_perturb
use mod_save, only : save_for_restart
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (fluxes_t), intent (in out) :: fluxes
type (obs_t), intent (in out) :: obs
type (states_t), intent (in out) :: states
type (covar_t), intent (in out) :: covar
character(max_path_len) :: filename
character(max_path_len) :: rowfmt
character(len=4) :: adate
logical :: lexist
real, dimension(:,:,:), allocatable :: flx
real, dimension(nbox,ntstate) :: flx_box
real, dimension(nbox,ntstate) :: err_box
real, dimension(nbox) :: xerr
real, dimension(nbox,nbox) :: corr
integer :: yyyymmdd, hhmiss
real :: area
real(kind=8) :: jd
integer :: ix, jy, ix0, jy0
integer :: n, nb, nread, num
integer :: ierr
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (fluxes_t), intent (in out) :: fluxes
type (obs_t), intent (in out) :: obs
type (states_t), intent (in out) :: states
type (covar_t), intent (in out) :: covar
character(max_path_len) :: filename
character(max_path_len) :: rowfmt
character(len=4) :: adate
logical :: lexist
real, dimension(:,:,:), allocatable :: flx
real, dimension(nbox,ntstate) :: flx_box
real, dimension(nbox,ntstate) :: err_box
real, dimension(nbox) :: xerr
real, dimension(nbox,nbox) :: corr
integer :: yyyymmdd, hhmiss
real :: area
real(kind=8) :: jd
integer :: ix, jy, ix0, jy0
integer :: n, nb, nread, num
integer :: ierr
real(kind=8), dimension(:), allocatable :: jdate
character(len=3), dimension(:), allocatable:: recs
real :: conc
character(len=20) :: dump
! global prior fluxes
! -------------------
......@@ -369,7 +374,10 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! ----------------------
! initial concentrations from termination of trajectories
call init_cini(files, obs)
! only do for a fresh run
if ( config%restart.eq.0 ) then
call init_cini(files, obs)
endif
! add random perturbation for MC
! ------------------------------
......@@ -382,6 +390,39 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
call perturb_obs(config, obs)
endif
! save datatype obs for restart
! -----------------------------
if ( config%restart.eq.1 ) then
! read datatype obs from file
write(logid,*) 'WARNING init_ghg: reading cini from previous run'
open(100,file=trim(files%path_output)//'obsdatatype.txt',status='old',action='read',iostat=ierr)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR init_ghg: cannot open file obsdatatype.txt'
stop
endif
allocate( recs(nobs) )
allocate( jdate(nobs) )
write(rowfmt,'(A,I6,A)') '(A3,1X,F14.6,1X,F10.4,1X,',ncini,'(F10.4,1X))'
read(100,*) dump
do n = 1, nobs
read(100,rowfmt,iostat=ierr) recs(n), jdate(n), conc, obs%cini(n,:)
if ( ierr.ne.0 ) then
write(logid,*) 'ERROR init_ghg: reading obsdatatype.txt'
exit
endif
end do
! check consistency
if ( recs(1).ne.obs%recs(1).or.recs(nobs).ne.obs%recs(nobs).or.jdate(1).ne.obs%jdate(1).or.jdate(nobs).ne.obs%jdate(nobs) ) then
write(logid,*) 'ERROR init_ghg: data in obsdatatype.txt inconsistent with run'
stop
endif
deallocate( recs )
deallocate( jdate )
else
! save datatype obs
call save_for_restart(files, obs)
endif
end subroutine init_ghg
......
......@@ -172,6 +172,9 @@ subroutine m1qn3_interface(iter, grad, files, config, fluxes, obs, states, covar
close(100)
write(logid,*) 'Cost read for iteration '//aiter//': ',cost
f = cost
! transform x to physical space (only needed if lastiter = maxiter)
px = chi2phi(config, states, covar, x)
states%px = real(px(:,1),kind=4)
else
! transform x to physical space
px = chi2phi(config, states, covar, x)
......
......@@ -161,6 +161,7 @@ module mod_obs
obs%delta(:) = 0.
obs%bkg(:) = 0.
obs%cini(:,:) = 0.
obs%cinipos(:) = 0.
obs%measerr(:) = 0.
obs%err(:) = 0.
obs%ocn(:) = 0.
......@@ -189,6 +190,7 @@ module mod_obs
if ( allocated(obs%delta) ) deallocate(obs%delta)
if ( allocated(obs%bkg) ) deallocate(obs%bkg)
if ( allocated(obs%cini) ) deallocate(obs%cini)
if ( allocated(obs%cinipos) ) deallocate(obs%cinipos)
if ( allocated(obs%measerr) ) deallocate(obs%measerr)
if ( allocated(obs%err) ) deallocate(obs%err)
if ( allocated(obs%ocn) ) deallocate(obs%ocn)
......
......@@ -35,7 +35,7 @@ module mod_save
implicit none
private
public :: save_state, save_co2, save_ghg, save_obs, save_nee_hourly, save_scalar_cini, check
public :: save_state, save_co2, save_ghg, save_obs, save_nee_hourly, save_scalar_cini, save_for_restart, check
contains
......@@ -716,6 +716,38 @@ module mod_save
end subroutine save_scalar_cini
! --------------------------------------------------
! save_for_restart
! --------------------------------------------------
!> save_for_restart
!!
!! Purpose: Saves the datatype obs for
!! eventual restart
!!
! --------------------------------------------------
subroutine save_for_restart(files, obs)
use mod_settings
use mod_obs
use mod_var
type (files_t), intent(in) :: files
type (obs_t), intent(in) :: obs
character(len=max_path_len) :: rowfmt
integer :: i, ierr
! write to file
open(100,file=trim(files%path_output)//'obsdatatype.txt',status='replace',action='write',iostat=ierr)
write(rowfmt,'(A,I6,A)') '(A3,1X,F14.6,1X,F10.4,1X,',ncini,'(F10.4,1X))'
write(100,*) 'recs jdate conc cini'
do i = 1, nobs
write(100,rowfmt) obs%recs(i), obs%jdate(i), obs%conc(i), obs%cini(i,:)
end do
close(100)
end subroutine save_for_restart
! --------------------------------------------------
! check
! --------------------------------------------------
......
Markdown is supported
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