Index: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcflx.F90
===================================================================
--- NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcflx.F90	(revision 13483)
+++ NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcflx.F90	(revision 13484)
@@ -29,5 +29,4 @@
    PUBLIC sbc_flx       ! routine called by step.F90
 
-   INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read 
    INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file
    INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file
@@ -35,4 +34,7 @@
    INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file
    INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
+ !!INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux flux
+   INTEGER , PARAMETER ::   jpfld   = 5 !! 6 ! maximum number of files to read 
+
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
 
@@ -59,4 +61,5 @@
       !!                   net downward radiative flux            qsr   (watt/m2)
       !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s)
+      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/m2/s)
       !!
       !!      CAUTION :  - never mask the surface stress fields
@@ -71,5 +74,5 @@
       !!              - emp         upward mass flux (evap. - precip.)
       !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero
-      !!                            if ice is present
+      !!                            if ice
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
@@ -85,6 +88,6 @@
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
-      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read
-      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp
+      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx ! informations about the fields to be read
+      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx
       !!---------------------------------------------------------------------
       !
@@ -107,5 +110,5 @@
          slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau
          slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr 
-         slf_i(jp_emp ) = sn_emp
+         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx
          !
          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
@@ -120,6 +123,4 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' )
          !
-         sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present)
-         !
       ENDIF
 
@@ -128,35 +129,24 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
 
-         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
-         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
+         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)  ! modify now Qsr to include the diurnal cycle
+         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
          ENDIF
          DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
             DO ji = 1, jpi
-               utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
-               vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
-               qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
-               emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
+               utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1)
+               vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1)
+               qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
+               emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1)
+               !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1) 
             END DO
          END DO
          !                                                        ! add to qns the heat due to e-p
-         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
+         !clem: I do not think it is needed
+         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
          !
-         qns(:,:) = qns(:,:) * tmask(:,:,1)
-         emp(:,:) = emp(:,:) * tmask(:,:,1)
+         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x) 
+         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, &
+            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  )
          !
-         !                                                        ! module of wind stress and wind speed at T-point
-         zcoef = 1. / ( zrhoa * zcdrag )
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vect. opt.
-               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-               taum(ji,jj) = zmod
-               wndm(ji,jj) = SQRT( zmod * zcoef )
-            END DO
-         END DO
-         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
-         CALL lbc_lnk( 'sbcflx', taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( 'sbcflx', wndm(:,:), 'T', 1. )
-
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
             WRITE(numout,*) 
@@ -172,4 +162,19 @@
          !
       ENDIF
+      !                                                           ! module of wind stress and wind speed at T-point
+      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
+      zcoef = 1. / ( zrhoa * zcdrag )
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vect. opt.
+            ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) )
+            zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) ) 
+            zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1)
+            taum(ji,jj) = zmod
+            wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used?
+         END DO
+      END DO
+      !
+      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp )
+      !
       !
    END SUBROUTINE sbc_flx
