Skip to content
Snippets Groups Projects
Commit 71a38746 authored by ronesy's avatar ronesy
Browse files

correction for prep_regions for satellite observations

parent 03aa365c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
prep_regions/get_surfinf.f90
+ 77
− 44
View file @ 71a38746
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
!! Inputs !! Inputs
!! files - files data structure !! files - files data structure
!! recs - receptor IDs for each observation !! recs - receptor IDs for each observation
!! cdryair - total column dry air concentration (mol/m2)
!! obstimes - observation time stamps in julian days !! obstimes - observation time stamps in julian days
!! avetimes - observation averaging time period in julian days !! avetimes - observation averaging time period in julian days
!! surfinf - total surface influence (zero matrix on input) !! surfinf - total surface influence (zero matrix on input)
...@@ -36,7 +37,7 @@ ...@@ -36,7 +37,7 @@
!! !!
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) subroutine get_surfinf(files, config, recs, cdryair, obstimes, avetimes, surfinf)
use mod_var use mod_var
use mod_dates use mod_dates
...@@ -50,12 +51,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -50,12 +51,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
type (config_t), intent (in) :: config type (config_t), intent (in) :: config
real(kind=8), dimension(maxobs), intent (in) :: obstimes real(kind=8), dimension(maxobs), intent (in) :: obstimes
real(kind=8), dimension(maxobs), intent (in) :: avetimes real(kind=8), dimension(maxobs), intent (in) :: avetimes
character(4), dimension(maxobs), intent (in) :: recs character(recname_len), dimension(maxobs), intent (in) :: recs
real, dimension(maxobs), intent (in) :: cdryair
real, dimension(nxregrid,nyregrid), intent (in out) :: surfinf real, dimension(nxregrid,nyregrid), intent (in out) :: surfinf
character(len=max_path_len) :: path_flexrec, filename, filedates character(len=max_path_len) :: path_flexrec, filename, filedates
character(len=max_name_len) :: varname character(len=max_name_len) :: varname
character(len=6) :: adate character(len=6) :: adate, anretr
character(len=8) :: areldate, areltime character(len=8) :: areldate, areltime
logical :: lexist logical :: lexist
integer :: ierr integer :: ierr
...@@ -69,9 +71,10 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -69,9 +71,10 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
real :: bndx, bndy, delx, dely real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift, new_grid_time, standard_grid_time integer :: numx, numy, xshift, new_grid_time, standard_grid_time
integer :: numpoint, release_nr, nr_footprints, ngrid integer :: numpoint, release_nr, nr_footprints, ngrid, nretr
integer :: ibdate, ibtime integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases real(kind=8), dimension(maxpoint,2) :: releases
logical :: lsatellite
! indices to regional domain ! indices to regional domain
! -------------------------- ! --------------------------
...@@ -88,6 +91,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -88,6 +91,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
prevmonth = 0 prevmonth = 0
grid(:,:,:) = 0. grid(:,:,:) = 0.
lsatellite = .false.
do i = 1, nobs do i = 1, nobs
...@@ -96,6 +100,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -96,6 +100,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
grid_tmp(:,:,:) = 0. grid_tmp(:,:,:) = 0.
nr_footprints = 0 nr_footprints = 0
! check if observation is a satellite retrieval
read(recs(i),*,iostat=ierr) nretr
print*, 'recs(i), nretr, ierr = ',recs(i), nretr, ierr
if ( ierr.eq.0 ) lsatellite = .true.
print*, 'lsatellite = ',lsatellite
print*, 'nretr = ',nretr
! define file name ! define file name
call caldate(obstimes(i), jjjjmmdd, hhmiss) call caldate(obstimes(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100 write(adate,'(I6.6)') jjjjmmdd/100
...@@ -103,9 +114,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -103,9 +114,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
write(areltime,'(I6.6)') hhmiss write(areltime,'(I6.6)') hhmiss
month = jjjjmmdd/100 month = jjjjmmdd/100
print*, 'Get footprints for observation nr:', i, recs(i), obstimes(i), jjjjmmdd, hhmiss if ( lsatellite ) then
path_flexrec = trim(files%path_flexpart)//trim(recs(i))//'/'//trim(adate)//'/' path_flexrec = trim(files%path_flexsat)//trim(areldate)//'/'
else
path_flexrec = trim(files%path_flexpart)//trim(recs(i))//'/'//trim(adate)//'/'
endif
filename = trim(path_flexrec)//'header' filename = trim(path_flexrec)//'header'
print*, 'path_flexrec = ',path_flexrec
! get header and release information related to current observation location ! get header and release information related to current observation location
! (if not read and same for previous observation) ! (if not read and same for previous observation)
...@@ -120,7 +135,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -120,7 +135,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
go to 10 go to 10
endif endif
else else
if(trim(recs(i)).ne.trim(recs(i-1)))then if ( trim(recs(i)).ne.trim(recs(i-1)).and.(.not.lsatellite) ) then
inquire ( file=trim(filename),exist=lexist ) inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header' print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
...@@ -133,21 +148,30 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -133,21 +148,30 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
endif endif
endif endif
! find the first release number corresponding to the current observation ! define grid filenames
release_nr=1 if ( lsatellite ) then
do while ( int(releases(release_nr,1)*1e5,kind=8).lt.int(obstimes(i)*1e5,kind=8) .and. release_nr.lt.numpoint ) ! satellites
release_nr = release_nr+1 write(anretr,'(I6.6)') nretr
end do if (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//'_'//trim(anretr)
if (config%nested) filename = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//'_'//trim(anretr)
! define first filename else
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss) ! not satellites
! round seconds ! find the first release number corresponding to the current observation
hhmiss = (hhmiss/100)*100 release_nr=1
write(areldate,'(I8.8)') jjjjmmdd do while ( int(releases(release_nr,1)*1e5,kind=8).lt.int(obstimes(i)*1e5,kind=8) .and. release_nr.lt.numpoint )
write(areltime,'(I6.6)') hhmiss release_nr = release_nr+1
if (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001' end do
if (config%nested) filename = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001' ! 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
if (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
if (config%nested) filename = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
endif
filedates = trim(path_flexrec)//'dates' filedates = trim(path_flexrec)//'dates'
print*, 'filename = ',filename
! read flexpart footprint ! read flexpart footprint
inquire ( file=trim(filename),exist=lexist ) inquire ( file=trim(filename),exist=lexist )
...@@ -156,35 +180,44 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf) ...@@ -156,35 +180,44 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
gtime_tmp = gtime ! Save the output date of each field in footprint grid gtime_tmp = gtime ! Save the output date of each field in footprint grid
grid_tmp = grid ! Save the footprint to temporary grid for summation grid_tmp = grid ! Save the footprint to temporary grid for summation
nr_footprints = 1 ! Start counting number of footprints to be averaged nr_footprints = 1 ! Start counting number of footprints to be averaged
! current release_nr was already read, start with next ! current release_nr was already read, start with next
release_nr = release_nr+1 release_nr = release_nr+1
do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obstimes(i)+avetimes(i))*1e5,kind=8) .and. release_nr.lt.numpoint ) if ( .not.lsatellite ) then
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss) ! check if other footprints need to be read for this observation
write(adate,'(I6.6)') jjjjmmdd/100 do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obstimes(i)+avetimes(i))*1e5,kind=8) &
write(areldate,'(I8.8)') jjjjmmdd .and. release_nr.lt.numpoint )
hhmiss=int(hhmiss/100)*100 call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
write(areltime,'(I6.6)') hhmiss write(adate,'(I6.6)') jjjjmmdd/100
path_flexrec = trim(files%path_flexpart)//trim(recs(i))//'/'//trim(adate)//'/' write(areldate,'(I8.8)') jjjjmmdd
if (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001' hhmiss=int(hhmiss/100)*100
if (config%nested) filename = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001' write(areltime,'(I6.6)') hhmiss
filedates = trim(path_flexrec)//'dates' path_flexrec = trim(files%path_flexpart)//trim(recs(i))//'/'//trim(adate)//'/'
call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime) if (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
nr_footprints=nr_footprints+1 if (config%nested) filename = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
release_nr = release_nr+1 filedates = trim(path_flexrec)//'dates'
! sum the grid with previous footprints call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
! check for corresponding output date in gtime nr_footprints=nr_footprints+1
do standard_grid_time = 1, maxngrid release_nr = release_nr+1
do new_grid_time = 1, maxngrid ! sum the grid with previous footprints
if ( gtime(new_grid_time) .eq. gtime_tmp(standard_grid_time)) then ! check for corresponding output date in gtime
grid_tmp(:,:,standard_grid_time) = grid_tmp(:,:,standard_grid_time)+grid(:,:,new_grid_time) do standard_grid_time = 1, maxngrid
endif do new_grid_time = 1, maxngrid
if ( gtime(new_grid_time) .eq. gtime_tmp(standard_grid_time)) then
grid_tmp(:,:,standard_grid_time) = grid_tmp(:,:,standard_grid_time)+grid(:,:,new_grid_time)
endif
end do
end do end do
end do end do ! finished summing footprints to get average
end do ! finished summing footprints to get average endif ! lsatellite
! (needed if combine satellite and ground-based footprints for grid definition)
print*, 'Total nr footprints for this observation:', nr_footprints print*, 'Total nr footprints for this observation:', nr_footprints
grid = grid_tmp/real(nr_footprints) grid = grid_tmp/real(nr_footprints)
! convert s.m3/kg to s.m2/kg ! convert s.m3/kg to s.m2/kg
grid = grid/outheight(1) grid = grid/outheight(1)
if ( lsatellite ) then
! convert from column SRR units to that consistent with mole fractions
grid = grid/cdryair(i)
endif
! add to total surface influence ! add to total surface influence
surfinf(:,:) = surfinf(:,:) + sum(grid(ix1:ix2,jy1:jy2,:),dim=3) surfinf(:,:) = surfinf(:,:) + sum(grid(ix1:ix2,jy1:jy2,:),dim=3)
prevmonth = month prevmonth = month
......
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