Commit ddf4d9f1 authored by ronesy's avatar ronesy
Browse files

Correction for FP resolutions <1 degree for zoom option and saving of surfinf in ncdf file

parent 12b949dc
......@@ -42,6 +42,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
use mod_dates
use mod_settings
use mod_flexpart
use mod_ncdf
implicit none
......@@ -53,6 +54,7 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
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=8) :: areldate, areltime
logical :: lexist
......@@ -196,5 +198,10 @@ subroutine get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
end do
! write surfinf to ncdf file
filename = trim(files%path_output)//'surfinf.nc'
varname = 'surfinf'
call write_ncdf(files, filename, varname, surfinf)
end subroutine get_surfinf
......@@ -71,7 +71,7 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
real, dimension(nxregrid,nyregrid), intent(in) :: surfinf
real, dimension(nxregrid,nyregrid), intent(in out) :: nbox_xy
real, parameter :: sensitivity_min = 1.e15
real, parameter :: sensitivity_min = 1.e15 ! default 1.e15
integer, parameter :: maxbox = 50000
real :: flux_av
integer :: n, i, ix, jy, iix, jjy, i1, j1, ii, jj
......@@ -100,7 +100,12 @@ subroutine map_grid(config, flux, surfinf, nbox_xy)
lat_lump(:) = 0.
zoomflag = .false.
! mean flux
flux_av = sum(flux)/float(nxregrid*nyregrid)
if ( config%spec.eq.'co2' ) then
! use absolute value of NEE
flux_av = sum(abs(flux))/float(nxregrid*nyregrid)
else
flux_av = sum(flux)/float(nxregrid*nyregrid)
endif
! increase resolution stepwise from coarsest to finest grid
! first iteration at resolution equal to product of mapdec
......@@ -119,14 +124,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)
lon_lump(nfx) = rllx + real((nfx-1)*nres)*dx
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)
lat_lump(nfy) = rlly + real((nfy-1)*nres)*dy
sens_lump(nfx,nfy) = 0.
contrib_lump(nfx,nfy) = 0.
water(nfx,nfy) = 0.
......@@ -186,22 +191,18 @@ 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 ( 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
if ( lon_lump(iix).lt.zllx.or.(lon_lump(iix)+real(nres)*dx).gt.zurx.or. &
lat_lump(jjy).lt.zlly.or.(lat_lump(jjy)+real(nres)*dy).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
......
......@@ -23,26 +23,25 @@ module mod_ncdf
! read_ncdf
! --------------------------------------------------
subroutine read_ncdf(files, config, flux)
subroutine read_ncdf(files, config, filename, varname, datain)
implicit none
type(files_t), intent(in) :: files
type(config_t), intent(in) :: config
real, dimension(nxregrid,nyregrid), intent(in out) :: flux
character(len=max_path_len), intent(in) :: filename
character(len=max_name_len), intent(in) :: varname
real, dimension(nxregrid,nyregrid), intent(in out) :: datain
character(len=max_path_len) :: filename
character(len=max_name_len) :: varname, lonname, latname, timename
character(len=max_name_len) :: lonname, latname, timename
logical :: lexist
integer :: ncid, lon_dimid, lat_dimid, time_dimid
integer :: ncid, lon_dimid, lat_dimid, time_dimid, ncerr
integer :: varid, lonid, latid
integer :: nlon, nlat, nt
integer :: ix, jy
real, dimension(:), allocatable :: lon, lat
real, dimension(:,:,:), allocatable :: flx
filename = trim(files%path_prior)//trim(files%filename_flx)
varname = files%varname_flx
lonname = files%lonname_flx
latname = files%latname_flx
timename = files%timename_flx
......@@ -58,11 +57,16 @@ module mod_ncdf
call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
call check( nf90_inq_dimid(ncid,trim(lonname),lon_dimid) )
call check( nf90_inq_dimid(ncid,trim(latname),lat_dimid) )
call check( nf90_inq_dimid(ncid,trim(timename),time_dimid) )
ncerr = nf90_inq_dimid(ncid,trim(timename),time_dimid)
if ( ncerr.lt.0 ) write(*,*) 'WARNING: file has no time dimension'
call check( nf90_inquire_dimension(ncid,lon_dimid,len=nlon) )
call check( nf90_inquire_dimension(ncid,lat_dimid,len=nlat) )
call check( nf90_inquire_dimension(ncid,time_dimid,len=nt) )
if ( ncerr.eq.0 ) then
call check( nf90_inquire_dimension(ncid,time_dimid,len=nt) )
else
nt = 1
endif
print*, 'read_ncdf: nlon, nlat, nt =', nlon, nlat, nt
......@@ -91,13 +95,21 @@ module mod_ncdf
! read flux
call check( nf90_inq_varid(ncid,trim(varname),varid) )
call check( nf90_get_var(ncid,varid,flx,start=(/ix,jy,1/),count=(/nxregrid,nyregrid,nt/)) )
if ( ncerr.eq.0 ) then
call check( nf90_get_var(ncid,varid,flx,start=(/ix,jy,1/),count=(/nxregrid,nyregrid,nt/)) )
else
call check( nf90_get_var(ncid,varid,flx(:,:,1),start=(/ix,jy/),count=(/nxregrid,nyregrid/)) )
endif
call check( nf90_close(ncid) )
deallocate(lon)
deallocate(lat)
flux = sum(flx,dim=3)/real(nt)
if ( ncerr.eq.0 ) then
datain = sum(flx,dim=3)/real(nt)
else
datain = flx(:,:,1)
endif
deallocate(flx)
......@@ -107,14 +119,16 @@ module mod_ncdf
! write_ncdf
! --------------------------------------------------
subroutine write_ncdf(files, nbox_xy)
subroutine write_ncdf(files, filename, varname, datain)
implicit none
type(files_t), intent(in) :: files
real, dimension(nxregrid,nyregrid), intent(in) :: nbox_xy
character(len=max_path_len), intent(in) :: filename
character(len=max_name_len), intent(in) :: varname
real, dimension(nxregrid,nyregrid), intent(in) :: datain
character(len=max_name_len) :: varname, lonname, latname
character(len=max_name_len) :: lonname, latname
character(len=*), parameter :: varunit = 'none'
character(len=*), parameter :: dimunit = 'degrees'
character(len=*), parameter :: units = 'unit'
......@@ -125,9 +139,8 @@ module mod_ncdf
lonname = files%lonname_regs
latname = files%latname_regs
varname = files%varname_regs
call check( nf90_create(trim(files%file_regions), nf90_clobber, ncid) )
call check( nf90_create(trim(filename), nf90_clobber, ncid) )
call check( nf90_def_dim(ncid, lonname, nxregrid, lon_dimid) )
call check( nf90_def_var(ncid, lonname, nf90_real, lon_dimid, lon_varid) )
......@@ -144,7 +157,7 @@ module mod_ncdf
call check( nf90_put_var(ncid, lon_varid, reg_lon) )
call check( nf90_put_var(ncid, lat_varid, reg_lat) )
call check( nf90_put_var(ncid, varid, nbox_xy) )
call check( nf90_put_var(ncid, varid, datain) )
call check( nf90_close(ncid) )
......
......@@ -40,12 +40,15 @@ program prep_regions
character(max_path_len) :: settings_files
character(max_path_len) :: settings_config
character(max_path_len) :: settings_regions
character(max_path_len) :: filename
character(max_name_len) :: varname
real, dimension(:,:), allocatable :: surfinf
real, dimension(:,:), allocatable :: flux, mask
real(kind=8), dimension(:), allocatable :: obstimes, avetimes
character(4), dimension(:), allocatable :: recs
real, dimension(:,:), allocatable :: nbox_xy
integer :: ierr,ix,jy
logical :: lexist
character(len=20) :: rowfmt
......@@ -86,7 +89,9 @@ program prep_regions
! read prior fluxes
allocate( flux(nxregrid,nyregrid), stat=ierr )
if ( ierr.ne.0 ) stop 'ERROR: not enough memory'
call read_ncdf(files, config, flux)
filename = trim(files%path_prior)//trim(files%filename_flx)
varname = files%varname_flx
call read_ncdf(files, config, filename, varname, flux)
! read observations
allocate( obstimes(maxobs), stat=ierr )
......@@ -97,20 +102,28 @@ program prep_regions
avetimes(:) = 0d0
call read_obs(files, recs, obstimes, avetimes)
print*, 'recs: ',recs(1)
! mean surface influence (i.e. flux sensitivity)
allocate( surfinf(nxregrid,nyregrid), stat=ierr )
if ( ierr.ne.0 ) stop 'ERROR: not enough memory'
surfinf(:,:) = 0.
call get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
inquire ( file = trim(files%path_output)//'surfinf.nc', exist=lexist )
if ( lexist ) then
! read surfinf from file
filename = trim(files%path_output)//'surfinf.nc'
varname = 'surfinf'
call read_ncdf(files, config, filename, varname, surfinf)
else
call get_surfinf(files, config, recs, obstimes, avetimes, surfinf)
endif
! calculate regions
allocate( nbox_xy(nxregrid,nyregrid), stat=ierr )
call map_grid(config, flux, surfinf, nbox_xy)
! write regions file
call write_ncdf(files, nbox_xy)
filename = files%file_regions
varname = 'regions'
call write_ncdf(files, filename, varname, nbox_xy)
! write lsm
allocate( mask(nxregrid,nyregrid), stat=ierr)
......
......@@ -52,7 +52,7 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
real, dimension(16) :: temp
integer :: ierr
integer :: cnt
integer :: narg, n
integer :: narg, n, reclen
integer :: nf, nfiles
integer :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss
real(kind=8) :: jdate, avetime
......@@ -89,9 +89,16 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
do nf = 1, nfiles
! check this file belongs to receptor in reclist
recs = filelist(nf)(1:3)
! added option for 4 character length names
if (filelist(nf)(4:4).eq.'_' )then
recs = filelist(nf)(1:3)
reclen=3
else
recs = filelist(nf)(1:4)
reclen=4
endif
print*, 'read_obs: rec = ',recs
if ( .not.( any(recname(:)(1:3).eq.to_upper(recs)) ) ) go to 10
if ( .not.( any(recname(:)(1:4).eq.to_upper(recs)) ) .and. .not.( any(recname(:)(1:3).eq.to_upper(recs)) )) go to 10
! open input file
open(100,file=trim(files%path_obs)//trim(filelist(nf)),action='read',status='old',iostat=ierr)
......@@ -108,12 +115,22 @@ subroutine read_obs(files, reclist, obstimes, avetimes)
! read data
read_loop: do
! if obs file contains avetime
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
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
else
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 obs file doesn't contain avetime (backwards compatability)
! avetime = 0d0
! 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
! 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
! else
! 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
! endif
if ( ierr.gt.0 ) exit read_loop
if ( jdate.ge.(juldatef+1d0) ) exit read_loop
if ( jdate.lt.juldatei ) cycle read_loop
......
......@@ -33,6 +33,7 @@ subroutine read_reclist(filename)
character(len=200) :: line
integer :: ierr
integer :: cnt
real, dimension(:), allocatable :: reclat, reclon, recalt
! count number of receptors
......@@ -57,13 +58,18 @@ subroutine read_reclist(filename)
! read receptors
allocate ( recname(nrec) )
allocate ( reclat(nrec) )
allocate ( reclon(nrec) )
allocate ( recalt(nrec) )
open(100,file=trim(filename),action='read',status='old',iostat=ierr)
write(*,*) 'Receptors: '
do cnt = 1, nrec
read(100,*) recname(cnt)
write(*,*) recname(cnt)
read(100,fmt='(A4,1X,F6.2,1X,F6.2,1X,F7.2)') recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
write(*,*) recname(cnt), reclat(cnt), reclon(cnt), recalt(cnt)
! read(100,*) recname(cnt)
! write(*,*) recname(cnt)
end do
close(100)
......
Supports Markdown
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