Correction of start interval for long averaging intervals and limit of last...

Correction of start interval for long averaging intervals and limit of last release time to end of month

prep_flexpart/process_obs.f90

@@ -69,7 +69,7 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
   integer                                               :: jjjjmmdd, hhmiss, yyyy, mm, dd, hh, mi, ss
   integer                                               :: eomday, hloc
   real                                                  :: conc, err
-  real(kind=8)                                          :: jdate, jdatestart, jdateend, jdt
+  real(kind=8)                                          :: jdate, jdatestart, jdateend, jdt, jdsom, jdeom
   real(kind=8), dimension(nobs,6)                       :: buffer
   real, dimension(nobs)                                 :: work
 
@@ -79,27 +79,43 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
     n = 1
     i = 0
     jdate = jdobs(n)
-!    do while( jdate.lt.jdobs(nobs) )
+    call caldate(jdate, jjjjmmdd, hhmiss)
+    eomday = calceomday(jjjjmmdd/100)
+    jdsom = juldate((jjjjmmdd/100)*100+1, 0)
+    jdeom = jdsom + real(eomday,kind=8)
+    print*, 'jdsom: ',jdsom
+    print*, 'jdeom: ',jdeom
+    ! loop over all observations in current month
     do while( jdate.le.jdobs(nobs) )
       jdt = dint(jdate*1.e6)/1.e6
       if ( settings%intaverage.ge.6. ) then
-        ! if long averaging interval start at 00:00
-        jdatestart = dint(jdate)
+        if ( n.eq.1 ) then
+          ! if long averaging interval start at 00:00
+          jdatestart = dint(jdate)
+        else
+          ! continue regular intervals (start at end of last interval)
+          jdatestart = jdateend
+        endif
       else
-        jdatestart = jdate
+        ! otherwise start at current timestamp
+        jdatestart = jdt
       endif
       if ( settings%intaverage.gt.0. ) then
+        ! averaging interval extends up to jdateend
         jdateend = dint((jdatestart+real(settings%intaverage/24.,kind=8))*1.e6)/1.e6
       else
         jdateend = dint((jdatestart+real(freq/24.,kind=8))*1.e6)/1.e6
       endif
+      ! if interval crosses month set to min of current month
+      jdateend = min(jdateend,jdeom)
+      ! check that not beyond end of run time
       jdateend = min(jdateend,jdatef)
       print*, 'jdate: ',jdate
       print*, 'jdatestart: ',jdatestart
       print*, 'jdateend: ',jdateend
       cnt = 0
       buffer(:,:) = 0.
-!      do while( (jdt.ge.jdatestart).and.(jdt.lt.min(jdobs(nobs),jdateend)) )
+      ! loop over observations in averaging interval (fill buffer)
       do while( (jdt.ge.jdatestart).and.(jdt.lt.jdateend) )
         call caldate(jdt, jjjjmmdd, hhmiss)
         print*, 'jdt:',jdt
@@ -161,12 +177,10 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
         if( cnt.gt.0 ) then
           i = i + 1
           ! round jdi to lsync frequency to match flexpart
-!          call caldate(buffer(1,1), jjjjmmdd, hhmiss)
           call caldate(jdatestart, jjjjmmdd, hhmiss)
           hhmiss = int(real(hhmiss)/real((lsync/60)*100))*lsync*100/60
           obs%jdi(i) = juldate(jjjjmmdd, hhmiss)
           ! round jdf to lsync frequency to match flexpart
-!          call caldate((buffer(cnt,1)+real(freq/24.,kind=8)), jjjjmmdd, hhmiss)
           call caldate(jdateend, jjjjmmdd, hhmiss)
           hhmiss = int(real(hhmiss)/real((lsync/60)*100))*lsync*100/60
           obs%jdf(i) = juldate(jjjjmmdd, hhmiss)
@@ -208,7 +222,7 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
 20  continue
     nobs = i
     print*, 'Number observations after selection/averaging: ',nobs
-  endif
+  endif ! averaging and/or selection
 
   ! if no selection and no averaging
   if( settings%lselect.eq.0.and.settings%laverage.eq.0) then

prep_flexpart/read_wdcgg.f90

@@ -149,7 +149,7 @@ subroutine read_wdcgg(settings, jd, nr, nobs, obs)
     jdateutc = jdate - dint(real(lon,kind=8)*24d0/360d0)/24d0
     if ( conc.le.-999. ) cycle read_loop
     if ( nflag.ne.0. ) cycle read_loop
-    if ( (jdate.lt.jd).or.(jdate.ge.min(jreldatef,(jd+real(eomday)))) ) cycle read_loop
+    if ( (jdateutc.lt.jd).or.(jdate.ge.min(jreldatef,(jd+real(eomday)))) ) cycle read_loop
     cnt = cnt + 1
     jdobs(cnt) = jdateutc
     latobs(cnt) = lat