Commit 9fe5d46a authored by ronesy's avatar ronesy
Browse files

Updates to be able to use footprints at regular release times and changes to...

Updates to be able to use footprints at regular release times and changes to make global grid map with a zoom region
parent 6b9570df
# ==================================================
#
# REGION SETTINGS
#
# comment lines start with '#'
# each parameter line starts with 'parameter name:'
#
# ==================================================
# Define map globally with "zoom" (true) or regionally (false)
usezoom: .true.
# Use prior fluxes
# weight surface influnce by annual fluxes (true/false)
useflux: .true.
# Use land use
# aggregate gridcells with > 99% water (true/false)
uselanduse: .true.
# Cut-off fractile
# upper fractile of the sensitivity used to define threshold
cutoff: 0.2
# Number of resolution steps (integer)
nmap: 4
# Resolution steps (integer vector of length nmap)
# e.g. if flexpart resolution is 0.5 degree then: mapdec = 1,2,2
# will give 3 resolution steps:
# 1) 2 x 2 x 1 x 0.5 = 2 degrees
# 2) 2 x 1 x 0.5 = 1 degree
# 3) 1 x 0.5 = 0.5 degree
mapdec: 1,2,2,2
# If usezoom is true
# zllx = lon of western edge of zoom domain
# zlly = lat of southern edge zoom domain
# zurx = lon of eastern edge of zoom domain
# zury = lat of northern edge of zoom domain
zllx: -12
zlly: 38
zurx: 28
zury: 70
......@@ -27,6 +27,7 @@
!! files - files data structure
!! recs - receptor IDs for each observation
!! obstimes - observation time stamps in julian days
!! avetimes - observation averaging time period in julian days
!! surfinf - total surface influence (zero matrix on input)
!!
!! Externals
......@@ -35,7 +36,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine get_surfinf(files, config, recs, obstimes, surfinf)
subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
use mod_var
use mod_dates
......@@ -47,7 +48,8 @@ subroutine get_surfinf(files, config, recs, obstimes, surfinf)
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
real(kind=8), dimension(maxobs), intent (in) :: obstimes
character(3), dimension(maxobs), intent (in) :: recs
real(kind=8), dimension(maxobs), intent (in) :: avetimes
character(4), dimension(maxobs), intent (in) :: recs
real, dimension(nxregrid,nyregrid), intent (in out) :: surfinf
character(len=max_path_len) :: path_flexrec, filename, filedates
......@@ -55,11 +57,20 @@ subroutine get_surfinf(files, config, recs, obstimes, surfinf)
character(len=8) :: areldate, areltime
logical :: lexist
integer :: ierr
integer :: jjjjmmdd, hhmiss
real, dimension(nxgrid,nygrid,ngrid) :: grid
real(kind=8), dimension(ngrid) :: gtime
integer :: jjjjmmdd, hhmiss, month, prevmonth
integer :: i, ix1, ix2, jy1, jy2
real, dimension(nxgrid,nygrid,maxngrid) :: grid ! flexpart flux sensitivity (or footprints)
real, dimension(nxgrid,nygrid,maxngrid) :: grid_tmp ! flexpart flux sensitivity (or footprints), temporary to later calculate average
real(kind=8), dimension(maxngrid) :: gtime ! flux sensitivity time stamp (julian days)
real(kind=8), dimension(maxngrid) :: gtime_tmp ! flux sensitivity time stamp (julian days) of first footprint per observation
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift, new_grid_time, standard_grid_time
integer :: numpoint, release_nr, nr_footprints, ngrid
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
! indices to regional domain
! --------------------------
......@@ -73,36 +84,113 @@ subroutine get_surfinf(files, config, recs, obstimes, surfinf)
! loop over obs times
! -------------------
prevmonth = 0
grid(:,:,:) = 0.
do i = 1, nobs
! reset values
gtime_tmp(:) = 0d0
grid_tmp(:,:,:) = 0.
nr_footprints = 0
! define file name
call caldate(obstimes(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
print*, 'jjjjmmdd = ',jjjjmmdd
print*, 'hhmiss = ',hhmiss
print*, 'recs(i) = ',trim(recs(i))
print*, 'obstimes(i) = ',obstimes(i)
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 (.not.config%nested) filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
filename = trim(path_flexrec)//'header'
! get header and release information related to current observation location
! (if not read and same for previous observation)
if ( i.eq.1 .or. month.ne.prevmonth ) then
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
else
print*, 'WARNING: cannot find ',filename
go to 10
endif
else
if(trim(recs(i)).ne.trim(recs(i-1)))then
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
else
print*, 'WARNING: cannot find ',filename
go to 10
endif
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'
filedates = trim(path_flexrec)//'dates'
print*, 'get_surfinf: grid filename = ',filename
! read flexpart footprint
inquire ( file=trim(filename),exist=lexist )
if ( lexist ) then
write(*,*) 'Reading flexpart :'//trim(filename)
call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, grid, gtime)
call read_grid(filename, filedates, obstimes(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
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
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
end do
end do
end do ! finished summing footprints to get average
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)
! add to total surface influence
surfinf(:,:) = surfinf(:,:) + sum(grid(ix1:ix2,jy1:jy2,:),dim=3)
prevmonth = month
else
write(*,*) 'WARNING: cannot find '//trim(filename)
! file doesn't exist
print*, 'WARNING: cannot find '//trim(filename)
go to 10
endif
surfinf(:,:) = surfinf(:,:) + sum(grid(ix1:ix2,jy1:jy2,:),dim=3)
endif
10 continue
......
......@@ -53,9 +53,9 @@ subroutine initialize(files, config)
character(4) :: ayear
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real(kind=8), dimension(maxpoint,2) :: releases
real :: bndx, bndy, delx, dely
integer :: numx, numy
integer :: numx, numy, xshift, ngrid
integer :: i, ierr
! read receptor list file
......@@ -92,9 +92,10 @@ subroutine initialize(files, config)
print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
numx, numy, bndx, bndy, delx, dely, xshift)
nxgrid = numx
nygrid = numy
nxshift = xshift
llx = bndx
lly = bndy
dx = delx
......@@ -104,10 +105,10 @@ subroutine initialize(files, config)
if ( .not.allocated(glon) ) allocate( glon(nxgrid) )
if ( .not.allocated(glat) ) allocate( glat(nygrid) )
do i = 1, nxgrid
glon(i) = llx + (i-1)*dx
glon(i) = llx + (i-1)*dx + 0.5*dx
end do
do i = 1, nygrid
glat(i) = lly + (i-1)*dy
glat(i) = lly + (i-1)*dy + 0.5*dy
end do
print*, 'initialize: glon = ',glon
print*, 'initialize: glat = ',glat
......@@ -155,16 +156,32 @@ subroutine initialize(files, config)
! initialize grid variables
! -------------------------
rllx = config%w_edge_lon
rlly = config%s_edge_lat
rurx = config%e_edge_lon
rury = config%n_edge_lat
rdx = config%xres
rdy = config%yres
nxregrid = int((rurx - rllx)/rdx)
nyregrid = int((rury - rlly)/rdy)
if ( usezoom ) then
! map defined globally - use input from header file
rllx = llx
rlly = lly
rurx = llx + nxgrid*dx
rury = lly + nygrid*dy
rdx = dx
rdy = dy
nxregrid = nxgrid
nyregrid = nygrid
else
! aggregation calculated over regional domain only
rllx = config%w_edge_lon
rlly = config%s_edge_lat
rurx = config%e_edge_lon
rury = config%n_edge_lat
rdx = config%xres
rdy = config%yres
nxregrid = int((rurx - rllx)/rdx)
nyregrid = int((rury - rlly)/rdy)
endif
ngrid = trajdays*ndgrid
print*, 'initialize: ngrid = ',ngrid
maxngrid = trajdays*ndgrid+2
print*, 'initialize: maxngrid = ',maxngrid
! check region boundaries
! -----------------------
......@@ -189,6 +206,25 @@ subroutine initialize(files, config)
write(*,*) 'ERROR: flexpart dimension does not match domain'
stop
endif
if ( usezoom ) then
! check zoom domain within global domain
if( zllx.lt.llx ) then
write(*,*) 'ERROR: zoom western boundary outside flexpart domain'
stop
endif
if( zlly.lt.lly ) then
write(*,*) 'ERROR: zoom southern boundary outside flexpart domain'
stop
endif
if( zurx.gt.(llx+nxgrid*dx) ) then
write(*,*) 'ERROR: zoom eastern boundary outside flexpart domain'
stop
endif
if( zury.gt.(lly+nygrid*dy) ) then
write(*,*) 'ERROR: zoom northern boundary outside flexpart domain'
stop
endif
endif
! initialize regional lon and lat dimensions
if ( .not.allocated(reg_lon) ) allocate( reg_lon(nxregrid) )
......@@ -216,6 +252,13 @@ subroutine initialize(files, config)
write(*,*) 'lly: ',lly
write(*,*) 'dx: ',dx
write(*,*) 'dy: ',dy
if ( usezoom ) then
write(*,*) 'Zoom domain:'
write(*,*) 'zllx: ',zllx
write(*,*) 'zlly: ',zlly
write(*,*) 'zurx: ',zurx
write(*,*) 'zury: ',zury
endif
write(*,*) 'nxgrid: ',nxgrid
write(*,*) 'nygrid: ',nygrid
write(*,*) 'nzgrid: ',nzgrid
......
......@@ -15,8 +15,8 @@ SRCS = mod_var.f90 \
mod_settings.f90 \
mod_strings.f90 \
mod_dates.f90 \
mod_flexpart.f90 \
mod_ncdf.f90 \
mod_flexpart.f90 \
sort.f90 \
read_reclist.f90 \
initialize.f90 \
......
......@@ -33,6 +33,8 @@
!! given positive numbers.
!!
!! Options (defined in SETTINGS):
!! usezoom - (true/false) global map with fine resolution only over region
!! of interest (i.e. nested region in SETTINGS_config)
!! useflux - (true/false) weight surface influence by annual fluxes
!! uselanduse - (true/false) avoid gridcells with > 99% water
!! nmap - (integer) number of resolution steps to use
......@@ -80,7 +82,10 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
real, dimension(2*maxbox) :: x_map, y_map, sens_map, water_box
real, dimension(nxregrid,nyregrid) :: sens_lump, contrib_lump, water
real, dimension(nxregrid*nyregrid) :: sort_contrib, sort_contrib_inv
logical, dimension(nxregrid,nyregrid) :: mark
logical, dimension(nxregrid,nyregrid) :: mark
logical :: zoomflag
real, dimension(nxregrid) :: lon_lump
real, dimension(nyregrid) :: lat_lump
real :: cut_contrib
! initialize some variables
......@@ -91,6 +96,9 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
sort_contrib(:) = 0.
water_box(:) = 0.
water(:,:) = 0.
lon_lump(:) = 0.
lat_lump(:) = 0.
zoomflag = .false.
! mean flux
flux_av = sum(flux)/float(nxregrid*nyregrid)
......@@ -111,12 +119,14 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
do ix = 1, numx, nres
nfx = nfx + 1
nfy = 0
lon_lump(nfx) = rllx + real((nfx-1)*nres)
numy = (nyregrid/nres)*nres
if (numy.lt.nyregrid) then
numy = numy + nres
endif
do jy = 1, numy, nres
nfy = nfy + 1
lat_lump(nfy) = rlly + real((nfy-1)*nres)
sens_lump(nfx,nfy) = 0.
contrib_lump(nfx,nfy) = 0.
water(nfx,nfy) = 0.
......@@ -176,12 +186,26 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
sort_contrib = sort_contrib_inv
cut_contrib = sort_contrib(nint(cutoff*float(nsort))+1)
print*, 'lon_lump: ',lon_lump
print*, 'lat_lump: ',lat_lump
! for all lumped boxes below threshold, mark lumped box, which will
! subsequently not be subdivided any further
do iix = 1, nfx
do jjy = 1, nfy
if ((.not.mark(iix,jjy)).and.((contrib_lump(iix,jjy).le.cut_contrib)).or.(water(iix,jjy).gt.0.99)) then
if ((water(iix,jjy).gt.0.99).or.(sens_lump(iix,jjy).lt.sensitivity_min)) then
if ( usezoom ) then
if ( lon_lump(iix).lt.zllx.or.(lon_lump(iix)+real(nres)).gt.zurx.or. &
lat_lump(jjy).lt.zlly.or.(lat_lump(jjy)+real(nres)).gt.zury ) then
zoomflag = .true.
else
zoomflag = .false.
endif
endif
if(nres==8) print*, 'zoomflag = ',zoomflag
if ( (.not.mark(iix,jjy)).and.((contrib_lump(iix,jjy).le.cut_contrib).or.zoomflag.or.(water(iix,jjy).gt.0.99)) ) then
if ( zoomflag.or.(water(iix,jjy).gt.0.99).or.(sens_lump(iix,jjy).lt.sensitivity_min) ) then
! if ((.not.mark(iix,jjy)).and.((contrib_lump(iix,jjy).le.cut_contrib)).or.(water(iix,jjy).gt.0.99)) then
! if ((water(iix,jjy).gt.0.99).or.(sens_lump(iix,jjy).lt.sensitivity_min)) then
nbox_invert = nbox_invert + 1
if (nbox_invert.gt.maxbox) then
write(*,*) 'ERROR: maxbox too small'
......@@ -201,7 +225,7 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
nbox_xy(ii,jj) = nb
end do
end do
endif
endif
endif
end do ! nfy
end do ! nfx
......
!---------------------------------------------------------------------------------------
! PREP_REGIONS: mod_flexpart
! FLEXINVERT: mod_flexpart
!---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
......@@ -19,7 +19,7 @@
!
!> mod_flexpart
!!
!! Purpose: Module for interfacing with flexpart output.
!! Purpose: Module for interfacing with flexpart output
!!
!---------------------------------------------------------------------------------------
......@@ -28,11 +28,13 @@ module mod_flexpart
use mod_var
use mod_dates
use netcdf
use mod_ncdf, only : check
implicit none
private
public :: read_header, read_grid, read_init, get_nread, read_dates
public :: read_header, read_grid, read_grid_ncdf, read_init, read_factor, get_nread, read_dates
contains
......@@ -47,33 +49,33 @@ module mod_flexpart
! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
numx, numy, bndx, bndy, delx, dely, xshift)
use mod_var
implicit none
character(len=max_path_len), intent(in) :: filename
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
integer, dimension(maxpoint), intent(in out) :: releases
integer, intent(in out) :: numx, numy
real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion
logical :: lexist
integer :: nageclass, imethod, nspec
integer :: loutstep, loutaver, loutsample
integer :: jjjjmmdd, ihmmss
real :: xp2, yp2
character(len=10), dimension(maxspec) :: species
integer, dimension(:), allocatable :: ireleasestart, ireleaseend
integer, dimension(:), allocatable :: npart, nkind
character(len=45), dimension(:), allocatable :: compoint
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
integer :: n, i, j, ierr
character(len=max_path_len), intent(in) :: filename
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2), intent(in out) :: releases
integer, intent(in out) :: numx, numy, xshift
real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion
logical :: lexist
integer :: nageclass, imethod, nspec
integer :: loutstep, loutaver, loutsample
integer :: jjjjmmdd, ihmmss
real :: xp2, yp2
character(len=10), dimension(maxspec) :: species
integer, dimension(:,:), allocatable :: ireleasestart, ireleaseend
integer, dimension(:), allocatable :: npart, nkind
character(len=45), dimension(:), allocatable :: compoint
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
integer :: n, i, j, ierr
inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then
......@@ -100,8 +102,8 @@ module mod_flexpart
end do
read(100) numpoint
allocate( ireleasestart(numpoint) )
allocate( ireleaseend(numpoint) )
allocate( ireleasestart(numpoint,1) )
allocate( ireleaseend(numpoint,1) )
allocate( xpoint(numpoint) )
allocate( ypoint(numpoint) )
allocate( zpoint1(numpoint) )
......@@ -112,7 +114,7 @@ module mod_flexpart
allocate( compoint(numpoint) )
do i = 1, numpoint
read(100) ireleasestart(i), ireleaseend(i)
read(100) ireleasestart(i,1), ireleaseend(i,1)
read(100) xpoint(i), ypoint(i), xp2, yp2, zpoint1(i), zpoint2(i)
read(100) npart(i), nkind(i)
read(100) compoint(i)
......@@ -126,18 +128,28 @@ module mod_flexpart
read(100) nageclass,(lage(i), i = 1, nageclass)
close(100)
releases(1:numpoint) = ireleasestart(:)
releases(1:numpoint,1) = juldate(ibdate,ibtime)+real(ireleasestart(:,1)/real(3600*24),kind=8)
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
write(*,*) 'numpoint:', numpoint
write(*,*) 'RELEASES start, first 26 values:', releases(1:26,1)
write(*,*) 'RELEASES end, first 26 values:', releases(1:26,2)
ageclass = lage(1)
trajdays = ageclass/3600/24
ndgrid = abs(24*3600/loutstep)
if( real(numx)*delx.eq.360. ) then
! global domain - check offset relative to 180W
nxshift = int((bndx + 180.)/delx)
xshift = int((bndx + 180.)/delx)
bndx = -180.
else
! not global domain
nxshift = 0
xshift = 0
endif
deallocate(ireleasestart)
......@@ -163,33 +175,35 @@ module mod_flexpart
!!
! --------------------------------------------------
subroutine read_grid(filename, filedates, jd, numx, numy, grid, gtime)
subroutine read_grid(filename, filedates, jd, numx, numy, xshift, ngrid, grid, gtime)
use mod_var
implicit none
character(len=max_path_len), intent(in) :: filename, filedates