[QUESTION]How does obspack_mod.F90 sample?

[helpyuan]

Ask a question about GEOS-Chem:

Hi all,

I wonder what is the sampling method in obspack_mod.F90 module? I did not find the keywords about interpolation in obspack_mod.F90, and I found that the subroutine "ObsPack_Get_Indices" seemed to select the grid point closest to the longitude, latitude, and height of the site. So is this the way to sample?

In addition, I have a question. In the subroutine "ObsPack_Get_Indices", why do I and J add 1.5d0 ?

I = INT( ( State_Diag%ObsPack_Longitude(iObs) + 180.0_f8                 &
                                                  - ( I0 * State_Grid%DX ) ) &
                                                  / State_Grid%DX + 1.5d0  )

J = INT( ( State_Diag%ObsPack_Latitude(iObs)  +  90.0_f8                 &
                                                  - ( J0 * State_Grid%DY ) ) &
                                                  / State_Grid%DY + 1.5d0  )

I appreciate any help.

[yantosca]

Thanks for writing @helpyuan. In the ObsPack source code the sampling interval is construced symmetrically around the ObsPack observation time:

       !====================================================================
       ! Construct the sampling intervals from the center time
       !====================================================================

       ! Pick the start and end time of the averaging interval
       ! depending on the averaging strategy listed in the file
       SELECT CASE ( State_Diag%ObsPack_Strategy(N) )

          !------------------
          ! DO-NOT-SAMPLE
          ! set start > end
          !------------------
          CASE( 0 )

             IF ( Window_Start_Time ) THEN
                State_Diag%ObsPack_Ival_Center(N) =                          &
                     State_Diag%ObsPack_Ival_Start(N)  - 1.0_f8
             ELSE
                State_Diag%ObsPack_Ival_Start(N) =                           &
                     State_Diag%ObsPack_Ival_Center(N) + 1.0_f8
             ENDIF

             State_Diag%ObsPack_Ival_End(N) =                                &
                  State_Diag%ObsPack_Ival_Center(N)    - 1.0_f8

          !------------------
          ! 4-hour window
          !------------------
          CASE( 1 )

             IF ( Window_Start_Time ) THEN
                State_Diag%ObsPack_Ival_Center(N) =                          &
                     State_Diag%ObsPack_Ival_Start(N)  + 7200.0_f8
             ELSE
                State_Diag%ObsPack_Ival_Start(N) =                           &
                     State_Diag%ObsPack_Ival_Center(N) - 7200.0_f8
             ENDIF

             State_Diag%ObsPack_Ival_End(N) =                                &
                  State_Diag%ObsPack_Ival_Center(N)    + 7200.0_f8

          !------------------
          ! 1-hour window
          !------------------
          CASE( 2 )

             IF ( Window_Start_Time ) THEN
                State_Diag%ObsPack_Ival_Center(N) =                          &
                     State_Diag%ObsPack_Ival_Start(N)  + 1800.0_f8
             ELSE
                State_Diag%ObsPack_Ival_Start(N) =                           &
                     State_Diag%ObsPack_Ival_Center(N) - 1800.0_f8
             ENDIF

             State_Diag%ObsPack_Ival_End(N) =                                &
                  State_Diag%ObsPack_Ival_center(N)    + 1800.0_f8

          !------------------
          ! 90-minute window
          !------------------
          CASE( 3 )

             IF ( Window_Start_Time ) THEN
                State_Diag%ObsPack_Ival_Center(N) =                           &
                     State_Diag%ObsPack_Ival_Start(N)  + 2700.0_f8
             ELSE

                State_Diag%ObsPack_Ival_Start(N) =                           &
                     State_Diag%ObsPack_Ival_Center(N) - 2700.0_f8
             ENDIF

             State_Diag%ObsPack_Ival_End(N) =                                &
                  State_Diag%ObsPack_Ival_Center(N)    + 2700.0_f8

          !---------------------
          ! Instaneous Sampling
          !---------------------
          CASE( 4 )

             IF ( Window_Start_Time ) THEN
                State_Diag%ObsPack_Ival_Center(N) =                          &
                     State_Diag%ObsPack_Ival_Start(N)
             ELSE
                State_Diag%ObsPack_Ival_Start(N) =                           &
                     State_Diag%ObsPack_Ival_Center(N)
             ENDIF

             State_Diag%ObsPack_Ival_End(N) =                                &
                  State_Diag%ObsPack_Ival_Center(N)

          !------------------
          ! Exit w/ error
          !------------------
          CASE DEFAULT
             ErrMsg = 'Observation with ObsPack_id: '                     // &
                      TRIM( ADJUSTL( State_Diag%ObsPack_id(N) ) )         // &
                      ' has an unknown or invalid sampling strategy!'
             CALL GC_Error( ErrMsg, RC, ThisLoc )
             RETURN

       END SELECT

We assume that the time listed is the start of the averaging interval (unless specified otherwise in the ObsPack netCDF file... some older versions of the ObsPack data have the reported time in the file being the center of the averaging interval, not the start. So the GEOS-Chem Obspack diagnostic can handle both instances).

If you look at the 1-hour example, then the

center of avg interval = start of avg interval + dt/2 
                       = start of avg interval + 1800s
end of avg interval    = start of avg interval + dt
                       = start of avg interval + 3600s

[yantosca]

I think this is why the 1.5 is used: I know that INT( x + 0.5 ) is a method of rounding X to the nearest integer. I also think the 1.0 comes in because the result of the equation starts from 0 and you need to add 1 because Fortran arrays have indices starting at 1.

[yantosca]

Also tagging @SaptSinha and @Jourdan-He

[helpyuan]

Thank for your patience @yantosca .

I have understood the sampling in time, but I am still curious about how to sample in space. It is spatial interpolateion? Or another way?

[yantosca]

I believe it finds the grid box closest to the observation.

[helpyuan]

Thanks @yantosca !

[yantosca]

You are welcome @helpyuan. I will close out this issue now.