Commit cf142160 authored by ronesy's avatar ronesy
Browse files

CGZ's changes to use hourly FP releases

parent 240438bb
!---------------------------------------------------------------------------------------
! FLEXINVERT: average_fp
!---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! FLEXINVERT is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> average_fp
!!
!! Purpose: Averages footprints over multiple releases. For instance;
!! averaging 24 hourly releases to get the footprint that is representative
!! for observations averaged over 1 day
!! Note: this also requires that the file with observations includes an
!! averaging or total time per observation. Added to allow more flexibility
!! using the same footprints for different types of simulations and projects.
!!
!! Interface:
!!
!! Inputs
!! files - file data structure
!! config - configuration settings data structure
!! obs - observation data structure
!!
!! Externals
!! caldate
!! get_nread
!! read_dates
!! read_header
!! read_grid
!! save_average_fp
!!
!---------------------------------------------------------------------------------------
subroutine average_fp(files, config, obs)
use mod_flexpart
use mod_settings
use mod_var
use mod_dates
use mod_obs
use mod_tools
use mod_save, only : save_average_fp
implicit none
type (files_t), intent (in) :: files
type (config_t), intent (in) :: config
type (obs_t), intent (in out) :: obs
character(len=max_path_len) :: path_flexrec, subdir_fp, filename, filenest, filename_out
character(len=max_path_len) :: filedates
character(len=6) :: adate
character(len=8) :: areldate, areltime
logical :: lexist
integer :: i
integer :: jjjjmmdd, hhmiss
integer :: month, prevmonth
integer :: countrate
integer :: tread1, tread2
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, dimension(nnxgrid,nnygrid,maxngrid):: gridnest ! flexpart nest flux sensitivity (or footprints)
real, dimension(nnxgrid,nnygrid,maxngrid):: gridnest_tmp ! flexpart nest 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
integer :: ngrid
real :: bndx, bndy, delx, dely
integer :: numx, numy, xshift, new_grid_time, standard_grid_time
integer :: numpoint, release_nr, nr_footprints
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
! loop over observations
! ----------------------
prevmonth = 0
do i = 1, nobs
! reset values
gtime_tmp(:) = 0d0
grid_tmp(:,:,:) = 0.
gridnest_tmp(:,:,:) = 0.
nr_footprints = 0
! Check month of observation
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
write(logid,*) 'Get footprints for observation nr:', i, obs%jdate(i), jjjjmmdd, hhmiss
month = jjjjmmdd/100
write(adate,'(I6.6)') jjjjmmdd/100
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(i))//'/'//trim(adate)//'/'
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
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(i)).ne.trim(obs%recs(i-1)))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely, xshift)
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(obs%jdate(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
filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
filenest = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
filedates = trim(path_flexrec)//'dates'
! read flexpart footprints
!-------------------------
inquire ( file=trim(filename), exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart :'//trim(filename)
call system_clock(tread1, countrate)
! Read first grid_time file from start of release
call read_grid(filename, filedates, obs%jdate(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
! Same for nested grid if applicable
if ( config%nested ) then
! read nested footprints
inquire ( file=trim(filenest),exist=lexist )
if ( lexist ) then
write(logid,*) 'Reading flexpart nest :'//trim(filenest)
call read_grid(filenest, filedates, obs%jdate(i), nnxgrid, nnygrid, nnxshift, ngrid, gridnest, gtime)
else
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
endif
gridnest_tmp = gridnest ! Save the footprint to temporary grid for summation
! Current release_nr was already read, start with next
release_nr = release_nr+1
write(logid,*) 'average_fp: obs%jdate(i)+obs%avetime(i) = ',obs%jdate(i)+obs%avetime(i)
write(logid,*) 'average_fp: releases(release_nr,2) = ', releases(release_nr,2)
do while ( int(releases(release_nr,2)*1e5,kind=8).le.int((obs%jdate(i)+obs%avetime(i))*1e5,kind=8) .and. release_nr.lt.numpoint )
call caldate(releases(release_nr,1), jjjjmmdd, hhmiss)
month = jjjjmmdd/100
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(obs%recs(i))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_time_'//trim(areldate)//trim(areltime)//'_001'
filenest = trim(path_flexrec)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001'
filedates = trim(path_flexrec)//'dates'
write(logid,*) 'Reading flexpart:'//trim(filename)
call read_grid(filename, filedates, obs%jdate(i), nxgrid, nygrid, nxshift, ngrid, grid, gtime)
nr_footprints=nr_footprints+1
! Same for nested grid if applicable
if ( config%nested ) then
write(logid,*) 'Reading flexpart nest:'//trim(filenest)
call read_grid(filenest, filedates, obs%jdate(i), nnxgrid, nnygrid, nnxshift, ngrid, gridnest, gtime)
endif
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)
if ( config%nested ) then
gridnest_tmp(:,:,standard_grid_time)=gridnest_tmp(:,:,standard_grid_time)+gridnest(:,:,new_grid_time)
endif
endif
end do
end do
end do ! Finished summing footprints to get average
! Replace grid by average grid and save
! -------------------------------------
write(logid,*) 'Total nr footprints for this observation:', nr_footprints
grid = grid_tmp/real(nr_footprints)
gridnest = gridnest_tmp/real(nr_footprints)
gtime = gtime_tmp
! Save as netcdf
! output filename same as start date and time of observation
call caldate(obs%jdate(i), jjjjmmdd, hhmiss)
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
! create subdirectory in path_output for files
subdir_fp = trim(obs%recs(i))//'/'//trim(adate)//'/'
call system('mkdir -p '//trim(files%path_output)//trim(obs%recs(i)))
call system('mkdir -p '//trim(files%path_output)//trim(subdir_fp))
filename_out = trim(subdir_fp)//'grid_time_'//trim(areldate)//trim(areltime)//'_001.nc'
call save_average_fp(files, filename_out, nxgrid, nygrid, maxngrid, ngrid, grid, gbl_lon, gbl_lat, gtime)
! If applicable same for nested grid
if ( config%nested ) then
filename_out = trim(subdir_fp)//'grid_time_nest_'//trim(areldate)//trim(areltime)//'_001.nc'
call save_average_fp(files, filename_out, nnxgrid, nnygrid, maxngrid, ngrid, gridnest, nest_lon, nest_lat, gtime)
endif
call system_clock(tread2, countrate)
write(logid,*) 'average_fp: time to average grid (s) = ',(tread2-tread1)/countrate
prevmonth = month
! If filename not found
else
write(logid,*) 'WARNING: cannot find '//trim(filename)
go to 10
endif
10 continue
end do ! loop over observations
end subroutine average_fp
!---------------------------------------------------------------------------------------
! FLEXINVERT: get_initfpctm
!---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! FLEXINVERT is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> get_initfpctm
!!
!! Purpose: Reads the initial concentration files and interpolates these to
!! the flexpart grid
!!
!! Interface:
!!
!! Inputs
!! filename - name of 4D concentration file
!! varname - variable name
!! concini - concentrations on flexpart grid (zero array on input)
!!
!! Externals
!! convertgrid
!! interpz
!! check
!!
!---------------------------------------------------------------------------------------
subroutine get_initfpctm(filename, varname, concini)
use mod_var
use mod_save, only : check
use netcdf
implicit none
character(len=max_path_len), intent(in) :: filename
character(len=max_name_len), intent(in) :: varname
real, dimension(nxgrid,nygrid,nzgrid), intent(in out) :: concini
character(len=*), parameter :: lonname='longitude'
character(len=*), parameter :: latname='latitude'
character(len=*), parameter :: prsname='pressure'
character(len=*), parameter :: levname='lev'
character(len=*), parameter :: timename='time'
integer :: ncid, xid, yid, zid, tid, varid
integer :: kz,jy, ix
integer :: nx, ny, nz, nt
integer :: nx_out, ny_out
logical :: midpoints
real, dimension(:), allocatable :: lon, lat, prs, geoh
real, dimension(nxgrid) :: lon_out
real, dimension(nygrid) :: lat_out
real, dimension(:,:,:), allocatable :: work, workout, concout
real, dimension(:,:,:,:), allocatable :: conc
integer :: status
logical :: preslevel
! initialize
! ----------
! don't pass globals to subroutines
lon_out = gbl_lon
lat_out = gbl_lat
nx_out = nxgrid
ny_out = nygrid
! read netcdf file
! ----------------
write(logid,*) 'Reading file: '//trim(filename)
preslevel=.true. !CGZ: defaults for noaa interpolated sfc fields
midpoints = .false.
! read dimensions
call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
call check( nf90_inq_dimid(ncid,lonname,xid) )
call check( nf90_inq_dimid(ncid,latname,yid) )
!The data may be on pressure or z-levels, check this (info lost when using "check")
status = nf90_inq_dimid(ncid,prsname,zid)
if (status.eq.-46) then
!Error code -46 is given when the dimension name does not exist, try if level rather than pressure is available.
print*, 'Pressure not available, test if initial concentration on z-level.'
status = nf90_inq_dimid(ncid,levname,zid)
print*, status, trim(filename), levname
call check( nf90_inq_dimid(ncid,levname,zid) )
preslevel=.false.
midpoints=.true. !DEFAULT FOR FLEXPART CTM NUDGING FIELDS
print*, 'Read z-lev dim'
else
call check( nf90_inq_dimid(ncid,prsname,zid) )
endif
!Time does not have to be available, check this here (info lost when using "check")
status=nf90_inq_dimid(ncid,timename,tid)
if (status.eq.-46) then
!No time variable, assume it's a mean value
nt=1
else
call check( nf90_inq_dimid(ncid,timename,tid) ) !Make sure to write error if it was not -46
call check( nf90_inquire_dimension(ncid,tid,len=nt) ) !Get size of time dimension
endif
call check( nf90_inquire_dimension(ncid,xid,len=nx) )
! print*, 'Read xid dim'
call check( nf90_inquire_dimension(ncid,yid,len=ny) )
! print*, 'Read yid dim'
call check( nf90_inquire_dimension(ncid,zid,len=nz) )
! print*, 'Read zid dim'
! mean so time dimension should be 1
if( nt.ne.1 ) then
write(logid,*) 'ERROR get_initfpctm: time dimension should have length of 1 in '//trim(filename)
stop
endif
allocate( lon(nx) )
allocate( lat(ny) )
allocate( prs(nz) )
allocate( conc(nx,ny,nz,nt) )
allocate( work(nx,ny,1) )
allocate( workout(nxgrid,nygrid,1) )
allocate( concout(nxgrid,nygrid,nz) )
allocate( geoh(nz) )
! read variables
call check( nf90_inq_varid(ncid,lonname,xid) )
call check( nf90_get_var(ncid,xid,lon) )
call check( nf90_inq_varid(ncid,latname,yid) )
call check( nf90_get_var(ncid,yid,lat) )
if(preslevel)then
call check( nf90_inq_varid(ncid,prsname,zid) )
call check( nf90_get_var(ncid,zid,prs) )
else
call check( nf90_inq_varid(ncid,levname,zid) )
call check( nf90_get_var(ncid,zid,geoh) ) !Save height as geopotential height
endif
call check( nf90_inq_varid(ncid,trim(varname),varid) )
call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
call check( nf90_close(ncid) )
! calculate geopotential height for each layer
if(preslevel)then
do kz = 1, nz
geoh(kz) = (gasc*ts/gc/mmol)*log(psurf/prs(kz))
enddo
else
!print*, 'CGZ, print debug, height lev:', geoh(:)
endif
! print*, 'get_initfpctm: nx, ny, nz = ',nx,ny,nz
! print*, 'get_initfpctm: lon = ',lon
! print*, 'get_initfpctm: lat = ',lat
! interpolate to flexpart grid
! ----------------------------
! interpolate horizontally
do kz = 1, nz
work(:,:,1) = conc(:,:,kz,1)
! gbl_lon and gbl_lat are lon and lat of global flexpart domain
call convertgrid(midpoints, nx, ny, lon, lat, work, nx_out, ny_out, lon_out, lat_out, workout)
concout(:,:,kz) = workout(:,:,1)
end do
! interpolate vertically
call interpz(nz, geoh, concout, concini)
!write(logid,*)'CGZ Writing grid to file:', trim(filename)//'_to_concini.txt'
!write(logid,*)'CGZ grid dimension:', nxgrid, nygrid, nzgrid
!write(logid,*)'CGZ min value grid:', MINVAL(concout(:,:,1)), MINVAL(concout(:,:,:))
!open(100,file=trim(filename)//'_to_concini.txt',action='write',status='replace')
!do kz = 1, nzgrid
! do jy = 1, nygrid
! do ix = 1, nxgrid
! write(100,*) ix,jy,kz, concini(ix,jy,kz)
! end do
! end do
!end do
!close(100)
deallocate( lon )
deallocate( lat )
deallocate( prs )
deallocate( conc )
deallocate( work )
deallocate( workout )
deallocate( concout )
deallocate( geoh )
end subroutine get_initfpctm
!---------------------------------------------------------------------------------------
! FLEXINVERT: init_fpctm
!---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! FLEXINVERT is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> init_fpctm
!!
!! Purpose: Calculates the initial concentration for each observation using the
!! 4D initial concentration fields and the flexpart grid_initial files.
!!
!! Interface:
!!
!! Inputs
!! files - file data structure
!! obs - observation data structure
!!
!! Externals
!! caldate
!! read_init
!! get_initfpctm
!!
!---------------------------------------------------------------------------------------
subroutine init_fpctm(files, obs)
use mod_flexpart
use mod_settings
use mod_strings
use mod_var
use mod_dates
use mod_obs
use mod_tools
implicit none
type (files_t), intent (in) :: files
type (obs_t), intent (in out) :: obs
character(max_path_len) :: path_flexrec, filename
character(max_name_len) :: varname
character(len=8) :: adate, areldate, areltime, ayear
logical :: lexist
integer :: lastjjjjmm, prejjjjmm, projjjjmm
integer :: jjjjmm, jjjjmmdd, hhmiss, mm, eomday
integer :: ix, jy, kz, i, n
real :: conc
real(kind=8) :: jdmid, jdpre, jdlast
real(kind=8), dimension(nobs) :: jdates
real, dimension(nobs,ncini) :: cini
integer, dimension(nobs) :: ind
real, dimension(nxgrid,nygrid,nzgrid) :: gridini
real, dimension(nxgrid,nygrid,nzgrid) :: concini, preconcini, proconcini
! loop over observations
! ----------------------
lastjjjjmm = 0
concini(:,:,:) = 0.
preconcini(:,:,:) = 0.
proconcini(:,:,:) = 0.
cini(:,:) = 0.
! sort obs chronologically for efficiency
jdates = obs%jdate(:)
!print*, 'before sort, init_cini: ind = ',ind
!print*, 'before sort, init_cini: jdates = ',jdates
!print*, 'before sort, init_cini: obs%jdate = ',obs%jdate
call sort(nobs,jdates,ind)
!print*, 'after sort, init_cini: ind = ',ind
!print*, 'after sort,init_cini: jdates = ',jdates
!print*, 'after sort, init_cini: obs%jdate = ',obs%jdate
call flush()
do i = 1, nobs
! define file name
call caldate(jdates(i), jjjjmmdd, hhmiss)
! round hhmiss to nearest 10 min
! print*, 'init_fpctm: hhmiss = ',hhmiss
! hhmiss = nint(real(hhmiss)/1000.)*1000
! print*, 'init_fpctm: after rounding hhmiss = ',hhmiss
write(adate,'(I6.6)') jjjjmmdd/100
write(areldate,'(I8.8)') jjjjmmdd
write(areltime,'(I6.6)') hhmiss
path_flexrec = trim(files%path_flexpart)//trim(obs%recs(ind(i)))//'/'//trim(adate)//'/'
filename = trim(path_flexrec)//'grid_initial_'//trim(areldate)//trim(areltime)//'_001'
print*, 'init_cini: filename = ',filename
! read flexpart initial concentrations weighting
inquire(file=trim(filename),exist=lexist)
if(lexist) then
call read_init(filename, gridini)
! convert from ppt to fraction of one
gridini = gridini*1.e-12
else
write(logid,*) 'WARNING: cannot find ',filename
go to 10
endif
! calculate termination time of trajectory
call caldate(jdates(i) - trajdays, jjjjmmdd, hhmiss)
jjjjmm = jjjjmmdd/100
!write(logid, *) 'init_fpctm: current month BG:', jjjjmm, ' for obs:', obs%jdate(i)
call flush()
! preceding and proceding months
mm = jjjjmm - (jjjjmm/100)*100
if( mm.gt.1 ) then
prejjjjmm = jjjjmm - 1
else