New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ldfslp.F90 in branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 2840

Last change on this file since 2840 was 2840, checked in by agn, 13 years ago

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