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 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 13 years ago

set proper svn properties to all files...

  • Property svn:keywords set to Id
File size: 50.7 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   !!----------------------------------------------------------------------
12#if   defined key_ldfslp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
15   !!----------------------------------------------------------------------
16   !!   ldf_slp      : compute the slopes of neutral surface
17   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface
18   !!   ldf_slp_init : initialization of the slopes computation
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE ldftra_oce
23   USE ldfdyn_oce
24   USE phycst          ! physical constants
25   USE zdfmxl          ! mixed layer depth
26   USE eosbn2
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE in_out_manager  ! I/O manager
29   USE prtctl          ! Print control
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   ldf_slp         ! routine called by step.F90
35   PUBLIC   ldf_slp_init    ! routine called by opa.F90
36   PUBLIC ldf_slp_grif   !              "
37
38   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag
39   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points
40   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points
41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   wslp2                !: wslp**2 from Griffies quarter cells
42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   alpha, beta          !: alpha,beta at T points
43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd
44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psix_eiv
45   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psiy_eiv
46   
47   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt
48   REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
49   REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
50
51   !! * Substitutions
52#  include "domzgr_substitute.h90"
53#  include "ldftra_substitute.h90"
54#  include "ldfeiv_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
58   !! $Id$
59   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt)
60   !!----------------------------------------------------------------------
61
62CONTAINS
63
64   SUBROUTINE ldf_slp( kt, prd, pn2 )
65      !!----------------------------------------------------------------------
66      !!                 ***  ROUTINE ldf_slp  ***
67      !!
68      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
69      !!              surfaces referenced locally) ('key_traldfiso').
70      !!
71      !! ** Method  :   The slope in the i-direction is computed at U- and
72      !!      W-points (uslp, wslpi) and the slope in the j-direction is
73      !!      computed at V- and W-points (vslp, wslpj).
74      !!      They are bounded by 1/100 over the whole ocean, and within the
75      !!      surface layer they are bounded by the distance to the surface
76      !!      ( slope<= depth/l  where l is the length scale of horizontal
77      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
78      !!      of 10cm/s)
79      !!        A horizontal shapiro filter is applied to the slopes
80      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
81      !!      the slope of the model level surface.
82      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
83      !!      [slopes already set to zero at level 1, and to zero or the ocean
84      !!      bottom slope (ln_sco=T) at level jpk in inildf]
85      !!
86      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
87      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
88      !!----------------------------------------------------------------------
89      USE oce , zgru  => ua   ! use ua as workspace
90      USE oce , zgrv  => va   ! use va as workspace
91      USE oce , zwy   => ta   ! use ta as workspace
92      USE oce , zwz   => sa   ! use sa as workspace
93      !!
94      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
95      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   prd   ! in situ density
96      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
97      !!
98      INTEGER  ::   ji , jj , jk    ! dummy loop indices
99      INTEGER  ::   ii0, ii1, iku   ! temporary integer
100      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
101      REAL(wp) ::   zeps, zmg, zm05g, zalpha        ! temporary scalars
102      REAL(wp) ::   zcoef1, zcoef2, zcoef3          !    -         -
103      REAL(wp) ::   zcofu , zcofv , zcofw           !    -         -
104      REAL(wp) ::   zau, zbu, zai, zbi, z1u, z1wu   !    -         -
105      REAL(wp) ::   zav, zbv, zaj, zbj, z1v, z1wv   !
106      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww     ! 3D workspace
107      !!----------------------------------------------------------------------
108     
109      zeps  =  1.e-20                             ! Local constant initialization
110      zmg   = -1.0 / grav
111      zm05g = -0.5 / grav
112      !
113      zww(:,:,:) = 0.e0
114      zwz(:,:,:) = 0.e0
115      !                                           ! horizontal density gradient computation
116      DO jk = 1, jpk
117         DO jj = 1, jpjm1
118            DO ji = 1, fs_jpim1   ! vector opt.
119               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
120               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
121            END DO
122         END DO
123      END DO
124      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
125# if defined key_vectopt_loop 
126         DO jj = 1, 1
127            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
128# else
129         DO jj = 1, jpjm1
130            DO ji = 1, jpim1
131# endif
132               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1               ! last ocean level
133               ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
134               zgru(ji,jj,iku) = gru(ji,jj) 
135               zgrv(ji,jj,ikv) = grv(ji,jj)               
136            END DO
137         END DO
138      ENDIF
139     
140      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
141
142     
143      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
144      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
145      !               
146      !                                           !* Local vertical density gradient evaluated from N^2
147      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
148         DO jj = 1, jpj
149            DO ji = 1, jpi
150               zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   &
151                  &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. )
152            END DO
153         END DO
154      END DO
155      DO jk = 2, jpkm1                            !* Slopes at u and v points
156         DO jj = 2, jpjm1
157            DO ji = fs_2, fs_jpim1   ! vector opt.
158               ! horizontal and vertical density gradient at u- and v-points
159               zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk)
160               zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk)
161               zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj  ,jk) )
162               zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji  ,jj+1,jk) )
163               ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
164               !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
165               zbu = MIN( zbu, -100.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) )
166               zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) )
167               ! uslp and vslp output in zwz and zww, resp.
168               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
169               zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           &
170                  &                    + zalpha  * uslpml(ji,jj)                                                  &
171                  &                              * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
172                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
173               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
174               zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           &
175                  &                    + zalpha  * vslpml(ji,jj)                                                  &
176                  &                              * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
177                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
178            END DO
179         END DO
180      END DO
181      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
182      !
183      zcofu = 1. / 16.                            !* horizontal Shapiro filter
184      zcofv = 1. / 16.
185      DO jk = 2, jpkm1
186         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
187            DO ji = 2, jpim1 
188               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
189                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
190                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
191                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
192                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
193               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
194                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
195                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
196                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
197                  &                       + 4.*  zww(ji,jj    ,jk)                       )
198            END DO
199         END DO
200         DO jj = 3, jpj-2                               ! other rows
201            DO ji = fs_2, fs_jpim1   ! vector opt.
202               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
203                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
204                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
205                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
206                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
207               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
208                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
209                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
210                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
211                  &                       + 4.*  zww(ji,jj    ,jk)                       )
212            END DO
213         END DO
214         !                                        !* decrease along coastal boundaries
215         DO jj = 2, jpjm1
216            DO ji = fs_2, fs_jpim1   ! vector opt.
217               z1u  = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5
218               z1v  = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5
219               z1wu = ( umask(ji,jj,jk)   + umask(ji,jj,jk+1) )*.5
220               z1wv = ( vmask(ji,jj,jk)   + vmask(ji,jj,jk+1) )*.5
221               uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu
222               vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv
223            END DO
224         END DO
225      END DO
226
227
228      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
229      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
230      !               
231      !                                           !* Local vertical density gradient evaluated from N^2
232      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
233         DO jj = 1, jpj
234            DO ji = 1, jpi
235               zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
236            END DO
237         END DO
238      END DO
239      DO jk = 2, jpkm1                            !* Slopes at w point
240         DO jj = 2, jpjm1
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               !                                        ! horizontal density i-gradient at w-points
243               zcoef1 = MAX( zeps, umask(ji-1,jj,jk  )+umask(ji,jj,jk  )    &
244                  &               +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )
245               zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
246               zai = zcoef1 * (  zgru(ji  ,jj,jk  ) + zgru(ji  ,jj,jk-1)   &
247                  &            + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk  ) ) * tmask (ji,jj,jk)
248               !                                        ! horizontal density j-gradient at w-points
249               zcoef2 = MAX( zeps, vmask(ji,jj-1,jk  )+vmask(ji,jj,jk-1)   &
250                  &               +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk  ) )
251               zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
252               zaj = zcoef2 * (  zgrv(ji,jj  ,jk  ) + zgrv(ji,jj  ,jk-1)   &
253                  &            + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk  ) ) * tmask (ji,jj,jk)
254               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
255               !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
256               zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) )
257               zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) )
258               !                                        ! wslpi and wslpj output in zwz and zww, resp.
259               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
260               zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
261               zwz(ji,jj,jk) = (     zai / ( zbi - zeps)  * ( 1. - zalpha )   &
262                  &             + zcoef3 * wslpiml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
263               zww(ji,jj,jk) = (     zaj / ( zbj - zeps)  * ( 1. - zalpha )   &
264                  &             + zcoef3 * wslpjml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
265            END DO
266         END DO
267      END DO
268      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
269      !
270      !                                           !* horizontal Shapiro filter
271      DO jk = 2, jpkm1
272         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
273            DO ji = 2, jpim1
274               zcofw = tmask(ji,jj,jk) / 16.
275               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
276                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
277                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
278                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
279                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
280
281               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
282                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
283                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
284                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
285                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
286            END DO
287         END DO 
288         DO jj = 3, jpj-2                               ! other rows
289            DO ji = fs_2, fs_jpim1   ! vector opt.
290               zcofw = tmask(ji,jj,jk) / 16.
291               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
292                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
293                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
294                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
295                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
296
297               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
298                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
299                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
300                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
301                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
302            END DO
303         END DO
304         !                                        !* decrease along coastal boundaries
305         DO jj = 2, jpjm1
306            DO ji = fs_2, fs_jpim1   ! vector opt.
307               z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5
308               z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5
309               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v
310               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v
311            END DO
312         END DO
313      END DO
314         
315      ! III.  Specific grid points     
316      ! ===========================
317      !               
318      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
319         !                                                    ! Gibraltar Strait
320         ij0 =  50   ;   ij1 =  53
321         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
322         ij0 =  51   ;   ij1 =  53
323         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
324         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
325         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
326         !
327         !                                                    ! Mediterrannean Sea
328         ij0 =  49   ;   ij1 =  56
329         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
330         ij0 =  50   ;   ij1 =  56
331         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
332         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
333         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
334      ENDIF
335
336     
337      ! IV. Lateral boundary conditions
338      ! ===============================
339      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
340      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
341
342
343      IF(ln_ctl) THEN
344         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
345         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
346      ENDIF
347      !
348   END SUBROUTINE ldf_slp
349
350   SUBROUTINE ldf_slp_grif ( kt )
351     !!----------------------------------------------------------------------
352     !!                 ***  ROUTINE ldf_slp_grif  ***
353     !!
354     !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
355     !!      of iso-pycnal surfaces referenced locally) ('key_traldfiso')
356     !!      at W-points using the Griffies quarter-cells.  Also calculates
357     !!      alpha and beta at T-points for use in the Griffies isopycnal
358     !!      scheme.
359     !!
360     !! ** Method  :   The slope in the i-direction is computed at U- and
361     !!      W-points (uslp, wslpi) and the slope in the j-direction is
362     !!      computed at V- and W-points (vslp, wslpj).
363     !!
364     !! ** Action : - alpha, beta
365     !!               wslp2 squared slope of neutral surfaces at w-points.
366     !!
367     !! History :
368     !!   9.0  !  06-10  (C. Harris)  New subroutine
369     !!----------------------------------------------------------------------
370     !! * Modules used
371     USE oce            , zdit  => ua,  &  ! use ua as workspace
372          zdis  => va,  &  ! use va as workspace
373          zdjt  => ta,  &  ! use ta as workspace
374          zdjs  => sa      ! use sa as workspace
375     !! * Arguments
376     INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
377
378     !! * Local declarations
379     INTEGER  ::   ji, jj, jk, ip, jp, kp  ! dummy loop indices
380     INTEGER  ::   iku, ikv                ! temporary integer
381     REAL(wp) ::   &
382          zt, zs, zh, zt2, zsp5, zp1t1,   &  ! temporary scalars
383          zdenr, zrhotmp, zdndt, zdddt,   &  !     "        "
384          zdnds, zddds, znum, zden,       &  !     "        "
385          zslope, za_sxe, zslopec, zdsloper,&  !     "        "
386          zfact, zepsln, zatempw,zatempu,zatempv, &   !     "        "
387          ze1ur,ze2vr,ze3wr,zdxt,zdxs,zdyt,zdys,zdzt,zdzs,zvolf,&
388          zr_slpmax,zdxrho,zdyrho,zabs_dzrho
389     REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) ::   &
390          zsx,zsy
391     REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   &
392          zsx_ml_base,zsy_ml_base
393     REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
394          zdkt,zdks
395     REAL(wp), DIMENSION(jpi,jpj) ::   &
396          zr_ml_basew
397     !!----------------------------------------------------------------------
398
399     ! 0. Local constant initialization
400     ! --------------------------------
401     zr_slpmax = 1.0_wp/slpmax
402
403     ! zslopec=0.004
404     ! zdsloper=1000.0
405     zepsln=1e-25
406
407     SELECT CASE ( nn_eos )
408
409     CASE ( 0 )               ! Jackett and McDougall (1994) formulation
410
411        !                                                ! ===============
412        DO jk = 1, jpk                                   ! Horizontal slab
413           !                                             ! ===============
414           DO jj = 1, jpjm1
415              DO ji = 1, fs_jpim1
416                 zt = tb(ji,jj,jk)          ! potential temperature
417                 zs = sb(ji,jj,jk) - 35.0   ! salinity anomaly (s-35)
418                 zh = fsdept(ji,jj,jk)      ! depth in meters
419
420                 beta(ji,jj,jk) = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &
421                      &                            - 0.301985e-05 ) * zt      &
422                      &   + 0.785567e-03                                      &
423                      &   + (     0.515032e-08 * zs                           &
424                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
425                      &   +(  (   0.121551e-17 * zh                           &
426                      &         - 0.602281e-15 * zs                           &
427                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
428                      &                             + 0.408195e-10   * zs     &
429                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
430                      &                             - 0.121555e-07 ) * zh
431
432                 alpha(ji,jj,jk) = - beta(ji,jj,jk) *                       &
433                      &     (((( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &
434                      &                               - 0.203814e-03 ) * zt   &
435                      &                               + 0.170907e-01 ) * zt   &
436                      &   + 0.665157e-01                                      &
437                      &   +     ( - 0.678662e-05 * zs                         &
438                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
439                      &   +   ( ( - 0.302285e-13 * zh                         &
440                      &           - 0.251520e-11 * zs                         &
441                      &           + 0.512857e-12 * zt * zt           ) * zh   &
442                      &           - 0.164759e-06 * zs                         &
443                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
444                      &                               + 0.380374e-04 ) * zh)
445
446              ENDDO
447           ENDDO
448        ENDDO
449
450     CASE ( 1 )
451
452        alpha(:,:,:)=-rn_alpha
453        beta(:,:,:)=0.0
454
455     CASE ( 2 )
456
457        alpha(:,:,:)=-rn_alpha
458        beta (:,:,:)=rn_beta
459
460     CASE ( 3 )
461
462        DO jk =1, jpk
463           DO jj = 1, jpjm1
464              DO ji = 1, fs_jpim1
465                 zt = tb(ji,jj,jk)
466                 zs = sb(ji,jj,jk)
467                 zh = fsdept(ji,jj,jk)
468                 zt2 = zt**2
469                 zsp5 = sqrt(ABS(zs))
470                 zp1t1=zh*zt
471                 znum=9.99843699e+02+zt*(7.35212840e+00+zt*(-5.45928211e-02+3.98476704e-04*zt)) &
472                      +zs*(2.96938239e+00-7.23268813e-03*zt+2.12382341e-03*zs)  &
473                      +zh*(1.04004591e-02+1.03970529e-07*zt2+5.18761880e-06*zs+ &
474                      zh*(-3.24041825e-08-1.23869360e-11*zt2))
475                 zden=1.00000000e+00+zt*(7.28606739e-03+zt*(-4.60835542e-05+zt*(3.68390573e-07+zt*1.80809186e-10))) &
476                      +zs*(2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(4.76534122e-06+1.63410736e-09*zt2)) &
477                      + zh*(5.30848875e-06+zh*zt*(-3.03175128e-16*zt2-1.27934137e-17*zh))
478                 zdenr=1.0/zden
479                 zrhotmp=znum*zdenr
480                 zdndt=7.35212840e+00+zt*(-1.091856422e-01+1.195430112e-03*zt)-7.23268813e-03*zs &
481                      +zp1t1*(2.07941058e-07-2.4773872e-11*zh)
482                 zdddt=7.28606739e-03+zt*(-9.21671084e-05+zt*(1.105171719e-06+7.23236744e-10*zt))  &
483                      +zs*(-9.27062484e-06+zt*(-5.35030929e-10*zt+3.26821472e-09*zsp5))  &
484                      +zh*zh*(-9.09525384e-16*zt2-1.27934137e-17*zh)
485                 zdnds=2.96938239e+00-7.23268813e-03*zt+2*2.12382341e-03*zs+5.18761880e-06*zh
486                 zddds=2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(7.14801183e-06+2.45116104e-09*zt2)
487                 alpha(ji,jj,jk)=(zdndt-zrhotmp*zdddt)*zdenr
488                 beta(ji,jj,jk)=zdenr*(zdnds-zrhotmp*zddds)
489
490              END DO
491           END DO
492        END DO
493
494     CASE DEFAULT
495
496        IF(lwp) WRITE(numout,cform_err)
497        IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
498        nstop = nstop + 1
499
500     END SELECT
501
502     CALL lbc_lnk( alpha, 'T', 1. )
503     CALL lbc_lnk( beta, 'T', 1. )
504
505
506     ! ---------------------------------------
507     ! 1. Horizontal tracer gradients at T-level jk
508     ! ---------------------------------------
509     DO jk = 1, jpkm1
510        DO jj = 1, jpjm1
511           DO ji = 1, fs_jpim1   ! vector opt.
512              ! i-gradient of T and S at jj
513              zdit (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk)
514              zdis (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk)
515              ! j-gradient of T and S at jj
516              zdjt (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk)
517              zdjs (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk)
518           END DO
519        END DO
520     END DO
521
522     IF( ln_zps ) THEN      ! partial steps correction at the last level
523# if defined key_vectopt_loop   &&   ! defined key_mpp_omp
524     jj = 1
525        DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
526# else
527           DO jj = 1, jpjm1
528              DO ji = 1, jpim1
529# endif
530                 ! last ocean level
531                 iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
532                 ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
533                 ! i-gradient of T and S
534                 zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem)
535                 zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal)
536                 ! j-gradient of T and S
537                 zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem)
538                 zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal)
539# if ! defined key_vectopt_loop   ||   defined key_mpp_omp
540              END DO
541# endif
542           END DO
543        ENDIF
544
545        ! ---------------------------------------
546        ! 1. Vertical tracer gradient at w-level jk
547        ! ---------------------------------------
548        DO jk = 2, jpk
549           zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk)
550           zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk)
551        END DO
552
553        zdkt(:,:,1) = 0.0
554        zdks(:,:,1) = 0.0
555
556        ! ---------------------------------------
557        ! Depth of the w-point below ML base
558        ! ---------------------------------------
559        DO jj = 1, jpj
560           DO ji = 1, jpi
561              jk = nmln(ji,jj)
562              zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1)
563           END DO
564        END DO
565
566
567        wslp2(:,:,:)=0.0
568        tfw(:,:,:) = 0.0
569        sfw(:,:,:) = 0.0
570        ftu(:,:,:) = 0.0
571        fsu(:,:,:) = 0.0
572        ftv(:,:,:) = 0.0
573        fsv(:,:,:) = 0.0
574        ftud(:,:,:) = 0.0
575        fsud(:,:,:) = 0.0
576        ftvd(:,:,:) = 0.0
577        fsvd(:,:,:) = 0.0
578        psix_eiv(:,:,:) = 0.0
579        psiy_eiv(:,:,:) = 0.0
580
581        ! ----------------------------------------------------------------------
582        ! x-z plane
583        ! ----------------------------------------------------------------------
584
585        ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1)
586        DO ip=0,1
587           DO kp=0,1
588
589              DO jk = 1, jpkm1
590                 DO jj = 1, jpjm1
591                    DO ji = 1, fs_jpim1   ! vector opt.
592
593                       ze1ur=1.0/e1u(ji,jj)
594                       zdxt=zdit(ji,jj,jk)*ze1ur
595                       zdxs=zdis(ji,jj,jk)*ze1ur
596
597                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp)
598                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr
599                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr
600                       ! Calculate the horizontal and vertical density gradient
601                       zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs
602                       zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln
603                       ! Limit by slpmax, and mask by psi-point
604                       zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) &
605                            & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho))
606                    END DO
607                 END DO
608              END DO
609
610           END DO
611        END DO
612
613        ! calculate slope of triad immediately below mixed-layer base
614        DO ip=0,1
615           DO kp=0,1
616              DO jj = 1, jpjm1
617                 DO ji = 1, fs_jpim1
618                    jk = nmln(ji+ip,jj)
619                    zsx_ml_base(ji+ip,jj,1-ip,kp)=zsx(ji+ip,jj,jk+1-kp,1-ip,kp)
620                 END DO
621              END DO
622           END DO
623        END DO
624
625        ! Below ML use limited zsx as is
626        ! Inside ML replace by linearly reducing zsx_ml_base towards surface
627        DO ip=0,1
628           DO kp=0,1
629
630              DO jk = 1, jpkm1
631
632                 DO jj = 1, jpjm1
633
634                    DO ji = 1, fs_jpim1   ! vector opt.
635                       ! k index of uppermost point(s) of triad is jk+kp-1
636                       ! must be .ge. nmln(ji,jj) for zfact=1.
637                       !                   otherwise  zfact=0.
638                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj))
639                       zsx(ji+ip,jj,jk,1-ip,kp) = zfact*zsx(ji+ip,jj,jk,1-ip,kp) + &
640                            & (1.0_wp-zfact)*(fsdepw(ji+ip,jj,jk+kp)*zr_ml_basew(ji+ip,jj))*zsx_ml_base(ji+ip,jj,1-ip,kp) 
641                    END DO
642
643                 END DO
644
645              END DO
646           END DO
647        END DO
648
649        ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points
650        DO ip=0,1
651           DO kp=0,1
652
653              DO jk = 1, jpkm1
654
655                 DO jj = 1, jpjm1
656
657                    DO ji = 1, fs_jpim1   ! vector opt.
658
659                       ze1ur=1.0/e1u(ji,jj)
660                       zdxt=zdit(ji,jj,jk)*ze1ur
661                       zdxs=zdis(ji,jj,jk)*ze1ur
662
663                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp)
664                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr
665                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr
666                       zslope=zsx(ji+ip,jj,jk,1-ip,kp)
667
668                       zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk)
669
670                       ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur
671                       fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur
672                       ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur
673                       fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur
674                       tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt
675                       sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs
676                       wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ &
677                            & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2
678                       psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zslope
679
680                    END DO
681                 END DO
682
683              END DO
684           END DO
685        END DO
686
687        ! ----------------------------------------------------------------------
688        ! y-z plane
689        ! ----------------------------------------------------------------------
690
691        ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1)
692        DO jp=0,1
693           DO kp=0,1
694
695              DO jk = 1, jpkm1
696                 DO jj = 1, jpjm1
697                    DO ji = 1, fs_jpim1   ! vector opt.
698
699                       ze2vr=1.0/e2v(ji,jj)
700                       zdyt=zdjt(ji,jj,jk)*ze2vr
701                       zdys=zdjs(ji,jj,jk)*ze2vr
702
703                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp)
704                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr
705                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr
706                       ! Calculate the horizontal and vertical density gradient
707                       zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys
708                       zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln
709                       ! Limit by slpmax, and mask by psi-point
710                       zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) &
711                            & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho))
712                    END DO
713                 END DO
714              END DO
715
716           END DO
717        END DO
718
719        ! calculate slope of triad immediately below mixed-layer base
720        DO jp=0,1
721           DO kp=0,1
722              DO jj = 1, jpjm1
723                 DO ji = 1, fs_jpim1
724                    jk = nmln(ji,jj+jp)
725                    zsy_ml_base(ji,jj+jp,1-jp,kp)=zsy(ji,jj+jp,jk+1-kp,1-jp,kp)
726                 END DO
727              END DO
728           END DO
729        END DO
730
731        ! Below ML use limited zsy as is
732        ! Inside ML replace by linearly reducing zsy_ml_base towards surface
733        DO jp=0,1
734           DO kp=0,1
735
736              DO jk = 1, jpkm1
737
738                 DO jj = 1, jpjm1
739
740                    DO ji = 1, fs_jpim1   ! vector opt.
741                       ! k index of uppermost point(s) of triad is jk+kp-1
742                       ! must be .ge. nmln(ji,jj) for zfact=1.
743                       !                   otherwise  zfact=0.
744                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp))
745                       zsy(ji,jj+jp,jk,1-jp,kp) = zfact*zsy(ji,jj+jp,jk,1-jp,kp) + &
746                            & (1.0_wp-zfact)*(fsdepw(ji,jj+jp,jk+kp)*zr_ml_basew(ji,jj+jp))*zsy_ml_base(ji,jj+jp,1-jp,kp) 
747                    END DO
748
749                 END DO
750
751              END DO
752           END DO
753        END DO
754
755        ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points
756        DO jp=0,1
757           DO kp=0,1
758
759              DO jk = 1, jpkm1
760
761                 DO jj = 1, jpjm1
762
763                    DO ji = 1, fs_jpim1   ! vector opt.
764
765                       ze2vr=1.0/e2v(ji,jj)
766                       zdyt=zdjt(ji,jj,jk)*ze2vr
767                       zdys=zdjs(ji,jj,jk)*ze2vr
768
769                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp)
770                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr
771                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr
772                       zslope=zsy(ji,jj+jp,jk,1-jp,kp)
773
774                       zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk)
775
776                       ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr
777                       fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr
778                       ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr
779                       fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr
780                       tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt
781                       sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys
782                       wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ &
783                            & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2
784                       psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zslope
785
786                    END DO
787                 END DO
788
789              END DO
790           END DO
791        END DO
792
793        tfw(:,:,1)=0.0
794        sfw(:,:,1)=0.0
795        wslp2(:,:,1)=0.0
796
797        CALL lbc_lnk( wslp2, 'W', 1. )
798        CALL lbc_lnk( tfw, 'W', 1. )
799        CALL lbc_lnk( sfw, 'W', 1. )
800        CALL lbc_lnk( ftu, 'U', -1. )
801        CALL lbc_lnk( fsu, 'U', -1. )
802        CALL lbc_lnk( ftv, 'V', -1. )
803        CALL lbc_lnk( fsv, 'V', -1. )
804        CALL lbc_lnk( ftud, 'U', -1. )
805        CALL lbc_lnk( fsud, 'U', -1. )
806        CALL lbc_lnk( ftvd, 'V', -1. )
807        CALL lbc_lnk( fsvd, 'V', -1. )
808        CALL lbc_lnk( psix_eiv, 'U', -1. )
809        CALL lbc_lnk( psiy_eiv, 'V', -1. )
810
811
812      END SUBROUTINE ldf_slp_grif
813
814
815   SUBROUTINE ldf_slp_mxl( prd, pn2 )
816      !!----------------------------------------------------------------------
817      !!                  ***  ROUTINE ldf_slp_mxl  ***
818      !!
819      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
820      !!              the mixed layer.
821      !!
822      !! ** Method :
823      !!      The slope in the i-direction is computed at u- and w-points
824      !!   (uslp, wslpi) and the slope in the j-direction is computed at
825      !!   v- and w-points (vslp, wslpj).
826      !!   They are bounded by 1/100 over the whole ocean, and within the
827      !!   surface layer they are bounded by the distance to the surface
828      !!   ( slope<= depth/l  where l is the length scale of horizontal
829      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
830      !!   of 10cm/s)
831      !!
832      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
833      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
834      !!----------------------------------------------------------------------
835      USE oce , zgru  => ua   ! ua, va used as workspace and set to hor.
836      USE oce , zgrv  => va   ! density gradient in ldf_slp
837      !!
838      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density
839      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
840      !!
841      INTEGER  ::   ji, jj, jk           ! dummy loop indices
842      INTEGER  ::   ik, ikm1             ! temporary integers
843      REAL(wp) ::   zeps, zmg, zm05g     ! temporary scalars
844      REAL(wp) ::   zcoef1, zcoef2       !    -         -
845      REAL(wp) ::   zau, zbu, zai, zbi   !    -         -
846      REAL(wp) ::   zav, zbv, zaj, zbj   !    -         -
847      REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace
848      !!----------------------------------------------------------------------
849
850      zeps  =  1.e-20                    ! Local constant initialization
851      zmg   = -1.0 / grav
852      zm05g = -0.5 / grav
853      !
854      uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
855      vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0
856      wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0
857      wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
858
859      !                                  ! surface mixed layer mask
860      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
861# if defined key_vectopt_loop
862         DO jj = 1, 1
863            DO ji = 1, jpij   ! vector opt. (forced unrolling)
864# else
865         DO jj = 1, jpj
866            DO ji = 1, jpi
867# endif
868               ik = nmln(ji,jj) - 1
869               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0
870               ELSE                  ;   omlmask(ji,jj,jk) = 0.e0
871               ENDIF
872            END DO
873         END DO
874      END DO
875
876
877      ! Slopes of isopycnal surfaces just before bottom of mixed layer
878      ! --------------------------------------------------------------
879      ! The slope are computed as in the 3D case.
880      ! A key point here is the definition of the mixed layer at u- and v-points.
881      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
882      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
883      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
884      ! induce velocity field near the base of the mixed layer.
885      !-----------------------------------------------------------------------
886      !
887      zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2)
888      zwy(jpi,:) = 0.e0
889# if defined key_vectopt_loop
890      DO jj = 1, 1
891         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
892# else
893      DO jj = 1, jpjm1
894         DO ji = 1, jpim1
895# endif
896            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation
897            ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values
898            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
899               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
900         END DO
901      END DO
902      CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change
903
904      !                            !* Slope at u points
905# if defined key_vectopt_loop
906      DO jj = 1, 1
907         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
908# else
909      DO jj = 2, jpjm1
910         DO ji = 2, jpim1
911# endif
912            ! horizontal and vertical density gradient at u-points
913            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
914            ik = MIN( ik, jpkm1 )
915            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
916            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
917            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
918            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
919            zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) )
920            ! uslpml
921            uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
922         END DO
923      END DO
924      CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change)
925
926      !                            !* vertical density gradient for v-slope (from N^2)
927# if defined key_vectopt_loop
928      DO jj = 1, 1
929         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
930# else
931      DO jj = 1, jpjm1
932         DO ji = 1, jpim1
933# endif
934            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
935            ik = MIN( ik, jpkm1 )
936            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
937               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
938         END DO
939      END DO
940      CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change
941
942      !                            !* Slope at v points
943# if defined key_vectopt_loop
944      DO jj = 1, 1
945         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
946# else
947      DO jj = 2, jpjm1
948         DO ji = 2, jpim1
949# endif
950            ! horizontal and vertical density gradient at v-points
951            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
952            ik = MIN( ik,jpkm1 )
953            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
954            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
955            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
956            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
957            zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) )
958            ! vslpml
959            vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
960         END DO
961      END DO
962      CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change)
963
964
965      !                            !* vertical density gradient for w-slope (from N^2)
966# if defined key_vectopt_loop
967      DO jj = 1, 1
968         DO ji = 1, jpij   ! vector opt. (forced unrolling)
969# else
970      DO jj = 1, jpj
971         DO ji = 1, jpi
972# endif
973            ik = nmln(ji,jj) + 1
974            ik = MIN( ik, jpk )
975            ikm1 = MAX ( 1, ik-1)
976            zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
977               &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
978         END DO
979      END DO
980
981      !                            !* Slopes at w points
982# if defined key_vectopt_loop
983      DO jj = 1, 1
984         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
985# else
986      DO jj = 2, jpjm1
987         DO ji = 2, jpim1
988# endif
989            ik = nmln(ji,jj) + 1
990            ik = MIN( ik, jpk )
991            ikm1 = MAX ( 1, ik-1 )
992            ! horizontal density i-gradient at w-points
993            zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   &
994               &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
995            zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
996            zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   &
997               &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik)
998            ! horizontal density j-gradient at w-points
999            zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    &
1000               &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) )
1001            zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
1002            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
1003               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
1004            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
1005            !                   static instability:
1006            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
1007            zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
1008            zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) )
1009            ! wslpiml and wslpjml
1010            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
1011            wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
1012         END DO
1013      END DO
1014      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions
1015      !
1016   END SUBROUTINE ldf_slp_mxl
1017
1018
1019   SUBROUTINE ldf_slp_init
1020      !!----------------------------------------------------------------------
1021      !!                  ***  ROUTINE ldf_slp_init  ***
1022      !!
1023      !! ** Purpose :   Initialization for the isopycnal slopes computation
1024      !!
1025      !! ** Method  :   read the nammbf namelist and check the parameter
1026      !!      values called by tra_dmp at the first timestep (nit000)
1027      !!----------------------------------------------------------------------
1028      INTEGER ::   ji, jj, jk   ! dummy loop indices
1029      !!----------------------------------------------------------------------
1030     
1031      IF(lwp) THEN   
1032         WRITE(numout,*)
1033         WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
1034         WRITE(numout,*) '~~~~~~~'
1035      ENDIF
1036
1037      ! Direction of lateral diffusion (tracers and/or momentum)
1038      ! ------------------------------
1039      ! set the slope to zero (even in s-coordinates)
1040
1041      uslp (:,:,:) = 0.e0
1042      vslp (:,:,:) = 0.e0
1043      wslpi(:,:,:) = 0.e0
1044      wslpj(:,:,:) = 0.e0
1045
1046      uslpml (:,:) = 0.e0
1047      vslpml (:,:) = 0.e0
1048      wslpiml(:,:) = 0.e0
1049      wslpjml(:,:) = 0.e0
1050
1051      IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
1052         IF(lwp) THEN
1053            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
1054         ENDIF
1055
1056         ! geopotential diffusion in s-coordinates on tracers and/or momentum
1057         ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
1058         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
1059
1060         ! set the slope of diffusion to the slope of s-surfaces
1061         !      ( c a u t i o n : minus sign as fsdep has positive value )
1062         DO jk = 1, jpk
1063            DO jj = 2, jpjm1
1064               DO ji = fs_2, fs_jpim1   ! vector opt.
1065                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
1066                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
1067                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
1068                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
1069               END DO
1070            END DO
1071         END DO
1072         ! Lateral boundary conditions on the slopes
1073         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
1074         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
1075      ENDIF
1076      !
1077   END SUBROUTINE ldf_slp_init
1078
1079#else
1080   !!------------------------------------------------------------------------
1081   !!   Dummy module :                 NO Rotation of lateral mixing tensor
1082   !!------------------------------------------------------------------------
1083   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
1084CONTAINS
1085   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
1086      INTEGER, INTENT(in) :: kt 
1087      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
1088      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
1089   END SUBROUTINE ldf_slp
1090   SUBROUTINE ldf_slp_init       ! Dummy routine
1091   END SUBROUTINE ldf_slp_init
1092#endif
1093
1094   !!======================================================================
1095END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.