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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 2371

Last change on this file since 2371 was 2371, checked in by acc, 13 years ago

nemo_v3_3_beta. Improvement of the Griffies operator code (#680). Some aspects are still undergoing scientific assessment but the code structure is ready for release. Dynamical allocation of the large triad arrays has been introduced to reduce memory when the new operator is not in use.

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