EnKF/sanity_check/p_sanity_check.F90

PROGRAM sanity_check_LIM3
! Checks if a LIM3 restart is physically consistent and outputs the updated
! version. Strongly based on limupdate.F90, and on sanity_check_LIM2
! (c) Original from Chris König Beatty
! (c) Re-written by Francois Massonnet, November 2012,
! (c) Refreshed for NEMO3.6, Francois Massonnet, April 2016
! francois.massonnet@uclouvain.be
!
! 
! TODO: create subroutines to automaticate the extraction of NetCDF!

! **************************************************************************************************
!                                           PLAN of the program
!
!
! 0. Header
! 1. Check existence of NetCDF files, grab command line args
!    1) Ice analysis file
!    2) Ice forecast file
!    3) Oce analysis file
!    4) Oce forecast file
!    5) Mask file
!    6) Mesh file          
! 2. Get dimensions, nb ice categories, nb ice/snow layers, allocate, get ice thickness bounds
! 3. Load ice+oce and mask+mesh variables from files
! 4. Ice analysis sanity check
!    1) Regularize ice concentration  and adapt other variables
!       accordingly
!    2) Manually update the snow and ice heat content according to volume
!       changes
!    3) Rebin categories with thickness out of bounds (CALL itd_reb() )
!       Note that itd_reb() calls shift_ice if a shift needs be done        
!    4) Several ice corrections
!       For example, set volume to zero where concentration is zero.         
!    5) Shrink ice (CALL shrink_ice() )
!       If sum(a_i) > 1 or < 0 or some a_i > 1 or < 0, then ice is shrunk 
!       or reset to zero.       
!    6) Rebin categories with thickness out of bounds (CALL itd_reb() )
!       We do it once again as thickness may have changed 
!       during processes 2) and 3)     
!    7) Final concentration check (CALL conc_check() )
!       The program stops if > 1 or < 0 occurs in sum or individual categories
!       This can ultimately cause trouble in the code as a term (1-sum(a_i))^0.6
!       occurs in a routine => negative power 0.6 is very baaaad
!       If ice concentration is just above 1 or just below zero (numerics) then
!       the program resets to zero or one.
! 5. Oce analysis sanity check
!    1) Compute the difference between forecast and analyzed sea salinity
!       If larger than XXX PSU, bound variations by that value
!       Also update the sea surface salinity variable accordingly 
!    2) Same for sea temperature; replace analysed temperature by forecast if
!    less than -2°C
!    3) Limit ocean velocities to [-2,2] m/s
! 6. Record the files. The original file is copied and the variables are dumped
!    in it
!    1) Ice analysis file
!    2) Oce analysis file      
!***************************************************************************************************
!

! 0. Header
! ---------
USE netcdf
USE my_variables
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: jl, ji, jj, jk
!

! 1. Grab Command Line Arguments
! ------------------------------
IF (iargc()==4) THEN
  CALL getarg(1, cfile_analysis_ice)
  CALL getarg(2, cfile_forecast_ice)
  CALL getarg(3, cfile_analysis_oce)
  CALL getarg(4, cfile_forecast_oce)
  
  ! Check if files exist
  ! 1) Ice analysis file
  INQUIRE(FILE=TRIM(cfile_analysis_ice), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cfile_analysis_ice)
    STOP
  END IF
  
  ! 2) Ice forecast file
  INQUIRE(FILE=TRIM(cfile_forecast_ice), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cfile_forecast_ice)
    STOP
  END IF

  ! 3) Oce analysis file
  INQUIRE(FILE=TRIM(cfile_analysis_oce), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cfile_analysis_oce)
    STOP
  END IF

  ! 4) Oce forecast file
  INQUIRE(FILE=TRIM(cfile_forecast_oce), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cfile_forecast_oce)
    STOP
  END IF

  ! 5) Mask file
  INQUIRE(FILE=TRIM(cmaskfile), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cmaskfile)
    STOP
  END IF

  ! 6) Mesh file
  INQUIRE(FILE=TRIM(cmeshfile), EXIST=l_ex)
  IF ( .NOT. l_ex ) THEN
    WRITE(*,*) "(sanity_check_lim3) File does not exist: " //&
    &TRIM(cmeshfile)
    STOP
  END IF

  ! Everything went OK:
  WRITE(*,*) "(sanity_check_lim3) Starting ..."
ELSE
  ! Write info
  WRITE(*,*)
  WRITE(*,*) " sanity_check_LIM3 needs arguments: "
  WRITE(*,*) " -analysis_file_ice "
  WRITE(*,*) " -forecast_file_ice "
  WRITE(*,*) " -analysis_file_oce "
  WRITE(*,*) " -forecast_file_oce "
  WRITE(*,*) " Checks NEMO-LIM3 ice and ocean analyses restarst (netcdf) file for sanity and fixes&
             &them if necessary."
  WRITE(*,*)
  WRITE(*,*) " Sanity means for now:"
  WRITE(*,*) " Strongly follow limupdate.F90"
  WRITE(*,*) " Files mask.nc and mesh_hgr.nc need to be in the current directory"
  WRITE(*,*)
  WRITE(*,*) " Hope to see you again soon."
  WRITE(*,*)
  WRITE(*,*) " Chris König Beatty "
  WRITE(*,*) " Francois Massonnet -- francois.massonnet@uclouvain.be"
  WRITE(*,*) " Last update: 2013"
  WRITE(*,*) " Last update: 2016 (to work with NEMO3.6)"

  STOP "(sanity_check): Stopped."
END IF ! iargc

! 2. Get dimensions, and allocate the variables
! ---------------------------------------------
CALL get_dimensions()
! Get number of ice categories
CALL get_nb_cat()
! Get number of ice layers
CALL get_nb_il()

! Get number of snow layers
CALL get_nb_sl()

! Allocate variables
CALL allocate_variables()


! Register ice thickness limits
! See Clement Rousset LIM3 paper or LIM3 doc:
! http://www.geosci-model-dev.net/8/2991/2015/gmd-8-2991-2015.pdf
! or the routine sbcice_lim.F90
!
! 
! 
! 
zalpha = 0.05             ! exponent of the transform function
rn_himean = 2.0           ! the expected mean thickness over the domain
zhmax  = 3.*rn_himean

DO jl = 1, jpl
  znum = jpl * ( zhmax+1 )**zalpha
  zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl
  hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
END DO

! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl)
hi_max(jpl) = 99._wp
WRITE(*,*), "Ice categories boundaries [m] are", hi_max


! 3. Load variables
! -----------------
CALL load_variables()
WRITE(*,*), "babysit after load_variables:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)

! 4. Ice analysis sanity check
! ----------------------------
! 1) Regularize ice concentration and other variables based on that
CALL regularize()


WRITE(*,*), "babysit after regularize:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "

! 2) Heat content update

! F. Massonnet - test CALL update_hc()
! CALL update_hc()

! 3) Rebin categories with thickness out of bounds
! The code here follows closely the contents of subroutine limi_itd_th_reb

CALL itd_reb()
WRITE(*,*), "babysit after itd_reb #1:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "

! 4) Several ice corrections

CALL several_ice_corrections()

WRITE(*,*), "babysit several ice corrections:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "

! 5) Shrink ice

CALL shrink_ice()

WRITE(*,*), "babysit after shrink_ice:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "

! 6) Rebin categories with thickness out of bounds

CALL itd_reb()

WRITE(*,*), "babysit after itd_reb #2:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "

! 7) Final check. 
! Stops if total conc > 1 or < 0, same for inidividual conc
! If exceeds 1 or is below zero for numerical reasons, reset. 

CALL conc_check()

WRITE(*,*), "babysit after conc check:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,1,2)
WRITE(*,*), " "


! 5. Ocean analysis sanity check
! ------------------------------
! 1) Salinity
CALL salinity_check()

! 2) Temperature
CALL temperature_check()

! 3) Velocity
CALL velocity_check()
! 6. Record NetCDF
! ----------------
! 1) Ice variables
CALL record_ice()
WRITE(*,*), "babysit after record ice:"
WRITE(*,*), "e_s_l01 i j", e_s(jiindx,jjindx,1,1)
WRITE(*,*), "e_s_l02 i j", e_s(jiindx,jjindx,2,1)
WRITE(*,*), "found negative val"
DO jk = 1, nlay_s
 DO jl = 1, jpl
  DO ji = 1, jpi
    DO jj = 1 , jpj
      IF ( e_s(ji,jj,jl,jk) < rzero ) THEN
       WRITE(*,*), "e_s < 0.0 in i=",ji," j=",jj
      ENDIF
    END DO ! jj
  END DO ! ji
 END DO ! jl
END DO ! jk

! 2) Oce variables
CALL record_oce()

END PROGRAM sanity_check_LIM3

! SUBROUTINES
! ¨¨¨¨¨¨¨¨¨¨¨

SUBROUTINE itd_reb()
  ! --------------------------------------------------------
  ! This routine is strongly inspired from lim_itd_th_reb
  ! When in situ thickness exceeds bounds, it transfers ice
  ! to neighbouring categories
  ! This routine calls shiftice() defined below
  ! --------------------------------------------------------
  USE my_variables
  IMPLICIT NONE
  !  Dummy variables
  INTEGER                 :: ji, jj, jl

  WRITE(*,*) 'Entering itd_reb()'
  ! 4.2.1 Compute in situ ice thickness in the categories (if there's ice)
DO jl = 1, jpl
  DO ji = 1, jpi
    DO jj = 1 , jpj
      IF ( a_i(ji,jj,jl) > epsi10 ) THEN
        ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
      ELSE
        ht_i(ji,jj,jl) = rzero
      ENDIF
    END DO ! jj
  END DO ! ji
END DO ! jl

! Print data at particular location before rebinning
!WRITE(*,*) 'Before rebinning: '
!WRITE(*,*) 'a_i : ',a_i(jiindx,jjindx,:)
!WRITE(*,*) 'ht_i : ',ht_i(jiindx,jjindx,:)
!WRITE(*,*) 'Total concentration: ', SUM(a_i(jiindx,jjindx,:))
!WRITE(*,*) 'Total volume       : ', SUM(v_i(jiindx,jjindx,:))

! 4.2.2 Make sure thickness of first category is at least hi_max_typ
! Not sure to understand what to do here

! 4.2.3. If a category is not in bounds, shift the entire area, volume and
! energy to the neighbouring category

! Initialize shift arrays
DO jl = 1, jpl
  idonor(:,:,jl) = 0
  ! idonor is the index of the category that is giving ice with respect to 
  ! the running category.
  ! Example: we are looping over categories. When jl = 3 (third category),
  ! we notice that this category has way too much ice. Then idonor(ji,jj,3)
  ! will be 3. Example 2 (jl=3): The fourth category has too much ice
  ! Then idonor(ji,jj,3) = 4
  zdaice(:,:,jl) = 0
  zdvice(:,:,jl) = 0
END DO

! 4.2.3.1 Move thin categories up
! Strategy:     -loop over all categories up to the last-but-one
!               -identify thicknesses in the current category  that are too big
!               -if applicable, transfer all volume, area and energy to the
!               jpl + 1 category
DO jl = 1, jpl - 1
  zshiftflag = 0
  DO ji = 1, jpi
    DO jj= 1, jpj
      IF ( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl)  ) THEN
        IF ( ldebug ) THEN
          WRITE(*,*), "Found excessive ice thickness in category",jl
          WRITE(*,*), ht_i(ji,jj,jl), "larger than" , hi_max(jl)  
          WRITE(*,*), "At grid point" , ji, jj 
        END IF
        zshiftflag        = 1     ! There's at least one point in the domain
                                 ! with too thick a sea ice in this category
        idonor(ji,jj,jl) = jl    ! The running category is then donor
        zdaice(ji,jj,jl) = a_i(ji,jj,jl) ! Amount of ice to be transferred
        zdvice(ji,jj,jl) = v_i(ji,jj,jl) ! (everything)
      END IF
    END DO! jj
  END DO ! ji

  IF (zshiftflag == 1) THEN ! this is the case if there was at least one
                            ! excessive thickness in the running category
                            ! somewhere on the domain
    CALL shiftice()
    idonor(:,:,jl) = 0        ! Reset before we move to next category
    zdaice(:,:,jl) = rzero
    zdvice(:,:,jl) = rzero
  END IF
END DO ! jl, categories
! 4.2.3.2 Move thick categories down
! Strategy:     -loop over all categories starting from last-but-one 
!                to the first one
!               -identify if the thickness of the category above is smaller
!                than the upper limit for the running category
!               -if so, transfer everything to the running category

DO jl = jpl-1, 1, -1    ! first, last, step
  zshiftflag = 0
  DO ji = 1, jpi
    DO jj= 1, jpj
      IF ( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl)  ) THEN
        IF ( ldebug ) THEN
          WRITE(*,*), "Found too small ice thickness in category",jl+1
          WRITE(*,*), ht_i(ji,jj,jl+1), "smaller than" , hi_max(jl)
        END IF
        zshiftflag        = 1     ! There's at least one point in the domain
                                  ! with too thin a sea ice in this category
        idonor(ji,jj,jl) = jl+1   ! The jl+1  category is then donor
        zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) ! Amount of ice to be transferred
        zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) ! (everything)
      END IF
    END DO! jj
  END DO ! ji
  IF (zshiftflag == 1) THEN ! this is the case if there was at least one
                            ! too small thickness in the jl+1 category
    CALL shiftice()
    idonor(:,:,jl) = rzero        ! Reset before we move to next category
    zdaice(:,:,jl) = rzero
    zdvice(:,:,jl) = rzero
  END IF
END DO ! jl, categories
! 4.2.3.3 Conservation check
! Print data at particular location after rebinning
!WRITE(*,*) 'After rebinning: '
!WRITE(*,*) 'a_i : ',a_i(jiindx,jjindx,:)
!WRITE(*,*) 'ht_i : ',ht_i(jiindx,jjindx,:)
!WRITE(*,*) 'Total concentration: ', SUM(a_i(jiindx,jjindx,:))
!WRITE(*,*) 'Total volume       : ', SUM(v_i(jiindx,jjindx,:))


WRITE(*,*) 'Leaving itd_reb()'
END SUBROUTINE itd_reb

SUBROUTINE shiftice()

  !-------------------------------------------------------------
  ! This routine is (strongly) inspired by lim_itd_shiftice
  ! It re-arranges ice thickness in the categories
  ! 
  ! This routine can potentially be called 2*4 times,
  ! for the first jpl-1 categories to upgrade their too thick ice
  ! and for the jpl-1 last categories to downgrade their too thin ice
  !
  !
  ! Francois Massonnet UCL 2012 francois.massonnet@uclouvain.be
  !-------------------------------------------------------------

  ! Variable declaration and importation

  USE my_variables
  
  IMPLICIT NONE
  LOGICAL ::  zdaice_negative, zdvice_negative, zdaice_greater_aicen,&
                &zdvice_greater_vicen
   
  ! Number of cells with ice to transfer
  REAL(wp),  DIMENSION(jpi,jpj)      :: zworka
  REAL(wp),  DIMENSION(jpi,jpj,jpl)  :: zaTsfn
  REAL(wp)                                  :: zdvsnow, zdesnow, zdo_aice
                                          ! Volume of snow transferred,
                                          ! Snow heat, ice age
                                          ! local variables
  REAL(wp)                                  :: zdsm_vice, zdaTsf, zdeice
                                          ! ice age, surface temperature,
                                          ! ice heat content
                                          ! local variables
  !  Dummy variables
  INTEGER                 :: ji, jj, jl, jk, jl1, jl2
  ! End definitions
  ! ----------------------------------------------------------------------

  ! Welcome the user
  WRITE(*,*) 'Entering SUBROUTINE shiftice'
  ! Define surface temperature times concentration
  ! This has more dimensions like energy. It will be used 
  ! when surface temperature will be "transferred" after rebinning
  DO jl = 1 , jpl
    zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
  END DO  
  

  ! Subroutine

  DO jl = 1 , jpl - 1
    zdaice_negative = .FALSE.
    zdvice_negative = .FALSE.
    zdaice_greater_aicen = .FALSE.
    zdvice_greater_vicen = .FALSE.

    DO ji = 1 , jpi 
      DO jj = 1 , jpj
   
        !---------------------------------------------------------------------
        ! 1) Check for daice or dvice out of bounds and reset if necessary
        !---------------------------------------------------------------------

        IF ( idonor(ji,jj,jl) .GT. 0 ) THEN ! If the running category is giving
          jl1=idonor(ji,jj,jl)              ! record the donor category
          !WRITE(*,*), "At grid point ", ji, jj, "category", jl1, "is giving ice"
         
          ! Tackle the cases with problems. Normally, zdaice and zdvice should
          ! be positive, but ...
          ! Check for negative transfers of ice given in input
 
          ! 1. Ice area
          IF (zdaice(ji,jj,jl) .LT. rzero ) THEN
            IF (zdaice(ji,jj,jl) .GT. -epsi10 ) THEN
              WRITE(*,*)," zdaice is negative but artificial"
              IF ( ( jl1 .EQ. jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )&
                   &.OR.&
                   &( jl1 .EQ. jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )) THEN
                ! You are here if one of these 2 conditions are verified
                ! 1) The running category is the donor and has too thick a
                ! sea ice
                ! 2) The running category is one category below the donor,
                ! which has too thin a sea ice
                !
                ! If you still with me, it means that a transfer needs to be
                ! done but the amount of ice is negative due to roundoff
                ! error or something. Let us reset zdaice and zdvice to the 
                ! value of a_i and v_i in the category for the transfer
                zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
                zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
               
              ELSE ! None of the two conditions above is verified
                   ! That is, 
                   ! 1) The running category is not the donor or has a correct
                   ! ice
                   ! AND
                   ! 2) The category above the running one is not the donor, or
                   ! the ice is correct in this above category
                   ! 
                   ! Since everything is fine, nothing should be done
                zdaice(ji,jj,jl) = rzero
                zdvice(ji,jj,jl) = rzero 
              END IF
            ELSE ! zdaice < - eps10
              WRITE(*,*) "zdaice is really negative"
              zdaice_negative = .TRUE.
            END IF
          END IF ! zdaice < 0        

          ! 2. Repeat for ice volume
          IF (zdvice(ji,jj,jl) .LT. rzero ) THEN
            IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN
              WRITE(*,*)," zdvice is negative but artificial"
              IF ( ( jl1 .EQ. jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )&
                   &.OR.&
                   &( jl1 .EQ. jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )) THEN
                ! You are here if one of these 2 conditions are verified
                ! 1) The running category is the donor and has too thick a
                ! sea ice
                ! 2) The running category is one category below the donor,
                ! which has too thin a sea ice
                !
                ! If you still with me, it means that a transfer needs to be
                ! done but the amount of ice is negative due to roundoff
                ! error or something. Let us reset zdaice and zdvice to the 
                ! value of a_i and v_i in the category for the transfer
                zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
                zdvice(ji,jj,jl) = v_i(ji,jj,jl1)

              ELSE ! None of the two conditions above is verified
                   ! That is, 
                   ! 1) The running category is not the donor or has a correct
                   ! ice
                   ! AND
                   ! 2) The category above the running one is not the donor, or
                   ! the ice is correct in this above category
                   ! 
                   ! Since everything is fine, nothing should be done
                zdaice(ji,jj,jl) = rzero
                zdvice(ji,jj,jl) = rzero
              END IF
            ELSE ! zdvice < - eps10
              WRITE(*,*) "zdvice is really negative"
              zdvice_negative = .TRUE.
            END IF
          END IF ! zdvice < 0        

          ! 3.A. If the area to be transferred equals the area in the running
          ! category, then just update it to the exact value
          ! (i.e. round it )
          IF   ( zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10  .AND.&
             &   zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1) + epsi10  ) THEN
            ! if concentration to be transferred is "equal"
            ! to the concentration of the donor. This is obviously the case
            ! if the running category is the donor 
            zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
            zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
          ELSE
            ! Otherwise, set the switch to true
            zdaice_greater_aicen = .TRUE.  
          END IF !  zdaice "=" a_i

          ! 3.B. (Repeat for volume) 
          ! If the volume to be transferred equals the volume in the running
          ! category, then keep it as is          
          IF   ( zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1) - epsi10  .AND.&
             &   zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1) + epsi10  ) THEN
            ! if volume to be transferred is "equal"
            ! to the volume of the donor. This is obviously the case
            ! if the running category is the donor 
            zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
            zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
          ELSE
            zdvice_greater_vicen = .TRUE.
          END IF !  zdaice "=" a_i

        ELSE   ! if idonor
          ! Nothing happens since this category is not giving ice
        END IF ! if idonor

      END DO ! jj
    END DO ! ji
  END DO ! jl, category


  !-------------------------------------------------
  ! 2) Transfer volume and energy between categories
  !-------------------------------------------------
  DO jl = 1, jpl - 1
    DO ji = 1, jpi
      DO jj = 1, jpj
        IF (zdaice(ji,jj,jl) .GT. rzero ) THEN
          jl1 = idonor(ji,jj,jl)
          ! Define proportionality factor.
          ! zworka will be the ratio between transferred volume of ice and
          ! actual volume of ice in the category. Auxiliary quantities (snow volume, snow
          ! heat content, ice age, salinity, etc.) will be transferred following
          ! this ratio. Indeed, if out of 4 m of ice, 1 m is transferred, then
          ! 1/4 of the snow volume will be transferred, too.
          zindb = MAX( rzero , SIGN( rone , v_i(ji,jj,jl1) - epsi10 )  ) 
          ! Thus zindb is equal to 1 if the donor has positive and not artificially
          ! close to zero ice volume to give, to zero otherwise
          zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10  ) *&
                &zindb
          ! zworka is zero where the donor has no ice, otherwise
          ! equal to zdvice/vice of the current category


          ! We have a donor, but who is the benefiter? who will receive ice?
          IF ( jl1 == jl ) THEN ! If the donor is the current category, then
            jl2=jl1+1           ! the receiver is the one above
          ELSE                  ! Otherwise (the donor is the above categ)
            jl2=jl              ! then it's the current category!
          END IF

          ! -----------------------
          ! Go for the transfers!!!
          ! -----------------------
          
          ! A. Ice areas
          ! ------------
          a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - zdaice(ji,jj,jl)
          a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl)
          ! The donor loses , the receiver gains

          ! B. Ice volumes
          ! --------------
          v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - zdvice(ji,jj,jl)
          v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl)

          ! C. Snow volumes
          ! --------------- 
          zdvsnow = v_s(ji,jj,jl1) * zworka(ji,jj) 
          v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - zdvsnow
          v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + zdvsnow

          ! D. Snow heat content
          ! --------------------
          IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN
           WRITE(*,*) 'Before itd_reb e_s_l01:' , e_s(ji,jj,1,jl)
         END IF

          DO jk = 1 , nlay_s
            zdesnow = e_s(ji,jj,jk,jl1) * zworka(ji,jj) 
            e_s(ji,jj,jk,jl1) = e_s(ji,jj,jk,jl1) - zdesnow
            e_s(ji,jj,jk,jl2) = e_s(ji,jj,jk,jl2) + zdesnow
          IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN
           WRITE(*,*) 'Res shrink jl1:' , e_s(ji,jj,jk,jl1)
           WRITE(*,*) 'Res shrink jl2:' , e_s(ji,jj,jk,jl2)
         END IF
         END DO

         IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN
           WRITE(*,*) 'after itd_reb e_s_l01:' , e_s(ji,jj,1,jl)
         END IF

          ! E. Ice age
          ! ----------
          ! Ice age is defined as areal. If 1/2 of the area
          ! is transferred then 1/2 of the age too
          ! Wat een cromme definitie!
          zdo_aice = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl)
          oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - zdo_aice
          oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + zdo_aice

          ! F. Ice salinity
          ! ---------------
          zdsm_vice = smv_i(ji,jj,jl1) * zworka(ji,jj)
          smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - zdsm_vice
          smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + zdsm_vice

          ! G. Surface temperature
          ! ----------------------
          zdaTsf = t_su(ji,jj,jl1) * zdaice(ji,jj,jl)
          zaTsfn(ji,jj,jl1) = zaTsfn(ji,jj,jl1) - zdaTsf
          zaTsfn(ji,jj,jl2) = zaTsfn(ji,jj,jl2) + zdaTsf

          ! H. Ice heat content
          ! -------------------
          DO jk = 1 , nlay_i
            zdeice = e_i(ji,jj,jk,jl1) * zworka(ji,jj)
            e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - zdeice
            e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + zdeice
          END DO
        ELSE
          ! WRITE(*,*),"Nothing to transfer"
        END IF
      END DO ! jpj
    END DO ! jpi
  END DO ! jl, category

  ! ---------------------------------------
  ! 3) Update ice thickness and temperature
  ! ---------------------------------------
  DO jl = 1 , jpl
    DO ji = 1 , jpi
      DO jj = 1 , jpj
        IF ( a_i(ji,jj,jl) > epsi10 ) THEN
          ht_i(ji,jj,jl) = v_i(ji,jj,jl)    / a_i(ji,jj,jl)
          t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)
        ELSE
          ht_i(ji,jj,jl) = rzero
          t_su(ji,jj,jl) = rtt  ! If no ice then set "ice" temperature to
                                ! freezing point
        END IF
      END DO ! jj
    END DO !ji
  END DO ! jl  
  WRITE(*,*) 'Leaving SUBROUTINE shiftice'
END SUBROUTINE shiftice

SUBROUTINE conc_check()
  USE my_variables
  IMPLICIT NONE
  
  !  Dummy variables
  INTEGER                 :: ji, jj, jl


  DO ji=1,jpi
  DO jj = 1,jpj

    DO jl=1,jpl
      ! 1. Individual concentrations greater than zero
      IF ( a_i(ji,jj,jl) .LT. rzero  ) THEN
        IF (a_i(ji,jj,jl) .LT. -epsi10) THEN ! We have a "true" negative conc
          WRITE(*,*) 'ERROR: final check: a_i negative at ',ji,jj
          WRITE(*,*) 'for category ',jl
          WRITE(*,*) 'Value: ', a_i(ji,jj,jl)
          WRITE(*,*) 'ABORT'
          STOP
        ELSE                             ! We have a fake negative conc
          a_i(ji,jj,jl) = rzero
        ENDIF
      ENDIF

      ! 2. Individual concentrations less than one
      IF ( a_i(ji,jj,jl) .GT. zamax  ) THEN
        IF (a_i(ji,jj,jl) - zamax .GT. epsi10) THEN  ! We have a "true" positive conc
          WRITE(*,*) 'ERROR: final check: a_i more than zamax at ',ji,jj
          WRITE(*,*) 'for category ',jl
          WRITE(*,*) 'Value: ', a_i(ji,jj,jl)
          WRITE(*,*) 'ABORT'
          STOP
        ELSE                             ! We have a fake more than one conc
          a_i(ji,jj,jl) = zamax
        ENDIF
      ENDIF
    END DO ! jl  

    ! 3. Total concentration greater than zero
    IF (  SUM(a_i(ji,jj,:)) .LT. rzero  ) THEN
      IF (SUM(a_i(ji,jj,:)) .LT. -epsi10) THEN  ! We have a "true" negative total conc
        WRITE(*,*) 'ERROR: final check: at_i negative at ',ji,jj
        WRITE(*,*) 'Value: ', SUM(a_i(ji,jj,:))
        WRITE(*,*) 'ABORT'
        STOP
      ELSE                             ! We have a fake negative conc
        ! Let's reset all categories to zero
        DO jl = 1, jpl
          a_i(ji,jj,jl)=rzero
        END DO
      ENDIF
    ENDIF

    ! 4. Total concentration less than one
    
    IF (  SUM(a_i(ji,jj,:)) .GT. zamax  ) THEN
      IF (SUM(a_i(ji,jj,:)) - zamax .GT. epsi10) THEN  ! We have a "true"positive conc
        WRITE(*,*) 'ERROR: final check: at_i more than one at ',ji,jj
        WRITE(*,*) 'Value: ', SUM(a_i(ji,jj,:))
        WRITE(*,*) 'Individual: ', a_i(ji,jj,:)
        WRITE(*,*) 'ABORT'
        STOP
      ELSE                             ! We have a fake more than one conc
        ! Define the excess, which is by definition negligible
        zda_i = SUM(a_i(ji,jj,:)) - zamax ! (positive)
        ! Give this excess to the first category than can accept it, i.e. that
        ! has less than zamax - zda_i
        l_adjust = .TRUE.
        DO jl = 1, jpl
          IF ( ( a_i(ji,jj,jl) .GT. zda_i ) .AND. l_adjust ) THEN
            a_i(ji,jj,jl) = a_i(ji,jj,jl) - zda_i
            l_adjust = .FALSE.
          END IF
        ENDDO ! jl
        IF ( l_adjust ) THEN
        WRITE(*,*) 'It was not possible to give the excess of ice.'
            WRITE(*,*) 'At grid point ',ji,jj
            WRITE(*,*) 'Category', jl
            WRITE(*,*) 'The excess is ', zda_i
            WRITE(*,*) 'The conc. in categories are ', a_i(ji,jj,:)
            WRITE(*,*) ' ABORT'
            STOP
        END IF
      ENDIF
    ENDIF

  END DO !jj
END DO !ji  

END SUBROUTINE conc_check

SUBROUTINE get_dimensions()
USE my_variables
USE netcdf
IMPLICIT NONE

!  Dummy variables

WRITE(*,*), 'Entering get_dimensions()'

! 2.A Get dimensions
! ------------------
! Open the oce restart file
ierr = nf90_open(TRIM(cfile_analysis_oce),nf90_Write,incid_oce_an_in)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Opening. Abort"
  WRITE(*,*), "File :" // TRIM(cfile_analysis_oce)
  STOP
END IF

! Inquire variable id (here the variable does not matter!)
ierr = nf90_inq_varid(incid_oce_an_in, "sn", ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Variable ID Inquiring. Abort"
  WRITE(*,*), "File :"//TRIM(cfile_analysis_oce)
  WRITE(*,*), "Variable: sn"
  STOP
END IF

! Inquire variable
ierr = nf90_inquire_variable(incid_oce_an_in, ivarid, dimids=idimid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Variable Inquiring. Abort"
  WRITE(*,*), "File :"//TRIM(cfile_analysis_oce)
  WRITE(*,*), "Variable: sn"
  STOP
END IF


! Get dimensions
ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(1), len=jpi)
ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(2), len=jpj)
ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(3), len=jpk)

IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF dimensions Inquiring. Abort"
  WRITE(*,*), "File :"//TRIM(cfile_analysis_oce)
  STOP
END IF

WRITE(*,*), "The model horizontal dimensions are" , jpi, "by",jpj
WRITE(*,*), "The model vertical   dimension   is" , jpk

WRITE(*,*), 'Leaving get_dimensions()'

END SUBROUTINE get_dimensions

SUBROUTINE get_nb_cat()
  USE my_variables
  USE netcdf
  IMPLICIT NONE
  !  Dummy variables
  INTEGER                 :: jl

  WRITE(*,*) 'Entering get_nb_cat()'

  ! Open file
  ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_SI3) Error NetCDF Ice file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  ! Get number of ice categories
  ! (2024-NEMOv4.2.1-SI3) Categories are now a dimension of variable, changes in consequence (Alison Delhasse)
  
  !cvarroot="a_i_htc"
  !jl=1
  !WRITE(cinteger,"(i1)") jl
  !ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger), ivarid)
  !jl=jl+1
  !
  !DO WHILE (ierr == nf90_noerr) 
  !  WRITE(cinteger,"(i1)") jl
  !  !WRITE(*,*),"Checking existence of " // TRIM(cconcroot)//TRIM(cinteger)
  !  ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//TRIM(cinteger), ivarid)
  !  jl=jl+1
  !ENDDO
  !
  !jpl=jl-2
  !WRITE(*,*), "There are ", jpl, "ice categories"

  ! Inquire variable id (here the variable does not matter!)
  ierr = nf90_inq_varid(incid_ice_an_in, "a_i", ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
     WRITE(*,*), "(sanity_check_SI3) Error NetCDF Variable ID Inquiring. Abort"
     WRITE(*,*), "File :"//TRIM(cfile_analysis_ice)
     WRITE(*,*), "Variable: a_i"
     STOP
  END IF

  ! Inquire variable
  ierr = nf90_inquire_variable(incid_ice_an_in, ivarid, dimids=idimid)
  IF (ierr .NE. nf90_noerr ) THEN
     WRITE(*,*), "(sanity_check_SI3) Error NetCDF Variable Inquiring. Abort"
     WRITE(*,*), "File :"//TRIM(cfile_analysis_ice)
     WRITE(*,*), "Variable: a_i"
     STOP
  END IF


  ! Get dimensions
  ierr = nf90_inquire_dimension(incid_ice_an_in, idimid(3), len=jpl)
  !ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(2), len=jpj)
  !ierr = nf90_inquire_dimension(incid_oce_an_in, idimid(3), len=jpk)

  IF (ierr .NE. nf90_noerr ) THEN
     WRITE(*,*), "(sanity_check_SI3) Error NetCDF dimensions Inquiring. Abort"
     WRITE(*,*), "File :"//TRIM(cfile_analysis_ice)
     STOP
  END IF


  WRITE(*,*), "There are ", jpl, "ice categories"

  ! Close
  ierr = nf90_close(incid_ice_an_in);
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  WRITE(*,*) 'Leaving get_nb_cat()'
  
END SUBROUTINE get_nb_cat

SUBROUTINE get_nb_il()
  USE my_variables
  USE netcdf
  IMPLICIT NONE

  !  Dummy variables
  INTEGER                 ::  jk
  
  WRITE(*,*) 'Entering get_nb_il()'
  ! Open file
  ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  ! Changes in variables names in SI3 (2024) - changes in the subroutines Alison Delhasse
  ! cvarroot="tempt_il"
  cvarroot="e_i"
  jk=1
  WRITE(cinteger,"(i1)") jk
  ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//'_l0'//TRIM(cinteger), ivarid)
  jk=jk+1
  DO WHILE (ierr == nf90_noerr)
    WRITE(cinteger,"(i1)") jk
    ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//'_l0'//TRIM(cinteger), ivarid)
    jk=jk+1
  ENDDO
  nlay_i = jk-2
  WRITE(*,*), "There are ", nlay_i, "ice layers"
 
  ! Close
  ierr = nf90_close(incid_ice_an_in);
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  WRITE(*,*) 'Leaving get_nb_il()'
END SUBROUTINE get_nb_il

SUBROUTINE get_nb_sl()
  USE my_variables
  USE netcdf
  IMPLICIT NONE

  !  Dummy variables
  INTEGER                 :: jk

  WRITE(*,*) 'Entering get_nb_sl()'

  ! Open file
  ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  ! Changes in variables names in SI3 (2024) - changes in the subroutines Alison Delhasse
  cvarroot="e_s"
  jk=1
  WRITE(cinteger,"(i1)") jk
  ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//'_l0'//TRIM(cinteger), ivarid)
  jk=jk+1
  DO WHILE (ierr == nf90_noerr)
    WRITE(cinteger,"(i1)") jk
    ierr = nf90_inq_varid(incid_ice_an_in, TRIM(cvarroot)//'_l0'//TRIM(cinteger),ivarid)
    jk=jk+1
  ENDDO
  nlay_s = jk-2
  WRITE(*,*), "There are ", nlay_s, "snow layers"

  ! Close
  ierr = nf90_close(incid_ice_an_in);
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    STOP
  END IF

  WRITE(*,*) 'Leaving get_nb_sl()'
END SUBROUTINE get_nb_sl

SUBROUTINE allocate_variables()
  USE my_variables
  IMPLICIT NONE

WRITE(*,*), 'Entering allocate_variables()'

ALLOCATE( ilandmask(jpi, jpj)              )
ALLOCATE( a_i(      jpi, jpj,       jpl)  ,&
         &v_i(      jpi, jpj,       jpl)  ,&
         &v_i_fc(   jpi, jpj,       jpl)  ,&
         &v_s(      jpi, jpj,       jpl)  ,&
         &v_s_fc(   jpi, jpj,       jpl)  ,&
         &oa_i (    jpi, jpj,       jpl)  ,&
         &smv_i(    jpi, jpj,       jpl)  ,&
         &t_su(     jpi, jpj,       jpl)   )
ALLOCATE( ht_i(     jpi, jpj,       jpl)   )
ALLOCATE( e_i(      jpi, jpj,nlay_i,jpl)   )
ALLOCATE( e_s(      jpi, jpj,nlay_s,jpl)   )
ALLOCATE( sss_m_an( jpi, jpj           )  ,&
          sss_m_fc( jpi, jpj           )  ,&
          sst_m_an( jpi, jpj           )  ,&
          sst_m_fc( jpi, jpj           )   )
ALLOCATE( sn_an   ( jpi, jpj,jpk       )  ,&
          sn_fc   ( jpi, jpj,jpk       )  ,&
          tn_an   ( jpi, jpj,jpk       )  ,&
          tn_fc   ( jpi, jpj,jpk       )  ,&
          un_an   ( jpi, jpj,jpk       )  ,&
          un_fc   ( jpi, jpj,jpk       )  ,&
          ub_an   ( jpi, jpj,jpk       )  ,&
          ub_fc   ( jpi, jpj,jpk       )  ,&
          vn_an   ( jpi, jpj,jpk       )  ,&
          vn_fc   ( jpi, jpj,jpk       )  ,&
          vb_an   ( jpi, jpj,jpk       )  ,&
          vb_fc   ( jpi, jpj,jpk       )   )
ALLOCATE( ssu_m_an( jpi, jpj           )  ,&
          ssu_m_fc( jpi, jpj           )  ,&
          ssv_m_an( jpi, jpj           )  ,&
          ssv_m_fc( jpi, jpj           )   )
ALLOCATE( hi_max(                   jpl)   )
ALLOCATE( idonor(   jpi, jpj,       jpl)   )
ALLOCATE( internal_melt(jpi, jpj,   jpl)   )
ALLOCATE( zdaice(   jpi, jpj,       jpl)  ,&
         &zdvice(   jpi, jpj,       jpl)   )
ALLOCATE( ze1t(     jpi, jpj           )  ,&
         &ze2t(     jpi, jpj           )  ,&
         &zcellarea(jpi, jpj           )   )
ALLOCATE(zheat_res( jpi, jpj           )  ,&
         &zdmsnif(  jpi, jpj           )   )
ALLOCATE(at_i(      jpi, jpj           )  ,&
         snwice_mass(jpi,jpj           )  ,&
         snwice_mass_b(jpi,jpj         )   )
  
WRITE(*,*) 'Leaving allocate_variables()'
END SUBROUTINE allocate_variables

SUBROUTINE load_variables()

USE my_variables
USE netcdf
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: jl, jk

WRITE(*,*) 'Entering load_variables()'
! 3.A Mask
! --------
! Open
ierr = nf90_open(TRIM(cmaskfile),nf90_NoWrite,incid_mask)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Mask file. Abort"
  WRITE(*,*), TRIM(cmaskfile)
  STOP
END IF

! Request variable ID
ierr = nf90_inq_varid(incid_mask, cmaskvar, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort"
  WRITE(*,*), TRIM(cmaskfile)
  STOP
END IF

! Get the NetCDF variable and put it in the Fortran variable
ierr = nf90_get_var(incid_mask, ivarid, ilandmask)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort"
  WRITE(*,*), TRIM(cmaskfile)
  STOP
END IF

! Close mask file
ierr = nf90_close(incid_mask)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort"
  WRITE(*,*), TRIM(cmaskfile)
  STOP
END IF

! 3.B Mesh file
! -------------
! Open
ierr = nf90_open(TRIM(cmeshfile),nf90_NoWrite,incid_mesh)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Mask file. Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

WRITE(*,*), "Mesh loaded"

! Request variable ID
ierr = nf90_inq_varid(incid_mesh, ce1tvar, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

! Get the NetCDF variable and put it in the Fortran variable
ierr = nf90_get_var(incid_mesh, ivarid, ze1t)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

! Repeat for e2t
! Request variable ID
ierr = nf90_inq_varid(incid_mesh, ce2tvar, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Inquire varID . Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

! Get the NetCDF variable and put it in the Fortran variable
ierr = nf90_get_var(incid_mesh, ivarid, ze2t)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF getVar varID . Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

! Close mesh file
ierr = nf90_close(incid_mesh)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort"
  WRITE(*,*), TRIM(cmeshfile)
  STOP
END IF

! Computing zcellarea
zcellarea(:,:) = ze1t(:,:) * ze2t(:,:)
! 3.C Ocean variables
! -------------------
!
! 3.C.1 Open analysis
ierr = nf90_open(trim(cfile_analysis_oce),nf90_NoWrite,incid_oce_an_in)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF


! The following lines commented out because SSS does not appear anymore
! as restart variable (2020)

!! 3.C.1.A Sea surface salinity
!! Request variable ID
!ierr = nf90_inq_varid(incid_oce_an_in, csss_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!
!! Get variable
!ierr = nf90_get_var(incid_oce_an_in, ivarid, sss_m_an)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!
! 3.C.1.B Salinity
! Request variable ID
ierr = nf90_inq_varid(incid_oce_an_in, csn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, sn_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! The following lines commented out because SST does not longer
! appear as restart file (2020)
!! 3.C.1.C Sea surface temperature
!! Request variable ID
!ierr = nf90_inq_varid(incid_oce_an_in, csst_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!
!! Get variable
!ierr = nf90_get_var(incid_oce_an_in, ivarid, sst_m_an)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF

! 3.C.1.D. Temperature
ierr = nf90_inq_varid(incid_oce_an_in, ctn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, tn_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! 3.C.1.E U-velocity ("un")
ierr = nf90_inq_varid(incid_oce_an_in, cun, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, un_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! 3.C.1.F U-velocity ("ub")
ierr = nf90_inq_varid(incid_oce_an_in, cub, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, ub_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! 3.C.1.G V-velocity ("vn")
ierr = nf90_inq_varid(incid_oce_an_in, cvn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, vn_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! 3.C.1.H V-velocity ("vb")
ierr = nf90_inq_varid(incid_oce_an_in, cvb, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF
  
! Get variable
ierr = nf90_get_var(incid_oce_an_in, ivarid, vb_an)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! The following lines commented out because surface velocities no longer
! appear as restart variables (2020)

!! 3.C.1.I SSU-velocity ("ssu_m")
!ierr = nf90_inq_varid(incid_oce_an_in, cssu_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!  
!! Get variable
!ierr = nf90_get_var(incid_oce_an_in, ivarid, ssu_m_an)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!
!! 3.C.1.J SSV-velocity ("ssv_m")
!ierr = nf90_inq_varid(incid_oce_an_in, cssv_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF
!  
!! Get variable
!ierr = nf90_get_var(incid_oce_an_in, ivarid, ssv_m_an)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_analysis_oce)
!  STOP
!END IF

! Close analysis
ierr = nf90_close(incid_oce_an_in);
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean close file. Abort"
  WRITE(*,*), TRIM(cfile_analysis_oce)
  STOP
END IF

! 3.C.2. Open forecast
ierr = nf90_open(trim(cfile_forecast_oce),nf90_NoWrite,incid_oce_fc_in)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF


! Following lines commented out as the variable no longer appears 
! in the restart files (2020

!! 3.C.2.A Sea surface salinity
!! Request variable ID
!ierr = nf90_inq_varid(incid_oce_fc_in, csss_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF
!
!! Get variable
!ierr = nf90_get_var(incid_oce_fc_in, ivarid, sss_m_fc)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF

! 3.C.2.B. Salinity
! Request variable ID
ierr = nf90_inq_varid(incid_oce_fc_in, csn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, sn_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Following lines commented out as the variable no longer appears
! in the restart files (2020)

!! 3.C.2.C Sea surface temperature
!! Request variable ID
!ierr = nf90_inq_varid(incid_oce_fc_in, csst_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF
!
!! Get variable
!ierr = nf90_get_var(incid_oce_fc_in, ivarid, sst_m_fc)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF

! 3.C.2.D. Temperature
! Request variable ID
ierr = nf90_inq_varid(incid_oce_fc_in, ctn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, tn_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! 3.C.2.E U-velocity ("un")
ierr = nf90_inq_varid(incid_oce_fc_in, cun, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, un_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! 3.C.2.F U-velocity ("ub")
ierr = nf90_inq_varid(incid_oce_fc_in, cub, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, ub_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! 3.C.2.G V-velocity ("vn")
ierr = nf90_inq_varid(incid_oce_fc_in, cvn, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, vn_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! 3.C.2.H V-velocity ("vb")
ierr = nf90_inq_varid(incid_oce_fc_in, cvb, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF
  
! Get variable
ierr = nf90_get_var(incid_oce_fc_in, ivarid, vb_fc)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF

! Following lines commented out as the variable no longer appears
! in the restart files (2020)

!! 3.C.2.I SSU-velocity ("ssu_m")
!ierr = nf90_inq_varid(incid_oce_fc_in, cssu_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF
!  
!! Get variable
!ierr = nf90_get_var(incid_oce_fc_in, ivarid, ssu_m_fc)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF

!! 3.C.2.J SSV-velocity ("ssv_m")
!ierr = nf90_inq_varid(incid_oce_fc_in, cssv_m, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF
!  
! Get variable
!ierr = nf90_get_var(incid_oce_fc_in, ivarid, ssv_m_fc)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file get. Abort"
!  WRITE(*,*), TRIM(cfile_forecast_oce)
!  STOP
!END IF



! Close forecast
ierr = nf90_close(incid_oce_fc_in);
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean close file. Abort"
  WRITE(*,*), TRIM(cfile_forecast_oce)
  STOP
END IF


WRITE(*,*), "Ocean variables loaded"

! 3.D Ice variables
! -----------------

!#################!
!# Open forecast #
!#################!

ierr = nf90_open(trim(cfile_forecast_ice),nf90_NoWrite,incid_ice_fc_in)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort"
  WRITE(*,*), TRIM(cfile_forecast_ice)
  STOP
END IF


! (2024) ice category is now a dimension of the variables. No need to loop with the cat index to load variables. 
!DO jl = 1 , jpl

  ! Read ice volume forecast
  cvarroot='v_i'
  cvarname = TRIM(cvarroot)

  ! Inquire variable ID
  ierr = nf90_inq_varid(incid_ice_fc_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_forecast_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Read variable
  ierr = nf90_get_var(incid_ice_fc_in, ivarid, v_i_fc)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_forecast_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Read snow volume forecast
  cvarroot='v_s'
  cvarname = TRIM(cvarroot)

  ! Inquire variable ID
  ierr = nf90_inq_varid(incid_ice_fc_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_forecast_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Read variable
  ierr = nf90_get_var(incid_ice_fc_in, ivarid, v_s_fc)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_forecast_ice)
    WRITE(*,*), cvarname
    STOP
  END IF
!END DO ! jl

! Close forecast
ierr = nf90_close(incid_ice_fc_in);
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort"
  WRITE(*,*), TRIM(cfile_forecast_ice)
  STOP
END IF

WRITE(*,*), "Ice forecast variables loaded"

!!!!---------------------------------------------------------------------!!!

!#################!
!# Open analysis #!
!#################!

ierr = nf90_open(trim(cfile_analysis_ice),nf90_NoWrite,incid_ice_an_in)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  STOP
END IF

! Time counter
! Request variable id
cvarname = 'time_counter'
ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_ice_an_in, ivarid, ztime_counter )
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF


! Strategy: Loop over categories, and for specific variables over layers
! (2024) Categories as dim now!
!DO jl = 1 , jpl
  ! 3.D.1. Ice concentration
  cvarroot='a_i'
  cvarname = TRIM(cvarroot)
!  WRITE(*,*), "Working with variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Get variable
  ierr = nf90_get_var(incid_ice_an_in, ivarid, a_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! 3.D.2. Ice volume per surface area
  cvarroot='v_i'
  cvarname = TRIM(cvarroot)
!  WRITE(*,*), "Working with variable " // cvarname

  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Get variable
  ierr = nf90_get_var(incid_ice_an_in, ivarid, v_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! 3.D.3. Snow volume per surface area
  cvarroot='v_s'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Working with variable " // cvarname

  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Get variable
  ierr = nf90_get_var(incid_ice_an_in, ivarid, v_s)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! 3.D.4. Ice age
  cvarroot='oa_i'
  cvarname = TRIM(cvarroot)
!  WRITE(*,*), "Working with variable " // cvarname
    
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Get variable
  ierr = nf90_get_var(incid_ice_an_in, ivarid, oa_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! 3.D.5. Ice salinity
  cvarroot='sv_i'
  cvarname = TRIM(cvarroot) 
!  WRITE(*,*), "Working with variable " // cvarname

  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! Get variable
  ierr = nf90_get_var(incid_ice_an_in, ivarid, smv_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
    WRITE(*,*), TRIM(cfile_analysis_ice)
    WRITE(*,*), cvarname
    STOP
  END IF

  ! 3.D.5. Ice surface temperature
  !cvarroot='t_su_htc'
  !WRITE(cinteger,'(i1)') jl
  !cvarname = TRIM(cvarroot) // TRIM(cinteger)
! ! WRITE(*,*), "Working with variable " // cvarname

  !! Request variable ID
  !ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
  !IF (ierr .NE. nf90_noerr ) THEN
  !  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
  !  WRITE(*,*), TRIM(cfile_analysis_ice)
  !  WRITE(*,*), cvarname
  !  STOP
  !END IF

  !! Get variable
  !ierr = nf90_get_var(incid_ice_an_in, ivarid, t_su(:,:,jl) )
  !IF (ierr .NE. nf90_noerr ) THEN
  !  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
  !  WRITE(*,*), TRIM(cfile_analysis_ice)
  !  WRITE(*,*), cvarname
  !  STOP
  !END IF

cvarroot='t_su'
cvarname = TRIM(cvarroot)
 WRITE(*,*), "Working with variable " // cvarname

! Request variable ID
ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_ice_an_in, ivarid, t_su)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

  ! 3.D.X Variables that are defined on several ice layers
  DO jk = 1 , nlay_i
    ! Ice enthalpy in layers
    cvarroot='e_i_l0'
    WRITE(cinteger2,'(i1)') jk 
    cvarname = TRIM(cvarroot)//TRIM(cinteger2)
!    WRITE(*,*),"Working with variable " // cvarname
    
    ! Request variable ID
    ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
      WRITE(*,*), TRIM(cfile_analysis_ice)
      WRITE(*,*), cvarname
      STOP
    END IF

    ! Get variable
    ierr = nf90_get_var(incid_ice_an_in, ivarid, e_i(:,:,jk,:) )
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
      WRITE(*,*), TRIM(cfile_analysis_ice)
      WRITE(*,*), cvarname
      STOP
    END IF 
  ENDDO ! jk, layers 

  DO jk = 1 , nlay_s
    ! Snow enthalpy in layers
    cvarroot='e_s_l0'
    WRITE(cinteger2,'(i1)') jk
    cvarname = TRIM(cvarroot)// TRIM(cinteger2)
!    WRITE(*,*),"Working with variable " // cvarname

    ! Request variable ID
    ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
      WRITE(*,*), TRIM(cfile_analysis_ice)
      WRITE(*,*), cvarname
      STOP
    END IF

    ! Get variable
    ierr = nf90_get_var(incid_ice_an_in, ivarid, e_s(:,:,jk,:) )
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
      WRITE(*,*), TRIM(cfile_analysis_ice)
      WRITE(*,*), cvarname
      STOP
    END IF
  ENDDO ! jk, layers
!ENDDO ! jl, categories

! Load data that don't depend on categories

WRITE(*,*), "Ice analyse 3-D variables loaded (category and layers)"

! Snow ice mass
cvarname = 'snwice_mass'
! Request variable ID
ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_ice_an_in, ivarid, snwice_mass(:,:) )
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Snow ice mass_b
cvarname = 'snwice_mass_b'
! Request variable ID
ierr = nf90_inq_varid(incid_ice_an_in, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Get variable
ierr = nf90_get_var(incid_ice_an_in, ivarid, snwice_mass_b(:,:) )
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice file get. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  WRITE(*,*), cvarname
  STOP
END IF

! Close ice analysis file
ierr = nf90_close(incid_ice_an_in);
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ice close file. Abort"
  WRITE(*,*), TRIM(cfile_analysis_ice)
  STOP
END IF

WRITE(*,*), "Ice forecast without category dim variables loaded"

WRITE(*,*) 'Leaving load_variables'
END SUBROUTINE load_variables

SUBROUTINE several_ice_corrections()
  USE my_variables
  IMPLICIT NONE

  !  Dummy variables
  INTEGER                 :: ji, jj, jl, jk

  WRITE(*,*), 'Entering several_ice_corrections()'

  DO ji = 1 , jpi
    DO jj = 1 , jpj
      DO jl = 1 , jpl
        ! Define switches (masks)
        !zindb =  MAX( rzero , SIGN( zamax , a_i(ji,jj,jl) - epsi06) )
        !zindsn = MAX(rzero ,  SIGN( zamax , v_s(ji,jj,jl) - epsi06) )
        !zindic = MAX(rzero ,  SIGN( zamax , v_i(ji,jj,jl) - epsi04) )
        zindb =  MAX( rzero , SIGN( rone , a_i(ji,jj,jl) - epsi06) )
        zindsn = MAX(rzero ,  SIGN( rone , v_s(ji,jj,jl) - epsi06) )
        zindic = MAX(rzero ,  SIGN( rone , v_i(ji,jj,jl) - epsi04) )

        zindb = zindb * zindic ! Mask where there is  conc. AND volume 

        ! A. Corrections to ice age
        ![AD] test without modifing age, so analysis = forecast
        !IF ( ( oa_i(ji,jj,jl) - rone ) * 86400.0 .GT. ( rdt_ice*ztime_counter*a_i(ji,jj,jl) ) )&
        !&THEN
          !WRITE(*,*) 'Correction on ice age'
        !  oa_i(ji,jj,jl) = rdt_ice * ztime_counter / 86400.0 * a_i(ji,jj,jl)
        !END IF
        !oa_i(ji,jj,jl) = zindb*oa_i(ji,jj,jl)
        
        ! B. Corrections to snow thickness
        v_s(ji,jj,jl) = zindsn*zindb*v_s(ji,jj,jl)

        ! C. Corrections to ice thickness
        v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl) , rzero ) 
        v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)
        

        ! D. Snow transformed to ice  if original ice cover disappears
        zindg = REAL(ilandmask(ji,jj) ) * MAX (rzero, SIGN( rone , - v_i (ji,jj,jl)  ) )
        ! (is it possible to have zindg not zero??)
        v_i(ji,jj,jl) = v_i(ji,jj,jl) + zindg * zrhosn * v_s(ji,jj,jl) / zrho0
        ! The last term of the above equation is the water volume equivalent to the snow
        ! volume I guess
        v_s(ji,jj,jl) = (rone - zindg ) * v_s(ji,jj,jl) + &
        & zindg * v_i(ji,jj,jl) * (zrho0 - zrhoic ) / zrhosn

        ! E. Correction to ice concentration
        a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX(a_i(ji,jj,jl), epsi06)

        ! F. Correction to snow heat contenti
        ![AD] Add loop on layer also for e_s to test! 
        IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN
          WRITE(*,*) 'Before i=307 j=307 e_s_l01:' , e_s(ji,jj,1,jl)  
        END IF

        DO jk = 1 , nlay_s
           e_s(ji,jj,jk,jl) = zindsn * &
           & ( MIN ( MAX ( rzero, e_s(ji,jj,jk,jl) ), zbigvalue ) ) + &
           & ( rone - zindsn ) * rzero
        END DO

         IF ( ldebug .AND. ( ji == jiindx ) .AND. ( jj == jjindx ) ) THEN
           WRITE(*,*) 'After i=307 j=307 e_s_l010:' , e_s(ji,jj,1,jl)
         END IF

        ! G. Correction to ice heat content
        DO jk = 1 , nlay_i
          e_i(ji,jj,jk,jl) = zindic * &
          &( MIN( MAX( rzero, e_i(ji,jj,jk,jl) ), zbigvalue ) ) &
          & + (rone - zindic ) * rzero  
        END DO ! nlay_i

        ! H. Correction to ice salinity
        IF ( (num_sal .EQ. 2) .OR. (num_sal .EQ. 4) ) THEN
          ! If we are in salinity profile mode
          ! Salinity stays in bounds
          smv_i(ji,jj,jl) = MAX( MIN( (zrhoic-zrhosn)/zrhoic * &
          & sss_m_an(ji,jj) , smv_i(ji,jj,jl) ) , & 
            0.1 * v_i(ji,jj,jl) )
        ENDIF 
         
        ! I. Thickness of first category above 5cm
        IF ( jl == 1) THEN
          ht_i(ji,jj,jl) = zindb * v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi06)
          zh             = MAX(rone , zindb*zhiclim/ & 
                           & MAX(ht_i(ji,jj,jl),epsi20 ) )
          ! This is a correction factor equal to zhiclim/h_insitu, i.e. the ratio
          ! between minimal and actual in situ thickness. 
        
          ! Something to do for the snow
 
          ! The ice concentration is shrunk so that the ice thickness is at least
          ! zhiclim
          a_i(ji,jj,jl) = a_i(ji,jj,jl) / zh
        END IF
      
    END DO ! jl

    ! Reset snowice to zero if necessary
    IF ( (snwice_mass(ji,jj) .LT. rzero) .OR. (snwice_mass_b(ji,jj) .LT. rzero) ) THEN
      snwice_mass(ji,jj) = rzero
      snwice_mass_b(ji,jj)=rzero
    END IF
  END DO ! jj
END DO ! ji

  WRITE(*,*), 'Leaving several_ice_corrections()'
END SUBROUTINE several_ice_corrections

SUBROUTINE shrink_ice()
USE my_variables
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: ji, jj, jl

WRITE(*,*), 'Entering shrink_ice()'
!  Total concentration cannot exceed zamax
at_i(:,:) = rzero

DO jl = 1 , jpl
  at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
END DO

DO ji = 1 , jpi
  DO jj = 1 , jpj
    ! Define the excessive concentration
    zda_ex = MAX( at_i(ji,jj) - zamax , rzero )
    ! (i.e. rzero, except if excess)
    DO jl = 1 , jpl
      zindb = MAX( rzero , SIGN( rone , v_i(ji,jj,jl) ) )
      ! (zindb is a mask)
      ! Spread the excess over the different categories with weight equal
      ! to concentration in each of them
      zda_i = a_i(ji,jj,jl) * zda_ex / MAX(at_i(ji,jj),epsi06) * zindb
      ! Correction of limupdate after Antoine Barthélemy
      ! Indeed the volume should not be changed.
      !! zdv_i = v_i(ji,jj,jl) * zda_i  / MAX(at_i(ji,jj),epsi06)
      ! We take out the excess of ice and put it as volume
      a_i(ji,jj,jl) = a_i(ji,jj,jl) - zda_i
      !! v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdv_i 
    END DO
  END DO
END DO
WRITE(*,*), 'Leaving shrink_ice()'

END SUBROUTINE shrink_ice

SUBROUTINE record_ice()
USE my_variables
USE netcdf
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: jl, jk

WRITE(*,*) 'Entering record_ice()'
WRITE(*,*) 'Recording the NetCDF'


! 8.1 Record ice variables
! We copy the input file and store everything into this copy
CALL system('cp -f '//trim(cfile_analysis_ice)//' '//trim(cfileout_ice) )

ierr = nf90_open(trim(cfileout_ice),nf90_Write,incid_ice_an_out)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF

! Loop over categories
! (2024)  Categories as dimension now, no need for the loop anymore. 
!         I also change the var names to correspond with the current nemo restarts. 
!DO jl = 1 , jpl
  ! 8.1.1 a_i
  !~~~~~~~~
  cvarroot='a_i'
  !WRITE(cinteger,'(i1)') jl
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, a_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! 8.1.2 v_i
  !~~~~~~~~
  cvarroot='v_i'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, v_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF
  
  ! 8.1.3 v_s
  !~~~~~~~~
  cvarroot='v_s'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, v_s)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF
  
  ! 8.1.4 smv_i
  !~~~~~~~~~~
  cvarroot='sv_i'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, smv_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. smv_i put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! 8.1.5 oa_i
  ! ~~~~~~~~
  cvarroot='oa_i'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, oa_i)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! 8.1.6 t_su
  ! ~~~~~~~~
  cvarroot='t_su'
  cvarname = TRIM(cvarroot)
  WRITE(*,*), "Recording variable " // cvarname
  ! Request variable ID
  ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! Put variable 
  ierr = nf90_put_var(incid_ice_an_out, ivarid, t_su)
  IF (ierr .NE. nf90_noerr ) THEN
    WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
    WRITE(*,*), TRIM(cfileout_ice)
    STOP
  END IF

  ! 8.1.7 Ice enthalpy (layers)
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~
  cvarroot='e_i_l0'
  DO jk = 1 , nlay_i
    WRITE(cinteger2,'(i1)') jk
    cvarname = TRIM(cvarroot) // TRIM(cinteger2)
    WRITE(*,*),"Recording variable " // cvarname
    ! Request variable ID
    ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
      WRITE(*,*), TRIM(cfileout_ice)
      STOP
    END IF

    ! Put variable 
    ierr = nf90_put_var(incid_ice_an_out, ivarid, e_i(:,:,jk,:))
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. e_i put. Abort"
      WRITE(*,*), TRIM(cfileout_ice)
      STOP
    END IF
  END DO !jk

  ! 8.1.8 Snow enthalpy (layers)
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~
  cvarroot='e_s_l0'
  DO jk = 1 , nlay_s
  WRITE(cinteger2,'(i1)') jk  
  cvarname = TRIM(cvarroot) // TRIM(cinteger2)
    WRITE(*,*),"Recording variable " // cvarname
    ! Request variable ID
    ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
      WRITE(*,*), TRIM(cfileout_ice)
      STOP
    END IF

    ! Put variable
    ierr = nf90_put_var(incid_ice_an_out, ivarid, e_s(:,:,jk,:))
  
    IF (ierr .NE. nf90_noerr ) THEN
      WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. e_s put. Abort"
      WRITE(*,*), TRIM(cfileout_ice)
      STOP
    END IF
  END DO !jk

!END DO ! jl, categories

! Put variables that don't depend on categories

! Snow ice mass
! Request variable id
cvarname = 'snwice_mass'
ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF

! Put variable
ierr = nf90_put_var(incid_ice_an_out, ivarid, snwice_mass(:,:))

IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF

cvarname = 'snwice_mass_b'
ierr = nf90_inq_varid(incid_ice_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF

! Put variable
ierr = nf90_put_var(incid_ice_an_out, ivarid, snwice_mass_b(:,:))

IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF ice file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF



! Close the NetCDF file
ierr = nf90_close(incid_ice_an_out)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort"
  WRITE(*,*), TRIM(cfileout_ice)
  STOP
END IF

WRITE(*,*) 'Leaving record_ice()'
END SUBROUTINE record_ice

SUBROUTINE record_oce()
USE my_variables
USE netcdf
IMPLICIT NONE
WRITE(*,*) 'Entering record_oce()'

! Record oce variables
! We copy the input file and store everything into this copy
CALL system('cp -f '//trim(cfile_analysis_oce)//' '//trim(cfileout_oce) )

! Open the file
ierr = nf90_open(trim(cfileout_oce),nf90_Write,incid_oce_an_out)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF Ocean file. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! A. sn (salinity)
cvarname= csn
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)

IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, sn_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Following lines commented out as the variable no longer appears
! in the restart files (2020)

!! B. sss_m (sea surface salinity)
!cvarname= csss_m
!ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF

!! Put variable 
!ierr = nf90_put_var(incid_oce_an_out, ivarid, sss_m_an(:,:))
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF

! C. tn (temperature)
cvarname= ctn
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)

IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, tn_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Following lines commented out as the variable no longer appears
! in the restart files (2020)

!! D. sst_m (sea surface temperature)
!cvarname= csst_m
!ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF

!! Put variable 
!ierr = nf90_put_var(incid_oce_an_out, ivarid, sst_m_an(:,:))
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF

! E. un (sea velocity, "un")
cvarname= cun
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, un_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! F. ub (sea velocity, "ub")
cvarname= cub
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, ub_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF


! G. vn (sea velocity, "vn")
cvarname= cvn
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, vn_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! H. vb (sea velocity, "vb")
cvarname= cvb
ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Put variable 
ierr = nf90_put_var(incid_oce_an_out, ivarid, vb_an(:,:,:))
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF

! Following lines commented out as the variable no longer appears
! in the restart files (2020)

!! I. ssu_m (sea surface velocity, "ssu_m")
!cvarname= cssu_m
!ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF
!
!! Put variable 
!ierr = nf90_put_var(incid_oce_an_out, ivarid, ssu_m_an(:,:))
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF
!
!! J. ssv_m (sea surface velocity, "ssv_m")
!cvarname= cssv_m
!ierr = nf90_inq_varid(incid_oce_an_out, cvarname, ivarid)
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. inquiry. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF
!
!! Put variable 
!ierr = nf90_put_var(incid_oce_an_out, ivarid, ssv_m_an(:,:))
!IF (ierr .NE. nf90_noerr ) THEN
!  WRITE(*,*), "(sanity_check_LIM3) Error NetCDF oce file var. put. Abort"
!  WRITE(*,*), TRIM(cfileout_oce)
!  STOP
!END IF

! Close file
ierr = nf90_close(incid_oce_an_out)
IF (ierr .NE. nf90_noerr ) THEN
  WRITE(*,*), "(sanity_check_LIM3) Error Closing NetCDF . Abort"
  WRITE(*,*), TRIM(cfileout_oce)
  STOP
END IF


WRITE(*,*) 'Leaving record_oce'
END SUBROUTINE record_oce

SUBROUTINE salinity_check()
USE my_variables
USE netcdf
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: ji, jj, jk

REAL                    :: zmaxsaldiff=9999.0


WRITE(*,*) 'Entering salinity_check() '

DO ji=1,jpi
  DO jj=1,jpj
    DO jk=1,jpk
      zsaldiff=sn_an(ji,jj,jk) - sn_fc(ji,jj,jk) 
      IF (ABS(zsaldiff) .GT. zmaxsaldiff  ) THEN

        IF (ABS(zsaldiff) .GT. 100.0 * zmaxsaldiff  ) THEN
          WRITE(*,*) "Very large salinity variation"
          WRITE(*,*) "at point ji,jj,jk = ", ji,jj,jk
          WRITE(*,*) "sn_an is ",sn_an(ji,jj,jk)
          WRITE(*,*) "sn_fc is ",sn_fc(ji,jj,jk)
          WRITE(*,*) "diff is", zsaldiff
        END IF
        !WRITE(*,*) "Salinity difference at ji,jj,jk=",ji,jj,jk
        !WRITE(*,*) "Difference is ",ABS(zsaldiff)
        !WRITE(*,*) "sn_an is ",sn_an(ji,jj,jk)
        !WRITE(*,*) "sn_fc is ",sn_fc(ji,jj,jk)
          
        sn_an(ji,jj,jk) = sn_fc(ji,jj,jk) + SIGN(zmaxsaldiff,zsaldiff)
        !WRITE(*,*) "sn_an is now ",sn_an(ji,jj,jk)
      END IF
    END DO ! jk

      zsaldiff=(sss_m_an(ji,jj) - sss_m_fc(ji,jj)) / REAL( nn_fsbc - 1 )

      IF (ABS(zsaldiff) .GT. zmaxsaldiff  ) THEN
        IF (ABS(zsaldiff) .GT. 100.0 * zmaxsaldiff  ) THEN
          WRITE(*,*) "Very large SSS variation  at ji,jj  = ",ji,jj
          WRITE(*,*) "sss_m_an is ", sss_m_an(ji,jj)
          WRITE(*,*) "sss_m_fc is ", sss_m_fc(ji,jj)
          WRITE(*,*) "diff is"     , zsaldiff
        END IF
        !WRITE(*,*) "Sea surface salinity difference at ji,jj=",ji,jj
        !WRITE(*,*) "Difference is ",ABS(zsaldiff)
        !WRITE(*,*) "sss_m_an is ",sss_m_an(ji,jj) / REAL( nn_fsbc - 1 )
        !WRITE(*,*) "sss_m_fc is ",sss_m_fc(ji,jj) / REAL( nn_fsbc - 1 )

        sss_m_an(ji,jj) = sss_m_fc(ji,jj) + SIGN(REAL( nn_fsbc - 1 ) * zmaxsaldiff , zsaldiff)

        !WRITE(*,*) "sss_m_an is now ",sss_m_an(ji,jj) / REAL( nn_fsbc - 1 )

        ! Important note
        ! The reason for the nn_fsbc - 1 factor is the following. In the
        ! restarts, the variable sss_m is (nn_fsbc -1) times the first layer of
        ! the sea salinity, where nn_fsbc is the frequency call of the ice to
        ! the ocean. As I understood, this is intended to facilitate the switch
        ! from one frequency call to another. The point is that the variable
        ! sss_m is the cumulative sea surface salinity over the (nn_fsbc-1) time
        ! steps. Stated the other way around, it's nn_fsbc-1 times its mean
        ! value. Hence here we divide sss_m by (nn_fsbc-1), ensure that
        ! variations in salinity are not too strong, and give back the corrected
        ! quantity multiplied by (nn_fsbc-1)
      END IF
  END DO
END DO
WRITE(*,*) 'Leaving salinity_check()'
END SUBROUTINE salinity_check

SUBROUTINE temperature_check()
USE my_variables
USE netcdf
IMPLICIT NONE
!  Dummy variables
INTEGER                 :: ji, jj, jk

REAL                    :: zmaxtempdiff=100.0
REAL                    :: ztempmin           ! Minimum temperature for surf sea water (fct of salinity)

WRITE(*,*) 'Entering temperature_check() '

! 3D- temperature
DO ji=1,jpi
  DO jj=1,jpj
    DO jk=1,jpk
      ztempdiff=tn_an(ji,jj,jk) - tn_fc(ji,jj,jk) 
      IF (ABS(ztempdiff) .GT. zmaxtempdiff  ) THEN

        IF (ABS(ztempdiff) .GT. 100.0 * zmaxtempdiff  ) THEN
          WRITE(*,*) "Very large temperature variation detected!!"
          WRITE(*,*) "at point ji,jj,jk = ", ji,jj,jk
          WRITE(*,*) "tn_an is ",tn_an(ji,jj,jk)
          WRITE(*,*) "tn_fc is ",tn_fc(ji,jj,jk)
        END IF
        !WRITE(*,*) "Temperature difference at ji,jj,jk=",ji,jj,jk
        !WRITE(*,*) "Difference is ",ABS(ztempdiff)
        !WRITE(*,*) "tn_an is ",tn_an(ji,jj,jk)
        !WRITE(*,*) "tn_fc is ",tn_fc(ji,jj,jk)
          
        tn_an(ji,jj,jk) = tn_fc(ji,jj,jk) + SIGN(zmaxtempdiff,ztempdiff)
        
        !WRITE(*,*) "tn_an is now ",tn_an(ji,jj,jk)
      END IF ! if diff is large

    END DO ! jk

    DO jk =1,1
      ! Reset surface temperature to freezing point if below freezing point
      ! This depends on the formulation chosen in the namelist (nn_eos).
      
      ! The formula comes from the NEMO code (routine eosbn2.f90)
      !
      ! In the case nn_eos = 0 (UNESCO formulation)
      ! ztempmin=  ( - 0.0575_wp + 1.710523e-3_wp * SQRT(  sn_an(ji,jj,jk) )   &
      !       &                     - 2.154996e-4_wp *    sn_an(ji,jj,jk)   ) * sn_an(ji,jj,jk) 
      ! 
      ! In the case nn_eos = -1 or 1 (TEOS-10 formulation)
      r1_S0 = 0.875_wp/35.16504_wp
 
      zs = sqrt( abs(sn_an(ji,jj,jk)) * r1_S0)

      ztempmin = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp

      ztempmin = ztempmin * sn_an(ji,jj,jk)
      
      ! This is the freezing point of sea water at the surface.
      IF ( tn_an(ji,jj,jk) .LT. ztempmin ) THEN
        WRITE(*,*) 'At grid point ', ji,jj,jk
        WRITE(*,*) 'tn_an reset from', tn_an(ji,jj,jk),'to ', ztempmin
        tn_an(ji,jj,jk) = ztempmin
      END IF

    END DO ! jk

! SST 
      ztempdiff=(sst_m_an(ji,jj) - sst_m_fc(ji,jj)) / REAL( nn_fsbc - 1 )

      IF (ABS(ztempdiff) .GT. zmaxtempdiff  ) THEN
        IF (ABS(ztempdiff) .GT. 100.0 * zmaxtempdiff  ) THEN
          WRITE(*,*) "Very large SST variation deteced at ji,jj  = ",ji,jj
          WRITE(*,*) "sst_m_an is ", sst_m_an(ji,jj)
          WRITE(*,*) "sst_m_fc is ", sst_m_fc(ji,jj)
          WRITE(*,*) "diff is"     , ztempdiff
        END IF
        !WRITE(*,*) "Sea surface temperature difference at ji,jj=",ji,jj
        !WRITE(*,*) "Difference is ",ABS(ztempdiff)
        !WRITE(*,*) "sst_m_an is ",sst_m_an(ji,jj) / REAL( nn_fsbc - 1 )
        !WRITE(*,*) "sst_m_fc is ",sst_m_fc(ji,jj) / REAL( nn_fsbc - 1 )

        sst_m_an(ji,jj) = sst_m_fc(ji,jj) + SIGN(REAL( nn_fsbc - 1 ) * zmaxtempdiff , ztempdiff)

        !WRITE(*,*) "sst_m_an is now ",sst_m_an(ji,jj) / REAL( nn_fsbc - 1 )

        ! Important note
        ! The reason for the nn_fsbc - 1 factor is the following. In the
        ! restarts, the variable sst_m is (nn_fsbc -1) times the first layer of
        ! the sea temperature, where nn_fsbc is the frequency call of the ice to
        ! the ocean. As I understood, this is intended to facilitate the switch
        ! from one frequency call to another. The point is that the variable
        ! sst_m is the cumulative sea surface temperature over the (nn_fsbc-1) time
        ! steps. Stated the other way around, it's nn_fsbc-1 times its mean
        ! value. Hence here we divide sst_m by (nn_fsbc-1), ensure that
        ! variations in temperature are not too strong, and give back the corrected
        ! quantity multiplied by (nn_fsbc-1)
      END IF
   
      
      ! For nn_eos = 0
      ! ztempmin=  ( - 0.0575_wp + 1.710523e-3_wp * SQRT(  sn_an(ji,jj,1) )   &
      !        &                     - 2.154996e-4_wp *    sn_an(ji,jj,1)   ) * sn_an(ji,jj,1) 

      ! For nn_eos = -1 or 1
      r1_S0 = 0.875_wp/35.16504_wp

      zs = sqrt( abs(sn_an(ji,jj,1)) * r1_S0)

      ztempmin = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp

      ztempmin = ztempmin * sn_an(ji,jj,1)
       

      IF ( ( sst_m_an(ji,jj) / REAL( nn_fsbc - 1 ) ) .LT. ztempmin ) THEN
        WRITE(*,*), 'At grid point ', ji,jj
        WRITE(*,*), 'SST reset from ',sst_m_an(ji,jj)/REAL( nn_fsbc - 1 ),' to ' ,ztempmin
        sst_m_an(ji,jj) = ztempmin * REAL (nn_fsbc - 1) 
      END IF

  END DO ! jj
END DO ! ji
WRITE(*,*) 'Leaving temperature_check()'
END SUBROUTINE temperature_check

SUBROUTINE  update_hc()
  USE my_variables
  IMPLICIT NONE

  INTEGER                 :: ji, jj, jk, jl
  REAL(wp) :: zratio, ztmelts
  REAL(wp), PARAMETER :: t_init = 270.0  ! Initial temperature of ice when it is forming. Inspired from limistate.
  REAL(wp)            :: zhc             ! Heat content (in ice or snow)
  REAL(wp), PARAMETER :: zunit_fac= 1.0e9! This 1.0e9 is because the e_i and e_s variables are saved in Giga Joules / m2 in
                                         ! the restart (but in Joules/m2 in the code), I guess because the restart cannot take
                                         ! numbers with large exponents.
                                         ! For info: the energy to melt 1 meter of ice at 0°C is
                                         ! 334 000 [J / kg] * 1 [m] * 1000 [kg/m3] = 0.334 x 10^9 J / m2
                                         

  WRITE(*,*) 'Entering update_hc()'
  
  DO jl = 1, jpl
    DO ji = 1, jpi
      DO jj = 1, jpj
        ! Ice layers
        ! Case 1: there was ice in the forecast
        IF (v_i_fc(ji,jj,jl) .GT. epsi10) THEN
            zratio = v_i(ji,jj,jl) / v_i_fc(ji,jj,jl)
            ! Update the ice heat content by that amount in each layer
            DO jk = 1, nlay_i
              e_i(ji,jj,jk,jl) = zratio * e_i(ji,jj,jk,jl)

            END DO ! jk

        ! Case 2: there was no ice in the forecast
        ELSEIF (v_i(ji,jj,jl) .GT. epsi06 ) THEN
           ztmelts = - tmut * smv_i(ji,jj,jl) + rtt
             zhc          = zrhoic *   (&                               ! [kg/m3]
                                        &zcpic * (ztmelts - t_init)&   ! [J/(kg.K) ] * [K ] = J/kg
                                        &+zlfus* (1- (ztmelts - rtt)/(t_init-rtt) )&  ! [J/kg]*[] 
                                        &-zrcp * (ztmelts - rtt)& ! [J/(kg.K)] * [K]
                                        &)
             ! zhc is in J/m3
             ! The amount of heat in each layer is that divided by the number of
             ! layers, multiplied by the sea ice volume (v_i*cell area) and by
             ! unit_fac (see LIM3 code) to get the units in GigaJoule
             DO jk = 1,nlay_i
                e_i(ji,jj,jk,jl) = zhc * v_i(ji,jj,jl) * zcellarea(ji,jj) / nlay_i / zunit_fac
             END DO
             !WRITE(*,*) "Ice was created by the filter at point ", ji,jj
             !WRITE(*,*) "Ice volume in forecast was",v_i_fc(ji,jj,jl)
             !WRITE(*,*) "Ice volume in analysis is ",v_i(ji,jj,jl)
             !WRITE(*,*) "In category ", jl
             !WRITE(*,*) "New e_i is ",e_i(ji,jj,:,jl)
        END IF
             !IF ( ji ==127 .AND. jj == 124 .AND. jl==5 ) THEN
             !   WRITE(*,*) "ji = ",ji,"jj = ",jj, "jl =", jl
             !   WRITE(*,*) "e_i: ", e_i(ji,jj,:,jl)
             !   WRITE(*,*) "zratio: ",zratio
             !   WRITE(*,*) "v_i", v_i(ji,jj,jl)
             !   WRITE(*,*) "v_i_fc: ", v_i_fc(ji,jj,jl)
             !   WRITE(*,*) "Stop because requested so"
             !   !STOP
             !END IF
 
        ! Snow layer
        ! Case 1: there was snow in the forecast
        IF (v_s_fc(ji,jj,jl) .GT. epsi06) THEN
            zratio = v_s(ji,jj,jl) / v_s_fc(ji,jj,jl)
            ! Update the ice heat content by that amount in each layer
            DO jk = 1, nlay_s
              e_s(ji,jj,jk,jl) = zratio * e_s(ji,jj,jk,jl)
            END DO ! jk

        ! Case 2: there was no snow in the forecast
        ELSEIF (v_s(ji,jj,jl) .GT. epsi06) THEN
             zhc          = zrhosn *   (&                               ! [kg/m3]
                                        &zcpic * (rtt - t_init)&   ! [J/(kg.K) ] * [K ] = J/kg
                                        &+zlfus& ! [J/(kg)] 
                                        &) 
             ! zhc is in J/m3
             ! The amount of heat in each layer is that divided by the number of
             ! layers, multiplied by the snow volume (v_s*cell area)
             DO jk = 1,nlay_s
                e_s(ji,jj,jk,jl) = zhc * v_s(ji,jj,jl) * zcellarea(ji,jj) / nlay_s / zunit_fac
             END DO
        END IF
      END DO! jj
    END DO ! ji
  END DO !jl

  WRITE(*,*) 'Leaving update_hc()'
END SUBROUTINE update_hc

SUBROUTINE regularize()
  USE my_variables
  IMPLICIT NONE

  INTEGER                 :: ji, jj, jk, jl

  ! Resets < 0 concentrations to 0
  ! Resets variables to zero where no ice concentration

  DO ji=1,jpi
    DO jj=1,jpj
      DO jl =1,jpl
        IF ( a_i(ji,jj,jl) .LT. rzero ) THEN
                a_i(  ji,jj,jl) = rzero
                v_i(  ji,jj,jl) = rzero
                smv_i(ji,jj,jl) = rzero
                v_s(  ji,jj,jl) = rzero
                oa_i( ji,jj,jl) = rzero
                DO jk=1,nlay_i
                        e_i(ji,jj,jk,jl) = rzero
                END DO
                DO jk=1,nlay_s
                        e_s(ji,jj,jk,jl) = rzero
                END DO
        END IF
      END DO
    
      ! Variables that do not depend on categories
      IF (  SUM(a_i(ji,jj,:)) .LT. rzero  ) THEN
        snwice_mass(ji,jj) = rzero
        snwice_mass_b(ji,jj) = rzero
      END IF

    END DO!jj
  END DO ! ji
END SUBROUTINE regularize

SUBROUTINE velocity_check()
  USE my_variables
  IMPLICIT NONE
  INTEGER                :: ji,jj,jk
  
  REAL                   :: zmaxvel = 2.0
  REAL                   :: zmaxsurfvel= 6.0  
  ! Resets ocean velocities to [-2,2] ms-1
  ! Resets sea surface velocities to [-4,4] ms-1
  DO ji=1,jpi
    DO jj = 1,jpj
      DO jk=1,jpk
     
      IF ( ABS(un_an(ji,jj,jk)) .GT. zmaxvel ) THEN
        WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk
        WRITE(*,*) "un_an is ",un_an(ji,jj,jk)
        !un_an(ji,jj,jk) =  un_fc(ji,jj,jk) !zmaxvel * SIGN( rone , un_an(ji,jj,jk) )
        un_an(ji,jj,jk) =   zmaxvel * SIGN( rone , un_an(ji,jj,jk) )
        WRITE(*,*) "un_an reset to",un_an(ji,jj,jk)
      END IF

      IF ( ABS(ub_an(ji,jj,jk)) .GT. zmaxvel ) THEN
        WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk
        WRITE(*,*) "ub_an is ",ub_an(ji,jj,jk)
        !ub_an(ji,jj,jk) =  ub_fc(ji,jj,jk) !zmaxvel * SIGN( rone , ub_an(ji,jj,jk) )
        ub_an(ji,jj,jk) =  zmaxvel * SIGN( rone , ub_an(ji,jj,jk) )
        WRITE(*,*) "ub_an reset to",ub_an(ji,jj,jk)
      END IF

      IF ( ABS(vn_an(ji,jj,jk)) .GT. zmaxvel ) THEN
        WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk
        WRITE(*,*) "vn_an is ",vn_an(ji,jj,jk)
        !vn_an(ji,jj,jk) =  vn_fc(ji,jj,jk) !zmaxvel * SIGN( rone , vn_an(ji,jj,jk) )
        vn_an(ji,jj,jk) =  zmaxvel * SIGN( rone , vn_an(ji,jj,jk) )
        WRITE(*,*) "vn_an reset to",vn_an(ji,jj,jk)
      END IF

      IF ( ABS(vb_an(ji,jj,jk)) .GT. zmaxvel ) THEN
        WRITE(*,*) "Fast ocean found at ji,jj,jk=",ji,jj,jk
        WRITE(*,*) "vb_an is ",vb_an(ji,jj,jk)
        !vb_an(ji,jj,jk) =  vb_fc(ji,jj,jk) !zmaxvel * SIGN( rone , vb_an(ji,jj,jk) )
        vb_an(ji,jj,jk) =  zmaxvel * SIGN( rone , vb_an(ji,jj,jk) )
        WRITE(*,*) "vb_an reset to",vb_an(ji,jj,jk)
      END IF

      END DO !jk

      ! Surface velocity
      IF (ABS(ssu_m_an(ji,jj)) .GT. zmaxsurfvel ) THEN
        WRITE(*,*) "Fast sea surface velocity found at ji,jj=" , ji,jj
        WRITE(*,*) "ssu_m_an is ",ssu_m_an(ji,jj)
        !ssu_m_an(ji,jj) = ssu_m_fc(ji,jj) !zmaxsurfvel * SIGN(rone, ssu_m_an(ji,jj) )
        ssu_m_an(ji,jj) = zmaxsurfvel * SIGN(rone, ssu_m_an(ji,jj) )
        WRITE(*,*) "ssu_m_an reset to", ssu_m_an(ji,jj)
      END IF 

      IF (ABS(ssv_m_an(ji,jj)) .GT. zmaxsurfvel ) THEN
        WRITE(*,*) "Fast sea surface velocity found at ji,jj=" , ji,jj
        WRITE(*,*) "ssv_m_an is ",ssv_m_an(ji,jj)
        !ssv_m_an(ji,jj) = ssv_m_fc(ji,jj) !zmaxsurfvel * SIGN(rone, ssv_m_an(ji,jj) )
        ssv_m_an(ji,jj) = zmaxsurfvel * SIGN(rone, ssv_m_an(ji,jj) )
        WRITE(*,*) "ssv_m_an reset to", ssv_m_an(ji,jj)
      END IF

    END DO
  END DO
END SUBROUTINE velocity_check