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.
ldfslp.F90 on Ticket #435 – Attachment – NEMO

Ticket #435: ldfslp.F90

File ldfslp.F90, 47.0 KB (added by gm, 13 years ago)

ldfslp for v3.3beta - suppress Griffies operator for v3.2

Line 
1MODULE ldfslp
2   !!======================================================================
3   !!                       ***  MODULE  ldfslp  ***
4   !! Ocean physics: slopes of neutral surfaces
5   !!======================================================================
6   !! History :  OPA  ! 1994-12  (G. Madec, M. Imbard)  Original code
7   !!            8.0  ! 1997-06  (G. Madec)  optimization, lbc
8   !!            8.1  ! 1999-10  (A. Jouzeau)  NEW profile in the mixed layer
9   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
10   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
11   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator
12   !!----------------------------------------------------------------------
13#if   defined key_ldfslp   ||   defined key_esopa
14   !!----------------------------------------------------------------------
15   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
16   !!----------------------------------------------------------------------
17   !!   ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)
18   !!   ldf_slp      : calculates the slopes of neutral surface   (Madec operator)
19   !!   ldf_slp_mxl  : calculates the slopes at the base of the mixed layer (Madec operator)
20   !!   ldf_slp_init : initialization of the slopes computation
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE ldftra_oce      ! lateral diffusion: traceur
25   USE ldfdyn_oce      ! lateral diffusion: dynamics
26   USE phycst          ! physical constants
27   USE zdfmxl          ! mixed layer depth
28   USE eosbn2          ! equation of states
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE in_out_manager  ! I/O manager
31   USE prtctl          ! Print control
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   ldf_slp        ! routine called by step.F90
37   PUBLIC   ldf_slp_grif   ! routine called by step.F90
38   PUBLIC   ldf_slp_init   ! routine called by opa.F90
39
40   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
41   !                                                                             !! Madec operator
42   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   uslp, wslpi          !: i_slope at U- and W-points
43   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   vslp, wslpj          !: j-slope at V- and W-points
44   !                                                                             !! Griffies operator
45   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   wslp2                !: wslp**2 from Griffies quarter cells
46   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials
47   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate
48
49   !                                                              !! Madec operator
50   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   omlmask           ! mask of the surface mixed layer at T-pt
51   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
52   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
53   INTEGER , DIMENSION(:,:)  , ALLOCATABLE ::   mbku  , mbkv      ! vertical index of the (upper) bottom ocean U/V-level
54
55   REAL(wp) ::   repsln = 1.e-25_wp       ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
56
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59#  include "ldftra_substitute.h90"
60#  include "ldfeiv_substitute.h90"
61#  include "vectopt_loop_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
64   !! $Id: ldfslp.F90 2371 2010-11-10 15:38:27Z acc $
65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE ldf_slp( kt, prd, pn2 )
70      !!----------------------------------------------------------------------
71      !!                 ***  ROUTINE ldf_slp  ***
72      !!
73      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
74      !!              surfaces referenced locally) (ln_traldf_iso=T).
75      !!
76      !! ** Method  :   The slope in the i-direction is computed at U- and
77      !!      W-points (uslp, wslpi) and the slope in the j-direction is
78      !!      computed at V- and W-points (vslp, wslpj).
79      !!      They are bounded by 1/100 over the whole ocean, and within the
80      !!      surface layer they are bounded by the distance to the surface
81      !!      ( slope<= depth/l  where l is the length scale of horizontal
82      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
83      !!      of 10cm/s)
84      !!        A horizontal shapiro filter is applied to the slopes
85      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
86      !!      the slope of the model level surface.
87      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
88      !!      [slopes already set to zero at level 1, and to zero or the ocean
89      !!      bottom slope (ln_sco=T) at level jpk in inildf]
90      !!
91      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
92      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
93      !!----------------------------------------------------------------------
94      USE oce , zgru  => ua   ! use ua as workspace
95      USE oce , zgrv  => va   ! use va as workspace
96      USE oce , zww   => ta   ! use ta as workspace
97      USE oce , zwz   => sa   ! use sa as workspace
98      !!
99      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
100      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   prd   ! in situ density
101      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
102      !!
103      INTEGER  ::   ji , jj , jk    ! dummy loop indices
104      INTEGER  ::   ii0, ii1, iku   ! local integer
105      INTEGER  ::   ij0, ij1, ikv   ! local integer
106      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zalpha   ! local scalars
107      REAL(wp) ::   zci, zau, zbu, zai, zbi, zbw         !   -      -
108      REAL(wp) ::   zcj, zav, zbv, zaj, zbj, zck         !   -      -
109      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzr     ! 3D workspace
110      !!----------------------------------------------------------------------
111     
112      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
113      z1_16  =  1.0_wp / 16._wp
114      zm1_g  = -1.0_wp / grav
115      zm1_2g = -0.5_wp / grav
116      !
117      zww(:,:,:) = 0.e0
118      zwz(:,:,:) = 0.e0
119      !
120      !                          !==   i- & j-gradient of density   ==!
121      DO jk = 1, jpk
122         DO jj = 1, jpjm1
123            DO ji = 1, fs_jpim1   ! vector opt.
124               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
125               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
126            END DO
127         END DO
128      END DO
129      IF( ln_zps ) THEN                ! partial steps correction at the bottom ocean level
130# if defined key_vectopt_loop 
131         DO jj = 1, 1
132            DO ji = 1, jpij-jpi   ! Vect. Opt. (forced unrolling)
133# else
134         DO jj = 1, jpjm1
135            DO ji = 1, jpim1
136# endif
137               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1      ! last ocean level
138               ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
139               zgru(ji,jj,iku) = gru(ji,jj) 
140               zgrv(ji,jj,ikv) = grv(ji,jj)               
141            END DO
142         END DO
143      ENDIF
144      !
145      DO jk = 2, jpkm1           !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
146         DO jj = 1, jpj   
147            DO ji = 1, jpi                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
148               !                          !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
149               !                          !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
150               !                          !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
151               !                          ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
152               zdzr(ji,jj,jk) = zm1_g * ( prd(ji,jj,jk) + 1._wp )              &
153                  &                   * ( pn2(ji,jj,jk) + pn2(ji,jj,jk+1) ) * ( 1._wp - 0.5_wp * tmask(ji,jj,jk+1) )
154            END DO
155         END DO
156      END DO
157     
158      !                          !==   Slopes just below the mixed layer   ==!
159      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
160
161     
162      ! I.  slopes at u and v point       | uslp = d/di( prd ) / d/dz( prd )
163      ! ===========================       | vslp = d/dj( prd ) / d/dz( prd )
164      !               
165      DO jk = 2, jpkm1           !==   Slopes at u- & v-points   ==!
166         DO jj = 2, jpjm1
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168               !                                      ! horizontal and vertical density gradient at u- and v-points
169               zau = zgru(ji,jj,jk) / e1u(ji,jj)
170               zav = zgrv(ji,jj,jk) / e2v(ji,jj)
171               zbu = 0.5 * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
172               zbv = 0.5 * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
173               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
174               !                                        kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
175               zbu = MIN(  zbu, -100.* ABS( zau ) , -7.e+3/fse3u(ji,jj,jk)* ABS( zau )  )
176               zbv = MIN(  zbv, -100.* ABS( zav ) , -7.e+3/fse3v(ji,jj,jk)* ABS( zav )  )
177               !                                      ! uslp and vslp output in zwz and zww, resp.
178               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
179               zwz(ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           &
180                  &                   + zalpha  * uslpml(ji,jj)                                                  &
181                  &                             * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
182                  &                             / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
183               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
184               zww(ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           &
185                  &                   + zalpha  * vslpml(ji,jj)                                                  &
186                  &                             * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
187                  &                             / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
188!!gm  modif to suppress omlmask....
189!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
190!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
191!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
192!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
193!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
194!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
195!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
196!!gm end modif
197            END DO
198         END DO
199      END DO
200      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
201      !
202      !                                            !* horizontal Shapiro filter
203      DO jk = 2, jpkm1
204         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
205            DO ji = 2, jpim1   ! NO vector opt.
206               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
207                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
208                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
209                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
210                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
211               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
212                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
213                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
214                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
215                  &                       + 4.*  zww(ji,jj    ,jk)                       )
216            END DO
217         END DO
218         DO jj = 3, jpj-2                               ! other rows
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
221                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
222                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
223                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
224                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
225               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
226                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
227                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
228                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
229                  &                       + 4.*  zww(ji,jj    ,jk)                       )
230            END DO
231         END DO
232         !               
233         DO jj = 2, jpjm1                         !* decrease along coastal boundaries
234            DO ji = fs_2, fs_jpim1   ! vector opt.
235               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5   &
236                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5
237               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5   &
238                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5
239            END DO
240         END DO
241      END DO
242
243
244      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
245      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
246      !               
247      DO jk = 2, jpkm1
248         DO jj = 2, jpjm1
249            DO ji = fs_2, fs_jpim1   ! Vect. Opt.
250               !                                  !* Local vertical density gradient evaluated from N^2
251               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
252               !
253               !                                  !* Slopes at w point
254               !                                        ! i- & j-gradient of density at w-points
255               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           &
256                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj)
257               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           &
258                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj)
259               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
260                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
261               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
262                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
263               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
264               !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
265               zbi = MIN( zbw ,- 100.* ABS( zai ) , -7.e+3/fse3w(ji,jj,jk)* ABS( zai )  )
266               zbj = MIN( zbw , -100.* ABS( zaj ) , -7.e+3/fse3w(ji,jj,jk)* ABS( zaj )  )
267               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
268               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
269               zck    = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
270               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1. - zalpha ) + zck * wslpiml(ji,jj) * zalpha  ) * tmask(ji,jj,jk)
271               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1. - zalpha ) + zck * wslpjml(ji,jj) * zalpha  ) * tmask(ji,jj,jk)
272
273!!gm  modif to suppress omlmask....
274!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
275!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
276!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
277!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
278!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
279!!gm end modif
280            END DO
281         END DO
282      END DO
283      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
284      !
285      !                                           !* horizontal Shapiro filter
286      DO jk = 2, jpkm1
287         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
288            DO ji = 2, jpim1
289               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
290                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
291                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
292                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
293                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
294                  !
295               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
296                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
297                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
298                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
299                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
300            END DO
301         END DO 
302         DO jj = 3, jpj-2                               ! other rows
303            DO ji = fs_2, fs_jpim1   ! vector opt.
304               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
305                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
306                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
307                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
308                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
309
310               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
311                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
312                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
313                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
314                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
315            END DO
316         END DO
317         !                                        !* decrease along coastal boundaries
318         DO jj = 2, jpjm1
319            DO ji = fs_2, fs_jpim1   ! vector opt.
320               zck = (  umask(ji,jj,jk) + umask(ji-1,jj,jk)  &
321                  &   + vmask(ji,jj,jk) + vmask(ji,jj-1,jk)  ) * 0.25
322               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 
323               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
324            END DO
325         END DO
326      END DO
327         
328      ! III.  Specific grid points     
329      ! ===========================
330      !               
331      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
332         !                                                    ! Gibraltar Strait
333         ij0 =  50   ;   ij1 =  53
334         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
335         ij0 =  51   ;   ij1 =  53
336         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
337         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
338         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
339         !
340         !                                                    ! Mediterrannean Sea
341         ij0 =  49   ;   ij1 =  56
342         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
343         ij0 =  50   ;   ij1 =  56
344         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
345         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
346         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
347      ENDIF
348
349      ! IV. Lateral boundary conditions
350      ! ===============================
351      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
352      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
353
354      IF(ln_ctl) THEN
355         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
356         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
357      ENDIF
358      !
359   END SUBROUTINE ldf_slp
360   
361
362   SUBROUTINE ldf_slp_grif ( kt )
363      !!----------------------------------------------------------------------
364      !!                 ***  ROUTINE ldf_slp_grif  ***
365      !!
366      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
367      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T)
368      !!      at W-points using the Griffies quarter-cells. 
369      !!
370      !! ** Method  :   calculates alpha and beta at T-points
371      !!
372      !! ** Action : - triadi_g, triadj_g   T-pts i- and j-slope triads relative to geopot. (used for eiv)
373      !!             - triadi , triadj    T-pts i- and j-slope triads relative to model-coordinate
374      !!             - wslp2              squared slope of neutral surfaces at w-points.
375      !!----------------------------------------------------------------------
376      USE oce,   zalbet  => ua   ! use ua as workspace
377      !!
378      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
379      !!
380      INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices
381      INTEGER  ::   iku, ikv                    ! local integer
382      REAL(wp) ::   zfacti, zfactj              ! local scalars
383      REAL(wp) ::   zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj
384      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim
385      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim
386      REAL(wp) ::   zdzrho_raw
387      REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) ::   zdzrho, zdyrho, zdxrho     ! Horizontal and vertical density gradients
388      REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   zti_mlb, ztj_mlb
389      REAL(wp), DIMENSION(jpi,jpj)         ::   z1_mlbw
390      !!----------------------------------------------------------------------
391
392      !--------------------------------!
393      !  Some preliminary calculation  !
394      !--------------------------------!
395      !
396      CALL eos_alpbet( tsb, zalbet )         !==  before local thermal/haline expension ratio at T-points  ==!
397      !
398      DO jl = 0, 1                           !==  unmasked before density i- j-, k-gradients  ==!
399         !
400         ip = jl   ;   jp = jl   ;   kp = jl       ! guaranteed nonzero gradients ( absolute value larger than repsln)
401         DO jk = 1, jpkm1                          ! done each pair of triad
402            DO jj = 1, jpjm1                       ! NB: not masked ==>  a minimum value is set
403               DO ji = 1, fs_jpim1   ! vector opt.
404                  zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point
405                  zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
406                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point
407                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
408                  IF( jk+kp /= 1 ) THEN                                       ! k-gradient of T & S a jk+kp
409                     zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) )
410                     zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) )
411                  ELSE
412                     zdkt = 0._wp                                             ! 1st level gradient set to zero
413                     zdks = 0._wp
414                  ENDIF
415                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zdis ) / e1u(ji,jj)
416                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zdjs ) / e2v(ji,jj)
417                  zdzrho_raw = ( - zalbet(ji   ,jj   ,jk) * zdkt + zdks ) / fse3w(ji,jj,jk+kp)
418                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MIN(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
419                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MIN(   repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
420                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )    ! force zdzrho >= repsln > 0
421               END DO
422            END DO
423         END DO
424         !
425         IF( ln_zps ) THEN                         ! partial steps: correction of i- & j-grad at ocean last levels
426# if defined key_vectopt_loop
427            DO jj = 1, 1
428               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
429# else
430            DO jj = 1, jpjm1
431               DO ji = 1, jpim1
432# endif
433                  iku  = mbku(ji,jj)          ;   ikv  = mbkv(ji,jj)             ! last ocean level (u- & v-points)
434                  zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature
435                  zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity
436                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zdis ) / e1u(ji,jj)
437                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zdjs ) / e2v(ji,jj)
438                  zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MIN( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign
439                  zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MIN( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
440               END DO
441            END DO
442         ENDIF
443         !
444      END DO
445      !
446      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==!
447         DO ji = 1, jpi                            ! i.e. 1 / (hmld+e3t(nmln))  where hmld=depw(nmln)
448            jk = MIN( nmln(ji,jj), mbathy(ji,jj) - 1 ) + 1     ! MIN in case ML depth is the ocean depth
449            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)
450         END DO
451      END DO
452      !
453      !                                      !==  intialisations to zero  ==!
454      !
455      wslp2  (:,:,:)     = 0._wp                                        ! wslp2 will be cumulated 3D field set to zero
456      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
457      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
458!!gm _iso set to zero missing
459      triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero
460      triadj  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp
461     
462      !-------------------------------------!
463      !  Triads just below the Mixed Layer  !
464      !-------------------------------------!
465      !
466      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
467         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED
468            DO jj = 1, jpjm1
469               DO ji = 1, fs_jpim1
470                  ip = jl   ;   jp = jl
471                  jk = MIN( nmln(ji+ip,jj), mbathy(ji+ip,jj) - 1 ) + 1     ! ML level+1 (MIN in case ML depth is the ocean depth)
472                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
473                     &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk)
474                  jk = MIN( nmln(ji,jj+jp), mbathy(ji,jj+jp) - 1 ) + 1
475                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
476                     &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
477                  zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
478                  ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
479               END DO
480            END DO
481         END DO
482      END DO
483
484      !-------------------------------------!
485      !  Triads with surface limits         !
486      !-------------------------------------!
487      !
488      DO kp = 0, 1                        ! k-index of triads
489         DO jl = 0, 1
490            ip = jl   ;   jp = jl         ! i- and j-indices of triads (i-k and j-k planes)
491            DO jk = 1, jpkm1
492               DO jj = 1, jpjm1
493                  DO ji = 1, fs_jpim1   ! vector opt.
494                     !
495                     ! Calculate slope relative to geopotentials used for GM skew fluxes
496                     ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth)
497                     ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point
498                     ! masked by umask taken at the level of dz(rho)
499                     !
500                     ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti)
501                     !
502                     zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked
503                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp)
504                     zti_coord = ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
505                     ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)                  ! unmasked
506                     zti_g_raw = zti_raw + zti_coord      ! ref to geopot surfaces
507                     ztj_g_raw = ztj_raw + ztj_coord
508                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
509                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
510                     !
511                     ! Below  ML use limited zti_g as is
512                     ! Inside ML replace by linearly reducing sx_mlb towards surface
513                     !
514                     zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp )  ! k index of uppermost point(s) of triad is jk+kp-1
515                     zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp )  ! must be .ge. nmln(ji,jj) for zfact=1
516                     !                                                          !                   otherwise  zfact=0
517                     zti_g_lim =           zfacti   * zti_g_lim                       &
518                        &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
519                        &                           * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj)
520                     ztj_g_lim =           zfactj   * ztj_g_lim                       &
521                        &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
522                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp)
523                     !
524                     triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp)                     ! masked
525                     triadj_g(ji   ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)
526                     !
527                     ! Get coefficients of isoneutral diffusion tensor
528                     ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients)
529                     ! 2. We require that isoneutral diffusion  gives no vertical buoyancy flux
530                     !     i.e. 33 term = (real slope* 31, 13 terms)
531                     ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms
532                     ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2
533                     !
534                     zti_lim  = zti_g_lim - zti_coord    ! remove coordinate slope => relative to coordinate surfaces
535                     ztj_lim  = ztj_g_lim - ztj_coord                                 
536                     zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp)      ! square of limited slopes     ! masked <<==
537                     ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp)
538                     !
539                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   )
540                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   )
541                     zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp)
542                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
543                     !
544                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim2 / zti_raw                                   ! masked
545                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw
546                     !
547!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
548                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2      ! masked
549                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2
550                  END DO
551               END DO
552            END DO
553         END DO
554      END DO
555      !
556      wslp2(:,:,1) = 0._wp                ! force the surface wslp to zero
557     
558      CALL lbc_lnk( wslp2, 'W', 1. )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked
559      !
560   END SUBROUTINE ldf_slp_grif
561
562
563   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )
564      !!----------------------------------------------------------------------
565      !!                  ***  ROUTINE ldf_slp_mxl  ***
566      !!
567      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
568      !!              the mixed layer.
569      !!
570      !! ** Method  :   The slope in the i-direction is computed at u- & w-points
571      !!              (uslp, wslpi) and the slope in the j-direction is computed at
572      !!              v- and w-points (vslp, wslpj).
573      !!   They are bounded by 1/100 over the whole ocean, and within the
574      !!   surface layer they are bounded by the distance to the surface
575      !!   ( slope<= depth/l  where l is the length scale of horizontal
576      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
577      !!   of 10cm/s)
578      !!
579      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
580      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
581      !!----------------------------------------------------------------------
582      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd            ! in situ density
583      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.)
584      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
585      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
586      !!
587      INTEGER  ::   ji , jj , jk         ! dummy loop indices
588      INTEGER  ::   iku, ikv, ik, ikm1   ! local integers
589      REAL(wp) ::   zm1_g , zau, zbu, zai, zbi, zci, zbw    ! local scalars
590      REAL(wp) ::   zm1_2g, zav, zbv, zaj, zbj, zcj, zeps   !   -      -
591      !!----------------------------------------------------------------------
592
593      zeps   =  1.e-20                   ! Local constant initialization
594      zm1_g  = -1.0 / grav
595      zm1_2g = -0.5 / grav
596      !
597      !                                  ! surface mixed layer mask
598      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
599# if defined key_vectopt_loop
600         DO jj = 1, 1
601            DO ji = 1, jpij   ! vector opt. (forced unrolling)
602# else
603         DO jj = 1, jpj
604            DO ji = 1, jpi
605# endif
606               ik = nmln(ji,jj) - 1
607               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0
608               ELSE                  ;   omlmask(ji,jj,jk) = 0.e0
609               ENDIF
610            END DO
611         END DO
612      END DO
613
614
615      ! Slopes of isopycnal surfaces just before bottom of mixed layer
616      ! --------------------------------------------------------------
617      ! The slope are computed as in the 3D case.
618      ! A key point here is the definition of the mixed layer at u- and v-points.
619      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
620      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
621      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
622      ! induce velocity field near the base of the mixed layer.
623      !-----------------------------------------------------------------------
624      !
625# if defined key_vectopt_loop
626      DO jj = 1, 1
627         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
628# else
629      DO jj = 2, jpjm1
630         DO ji = 2, jpim1
631# endif
632            !                    !==   Slope at u- & v-points just below the Mixed Layer   ==!
633            !
634            !                          !- vertical density gradient for u- and v-slopes (from dzr at T-point)
635            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
636            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
637            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,ikv+1) )
638            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv+1) )
639            !                          !- horizontal density gradient at u- & v-points
640            zau = p_gru(ji,jj,iku) / e1u(ji,jj)
641            zav = p_grv(ji,jj,ikv) / e2v(ji,jj)
642            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
643            !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
644            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik)* ABS( zau )  )
645            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav )  )
646            !                          !- Slope at u- & v-points (uslpml, vslpml)
647            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik)
648            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik)
649            !
650            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
651            !
652            ik   = MIN( nmln(ji,jj) + 1, jpk )
653            ikm1 = MAX( 1, ik-1 )
654            !                          !- vertical density gradient for w-slope (from N^2)
655            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
656            !                          !- horizontal density i- & j-gradient at w-points
657            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           &
658               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj) 
659            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ikm1)           &
660               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ik  ) , zeps  ) * e2t(ji,jj)
661            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ikm1)           &
662               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ik  )  ) / zci  * tmask(ji,jj,ik)
663            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ikm1)           &
664               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ik  )  ) / zcj  * tmask(ji,jj,ik)
665            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
666            !                             kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
667            zbi = MIN(  zbw , -100.* ABS( zai ) , -7.e+3/fse3w(ji,jj,ik)* ABS( zai )  )
668            zbj = MIN(  zbw , -100.* ABS( zaj ) , -7.e+3/fse3w(ji,jj,ik)* ABS( zaj )  )
669            !                          !- i- & j-slope at w-points (wslpiml, wslpjml)
670            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
671            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
672         END DO
673      END DO
674!!gm this lbc_lnk should be useless....
675      CALL lbc_lnk( uslpml , 'U', -1. )   ;   CALL lbc_lnk( vslpml , 'V', -1. )     ! lateral boundary cond. (sign change)
676      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )     ! lateral boundary cond. (sign change)
677      !
678   END SUBROUTINE ldf_slp_mxl
679
680
681   SUBROUTINE ldf_slp_init
682      !!----------------------------------------------------------------------
683      !!                  ***  ROUTINE ldf_slp_init  ***
684      !!
685      !! ** Purpose :   Initialization for the isopycnal slopes computation
686      !!
687      !! ** Method  :   read the nammbf namelist and check the parameter
688      !!      values called by tra_dmp at the first timestep (nit000)
689      !!----------------------------------------------------------------------
690      INTEGER ::   ji, jj, jk   ! dummy loop indices
691      INTEGER ::   ierr         ! local integer
692      !!
693      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
694      !!----------------------------------------------------------------------
695     
696      IF(lwp) THEN   
697         WRITE(numout,*)
698         WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing'
699         WRITE(numout,*) '~~~~~~~~~~~~'
700      ENDIF
701     
702      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes
703         !                                   ! allocate the arrays
704         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) ,   &
705            &      triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1) ,                        &
706            &      mbku    (jpi,jpj)             , mbkv    (jpi,jpj)             , STAT=ierr )
707         IF( ierr > 0 ) THEN
708            CALL ctl_stop( 'ldf_slp_init : unable to allocate Griffies operator slope ' )   ;   RETURN
709         ENDIF
710         !
711         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
712         !
713         IF( ( ln_traldf_hor .AND. ln_dynldf_hor ) .AND. ln_sco )   &
714            &     CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )
715         !
716         DO jj = 1, jpjm1                    ! bottom k-index of u- (v-) level
717            DO ji = 1, jpim1
718               mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
719               mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
720            END DO
721         END DO
722         zmbk(:,:) = REAL( mbku(:,:) )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku(:,:) = MAX( INT( zmbk(:,:) ), 1 )
723         zmbk(:,:) = REAL( mbkv(:,:) )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv(:,:) = MAX( INT( zmbk(:,:) ), 1 )
724         !
725      ELSE                             ! Madec operator : slopes at u-, v-, and w-points
726         ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                &
727            &   omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj)   , vslpml(jpi,jpj)    , wslpiml(jpi,jpj)   , wslpjml(jpi,jpj) , STAT=ierr )
728         IF( ierr > 0 ) THEN
729            CALL ctl_stop( 'ldf_slp_init : unable to allocate Madec operator slope ' )   ;   RETURN
730         ENDIF
731
732         ! Direction of lateral diffusion (tracers and/or momentum)
733         ! ------------------------------
734         uslp (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp      ! set the slope to zero (even in s-coordinates)
735         vslp (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp
736         wslpi(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp
737         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
738
739!!gm I no longer understand this.....
740         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
741            IF(lwp) THEN
742               WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
743            ENDIF
744
745            ! geopotential diffusion in s-coordinates on tracers and/or momentum
746            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
747            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
748
749            ! set the slope of diffusion to the slope of s-surfaces
750            !      ( c a u t i o n : minus sign as fsdep has positive value )
751            DO jk = 1, jpk
752               DO jj = 2, jpjm1
753                  DO ji = fs_2, fs_jpim1   ! vector opt.
754                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
755                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
756                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
757                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
758                  END DO
759               END DO
760            END DO
761            ! Lateral boundary conditions on the slopes
762            CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
763            CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
764         ENDIF
765      ENDIF     
766      !
767   END SUBROUTINE ldf_slp_init
768
769#else
770   !!------------------------------------------------------------------------
771   !!   Dummy module :                 NO Rotation of lateral mixing tensor
772   !!------------------------------------------------------------------------
773   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
774CONTAINS
775   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
776      INTEGER, INTENT(in) :: kt 
777      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
778      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
779   END SUBROUTINE ldf_slp
780   SUBROUTINE ldf_slp_init       ! Dummy routine
781   END SUBROUTINE ldf_slp_init
782#endif
783
784   !!======================================================================
785END MODULE ldfslp