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

Changes to average/select grid_initial files if used regular releases

parent cf142160
......@@ -26,6 +26,7 @@
!!
!! Inputs
!! files - file data structure
!! config - configuration settings data structure
!! obs - observation data structure
!!
!! Externals
......@@ -35,7 +36,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine init_cini(files, obs)
subroutine init_cini(files, config, obs)
use mod_flexpart
use mod_settings
......@@ -48,6 +49,7 @@ subroutine init_cini(files, obs)
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (obs_t), intent (in out) :: obs
character(max_path_len) :: path_flexrec, filename
......@@ -57,14 +59,21 @@ subroutine init_cini(files, obs)
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
integer :: ix, jy, kz, i, n
integer :: month, prevmonth
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift
integer :: numpoint, release_nr, nr_init
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
! loop over observations
! ----------------------
......@@ -74,35 +83,108 @@ subroutine init_cini(files, obs)
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:,:) = 0.
prevmonth = 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
do i = 1, nobs
! define file name
! reset values
gridini(:,:,:) = 0.
nr_init = 0
! check month of observation
call caldate(jdates(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
month = jjjjmmdd/100
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
! read flexpart initial concentrations weighting
inquire(file=trim(filename),exist=lexist)
if(lexist) then
call read_init(filename, gridini)
! convert from ppt to fraction of one
gridini = gridini*1.e-12
if ( config%average_fp ) then
! match and, if needed, average grid_initial files to observation timestamp
! get header and release information related to current observation location
! (if not already read and same for previous observation)
filename = trim(path_flexrec)//'header'
if ( i.eq.1 .or. month.ne.prevmonth ) then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
endif
endif
! find the first release number corresponding to the current observation
release_nr=1
do while ( int(releases(release_nr,1)*1.e5,kind=8).lt.int(jdates(i)*1.e5,kind=8) .and. release_nr.lt.numpoint )
release_nr = release_nr+1
end do
! define first filename
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
! round seconds
hhmiss = (hhmiss/100)*100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
! print*, 'init_cini: filename = ',filename
! read first flexpart grid_init files
inquire ( file=trim(filename), exist=lexist )
if ( lexist ) then
call read_init(filename, gridini_tmp)
gridini = gridini + gridini_tmp
nr_init = 1
! current release_nr was already read, start with next
release_nr = release_nr + 1
do while ( int(releases(release_nr,2)*1.e5,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1.e5,kind=8) .and. release_nr.lt.numpoint )
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
hhmiss=int(hhmiss/100)*100
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
! print*, 'init_cini: filename next = ',filename
! sum the grid with previous initial grids
call read_init(filename, gridini_tmp)
gridini = gridini + gridini_tmp
release_nr = release_nr+1
nr_init = nr_init + 1
end do ! finished summing grid_initial
! convert from ppt to fraction of one
gridini = gridini*1.e-12
gridini = gridini/real(nr_init)
else
! filename not found
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
endif
! timestamps of grid_init files match observations
! define file name
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
! read flexpart grid_init file
inquire(file=trim(filename),exist=lexist)
if(lexist) then
call read_init(filename, gridini)
! convert from ppt to fraction of one
gridini = gridini*1.e-12
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
endif
endif
! calculate termination time of trajectory
call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
......@@ -142,7 +224,8 @@ subroutine init_cini(files, obs)
call get_initconc(filename, varname, proconcini)
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
proconcini = concini
endif
else
! get current concentration field
......@@ -179,7 +262,8 @@ subroutine init_cini(files, obs)
call get_initconc(filename, varname, proconcini)
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
proconcini = concini
endif
endif
endif
......@@ -210,16 +294,14 @@ subroutine init_cini(files, obs)
10 continue
lastjjjjmm = jjjjmm
prevmonth = month
end do
! reorder cini for observations
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: cini = ',cini
do i = 1, nobs
obs%cini(ind(i),:) = cini(i,:)
end do
! print*, 'init_cini: obs%cini = ',obs%cini
end subroutine init_cini
......
......@@ -646,7 +646,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
call init_cini(files, obs)
call init_cini(files, config, obs)
endif
! add random perturbation for MC
......
......@@ -26,6 +26,7 @@
!!
!! Inputs
!! files - file data structure
!! config - configuration settings data structure
!! obs - observation data structure
!!
!! Externals
......@@ -35,7 +36,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine init_fpctm(files, obs)
subroutine init_fpctm(files, config, obs)
use mod_flexpart
use mod_settings
......@@ -48,23 +49,31 @@ subroutine init_fpctm(files, obs)
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (obs_t), intent (in out) :: obs
character(max_path_len) :: path_flexrec, filename
character(max_name_len) :: varname
character(len=8) :: adate, areldate, areltime, ayear
character(len=8) :: adate, areldate, areltime
logical :: lexist
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
integer :: ix, jy, kz, i, n
integer :: month, prevmonth
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: gridini, gridini_tmp
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift
integer :: numpoint, release_nr, nr_init
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
! loop over observations
! ----------------------
......@@ -74,49 +83,113 @@ subroutine init_fpctm(files, obs)
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:,:) = 0.
prevmonth = 0
! sort obs chronologically for efficiency
jdates = 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
! define file name
! reset values
gridini(:,:,:) = 0.
nr_init = 0
! check month of observation
call caldate(jdates(i), jjjjmmdd, hhmiss)
! round hhmiss to nearest 10 min
! print*, 'init_fpctm: hhmiss = ',hhmiss
! hhmiss = nint(real(hhmiss)/1000.)*1000
! print*, 'init_fpctm: after rounding hhmiss = ',hhmiss
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
month = jjjjmmdd/100
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
print*, 'init_cini: pathflexrec = ',path_flexrec
if ( config%average_fp ) then
! match and, if needed, average grid_initial files to observation timestamp
! get header and release information related to current observation location
! (if not already read and same for previous observation)
filename = trim(path_flexrec)//'header'
if ( i.eq.1 .or. month.ne.prevmonth ) then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
endif
endif
! find the first release number corresponding to the current observation
release_nr=1
do while ( int(releases(release_nr,1)*1e5,kind=8).lt.int(jdates(i)*1e5,kind=8) .and. release_nr.lt.numpoint )
release_nr = release_nr+1
end do
! define first filename
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
! round seconds
hhmiss = (hhmiss/100)*100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
! read first flexpart grid_init files
inquire ( file=trim(filename), exist=lexist )
if ( lexist ) then
call read_init(filename, gridini_tmp)
gridini = gridini + gridini_tmp
nr_init = 1
! current release_nr was already read, start with next
release_nr = release_nr + 1
do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((jdates(i)+obs%avetime(ind(i)))*1e5,kind=8) .and. release_nr.lt.numpoint )
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
hhmiss=int(hhmiss/100)*100
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
! sum the grid with previous initial grids
call read_init(filename, gridini_tmp)
gridini = gridini + gridini_tmp
release_nr = release_nr+1
nr_init = nr_init + 1
end do ! finished summing grid_initial
! convert from ppt to fraction of one
gridini = gridini*1.e-12
gridini = gridini/real(nr_init)
else
! filename not found
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
! read flexpart initial concentrations weighting
inquire(file=trim(filename),exist=lexist)
if(lexist) then
call read_init(filename, gridini)
! convert from ppt to fraction of one
gridini = gridini*1.e-12
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
endif
! timestamps of grid_init files match observations
! define file name
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
! read flexpart grid_init file
inquire(file=trim(filename),exist=lexist)
if(lexist) then
call read_init(filename, gridini)
! convert from ppt to fraction of one
gridini = gridini*1.e-12
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
endif
endif
! calculate termination time of trajectory
call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
jjjjmm = jjjjmmdd/100
!write(logid, *) 'init_fpctm: current month BG:', jjjjmm, ' for obs:', obs%jdate(i)
call flush()
! preceding and proceding months
mm = jjjjmm - (jjjjmm/100)*100
if( mm.gt.1 ) then
......@@ -129,54 +202,41 @@ subroutine init_fpctm(files, obs)
else
projjjjmm = (jjjjmm/100 + 1)*100 + 1
endif
call flush()
!Original settings
jdmid = juldate(jjjjmm*100+14, 0) ! jdmid = juldate(jjjjmm*100+14, 0)
jdpre = juldate(prejjjjmm*100+14, 0) ! jdpre = juldate(prejjjjmm*100+14, 0)
if (lastjjjjmm.ne.0)then !CGZ, add if to prevent 0-date in first timestep and allow for an error in mod_dates.f90 if a 0-date is calculated
jdlast = juldate(lastjjjjmm*100+14, 0) !jdlast = juldate(lastjjjjmm*100+14, 0)
jdmid = juldate(jjjjmm*100+14, 0)
jdpre = juldate(prejjjjmm*100+14, 0)
! prevent 0 date in first time step (need to check this)
if( lastjjjjmm.ne.0 ) then
jdlast = juldate(lastjjjjmm*100+14, 0)
else
jdlast = 0
endif
eomday = calceomday(prejjjjmm)
! get new concentrations if not same month as previous observation
if( jjjjmm.ne.lastjjjjmm ) then
!write(logid,* ) 'init_fpctm: checking mid month:', jdmid, jdlast, 'switches:', lastjjjjmm.ne.0, jdmid.eq.(jdlast+dble(eomday))
if( (lastjjjjmm.ne.0).and.(jdmid.eq.(jdlast+dble(eomday))) ) then
! this is proceeding month
preconcini = concini
concini = proconcini
! get proceeding concentration field
write(adate,'(I6.6)') projjjjmm
write(ayear,'(I4.4)') projjjjmm/100
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
filename = str_replace(filename, 'YYYY', ayear)
call flush()
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
call get_initfpctm(filename, varname, proconcini)
else
write(logid,*) 'WARNING: cannot find ',filename
write(logid,*) 'WARNING; next background file missing, copy present into next to avoid 0 value'
proconcini=concini
!go to 10 !CGZ; switched this off because it will cause a obs%cini(i) of 0, but otherwise observation is still considered in simulation!
!Should choose to remove this observation, or stop simulation to ask for an appropriate concentration field?
write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
proconcini = concini
endif
else
! get current concentration field
write(adate,'(I6.6)') jjjjmm
write(ayear,'(I4.4)') jjjjmm/100
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
filename = str_replace(filename, 'YYYY', ayear)
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
......@@ -187,12 +247,8 @@ subroutine init_fpctm(files, obs)
endif
! get previous concentration field
write(adate,'(I6.6)') prejjjjmm
write(ayear,'(I4.4)') prejjjjmm/100
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
filename = str_replace(filename, 'YYYY', ayear)
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
......@@ -203,21 +259,16 @@ subroutine init_fpctm(files, obs)
endif
! get next concentration field
write(adate,'(I6.6)') projjjjmm
write(ayear,'(I4.4)') projjjjmm/100
filename = str_replace(files%file_initconc, 'YYYYMM', adate)
filename = trim(files%path_initconc)//trim(filename)
filename = str_replace(filename, 'YYYY', ayear)
inquire(file=trim(filename),exist=lexist)
if(lexist) then
varname = files%varname_init
call get_initfpctm(filename, varname, proconcini)
else
write(logid,*) 'WARNING: cannot find ',filename
write(logid,*) 'CGZ WARNING; next background file missing, copy present into next to avoid 0 value'
write(logid,*) 'WARNING: next background file missing, copy present into next to avoid 0 value'
proconcini=concini
!go to 10 !CGZ; switched this off because it will cause a obs%cini(i) of 0, but otherwise observation is still considered in simulation!
!Should choose to remove this observation, or stop simulation to ask for an appropriate concentration field?
endif
endif
endif
......@@ -230,7 +281,6 @@ subroutine init_fpctm(files, obs)
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
......@@ -243,29 +293,19 @@ subroutine init_fpctm(files, obs)
cini(i,n) = cini(i,n) + gridini(ix,jy,kz)*conc
end do
end do
!cgz: check cini
!if(obs%cini(i).lt.0.1) then
! write(logid,*) 'CGZ: zero concentration for months:', prejjjjmm, jjjjmm, projjjjmm, obs%jdate(i)
! write(logid,*) ( (obs%jdate(i) - trajdays).lt.jdmid )
! write(logid,*) sum(preconcini), sum(concini), sum(proconcini), sum(gridini)
!endif
end do
10 continue
lastjjjjmm = jjjjmm
! write(logid,*) 'init_fpctm: finished obs: nr=', i, ', obs=', obs%conc(ind(i)), ', cini=', sum(cini(i,))
prevmonth = month
end do
! reorder cini for observations
! print*, 'init_cini: ind = ',ind
! print*, 'init_cini: cini = ',cini
do i = 1, nobs
obs%cini(ind(i),:) = cini(i,:)
end do
! print*, 'init_cini: obs%cini = ',obs%cini
end subroutine init_fpctm
......
......@@ -376,8 +376,8 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
if ( config%restart.eq.0 ) then
! call init_cini(files, obs)
call init_fpctm(files, obs)
call init_cini(files, config, obs)
! call init_fpctm(files, config, obs)
endif
! add random perturbation for MC
......
......@@ -102,32 +102,26 @@ subroutine read_obs(config, files)
! read data
read_loop: do
! for backwards compatibility use flag to indicate if column "avetime" was included in obs file
if ( config%average_fp ) then
! if avetime included
if (reclen.eq.4) then
read(100,fmt='(A4,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
elseif (reclen.eq.3) then
read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
endif
else
! if avetime not included
avetime = 0d0
if (reclen.eq.4) then
read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, conc, err, num
elseif (reclen.eq.3) then
read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, conc, err, num
endif
! if avetime included in obs file
if (reclen.eq.4) then
read(100,fmt='(A4,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
elseif (reclen.eq.3) then
read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
endif
! if avetime not included (backwards compatability with observations prepared before update)