MODULE obs_readmdt
   !!======================================================================
   !!                       ***  MODULE obs_readmdt  ***
   !! Observation diagnostics: Read the MDT for SLA data (skeleton for now)
   !!======================================================================

   !!----------------------------------------------------------------------
   !!   obs_rea_mdt : Driver for reading MDT
   !!----------------------------------------------------------------------

   !! * Modules used   
   USE par_kind, ONLY : &       ! Precision variables
      & wp, &
      & dp, &
      & sp
   USE par_oce, ONLY : &        ! Domain parameters
      & jpi, &
      & jpj, &
      & jpim1
   USE in_out_manager, ONLY : & ! I/O manager
      & lwp,    &
      & numout 
   USE obs_surf_def             ! Surface observation definitions
   USE dom_oce, ONLY : &        ! Domain variables
      & tmask, &
      & tmask_i, &
      & e1t,   &
      & e2t,   &
      & gphit, &
      & glamt
   USE obs_const, ONLY : &
      & obfillflt              ! Fillvalue
   USE oce, ONLY : &           ! Model variables
      & sshn
   USE obs_inter_sup           ! Interpolation support routines
   USE obs_inter_h2d           ! 2D interpolation
   USE obs_utils               ! Various observation tools
   USE lib_mpp, only: &        ! MPP routines
      & lk_mpp, &
      & mpp_sum
   USE iom_nf90   
   USE netcdf                   ! NetCDF library

   IMPLICIT NONE

   !! * Routine accessibility
   PRIVATE

   INTEGER, PUBLIC :: nmsshc = 1        ! MDT correction scheme
   REAL(wp), PUBLIC :: mdtcorr = 1.61   ! User specified MDT correction
   REAL(wp), PUBLIC :: mdtcutoff = 65.0  ! MDT cutoff for computed correction
   PUBLIC obs_rea_mdt     ! Read the MDT
   PUBLIC obs_offset_mdt  ! Remove the offset between the model MDT and the 
                          ! used one

CONTAINS

   SUBROUTINE obs_rea_mdt( kslano, sladata, k2dint )
      !!---------------------------------------------------------------------
      !!
      !!                   *** ROUTINE obs_rea_mdt ***
      !!
      !! ** Purpose : Read from file the MDT data (skeleton)
      !!
      !! ** Method  : 
      !!
      !! ** Action  : 
      !!
      !! References :
      !!
      !! History :  
      !!      ! :  2007-03 (K. Mogensen) Initial skeleton version
      !!      ! :  2007-04 (E. Remy) migration and improvement from OPAVAR
      !!----------------------------------------------------------------------
      !! * Modules used
      USE iom

      !! * Arguments
      INTEGER, INTENT(IN) :: kslano          ! Number of SLA Products
      TYPE(obs_surf), DIMENSION(kslano), INTENT(INOUT) :: &
         & sladata       ! SLA data
      INTEGER, INTENT(IN) :: k2dint

      !! * Local declarations

      CHARACTER(LEN=12), PARAMETER :: &
         & cpname = 'obs_rea_mdt'
      CHARACTER(LEN=20), PARAMETER :: &
         & mdtname = 'slaReferenceLevel.nc'

      INTEGER :: jslano      ! Data set loop variable
      INTEGER :: jobs        ! Obs loop variable
      INTEGER :: jpimdt      ! Number of grid point in latitude for the MDT
      INTEGER :: jpjmdt      ! Number of grid point in longitude for the MDT
      INTEGER :: iico        ! Grid point indicies
      INTEGER :: ijco 
      INTEGER :: i_nx_id     ! Index to read the NetCDF file
      INTEGER :: i_ny_id     ! 
      INTEGER :: i_file_id   ! 
      INTEGER :: i_var_id
      INTEGER :: i_stat

      REAL(wp), DIMENSION(jpi,jpj) :: & 
         & z_mdt,       &  ! Array to store the MDT values
         & mdtmask         ! Array to store the mask for the MDT
      REAL(wp), DIMENSION(1) :: &
         & zext, &
         & zobsmask
      REAL(wp), DIMENSION(2,2,1) :: &
         & zweig
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
         & zmask, &
         & zmdtl, &
         & zglam, &
         & zgphi
         
      REAL(wp) :: zlam
      REAL(wp) :: zphi
      REAL(wp) :: zfill
      REAL(sp) :: zinfill
      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
         & igrdi, &
         & igrdj
      INTEGER :: nummdt

      IF(lwp)WRITE(numout,*) 
      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : '
      IF(lwp)WRITE(numout,*) ' ------------- '
      IF(lwp)WRITE(numout,*) '   Read MDT for referencing altimeter', &
         &                   '   anomalies'

      ! Open the file
      
      CALL iom_open( mdtname, nummdt )
      
      ! Get the MDT data

      CALL iom_get( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 )

      ! Close the file

      CALL iom_close(nummdt)     
      
      ! Read in the fill value
      zinfill = 0.0
      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
      zfill = zinfill
      i_stat = nf90_close( nummdt )

! setup mask based on tmask and MDT mask
! set mask to 0 where the MDT is set to fillvalue

      WHERE(z_mdt(:,:) /= zfill)
         mdtmask(:,:)=tmask(:,:,1)
      ELSEWHERE
         mdtmask(:,:)=0
      END WHERE

      ! Remove the offset between the MDT used with the sla and the model MDT

      IF ( nmsshc == 1 .OR. nmsshc == 2 ) CALL obs_offset_mdt( z_mdt, zfill )

      ! Intepolate the MDT already on the model grid at the observation point
  
      DO jslano = 1, kslano

         ALLOCATE( &
            & igrdi(2,2,sladata(jslano)%nsurf), &
            & igrdj(2,2,sladata(jslano)%nsurf), &
            & zglam(2,2,sladata(jslano)%nsurf), &
            & zgphi(2,2,sladata(jslano)%nsurf), &
            & zmask(2,2,sladata(jslano)%nsurf), &
            & zmdtl(2,2,sladata(jslano)%nsurf)  &
            & )
         
         DO jobs = 1, sladata(jslano)%nsurf

            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)

         END DO

         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
            &                  igrdi, igrdj, glamt, zglam )
         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
            &                  igrdi, igrdj, gphit, zgphi )
         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
            &                  igrdi, igrdj, mdtmask, zmask )
         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
            &                  igrdi, igrdj, z_mdt, zmdtl )

         DO jobs = 1, sladata(jslano)%nsurf
            
            zlam = sladata(jslano)%rlam(jobs)
            zphi = sladata(jslano)%rphi(jobs)

            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
               &                   zmask(:,:,jobs), zweig, zobsmask )
            
            CALL obs_int_h2d( 1, 1,      &
                   &              zweig, zmdtl(:,:,jobs),  zext )
 
            sladata(jslano)%rext(jobs,2) = zext(1)

! mark any masked data with a QC flag
            IF ( zobsmask(1) == 0 ) sladata(jslano)%nqc(jobs) = 11

         END DO
         
         DEALLOCATE( &
            & igrdi, &
            & igrdj, &
            & zglam, &
            & zgphi, &
            & zmask, &
            & zmdtl  &
            & )

      END DO

   END SUBROUTINE obs_rea_mdt

   SUBROUTINE obs_offset_mdt( mdt, zfill )

      !!---------------------------------------------------------------------
      !!
      !!                   *** ROUTINE obs_offset_mdt ***
      !!
      !! ** Purpose : Compute a correction term for the MDT on the model grid
      !!             !!!!! IF it is on the model grid
      !!
      !! ** Method  : Compute the mean difference between the model and the 
      !!              used MDT and remove the offset. 
      !!
      !! ** Action  : 
      !!
      !! References :
      !!
      !! History :  
      !!      ! :  2007-04 (E. Remy) migration from OPAVAR
      !!----------------------------------------------------------------------
      !! * Modules used

      !! * Arguments
      REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
         & mdt           ! MDT used on the model grid
      REAL(wp), INTENT(IN) :: zfill 

      !! * Local declarations 
      REAL(wp) :: zdxdy
      REAL(wp) :: zarea
      REAL(wp) :: zeta1
      REAL(wp) :: zeta2
      REAL(wp) :: zcorr_mdt  
      REAL(wp) :: zcorr_bcketa
      REAL(wp) :: zcorr
      REAL(wp), DIMENSION(jpi,jpj) :: zpromsk
      INTEGER :: jj
      INTEGER :: ji
      CHARACTER(LEN=14), PARAMETER :: &
         & cpname = 'obs_offset_mdt'
   
      !  Initialize the local mask, for domain projection 
      !  Also exclude mdt points which are set to missing data

      DO ji = 1, jpi
        DO jj = 1, jpj
           zpromsk(ji,jj) = tmask_i(ji,jj)
           IF (    ( gphit(ji,jj) .GT.  mdtcutoff ) &
              &.OR.( gphit(ji,jj) .LT. -mdtcutoff ) &
              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
              &        zpromsk(ji,jj) = 0.0
        END DO
      END DO 

      ! Compute MSSH mean over [0,360] x [-mdtcutoff,mdtcutoff]

      zarea = 0.0
      zeta1 = 0.0
      zeta2 = 0.0

      DO jj = 1, jpj
         DO ji = 1, jpi
          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
          zarea = zarea + zdxdy
          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
        END DO      
      END DO

      IF( lk_mpp) CALL mpp_sum( zeta1 )
      IF( lk_mpp) CALL mpp_sum( zeta2 )
      IF( lk_mpp) CALL mpp_sum( zarea )
      
      zcorr_mdt = zeta1 / zarea
      zcorr_bcketa  = zeta2 / zarea

      !  Define correction term

      zcorr = zcorr_mdt - zcorr_bcketa

      !  Correct spatial mean of the MSSH

      IF ( nmsshc == 1 ) mdt(:,:) = mdt(:,:) - zcorr  

      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT

      IF ( nmsshc == 2 ) mdt(:,:) = mdt(:,:) - mdtcorr

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) ' obs_readmdt : mdtcutoff     = ', mdtcutoff
         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
         WRITE(numout,*) '               zcorr         = ', zcorr
         WRITE(numout,*) '               nmsshc        = ', nmsshc
         WRITE(numout,*) 
      ENDIF

      IF ( nmsshc == 0 ) WRITE(numout,*) &
         &                 '           MSSH correction is not applied'
      IF ( nmsshc == 1 ) WRITE(numout,*) &
         &                 '           MSSH correction is applied'
      IF ( nmsshc == 2 ) WRITE(numout,*) &
         &                 '           User defined MSSH correction' 


   END SUBROUTINE obs_offset_mdt
 
END MODULE obs_readmdt