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_crs.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp_crs.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by cbricaud, 9 years ago

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