Commit 3d3fb4db authored by cgz's avatar cgz
Browse files

Merge branche 'master' of git.nilu.no:flexpart/flexinvertplus

parents a7d020ca 9fe5d46a
# ==================================================
#
# 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