source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

  • Property svn:keywords set to Id
File size: 54.3 KB
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     1.0  ! 2002-10  (G. Madec)  Free form, F90
10   !!             -   ! 2005-10  (A. Beckmann)  correction for s-coordinates
11   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator
12   !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML
13   !!----------------------------------------------------------------------
14#if   defined key_ldfslp   ||   defined key_esopa
15   !!----------------------------------------------------------------------
16   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
17   !!----------------------------------------------------------------------
18   !!   ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)
19   !!   ldf_slp      : calculates the slopes of neutral surface   (Madec operator)
20   !!   ldf_slp_mxl  : calculates the slopes at the base of the mixed layer (Madec operator)
21   !!   ldf_slp_init : initialization of the slopes computation
22   !!----------------------------------------------------------------------
23   USE timing, ONLY: timing_start, timing_stop
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE ldftra_oce      ! lateral diffusion: traceur
27   USE ldfdyn_oce      ! lateral diffusion: dynamics
28   USE phycst          ! physical constants
29   USE zdfmxl          ! mixed layer depth
30   USE eosbn2          ! equation of states
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE in_out_manager  ! I/O manager
33   USE prtctl          ! Print control
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   ldf_slp        ! routine called by step.F90
39   PUBLIC   ldf_slp_grif   ! routine called by step.F90
40   PUBLIC   ldf_slp_init   ! routine called by opa.F90
41
42   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
43   !                                                                             !! Madec operator
44   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   uslp, wslpi          !: i_slope at U- and W-points
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   vslp, wslpj          !: j-slope at V- and W-points
47   !                                                                !! Griffies operator
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wslp2                !: wslp**2 from Griffies quarter cells
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   &
50             &   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   &
52             &   triadi  , triadj     !: isoneutral slopes relative to model-coordinate
53
54   !                                                              !! Madec operator
55   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   omlmask           ! mask of the surface mixed layer at T-pt
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
58   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
59
60   REAL(wp) ::   repsln = 1.e-25_wp       ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
61
62   ! Workspace arrays for ldf_slp_grif. These could be replaced by several 3D and 2D workspace
63   ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace
64   ! arrays can't be used here because of the zero-indexing of some of the ranks. ARPDBG.
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   zdzrho , zdyrho, zdxrho     ! Horizontal and vertical density gradients
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only
67
68   !! * Control permutation of array indices
69#  include "ldfslp_ftrans.h90"
70!FTRANS zdxrho :I :I :z :
71!FTRANS zdyrho :I :I :z :
72!FTRANS zdzrho :I :I :z :
73#  include "oce_ftrans.h90"
74#  include "dom_oce_ftrans.h90"
75#  include "ldftra_oce_ftrans.h90"
76#  include "ldfdyn_oce_ftrans.h90"
77
78   !! * Substitutions
79#  include "domzgr_substitute.h90"
80#  include "ldftra_substitute.h90"
81#  include "ldfeiv_substitute.h90"
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
85   !! $Id$
86   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
87   !!----------------------------------------------------------------------
88CONTAINS
89
90   INTEGER FUNCTION ldf_slp_alloc()
91      !!----------------------------------------------------------------------
92      !!              ***  FUNCTION ldf_slp_alloc  ***
93      !!----------------------------------------------------------------------
94      !
95      ALLOCATE( zdxrho (jpi,jpj,jpk,0:1) , zti_mlb(jpi,jpj,0:1,0:1) ,     &
96         &      zdyrho (jpi,jpj,jpk,0:1) , ztj_mlb(jpi,jpj,0:1,0:1) ,     &
97         &      zdzrho (jpi,jpj,jpk,0:1)                            , STAT=ldf_slp_alloc )
98         !
99      IF( lk_mpp             )   CALL mpp_sum ( ldf_slp_alloc )
100      IF( ldf_slp_alloc /= 0 )   CALL ctl_warn('ldf_slp_alloc : failed to allocate arrays.')
101      !
102   END FUNCTION ldf_slp_alloc
103
104
105   SUBROUTINE ldf_slp( kt, prd, pn2 )
106      !!----------------------------------------------------------------------
107      !!                 ***  ROUTINE ldf_slp  ***
108      !!
109      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
110      !!              surfaces referenced locally) (ln_traldf_iso=T).
111      !!
112      !! ** Method  :   The slope in the i-direction is computed at U- and
113      !!      W-points (uslp, wslpi) and the slope in the j-direction is
114      !!      computed at V- and W-points (vslp, wslpj).
115      !!      They are bounded by 1/100 over the whole ocean, and within the
116      !!      surface layer they are bounded by the distance to the surface
117      !!      ( slope<= depth/l  where l is the length scale of horizontal
118      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
119      !!      of 10cm/s)
120      !!        A horizontal shapiro filter is applied to the slopes
121      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
122      !!      the slope of the model level surface.
123      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
124      !!      [slopes already set to zero at level 1, and to zero or the ocean
125      !!      bottom slope (ln_sco=T) at level jpk in inildf]
126      !!
127      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
128      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
129      !!----------------------------------------------------------------------
130      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
131      USE oce     , ONLY:   zgru => ua       , zww => va   ! (ua,va) used as workspace
132      USE oce     , ONLY:   zgrv => ta       , zwz => sa   ! (ta,sa) used as workspace
133      USE wrk_nemo, ONLY:   zdzr => wrk_3d_1               ! 3D workspace
134      !! DCSE_NEMO: need additional directives for renamed module variables
135!FTRANS zgru :I :I :z
136!FTRANS zww  :I :I :z
137!FTRANS zgrv :I :I :z
138!FTRANS zwz  :I :I :z
139!FTRANS zdzr :I :I :z
140
141      !!
142      INTEGER , INTENT(in)                   ::   kt    ! ocean time-step index
143!FTRANS prd :I :I :z
144!FTRANS pn2 :I :I :z
145      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   prd   ! in situ density
146      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
147      !!
148      INTEGER  ::   ji , jj , jk    ! dummy loop indices
149      INTEGER  ::   ii0, ii1, iku   ! temporary integer
150      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
151      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16     ! local scalars
152      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
153      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
154      REAL(wp) ::   zck, zfk,      zbw             !   -      -
155      !!----------------------------------------------------------------------
156
157      CALL timing_start('ldf_slp')
158
159      IF( wrk_in_use(3, 1) ) THEN
160         CALL ctl_stop('ldf_slp: requested workspace arrays are unavailable')   ;   RETURN
161      ENDIF
162
163      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
164      z1_16  =  1.0_wp / 16._wp
165      zm1_g  = -1.0_wp / grav
166      zm1_2g = -0.5_wp / grav
167      !
168      zww(:,:,:) = 0._wp
169      zwz(:,:,:) = 0._wp
170      !
171#if defined key_z_first
172      DO jj = 1, jpjm1           !==   i- & j-gradient of density   ==!
173         DO ji = 1, jpim1
174            DO jk = 1, jpkf
175#else
176      DO jk = 1, jpkf            !==   i- & j-gradient of density   ==!
177         DO jj = 1, jpjm1
178            DO ji = 1, fs_jpim1   ! vector opt.
179#endif
180               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
181               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
182            END DO
183         END DO
184      END DO
185      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
186!! DCSE_NEMO: Attention! key_vectopt_loop will break key_z_first
187# if ( defined key_vectopt_loop ) && ! ( defined key_z_first )
188         DO jj = 1, 1
189            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
190# else
191         DO jj = 1, jpjm1
192            DO ji = 1, jpim1
193# endif
194               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
195               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)               
196            END DO
197         END DO
198      ENDIF
199      !
200#if defined key_z_first
201      DO jj = 1, jpj
202         DO ji = 1, jpi
203            zdzr(ji,jj,1) = 0._wp   !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
204            DO jk = 2, jpkfm1
205               zdzr(ji,jj,jk) = zm1_g * ( prd(ji,jj,jk) + 1._wp )              &
206                  &                   * ( pn2(ji,jj,jk) + pn2(ji,jj,jk+1) ) * ( 1._wp - 0.5_wp * tmask(ji,jj,jk+1) )
207            END DO
208         END DO
209      END DO
210#else
211      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
212      DO jk = 2, jpkfm1
213         !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
214         !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
215         !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
216         !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
217         !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
218         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
219            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
220      END DO
221#endif
222      !
223      !                          !==   Slopes just below the mixed layer   ==!
224      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
225
226     
227      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
228      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
229      !               
230#if defined key_z_first
231      DO jj = 2, jpjm1               !* Slopes at u and v points
232         DO ji = 2, jpim1
233            DO jk = 2, jpkfm1
234#else
235      DO jk = 2, jpkfm1                            !* Slopes at u and v points
236         DO jj = 2, jpjm1
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238#endif
239               !                                      ! horizontal and vertical density gradient at u- and v-points
240               zau = zgru(ji,jj,jk) / e1u(ji,jj)
241               zav = zgrv(ji,jj,jk) / e2v(ji,jj)
242               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
243               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
244               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
245               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
246               zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  )
247               zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
248               !                                      ! uslp and vslp output in zwz and zww, resp.
249               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
250               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
251               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
252                  &                   + zfi  * uslpml(ji,jj)                                                     &
253                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
254                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
255               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
256                  &                   + zfj  * vslpml(ji,jj)                                                     &
257                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
258                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
259!!gm  modif to suppress omlmask.... (as in Griffies case)
260!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
261!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
262!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
263!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
264!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
265!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
266!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
267!!gm end modif
268            END DO
269         END DO
270      END DO
271      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
272      !
273      !                                            !* horizontal Shapiro filter
274#if defined key_z_first
275      DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
276         DO ji = 2, jpim1 
277            DO jk = 2, jpkfm1
278#else
279      DO jk = 2, jpkfm1
280         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
281            DO ji = 2, jpim1 
282#endif
283               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
284                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
285                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
286                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
287                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
288               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
289                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
290                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
291                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
292                  &                       + 4.*  zww(ji,jj    ,jk)                       )
293            END DO
294         END DO
295#if defined key_z_first
296      END DO
297      DO jj = 3, jpj-2                               ! other rows
298         DO ji = 2, jpim1
299            DO jk = 2, jpkfm1
300#else
301         DO jj = 3, jpj-2                               ! other rows
302            DO ji = fs_2, fs_jpim1   ! vector opt.
303#endif
304               uslp(ji,jj,jk) = z1_16 * (        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)                       )
309               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
310                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
311                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
312                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
313                  &                       + 4.*  zww(ji,jj    ,jk)                       )
314            END DO
315         END DO
316#if defined key_z_first
317      END DO
318      !                                           !* decrease along coastal boundaries
319      DO jj = 2, jpjm1
320         DO ji = 2, jpim1
321            DO jk = 2, jpkfm1
322#else
323         !                                        !* decrease along coastal boundaries
324         DO jj = 2, jpjm1
325            DO ji = fs_2, fs_jpim1   ! vector opt.
326#endif
327               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
328                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp
329               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
330                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp
331            END DO
332         END DO
333      END DO
334
335
336      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
337      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
338      !               
339#if defined key_z_first
340      DO jj = 2, jpjm1
341         DO ji = 2, jpim1
342            DO jk = 2, jpkfm1
343#else
344      DO jk = 2, jpkfm1
345         DO jj = 2, jpjm1
346            DO ji = fs_2, fs_jpim1   ! vector opt.
347#endif
348               !                                  !* Local vertical density gradient evaluated from N^2
349               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
350               !                                  !* Slopes at w point
351               !                                        ! i- & j-gradient of density at w-points
352               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           &
353                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj)
354               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           &
355                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj)
356               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
357                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk)
358               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
359                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk)
360               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
361               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
362               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  )
363               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
364               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
365               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
366               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
367               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
368               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
369
370!!gm  modif to suppress omlmask....  (as in Griffies operator)
371!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
372!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
373!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
374!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
375!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
376!!gm end modif
377            END DO
378         END DO
379      END DO
380      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions
381      !
382      !                                           !* horizontal Shapiro filter
383#if defined key_z_first
384      DO jj = 2, jpjm1, MAX(1, jpj-3)                           ! rows jj=2 and =jpjm1 only
385         DO ji = 2, jpim1
386            DO jk = 2, jpkfm1
387#else
388      DO jk = 2, jpkfm1
389         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
390            DO ji = 2, jpim1
391#endif
392               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
393                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
394                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
395                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
396                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
397
398               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
399                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
400                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
401                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
402                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
403            END DO
404         END DO 
405#if defined key_z_first
406      END DO
407      DO jj = 3, jpj-2                                  ! other rows
408         DO ji = 2, jpim1
409            DO jk = 2, jpkfm1
410#else
411         DO jj = 3, jpj-2                               ! other rows
412            DO ji = fs_2, fs_jpim1   ! vector opt.
413#endif
414               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
415                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
416                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
417                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
418                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
419
420               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
421                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
422                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
423                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
424                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk)
425            END DO
426         END DO
427#if defined key_z_first
428      END DO
429      !                                           !* decrease along coastal boundaries
430      DO jj = 2, jpjm1
431         DO ji = 2, jpim1
432            DO jk = 2, jpkfm1
433#else
434         !                                        !* decrease along coastal boundaries
435         DO jj = 2, jpjm1
436            DO ji = fs_2, fs_jpim1   ! vector opt.
437#endif
438               zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   &
439                  &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
440               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck
441               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
442            END DO
443         END DO
444      END DO
445         
446      ! III.  Specific grid points     
447      ! ===========================
448      !               
449      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
450         !                                                    ! Gibraltar Strait
451         ij0 =  50   ;   ij1 =  53
452         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
453         ij0 =  51   ;   ij1 =  53
454         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
455         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
456         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
457         !
458         !                                                    ! Mediterranean Sea
459         ij0 =  49   ;   ij1 =  56
460         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
461         ij0 =  50   ;   ij1 =  56
462         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
463         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
464         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , :jpkf ) = 0._wp
465      ENDIF
466
467      ! IV. Lateral boundary conditions
468      ! ===============================
469      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
470      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
471
472
473      IF(ln_ctl) THEN
474         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
475         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
476      ENDIF
477      !
478      IF( wrk_not_released(3, 1) )   CALL ctl_stop('ldf_slp: failed to release workspace arrays')
479      !
480      CALL timing_stop('ldf_slp', 'section')
481
482   END SUBROUTINE ldf_slp
483   
484   !! * Reset control of array index permutation
485!FTRANS CLEAR
486#  include "ldfslp_ftrans.h90"
487!FTRANS zdxrho :I :I :z :
488!FTRANS zdyrho :I :I :z :
489!FTRANS zdzrho :I :I :z :
490#  include "oce_ftrans.h90"
491#  include "dom_oce_ftrans.h90"
492#  include "ldftra_oce_ftrans.h90"
493#  include "ldfdyn_oce_ftrans.h90"
494
495   SUBROUTINE ldf_slp_grif ( kt )
496      !!----------------------------------------------------------------------
497      !!                 ***  ROUTINE ldf_slp_grif  ***
498      !!
499      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
500      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T)
501      !!      at W-points using the Griffies quarter-cells. 
502      !!
503      !! ** Method  :   calculates alpha and beta at T-points
504      !!
505      !! ** Action : - triadi_g, triadj_g   T-pts i- and j-slope triads relative to geopot. (used for eiv)
506      !!             - triadi , triadj    T-pts i- and j-slope triads relative to model-coordinate
507      !!             - wslp2              squared slope of neutral surfaces at w-points.
508      !!----------------------------------------------------------------------
509      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
510      USE oce     , ONLY:   zdit    => ua       , zdis   => va         ! (ua,va) used as workspace
511      USE oce     , ONLY:   zdjt    => ta       , zdjs   => sa         ! (ta,sa) used as workspace
512      USE wrk_nemo, ONLY:   zdkt    => wrk_3d_2 , zdks   => wrk_3d_3   ! 3D workspace
513      USE wrk_nemo, ONLY:   zalpha  => wrk_3d_4 , zbeta => wrk_3d_5    ! alpha, beta at T points, at depth fsgdept
514      USE wrk_nemo, ONLY:   z1_mlbw => wrk_2d_1
515      !! DCSE_NEMO: need additional directives for renamed module variables
516!FTRANS zdit   :I :I :z
517!FTRANS zdis   :I :I :z
518!FTRANS zdjt   :I :I :z
519!FTRANS zdjs   :I :I :z
520!FTRANS zdkt   :I :I :z
521!FTRANS zdks   :I :I :z
522!FTRANS zalpha :I :I :z
523!FTRANS zbeta  :I :I :z
524      !
525      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
526      !
527      INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices
528      INTEGER  ::   iku, ikv                                  ! local integer
529      REAL(wp) ::   zfacti, zfactj, zatempw,zatempu,zatempv   ! local scalars
530      REAL(wp) ::   zbu, zbv, zbti, zbtj                      !   -      -
531      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim
532      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim
533      REAL(wp) ::   zdzrho_raw
534      !!----------------------------------------------------------------------
535
536      IF( wrk_in_use(3, 2,3,4,5) .OR. wrk_in_use(2, 1) )THEN
537         CALL ctl_stop('ldf_slp_grif: requested workspace arrays are unavailable')   ;   RETURN
538      ENDIF
539
540      !--------------------------------!
541      !  Some preliminary calculation  !
542      !--------------------------------!
543      !
544      CALL eos_alpbet( tsb, zalpha, zbeta )     !==  before thermal and haline expension coeff. at T-points  ==!
545      !
546#if defined key_z_first
547      DO jj = 1, jpjm1
548         DO ji = 1, jpim1
549            DO jk = 1, jpkfm1                    !==  before lateral T & S gradients at T-level jk  ==!
550#else
551      DO jk = 1, jpkfm1                          !==  before lateral T & S gradients at T-level jk  ==!
552         DO jj = 1, jpjm1
553            DO ji = 1, fs_jpim1   ! vector opt.
554#endif
555               zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk)   ! i-gradient of T and S at jj
556               zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk)
557               zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk)   ! j-gradient of T and S at jj
558               zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk)
559            END DO
560         END DO
561      END DO
562      IF( ln_zps ) THEN                               ! partial steps: correction at the last level
563# if ( defined key_vectopt_loop ) && ! ( defined key_z_first )
564         DO jj = 1, 1
565            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
566# else
567         DO jj = 1, jpjm1
568            DO ji = 1, jpim1
569# endif
570               zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem)                           ! i-gradient of T and S
571               zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal)
572               zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem)                           ! j-gradient of T and S
573               zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal)
574            END DO
575         END DO
576      ENDIF
577      !
578#if defined key_z_first
579      DO jj = 1, jpj
580         DO ji = 1, jpi
581            zdkt(ji,jj,1) = 0._wp            !==  before vertical T & S gradient at w-level  ==!
582            zdks(ji,jj,1) = 0._wp
583            DO jk = 2, jpkf
584               zdkt(ji,jj,jk) = ( tb(ji,jj,jk-1) - tb(ji,jj,jk) ) * tmask(ji,jj,jk)
585               zdks(ji,jj,jk) = ( sb(ji,jj,jk-1) - sb(ji,jj,jk) ) * tmask(ji,jj,jk)
586            END DO
587         END DO
588      END DO
589#else
590      zdkt(:,:,1) = 0._wp                    !==  before vertical T & S gradient at w-level  ==!
591      zdks(:,:,1) = 0._wp
592      DO jk = 2, jpkf
593         zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk)
594         zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk)
595      END DO
596#endif
597      !
598      !
599      DO jl = 0, 1                           !==  density i-, j-, and k-gradients  ==!
600         ip = jl   ;   jp = jl         ! guaranteed nonzero gradients ( absolute value larger than repsln)
601#if defined key_z_first
602         DO jj = 1, jpjm1                          ! NB: not masked due to the minimum value set
603            DO ji = 1, jpim1
604               DO jk = 1, jpkfm1                    ! done each pair of triad
605#else
606         DO jk = 1, jpkfm1                          ! done each pair of triad
607            DO jj = 1, jpjm1                       ! NB: not masked due to the minimum value set
608               DO ji = 1, fs_jpim1   ! vector opt.
609#endif
610                  zdxrho_raw = ( zalpha(ji+ip,jj   ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj   ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj)
611                  zdyrho_raw = ( zalpha(ji   ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji   ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj)
612                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )    ! keep the sign
613                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX(   repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
614               END DO
615            END DO
616         END DO
617      END DO
618     DO kp = 0, 1                           !==  density i-, j-, and k-gradients  ==!
619#if defined key_z_first
620         DO jj = 1, jpj                         ! NB: not masked due to the minimum value set
621            DO ji = 1, jpi
622               DO jk = 1, jpkfm1                    ! done each pair of triad
623#else
624         DO jk = 1, jpkfm1                          ! done each pair of triad
625            DO jj = 1, jpj                       ! NB: not masked due to the minimum value set
626               DO ji = 1, jpi   ! vector opt.
627#endif
628                  zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   &
629                     &       / fse3w(ji,jj,jk+kp)
630                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )                    ! force zdzrho >= repsln
631               END DO
632            END DO
633         END DO
634      END DO
635      !
636      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==!
637         DO ji = 1, jpi
638            jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
639            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)
640         END DO
641      END DO
642      !
643      !                                      !==  intialisations to zero  ==!
644      !
645      wslp2  (:,:,:)     = 0._wp                                           ! wslp2 will be cumulated 3D field set to zero
646      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp      ! set surface and bottom slope to zero
647      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
648!!gm _iso set to zero missing
649      triadi (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero
650      triadj (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp
651     
652      !-------------------------------------!
653      !  Triads just below the Mixed Layer  !
654      !-------------------------------------!
655      !
656      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
657         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED
658            DO jj = 1, jpjm1
659               DO ji = 1, fs_jpim1
660                  ip = jl   ;   jp = jl
661                  jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1         ! ML level+1 (MIN in case ML depth is the ocean depth)
662                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
663                     &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk)
664                  jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1
665                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
666                     &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
667                  zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
668                  ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
669               END DO
670            END DO
671         END DO
672      END DO
673
674      !-------------------------------------!
675      !  Triads with surface limits         !
676      !-------------------------------------!
677      !
678      DO kp = 0, 1                        ! k-index of triads
679         DO jl = 0, 1
680            ip = jl   ;   jp = jl         ! i- and j-indices of triads (i-k and j-k planes)
681#if defined key_z_first
682            DO jj = 1, jpjm1
683               DO ji = 1, jpim1
684                  DO jk = 1, jpkfm1
685#else
686            DO jk = 1, jpkfm1
687               DO jj = 1, jpjm1
688                  DO ji = 1, fs_jpim1   ! vector opt.
689#endif
690                     !
691                     ! Calculate slope relative to geopotentials used for GM skew fluxes
692                     ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth)
693                     ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point
694                     ! masked by umask taken at the level of dz(rho)
695                     !
696                     ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti)
697                     !
698                     zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked
699                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp)
700                     zti_coord = ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
701                     ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
702                  ! unmasked
703                     zti_g_raw = zti_raw + zti_coord      ! ref to geopot surfaces
704                     ztj_g_raw = ztj_raw + ztj_coord
705                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
706                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
707                     !
708                     ! Below  ML use limited zti_g as is
709                     ! Inside ML replace by linearly reducing sx_mlb towards surface
710                     !
711                     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
712                     zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp )  ! must be .ge. nmln(ji,jj) for zfact=1
713                     !                                                          !                   otherwise  zfact=0
714                     zti_g_lim =           zfacti   * zti_g_lim                       &
715                        &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
716                        &                           * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj)
717                     ztj_g_lim =           zfactj   * ztj_g_lim                       &
718                        &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
719                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp)                   ! masked
720                     !
721                     triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp)
722                     triadj_g(ji   ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)
723                     !
724                     ! Get coefficients of isoneutral diffusion tensor
725                     ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients)
726                     ! 2. We require that isoneutral diffusion  gives no vertical buoyancy flux
727                     !     i.e. 33 term = (real slope* 31, 13 terms)
728                     ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms
729                     ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2
730                     !
731                     zti_lim  = zti_g_lim - zti_coord    ! remove the coordinate slope  ==> relative to coordinate surfaces
732                     ztj_lim  = ztj_g_lim - ztj_coord                                 
733                     zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp)      ! square of limited slopes            ! masked <<==
734                     ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp)
735                     !
736                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   )
737                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   )
738                     zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp)
739                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
740                     !
741                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim2 / zti_raw                                          ! masked
742                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw
743                     !
744!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
745                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2             ! masked
746                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2
747                  END DO
748               END DO
749            END DO
750         END DO
751      END DO
752      !
753      wslp2(:,:,1) = 0._wp                ! force the surface wslp to zero
754     
755      CALL lbc_lnk( wslp2, 'W', 1. )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked
756      !
757      IF( wrk_not_released(3, 2,3,4,5) .OR.   &
758          wrk_not_released(2, 1)       )   CALL ctl_stop('ldf_slp_grif: failed to release workspace arrays')
759      !
760   END SUBROUTINE ldf_slp_grif
761
762   !! * Reset control of array index permutation
763!FTRANS CLEAR
764#  include "ldfslp_ftrans.h90"
765!FTRANS zdxrho :I :I :z :
766!FTRANS zdyrho :I :I :z :
767!FTRANS zdzrho :I :I :z :
768#  include "oce_ftrans.h90"
769#  include "dom_oce_ftrans.h90"
770#  include "ldftra_oce_ftrans.h90"
771#  include "ldfdyn_oce_ftrans.h90"
772
773   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )
774      !!----------------------------------------------------------------------
775      !!                  ***  ROUTINE ldf_slp_mxl  ***
776      !!
777      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
778      !!              the mixed layer.
779      !!
780      !! ** Method  :   The slope in the i-direction is computed at u- & w-points
781      !!              (uslpml, wslpiml) and the slope in the j-direction is computed
782      !!              at v- and w-points (vslpml, wslpjml) with the same bounds as
783      !!              in ldf_slp.
784      !!
785      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces
786      !!                vslpml, wslpjml    just below the mixed layer
787      !!                omlmask         :  mixed layer mask
788      !!----------------------------------------------------------------------
789!FTRANS prd   :I :I :z
790!FTRANS pn2   :I :I :z
791!FTRANS p_gru :I :I :z
792!FTRANS p_grv :I :I :z
793!FTRANS p_dzr :I :I :z
794      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   prd            ! in situ density
795      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.)
796      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
797      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
798      !!
799      INTEGER  ::   ji , jj , jk         ! dummy loop indices
800      INTEGER  ::   iku, ikv, ik, ikm1   ! local integers
801      REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars
802      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
803      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
804      REAL(wp) ::   zck, zfk,      zbw             !   -      -
805      !!----------------------------------------------------------------------
806
807      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
808      zm1_g  = -1.0_wp / grav
809      zm1_2g = -0.5_wp / grav
810      !
811      uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp
812      vslpml (1,:) = 0._wp      ;      vslpml (jpi,:) = 0._wp
813      wslpiml(1,:) = 0._wp      ;      wslpiml(jpi,:) = 0._wp
814      wslpjml(1,:) = 0._wp      ;      wslpjml(jpi,:) = 0._wp
815      !
816      !                          !==   surface mixed layer mask   !
817      DO jk = 1, jpkf                ! =1 inside the mixed layer, =0 otherwise
818# if ( defined key_vectopt_loop ) && ! ( defined key_z_first )
819         DO jj = 1, 1
820            DO ji = 1, jpij   ! vector opt. (forced unrolling)
821# else
822         DO jj = 1, jpj
823            DO ji = 1, jpi
824# endif
825               ik = nmln(ji,jj) - 1
826               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
827               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
828               ENDIF
829            END DO
830         END DO
831      END DO
832
833
834      ! Slopes of isopycnal surfaces just before bottom of mixed layer
835      ! --------------------------------------------------------------
836      ! The slope are computed as in the 3D case.
837      ! A key point here is the definition of the mixed layer at u- and v-points.
838      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
839      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
840      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
841      ! induce velocity field near the base of the mixed layer.
842      !-----------------------------------------------------------------------
843      !
844# if ( defined key_vectopt_loop ) && ! ( defined key_z_first )
845      DO jj = 1, 1
846         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
847# else
848      DO jj = 2, jpjm1
849         DO ji = 2, jpim1
850# endif
851            !                    !==   Slope at u- & v-points just below the Mixed Layer   ==!
852            !
853            !                          !- vertical density gradient for u- and v-slopes (from dzr at T-point)
854            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
855            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
856            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
857            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
858            !                          !- horizontal density gradient at u- & v-points
859            zau = p_gru(ji,jj,iku) / e1u(ji,jj)
860            zav = p_grv(ji,jj,ikv) / e2v(ji,jj)
861            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
862            !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
863            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  )
864            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
865            !                          !- Slope at u- & v-points (uslpml, vslpml)
866            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
867            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
868            !
869            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
870            !
871            ik   = MIN( nmln(ji,jj) + 1, jpkf )
872            ikm1 = MAX( 1, ik-1 )
873            !                          !- vertical density gradient for w-slope (from N^2)
874            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
875            !                          !- horizontal density i- & j-gradient at w-points
876            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           &
877               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj) 
878            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           &
879               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj)
880            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           &
881               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik)
882            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
883               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask(ji,jj,ik)
884            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
885            !                             kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
886            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai )  )
887            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  )
888            !                          !- i- & j-slope at w-points (wslpiml, wslpjml)
889            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
890            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
891         END DO
892      END DO
893!!gm this lbc_lnk should be useless....
894      CALL lbc_lnk( uslpml , 'U', -1. )   ;   CALL lbc_lnk( vslpml , 'V', -1. )   ! lateral boundary cond. (sign change)
895      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )   ! lateral boundary conditions
896      !
897   END SUBROUTINE ldf_slp_mxl
898
899   !! * Reset control of array index permutation
900!FTRANS CLEAR
901#  include "ldfslp_ftrans.h90"
902!FTRANS zdxrho :I :I :z :
903!FTRANS zdyrho :I :I :z :
904!FTRANS zdzrho :I :I :z :
905#  include "oce_ftrans.h90"
906#  include "dom_oce_ftrans.h90"
907#  include "ldftra_oce_ftrans.h90"
908#  include "ldfdyn_oce_ftrans.h90"
909
910   SUBROUTINE ldf_slp_init
911      !!----------------------------------------------------------------------
912      !!                  ***  ROUTINE ldf_slp_init  ***
913      !!
914      !! ** Purpose :   Initialization for the isopycnal slopes computation
915      !!
916      !! ** Method  :   read the nammbf namelist and check the parameter
917      !!              values called by tra_dmp at the first timestep (nit000)
918      !!----------------------------------------------------------------------
919      INTEGER ::   ji, jj, jk   ! dummy loop indices
920      INTEGER ::   ierr         ! local integer
921      !!----------------------------------------------------------------------
922     
923      IF(lwp) THEN   
924         WRITE(numout,*)
925         WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing'
926         WRITE(numout,*) '~~~~~~~~~~~~'
927      ENDIF
928     
929      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes
930         ALLOCATE(triadi_g(jpi,jpj,jpkorig,0:1,0:1), triadj_g(jpi,jpj,jpkorig,0:1,0:1), &
931                  wslp2(jpi,jpj,jpkorig) , STAT=ierr)
932         ALLOCATE(triadi  (jpi,jpj,jpkorig,0:1,0:1), triadj  (jpi,jpj,jpkorig,0:1,0:1), &
933                  STAT=ierr)
934         IF( ierr > 0             )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
935         IF( ldf_slp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate workspace arrays' )
936         !
937         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
938         !
939         IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco )   &
940            CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )
941         !
942      ELSE                             ! Madec operator : slopes at u-, v-, and w-points
943         ALLOCATE(uslp(jpi,jpj,jpkorig), vslp(jpi,jpj,jpkorig),   &
944                  wslpi(jpi,jpj,jpkorig), wslpj(jpi,jpj,jpkorig), &
945                  omlmask(jpi,jpj,jpkorig), &
946                  uslpml(jpi,jpj), vslpml(jpi,jpj), wslpiml(jpi,jpj), wslpjml(jpi,jpj), STAT=ierr)
947         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
948
949         ! Direction of lateral diffusion (tracers and/or momentum)
950         ! ------------------------------
951         uslp (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp      ! set the slope to zero (even in s-coordinates)
952         vslp (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp
953         wslpi(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp
954         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
955
956!!gm I no longer understand this.....
957         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
958            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
959
960            ! geopotential diffusion in s-coordinates on tracers and/or momentum
961            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
962            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
963
964            ! set the slope of diffusion to the slope of s-surfaces
965            !      ( c a u t i o n : minus sign as fsdep has positive value )
966#if defined key_z_first
967            DO jj = 2, jpjm1
968               DO ji = 2, jpim1
969                  DO jk = 1, jpkf
970#else
971            DO jk = 1, jpkf
972               DO jj = 2, jpjm1
973                  DO ji = fs_2, fs_jpim1   ! vector opt.
974#endif
975                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
976                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
977                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
978                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
979                  END DO
980               END DO
981            END DO
982            CALL lbc_lnk( uslp , 'U', -1. )   ;   CALL lbc_lnk( vslp , 'V', -1. )      ! Lateral boundary conditions
983            CALL lbc_lnk( wslpi, 'W', -1. )   ;   CALL lbc_lnk( wslpj, 'W', -1. )
984         ENDIF
985      ENDIF
986      !
987   END SUBROUTINE ldf_slp_init
988
989#else
990   !!------------------------------------------------------------------------
991   !!   Dummy module :                 NO Rotation of lateral mixing tensor
992   !!------------------------------------------------------------------------
993   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
994CONTAINS
995   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
996      INTEGER, INTENT(in) :: kt 
997      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
998      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
999   END SUBROUTINE ldf_slp
1000   SUBROUTINE ldf_slp_grif( kt )        ! Dummy routine
1001      INTEGER, INTENT(in) :: kt
1002      WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt
1003   END SUBROUTINE ldf_slp_grif
1004   SUBROUTINE ldf_slp_init       ! Dummy routine
1005   END SUBROUTINE ldf_slp_init
1006#endif
1007
1008   !!======================================================================
1009END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.