Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/dtatsd.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 14923)
@@ -6,6 +6,6 @@
    !! History :  OPA  ! 1991-03  ()  Original code
    !!             -   ! 1992-07  (M. Imbard)
-   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT 
-   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
+   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
    !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
    !!            3.4  ! 2010-11  (G. Madec, C. Ethe) Merge of dtatem and dtasal + remove CPP keys
@@ -49,9 +49,9 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE dta_tsd_init  ***
-      !!                    
-      !! ** Purpose :   initialisation of T & S input data 
-      !! 
+      !!
+      !! ** Purpose :   initialisation of T & S input data
+      !!
       !! ** Method  : - Read namtsd namelist
-      !!              - allocates T & S data structure 
+      !!              - allocates T & S data structure
       !!----------------------------------------------------------------------
       LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used
@@ -77,5 +77,5 @@
 
       IF( PRESENT( ld_tradmp ) )   ln_tsd_dmp = .TRUE.     ! forces the initialization when tradmp is used
-      
+
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
@@ -114,5 +114,4 @@
             CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN
          ENDIF
-         !
          !                         ! fill sf_tsd with sn_tem & sn_sal and control print
          slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal
@@ -150,9 +149,9 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE dta_tsd  ***
-      !!                    
+      !!
       !! ** Purpose :   provides T and S data at kt
-      !! 
+      !!
       !! ** Method  : - call fldread routine
-      !!              - ORCA_R2: add some hand made alteration to read data  
+      !!              - ORCA_R2: add some hand made alteration to read data
       !!              - 'key_orca_lev10' interpolates on 10 times more levels
       !!              - s- or mixed z-s coordinate: vertical interpolation on model mesh
@@ -162,6 +161,6 @@
       !! ** Action  :   ptsd   T-S data on medl mesh and interpolated at time-step kt
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt     ! ocean time-step
-      CHARACTER(LEN=3)                     , INTENT(in   ) ::   cddta  ! dmp or ini
+      INTEGER                          , INTENT(in   ) ::   kt     ! ocean time-step
+      CHARACTER(LEN=3)                 , INTENT(in   ) ::   cddta  ! dmp or ini
       REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data
       !
@@ -224,5 +223,5 @@
                      IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
                         zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
+                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi
                         zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
                      ENDIF
@@ -237,7 +236,7 @@
             ptsd(ji,jj,jpk,jp_sal) = 0._wp
          END_2D
-         ! 
+         !
       ELSE                                !==   z- or zps- coordinate   ==!
-         !                             
+         !
          DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
             ptsd(ji,jj,jk,jp_tem) = ptsd(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)    ! Mask
@@ -247,5 +246,5 @@
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
             DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-               ik = mbkt(ji,jj) 
+               ik = mbkt(ji,jj)
                IF( ik > 1 ) THEN
                   zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
@@ -255,5 +254,5 @@
                ik = mikt(ji,jj)
                IF( ik > 1 ) THEN
-                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
+                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
                   ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
                   ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 14923)
@@ -31,5 +31,5 @@
    !!   bn2           : compute the Brunt-Vaisala frequency
    !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature
-   !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio 
+   !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio
    !!   eos_rab_3d    : compute in situ thermal/haline expansion ratio
    !!   eos_rab_2d    : compute in situ thermal/haline expansion ratio for 2d fields
@@ -46,5 +46,5 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
    USE prtctl         ! Print control
    USE lbclnk         ! ocean lateral boundary conditions
@@ -63,5 +63,5 @@
    END INTERFACE
    !
-   INTERFACE eos_fzp 
+   INTERFACE eos_fzp
       MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
    END INTERFACE
@@ -91,15 +91,15 @@
 
    !                               !!!  simplified eos coefficients (default value: Vallis 2006)
-   REAL(wp), PUBLIC ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff. 
-   REAL(wp), PUBLIC ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff. 
-   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2        
-   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2        
-   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T  
-   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S  
-   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt  
-   
+   REAL(wp), PUBLIC ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff.
+   REAL(wp), PUBLIC ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff.
+   REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2
+   REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2
+   REAL(wp) ::   rn_mu1     = 1.4970e-4_wp     ! thermobaric coeff. in T
+   REAL(wp) ::   rn_mu2     = 1.1090e-5_wp     ! thermobaric coeff. in S
+   REAL(wp) ::   rn_nu      = 2.4341e-3_wp     ! cabbeling coeff. in theta*salt
+
    ! TEOS10/EOS80 parameters
    REAL(wp) ::   r1_S0, r1_T0, r1_Z0, rdeltaS
-   
+
    ! EOS parameters
    REAL(wp) ::   EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
@@ -119,6 +119,6 @@
    REAL(wp) ::   EOS022
    REAL(wp) ::   EOS003 , EOS103
-   REAL(wp) ::   EOS013 
-   
+   REAL(wp) ::   EOS013
+
    ! ALPHA parameters
    REAL(wp) ::   ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
@@ -135,5 +135,5 @@
    REAL(wp) ::   ALP012
    REAL(wp) ::   ALP003
-   
+
    ! BETA parameters
    REAL(wp) ::   BET000 , BET100 , BET200 , BET300 , BET400 , BET500
@@ -162,5 +162,5 @@
    REAL(wp) ::   PEN002 , PEN102
    REAL(wp) ::   PEN012
-   
+
    ! ALPHA_PEN parameters
    REAL(wp) ::   APE000 , APE100 , APE200 , APE300
@@ -241,5 +241,5 @@
       INTEGER                                 , INTENT(in   ) ::   ktts, ktrd, ktdep
       REAL(wp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
-      !                                                               ! 2 : salinity               [psu]
+      !                                                                  ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(A2D_T(ktrd) ,JPK     ), INTENT(  out) ::   prd   ! in situ density            [-]
       REAL(wp), DIMENSION(A2D_T(ktdep),JPK     ), INTENT(in   ) ::   pdep  ! depth                      [m]
@@ -301,5 +301,5 @@
                &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
                &  - rn_nu * zt * zs
-               !                                 
+               !
             prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
          END_3D
@@ -354,5 +354,5 @@
       INTEGER                                  , INTENT(in   ) ::   ktts, ktrd, ktrhop, ktdep
       REAL(wp), DIMENSION(A2D_T(ktts)  ,JPK,JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
-      !                                                                ! 2 : salinity               [psu]
+      !                                                                    ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(A2D_T(ktrd)  ,JPK     ), INTENT(  out) ::   prd    ! in situ density            [-]
       REAL(wp), DIMENSION(A2D_T(ktrhop),JPK     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
@@ -467,5 +467,5 @@
             END_3D
          ENDIF
-         
+
       CASE( np_seos )                !==  simplified EOS  ==!
          !
@@ -503,5 +503,6 @@
       END SELECT
       !
-      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', &
+         &                                  tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
       !
       IF( ln_timing )   CALL timing_stop('eos-pot')
@@ -534,5 +535,5 @@
       INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktrd
       REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
-      !                                                           ! 2 : salinity               [psu]
+      !                                                             ! 2 : salinity               [psu]
       REAL(wp), DIMENSION(A2D_T(ktdep)    ), INTENT(in   ) ::   pdep  ! depth                      [m]
       REAL(wp), DIMENSION(A2D_T(ktrd)     ), INTENT(  out) ::   prd   ! in situ density
@@ -666,21 +667,21 @@
          !
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-               !
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-               ztm = tmask(ji,jj,1)                                         ! tmask
-               !
-               zn0 = (((((EOS060*zt   &
-                  &   + EOS150*zs+EOS050)*zt   &
-                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                  !
-               !
-               prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
-               !
-            END_2D
+            !
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,1)                                         ! tmask
+            !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            !
+            prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
+            !
+         END_2D
 
       CASE( np_seos )                !==  simplified EOS  ==!
@@ -713,4 +714,5 @@
          !
       END SELECT
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prhop, clinfo1=' pot: ', kdim=1 )
       !
       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prhop, clinfo1=' eos-pot: ' )
@@ -741,5 +743,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
+      INTEGER                                , INTENT(in   ) ::   Kmm   ! time level index
       INTEGER                                , INTENT(in   ) ::   ktts, ktab
       REAL(wp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in   ) ::   pts   ! pot. temperature & salinity
@@ -873,5 +875,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
+      INTEGER                            , INTENT(in   ) ::   Kmm   ! time level index
       INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktab
       REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! pot. temperature & salinity
@@ -1113,5 +1115,5 @@
       !!                  ***  ROUTINE bn2  ***
       !!
-      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the 
+      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the
       !!                time-step of the input arguments
       !!
@@ -1120,10 +1122,10 @@
       !!      N.B. N^2 is set one for all to zero at jk=1 in istate module.
       !!
-      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point 
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
+      !! ** Action  :   pn2 : square of the brunt-vaisala frequency at w-point
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                                , INTENT(in   ) ::  Kmm   ! time level index
       INTEGER                                , INTENT(in   ) ::  ktab, ktn2
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
+      REAL(wp), DIMENSION(jpi,jpj,  jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
       REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
       REAL(wp), DIMENSION(A2D_T(ktn2),JPK     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2]
@@ -1137,7 +1139,7 @@
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpkm1 )      ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90
          zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
-            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
-            !
-         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
+            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )
+            !
+         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw
          zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
          !
@@ -1211,5 +1213,5 @@
 
 
-   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep )
+   SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
       !!
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
@@ -1234,6 +1236,6 @@
       !!----------------------------------------------------------------------
       INTEGER                       , INTENT(in   )           ::   kttf
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
+      REAL(wp), DIMENSION(jpi,jpj)  , INTENT(in   )           ::   psal   ! salinity   [psu]
+      REAL(wp), DIMENSION(jpi,jpj)  , INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m]
       REAL(wp), DIMENSION(A2D_T(kttf)), INTENT(out  )           ::   ptf    ! freezing temperature [Celsius]
       !
@@ -1267,5 +1269,5 @@
          CALL ctl_stop( 'eos_fzp_2d:', ctmp1 )
          !
-      END SELECT      
+      END SELECT
       !
   END SUBROUTINE eos_fzp_2d_t
@@ -1324,8 +1326,8 @@
       !! ** Purpose :   Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
       !!
-      !! ** Method  :   PE is defined analytically as the vertical 
+      !! ** Method  :   PE is defined analytically as the vertical
       !!                   primitive of EOS times -g integrated between 0 and z>0.
       !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd
-      !!                                                      = 1/z * /int_0^z rd dz - rd 
+      !!                                                      = 1/z * /int_0^z rd dz - rd
       !!                                where rd is the density anomaly (see eos_rhd function)
       !!                ab_pe are partial derivatives of PE anomaly with respect to T and S:
@@ -1391,5 +1393,5 @@
                !
             zn  = ( zn2 * zh + zn1 ) * zh + zn0
-            !                              
+            !
             pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
             !
@@ -1406,5 +1408,5 @@
                !
             zn  = ( zn2 * zh + zn1 ) * zh + zn0
-            !                              
+            !
             pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
             !
@@ -1504,5 +1506,5 @@
          IF(lwp) WRITE(numout,*) '   ==>>>   use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
          !
-         l_useCT = .TRUE.                          ! model temperature is Conservative temperature 
+         l_useCT = .TRUE.                          ! model temperature is Conservative temperature
          !
          rdeltaS = 32._wp
@@ -1885,5 +1887,5 @@
 
          r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
-         
+
          IF(lwp) THEN
             WRITE(numout,*)
@@ -1923,8 +1925,8 @@
       END SELECT
       !
-      rho0_rcp    = rho0 * rcp 
+      rho0_rcp    = rho0 * rcp
       r1_rho0     = 1._wp / rho0
       r1_rcp      = 1._wp / rcp
-      r1_rho0_rcp = 1._wp / rho0_rcp 
+      r1_rho0_rcp = 1._wp / rho0_rcp
       !
       IF(lwp) THEN
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isf_oce.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 14923)
@@ -1,6 +1,6 @@
 MODULE isf_oce
    !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  compute iceshelf melt and heat flux
+   !!                       ***  MODULE  isf_oce  ***
+   !! Ice shelves :  ice shelves variables defined in memory
    !!======================================================================
    !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
@@ -48,10 +48,11 @@
    !
    ! 0.3 -------- ice shelf cavity parametrised namelist parameter -------------
-   LOGICAL           , PUBLIC :: ln_isfpar_mlt   !: logical for the computation of melt inside the cavity
-   CHARACTER(LEN=256), PUBLIC :: cn_isfpar_mlt   !: melt formulation (cavity/param)
-   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_fwf   !: information about the isf melting file to be read
-   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmax  !: information about the grounding line depth file to be read
-   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmin  !: information about the calving   line depth file to be read
-   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_Leff  !: information about the effective length     file to be read
+   LOGICAL           , PUBLIC :: ln_isfpar_mlt      !: logical for the computation of melt inside the cavity
+   REAL(wp)          , PUBLIC :: rn_isfpar_bg03_gt0 !: temperature exchange coeficient [m/s]
+   CHARACTER(LEN=256), PUBLIC :: cn_isfpar_mlt      !: melt formulation (cavity/param)
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_fwf      !: information about the isf melting file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmax     !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmin     !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_Leff     !: information about the effective length     file to be read
    !
    ! 0.4 -------- coupling namelist parameter -------------
@@ -147,4 +148,5 @@
    END SUBROUTINE isf_alloc_par
 
+   
    SUBROUTINE isf_alloc_cav()
       !!---------------------------------------------------------------------
@@ -174,4 +176,5 @@
    END SUBROUTINE isf_alloc_cav
 
+   
    SUBROUTINE isf_alloc_cpl()
       !!---------------------------------------------------------------------
@@ -185,14 +188,14 @@
       ierr = 0
       !
-      ALLOCATE( risfcpl_ssh(jpi,jpj), risfcpl_tsc(jpi,jpj,jpk,jpts), risfcpl_vol(jpi,jpj,jpk), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      risfcpl_tsc(:,:,:,:) = 0.0 ; risfcpl_vol(:,:,:) = 0.0 ; risfcpl_ssh(:,:) = 0.0
-
-      IF ( ln_isfcpl_cons) THEN
-         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) ,risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )
+      ALLOCATE( risfcpl_ssh(jpi,jpj) , risfcpl_tsc(jpi,jpj,jpk,jpts) , risfcpl_vol(jpi,jpj,jpk) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      risfcpl_tsc(:,:,:,:) = 0._wp ; risfcpl_vol(:,:,:) = 0._wp ; risfcpl_ssh(:,:) = 0._wp
+
+      IF ( ln_isfcpl_cons ) THEN
+         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) , risfcpl_cons_ssh(jpi,jpj) , STAT=ialloc )
          ierr = ierr + ialloc
          !
-         risfcpl_cons_tsc(:,:,:,:) = 0.0 ; risfcpl_cons_vol(:,:,:) = 0.0 ; risfcpl_cons_ssh(:,:) = 0.0
+         risfcpl_cons_tsc(:,:,:,:) = 0._wp ; risfcpl_cons_vol(:,:,:) = 0._wp ; risfcpl_cons_ssh(:,:) = 0._wp
          !
       END IF
@@ -203,4 +206,5 @@
    END SUBROUTINE isf_alloc_cpl
 
+   
    SUBROUTINE isf_dealloc_cpl()
       !!---------------------------------------------------------------------
@@ -214,5 +218,5 @@
       ierr = 0
       !
-      DEALLOCATE( risfcpl_ssh, risfcpl_tsc, risfcpl_vol, STAT=ialloc )
+      DEALLOCATE( risfcpl_ssh , risfcpl_tsc , risfcpl_vol , STAT=ialloc )
       ierr = ierr + ialloc
       !
@@ -222,4 +226,5 @@
    END SUBROUTINE isf_dealloc_cpl
 
+   
    SUBROUTINE isf_alloc()
       !!---------------------------------------------------------------------
@@ -234,19 +239,19 @@
       ierr = 0       ! set to zero if no array to be allocated
       !
-      ALLOCATE(fwfisf_par(jpi,jpj)  , fwfisf_par_b(jpi,jpj), &
-         &     fwfisf_cav(jpi,jpj)  , fwfisf_cav_b(jpi,jpj), &
-         &     fwfisf_oasis(jpi,jpj),            STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
-      ierr = ierr + ialloc
-      !
-      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
-      ierr = ierr + ialloc
-      !
-      ALLOCATE( mskisf_cav(jpi,jpj), STAT=ialloc)
+      ALLOCATE( fwfisf_par  (jpi,jpj) , fwfisf_par_b(jpi,jpj) ,     &
+         &      fwfisf_cav  (jpi,jpj) , fwfisf_cav_b(jpi,jpj) ,     &
+         &      fwfisf_oasis(jpi,jpj)                         , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_cav(jpi,jpj) , STAT=ialloc )
       ierr = ierr + ialloc
       !
@@ -255,13 +260,13 @@
       !
       ! initalisation of fwf and tsc array to 0
-      risfload(:,:)       = 0.0_wp
-      fwfisf_oasis(:,:)   = 0.0_wp
-      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
-      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
-      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
-      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
-      !
-
+      risfload    (:,:)   = 0._wp
+      fwfisf_oasis(:,:)   = 0._wp
+      fwfisf_par  (:,:)   = 0._wp   ;   fwfisf_par_b  (:,:)   = 0._wp
+      fwfisf_cav  (:,:)   = 0._wp   ;   fwfisf_cav_b  (:,:)   = 0._wp
+      risf_cav_tsc(:,:,:) = 0._wp   ;   risf_cav_tsc_b(:,:,:) = 0._wp
+      risf_par_tsc(:,:,:) = 0._wp   ;   risf_par_tsc_b(:,:,:) = 0._wp
+      !
    END SUBROUTINE isf_alloc
-
+   
+   !!======================================================================
 END MODULE isf_oce
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 14923)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                       ***  MODULE  isfstp  ***
-   !! Surface module :  compute iceshelf load, melt and heat flux
+   !! Ice Shelves :  compute iceshelf load, melt and heat flux
    !!======================================================================
    !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
@@ -42,8 +42,7 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
  
-  SUBROUTINE isf_stp( kt, Kmm )
+   SUBROUTINE isf_stp( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isf_stp  ***
@@ -58,11 +57,12 @@
       !!              - compute fluxes
       !!              - write restart variables
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
-      !!----------------------------------------------------------------------
-      INTEGER :: jk                               ! loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t    ! e3t 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! ocean time step
+      INTEGER, INTENT(in) ::   Kmm   ! ocean time level index
+      !
+      INTEGER :: jk                              ! loop index
+#if defined key_qco
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t   ! 3D workspace
+#endif
       !!---------------------------------------------------------------------
       !
@@ -83,11 +83,15 @@
          ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
          rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
+#if defined key_qco
          DO jk = 1, jpk
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO 
-         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         CALL isf_tbl_lvl( ht(:,:), ze3t           , misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
+#else
+         CALL isf_tbl_lvl( ht(:,:),  e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav )
+#endif
          !
          ! 1.3: compute ice shelf melt
-         CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav)
+         CALL isf_cav( kt, Kmm, risf_cav_tsc, fwfisf_cav )
          !
       END IF
@@ -108,11 +112,15 @@
          ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
          rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
+#if defined key_qco
          DO jk = 1, jpk
             ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
          END DO
-         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         CALL isf_tbl_lvl( ht(:,:), ze3t           , misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+#else
+         CALL isf_tbl_lvl( ht(:,:),  e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
+#endif
          !
          ! 2.3: compute ice shelf melt
-         CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par)
+         CALL isf_par( kt, Kmm, risf_par_tsc, fwfisf_par )
          !
       END IF
@@ -128,5 +136,6 @@
    END SUBROUTINE isf_stp
 
-   SUBROUTINE isf_init(Kbb, Kmm, Kaa)
+   
+   SUBROUTINE isf_init( Kbb, Kmm, Kaa )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isfstp_init  ***
@@ -142,34 +151,25 @@
       !!              - call cav/param/isfcpl init routine
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa      ! ocean time level indices
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time level indices
+      !!----------------------------------------------------------------------
       !
       ! constrain: l_isfoasis need to be known
       !
-      ! Read namelist
-      CALL isf_nam()
-      !
-      ! Allocate public array
-      CALL isf_alloc()
-      !
-      ! check option compatibility
-      CALL isf_ctl()
-      !
-      ! compute ice shelf load
-      IF ( ln_isfcav ) CALL isf_load( Kmm, risfload )
+      CALL isf_nam()                                              ! Read namelist
+      !
+      CALL isf_alloc()                                            ! Allocate public array
+      !
+      CALL isf_ctl()                                              ! check option compatibility
+      !
+      IF( ln_isfcav ) CALL isf_load( Kmm, risfload )              ! compute ice shelf load
       !
       ! terminate routine now if no ice shelf melt formulation specify
-      IF ( ln_isf ) THEN
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation melt in the cavity
-         IF ( ln_isfcav_mlt ) CALL isf_cav_init()
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation parametrised melt
-         IF ( ln_isfpar_mlt ) CALL isf_par_init()
-         !
-         !---------------------------------------------------------------------------------------------------------------------
-         ! initialisation ice sheet coupling
-         IF( ln_isfcpl ) CALL isfcpl_init(Kbb, Kmm, Kaa)
+      IF( ln_isf ) THEN
+         !
+         IF( ln_isfcav_mlt )   CALL isf_cav_init()                ! initialisation melt in the cavity
+         !
+         IF( ln_isfpar_mlt )   CALL isf_par_init()                ! initialisation parametrised melt
+         !
+         IF( ln_isfcpl     )   CALL isfcpl_init( Kbb, Kmm, Kaa )  ! initialisation ice sheet coupling
          !
       END IF
@@ -177,4 +177,5 @@
   END SUBROUTINE isf_init
 
+  
   SUBROUTINE isf_ctl()
       !!---------------------------------------------------------------------
@@ -283,5 +284,6 @@
       END IF
    END SUBROUTINE isf_ctl
-   !
+
+   
    SUBROUTINE isf_nam
       !!---------------------------------------------------------------------
@@ -299,5 +301,6 @@
          &             sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff,                           &
          &             ln_isfcpl     , nn_drown      , ln_isfcpl_cons, ln_isfdebug, rn_vtide,    &
-         &             cn_isfload    , rn_isfload_T  , rn_isfload_S  , cn_isfdir
+         &             cn_isfload    , rn_isfload_T  , rn_isfload_S  , cn_isfdir  ,              &
+         &             rn_isfpar_bg03_gt0
       !!----------------------------------------------------------------------
       !
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/istate.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/istate.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/istate.F90	(revision 14923)
@@ -34,7 +34,8 @@
    USE lib_mpp         ! MPP library
    USE restart         ! restart
+
 #if defined key_agrif
+   USE agrif_oce       ! initial state interpolation
    USE agrif_oce_interp
-   USE agrif_oce
 #endif   
 
@@ -42,5 +43,5 @@
    PRIVATE
 
-   PUBLIC   istate_init   ! routine called by step.F90
+   PUBLIC   istate_init   ! routine called by nemogcm.F90
 
    !! * Substitutions
@@ -63,5 +64,5 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table for qco substitute
 !!gm see comment further down
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
@@ -73,15 +74,12 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
 
-!!gm  Why not include in the first call of dta_tsd ?  
-!!gm  probably associated with the use of internal damping...
        CALL dta_tsd_init        ! Initialisation of T & S input data
-!!gm to be moved in usrdef of C1D case
+
 !      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
-!!gm
 
-      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
-      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
-      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
+      rhd  (:,:,:      ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
+      rn2b (:,:,:      ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
+      ts   (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
+      rab_b(:,:,:,:    ) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
       uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
@@ -94,66 +92,51 @@
          ln_1st_euler = .true.                ! Set time-step indicator at nit000 (euler forward)
          CALL day_init 
-         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent
+         CALL agrif_istate_oce( Kbb, Kmm, Kaa )   ! Interp from parent
          !
-         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
-         ssh (:,:,Kmm)     = ssh(:,:,Kbb)
-         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
-         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)
+         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
+         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
       ELSE
 #endif
-      IF( ln_rstart ) THEN                    ! Restart from a file
-         !                                    ! -------------------
-         CALL rst_read( Kbb, Kmm )            ! Read the restart file
-         CALL day_init                        ! model calendar (using both namelist and restart infos)
-         !
-      ELSE                                    ! Start from rest
-         !                                    ! ---------------
-         numror = 0                           ! define numror = 0 -> no restart file to read
-         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
-         CALL day_init                        ! model calendar (using both namelist and restart infos)
-         !                                    ! Initialization of ocean to zero
-         !
-         IF( ln_tsd_init ) THEN               
-            CALL dta_tsd( nit000, 'ini', ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
+         IF( ln_rstart ) THEN                    ! Restart from a file
+            !                                    ! -------------------
+            CALL rst_read( Kbb, Kmm )            ! Read the restart file
+            CALL day_init                        ! model calendar (using both namelist and restart infos)
             !
-            uu (:,:,:,Kbb) = 0._wp
-            vv (:,:,:,Kbb) = 0._wp
+         ELSE                                    ! Start from rest
+            !                                    ! ---------------
+            numror = 0                           ! define numror = 0 -> no restart file to read
+            l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
+            CALL day_init                        ! model calendar (using both namelist and restart infos)
+            !                                    ! Initialization of ocean to zero
             !
-         ELSE                                 ! user defined initial T and S
-            DO jk = 1, jpk
-               zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
-            END DO
-            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb)  )         
-         ENDIF
-         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
-         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
-         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
-
-!!gm POTENTIAL BUG :
-!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
-!!             as well as gdept_ and gdepw_....   !!!!! 
-!!      ===>>>>   probably a call to domvvl initialisation here....
-
+            IF( ln_tsd_init ) THEN               
+               CALL dta_tsd( nit000, 'ini', ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
+               !
+               uu  (:,:,:,Kbb) = 0._wp               ! set the ocean at rest
+               vv  (:,:,:,Kbb) = 0._wp  
+               !
+            ELSE                                 ! user defined initial T and S
+               DO jk = 1, jpk
+                  zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
+               END DO
+               CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )         
+            ENDIF
+            ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+            uu    (:,:,:,Kmm) = uu   (:,:,:,Kbb)
+            vv    (:,:,:,Kmm) = vv   (:,:,:,Kbb)
 
          !
-!!gm to be moved in usrdef of C1D case
-!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-!            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
-!            CALL dta_uvd( nit000, zuvd )
-!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
-!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
-!            DEALLOCATE( zuvd )
-!         ENDIF
+!!gm ==>>>  to be moved in usrdef_istate of C1D case 
+         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
+            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
+            CALL dta_uvd( nit000, Kbb, zuvd )
+            uu(:,:,:,Kbb) = zuvd(:,:,:,1)  ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
+            vv(:,:,:,Kbb) = zuvd(:,:,:,2)  ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
+            DEALLOCATE( zuvd )
+         ENDIF
          !
-!!gm This is to be changed !!!!
-!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
-!         IF( .NOT.ln_linssh ) THEN
-!            DO jk = 1, jpk
-!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
-!            END DO
-!         ENDIF
-!!gm 
          ! 
-      ENDIF 
+         ENDIF 
 #if defined key_agrif
       ENDIF
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 14923)
@@ -24,4 +24,5 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! IOM
    USE lib_mpp        ! distribued memory computing library
    USE timing         ! Timing
@@ -34,7 +35,7 @@
    PUBLIC   sbc_fwb    ! routine called by step
 
-   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget
-   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year.
-   REAL(wp) ::   fwfold    ! fwfold to be suppressed
+   REAL(wp) ::   rn_fwb0   ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only)
+   REAL(wp) ::   a_fwb     ! annual domain averaged freshwater budget from the
+                           ! previous year
    REAL(wp) ::   area      ! global mean ocean surface (interior domain)
 
@@ -65,5 +66,5 @@
       INTEGER, INTENT( in ) ::   Kmm      ! ocean time level index
       !
-      INTEGER  ::   inum, ikty, iyear     ! local integers
+      INTEGER  ::   ios, inum, ikty       ! local integers
       REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp                ! local scalars
       REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread, zcoef          !   -      -
@@ -72,7 +73,14 @@
       REAL(wp)   ,DIMENSION(1) ::   z_fwfprv  
       COMPLEX(dp),DIMENSION(1) ::   y_fwfnow  
+      !
+      NAMELIST/namsbc_fwb/rn_fwb0
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 ) THEN
+         READ( numnam_ref, namsbc_fwb, IOSTAT = ios, ERR = 901 )
+901      IF( ios /= 0 ) CALL ctl_nam( ios, 'namsbc_fwb in reference namelist'     )
+         READ( numnam_cfg, namsbc_fwb, IOSTAT = ios, ERR = 902 )
+902      IF( ios >  0 ) CALL ctl_nam( ios, 'namsbc_fwb in configuration namelist' )
+         IF(lwm) WRITE( numond, namsbc_fwb )
          IF(lwp) THEN
             WRITE(numout,*)
@@ -80,7 +88,12 @@
             WRITE(numout,*) '~~~~~~~'
             IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero'
-            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget'
+            IF( kn_fwb == 4 )   WRITE(numout,*) '          instantaneously set to zero with heat and salt flux correction (ISOMIP+)'
             IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area'
-            IF( kn_fwb == 4 )   WRITE(numout,*) '          instantaneously set to zero with heat and salt flux correction (ISOMIP+)'
+            IF( kn_fwb == 2 ) THEN
+               WRITE(numout,*) '          adjusted from previous year budget'
+               WRITE(numout,*)
+               WRITE(numout,*) '   Namelist namsbc_fwb'
+               WRITE(numout,*) '      Initial freshwater adjustment flux [kg/m2/s] = ', rn_fwb0
+            END IF
          ENDIF
          !
@@ -111,4 +124,7 @@
             emp(:,:) = emp(:,:) - z_fwfprv(1)        * tmask(:,:,1)
             qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', zcoef * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1)        * tmask(:,:,1) )
          ENDIF
          !
@@ -131,47 +147,48 @@
             qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
             sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
-            !qns(:,:) = qns(:,:) + zcoef * ( -1.9 ) * tmask(:,:,1) ! (Eq. 35 AD2015) ! could be sst_m if we don't want any bouyancy fluxes
-            !sfx(:,:) = sfx(:,:) + z_fwf * ( 33.8 ) * tmask(:,:,1) ! (Eq. 36 AD2015) ! could be sss_m if we don't want any bouyancy fluxes
-            !qns(:,:) = qns(:,:) + zcoef * ( -1.0 ) * tmask(:,:,1) ! use for ISOMIP+ coupling sanity check (keep ssh cst while playing with cpl conservation option)
-            !sfx(:,:) = sfx(:,:) + z_fwf * ( 34.2 ) * tmask(:,:,1) ! use for ISOMIP+ coupling sanity check (keep ssh cst while playing with cpl conservation option)
-         ENDIF
-         !
-      CASE ( 2 )                             !==  fwf budget adjusted from the previous year  ==!
-         !
-         IF( kt == nit000 ) THEN                      ! initialisation
-            !                                         ! Read the corrective factor on precipitations (fwfold)
-            CALL ctl_opn( inum, 'EMPave_old.dat', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-            READ ( inum, "(24X,I8,2ES24.16)" ) iyear, a_fwb_b, a_fwb
-            CLOSE( inum )
-            fwfold = a_fwb                            ! current year freshwater budget correction
-            !                                         ! estimate from the previous year budget
+         ENDIF
+         !
+      CASE ( 2 )                             !==  fw adjustment based on fw budget at the end of the previous year  ==!
+         !
+         IF( kt == nit000 ) THEN                                                                    ! initialisation
+            !                                                                                       ! set the fw adjustment (a_fwb)
+            IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb',   ldstop = .FALSE. ) > 0 ) THEN        !    as read from restart file
+               IF(lwp) WRITE(numout,*) 'sbc_fwb : reading FW-budget adjustment from restart file'
+               CALL iom_get( numror, 'a_fwb',   a_fwb )
+            ELSE                                                                                    !    as specified in namelist
+               a_fwb = rn_fwb0
+            END IF
+            !
             IF(lwp)WRITE(numout,*)
-            IF(lwp)WRITE(numout,*)'sbc_fwb : year = ',iyear  , ' freshwater budget correction = ', fwfold
-            IF(lwp)WRITE(numout,*)'          year = ',iyear-1, ' freshwater budget read       = ', a_fwb
-            IF(lwp)WRITE(numout,*)'          year = ',iyear-2, ' freshwater budget read       = ', a_fwb_b
+            IF(lwp)WRITE(numout,*)'sbc_fwb : initial freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
+            !
          ENDIF   
-         !                                         ! Update fwfold if new year start
+         !                                         ! Update a_fwb if new year start
          ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!!
          IF( MOD( kt, ikty ) == 0 ) THEN
-            a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
+                                                      ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
             a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
-            fwfold =  a_fwb                           ! current year freshwater budget correction
-            !                                         ! estimate from the previous year budget
          ENDIF
          ! 
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
-            zcoef = fwfold * rcp
-            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1)
+            zcoef = a_fwb * rcp
+            emp(:,:) = emp(:,:) + a_fwb              * tmask(:,:,1)
             qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
-         ENDIF
-         !
-         IF( kt == nitend .AND. lwm ) THEN            ! save fwfold value in a file (only one required)
-            CALL ctl_opn( inum, 'EMPave.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
-            WRITE( inum, "(24X,I8,2ES24.16)" ) nyear, a_fwb_b, a_fwb
-            CLOSE( inum )
-         ENDIF
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zcoef * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -a_fwb              * tmask(:,:,1) )
+         ENDIF
+         ! Output restart information
+         IF( lrst_oce ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
+            IF(lwp) WRITE(numout,*) '~~~~'
+            CALL iom_rstput( kt, nitrst, numrow, 'a_fwb',   a_fwb )
+         END IF
+         !
+         IF( kt == nitend .AND. lwp ) WRITE(numout,*) 'sbc_fwb : final freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
          !
       CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
@@ -211,4 +228,7 @@
             qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction
             erp(:,:) = erp(:,:) + zerp_cor(:,:)
+            ! outputs
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zerp_cor(:,:) * rcp * sst_m(:,:) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -zerp_cor(:,:) )
             !
             IF( lwp ) THEN                   ! control print
Index: NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/tradmp.F90
===================================================================
--- NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 14922)
+++ NEMO/branches/2021/ticket2669_isf_fluxes/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 14923)
@@ -11,5 +11,5 @@
    !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules
    !!            3.2  ! 2009-08  (G. Madec, C. Talandier)  DOCTOR norm for namelist parameter
-   !!            3.3  ! 2010-06  (C. Ethe, G. Madec) merge TRA-TRC 
+   !!            3.3  ! 2010-06  (C. Ethe, G. Madec) merge TRA-TRC
    !!            3.4  ! 2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression of CPP keys
    !!            3.6  ! 2015-06  (T. Graham)  read restoring coefficient in a file
@@ -26,5 +26,5 @@
    USE c1d            ! 1D vertical configuration
    USE trd_oce        ! trends: ocean variables
-   USE trdtra         ! trends manager: tracers 
+   USE trdtra         ! trends manager: tracers
    USE zdf_oce        ! ocean: vertical physics
    USE phycst         ! physical constants
@@ -75,10 +75,10 @@
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_dmp  ***
-      !!                  
+      !!
       !! ** Purpose :   Compute the tracer trend due to a newtonian damping
       !!      of the tracer field towards given data field and add it to the
       !!      general tracer trends.
       !!
-      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed 
+      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
       !!      and add to the general tracer trends:
       !!                     ta = ta + resto * (t_dta - tb)
@@ -101,7 +101,7 @@
       IF( ln_timing )   CALL timing_start('tra_dmp')
       !
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 
+      IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN   !* Save ta and sa trends
+         ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )
+         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs)
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -140,9 +140,17 @@
       END SELECT
       !
+      ! outputs (clem trunk)
+      IF( iom_use('hflx_dmp_cea') )       &
+         &   CALL iom_put('hflx_dmp_cea', &
+         &   SUM( ( pts(:,:,:,jp_tem,Krhs) - ztrdts(:,:,:,jp_tem) ) * e3t(:,:,:,Kmm), dim=3 ) * rcp * rho0 ) ! W/m2
+      IF( iom_use('sflx_dmp_cea') )       &
+         &   CALL iom_put('sflx_dmp_cea', &
+         &   SUM( ( pts(:,:,:,jp_sal,Krhs) - ztrdts(:,:,:,jp_sal) ) * e3t(:,:,:,Kmm), dim=3 ) * rho0 )       ! g/m2/s
+      !
       IF( l_trdtra )   THEN       ! trend diagnostic
          ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
-         DEALLOCATE( ztrdts ) 
+         DEALLOCATE( ztrdts )
       ENDIF
       !                           ! Control print
@@ -158,10 +166,10 @@
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_dmp_init  ***
-      !! 
-      !! ** Purpose :   Initialization for the newtonian damping 
+      !!
+      !! ** Purpose :   Initialization for the newtonian damping
       !!
       !! ** Method  :   read the namtra_dmp namelist and check the parameters
       !!----------------------------------------------------------------------
-      INTEGER ::   ios, imask   ! local integers 
+      INTEGER ::   ios, imask   ! local integers
       !
       NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto
