#------------------------------------------------------------------------------
# Algorithm for Optimal Atmospheric Network Design
#------------------------------------------------------------------------------
#  optimize_network 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.
#
#  optimize_network 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.
#------------------------------------------------------------------------------
# Description
#------------------------------------------------------------------------------
#
#   This algorithm selects sites from a list of potential sites that maximise
#   the Information Content. The selected sites constitute an optimal network
#
#------------------------------------------------------------------------------
# Method
#------------------------------------------------------------------------------
#
#   The algorithm maximise the information content (IC) calculated for 
#   potential networks as calculated as:
#
#   IC = 0.5*ln(|B|) - 0.5*ln(|A|)
#
#   where B and A are the prior and posterior error covariance
#   matrices, respectively. Given that:
#
#   A = (H^TR^(-1)H + B^(-1))^(-1)
#
#   and 
#
#   IC = 0.5*ln(|B/A|)
#
#   then
#
#   IC = 0.5*ln(det(B)*det(H^TR^(-1)H + B^(-1)))
#
#   For computational efficiency the fact that 
#
#   ln(det(M)) = 2*tr(ln(L)) 
#
#   is used where M is any symmetric positive definite matrix and 
#   M = LL^T where L is the Cholesky Lower Triangle of M
#
#   Ref: Thompson and Pisso, A Flexible Algorithm for Network Design Based
#   on Information Theory, Atmos. Meas. Tech., 2022.
#
#------------------------------------------------------------------------------
# Usage
#------------------------------------------------------------------------------
# 
# Inputs:
#   - source-receptor-relationships (SRRs, or "footprints") for each potential
#     site (and existing sites if any) as NetCDF files
#   - prior fluxes as NetCDF files
#   - files of land-sea mask or other region mask (if used) as NetCDF files
#   - list of potential sites as ascii file
#   - list of existing sites (if used) as ascii file
#
# To run:
#   - edit network_config.r 
#   - specify the configuration file below
#
#------------------------------------------------------------------------------
# Written by Rona Thompson, Feb-2022
#------------------------------------------------------------------------------

# user configuration file
config_file<-'network_config.r'

#------------------------------------------------------------------------------

library(parallel)
library(MASS)
library(ncdf4)
source('area_grid.r')

#------------------------------------------------------------------------------

#--- initialization

# read configuration settings
source(config_file)

# number of sources
nsrc=length(ssig)
nsrc=max(1,nsrc)
print(paste('no. sources = ',nsrc,sep=''))

# read site list
stnlist=read.csv(filestn,header=TRUE,sep=',')
nsite=nrow(stnlist)
print(paste('no. sites = ',nsite,sep=''))
slat=stnlist[,'lat']
slon=stnlist[,'lon']
mrerr=stnlist[,'mrerr']
if(isomeas>0) isoerr=stnlist[,'isoerr']
stnlist=stnlist[,'site']

# read existing site list
if(existingTF){
  sitexist=read.csv(fileexist,header=FALSE,sep=',')
  sitexist=as.matrix(sitexist)
  nexist=length(sitexist)
  if(nexist>nsite){ 
    print('ERROR: number of existing sites exceeds total number of sites')
    stop()
  }
  selexist=match(sitexist,stnlist)
}else{
  nexist=0
}

# convert ssig 
ssig=ssig/1000
isoerr=isoerr/1000

#--- calculate B

# read prior flux
nc=nc_open(fileprior)
lon=ncvar_get(nc,lonname)
lat=ncvar_get(nc,latname)
time=ncvar_get(nc,timename)
dx=lon[2]-lon[1]
dy=lat[2]-lat[1]
ix1=which(abs(lon-0.5*dx-llx)==min(abs(lon-0.5*dx-llx)))
ix2=which(abs(lon+0.5*dx-urx)==min(abs(lon+0.5*dx-urx)))
jy1=which(abs(lat-0.5*dy-lly)==min(abs(lat-0.5*dy-lly)))
jy2=which(abs(lat+0.5*dy-ury)==min(abs(lat+0.5*dy-ury)))
lon=lon[ix1:ix2]
lat=lat[jy1:jy2]
nlon=length(lon)
nlat=length(lat)
ntime=length(time)

# loop over variables
vars=nc$var
nvar=length(vars)
# if using sources 
if(nsrc>1) flx=array(0,dim=c(nlon,nlat,nsrc,ntime))
if(nsrc==1) flx=array(0,dim=c(nlon,nlat,ntime))
n=0
for(nv in 1:nvar){
  if(vars[[nv]]$name!='source_names'){
    if(is.element(vars[[nv]]$name,var_exclude)) next
    n=n+1
    print(vars[[nv]]$name)
    tmp=ncvar_get(nc,vars[[nv]]$name)
    tmp=tmp[ix1:ix2,jy1:jy2,]
    tmp[is.na(tmp)]=0
    if(nsrc>1) flx[,,n,]=flx[,,n,]+tmp
    if(nsrc==1) flx=flx+tmp
  }
}
nc_close(nc)

# average fluxes to requested temporal resolution
if(ntime<nfield){
  print('ERROR: insufficient temporal resolution input fluxes')
  stop()
}   
navg=floor(ntime/nfield)
if(nsrc>1) flxavg=array(0,dim=c(nlon,nlat,nsrc,nfield))
if(nsrc==1) flxavg=array(0,dim=c(nlon,nlat,nfield))
if(navg>1){
  nn=1
  n=0
  while(nn<ntime){
    n=n+1
    if(nsrc>1) flxavg[,,,n]=apply(flx[,,,nn:(nn+navg-1)],c(1,2,3),mean)
    if(nsrc==1) flxavg[,,n]=apply(flx[,,nn:(nn+navg-1)],c(1,2),mean)
    nn=nn+navg 
  }
  flx=flxavg
}
flx=flx/numscale
flxerr_ll=flxerr_ll/numscale

# if spatial resolution differs from SRRs
approxy<-function(data) approx(data,x=lat,xout=lat_out,method="linear",rule=2,f=0.5)$y
approxx<-function(data) approx(data,x=lon,xout=lon_out,method="linear",rule=2,f=0.5)$y
if(dx!=xres|dy!=yres){
  lon_out=seq(llx+0.5*xres,urx-0.5*xres,xres)
  lat_out=seq(lly+0.5*yres,ury-0.5*yres,yres)
  nxout=length(lon_out)
  nyout=length(lat_out)
  if(nsrc>1) flx_out=array(0,dim=c(nxout,nyout,nsrc,nfield))
  if(nsrc==1) flx_out=array(0,dim=c(nxout,nyout,nfield))
  for(nf in 1:nfield){
    for(ns in 1:nsrc){
      if(nsrc>1) temp=flx[,,ns,nf]
      if(nsrc==1) temp=flx[,,nf]
      tmp=apply(temp,2,approxx)
      temp=apply(tmp,1,approxy)
      if(nsrc>1) flx_out[,,ns,nf]=t(temp)
      if(nsrc==1) flx_out[,,nf]=t(temp)
    }
  }
  area=matrix(0,length(lon),length(lat))
  for(jy in 1:length(lat)){
    area[,jy]=area_grid(latp=lat[jy],dx=dx,dy=dy)
  }
  area_out=matrix(0,length(lon_out),length(lat_out))
  for(jy in 1:length(lat_out)){
    area_out[,jy]=area_grid(latp=lat_out[jy],dx=xres,dy=yres)
  }
  if(nsrc>1) tot_in=sum(flx[,,1,1]*area)/1e9
  if(nsrc==1) tot_in=sum(flx[,,1]*area)/1e9
  if(nsrc>1) tot_out=sum(flx_out[,,1,1]*area_out)/1e9
  if(nsrc==1) tot_out=sum(flx_out[,,1]*area_out)/1e9
  print(paste('total in = ',tot_in,sep=''))
  print(paste('total out = ',tot_out,sep=''))
  lon=lon_out
  lat=lat_out
  flx=flx_out
  nlon=nxout
  nlat=nyout
}

# calculate prior uncertainty
if(nsrc>1) err=apply(flx,c(1,2,3),mean)*flxerr
if(nsrc==1) err=apply(flx,c(1,2),mean)*flxerr
err[err<flxerr_ll]=flxerr_ll

# read and apply region mask
if(length(filemask)>0){
  nc=nc_open(filemask)
  mlon=ncvar_get(nc,mlonname)
  mlat=ncvar_get(nc,mlatname)
  dx=mlon[2]-mlon[1]
  dy=mlat[2]-mlat[1]
  ix1=which(abs(mlon-0.5*dx-llx)==min(abs(mlon-0.5*dx-llx)))
  ix2=which(abs(mlon+0.5*dx-urx)==min(abs(mlon+0.5*dx-urx)))
  jy1=which(abs(mlat-0.5*dy-lly)==min(abs(mlat-0.5*dy-lly)))
  jy2=which(abs(mlat+0.5*dy-ury)==min(abs(mlat+0.5*dy-ury)))
  mask=ncvar_get(nc,maskname)
  mask=mask[ix1:ix2,jy1:jy2]
  mask[is.na(mask)]=0
  nc_close(nc)
  nerr=length(which(mask>0))
  gridoper=matrix(0,nrow=nerr,ncol=nlon*nlat)
  ii=0
  for(jy in 1:nlat){
    for(ix in 1:nlon){
      if(mask[ix,jy]>0){
        ii=ii+1
        gridoper[ii,((jy-1)*nlon+ix)]=1
      }
    }
  }
  if(nsrc==1){ 
    xerr=gridoper%*%as.vector(err)
  }else{
    xerr=vector('numeric',nerr*nsrc)
    for(ns in 1:nsrc){
      xerr[((ns-1)*nerr+1):(ns*nerr)]=gridoper%*%as.vector(err[,,ns])
    }
  }
}else{
  xerr=as.vector(err)
  nerr=nlon*nlat
}
xerr=abs(xerr)

# read lsm
nc=nc_open(filelsm)
mlon=ncvar_get(nc,mlonname)
mlat=ncvar_get(nc,mlatname)
dx=mlon[2]-mlon[1]
dy=mlat[2]-mlat[1]
ix1=which(abs(mlon-0.5*dx-llx)==min(abs(mlon-0.5*dx-llx)))
ix2=which(abs(mlon+0.5*dx-urx)==min(abs(mlon+0.5*dx-urx)))
jy1=which(abs(mlat-0.5*dy-lly)==min(abs(mlat-0.5*dy-lly)))
jy2=which(abs(mlat+0.5*dy-ury)==min(abs(mlat+0.5*dy-ury)))
lsm=ncvar_get(nc,maskname)
lsm[is.na(lsm)]=0
nc_close(nc)
lsm=lsm[ix1:ix2,jy1:jy2]

# calculate spatial correlation matrix
pih=pi/180
rearth=6.371e6
correl=matrix(0,nrow=nerr,ncol=nerr)
for(ind1 in 1:nerr){
  if(nerr!=(nlon*nlat)) ij1=which(gridoper[ind1,]==1)
  if(nerr==(nlon*nlat)) ij1=ind1
  jy=floor((ij1-1)/nlon)+1
  ix=ij1-(jy-1)*nlon
  lon1=lon[ix]
  lat1=lat[jy]
  if(lsm[ix,jy]==1){
    sigma1=sigma_land*1000
  }else{
    sigma1=sigma_ocean*1000
  }
  for(ind2 in ind1:nerr){
    if(nerr!=(nlon*nlat)) ij2=which(gridoper[ind2,]==1)
    if(nerr==(nlon*nlat)) ij2=ind2
    jy=floor((ij2-1)/nlon)+1
    ix=ij2-(jy-1)*nlon
    lon2=lon[ix]
    lat2=lat[jy]
    if(lsm[ix,jy]==1){
      sigma2=sigma_land*1000
    }else{
      sigma2=sigma_ocean*1000
    }
    if(ind1==ind2){ 
      ddx=0
    }else{
      dd=(1 - sin(lat2*pih)*sin(lat1*pih) - 
           cos(lat2*pih)*cos(lat1*pih)*cos((lon2-lon1)*pih))/2.
      dd=max(dd,0)
      ddx=rearth*2*asin(sqrt(dd))
    }
    cor=exp(-1*ddx/sigma1)
    # no correlation between land and ocean
    if(sigma1!=sigma2) cor=0
    if(cor<0.0001) cor=0
    # and to be sure
    if(ind1==ind2 ) cor=1
    correl[ind1,ind2]=cor
    correl[ind2,ind1]=cor
  }
}

# B matrix spatial correlation
B=matrix(0,nrow=nerr*nsrc,ncol=nerr*nsrc)
for(n in 1:nsrc){
  # assume no correlation between sources
  B[((n-1)*nerr+1):(n*nerr),((n-1)*nerr+1):(n*nerr)]=correl*xerr[((n-1)*nerr+1):(n*nerr)]%*%t(xerr[((n-1)*nerr+1):(n*nerr)])
}
sizeb=object.size(B)
print(sizeb,units="Mb")

# inverse of B
chB=chol(B)
invB=chol2inv(chB)
print('range chB')
print(range(chB))

#--- calculate H

# for potential sites
n=0
H=matrix(0,nrow=nsite,ncol=nerr*nfield*nsrc)
hours=seq(1,24,1)
monthdays=c(31,28,31,30,31,30,31,31,30,31,30,31)
for(stn in stnlist){
  n=n+1
  gridav=array(0,dim=c(nlon,nlat,nfield))
  tempfile=namegen
  tempfile=gsub('YYYY',year,tempfile)
  for(mm in 1:12){
    tempfile1=gsub('MM',sprintf('%02i',mm),tempfile)
    for(dd in 1:monthdays[mm]){
      if((length(seldays)>0)&(!is.element(dd,seldays))) next
      tempfile2=gsub('DD',sprintf('%02i',dd),tempfile1)
      for(hh in hours){
        if((length(selhour)>0)&(!is.element(hh,selhour))) next
        thefile=paste(stn,'/',gsub('HH',sprintf('%02i',(hh-1)),tempfile2),sep='')
        print(paste('reading ',thefile,sep=''))
        nc=nc_open(paste(pathfp,thefile,sep=''))
        if(length(seldays)>0) refday=seldays[1] else refday=1
        if(length(selhour)>0) refhour=selhour[1] else refhour=1
        if(mm==1&dd==refday&hh==refhour){
          glat=ncvar_get(nc,flatname)
          glon=ncvar_get(nc,flonname)
          dx=glon[2]-glon[1]
          dy=glat[2]-glat[1]
          glat=glat+0.5*dy
          glon=glon+0.5*dx
          ix1=which(abs(glon-0.5*dx-llx)==min(abs(glon-0.5*dx-llx)))
          ix2=which(abs(glon+0.5*dx-urx)==min(abs(glon+0.5*dx-urx)))
          jy1=which(abs(glat-0.5*dy-lly)==min(abs(glat-0.5*dy-lly)))
          jy2=which(abs(glat+0.5*dy-ury)==min(abs(glat+0.5*dy-ury)))
        }
        time=ncvar_get(nc,ftimename)
        ntime=length(time)
        tmp=ncvar_get(nc,fvarname,start=c(ix1,jy1,1),count=c((ix2-ix1+1),(jy2-jy1+1),ntime))
        nc_close(nc)
        nt=floor((mm-1)/navg)+1
        if(length(seldays)==0){ 
          gridav[,,nt]=gridav[,,nt]+apply(tmp,c(1,2),sum)/monthdays[mm]/navg
        }else{
          gridav[,,nt]=gridav[,,nt]+apply(tmp,c(1,2),sum)/length(seldays)/navg
        }
      } # hours
    } # days
  } # months
  for(nf in 1:nfield){
    for(ns in 1:nsrc){
      if(isomeas>0){
        # multiply by isotopic ratio
        if(nerr!=(nlon*nlat)) H[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=ssig[ns]*gridoper%*%as.vector(gridav[,,nf])
        if(nerr==(nlon*nlat)) H[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=ssig[ns]*as.vector(gridav[,,nf])
      }else{
        if(nerr!=(nlon*nlat)) H[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=gridoper%*%as.vector(gridav[,,nf])
        if(nerr==(nlon*nlat)) H[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=as.vector(gridav[,,nf])
      }
    }
  }
}

mmair=28.97
H=H*coef*conv*numscale*mmair/molmass/outheight
print('range of H')
print(range(H))

# for existing sites if any
if(nexist>0){
  n=0
  Hexist=matrix(0,nrow=nexist,ncol=nerr*nfield*nsrc)
  for(stn in sitexist){
    n=n+1
    gridav=array(0,dim=c(nlon,nlat,nfield))
    tempfile=namegen
    tempfile=gsub('YYYY',year,tempfile)
    for(mm in 1:12){
      nf=floor((mm-1)/navg)+1           
      tempfile1=gsub('MM',sprintf('%02i',mm),tempfile)
      for(dd in 1:monthdays[mm]){
        if((length(seldays_exist)>0)&(!is.element(dd,seldays_exist))) next
        tempfile2=gsub('DD',sprintf('%02i',dd),tempfile1)
        for(hh in hours){
          if((length(selhour_exist)>0)&(!is.element(hh,selhour_exist))) next
          thefile=paste(stn,'/',gsub('HH',sprintf('%02i',(hh-1)),tempfile2),sep='')
          print(paste('reading ',thefile,sep=''))
          nc=nc_open(paste(pathfp,thefile,sep=''))
          if(length(seldays_exist)>0) refday=seldays_exist[1] else refday=1
          if(length(selhour_exist)>0) refhour=selhour_exist[1] else refhour=1
          if(mm==1&dd==refday&hh==refhour){
            glat=ncvar_get(nc,flatname)
            glon=ncvar_get(nc,flonname)
            dx=glon[2]-glon[1]
            dy=glat[2]-glat[1]
            glat=glat+0.5*dy
            glon=glon+0.5*dx
            ix1=which(abs(glon-0.5*dx-llx)==min(abs(glon-0.5*dx-llx)))
            ix2=which(abs(glon+0.5*dx-urx)==min(abs(glon+0.5*dx-urx)))
            jy1=which(abs(glat-0.5*dy-lly)==min(abs(glat-0.5*dy-lly)))
            jy2=which(abs(glat+0.5*dy-ury)==min(abs(glat+0.5*dy-ury)))
          }
          time=ncvar_get(nc,ftimename)
          ntime=length(time)
          tmp=ncvar_get(nc,fvarname,start=c(ix1,jy1,1),count=c((ix2-ix1+1),(jy2-jy1+1),ntime))
          nc_close(nc)
          nt=floor((mm-1)/navg)+1
          if(length(seldays_exist)==0){
            gridav[,,nt]=gridav[,,nt]+apply(tmp,c(1,2),sum)/length(monthdays[mm])/navg
          }else{
            gridav[,,nt]=gridav[,,nt]+apply(tmp,c(1,2),sum)/length(seldays_exist)/navg
          }
        } # hours
      } # days
    } # months
    for(nf in 1:nfield){
      for(ns in 1:nsrc){
        if(isomeas==2){
          # multiply by isotopic ratio
          if(nerr!=(nlon*nlat)) Hexist[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=ssig[ns]*gridoper%*%as.vector(gridav[,,nf])
          if(nerr==(nlon*nlat)) Hexist[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=ssig[ns]*as.vector(gridav[,,nf])
        }else{
          if(nerr!=(nlon*nlat)) Hexist[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=gridoper%*%as.vector(gridav[,,nf])
          if(nerr==(nlon*nlat)) Hexist[n,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+ns*nerr)]=as.vector(gridav[,,nf])
        }
      }
    }
  }
  Hexist=Hexist*coef*conv*numscale*mmair/molmass/outheight
}
gc()

#--- calculate info content

if(calcinfoTF){

  # function calculating info content
  calcinfo<-function(ssite){
    Hsel=rbind(Hexist,H[ssite,])
    # observation uncertainty
    if(isomeas==0){
      Rsel=c(R,mrerr[ssite])
    }else if(isomeas>0){
      Rsel=c(R,isoerr[ssite])
    }
    iRsel=1/Rsel
    print('iRsel')
    print(iRsel)
    if(largeTF){
      # large problem - split by timestep with B defined for one timestep
      info=0
      for(nf in 1:nfield){
        # H^TR^-1H + B^-1
        work=t(matrix(Hsel[,((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr)],nrow(Hsel),nsrc*nerr))%*%diag(iRsel,nrow(Hsel),ncol=nrow(Hsel))%*%
                 matrix(Hsel[,((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr)],nrow(Hsel),nsrc*nerr) + invB
        chwork=chol(work)
        info = info + sum(diag(log(chB))) + sum(diag(log(chwork)))
      }
    }else{
      # manageable problem with B defined for one timestep
      # H^TR^-1H + B^-1
      work=t(Hsel)%*%diag(iRsel,nrow(Hsel),ncol=nrow(Hsel))%*%Hsel
      for(nf in 1:nfield){
        work[((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr),((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr)]=work[((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr),((nf-1)*nsrc*nerr+1):(nf*nsrc*nerr)]+invB
      }
      chwork=chol(work)
      info = nfield*sum(diag(log(chB))) + sum(diag(log(chwork)))
    }
    gc()
    print(paste('site, info = ',ssite,', ',info,sep=''))
    return(info)
  }
  
  # iterate over number sites to select
  sitepool=seq(1,nsite)
  optsites=vector('numeric',nsel)
  optinfo=vector('numeric',nsel)
  if(nexist==0) Hexist=NULL
  if(nexist==0) R=NULL
  if(nexist>0){
    if(isomeas<2) R=mrerr[selexist]
    if(isomeas==2) R=isoerr[selexist]
  }
  for(n in 1:nsel){
    # call calcinfo function on site pool using parallelization
    info<-mclapply(sitepool, calcinfo, mc.cores=ncores, mc.preschedule=TRUE)
    # equivalent but without parallelization
#    info<-lapply(sitepool, calcinfo)
    info=unlist(info)
    sel=which.max(info)
    ssel=sitepool[sel]
    print(paste('selected site: ',ssel,sep=''))
    sitepool=setdiff(sitepool,ssel)
    optsites[n]=ssel
    optinfo[n]=info[sel]
    Hexist=rbind(Hexist,H[ssel,])
    if(isomeas==0) R=c(R,mrerr[ssel])
    if(isomeas>0) R=c(R,isoerr[ssel])
  }

  # write output
  write.table(cbind(optsites,optinfo),file=paste(pathout,'info_opt_sites.txt',sep=''),append=FALSE,row.names=FALSE,col.names=FALSE,quote=FALSE)

} # calcinfoTF

#--- save matrices

x=ncdim_def('longitude','degrees',lon)
y=ncdim_def('latitude','degrees',lat)
t=ncdim_def('time','none',seq(1,nfield))
s=ncdim_def('source','none',seq(1,nsrc))

# mean footprint all sites
Hmean=apply(H,2,mean)
Hmean=Hmean/numscale
Hgrid=array(0,dim=c(nlon,nlat,nsrc,nfield))
for(nf in 1:nfield){
  for(ns in 1:nsrc){
    if(nerr!=(nlon*nlat)) work=Hmean[((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]%*%gridoper
    if(nerr==(nlon*nlat)) work=Hmean[((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]
    for(ny in 1:nlat){
      Hgrid[,ny,ns,nf]=work[((ny-1)*nlon+1):(ny*nlon)]
    }
  }
}
var=ncvar_def('grid',longname='mean SRR all sites',units='hm3/kg',prec='double',dim=list(x,y,s,t))
nc=nc_create(paste(pathout,'grid_mean_all_',year,'.nc',sep=''),var)
ncvar_put(nc,var,Hgrid)
nc_close(nc)

# if calcinfoTF is FALSE read info_opt_sites from previous run
if(!calcinfoTF){
  optsites=read.table(paste(pathout,'info_opt_sites.txt',sep=''),header=FALSE)
  optsites=optsites[,1]
}

# mean footprint optimal sites
Hopt=apply(H[optsites,],2,mean)
Hopt=Hopt/numscale
Hgrid=array(0,dim=c(nlon,nlat,nsrc,nfield))
for(nf in 1:nfield){
  for(ns in 1:nsrc){
    if(nerr!=(nlon*nlat)) work=Hopt[((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]%*%gridoper
    if(nerr==(nlon*nlat)) work=Hopt[((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]
    for(ny in 1:nlat){
      Hgrid[,ny,ns,nf]=work[((ny-1)*nlon+1):(ny*nlon)]
    }
  }
}
var=ncvar_def('grid',longname='mean SRR optimal sites',units='hm3/kg',prec='double',dim=list(x,y,s,t))
nc=nc_create(paste(pathout,'grid_mean_opt_',year,'.nc',sep=''),var)
ncvar_put(nc,var,Hgrid)
nc_close(nc)

# footprints all sites
Hgrid=array(0,dim=c(nsite,nlon,nlat,nsrc,nfield))
Hgrid=Hgrid/numscale
for(nn in 1:nsite){
  for(nf in 1:nfield){
    for(ns in 1:nsrc){
      if(nerr!=(nlon*nlat)) work=H[nn,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]%*%gridoper
      if(nerr==(nlon*nlat)) work=H[nn,((nf-1)*nsrc*nerr+(ns-1)*nerr+1):((nf-1)*nsrc*nerr+(ns*nerr))]
      for(ny in 1:nlat){
        Hgrid[nn,,ny,ns,nf]=work[((ny-1)*nlon+1):(ny*nlon)]
      }
    }
  }
}
n=ncdim_def('nsite','none',seq(1,nsite))
var=ncvar_def('grid',longname='mean SRR optimal sites',units='hm3/kg',prec='double',dim=list(n,x,y,s,t))
nc=nc_create(paste(pathout,'grid_all_',year,'.nc',sep=''),var)
ncvar_put(nc,var,Hgrid)
nc_close(nc)

# mean flux and uncertainty

if(nsrc>1) flx=apply(flx,c(1,2,4),sum)
flx=apply(flx,c(1,2),mean)
flx=flx*numscale
var=ncvar_def('flux',longname='total annual mean flux',units='kg/m2/h',prec='double',dim=list(x,y))
nc=nc_create(paste(pathout,'flux_mean_',year,'.nc',sep=''),var)
ncvar_put(nc,var,flx)
nc_close(nc)

if(nsrc>1) err=apply(err,c(1,2),sum)
err=abs(err)*numscale
var=ncvar_def('error',longname='total annual mea flux error',units='kg/m2/h',prec='double',dim=list(x,y))
nc=nc_create(paste(pathout,'fluxerr_mean_',year,'.nc',sep=''),var)
ncvar_put(nc,var,err)
nc_close(nc)