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 in branches/2016/dev_v3_6_STABLE_OMP/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2016/dev_v3_6_STABLE_OMP/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 6712

Last change on this file since 6712 was 6712, checked in by mocavero, 8 years ago

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