diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config index b44075636d9d7257c57dbd845e8b8a7990ef8d89..69cfb3a0e4f7a0906b55f6a4bf57ff1ce4dc4828 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 3e25473bc1dea09452c20c32fa56ef3917235a4e..1a5c46182150f99f6a2a1701eafdb851058d49d0 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 bab93305a1dc1292b1b91dc014ad9609175a46f0..624da68bbfdbdd6ac77f1b99d2c02694fb8cd415 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 5dd545e9274665105e652b473e3b320a85e72909..ca1acac1861a34cf75160b24f5fada67c06f17ce 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 77144ea3b866cda08d5b996f3b9246487fc7dd59..195199d1557f40a83c36224bb9ef83dcd8d0c859 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 47dc6c7a56cca46ca4945002f916d703f4fb8c78..e79273171330c43f5240d013f933a000b260d320 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