New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6043 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 – NEMO

Ignore:
Timestamp:
2015-12-14T10:27:28+01:00 (8 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5600 r6043  
    1111   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator 
    1212   !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML 
     13   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  add limiter on triad slopes 
    1314   !!---------------------------------------------------------------------- 
    14 #if   defined key_ldfslp   ||   defined key_esopa 
     15 
    1516   !!---------------------------------------------------------------------- 
    16    !!   'key_ldfslp'                      Rotation of lateral mixing tensor 
    17    !!---------------------------------------------------------------------- 
    18    !!   ldf_slp_grif  : calculates the triads of isoneutral slopes (Griffies operator) 
    1917   !!   ldf_slp       : calculates the slopes of neutral surface   (Madec operator) 
     18   !!   ldf_slp_triad : calculates the triads of isoneutral slopes (Griffies operator) 
    2019   !!   ldf_slp_mxl   : calculates the slopes at the base of the mixed layer (Madec operator) 
    2120   !!   ldf_slp_init  : initialization of the slopes computation 
     
    2322   USE oce            ! ocean dynamics and tracers 
    2423   USE dom_oce        ! ocean space and time domain 
    25    USE ldftra_oce     ! lateral diffusion: traceur 
    26    USE ldfdyn_oce     ! lateral diffusion: dynamics 
     24   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    2725   USE phycst         ! physical constants 
    2826   USE zdfmxl         ! mixed layer depth 
     
    3028   ! 
    3129   USE in_out_manager ! I/O manager 
     30   USE prtctl         ! Print control 
    3231   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    33    USE prtctl         ! Print control 
     32   USE lib_mpp        ! distribued memory computing library 
     33   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3434   USE wrk_nemo       ! work arrays 
    3535   USE timing         ! Timing 
    36    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3736 
    3837   IMPLICIT NONE 
    3938   PRIVATE 
    4039 
    41    PUBLIC   ldf_slp        ! routine called by step.F90 
    42    PUBLIC   ldf_slp_grif   ! routine called by step.F90 
    43    PUBLIC   ldf_slp_init   ! routine called by opa.F90 
    44  
    45    LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag 
    46    !                                                                             !! Madec operator 
    47    !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 
     40   PUBLIC   ldf_slp         ! routine called by step.F90 
     41   PUBLIC   ldf_slp_triad   ! routine called by step.F90 
     42   PUBLIC   ldf_slp_init    ! routine called by nemogcm.F90 
     43 
     44   LOGICAL , PUBLIC ::   l_ldfslp = .FALSE.     !: slopes flag 
     45 
     46   LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.       !: iso-neutral direction 
     47   LOGICAL , PUBLIC ::   ln_traldf_triad = .FALSE.      !: griffies triad scheme 
     48 
     49   LOGICAL , PUBLIC ::   ln_triad_iso    = .FALSE.      !: pure horizontal mixing in ML 
     50   LOGICAL , PUBLIC ::   ln_botmix_triad = .FALSE.      !: mixing on bottom 
     51   REAL(wp), PUBLIC ::   rn_sw_triad     = 1._wp        !: =1 switching triads ; =0 all four triads used  
     52   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp      !: slope limit 
     53 
     54   LOGICAL , PUBLIC ::   l_grad_zps = .FALSE.           !: special treatment for Horz Tgradients w partial steps (triad operator) 
     55    
     56   !                                                     !! Classic operator (Madec) 
    4857   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   uslp, wslpi          !: i_slope at U- and W-points 
    4958   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   vslp, wslpj          !: j-slope at V- and W-points 
    50    !                                                                !! Griffies operator 
     59   !                                                     !! triad operator (Griffies) 
    5160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wslp2                !: wslp**2 from Griffies quarter cells 
    5261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials 
    5362   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate 
    54  
    55    !                                                              !! Madec operator 
    56    !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 
     63   !                                                     !! both operators 
     64   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   ah_wslp2             !: ah * slope^2 at w-point 
     65   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   akz                  !: stabilizing vertical diffusivity 
     66    
     67   !                                                     !! Madec operator 
    5768   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   omlmask           ! mask of the surface mixed layer at T-pt 
    5869   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer 
     
    6374   !! * Substitutions 
    6475#  include "domzgr_substitute.h90" 
    65 #  include "ldftra_substitute.h90" 
    66 #  include "ldfeiv_substitute.h90" 
    6776#  include "vectopt_loop_substitute.h90" 
    6877   !!---------------------------------------------------------------------- 
    69    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     78   !! NEMO/OPA 4.0 , NEMO Consortium (2014) 
    7079   !! $Id$ 
    7180   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    105114      INTEGER  ::   ii0, ii1, iku   ! temporary integer 
    106115      INTEGER  ::   ij0, ij1, ikv   ! temporary integer 
    107       REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars 
     116      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 
    108117      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
    109118      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
    110119      REAL(wp) ::   zck, zfk,      zbw             !   -      - 
    111       REAL(wp) ::   zdepv, zdepu         !   -      - 
    112120      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww 
    113121      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr 
    114122      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zgru, zgrv 
    115       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhmlpu, zhmlpv 
    116123      !!---------------------------------------------------------------------- 
    117124      ! 
     
    119126      ! 
    120127      CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
    121       CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv ) 
    122  
    123       IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN  
    124       
    125          zeps   =  1.e-20_wp        !==   Local constant initialization   ==! 
    126          z1_16  =  1.0_wp / 16._wp 
    127          zm1_g  = -1.0_wp / grav 
    128          zm1_2g = -0.5_wp / grav 
    129          ! 
    130          zww(:,:,:) = 0._wp 
    131          zwz(:,:,:) = 0._wp 
    132          ! 
    133          DO jk = 1, jpk             !==   i- & j-gradient of density   ==! 
    134             DO jj = 1, jpjm1 
    135                DO ji = 1, fs_jpim1   ! vector opt. 
    136                   zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
    137                   zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
    138                END DO 
    139             END DO 
    140          END DO 
    141          IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level 
    142             DO jj = 1, jpjm1 
    143                DO ji = 1, jpim1 
    144                   zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    145                   zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
    146                END DO 
    147             END DO 
    148          ENDIF 
    149          IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
    150             DO jj = 1, jpjm1 
    151                DO ji = 1, jpim1 
    152                IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
    153                IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
    154                END DO 
    155             END DO 
    156          ENDIF 
    157          ! 
    158          !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
    159          ! interior value 
    160          DO jk = 2, jpkm1 
    161             !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
    162             !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
    163             !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1 
    164             !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2 
    165             !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster 
    166             zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              & 
    167                &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 
    168          END DO 
    169          ! surface initialisation  
    170          zdzr(:,:,1) = 0._wp  
    171          IF ( ln_isfcav ) THEN 
    172             ! if isf need to overwrite the interior value at at the first ocean point 
    173             DO jj = 1, jpjm1 
    174                DO ji = 1, jpim1 
    175                   zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
    176                END DO 
    177             END DO 
    178          END IF 
    179          ! 
    180          !                          !==   Slopes just below the mixed layer   ==! 
    181          CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml 
    182  
    183  
    184          ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd ) 
    185          ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    186          ! 
    187          IF ( ln_isfcav ) THEN 
    188             DO jj = 2, jpjm1 
    189                DO ji = fs_2, fs_jpim1   ! vector opt. 
    190                   IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
    191                   IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp) 
    192                   IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp) 
    193                   IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
    194                   IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp) 
    195                   IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp) 
    196                ENDDO 
    197             ENDDO 
    198          ELSE 
    199             DO jj = 2, jpjm1 
    200                DO ji = fs_2, fs_jpim1   ! vector opt. 
    201                   zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
    202                   zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
    203                ENDDO 
    204             ENDDO 
    205          END IF 
    206          DO jk = 2, jpkm1                            !* Slopes at u and v points 
    207             DO jj = 2, jpjm1 
    208                DO ji = fs_2, fs_jpim1   ! vector opt. 
    209                   !                                      ! horizontal and vertical density gradient at u- and v-points 
    210                   zau = zgru(ji,jj,jk) / e1u(ji,jj) 
    211                   zav = zgrv(ji,jj,jk) / e2v(ji,jj) 
    212                   zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) ) 
    213                   zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) ) 
    214                   !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    215                   !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    216                   zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  ) 
    217                   zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
    218                   !                                      ! uslp and vslp output in zwz and zww, resp. 
    219                   zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) ) 
    220                   zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) ) 
    221                   ! thickness of water column between surface and level k at u/v point 
    222                   zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              & 
    223                                    - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) ) 
    224                   zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              & 
    225                                    - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 
    226                   ! 
    227                   zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          & 
    228                      &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 
    229                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 
    230                   zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          & 
    231                      &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)  
    232                   zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 
    233                    
    234                   
     128 
     129      zeps   =  1.e-20_wp        !==   Local constant initialization   ==! 
     130      z1_16  =  1.0_wp / 16._wp 
     131      zm1_g  = -1.0_wp / grav 
     132      zm1_2g = -0.5_wp / grav 
     133      z1_slpmax = 1._wp / rn_slpmax 
     134      ! 
     135      zww(:,:,:) = 0._wp 
     136      zwz(:,:,:) = 0._wp 
     137      ! 
     138      DO jk = 1, jpk             !==   i- & j-gradient of density   ==! 
     139         DO jj = 1, jpjm1 
     140            DO ji = 1, fs_jpim1   ! vector opt. 
     141               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
     142               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
     143            END DO 
     144         END DO 
     145      END DO 
     146      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level 
     147         DO jj = 1, jpjm1 
     148            DO ji = 1, jpim1 
     149               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
     150               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     151            END DO 
     152         END DO 
     153      ENDIF 
     154      ! 
     155      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     156      DO jk = 2, jpkm1 
     157         !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
     158         !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
     159         !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1 
     160         !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2 
     161         !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster 
     162         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              & 
     163            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 
     164      END DO 
     165      ! 
     166      !                          !==   Slopes just below the mixed layer   ==! 
     167      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml 
     168 
     169 
     170      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd ) 
     171      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
     172      ! 
     173      DO jk = 2, jpkm1                            !* Slopes at u and v points 
     174         DO jj = 2, jpjm1 
     175            DO ji = fs_2, fs_jpim1   ! vector opt. 
     176               !                                      ! horizontal and vertical density gradient at u- and v-points 
     177               zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 
     178               zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 
     179               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) ) 
     180               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) ) 
     181               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
     182               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     183               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  ) 
     184               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
     185               !                                      ! uslp and vslp output in zwz and zww, resp. 
     186               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
     187               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
     188               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
     189                  &                   + zfi  * uslpml(ji,jj)                                                     & 
     190                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   & 
     191                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 
     192               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
     193                  &                   + zfj  * vslpml(ji,jj)                                                     & 
     194                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   & 
     195                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 
    235196!!gm  modif to suppress omlmask.... (as in Griffies case) 
    236 !                  !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0. 
    237 !                  zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 
    238 !                  zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 
    239 !                  zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 
    240 !                  zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 
    241 !                  zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 
    242 !                  zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 
     197!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0. 
     198!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 
     199!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 
     200!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 
     201!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 
     202!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 
     203!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 
    243204!!gm end modif 
    244                END DO 
    245             END DO 
    246          END DO 
    247          CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions 
    248          ! 
    249          !                                            !* horizontal Shapiro filter 
    250          DO jk = 2, jpkm1 
    251             DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
    252                DO ji = 2, jpim1 
    253                   uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
    254                      &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
    255                      &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
    256                      &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
    257                      &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
    258                   vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
    259                      &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
    260                      &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      & 
    261                      &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    & 
    262                      &                       + 4.*  zww(ji,jj    ,jk)                       ) 
    263                END DO 
    264             END DO 
    265             DO jj = 3, jpj-2                               ! other rows 
    266                DO ji = fs_2, fs_jpim1   ! vector opt. 
    267                   uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
    268                      &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
    269                      &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
    270                      &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
    271                      &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
    272                   vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
    273                      &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
    274                      &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)         & 
    275                      &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    & 
    276                      &                       + 4.*  zww(ji,jj    ,jk)                       ) 
    277                END DO 
    278             END DO 
    279             !                                        !* decrease along coastal boundaries 
    280             DO jj = 2, jpjm1 
    281                DO ji = fs_2, fs_jpim1   ! vector opt. 
    282                   uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    283                      &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
    284                      &                            *   umask(ji,jj,jk-1) 
    285                   vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    286                      &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   & 
    287                      &                            *   vmask(ji,jj,jk-1) 
    288                END DO 
    289             END DO 
    290          END DO 
    291  
    292  
    293          ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd ) 
    294          ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd ) 
    295          ! 
    296          DO jk = 2, jpkm1 
    297             DO jj = 2, jpjm1 
    298                DO ji = fs_2, fs_jpim1   ! vector opt. 
    299                   !                                  !* Local vertical density gradient evaluated from N^2 
    300                   zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 
    301                   !                                  !* Slopes at w point 
    302                   !                                        ! i- & j-gradient of density at w-points 
    303                   zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           & 
    304                      &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj) 
    305                   zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           & 
    306                      &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj) 
    307                   zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           & 
    308                      &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci 
    309                   zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           & 
    310                      &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj 
    311                   !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    312                   !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    313                   zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  ) 
    314                   zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
    315                   !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    316                   zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    317                   zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
    318                   zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
    319                      &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    320                   zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 
    321                      &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
     205            END DO 
     206         END DO 
     207      END DO 
     208      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions 
     209      ! 
     210      !                                            !* horizontal Shapiro filter 
     211      DO jk = 2, jpkm1 
     212         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
     213            DO ji = 2, jpim1 
     214               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
     215                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
     216                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
     217                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
     218                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
     219               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
     220                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
     221                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      & 
     222                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    & 
     223                  &                       + 4.*  zww(ji,jj    ,jk)                       ) 
     224            END DO 
     225         END DO 
     226         DO jj = 3, jpj-2                               ! other rows 
     227            DO ji = fs_2, fs_jpim1   ! vector opt. 
     228               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
     229                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
     230                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
     231                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
     232                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
     233               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
     234                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
     235                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      & 
     236                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    & 
     237                  &                       + 4.*  zww(ji,jj    ,jk)                       ) 
     238            END DO 
     239         END DO 
     240         !                                        !* decrease along coastal boundaries 
     241         DO jj = 2, jpjm1 
     242            DO ji = fs_2, fs_jpim1   ! vector opt. 
     243               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
     244                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp 
     245               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
     246                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp 
     247            END DO 
     248         END DO 
     249      END DO 
     250 
     251 
     252      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd ) 
     253      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd ) 
     254      ! 
     255      DO jk = 2, jpkm1 
     256         DO jj = 2, jpjm1 
     257            DO ji = fs_2, fs_jpim1   ! vector opt. 
     258               !                                  !* Local vertical density gradient evaluated from N^2 
     259               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 
     260               !                                  !* Slopes at w point 
     261               !                                        ! i- & j-gradient of density at w-points 
     262               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           & 
     263                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj) 
     264               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           & 
     265                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) * e2t(ji,jj) 
     266               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           & 
     267                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk) 
     268               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           & 
     269                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk) 
     270               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
     271               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     272               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  ) 
     273               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
     274               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
     275               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
     276               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 
     277               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
     278               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
    322279 
    323280!!gm  modif to suppress omlmask....  (as in Griffies operator) 
    324 !                  !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0. 
    325 !                  zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 
    326 !                  zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. ) 
    327 !                  zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 
    328 !                  zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 
     281!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0. 
     282!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 
     283!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. ) 
     284!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 
     285!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 
    329286!!gm end modif 
    330                END DO 
    331             END DO 
    332          END DO 
    333          CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions 
    334          ! 
    335          !                                           !* horizontal Shapiro filter 
    336          DO jk = 2, jpkm1 
    337             DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
    338                DO ji = 2, jpim1 
    339                   zcofw = tmask(ji,jj,jk) * z1_16 
    340                   wslpi(ji,jj,jk) = (          zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    341                        &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    342                        &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    343                        &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    344                        &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
    345  
    346                   wslpj(ji,jj,jk) = (          zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
    347                        &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
    348                        &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    349                        &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    350                        &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
    351                END DO 
    352             END DO 
    353             DO jj = 3, jpj-2                               ! other rows 
    354                DO ji = fs_2, fs_jpim1   ! vector opt. 
    355                   zcofw = tmask(ji,jj,jk) * z1_16 
    356                   wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    357                        &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    358                        &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    359                        &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    360                        &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
    361  
    362                   wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
    363                        &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
    364                        &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    365                        &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    366                        &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
    367                END DO 
    368             END DO 
    369             !                                        !* decrease along coastal boundaries 
    370             DO jj = 2, jpjm1 
    371                DO ji = fs_2, fs_jpim1   ! vector opt. 
    372                   zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
    373                      &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
    374                   wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 
    375                   wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 
    376                END DO 
    377             END DO 
    378          END DO 
    379  
    380          ! III.  Specific grid points 
    381          ! =========================== 
    382          ! 
    383          IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area 
    384             !                                                    ! Gibraltar Strait 
    385             ij0 =  50   ;   ij1 =  53 
    386             ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    387             ij0 =  51   ;   ij1 =  53 
    388             ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    389             ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    390             ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    391             ! 
    392             !                                                    ! Mediterrannean Sea 
    393             ij0 =  49   ;   ij1 =  56 
    394             ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    395             ij0 =  50   ;   ij1 =  56 
    396             ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    397             ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    398             ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    399          ENDIF 
    400  
    401  
    402          ! IV. Lateral boundary conditions 
    403          ! =============================== 
    404          CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
    405          CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
    406  
    407  
    408          IF(ln_ctl) THEN 
    409             CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk) 
    410             CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 
    411          ENDIF 
    412          ! 
    413  
    414       ELSEIF ( lk_vvl ) THEN  
    415   
    416          IF(lwp) THEN  
    417             WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'  
    418          ENDIF  
    419  
    420          ! geopotential diffusion in s-coordinates on tracers and/or momentum  
    421          ! The slopes of s-surfaces are computed at each time step due to vvl  
    422          ! The slopes for momentum diffusion are i- or j- averaged of those on tracers  
    423  
    424          ! set the slope of diffusion to the slope of s-surfaces  
    425          !      ( c a u t i o n : minus sign as fsdep has positive value )  
    426          DO jj = 2, jpjm1  
    427             DO ji = fs_2, fs_jpim1   ! vector opt.  
    428                uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1)  
    429                vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1)  
    430                wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5  
    431                wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5  
    432             END DO  
    433          END DO  
    434  
    435          DO jk = 2, jpk  
    436             DO jj = 2, jpjm1  
    437                DO ji = fs_2, fs_jpim1   ! vector opt.  
    438                   uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)  
    439                   vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
    440                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 
    441                     &                              * wmask(ji,jj,jk) * 0.5  
    442                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 
    443                     &                              * wmask(ji,jj,jk) * 0.5  
    444                END DO  
    445             END DO  
    446          END DO  
    447  
    448          ! Lateral boundary conditions on the slopes  
    449          CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )  
    450          CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )  
    451    
    452          if( kt == nit000 ) then  
    453             IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &  
    454                &                             ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj))  
    455          endif  
    456    
    457          IF(ln_ctl) THEN  
    458             CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)  
    459             CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)  
    460          ENDIF  
    461  
     287            END DO 
     288         END DO 
     289      END DO 
     290      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions 
     291      ! 
     292      !                                           !* horizontal Shapiro filter 
     293      DO jk = 2, jpkm1 
     294         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
     295            DO ji = 2, jpim1 
     296               zcofw = wmask(ji,jj,jk) * z1_16 
     297               wslpi(ji,jj,jk) = (         zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
     298                    &               +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
     299                    &               + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
     300                    &               +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
     301                    &               + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
     302 
     303               wslpj(ji,jj,jk) = (         zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
     304                    &               +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
     305                    &               + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
     306                    &               +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
     307                    &               + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
     308            END DO 
     309         END DO 
     310         DO jj = 3, jpj-2                               ! other rows 
     311            DO ji = fs_2, fs_jpim1   ! vector opt. 
     312               zcofw = wmask(ji,jj,jk) * z1_16 
     313               wslpi(ji,jj,jk) = (         zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
     314                    &               +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
     315                    &               + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
     316                    &               +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
     317                    &               + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
     318 
     319               wslpj(ji,jj,jk) = (         zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
     320                    &               +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
     321                    &               + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
     322                    &               +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
     323                    &               + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
     324            END DO 
     325         END DO 
     326         !                                        !* decrease in vicinity of topography 
     327         DO jj = 2, jpjm1 
     328            DO ji = fs_2, fs_jpim1   ! vector opt. 
     329               zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
     330                  &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
     331               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 
     332               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 
     333            END DO 
     334         END DO 
     335      END DO 
     336 
     337      ! IV. Lateral boundary conditions 
     338      ! =============================== 
     339      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
     340      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
     341 
     342 
     343      IF(ln_ctl) THEN 
     344         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk) 
     345         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 
    462346      ENDIF 
    463        
     347      ! 
    464348      CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
    465       CALL wrk_dealloc( jpi,jpj,     zhmlpu, zhmlpv) 
    466349      ! 
    467350      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp') 
     
    470353 
    471354 
    472    SUBROUTINE ldf_slp_grif ( kt ) 
    473       !!---------------------------------------------------------------------- 
    474       !!                 ***  ROUTINE ldf_slp_grif  *** 
     355   SUBROUTINE ldf_slp_triad ( kt ) 
     356      !!---------------------------------------------------------------------- 
     357      !!                 ***  ROUTINE ldf_slp_triad  *** 
    475358      !! 
    476359      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope 
    477       !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 
     360      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T) 
    478361      !!      at W-points using the Griffies quarter-cells. 
    479362      !! 
     
    490373      REAL(wp) ::   zfacti, zfactj              ! local scalars 
    491374      REAL(wp) ::   znot_thru_surface           ! local scalars 
    492       REAL(wp) ::   zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 
     375      REAL(wp) ::   zdit, zdis, zdkt, zbu, zbti, zisw 
     376      REAL(wp) ::   zdjt, zdjs, zdks, zbv, zbtj, zjsw 
    493377      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 
    494378      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 
    495379      REAL(wp) ::   zdzrho_raw 
     380      REAL(wp) ::   zbeta0, ze3_e1, ze3_e2 
    496381      REAL(wp), POINTER, DIMENSION(:,:)     ::   z1_mlbw 
     382      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalbet 
    497383      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients 
    498384      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only 
    499385      !!---------------------------------------------------------------------- 
    500386      ! 
    501       IF( nn_timing == 1 )  CALL timing_start('ldf_slp_grif') 
     387      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_triad') 
    502388      ! 
    503389      CALL wrk_alloc( jpi,jpj, z1_mlbw ) 
     390      CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 
    504391      CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    505392      CALL wrk_alloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
     
    519406                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point 
    520407                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    521                   zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) / e1u(ji,jj) 
    522                   zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj) 
    523                   zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    524                   zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     408                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 
     409                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 
     410                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN(  MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     411                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN(  MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
    525412               END DO 
    526413            END DO 
     
    533420                  zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature 
    534421                  zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity 
    535                   zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) / e1u(ji,jj) 
    536                   zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj) 
     422                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 
     423                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 
    537424                  zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    538425                  zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     
    554441                     zdks = 0._wp 
    555442                  ENDIF 
    556                   zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp) 
    557                   zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw )    ! force zdzrho >= repsln 
     443                  zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 
     444                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw )    ! force zdzrho >= repsln 
    558445                 END DO 
    559446            END DO 
     
    588475                  ! 
    589476                  jk = nmln(ji+ip,jj) + 1 
    590                   IF( jk .GT. mbkt(ji+ip,jj) ) THEN  !ML reaches bottom 
    591                     zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp 
    592                   ELSE 
    593                     ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 
    594                     zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      & 
    595                        &      - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk) 
    596                     zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
     477                  IF( jk > mbkt(ji+ip,jj) ) THEN   ! ML reaches bottom 
     478                     zti_mlb(ji+ip,jj   ,1-ip,kp) = 0.0_wp 
     479                  ELSE                              
     480                     ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 
     481                     zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      & 
     482                        &          - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk) 
     483                     ze3_e1    =  fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj)  
     484                     zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1  , ABS( zti_g_raw ) ), zti_g_raw ) 
    597485                  ENDIF 
    598486                  ! 
    599487                  jk = nmln(ji,jj+jp) + 1 
    600488                  IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom 
    601                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp 
     489                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp 
    602490                  ELSE 
    603                     ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      & 
    604                        &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
    605                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     491                     ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      & 
     492                        &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
     493                     ze3_e2    =  fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj) 
     494                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2  , ABS( ztj_g_raw ) ), ztj_g_raw ) 
    606495                  ENDIF 
    607496               END DO 
     
    632521                     zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked 
    633522                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp) 
    634  
     523                     ! 
    635524                     ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 
    636                      zti_coord = znot_thru_surface * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    637                      ztj_coord = znot_thru_surface * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)                  ! unmasked 
     525                     zti_coord = znot_thru_surface * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) 
     526                     ztj_coord = znot_thru_surface * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj)     ! unmasked 
    638527                     zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces 
    639528                     ztj_g_raw = ztj_raw - ztj_coord 
    640                      zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
    641                      ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     529                     ! additional limit required in bilaplacian case 
     530                     ze3_e1    = fse3w(ji+ip,jj   ,jk+kp) * r1_e1u(ji,jj) 
     531                     ze3_e2    = fse3w(ji   ,jj+jp,jk+kp) * r1_e2v(ji,jj) 
     532                     ! NB: hard coded factor 5 (can be a namelist parameter...) 
     533                     zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 
     534                     ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 
    642535                     ! 
    643536                     ! Below  ML use limited zti_g as is & mask 
     
    668561                     ! 
    669562                     IF( ln_triad_iso ) THEN 
    670                         zti_raw = zti_lim**2 / zti_raw 
    671                         ztj_raw = ztj_lim**2 / ztj_raw 
     563                        zti_raw = zti_lim*zti_lim / zti_raw 
     564                        ztj_raw = ztj_lim*ztj_lim / ztj_raw 
    672565                        zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 
    673566                        ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 
    674                         zti_lim =           zfacti   * zti_lim                       & 
    675                         &      + ( 1._wp - zfacti ) * zti_raw 
    676                         ztj_lim =           zfactj   * ztj_lim                       & 
    677                         &      + ( 1._wp - zfactj ) * ztj_raw 
     567                        zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 
     568                        ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 
    678569                     ENDIF 
    679                      triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim 
    680                      triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim 
    681                     ! 
    682                      zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   ) 
    683                      zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   ) 
    684                      zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 
    685                      zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 
    686                      ! 
    687                      !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked 
    688                      wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2      ! masked 
    689                      wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2 
     570                     !                                      ! switching triad scheme  
     571                     zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad    & 
     572                        &            * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) )  ) 
     573                     zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad    & 
     574                        &            * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) )  ) 
     575                     ! 
     576                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim * zisw 
     577                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 
     578                     ! 
     579                     zbu  = e1e2u(ji   ,jj   ) * fse3u(ji   ,jj   ,jk   ) 
     580                     zbv  = e1e2v(ji   ,jj   ) * fse3v(ji   ,jj   ,jk   ) 
     581                     zbti = e1e2t(ji+ip,jj   ) * fse3w(ji+ip,jj   ,jk+kp) 
     582                     zbtj = e1e2t(ji   ,jj+jp) * fse3w(ji   ,jj+jp,jk+kp) 
     583                     ! 
     584                     wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim      ! masked 
     585                     wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 
    690586                  END DO 
    691587               END DO 
     
    699595      ! 
    700596      CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 
     597      CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 
    701598      CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    702599      CALL wrk_dealloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
    703600      ! 
    704       IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_grif') 
    705       ! 
    706    END SUBROUTINE ldf_slp_grif 
     601      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_triad') 
     602      ! 
     603   END SUBROUTINE ldf_slp_triad 
    707604 
    708605 
     
    730627      INTEGER  ::   ji , jj , jk                   ! dummy loop indices 
    731628      INTEGER  ::   iku, ikv, ik, ikm1             ! local integers 
    732       REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars 
     629      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars 
    733630      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
    734631      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
     
    741638      zm1_g  = -1.0_wp / grav 
    742639      zm1_2g = -0.5_wp / grav 
     640      z1_slpmax = 1._wp / rn_slpmax 
    743641      ! 
    744642      uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp 
     
    752650            DO ji = 1, jpi 
    753651               ik = nmln(ji,jj) - 1 
    754                IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 
    755                   omlmask(ji,jj,jk) = 1._wp 
    756                ELSE 
    757                   omlmask(ji,jj,jk) = 0._wp 
     652               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
     653               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp 
    758654               ENDIF 
    759655            END DO 
     
    777673            ! 
    778674            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point) 
    779             iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
    780             ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
     675            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
     676            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
    781677            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) ) 
    782678            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) ) 
    783679            !                        !- horizontal density gradient at u- & v-points 
    784             zau = p_gru(ji,jj,iku) / e1u(ji,jj) 
    785             zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 
     680            zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 
     681            zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 
    786682            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    787683            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    788             zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  ) 
    789             zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  ) 
     684            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  ) 
     685            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  ) 
    790686            !                        !- Slope at u- & v-points (uslpml, vslpml) 
    791687            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 
     
    812708            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  ) 
    813709            !                        !- i- & j-slope at w-points (wslpiml, wslpjml) 
    814             wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik) 
    815             wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik) 
     710            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 
     711            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 
    816712         END DO 
    817713      END DO 
     
    831727      !! ** Purpose :   Initialization for the isopycnal slopes computation 
    832728      !! 
    833       !! ** Method  :   read the nammbf namelist and check the parameter 
    834       !!      values called by tra_dmp at the first timestep (nit000) 
     729      !! ** Method  :    
    835730      !!---------------------------------------------------------------------- 
    836731      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    845740         WRITE(numout,*) '~~~~~~~~~~~~' 
    846741      ENDIF 
    847  
    848       IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes 
    849          ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 
    850          ALLOCATE( triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1)                      , STAT=ierr ) 
    851          IF( ierr > 0             )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 
    852          ! 
     742      ! 
     743      ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr ) 
     744      IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' ) 
     745      ! 
     746      IF( ln_traldf_triad ) THEN        ! Griffies operator : triad of slopes 
     747         IF(lwp) WRITE(numout,*) '              Griffies (triad) operator initialisation' 
     748         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) ,     & 
     749            &      triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1) ,     & 
     750            &      wslp2   (jpi,jpj,jpk)                                         , STAT=ierr ) 
     751         IF( ierr > 0      )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 
    853752         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 
    854753         ! 
    855754      ELSE                             ! Madec operator : slopes at u-, v-, and w-points 
    856          ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                & 
    857             &   omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj)   , vslpml(jpi,jpj)    , wslpiml(jpi,jpj)   , wslpjml(jpi,jpj) , STAT=ierr ) 
     755         IF(lwp) WRITE(numout,*) '              Madec operator initialisation' 
     756         ALLOCATE( omlmask(jpi,jpj,jpk) ,                                                                        & 
     757            &      uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) ,     & 
     758            &      vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr ) 
    858759         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' ) 
    859760 
     
    865766         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp 
    866767 
    867          IF(ln_sco .AND.  (ln_traldf_hor .OR. ln_dynldf_hor )) THEN 
    868             IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
    869  
    870             ! geopotential diffusion in s-coordinates on tracers and/or momentum 
    871             ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 
    872             ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 
    873  
    874             ! set the slope of diffusion to the slope of s-surfaces 
    875             !      ( c a u t i o n : minus sign as fsdep has positive value ) 
    876             DO jk = 1, jpk 
    877                DO jj = 2, jpjm1 
    878                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    879                      uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
    880                      vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
    881                      wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
    882                      wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
    883                   END DO 
    884                END DO 
    885             END DO 
    886             CALL lbc_lnk( uslp , 'U', -1. )   ;   CALL lbc_lnk( vslp , 'V', -1. )      ! Lateral boundary conditions 
    887             CALL lbc_lnk( wslpi, 'W', -1. )   ;   CALL lbc_lnk( wslpj, 'W', -1. ) 
    888          ENDIF 
     768         !!gm I no longer understand this..... 
     769!!gm         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 
     770!            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
     771! 
     772!            ! geopotential diffusion in s-coordinates on tracers and/or momentum 
     773!            ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 
     774!            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 
     775! 
     776!            ! set the slope of diffusion to the slope of s-surfaces 
     777!            !      ( c a u t i o n : minus sign as fsdep has positive value ) 
     778!            DO jk = 1, jpk 
     779!               DO jj = 2, jpjm1 
     780!                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     781!                     uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
     782!                     vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
     783!                     wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 
     784!                     wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 
     785!                  END DO 
     786!               END DO 
     787!            END DO 
     788!            CALL lbc_lnk( uslp , 'U', -1. )   ;   CALL lbc_lnk( vslp , 'V', -1. )      ! Lateral boundary conditions 
     789!            CALL lbc_lnk( wslpi, 'W', -1. )   ;   CALL lbc_lnk( wslpj, 'W', -1. ) 
     790!!gm         ENDIF 
    889791      ENDIF 
    890792      ! 
     
    892794      ! 
    893795   END SUBROUTINE ldf_slp_init 
    894  
    895 #else 
    896    !!------------------------------------------------------------------------ 
    897    !!   Dummy module :                 NO Rotation of lateral mixing tensor 
    898    !!------------------------------------------------------------------------ 
    899    LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag 
    900 CONTAINS 
    901    SUBROUTINE ldf_slp( kt, prd, pn2 )   ! Dummy routine 
    902       INTEGER, INTENT(in) :: kt 
    903       REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 
    904       WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) 
    905    END SUBROUTINE ldf_slp 
    906    SUBROUTINE ldf_slp_grif( kt )        ! Dummy routine 
    907       INTEGER, INTENT(in) :: kt 
    908       WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 
    909    END SUBROUTINE ldf_slp_grif 
    910    SUBROUTINE ldf_slp_init              ! Dummy routine 
    911    END SUBROUTINE ldf_slp_init 
    912 #endif 
    913796 
    914797   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.