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_r12563_ASINTER-06_ABL_improvement/src/OCE/LDF – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/LDF/ldfslp.F90 @ 12808

Last change on this file since 12808 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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