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 NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LDF – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LDF/ldfslp.F90 @ 13807

Last change on this file since 13807 was 13807, checked in by mocavero, 4 years ago

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