From 71a3874653ee3eba15d6200da990091f9079679f Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Thu, 13 Jul 2023 16:53:36 +0200
Subject: [PATCH] correction for prep_regions for satellite observations
---
prep_regions/get_surfinf.f90 | 121 ++++++++++++++++++++++-------------
1 file changed, 77 insertions(+), 44 deletions(-)
diff --git a/prep_regions/get_surfinf.f90 b/prep_regions/get_surfinf.f90
index 1cd27f7..4b6eeb8 100644
--- a/prep_regions/get_surfinf.f90
+++ b/prep_regions/get_surfinf.f90
@@ -26,6 +26,7 @@
!! Inputs
!! files - files data structure
!! recs - receptor IDs for each observation
+!! cdryair - total column dry air concentration (mol/m2)
!! obstimes - observation time stamps in julian days
!! avetimes - observation averaging time period in julian days
!! surfinf - total surface influence (zero matrix on input)
@@ -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_dates
@@ -50,12 +51,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
type (config_t), intent (in) :: config
real(kind=8), dimension(maxobs), intent (in) :: obstimes
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
character(len=max_path_len) :: path_flexrec, filename, filedates
character(len=max_name_len) :: varname
- character(len=6) :: adate
+ character(len=6) :: adate, anretr
character(len=8) :: areldate, areltime
logical :: lexist
integer :: ierr
@@ -69,9 +71,10 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
real :: bndx, bndy, delx, dely
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
real(kind=8), dimension(maxpoint,2) :: releases
+ logical :: lsatellite
! indices to regional domain
! --------------------------
@@ -88,6 +91,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
prevmonth = 0
grid(:,:,:) = 0.
+ lsatellite = .false.
do i = 1, nobs
@@ -96,6 +100,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
grid_tmp(:,:,:) = 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
call caldate(obstimes(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
@@ -103,9 +114,13 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
write(areltime,'(I6.6)') hhmiss
month = jjjjmmdd/100
- print*, 'Get footprints for observation nr:', i, recs(i), obstimes(i), jjjjmmdd, hhmiss
- path_flexrec = trim(files%path_flexpart)//trim(recs(i))//'/'//trim(adate)//'/'
+ if ( lsatellite ) then
+ 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'
+ print*, 'path_flexrec = ',path_flexrec
! get header and release information related to current observation location
! (if not read and same for previous observation)
@@ -120,7 +135,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
go to 10
endif
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 )
if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
@@ -133,21 +148,30 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
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(obstimes(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
- 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'
+ ! define grid filenames
+ if ( lsatellite ) then
+ ! satellites
+ write(anretr,'(I6.6)') nretr
+ 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)
+ else
+ ! not satellites
+ ! 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(obstimes(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
+ 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'
+ print*, 'filename = ',filename
! read flexpart footprint
inquire ( file=trim(filename),exist=lexist )
@@ -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
grid_tmp = grid ! Save the footprint to temporary grid for summation
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
- do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obstimes(i)+avetimes(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(recs(i))//'/'//trim(adate)//'/'
- 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'
- filedates = trim(path_flexrec)//'dates'
- call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
- nr_footprints=nr_footprints+1
- release_nr = release_nr+1
- ! sum the grid with previous footprints
- ! check for corresponding output date in gtime
- do standard_grid_time = 1, maxngrid
- 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
+ if ( .not.lsatellite ) then
+ ! check if other footprints need to be read for this observation
+ do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obstimes(i)+avetimes(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(recs(i))//'/'//trim(adate)//'/'
+ 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'
+ filedates = trim(path_flexrec)//'dates'
+ call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
+ nr_footprints=nr_footprints+1
+ release_nr = release_nr+1
+ ! sum the grid with previous footprints
+ ! check for corresponding output date in gtime
+ do standard_grid_time = 1, maxngrid
+ 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 ! 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
grid = grid_tmp/real(nr_footprints)
! convert s.m3/kg to s.m2/kg
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
surfinf(:,:) = surfinf(:,:) + sum(grid(ix1:ix2,jy1:jy2,:),dim=3)
prevmonth = month
--
GitLab