MODULE obs_readsnowdepth
   !!======================================================================
   !!                       ***  MODULE obs_readsnowdepth  ***
   !! Observation diagnostics: Get the snow depth for freeboard conversion to thickness
   !!======================================================================
   !! History :      ! 2018-10 (E. Fiedler) Initial version
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   obs_rea_snowdepth    : Driver for reading MDT
   !!----------------------------------------------------------------------
   USE wrk_nemo         ! Memory Allocation
   USE par_kind         ! Precision variables
   USE par_oce          ! Domain parameters
   USE in_out_manager   ! I/O manager
   USE obs_surf_def     ! Surface observation definitions
   USE obs_inter_sup    ! Interpolation support routines
   USE obs_inter_h2d    ! 2D interpolation
   USE obs_utils        ! Various observation tools
   USE iom_nf90         ! IOM NetCDF
   USE netcdf           ! NetCDF library
   USE lib_mpp          ! MPP library
   USE dom_oce, ONLY : &                  ! Domain variables
      &                    tmask, tmask_i, e1t, e2t, gphit, glamt
   USE obs_const, ONLY :   obfillflt      ! Fillvalue
   USE oce      , ONLY :   sshn           ! Model variables



   IMPLICIT NONE
   PRIVATE
   
   PUBLIC   obs_rea_snowdepth     ! called by dia_obs

   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE obs_rea_snowdepth( fbddata, k2dint, thick_s )
      !!---------------------------------------------------------------------
      !!
      !!                   *** ROUTINE obs_rea_snowdepth ***
      !!
      !! ** Purpose : Get snowdepth from CICE
      !!
      !! ** Method  : 
      !!
      !! ** Action  : 
      !!----------------------------------------------------------------------
      USE iom
      !
      TYPE(obs_surf), INTENT(inout)  ::   fbddata   ! Sea ice freeboard data
      INTEGER       , INTENT(in)     ::   k2dint    ! ?
      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: thick_s ! Model snow depth
     
      CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_snowdepth'

      INTEGER ::   jobs                ! Obs loop variable
      INTEGER ::   jpi_thick_s, jpj_thick_s      ! Number of grid point in lat/lon for the snow depth
      INTEGER ::   iico, ijco          ! Grid point indices

      !
      REAL(wp), DIMENSION(1)     ::   zext, zobsmask
      REAL(wp), DIMENSION(2,2,1) ::   zweig
      !
      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, z_thick_s_l, zglam, zgphi
      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
      !
      REAL(wp), POINTER, DIMENSION(:,:) ::  z_thick_s, thick_s_mask
         
      REAL(wp) :: zlam, zphi, zfill, zinfill    ! local scalar
      !!----------------------------------------------------------------------

      CALL wrk_alloc(jpi,jpj,z_thick_s,thick_s_mask) 

      IF(lwp)WRITE(numout,*) 
      IF(lwp)WRITE(numout,*) ' obs_rea_snowdepth : Get model snow depth for freeboard conversion to sea ice thickness'
      IF(lwp)WRITE(numout,*) ' ------------- '
      CALL FLUSH(numout)

      ! Get ice thickness information
      z_thick_s = thick_s

      ! Setup mask based on tmask
      thick_s_mask(:,:) = tmask(:,:,1)

      ! Interpolate the snow depth already on the model grid at the observation point
  
      ALLOCATE( &
         & igrdi(2,2,fbddata%nsurf), &
         & igrdj(2,2,fbddata%nsurf), &
         & zglam(2,2,fbddata%nsurf), &
         & zgphi(2,2,fbddata%nsurf), &
         & zmask(2,2,fbddata%nsurf), &
         & z_thick_s_l(2,2,fbddata%nsurf)  &
         & )
         
      DO jobs = 1, fbddata%nsurf

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

      END DO

      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, glamt  , zglam )
      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, gphit  , zgphi )
      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, thick_s_mask, zmask )
      CALL obs_int_comm_2d( 2, 2, fbddata%nsurf, jpi, jpj, igrdi, igrdj, z_thick_s  , z_thick_s_l )

      DO jobs = 1, fbddata%nsurf
            
         zlam = fbddata%rlam(jobs)
         zphi = fbddata%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, z_thick_s_l(:,:,jobs),  zext )
 
         fbddata%rext(jobs,2) = zext(1)

! mark any masked data with a QC flag
         IF( zobsmask(1) == 0 )   fbddata%nqc(jobs) = IBSET(fbddata%nqc(jobs),15)

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

      CALL wrk_dealloc(jpi,jpj,z_thick_s,thick_s_mask) 
      IF(lwp)WRITE(numout,*) ' ------------- '
      !
   END SUBROUTINE obs_rea_snowdepth
   
   !!======================================================================
END MODULE obs_readsnowdepth