Commit 766b9644 authored by ronesy's avatar ronesy
Browse files

Removed optimize offsets option and changes to make more efficient

parent 7a1e2bbe
!---------------------------------------------------------------------------------------
! FLEXINVERT: calc_conc
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> calc_conc
!!
!! Purpose: Calculates fixed mixing ratios.
!!
!! Details: Calculates the fixed components of the mixing ratios
!!
!! Interface:
!!
!! Inputs
!! config - configuration settings data structure
!! fluxes - flux data structure
!! obs - observation data structure
!! ngrid - time dimension of hnest
!! hnest - transport operator nested domain
!! ind - index to low temporal resolution fluxes
!! iobs - index to obs
!!
!---------------------------------------------------------------------------------------
subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1, ix2, jy1, jy2)
use mod_settings
use mod_var
use mod_dates
use mod_fluxes
use mod_obs
implicit none
type (config_t), intent (in) :: config
type (fluxes_t), intent (in) :: fluxes
type (obs_t), intent (in out) :: obs
real(kind=8), dimension(maxngrid), intent (in) :: gtime
real, dimension(nxregrid,nyregrid,ngrid), intent (in) :: hnest
real, dimension(nxgrid,nygrid), intent (in) :: hbkg
integer, intent (in) :: ngrid, iobs, ix1, ix2, jy1, jy2
integer :: n, ns, ni, flag, ind, ilo, ihi
real, dimension(nxgrid,nygrid) :: flxbg
real :: bkgerr, ffferr
! background
! ----------
! index to global fluxes for start of footprint time
flag = 0
ni = 0
do while ( (flag.eq.0).and.(ni.lt.ntstate) )
ni = ni + 1
if ( gtime(1).lt.(fluxes%time(ni)+real(statres,kind=8)) ) flag = 1
enddo
ind = min(ni, ntstate)
! background fluxes
if ( trim(config%spec).eq.'co2' ) then
! CO2 species
flxbg = fluxes%flxnee(:,:,ind) + &
fluxes%flxff(:,:,ind) + &
fluxes%flxocn(:,:,ind)
else if ( trim(config%spec).eq.'ghg' ) then
! GHG species
flxbg = fluxes%flx(:,:,ind)
endif
! background mixing ratio
obs%bkg(iobs) = sum( hbkg(:,:) * flxbg(:,:) )
if (iobs.lt.10) print*, 'calc_conc: obs%bkg(iobs) = ',obs%bkg(iobs)
! background mixing ratio uncertainty
bkgerr = ( sum( hbkg(:,:) * abs(flxbg(:,:)) ) * config%flxerr )**2
! fixed contributions
! -------------------
if ( trim(config%spec).eq.'co2' ) then
! CO2 species
obs%nee(iobs) = 0.
obs%fff(iobs) = 0.
obs%ocn(iobs) = 0.
ffferr = 0.
else
! GHG species
obs%ghg(iobs) = 0.
endif
do n = 1, ngrid
! calculate indice to fluxes
flag = 0
ind = 0
do while ( (flag.eq.0).and.(ind.lt.ntstate) )
ind = ind + 1
if ( gtime(n).lt.(fluxes%time(ind)+real(statres,kind=8)) ) flag = 1
enddo
! index to lo-res fluxes
ilo = min(ind, ntstate)
if ( trim(config%spec).eq.'co2' ) then
! CO2 species
! index to hi-res fluxes (NEE and FFF only)
ihi = minloc(fluxes%timefp,dim=1,mask=abs(fluxes%timefp-gtime(n)).eq.minval(abs(fluxes%timefp-gtime(n))))
! change in mixing ratio
obs%nee(iobs) = obs%nee(iobs) + sum(hnest(:,:,n) * fluxes%flxnee_nest(:,:,ihi))
obs%fff(iobs) = obs%fff(iobs) + sum(hnest(:,:,n) * fluxes%flxff_nest(:,:,ihi))
ffferr = ffferr + (sum(hnest(:,:,n) * fluxes%flxff_nest(:,:,ihi) * config%ffferr))**2
if ( config%nested ) then
obs%ocn(iobs) = obs%ocn(iobs) + sum(hnest(:,:,n) * fluxes%flxocn_nest(:,:,ilo))
else
obs%ocn(iobs) = obs%ocn(iobs) + sum(hnest(:,:,n) * fluxes%flxocn(ix1:ix2,jy1:jy2,ilo))
endif
else
! GHG species
if ( config%nested ) then
! nested domain
obs%ghg(iobs) = obs%ghg(iobs) + sum(hnest(:,:,n) * fluxes%flx_nest(:,:,ilo))
else
! no nested domain
obs%ghg(iobs) = obs%ghg(iobs) + sum(hnest(:,:,n) * fluxes%flx(ix1:ix2,jy1:jy2,ilo))
endif
endif
end do ! ngrid
! observation space error
! -----------------------
if ( trim(config%spec).eq.'co2' ) then
obs%err(iobs) = sqrt(obs%measerr(iobs)**2 + bkgerr + ffferr)
else
obs%err(iobs) = sqrt(obs%measerr(iobs)**2 + bkgerr)
endif
end subroutine calc_conc
......@@ -89,7 +89,7 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
end do
end do
end do
covsum = covsum/real(ndt)
covsum = covsum/real(ndt)**2
! scale to total error of nested domain (unit Tg/y)
toterr = sqrt(covsum)*3600.*24.*365./1.e9/numscale
......
!---------------------------------------------------------------------------------------
! FLEXINVERT: get_initcams
!---------------------------------------------------------------------------------------
! 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_initcams
!!
!! 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_initcams(filename, varname, ntime, concini)
use mod_var
use mod_tools, only : sort
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
integer, intent(in) :: ntime
real, dimension(nxgrid,nygrid,nzgrid,ntime), intent(in out) :: concini
character(len=*), parameter :: lonname='longitude'
character(len=*), parameter :: latname='latitude'
character(len=*), parameter :: timename='time'
integer :: ncid, xid, yid, zid, tid, apid, bpid, varid
integer :: n, kz
integer :: nx, ny, nz, nt, nlev, nzmax
integer :: nx_out, ny_out
logical :: midpoints
real, dimension(:), allocatable :: lon, lat, lonin, latin
real(kind=8), dimension(:), allocatable :: lono, lato
real, dimension(:), allocatable :: ap, bp, pres, mpres, alt
integer, dimension(:), allocatable :: indx, indy
real, dimension(nxgrid) :: lon_out
real, dimension(nygrid) :: lat_out
real, dimension(:,:,:), allocatable :: work1, work2, work3, work4
real, dimension(:,:,:,:), allocatable :: conc, concin
! 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)
! 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) )
call check( nf90_inq_dimid(ncid,'hlevel',zid) )
call check( nf90_inq_dimid(ncid,timename,tid) )
call check( nf90_inquire_dimension(ncid,xid,len=nx) )
call check( nf90_inquire_dimension(ncid,yid,len=ny) )
call check( nf90_inquire_dimension(ncid,zid,len=nlev) )
call check( nf90_inquire_dimension(ncid,tid,len=nt) )
nz = nlev - 1
if ( ntime.ne.nt ) write(logid,*) 'ERROR get_initcams: inconsistent time dimension'
allocate( lon(nx) )
allocate( lat(ny) )
allocate( ap(nlev) )
allocate( bp(nlev) )
allocate( pres(nlev) )
allocate( mpres(nz) )
allocate( alt(nz) )
allocate( conc(nx,ny,nz,nt) )
! 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) )
call check( nf90_inq_varid(ncid,'ap',apid) )
call check( nf90_get_var(ncid,apid,ap) )
call check( nf90_inq_varid(ncid,'bp',bpid) )
call check( nf90_get_var(ncid,bpid,bp) )
call check( nf90_inq_varid(ncid,trim(varname),varid) )
call check( nf90_get_var(ncid,varid,conc(:,:,:,:)) )
call check( nf90_close(ncid) )
! adjust lat and lon of grid
allocate( indy(ny) )
allocate( indx(nx) )
allocate( lato(ny) )
allocate( lono(nx) )
lato = dble(lat)
lono = dble(lon)
call sort(ny, lato, indy)
call sort(nx, lono, indx)
conc = conc(indx,indy,:,:)
if ( lono(nx).ge.180. ) then
allocate( concin((nx-1),ny,nz,nt) )
allocate( lonin((nx-1)) )
concin = conc(1:(nx-1),:,:,:)
lonin = sngl(lono(1:(nx-1)))
else
allocate( concin(nx,ny,nz,nt) )
allocate( lonin(nx) )
concin = conc(:,:,:,:)
lonin = sngl(lono(:))
endif
allocate( latin(ny) )
if ( lato(ny).eq.90. ) then
do n = 1, ny
latin(n) = -90. + real(n - 1)*(180./96.)
end do
else
latin = sngl(lato)
endif
! print*, 'get_initcams: latin = ',latin
! print*, 'get_initcams: lonin = ',lonin
! calculate altitude
do n = 1, nlev
pres(n) = ap(n) + bp(n)*psurf
end do
! avoid pressure of zero
pres(nlev) = 2.
do n = 1, nz
mpres(n) = 0.5*(pres(n) + pres(n+1))
alt(n)=7500*log(psurf/pres(n+1))
end do
! print*, 'get_initcams: pres = ',pres
! print*, 'get_initcams: alt = ',alt
! interpolate to flexpart grid
! ----------------------------
if ( lono(nx).ge.180. ) then
allocate( work1((nx-1),ny,1) )
nx = nx - 1
else
allocate( work1(nx,ny,1) )
endif
! determine upper alt level of conc field for vertical interpolation
nzmax = minloc(alt, dim=1, mask=abs(alt-outheight(nzgrid)).eq.minval(abs(alt-outheight(nzgrid))))
nzmax = min(nz, nzmax+1)
! print*, 'get_initcams: nzmax = ',nzmax
allocate( work2(nxgrid,nygrid,1) )
allocate( work3(nxgrid,nygrid,nzmax) )
allocate( work4(nxgrid,nygrid,nzgrid) )
midpoints = .false.
do n = 1, nt
! interpolate horizontally
do kz = 1, nzmax
work1(:,:,1) = concin(:,:,kz,n)
! gbl_lon and gbl_lat are lon and lat of global flexpart domain
call convertgrid(midpoints, nx, ny, lonin, latin, work1, nx_out, ny_out, lon_out, lat_out, work2)
work3(:,:,kz) = work2(:,:,1)
end do
! interpolate vertically
call interpz(nzmax, alt(1:nzmax), work3, work4)
concini(:,:,:,n) = work4
end do
deallocate( lon )
deallocate( lat )
deallocate( lono )
deallocate( lato )
deallocate( lonin )
deallocate( latin )
deallocate( indx )
deallocate( indy )
deallocate( ap )
deallocate( bp )
deallocate( pres )
deallocate( mpres )
deallocate( alt )
deallocate( conc )
deallocate( concin )
deallocate( work1 )
deallocate( work2 )
deallocate( work3 )
deallocate( work4 )
end subroutine get_initcams
......@@ -54,7 +54,7 @@ subroutine get_initfpctm(filename, varname, concini)
character(len=*), parameter :: levname='lev'
character(len=*), parameter :: timename='time'
integer :: ncid, xid, yid, zid, tid, varid
integer :: kz,jy, ix
integer :: kz
integer :: nx, ny, nz, nt
integer :: nx_out, ny_out
logical :: midpoints
......@@ -79,7 +79,7 @@ subroutine get_initfpctm(filename, varname, concini)
! ----------------
write(logid,*) 'Reading file: '//trim(filename)
preslevel=.true. !CGZ: defaults for noaa interpolated sfc fields
preslevel = .true.
midpoints = .false.
! read dimensions
......@@ -87,38 +87,33 @@ subroutine get_initfpctm(filename, varname, concini)
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")
! 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'
! error code -46 is given when the dimension name does not exist, try if level rather than pressure is available.
write(logid,*) 'WARNING get_initfpctm: pressure not available, test if initial concentration on z-level.'
status = nf90_inq_dimid(ncid,levname,zid)
call check( nf90_inq_dimid(ncid,levname,zid) )
preslevel = .false.
midpoints = .true.
else
call check( nf90_inq_dimid(ncid,prsname,zid) )
call check( nf90_inq_dimid(ncid,prsname,zid) )
endif
!Time does not have to be available, check this here (info lost when using "check")
! 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
! 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
! write error if it was not -46 and get size of time dimension
call check( nf90_inq_dimid(ncid,timename,tid) )
call check( nf90_inquire_dimension(ncid,tid,len=nt) )
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
......@@ -145,7 +140,8 @@ subroutine get_initfpctm(filename, varname, concini)
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
! save height as geopotential height
call check( nf90_get_var(ncid,zid,geoh) )
endif
call check( nf90_inq_varid(ncid,trim(varname),varid) )
......@@ -154,45 +150,26 @@ subroutine get_initfpctm(filename, varname, concini)
call check( nf90_close(ncid) )
! calculate geopotential height for each layer
if(preslevel)then
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
! lon_out and lat_out 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 )
......
!---------------------------------------------------------------------------------------
! FLEXINVERT: get_initfpctm_hour
!---------------------------------------------------------------------------------------
! 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_hour
!!
!! 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_hour(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,nztime), intent(in out) :: concini
character(len=*), parameter :: lonname='longitude'
character(len=*), parameter :: latname='latitude'
character(len=*), parameter :: altname='lev'
character(len=*), parameter :: timename='Date'
integer :: ncid, xid, yid, zid, tid, varid
integer :: kz, it
integer :: nx, ny, nz, nt
integer :: nx_out, ny_out
logical :: midpoints
real, dimension(:), allocatable :: lon, lat, alt, time
real, dimension(nxgrid) :: lon_out
real, dimension(nygrid) :: lat_out
real, dimension(:,:,:), allocatable :: work, workout, workout2, concout
real, dimension(:,:,:,:), allocatable :: conc
! initialize
! ----------
! don't pass globals to subroutines
lon_out = gbl_lon
lat_out = gbl_lat
nx_out = nxgrid
ny_out = nygrid
! read netcdf file
! ----------------