From ea35c372b71668357eadf059bc2f2b713bb1406c Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Mon, 3 Mar 2025 15:59:27 +0100
Subject: [PATCH] Added possibility for different temporal resolution of input
ocean fluxes for global and regional domains and changed init_co2 for prior
fluxes with no leap years
---
settings/SETTINGS_co2_config | 1 +
source/init_co2.f90 | 68 +++++++++++++++++++++---------------
source/initialize.f90 | 1 +
source/mod_settings.f90 | 4 +++
source/mod_var.f90 | 1 +
source/read_obs.f90 | 1 -
6 files changed, 47 insertions(+), 29 deletions(-)
diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config
index b440756..69cfb3a 100644
--- a/settings/SETTINGS_co2_config
+++ b/settings/SETTINGS_co2_config
@@ -135,6 +135,7 @@ coeff_ff_reg: 1.
# other prior fluxes:
# nstep_ocn = time step of other fluxes (integer hours, for monthly data use 720)
nstep_ocn: 720
+nstep_ocn_reg: 720
# Measurement error: unit same as obs
# used if error in obs input <= zero
diff --git a/source/init_co2.f90 b/source/init_co2.f90
index 3e25473..1a5c461 100644
--- a/source/init_co2.f90
+++ b/source/init_co2.f90
@@ -87,7 +87,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
integer, dimension(ndt) :: cnt
integer :: nread, num, numocn
integer :: eomday, eomday_save
- integer :: jjjjmmdd, hhmiss, hh
+ integer :: jjjjmmdd, hhmiss, hh, jjjj, mm
real :: area
real(kind=8) :: jd, jdi
integer :: ix, jy, ix0, jy0
@@ -134,6 +134,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! number of fields to read for each flux time step
nread = int(statres)*24/nt_nee
+ if ( nread.eq.0 ) nread = 1
print*, 'init_co2: nee nread = ',nread
allocate( flx(nxgrid,nygrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
@@ -185,12 +186,6 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
stop
endif
write(logid,*) 'Reading fossil fuel :'//trim(filename)
- if ( nt_ff.lt.24 ) then
- ! assume timestamp is hours since 1-Jan each year
- jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
- else
- jdi = jd
- endif
call read_ncdf(filename, &
files%varname_ff, &
files%lonname_ff, &
@@ -198,7 +193,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
files%timename_ff,&
llx, lly, &
nxgrid, nygrid, &
- jdi, nread, num, &
+ jd, nread, num, &
flx)
print*, 'init_co2: ff num = ',num
! convert from input flux units to kg/m2/s
@@ -297,12 +292,20 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! loop over months reading all input fluxes
do while ( jd.le.juldatef )
call caldate(jd, jjjjmmdd, hhmiss)
- write(adate,'(I4.4)') jjjjmmdd/10000
- write(amonth,'(I2.2)') (jjjjmmdd-(jjjjmmdd/10000)*10000)/100
+ jjjj = jjjjmmdd/10000
+ mm = (jjjjmmdd-jjjj*10000)/100
+ write(adate,'(I4.4)') jjjj
+ write(amonth,'(I2.2)') mm
eomday = calceomday(jjjjmmdd/100)
eomday_save = eomday
! no leap years in input data
- if ( eomday.eq.29 ) eomday = 28
+ if ( (mod(jjjj,4).eq.0) ) then
+ if ( mm.eq.2 ) then
+ eomday = 28
+ else if ( mm.gt.2 ) then
+ jdi = jd - 1d0
+ endif
+ endif
print*, 'init_co2: ffreg eomday = ',eomday
! number of fields to read this month
nread = eomday*24/nt_ff_reg
@@ -321,12 +324,6 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
stop
endif
write(logid,*) 'Reading fossil fuel :'//trim(filename)
- if ( nt_ff_reg.lt.24 ) then
- ! assume timestamp is hours since 1-Jan each year
- jdi = (jd - juldate((jjjjmmdd/10000)*10000+101,0))*24d0 + juldate(19000101,0)
- else
- jdi = jd
- endif
call read_ncdf(filename, &
files%varnest_ff, &
files%lonnest_ff, &
@@ -398,7 +395,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
write(logid,*) 'Reading regional ocean fluxes'
! nested
fluxes%flxocn_nest(:,:,:) = 0.
- numocn = max(1,(nday*24)/nt_ocn)
+ numocn = max(1,(nday*24)/nt_ocn_reg)
print*, 'init_co2: ocn reg timesteps input to read = ',numocn
allocate( flxocn(nxregrid,nyregrid,numocn), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
@@ -411,9 +408,15 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
eomday = calceomday(jjjjmmdd/100)
eomday_save = eomday
! no leap years in input data
- if ( eomday.eq.29 ) eomday = 28
+ if ( (mod(jjjj,4).eq.0) ) then
+ if ( mm.eq.2 ) then
+ eomday = 28
+ else if ( mm.gt.2 ) then
+ jdi = jd - 1d0
+ endif
+ endif
! number of fields to read this month
- nread = eomday*24/nt_ocn
+ nread = eomday*24/nt_ocn_reg
if ( nread.eq.0 ) nread = 1
allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
if ( ierr.ne.0 ) stop 'ERROR init_co2: not enough memory'
@@ -435,12 +438,12 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
files%timenest_ocn,&
rllx, rlly, &
nxregrid, nyregrid, &
- jd, nread, num, &
+ jdi, nread, num, &
flx)
! timestamp in julian days
do i = 1, nread
if ( (n+i).gt.numocn ) exit
- timeocn(n+i) = jd + dble((i-1)*nt_ocn)/24d0
+ timeocn(n+i) = jd + dble((i-1)*nt_ocn_reg)/24d0
flxocn(:,:,n+i) = flx(:,:,i)
end do
deallocate( flx )
@@ -459,21 +462,21 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
! average/interpolate to state vector timestamp
allocate( datain(nxregrid,n) )
allocate( dataout(nxregrid,ntstate) )
- if ( statres.gt.real(nt_ocn/24) ) then
+ if ( statres.gt.real(nt_ocn_reg/24) ) then
! average flux
do jy = 1, nyregrid
datain = flxocn(:,jy,:)
call average(nxregrid, n, timeocn, datain, ntstate, fluxes%time, dataout)
fluxes%flxocn_nest(:,jy,:) = dataout(:,:)
end do
- else if ( statres.lt.real(nt_ocn/24) ) then
+ else if ( statres.lt.real(nt_ocn_reg/24) ) then
! interpolate flux
do jy = 1, nyregrid
datain = flxocn(:,jy,:)
call interp(nxregrid, n, timeocn, datain, ntstate, fluxes%time, dataout)
fluxes%flxocn_nest(:,jy,:) = dataout(:,:)
end do
- else if ( statres.eq.real(nt_ocn/24) ) then
+ else if ( statres.eq.real(nt_ocn_reg/24) ) then
fluxes%flxocn_nest(:,:,:) = flxocn(:,:,:)
endif
deallocate( datain )
@@ -498,12 +501,21 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
jd = juldatei
do while ( jd.le.juldatef )
call caldate(jd, jjjjmmdd, hhmiss)
- write(adate,'(I4.4)') jjjjmmdd/10000
+ jjjj = jjjjmmdd/10000
+ mm = (jjjjmmdd-jjjj*10000)/100
+ write(adate,'(I4.4)') jjjj
+ write(amonth,'(I2.2)') mm
eomday = calceomday(jjjjmmdd/100)
eomday_save = eomday
print*, 'init_co2: eomday = ',eomday
! no leap years in input data
- if ( eomday.eq.29 ) eomday = 28
+ if ( (mod(jjjj,4).eq.0) ) then
+ if ( mm.eq.2 ) then
+ eomday = 28
+ else if ( mm.gt.2 ) then
+ jdi = jd - 1d0
+ endif
+ endif
! number of fields to read this month
nread = eomday*nd_nee
if ( nread.eq.0 ) nread = 1
@@ -526,7 +538,7 @@ subroutine init_co2(files, config, fluxes, obs, states, covar)
files%timenest_nee,&
rllx, rlly, &
nxregrid, nyregrid, &
- jd, nread, num, &
+ jdi, nread, num, &
flx)
do i = 1, nread
if ( (n+i).gt.(nday*nd_nee) ) exit
diff --git a/source/initialize.f90 b/source/initialize.f90
index bab9330..624da68 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -260,6 +260,7 @@ subroutine initialize(files, config)
nt_ff = config%nstep_ff
nt_ff_reg = config%nstep_ff_reg
nt_ocn = config%nstep_ocn
+ nt_ocn_reg = config%nstep_ocn_reg
nt_flx = config%nstep_flx
! number of NEE values per day
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index 5dd545e..ca1acac 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -149,6 +149,7 @@ module mod_settings
integer :: nstep_ff ! timestep of global fossil fuel emissions in hours
integer :: nstep_ff_reg ! timestep of regional fossil fuel emissions in hours
integer :: nstep_ocn ! timestep of ocean prior fluxes in hours
+ integer :: nstep_ocn_reg ! timestep of regional ocean prior fluxes in hours
integer :: nstep_flx ! timestep of GHG prior fluxes in hours
real :: coeff_ff ! coefficient to convert from input global flux units to kg/m2/h
real :: coeff_ff_reg ! coefficient to convert from input regional flux units to kg/m2/h
@@ -820,6 +821,9 @@ module mod_settings
identifier = "nstep_ocn:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%nstep_ocn = int(cn)
+ identifier = "nstep_ocn_reg:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%nstep_ocn_reg = int(cn)
identifier = "nstep_flx:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%nstep_flx = int(cn)
diff --git a/source/mod_var.f90 b/source/mod_var.f90
index 77144ea..195199d 100644
--- a/source/mod_var.f90
+++ b/source/mod_var.f90
@@ -99,6 +99,7 @@ module mod_var
integer :: nt_ff ! timestep of global fossil fuel emissions (hours)
integer :: nt_ff_reg ! timestep of regional fossil fuel emissions (hours)
integer :: nt_ocn ! timestep of ocean prior fluxes (hours)
+ integer :: nt_ocn_reg ! timestep of regional ocean prior fluxes (hours)
integer :: nt_flx ! timestep of other prior fluxes (hours)
integer :: ntstate ! number of time steps in state vector
real :: statres ! time resolution of the state vector (for CO2 1 day) (days)
diff --git a/source/read_obs.f90 b/source/read_obs.f90
index 47dc6c7..e792731 100755
--- a/source/read_obs.f90
+++ b/source/read_obs.f90
@@ -117,7 +117,6 @@ subroutine read_obs(config, files)
! check if data is in inversion time interval
call parse_string(before,"_",args,narg)
adate = args(narg)
- print*, 'read_obs: adate = ',adate
read(adate,*) yyyymmdd
jdate = juldate(yyyymmdd, 0)
if ( jdate.lt.juldatei.or.jdate.gt.juldatef ) cycle
--
GitLab