ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/OBS/obs_oper.F90

MODULE obs_oper
   !!======================================================================
   !!                       ***  MODULE obs_oper  ***
   !! Observation diagnostics: Observation operators for various observation
   !!                          types
   !!======================================================================

   !!----------------------------------------------------------------------
   !!   obs_pro_opt :    Compute the model counterpart of temperature and
   !!                    salinity observations from profiles
   !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly
   !!                    observations
   !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature
   !!                    observations
   !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity
   !!                    observations
   !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration
   !!                    observations
   !!
   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
   !!                    components of velocity from observations.
   !!----------------------------------------------------------------------

   !! * Modules used   
   USE par_kind, ONLY : &         ! Precision variables
      & wp
   USE in_out_manager             ! I/O manager
   USE obs_inter_sup              ! Interpolation support
   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
      & obs_int_h2d, &
      & obs_int_h2d_init
   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
      & obs_int_z1d,    &
      & obs_int_z1d_spl
   USE obs_const,  ONLY :     &
      & obfillflt		  ! Fillvalue   
   USE dom_oce,       ONLY : &
      & glamt, glamu, glamv, &
      & gphit, gphiu, gphiv
   USE lib_mpp,       ONLY : &
      & ctl_warn, ctl_stop

   IMPLICIT NONE

   !! * Routine accessibility
   PRIVATE

   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
      &   obs_seaice_opt, &
      &   obs_vel_opt     ! Compute the model counterpart of velocity profile data

   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id: obs_oper.F90 4245 2013-11-19 11:19:21Z cetlod $
   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &
      &                    kdailyavtypes )
      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_pro_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of profiles
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : Linearly interpolate to each observation point using 
      !!              the model values at the corners of the surrounding grid box.
      !!
      !!    First, a vertical profile of horizontally interpolated model
      !!    now temperatures is computed at the obs (lon, lat) point.
      !!    Several horizontal interpolation schemes are available:
      !!        - distance-weighted (great circle) (k2dint = 0)
      !!        - distance-weighted (small angle)  (k2dint = 1)
      !!        - bilinear (geographical grid)     (k2dint = 2)
      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
      !!
      !!    Next, the vertical temperature profile is interpolated to the
      !!    data depth points. Two vertical interpolation schemes are
      !!    available:
      !!        - linear       (k1dint = 0)
      !!        - Cubic spline (k1dint = 1)
      !!
      !!    For the cubic spline the 2nd derivative of the interpolating 
      !!    polynomial is computed before entering the vertical interpolation 
      !!    routine.
      !!
      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
      !!    a daily mean model temperature field. So, we first compute
      !!    the mean, then interpolate only at the end of the day.
      !!
      !!    Note: the in situ temperature observations must be converted
      !!    to potential temperature (the model variable) prior to
      !!    assimilation. 
      !!??????????????????????????????????????????????????????????????
      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
      !!??????????????????????????????????????????????????????????????
      !!
      !! ** Action  :
      !!
      !! History :
      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
      !!      ! 06-03 (G. Smith) NEMOVAR migration
      !!      ! 06-10 (A. Weaver) Cleanup
      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
      !!      ! 07-03 (K. Mogensen) General handling of profiles
      !!-----------------------------------------------------------------------
  
      !! * Modules used
      USE obs_profiles_def ! Definition of storage space for profile obs.

      IMPLICIT NONE

      !! * Arguments
      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
      INTEGER, INTENT(IN) :: kt        ! Time step
      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kpk
      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
                                       !   (kit000-1 = restart time)
      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                    
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
         & ptn,    &    ! Model temperature field
         & psn,    &    ! Model salinity field
         & ptmask       ! Land-sea mask
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
         & pgdept       ! Model array of depth levels
      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
         & kdailyavtypes! Types for daily averages
      !! * Local declarations
      INTEGER ::   ji
      INTEGER ::   jj
      INTEGER ::   jk
      INTEGER ::   jobs
      INTEGER ::   inrc
      INTEGER ::   ipro
      INTEGER ::   idayend
      INTEGER ::   ista
      INTEGER ::   iend
      INTEGER ::   iobs
      INTEGER, DIMENSION(imaxavtypes) :: &
         & idailyavtypes
      REAL(KIND=wp) :: zlam
      REAL(KIND=wp) :: zphi
      REAL(KIND=wp) :: zdaystp
      REAL(KIND=wp), DIMENSION(kpk) :: &
         & zobsmask, &
         & zobsk,    &
         & zobs2k
      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
         & zweig
      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         & zmask, &
         & zintt, &
         & zints, &
         & zinmt, &
         & zinms
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zglam, &
         & zgphi
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdi, &
         & igrdj

      !------------------------------------------------------------------------
      ! Local initialization 
      !------------------------------------------------------------------------
      ! ... Record and data counters
      inrc = kt - kit000 + 2
      ipro = prodatqc%npstp(inrc)
 
      ! Daily average types
      IF ( PRESENT(kdailyavtypes) ) THEN
         idailyavtypes(:) = kdailyavtypes(:)
      ELSE
         idailyavtypes(:) = -1
      ENDIF

      ! Initialize daily mean for first timestep
      idayend = MOD( kt - kit000 + 1, kdaystp )

      ! Added kt == 0 test to catch restart case 
      IF ( idayend == 1 .OR. kt == 0) THEN
         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpi
                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
               END DO
            END DO
         END DO
      ENDIF

      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               ! Increment the temperature field for computing daily mean
               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
                  &                        + ptn(ji,jj,jk)
               ! Increment the salinity field for computing daily mean
               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
                  &                        + psn(ji,jj,jk)
            END DO
         END DO
      END DO
   
      ! Compute the daily mean at the end of day
      zdaystp = 1.0 / REAL( kdaystp )
      IF ( idayend == 0 ) THEN
         DO jk = 1, jpk
            DO jj = 1, jpj
               DO ji = 1, jpi
                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
                     &                        * zdaystp
                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
                  &                           * zdaystp
               END DO
            END DO
         END DO
      ENDIF

      ! Get the data for interpolation
      ALLOCATE( &
         & igrdi(2,2,ipro),      &
         & igrdj(2,2,ipro),      &
         & zglam(2,2,ipro),      &
         & zgphi(2,2,ipro),      &
         & zmask(2,2,kpk,ipro),  &
         & zintt(2,2,kpk,ipro),  &
         & zints(2,2,kpk,ipro)   &
         & )

      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
         iobs = jobs - prodatqc%nprofup
         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
      END DO

      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )

      ! At the end of the day also get interpolated means
      IF ( idayend == 0 ) THEN

         ALLOCATE( &
            & zinmt(2,2,kpk,ipro),  &
            & zinms(2,2,kpk,ipro)   &
            & )

         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
            &                  prodatqc%vdmean(:,:,:,1), zinmt )
         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
            &                  prodatqc%vdmean(:,:,:,2), zinms )

      ENDIF

      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro

         iobs = jobs - prodatqc%nprofup

         IF ( kt /= prodatqc%mstp(jobs) ) THEN
            
            IF(lwp) THEN
               WRITE(numout,*)
               WRITE(numout,*) ' E R R O R : Observation',              &
                  &            ' time step is not consistent with the', &
                  &            ' model time step'
               WRITE(numout,*) ' ========='
               WRITE(numout,*)
               WRITE(numout,*) ' Record  = ', jobs,                    &
                  &            ' kt      = ', kt,                      &
                  &            ' mstp    = ', prodatqc%mstp(jobs), &
                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
            ENDIF
            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
         ENDIF
         
         zlam = prodatqc%rlam(jobs)
         zphi = prodatqc%rphi(jobs)
         
         ! Horizontal weights and vertical mask

         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN

            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
               &                   zmask(:,:,:,iobs), zweig, zobsmask )

         ENDIF

         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN

            zobsk(:) = obfillflt

	    IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN

               IF ( idayend == 0 )  THEN
                  
                  ! Daily averaged moored buoy (MRB) data
                  
                  CALL obs_int_h2d( kpk, kpk,      &
                     &              zweig, zinmt(:,:,:,iobs), zobsk )
                  
                  
               ELSE
               
                  CALL ctl_stop( ' A nonzero' //     &
                     &           ' number of profile T BUOY data should' // &
                     &           ' only occur at the end of a given day' )

               ENDIF
	       
            ELSE 
               
               ! Point data

               CALL obs_int_h2d( kpk, kpk,      &
                  &              zweig, zintt(:,:,:,iobs), zobsk )

            ENDIF

            !-------------------------------------------------------------
            ! Compute vertical second-derivative of the interpolating 
            ! polynomial at obs points
            !-------------------------------------------------------------
            
            IF ( k1dint == 1 ) THEN
               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
                  &                  pgdept, zobsmask )
            ENDIF
            
            !-----------------------------------------------------------------
            !  Vertical interpolation to the observation point
            !-----------------------------------------------------------------
            ista = prodatqc%npvsta(jobs,1)
            iend = prodatqc%npvend(jobs,1)
            CALL obs_int_z1d( kpk,                &
               & prodatqc%var(1)%mvk(ista:iend),  &
               & k1dint, iend - ista + 1,         &
               & prodatqc%var(1)%vdep(ista:iend), &
               & zobsk, zobs2k,                   &
               & prodatqc%var(1)%vmod(ista:iend), &
               & pgdept, zobsmask )

         ENDIF

         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN

            zobsk(:) = obfillflt

            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN

               IF ( idayend == 0 )  THEN

                  ! Daily averaged moored buoy (MRB) data
                  
                  CALL obs_int_h2d( kpk, kpk,      &
                     &              zweig, zinms(:,:,:,iobs), zobsk )
                  
               ELSE

                  CALL ctl_stop( ' A nonzero' //     &
                     &           ' number of profile S BUOY data should' // &
                     &           ' only occur at the end of a given day' )

               ENDIF

            ELSE
               
               ! Point data

               CALL obs_int_h2d( kpk, kpk,      &
                  &              zweig, zints(:,:,:,iobs), zobsk )

            ENDIF


            !-------------------------------------------------------------
            ! Compute vertical second-derivative of the interpolating 
            ! polynomial at obs points
            !-------------------------------------------------------------
            
            IF ( k1dint == 1 ) THEN
               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
                  &                  pgdept, zobsmask )
            ENDIF
            
            !----------------------------------------------------------------
            !  Vertical interpolation to the observation point
            !----------------------------------------------------------------
            ista = prodatqc%npvsta(jobs,2)
            iend = prodatqc%npvend(jobs,2)
            CALL obs_int_z1d( kpk, &
               & prodatqc%var(2)%mvk(ista:iend),&
               & k1dint, iend - ista + 1, &
               & prodatqc%var(2)%vdep(ista:iend),&
               & zobsk, zobs2k, &
               & prodatqc%var(2)%vmod(ista:iend),&
               & pgdept, zobsmask )

         ENDIF

      END DO
 
      ! Deallocate the data for interpolation
      DEALLOCATE( &
         & igrdi, &
         & igrdj, &
         & zglam, &
         & zgphi, &
         & zmask, &
         & zintt, &
         & zints  &
         & )
      ! At the end of the day also get interpolated means
      IF ( idayend == 0 ) THEN
         DEALLOCATE( &
            & zinmt,  &
            & zinms   &
            & )
      ENDIF

      prodatqc%nprofup = prodatqc%nprofup + ipro 
      
   END SUBROUTINE obs_pro_opt

   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
      &                    psshn, psshmask, k2dint )
      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_sla_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of sea level anomaly
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : Linearly interpolate to each observation point using 
      !!              the model values at the corners of the surrounding grid box.
      !!
      !!    The now model SSH is first computed at the obs (lon, lat) point.
      !!
      !!    Several horizontal interpolation schemes are available:
      !!        - distance-weighted (great circle) (k2dint = 0)
      !!        - distance-weighted (small angle)  (k2dint = 1)
      !!        - bilinear (geographical grid)     (k2dint = 2)
      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
      !!  
      !!    The sea level anomaly at the observation points is then computed 
      !!    by removing a mean dynamic topography (defined at the obs. point).
      !!
      !! ** Action  :
      !!
      !! History :
      !!      ! 07-03 (A. Weaver)
      !!-----------------------------------------------------------------------
  
      !! * Modules used
      USE obs_surf_def  ! Definition of storage space for surface observations

      IMPLICIT NONE

      !! * Arguments
      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
      INTEGER, INTENT(IN) :: kt      ! Time step
      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
                                      !   (kit000-1 = restart time)
      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
         & psshn,  &    ! Model SSH field
         & psshmask     ! Land-sea mask
         
      !! * Local declarations
      INTEGER :: ji
      INTEGER :: jj
      INTEGER :: jobs
      INTEGER :: inrc
      INTEGER :: isla
      INTEGER :: iobs
      REAL(KIND=wp) :: zlam
      REAL(KIND=wp) :: zphi
      REAL(KIND=wp) :: zext(1), zobsmask(1)
      REAL(kind=wp), DIMENSION(2,2,1) :: &
         & zweig
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zmask, &
         & zsshl, &
         & zglam, &
         & zgphi
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdi, &
         & igrdj

      !------------------------------------------------------------------------
      ! Local initialization 
      !------------------------------------------------------------------------
      ! ... Record and data counters
      inrc = kt - kit000 + 2
      isla = sladatqc%nsstp(inrc)

      ! Get the data for interpolation

      ALLOCATE( &
         & igrdi(2,2,isla), &
         & igrdj(2,2,isla), &
         & zglam(2,2,isla), &
         & zgphi(2,2,isla), &
         & zmask(2,2,isla), &
         & zsshl(2,2,isla)  &
         & )
      
      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
         iobs = jobs - sladatqc%nsurfup
         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
         igrdj(1,2,iobs) = sladatqc%mj(jobs)
         igrdi(2,1,iobs) = sladatqc%mi(jobs)
         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
         igrdi(2,2,iobs) = sladatqc%mi(jobs)
         igrdj(2,2,iobs) = sladatqc%mj(jobs)
      END DO

      CALL obs_int_comm_2d( 2, 2, isla, &
         &                  igrdi, igrdj, glamt, zglam )
      CALL obs_int_comm_2d( 2, 2, isla, &
         &                  igrdi, igrdj, gphit, zgphi )
      CALL obs_int_comm_2d( 2, 2, isla, &
         &                  igrdi, igrdj, psshmask, zmask )
      CALL obs_int_comm_2d( 2, 2, isla, &
         &                  igrdi, igrdj, psshn, zsshl )

      ! Loop over observations

      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla

         iobs = jobs - sladatqc%nsurfup

         IF ( kt /= sladatqc%mstp(jobs) ) THEN
            
            IF(lwp) THEN
               WRITE(numout,*)
               WRITE(numout,*) ' E R R O R : Observation',              &
                  &            ' time step is not consistent with the', &
                  &            ' model time step'
               WRITE(numout,*) ' ========='
               WRITE(numout,*)
               WRITE(numout,*) ' Record  = ', jobs,                &
                  &            ' kt      = ', kt,                  &
                  &            ' mstp    = ', sladatqc%mstp(jobs), &
                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
            ENDIF
            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
            
         ENDIF
         
         zlam = sladatqc%rlam(jobs)
         zphi = sladatqc%rphi(jobs)

         ! Get weights to interpolate the model SSH to the observation point
         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
            &                   zmask(:,:,iobs), zweig, zobsmask )
         

         ! Interpolate the model SSH to the observation point
         CALL obs_int_h2d( 1, 1,      &
            &              zweig, zsshl(:,:,iobs),  zext )
         
         sladatqc%rext(jobs,1) = zext(1)
         ! ... Remove the MDT at the observation point
         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)

      END DO

      ! Deallocate the data for interpolation
      DEALLOCATE( &
         & igrdi, &
         & igrdj, &
         & zglam, &
         & zgphi, &
         & zmask, &
         & zsshl  &
         & )

      sladatqc%nsurfup = sladatqc%nsurfup + isla

   END SUBROUTINE obs_sla_opt

   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
      &                    psstn, psstmask, k2dint, ld_nightav )
      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_sst_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of surface temperature
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : Linearly interpolate to each observation point using 
      !!              the model values at the corners of the surrounding grid box.
      !!
      !!    The now model SST is first computed at the obs (lon, lat) point.
      !!
      !!    Several horizontal interpolation schemes are available:
      !!        - distance-weighted (great circle) (k2dint = 0)
      !!        - distance-weighted (small angle)  (k2dint = 1)
      !!        - bilinear (geographical grid)     (k2dint = 2)
      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
      !!
      !!
      !! ** Action  :
      !!
      !! History :
      !!        !  07-07  (S. Ricci ) : Original
      !!      
      !!-----------------------------------------------------------------------

      !! * Modules used
      USE obs_surf_def  ! Definition of storage space for surface observations
      USE sbcdcy

      IMPLICIT NONE

      !! * Arguments
      TYPE(obs_surf), INTENT(INOUT) :: &
         & sstdatqc     ! Subset of surface data not failing screening
      INTEGER, INTENT(IN) :: kt        ! Time step
      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
                                       !   (kit000-1 = restart time)
      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day  
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
         & psstn,  &    ! Model SST field
         & psstmask     ! Land-sea mask

      !! * Local declarations
      INTEGER :: ji
      INTEGER :: jj
      INTEGER :: jobs
      INTEGER :: inrc
      INTEGER :: isst
      INTEGER :: iobs
      INTEGER :: idayend
      REAL(KIND=wp) :: zlam
      REAL(KIND=wp) :: zphi
      REAL(KIND=wp) :: zext(1), zobsmask(1)
      REAL(KIND=wp) :: zdaystp
      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
         & icount_sstnight,      &
         & imask_night
      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
         & zintmp, &
         & zouttmp, & 
         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
      REAL(kind=wp), DIMENSION(2,2,1) :: &
         & zweig
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zmask, &
         & zsstl, &
         & zsstm, &
         & zglam, &
         & zgphi
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdi, &
         & igrdj
      LOGICAL, INTENT(IN) :: ld_nightav

      !-----------------------------------------------------------------------
      ! Local initialization 
      !-----------------------------------------------------------------------
      ! ... Record and data counters
      inrc = kt - kit000 + 2
      isst = sstdatqc%nsstp(inrc)

      IF ( ld_nightav ) THEN

      ! Initialize array for night mean

      IF ( kt .EQ. 0 ) THEN
         ALLOCATE ( icount_sstnight(kpi,kpj) )
         ALLOCATE ( imask_night(kpi,kpj) )
         ALLOCATE ( zintmp(kpi,kpj) )
         ALLOCATE ( zouttmp(kpi,kpj) )
         ALLOCATE ( zmeanday(kpi,kpj) )
         nday_qsr = -1   ! initialisation flag for nbc_dcy
      ENDIF

      ! Initialize daily mean for first timestep
      idayend = MOD( kt - kit000 + 1, kdaystp )

      ! Added kt == 0 test to catch restart case 
      IF ( idayend == 1 .OR. kt == 0) THEN
         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
         DO jj = 1, jpj
            DO ji = 1, jpi
               sstdatqc%vdmean(ji,jj) = 0.0
               zmeanday(ji,jj) = 0.0
               icount_sstnight(ji,jj) = 0
            END DO
         END DO
      ENDIF

      zintmp(:,:) = 0.0
      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
      imask_night(:,:) = INT( zouttmp(:,:) )

      DO jj = 1, jpj
         DO ji = 1, jpi
            ! Increment the temperature field for computing night mean and counter
            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  &
                   &                        + psstn(ji,jj)*imask_night(ji,jj)
            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj)
            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
         END DO
      END DO
   
      ! Compute the daily mean at the end of day

      zdaystp = 1.0 / REAL( kdaystp )

      IF ( idayend == 0 ) THEN 
         DO jj = 1, jpj
            DO ji = 1, jpi
               ! Test if "no night" point
               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
                    &                        / icount_sstnight(ji,jj) 
               ELSE
                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
               ENDIF
            END DO
         END DO
      ENDIF

      ENDIF

      ! Get the data for interpolation
      
      ALLOCATE( &
         & igrdi(2,2,isst), &
         & igrdj(2,2,isst), &
         & zglam(2,2,isst), &
         & zgphi(2,2,isst), &
         & zmask(2,2,isst), &
         & zsstl(2,2,isst)  &
         & )
      
      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
         iobs = jobs - sstdatqc%nsurfup
         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
      END DO
      
      CALL obs_int_comm_2d( 2, 2, isst, &
         &                  igrdi, igrdj, glamt, zglam )
      CALL obs_int_comm_2d( 2, 2, isst, &
         &                  igrdi, igrdj, gphit, zgphi )
      CALL obs_int_comm_2d( 2, 2, isst, &
         &                  igrdi, igrdj, psstmask, zmask )
      CALL obs_int_comm_2d( 2, 2, isst, &
         &                  igrdi, igrdj, psstn, zsstl )

      ! At the end of the day get interpolated means
      IF ( idayend == 0 .AND. ld_nightav ) THEN

         ALLOCATE( &
            & zsstm(2,2,isst)  &
            & )

         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
            &               sstdatqc%vdmean(:,:), zsstm )

      ENDIF

      ! Loop over observations

      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
         
         iobs = jobs - sstdatqc%nsurfup
         
         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
            
            IF(lwp) THEN
               WRITE(numout,*)
               WRITE(numout,*) ' E R R O R : Observation',              &
                  &            ' time step is not consistent with the', &
                  &            ' model time step'
               WRITE(numout,*) ' ========='
               WRITE(numout,*)
               WRITE(numout,*) ' Record  = ', jobs,                &
                  &            ' kt      = ', kt,                  &
                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
            ENDIF
            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
            
         ENDIF
         
         zlam = sstdatqc%rlam(jobs)
         zphi = sstdatqc%rphi(jobs)
         
         ! Get weights to interpolate the model SST to the observation point
         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
            &                   zmask(:,:,iobs), zweig, zobsmask )
            
         ! Interpolate the model SST to the observation point 

         IF ( ld_nightav ) THEN

           IF ( idayend == 0 )  THEN
               ! Daily averaged/diurnal cycle of SST  data
               CALL obs_int_h2d( 1, 1,      & 
                     &              zweig, zsstm(:,:,iobs), zext )
            ELSE 
               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
                     &           ' number of night SST data should' // &
                     &           ' only occur at the end of a given day' )
            ENDIF

         ELSE

            CALL obs_int_h2d( 1, 1,      &
            &              zweig, zsstl(:,:,iobs),  zext )

         ENDIF
         sstdatqc%rmod(jobs,1) = zext(1)
         
      END DO
      
      ! Deallocate the data for interpolation
      DEALLOCATE( &
         & igrdi, &
         & igrdj, &
         & zglam, &
         & zgphi, &
         & zmask, &
         & zsstl  &
         & )

      ! At the end of the day also get interpolated means
      IF ( idayend == 0 .AND. ld_nightav ) THEN
         DEALLOCATE( &
            & zsstm  &
            & )
      ENDIF
      
      sstdatqc%nsurfup = sstdatqc%nsurfup + isst

   END SUBROUTINE obs_sst_opt

   SUBROUTINE obs_sss_opt
      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_sss_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of sea surface salinity
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : 
      !!
      !! ** Action  :
      !!
      !! History :
      !!      ! ??-?? 
      !!-----------------------------------------------------------------------

      IMPLICIT NONE

   END SUBROUTINE obs_sss_opt

   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
      &                    pseaicen, pseaicemask, k2dint )

      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_seaice_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of surface temperature
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : Linearly interpolate to each observation point using 
      !!              the model values at the corners of the surrounding grid box.
      !!
      !!    The now model sea ice is first computed at the obs (lon, lat) point.
      !!
      !!    Several horizontal interpolation schemes are available:
      !!        - distance-weighted (great circle) (k2dint = 0)
      !!        - distance-weighted (small angle)  (k2dint = 1)
      !!        - bilinear (geographical grid)     (k2dint = 2)
      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
      !!
      !!
      !! ** Action  :
      !!
      !! History :
      !!        !  07-07  (S. Ricci ) : Original
      !!      
      !!-----------------------------------------------------------------------

      !! * Modules used
      USE obs_surf_def  ! Definition of storage space for surface observations

      IMPLICIT NONE

      !! * Arguments
      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
      INTEGER, INTENT(IN) :: kt       ! Time step
      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step 
                                      !   (kit000-1 = restart time)
      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
         & pseaicen,  &    ! Model sea ice field
         & pseaicemask     ! Land-sea mask
         
      !! * Local declarations
      INTEGER :: ji
      INTEGER :: jj
      INTEGER :: jobs
      INTEGER :: inrc
      INTEGER :: iseaice
      INTEGER :: iobs
       
      REAL(KIND=wp) :: zlam
      REAL(KIND=wp) :: zphi
      REAL(KIND=wp) :: zext(1), zobsmask(1)
      REAL(kind=wp), DIMENSION(2,2,1) :: &
         & zweig
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zmask, &
         & zseaicel, &
         & zglam, &
         & zgphi
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdi, &
         & igrdj

      !------------------------------------------------------------------------
      ! Local initialization 
      !------------------------------------------------------------------------
      ! ... Record and data counters
      inrc = kt - kit000 + 2
      iseaice = seaicedatqc%nsstp(inrc)

      ! Get the data for interpolation
      
      ALLOCATE( &
         & igrdi(2,2,iseaice), &
         & igrdj(2,2,iseaice), &
         & zglam(2,2,iseaice), &
         & zgphi(2,2,iseaice), &
         & zmask(2,2,iseaice), &
         & zseaicel(2,2,iseaice)  &
         & )
      
      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
         iobs = jobs - seaicedatqc%nsurfup
         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
      END DO
      
      CALL obs_int_comm_2d( 2, 2, iseaice, &
         &                  igrdi, igrdj, glamt, zglam )
      CALL obs_int_comm_2d( 2, 2, iseaice, &
         &                  igrdi, igrdj, gphit, zgphi )
      CALL obs_int_comm_2d( 2, 2, iseaice, &
         &                  igrdi, igrdj, pseaicemask, zmask )
      CALL obs_int_comm_2d( 2, 2, iseaice, &
         &                  igrdi, igrdj, pseaicen, zseaicel )
      
      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
         
         iobs = jobs - seaicedatqc%nsurfup
         
         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
            
            IF(lwp) THEN
               WRITE(numout,*)
               WRITE(numout,*) ' E R R O R : Observation',              &
                  &            ' time step is not consistent with the', &
                  &            ' model time step'
               WRITE(numout,*) ' ========='
               WRITE(numout,*)
               WRITE(numout,*) ' Record  = ', jobs,                &
                  &            ' kt      = ', kt,                  &
                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
            ENDIF
            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
            
         ENDIF
         
         zlam = seaicedatqc%rlam(jobs)
         zphi = seaicedatqc%rphi(jobs)
         
         ! Get weights to interpolate the model sea ice to the observation point
         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
            &                   zmask(:,:,iobs), zweig, zobsmask )
         
         ! ... Interpolate the model sea ice to the observation point
         CALL obs_int_h2d( 1, 1,      &
            &              zweig, zseaicel(:,:,iobs),  zext )
         
         seaicedatqc%rmod(jobs,1) = zext(1)
         
      END DO
      
      ! Deallocate the data for interpolation
      DEALLOCATE( &
         & igrdi,    &
         & igrdj,    &
         & zglam,    &
         & zgphi,    &
         & zmask,    &
         & zseaicel  &
         & )
      
      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice

   END SUBROUTINE obs_seaice_opt

   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
      &                    ld_dailyav )
      !!-----------------------------------------------------------------------
      !!
      !!                     ***  ROUTINE obs_vel_opt  ***
      !!
      !! ** Purpose : Compute the model counterpart of velocity profile
      !!              data by interpolating from the model grid to the 
      !!              observation point.
      !!
      !! ** Method  : Linearly interpolate zonal and meridional components of velocity 
      !!              to each observation point using the model values at the corners of 
      !!              the surrounding grid box. The model velocity components are on a 
      !!              staggered C- grid.
      !!
      !!    For velocity data from the TAO array, the model equivalent is
      !!    a daily mean velocity field. So, we first compute
      !!    the mean, then interpolate only at the end of the day.
      !!
      !! ** Action  :
      !!
      !! History :
      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
      !!-----------------------------------------------------------------------
    
      !! * Modules used
      USE obs_profiles_def ! Definition of storage space for profile obs.

      IMPLICIT NONE

      !! * Arguments
      TYPE(obs_prof), INTENT(INOUT) :: &
         & prodatqc        ! Subset of profile data not failing screening
      INTEGER, INTENT(IN) :: kt        ! Time step
      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
      INTEGER, INTENT(IN) :: kpj
      INTEGER, INTENT(IN) :: kpk 
      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
                                       !   (kit000-1 = restart time)
      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                    
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
         & pun,    &    ! Model zonal component of velocity
         & pvn,    &    ! Model meridional component of velocity
         & pumask, &    ! Land-sea mask
         & pvmask       ! Land-sea mask
      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
         & pgdept       ! Model array of depth levels
      LOGICAL, INTENT(IN) :: ld_dailyav
         
      !! * Local declarations
      INTEGER :: ji
      INTEGER :: jj
      INTEGER :: jk
      INTEGER :: jobs
      INTEGER :: inrc
      INTEGER :: ipro
      INTEGER :: idayend
      INTEGER :: ista
      INTEGER :: iend
      INTEGER :: iobs
      INTEGER, DIMENSION(imaxavtypes) :: &
         & idailyavtypes
      REAL(KIND=wp) :: zlam
      REAL(KIND=wp) :: zphi
      REAL(KIND=wp) :: zdaystp
      REAL(KIND=wp), DIMENSION(kpk) :: &
         & zobsmasku, &
         & zobsmaskv, &
         & zobsmask,  &
         & zobsk,     &
         & zobs2k
      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
         & zweigu,zweigv
      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
         & zumask, zvmask, &
         & zintu, &
         & zintv, &
         & zinmu, &
         & zinmv
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zglamu, zglamv, &
         & zgphiu, zgphiv
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdiu, &
         & igrdju, &
         & igrdiv, &
         & igrdjv

      !------------------------------------------------------------------------
      ! Local initialization 
      !------------------------------------------------------------------------
      ! ... Record and data counters
      inrc = kt - kit000 + 2
      ipro = prodatqc%npstp(inrc)

      ! Initialize daily mean for first timestep
      idayend = MOD( kt - kit000 + 1, kdaystp )

      ! Added kt == 0 test to catch restart case 
      IF ( idayend == 1 .OR. kt == 0) THEN
         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
         prodatqc%vdmean(:,:,:,1) = 0.0
         prodatqc%vdmean(:,:,:,2) = 0.0
      ENDIF

      ! Increment the zonal velocity field for computing daily mean
      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
      ! Increment the meridional velocity field for computing daily mean
      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
   
      ! Compute the daily mean at the end of day
      zdaystp = 1.0 / REAL( kdaystp )
      IF ( idayend == 0 ) THEN
         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
      ENDIF

      ! Get the data for interpolation
      ALLOCATE( &
         & igrdiu(2,2,ipro),      &
         & igrdju(2,2,ipro),      &
         & igrdiv(2,2,ipro),      &
         & igrdjv(2,2,ipro),      &
         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
         & zintu(2,2,kpk,ipro),  &
         & zintv(2,2,kpk,ipro)   &
         & )

      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
         iobs = jobs - prodatqc%nprofup
         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
      END DO

      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )

      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )

      ! At the end of the day also get interpolated means
      IF ( idayend == 0 ) THEN

         ALLOCATE( &
            & zinmu(2,2,kpk,ipro),  &
            & zinmv(2,2,kpk,ipro)   &
            & )

         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
            &                  prodatqc%vdmean(:,:,:,1), zinmu )
         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
            &                  prodatqc%vdmean(:,:,:,2), zinmv )

      ENDIF

! loop over observations

      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro

         iobs = jobs - prodatqc%nprofup

         IF ( kt /= prodatqc%mstp(jobs) ) THEN
            
            IF(lwp) THEN
               WRITE(numout,*)
               WRITE(numout,*) ' E R R O R : Observation',              &
                  &            ' time step is not consistent with the', &
                  &            ' model time step'
               WRITE(numout,*) ' ========='
               WRITE(numout,*)
               WRITE(numout,*) ' Record  = ', jobs,                    &
                  &            ' kt      = ', kt,                      &
                  &            ' mstp    = ', prodatqc%mstp(jobs), &
                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
            ENDIF
            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
         ENDIF
         
         zlam = prodatqc%rlam(jobs)
         zphi = prodatqc%rphi(jobs)

         ! Initialize observation masks

         zobsmasku(:) = 0.0
         zobsmaskv(:) = 0.0
         
         ! Horizontal weights and vertical mask

         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN

            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )

         ENDIF

         
         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN

            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )

         ENDIF

         ! Ensure that the vertical mask on u and v are consistent.

         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )

         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN

            zobsk(:) = obfillflt

	    IF ( ld_dailyav ) THEN

               IF ( idayend == 0 )  THEN
                  
                  ! Daily averaged data
                  
                  CALL obs_int_h2d( kpk, kpk,      &
                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
                  
                  
               ELSE
               
                  CALL ctl_stop( ' A nonzero' //     &
                     &           ' number of U profile data should' // &
                     &           ' only occur at the end of a given day' )

               ENDIF
	       
            ELSE 
               
               ! Point data

               CALL obs_int_h2d( kpk, kpk,      &
                  &              zweigu, zintu(:,:,:,iobs), zobsk )

            ENDIF

            !-------------------------------------------------------------
            ! Compute vertical second-derivative of the interpolating 
            ! polynomial at obs points
            !-------------------------------------------------------------
            
            IF ( k1dint == 1 ) THEN
               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
                  &                  pgdept, zobsmask )
            ENDIF
            
            !-----------------------------------------------------------------
            !  Vertical interpolation to the observation point
            !-----------------------------------------------------------------
            ista = prodatqc%npvsta(jobs,1)
            iend = prodatqc%npvend(jobs,1)
            CALL obs_int_z1d( kpk,                &
               & prodatqc%var(1)%mvk(ista:iend),  &
               & k1dint, iend - ista + 1,         &
               & prodatqc%var(1)%vdep(ista:iend), &
               & zobsk, zobs2k,                   &
               & prodatqc%var(1)%vmod(ista:iend), &
               & pgdept, zobsmask )

         ENDIF

         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN

            zobsk(:) = obfillflt

            IF ( ld_dailyav ) THEN

               IF ( idayend == 0 )  THEN

                  ! Daily averaged data
                  
                  CALL obs_int_h2d( kpk, kpk,      &
                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
                  
               ELSE

                  CALL ctl_stop( ' A nonzero' //     &
                     &           ' number of V profile data should' // &
                     &           ' only occur at the end of a given day' )

               ENDIF

            ELSE
               
               ! Point data

               CALL obs_int_h2d( kpk, kpk,      &
                  &              zweigv, zintv(:,:,:,iobs), zobsk )

            ENDIF


            !-------------------------------------------------------------
            ! Compute vertical second-derivative of the interpolating 
            ! polynomial at obs points
            !-------------------------------------------------------------
            
            IF ( k1dint == 1 ) THEN
               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
                  &                  pgdept, zobsmask )
            ENDIF
            
            !----------------------------------------------------------------
            !  Vertical interpolation to the observation point
            !----------------------------------------------------------------
            ista = prodatqc%npvsta(jobs,2)
            iend = prodatqc%npvend(jobs,2)
            CALL obs_int_z1d( kpk, &
               & prodatqc%var(2)%mvk(ista:iend),&
               & k1dint, iend - ista + 1, &
               & prodatqc%var(2)%vdep(ista:iend),&
               & zobsk, zobs2k, &
               & prodatqc%var(2)%vmod(ista:iend),&
               & pgdept, zobsmask )

         ENDIF

      END DO
 
      ! Deallocate the data for interpolation
      DEALLOCATE( &
         & igrdiu, &
         & igrdju, &
         & igrdiv, &
         & igrdjv, &
         & zglamu, zglamv, &
         & zgphiu, zgphiv, &
         & zumask, zvmask, &
         & zintu, &
         & zintv  &
         & )
      ! At the end of the day also get interpolated means
      IF ( idayend == 0 ) THEN
         DEALLOCATE( &
            & zinmu,  &
            & zinmv   &
            & )
      ENDIF

      prodatqc%nprofup = prodatqc%nprofup + ipro 
      
   END SUBROUTINE obs_vel_opt

END MODULE obs_oper