Commit 7ba83dd0 authored by ronesy's avatar ronesy
Browse files

Corrections to date/time for using averaging of footprints

parent 9fe5d46a
......@@ -121,7 +121,7 @@ subroutine average_fp(files, config, obs)
! Find the first release number corresponding to the current observation
release_nr=1
do while ( int(releases(release_nr,1)*1e5,kind=8).lt.int(obs%jdate(i)*1e5,kind=8) .and. release_nr.lt.numpoint )
do while ( int(releases(release_nr,1)*1e4,kind=8).lt.int(obs%jdate(i)*1e4,kind=8) .and. release_nr.lt.numpoint )
release_nr = release_nr+1
end do
......@@ -162,15 +162,14 @@ subroutine average_fp(files, config, obs)
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
gridnest_tmp = gridnest ! Save the footprint to temporary grid for summation
endif
gridnest_tmp = gridnest ! Save the footprint to temporary grid for summation
! Current release_nr was already read, start with next
release_nr = release_nr+1
write(logid,*) 'average_fp: obs%jdate(i)+obs%avetime(i) = ',obs%jdate(i)+obs%avetime(i)
write(logid,*) 'average_fp: releases(release_nr,2) = ', releases(release_nr,2)
do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obs%jdate(i)+obs%avetime(i))*1e5,kind=8) .and. release_nr.lt.numpoint )
do while ( int(releases(release_nr,2)*1e4,kind=8).le.int((obs%jdate(i)+obs%avetime(i))*1e4,kind=8) .and. release_nr.lt.numpoint )
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
month = jjjjmmdd/100
......
......@@ -121,7 +121,7 @@ subroutine init_cini(files, config, obs)
endif
! find the first release number corresponding to the current observation
release_nr=1
do while ( int(releases(release_nr,1)*1.e5,kind=8).lt.int(jdates(i)*1.e5,kind=8) .and. release_nr.lt.numpoint )
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
......@@ -140,7 +140,7 @@ subroutine init_cini(files, config, obs)
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.e5,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1.e5,kind=8) .and. release_nr.lt.numpoint )
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
......
......@@ -52,9 +52,9 @@ subroutine init_fpctm(files, config, obs)
type (config_t), intent (in) :: config
type (obs_t), intent (in out) :: obs
character(max_path_len) :: path_flexrec, filename
character(max_path_len) :: path_flexrec, pathname, filename
character(max_name_len) :: varname
character(len=8) :: adate, areldate, areltime
character(len=8) :: adate, ayear, areldate, areltime
logical :: lexist
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
......@@ -100,7 +100,7 @@ subroutine init_fpctm(files, config, obs)
write(adate,'(I6.6)') jjjjmmdd/100
month = jjjjmmdd/100
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
print*, 'init_cini: pathflexrec = ',path_flexrec
print*, 'init_fpctm: pathflexrec = ',path_flexrec
if ( config%average_fp ) then
......@@ -122,7 +122,7 @@ subroutine init_fpctm(files, config, obs)
endif
! find the first release number corresponding to the current observation
release_nr=1
do while ( int(releases(release_nr,1)*1e5,kind=8).lt.int(jdates(i)*1e5,kind=8) .and. release_nr.lt.numpoint )
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
......@@ -132,7 +132,7 @@ subroutine init_fpctm(files, config, obs)
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
print*, 'init_fpctm: filename = ',filename
! read first flexpart grid_init files
inquire ( file=trim(filename), exist=lexist )
if ( lexist ) then
......@@ -141,7 +141,7 @@ subroutine init_fpctm(files, config, obs)
nr_init = 1
! current release_nr was already read, start with next
release_nr = release_nr + 1
do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1e5,kind=8) .and. release_nr.lt.numpoint )
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
......@@ -172,7 +172,7 @@ subroutine init_fpctm(files, config, obs)
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
print*, 'init_fpctm: filename = ',filename
! read flexpart grid_init file
inquire(file=trim(filename),exist=lexist)
if(lexist) then
......@@ -221,8 +221,10 @@ subroutine init_fpctm(files, config, obs)
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(files%path_initconc)//trim(filename)
filename = trim(pathname)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
......@@ -235,8 +237,10 @@ subroutine init_fpctm(files, config, obs)
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(files%path_initconc)//trim(filename)
filename = trim(pathname)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
......@@ -247,8 +251,10 @@ subroutine init_fpctm(files, config, obs)
endif
! get previous concentration field
write(adate,'(I6.6)') prejjjjmm
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
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
......@@ -259,8 +265,10 @@ subroutine init_fpctm(files, config, obs)
endif
! get next concentration field
write(adate,'(I6.6)') projjjjmm
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
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
......
......@@ -132,10 +132,10 @@ module mod_flexpart
releases(1:numpoint,2) = juldate(ibdate,ibtime)+real(ireleaseend(:,1)/real(3600*24),kind=8)
releases(numpoint+1:maxpoint,:) = 0
! round dates
releases(1:numpoint,1) = dnint(releases(1:numpoint,1)*1.e5)
releases(1:numpoint,1) = releases(1:numpoint,1)/1.e5
releases(1:numpoint,2) = dnint(releases(1:numpoint,2)*1.e5)
releases(1:numpoint,2) = releases(1:numpoint,2)/1.e5
releases(1:numpoint,1) = dnint(releases(1:numpoint,1)*1.e6)
releases(1:numpoint,1) = releases(1:numpoint,1)/1.e6
releases(1:numpoint,2) = dnint(releases(1:numpoint,2)*1.e6)
releases(1:numpoint,2) = releases(1:numpoint,2)/1.e6
write(logid,*) 'numpoint:', numpoint
write(logid,*) 'RELEASES start, first 26 values:', releases(1:26,1)
write(logid,*) 'RELEASES end, first 26 values:', releases(1:26,2)
......
......@@ -222,7 +222,7 @@ module mod_obs
character(len=max_path_len) :: header
integer :: yyyymmdd, hhmiss
character(len=3) :: recs
character(len=4) :: recs
real(kind=8) :: jdate, avetime
real :: conc, err
integer :: ierr
......@@ -235,7 +235,7 @@ module mod_obs
open(200,file=trim(files%path_output)//'obsread.txt',action='read',status='old',iostat=ierr)
read(200,*) header
do i = 1, nobs
read(200,fmt='(A3,1X,I8,1X,I6,1X,D14.7,1X,F10.6,1X,F10.4,1X,F10.4)') &
read(200,fmt='(A4,1X,I8,1X,I6,1X,D14.7,1X,F10.6,1X,F10.4,1X,F10.4)') &
recs, yyyymmdd, hhmiss, jdate, avetime, conc, err
obs%recs(i) = recs
obs%jdate(i) = jdate
......
......@@ -120,8 +120,8 @@ subroutine read_obs(config, files)
! recs, yyyymmdd, hhmiss, jdate, conc, err, num
! endif
! round jdate
jdate = dnint(jdate*1.e5)/1.e5
avetime = dnint(avetime*1.e4)/1.e4
jdate = dnint(jdate*1.e6)/1.e6
avetime = dnint(avetime*1.e6)/1.e6
if ( ierr.ne.0 ) exit read_loop
if ( jdate.ge.(juldatef+1d0) ) exit read_loop
if ( jdate.lt.juldatei ) cycle read_loop
......
......@@ -145,10 +145,6 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(i))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
filenest = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
filedates = trim(path_flexrec)//'dates'
! read correction factor for dry dir (ratio rho_wet/rho_dry)
! if ( month.ne.prevmonth ) then
......@@ -174,62 +170,84 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! for nested fp run: assume dry air correction is negligible
! read flexpart footprints
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filename)
call system_clock(tread1, countrate)
! for backwards compatability add averaged footprints as option
if ( .not.config%average_fp ) then
call system_clock(tread1, countrate)
if ( .not.config%average_fp ) then
! footprint times match observation timestamps
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(i))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
filenest = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
filedates = trim(path_flexrec)//'dates'
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filename)
! read footprints from binary files
call read_grid(filename, filedates, obs%jdate(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
else
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
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'
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filename)
! read averaged footprints from ncdf files
subdir_fp = trim(obs%recs(i))//'/'//trim(adate)//'/'
filename = trim(files%path_output)//trim(subdir_fp)//'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc'
call read_grid_ncdf(filename, nxgrid, nygrid, ngrid, grid, gtime)
else
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
call system_clock(tread2, countrate)
! write(logid,*) 'simulate: time to read grid (s) = ',(tread2-tread1)/countrate
! correction for dry air (only if obs are mixing ratios)
! do n = 1, ngrid
! ind = minloc(ftime,dim=1,mask=abs(ftime-gtime(n)).eq.minval(abs(ftime-gtime(n))))
! grid(:,:,n) = grid(:,:,n)*factor(:,:,ind)
! end do
! convert from ppt to ppmv or ppbv
grid = grid*config%coeff*mmair/config%molarmass
! convert s.m3/kg to s.m2/kg
grid = grid/outheight(1)
! apply numerical scaling
grid = grid/numscale
else
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
call system_clock(tread2, countrate)
! write(logid,*) 'simulate: time to read grid (s) = ',(tread2-tread1)/countrate
! correction for dry air (only if obs are mixing ratios)
! do n = 1, ngrid
! ind = minloc(ftime,dim=1,mask=abs(ftime-gtime(n)).eq.minval(abs(ftime-gtime(n))))
! grid(:,:,n) = grid(:,:,n)*factor(:,:,ind)
! end do
! convert from ppt to ppmv or ppbv
grid = grid*config%coeff*mmair/config%molarmass
! convert s.m3/kg to s.m2/kg
grid = grid/outheight(1)
! apply numerical scaling
grid = grid/numscale
! read nested flexpart footprints
if ( config%nested ) then
! read nested footprints
inquire ( file=trim(filenest),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filenest)
! for backwards compatability add option for averaged footprints
if ( .not.config%average_fp ) then
if ( .not.config%average_fp ) then
! footprint times match observation timestamps
filenest = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
inquire ( file=trim(filenest),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filenest)
! read footprints from binary files
call read_grid(filenest, filedates, obs%jdate(i), nnxgrid, nnygrid, nnxshift, ngrid, gridnest, gtime)
else
write(logid,*) 'WARNING: cannot find '//trim(filenest)
go to 10
endif
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'
inquire ( file=trim(filenest),exist=lexist )
if ( lexist ) then
! read averaged footprints from ncdf files
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'
call read_grid_ncdf(filenest, nnxgrid, nnygrid, ngrid, gridnest, gtime)
else
write(logid,*) 'WARNING: cannot find '//trim(filenest)
go to 10
endif
! convert from ppt to ppmv or ppbv
gridnest = gridnest*config%coeff*mmair/config%molarmass
! convert s.m3/kg to s.m2/kg
gridnest = gridnest/outheight(1)
! apply numerical scaling
gridnest = gridnest/numscale
else
write(logid,*) 'WARNING: cannot find '//trim(filenest)
go to 10
endif
! convert from ppt to ppmv or ppbv
gridnest = gridnest*config%coeff*mmair/config%molarmass
! convert s.m3/kg to s.m2/kg
gridnest = gridnest/outheight(1)
! apply numerical scaling
gridnest = gridnest/numscale
endif
! calculate indices to state vector for footprint times
......
Supports Markdown
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