update to prep_syndata for different number of timesteps in original fluxes and analysis.nc

prep_syndata/get_flux.f90

@@ -46,7 +46,7 @@ subroutine get_flux(config, files, flx, state)
 
   type (config_t),                           intent(in)     :: config
   type (files_t),                            intent(in)     :: files
-  real, dimension(nxgrid,nygrid,nt_flx),     intent(in out) :: flx
+  real, dimension(nxgrid,nygrid,nflx),       intent(in out) :: flx
   real, dimension(nvar),                     intent(in out) :: state
 
   character(len=max_path_len)                               :: filename
@@ -62,7 +62,7 @@ subroutine get_flux(config, files, flx, state)
   real, dimension(:,:,:), allocatable                       :: tmp
   real, parameter                                           :: smallnum = 1.e-15
 
-  allocate( flx_time(nt_flx) )
+  allocate( flx_time(nflx) )
   flx_time(:) = 0.
 
   ! read prior fluxes
@@ -101,7 +101,7 @@ subroutine get_flux(config, files, flx, state)
   if ( config%spec.eq.'ghg' ) then
 
     ! GHG species
-    do while ( (jd.lt.juldatef).and.(n.lt.nt_flx) )
+    do while ( (jd.lt.juldatef).and.(n.lt.nflx) )
       n = n + 1
       call caldate(jd, jjjjmmdd, hhmiss)
       write(adate,'(I4.4)') jjjjmmdd/10000
@@ -143,7 +143,7 @@ subroutine get_flux(config, files, flx, state)
   else if ( config%spec.eq.'co2' ) then
 
     ! CO2 species
-    do while ( (jd.le.juldatef).and.(n.lt.nt_flx) )
+    do while ( (jd.le.juldatef).and.(n.lt.nflx) )
       n = n + 1
       call caldate(jd, jjjjmmdd, hhmiss)
       write(adate,'(I4.4)') jjjjmmdd/10000

prep_syndata/initialize.f90

@@ -132,15 +132,18 @@ subroutine initialize(files, config)
     flxtres = aint(365./real(config%nt_flx))
     ntstate = nint(real(nday)/statres)
     nt_flx = nday*nd_nee
+    nflx = nt_flx
   else
     ! statres is time resolution of state variables
     ! flxtres is temporal resolution of fluxes (same as state variables for non-CO2 tracers)
+    ! nflx is the number of input flux timesteps falling within the inversion interval
     ! ntstate is also number of state time invervals 
     statres = config%statres
-    flxtres = config%statres
+    flxtres = aint(365./real(config%nt_flx))
     ntstate = nint(real(nday)/statres)
     ndt = 1
     nt_flx = config%nt_flx
+    nflx = max(1, nday/int(365./real(config%nt_flx)))
   endif
   if ( ntstate.lt.1 ) then
     write(logid,*) 'ERROR initialize: number flux time steps < 1 -> increase time interval'
@@ -155,6 +158,7 @@ subroutine initialize(files, config)
   print*, 'flxtres: ',flxtres
   print*, 'ntstate: ',ntstate
   print*, 'nt_flx: ',nt_flx
+  print*, 'nflx: ',nflx
   print*, 'nday: ',nday
   print*, 'nyr: ',nyr
 

prep_syndata/mod_save.f90

@@ -60,7 +60,7 @@ module mod_save
 
     type (config_t),                       intent (in)     :: config
     type (files_t),                        intent (in)     :: files
-    real, dimension(nxgrid,nygrid,nt_flx), intent (in out) :: flx
+    real, dimension(nxgrid,nygrid,nflx), intent (in out) :: flx
     real, dimension(nvar),                 intent (in out) :: state
 
     character(len=max_name_len)                            :: filename, lonname, latname, timename, varname
@@ -75,8 +75,8 @@ module mod_save
     integer                                                :: d, n, nd, nn, nt, ii, i1, i2 
     integer                                                :: ix, jy, ix1, ix2, jy1, jy2
     real                                                   :: frac, hh, hutc
-    real, dimension(nxregrid,nyregrid,nt_flx)              :: flx_out
-    real, dimension(nt_flx)                                :: time
+    real, dimension(nxregrid,nyregrid,nflx)                :: flx_out
+    real, dimension(nflx)                                  :: time
     real, dimension(:,:,:), allocatable                    :: tmp
     real(kind=8)                                           :: jd
     integer                                                :: yyyymmdd, hhmiss
@@ -149,9 +149,10 @@ module mod_save
 
     if ( config%spec.eq.'ghg' ) then
       ! GHG species
-      jd = juldatei
-      n = 1
-      do while ( (jd.lt.juldatef).and.(n.lt.nt_flx) )
+!      jd = juldatei
+!      n = 1
+!      do while ( (jd.lt.juldatef).and.(n.lt.nflx) )
+      do n = 1, nflx
         nt = minloc(stat_time,dim=1,mask=abs(stat_time-flx_time(n)).eq.minval(abs(stat_time-flx_time(n))))
         nt = min(nt,ntstate)
         print*, 'stat_time(nt), flx_time(n) = ',stat_time(nt), flx_time(n)
@@ -208,14 +209,14 @@ module mod_save
             endif
           end do
         end do
-        jd = jd + statres
-        n = n + 1
       end do
     endif ! GHG
-    print*, 'mod_save: n = ',n
 
     ! replace fluxes in inversion domain with perturbed fluxes    
-    flx(ix1:ix2,jy1:jy2,1:n) = flx_out(:,:,1:n)
+    flx(ix1:ix2,jy1:jy2,:) = flx_out(:,:,:)
+
+!*** FOR SF6
+    where ( flx(:,:,:).lt.0. ) flx(:,:,:) = 0.
 
     ! undo numerical scaling
     flx = flx/numscale

prep_syndata/mod_var.f90

@@ -31,7 +31,7 @@ module mod_var
   integer, parameter :: max_path_len=200          ! max character length of paths and files
   integer, parameter :: max_name_len=50           ! max character length of variable names
   integer, parameter :: recname_len=3             ! length of receptor names
-  real, parameter    :: numscale=1.e20            ! numeric scaling factor
+  real, parameter    :: numscale=1.e12            ! numeric scaling factor
   integer, parameter :: maxpoint=1000             ! max number of releases per receptor per month
   integer, parameter :: maxspec=10                ! max number of species in a flexpart run
   integer, parameter :: maxlev=50                 ! max number of vertical levels in flexpart output
@@ -67,6 +67,7 @@ module mod_var
   integer            :: nyr                       ! number of years in inversion interval
   integer            :: nd_nee                    ! number of NEE values per day (in prior NEE file)
   integer            :: nt_flx                    ! number of flux values per year in prior flux files
+  integer            :: nflx                      ! number of input flux time steps falling within inversion interval
   integer            :: ntstate                   ! number of time steps in state vector
   real               :: statres                   ! time resolution of the state vector (for CO2 1 day) (days) 
   real               :: statres_hr                ! sub-daily time resolution of the state vector (hours) 

prep_syndata/prep_syndata.f90

@@ -60,8 +60,8 @@ program prep_syndata
     stop 'ERROR: need to specify SETTINGS_config file'
   endif
 
-  write(logid,*) 'Using: '//trim(settings_files)
-  write(logid,*) 'Using: '//trim(settings_config)
+  write(*,*) 'Using: '//trim(settings_files)
+  write(*,*) 'Using: '//trim(settings_config)
 
   ! initialization
   ! --------------
@@ -79,7 +79,7 @@ program prep_syndata
   call get_error_cov(files, config)
   
   ! get prior fluxes
-  allocate( flx(nxgrid,nygrid,nt_flx))
+  allocate( flx(nxgrid,nygrid,nflx))
   allocate( state(nvar) )
   flx(:,:,:) = 0.
   state(:) = 0.