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 trunk/NEMO/OPA_SRC/LDF – NEMO

source: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 1515

Last change on this file since 1515 was 1515, checked in by ctlod, 15 years ago

wrong sign in the lbc related to pn2 in the computation of the slope just below the mixed layer (step II), see ticket: #435

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 31.9 KB
Line 
1MODULE ldfslp
2   !!======================================================================
3   !!                       ***  MODULE  ldfslp  ***
4   !! Ocean physics: slopes of neutral surfaces
5   !!======================================================================
6   !! History :  OPA  ! 1994-12  (G. Madec, M. Imbard)  Original code
7   !!            8.0  ! 1997-06  (G. Madec)  optimization, lbc
8   !!            8.1  ! 1999-10  (A. Jouzeau)  NEW profile in the mixed layer
9   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
10   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
11   !!----------------------------------------------------------------------
12#if   defined key_ldfslp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
15   !!----------------------------------------------------------------------
16   !!   ldf_slp      : compute the slopes of neutral surface
17   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface
18   !!   ldf_slp_init : initialization of the slopes computation
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE ldftra_oce
23   USE ldfdyn_oce
24   USE phycst          ! physical constants
25   USE zdfmxl          ! mixed layer depth
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE in_out_manager  ! I/O manager
28   USE prtctl          ! Print control
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ldf_slp    ! routine called by step.F90
34
35   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag
36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points
37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points
38   
39   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt
40   REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
41   REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE ldf_slp( kt, prd, pn2 )
55      !!----------------------------------------------------------------------
56      !!                 ***  ROUTINE ldf_slp  ***
57      !!
58      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
59      !!              surfaces referenced locally) ('key_traldfiso').
60      !!
61      !! ** Method  :   The slope in the i-direction is computed at U- and
62      !!      W-points (uslp, wslpi) and the slope in the j-direction is
63      !!      computed at V- and W-points (vslp, wslpj).
64      !!      They are bounded by 1/100 over the whole ocean, and within the
65      !!      surface layer they are bounded by the distance to the surface
66      !!      ( slope<= depth/l  where l is the length scale of horizontal
67      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
68      !!      of 10cm/s)
69      !!        A horizontal shapiro filter is applied to the slopes
70      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
71      !!      the slope of the model level surface.
72      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
73      !!      [slopes already set to zero at level 1, and to zero or the ocean
74      !!      bottom slope (ln_sco=T) at level jpk in inildf]
75      !!
76      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
77      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
78      !!----------------------------------------------------------------------
79      USE oce , zgru  => ua   ! use ua as workspace
80      USE oce , zgrv  => va   ! use va as workspace
81      USE oce , zwy   => ta   ! use ta as workspace
82      USE oce , zwz   => sa   ! use sa as workspace
83      !!
84      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
85      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   prd   ! in situ density
86      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
87      !!
88      INTEGER  ::   ji , jj , jk    ! dummy loop indices
89      INTEGER  ::   ii0, ii1, iku   ! temporary integer
90      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
91      REAL(wp) ::   zeps, zmg, zm05g, zalpha        ! temporary scalars
92      REAL(wp) ::   zcoef1, zcoef2, zcoef3          !    -         -
93      REAL(wp) ::   zcofu , zcofv , zcofw           !    -         -
94      REAL(wp) ::   zau, zbu, zai, zbi, z1u, z1wu   !    -         -
95      REAL(wp) ::   zav, zbv, zaj, zbj, z1v, z1wv   !
96      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww     ! 3D workspace
97      !!----------------------------------------------------------------------
98     
99      IF( kt == nit000 )   CALL ldf_slp_init      ! initialization (first time-step only)
100
101      zeps  =  1.e-20                             ! Local constant initialization
102      zmg   = -1.0 / grav
103      zm05g = -0.5 / grav
104      !
105      zww(:,:,:) = 0.e0
106      zwz(:,:,:) = 0.e0
107      !                                           ! horizontal density gradient computation
108      DO jk = 1, jpk
109         DO jj = 1, jpjm1
110            DO ji = 1, fs_jpim1   ! vector opt.
111               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
112               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
113            END DO
114         END DO
115      END DO
116      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
117# if defined key_vectopt_loop 
118         DO jj = 1, 1
119            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
120# else
121         DO jj = 1, jpjm1
122            DO ji = 1, jpim1
123# endif
124               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1               ! last ocean level
125               ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
126               zgru(ji,jj,iku) = gru(ji,jj) 
127               zgrv(ji,jj,ikv) = grv(ji,jj)               
128            END DO
129         END DO
130      ENDIF
131     
132      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
133
134     
135      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
136      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
137      !               
138      !                                           !* Local vertical density gradient evaluated from N^2
139      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
140         DO jj = 1, jpj
141            DO ji = 1, jpi
142               zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   &
143                  &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. )
144            END DO
145         END DO
146      END DO
147      DO jk = 2, jpkm1                            !* Slopes at u and v points
148         DO jj = 2, jpjm1
149            DO ji = fs_2, fs_jpim1   ! vector opt.
150               ! horizontal and vertical density gradient at u- and v-points
151               zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk)
152               zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk)
153               zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj  ,jk) )
154               zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji  ,jj+1,jk) )
155               ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
156               !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
157               zbu = MIN( zbu, -100.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) )
158               zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) )
159               ! uslp and vslp output in zwz and zww, resp.
160               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
161               zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           &
162                  &                    + zalpha  * uslpml(ji,jj)                                                  &
163                  &                              * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
164                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
165               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
166               zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           &
167                  &                    + zalpha  * vslpml(ji,jj)                                                  &
168                  &                              * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
169                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
170            END DO
171         END DO
172      END DO
173      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
174      !
175      zcofu = 1. / 16.                            !* horizontal Shapiro filter
176      zcofv = 1. / 16.
177      DO jk = 2, jpkm1
178         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
179            DO ji = 2, jpim1 
180               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
181                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
182                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
183                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
184                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
185               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
186                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
187                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
188                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
189                  &                       + 4.*  zww(ji,jj    ,jk)                       )
190            END DO
191         END DO
192         DO jj = 3, jpj-2                               ! other rows
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
195                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
196                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
197                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
198                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
199               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
200                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
201                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
202                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
203                  &                       + 4.*  zww(ji,jj    ,jk)                       )
204            END DO
205         END DO
206         !                                        !* decrease along coastal boundaries
207         DO jj = 2, jpjm1
208            DO ji = fs_2, fs_jpim1   ! vector opt.
209               z1u  = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5
210               z1v  = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5
211               z1wu = ( umask(ji,jj,jk)   + umask(ji,jj,jk+1) )*.5
212               z1wv = ( vmask(ji,jj,jk)   + vmask(ji,jj,jk+1) )*.5
213               uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu
214               vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv
215            END DO
216         END DO
217      END DO
218
219
220      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
221      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
222      !               
223      !                                           !* Local vertical density gradient evaluated from N^2
224      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
225         DO jj = 1, jpj
226            DO ji = 1, jpi
227               zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
228            END DO
229         END DO
230      END DO
231      DO jk = 2, jpkm1                            !* Slopes at w point
232         DO jj = 2, jpjm1
233            DO ji = fs_2, fs_jpim1   ! vector opt.
234               !                                        ! horizontal density i-gradient at w-points
235               zcoef1 = MAX( zeps, umask(ji-1,jj,jk  )+umask(ji,jj,jk  )    &
236                  &               +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )
237               zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
238               zai = zcoef1 * (  zgru(ji  ,jj,jk  ) + zgru(ji  ,jj,jk-1)   &
239                  &            + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk  ) ) * tmask (ji,jj,jk)
240               !                                        ! horizontal density j-gradient at w-points
241               zcoef2 = MAX( zeps, vmask(ji,jj-1,jk  )+vmask(ji,jj,jk-1)   &
242                  &               +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk  ) )
243               zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
244               zaj = zcoef2 * (  zgrv(ji,jj  ,jk  ) + zgrv(ji,jj  ,jk-1)   &
245                  &            + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk  ) ) * tmask (ji,jj,jk)
246               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
247               !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
248               zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) )
249               zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) )
250               !                                        ! wslpi and wslpj output in zwz and zww, resp.
251               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
252               zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
253               zwz(ji,jj,jk) = (     zai / ( zbi - zeps)  * ( 1. - zalpha )   &
254                  &             + zcoef3 * wslpiml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
255               zww(ji,jj,jk) = (     zaj / ( zbj - zeps)  * ( 1. - zalpha )   &
256                  &             + zcoef3 * wslpjml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
257            END DO
258         END DO
259      END DO
260      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
261      !
262      !                                           !* horizontal Shapiro filter
263      DO jk = 2, jpkm1
264         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
265            DO ji = 2, jpim1
266               zcofw = tmask(ji,jj,jk) / 16.
267               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
268                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
269                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
270                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
271                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
272
273               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
274                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
275                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
276                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
277                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
278            END DO
279         END DO 
280         DO jj = 3, jpj-2                               ! other rows
281            DO ji = fs_2, fs_jpim1   ! vector opt.
282               zcofw = tmask(ji,jj,jk) / 16.
283               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
284                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
285                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
286                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
287                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
288
289               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
290                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
291                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
292                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
293                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
294            END DO
295         END DO
296         !                                        !* decrease along coastal boundaries
297         DO jj = 2, jpjm1
298            DO ji = fs_2, fs_jpim1   ! vector opt.
299               z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5
300               z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5
301               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v
302               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v
303            END DO
304         END DO
305      END DO
306         
307
308      ! III.  Specific grid points     
309      ! ===========================
310      !               
311      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
312         !                                                    ! Gibraltar Strait
313         ij0 =  50   ;   ij1 =  53
314         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
315         ij0 =  51   ;   ij1 =  53
316         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
317         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
318         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
319         !
320         !                                                    ! Mediterrannean Sea
321         ij0 =  49   ;   ij1 =  56
322         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
323         ij0 =  50   ;   ij1 =  56
324         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
325         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
326         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
327      ENDIF
328
329     
330      ! IV. Lateral boundary conditions
331      ! ===============================
332      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
333      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
334
335
336      IF(ln_ctl) THEN
337         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
338         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
339      ENDIF
340      !
341   END SUBROUTINE ldf_slp
342
343
344   SUBROUTINE ldf_slp_mxl( prd, pn2 )
345      !!----------------------------------------------------------------------
346      !!                  ***  ROUTINE ldf_slp_mxl  ***
347      !!
348      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
349      !!              the mixed layer.
350      !!
351      !! ** Method :
352      !!      The slope in the i-direction is computed at u- and w-points
353      !!   (uslp, wslpi) and the slope in the j-direction is computed at
354      !!   v- and w-points (vslp, wslpj).
355      !!   They are bounded by 1/100 over the whole ocean, and within the
356      !!   surface layer they are bounded by the distance to the surface
357      !!   ( slope<= depth/l  where l is the length scale of horizontal
358      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
359      !!   of 10cm/s)
360      !!
361      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
362      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
363      !!----------------------------------------------------------------------
364      USE oce , zgru  => ua   ! ua, va used as workspace and set to hor.
365      USE oce , zgrv  => va   ! density gradient in ldf_slp
366      !!
367      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density
368      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
369      !!
370      INTEGER  ::   ji, jj, jk           ! dummy loop indices
371      INTEGER  ::   ik, ikm1             ! temporary integers
372      REAL(wp) ::   zeps, zmg, zm05g     ! temporary scalars
373      REAL(wp) ::   zcoef1, zcoef2       !    -         -
374      REAL(wp) ::   zau, zbu, zai, zbi   !    -         -
375      REAL(wp) ::   zav, zbv, zaj, zbj   !    -         -
376      REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace
377      !!----------------------------------------------------------------------
378
379      zeps  =  1.e-20                    ! Local constant initialization
380      zmg   = -1.0 / grav
381      zm05g = -0.5 / grav
382      !
383      uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
384      vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0
385      wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0
386      wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
387
388      !                                  ! surface mixed layer mask
389      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
390# if defined key_vectopt_loop
391         DO jj = 1, 1
392            DO ji = 1, jpij   ! vector opt. (forced unrolling)
393# else
394         DO jj = 1, jpj
395            DO ji = 1, jpi
396# endif
397               ik = nmln(ji,jj) - 1
398               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0
399               ELSE                  ;   omlmask(ji,jj,jk) = 0.e0
400               ENDIF
401            END DO
402         END DO
403      END DO
404
405
406      ! Slopes of isopycnal surfaces just before bottom of mixed layer
407      ! --------------------------------------------------------------
408      ! The slope are computed as in the 3D case.
409      ! A key point here is the definition of the mixed layer at u- and v-points.
410      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
411      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
412      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
413      ! induce velocity field near the base of the mixed layer.
414      !-----------------------------------------------------------------------
415      !
416      zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2)
417      zwy(jpi,:) = 0.e0
418# if defined key_vectopt_loop
419      DO jj = 1, 1
420         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
421# else
422      DO jj = 1, jpjm1
423         DO ji = 1, jpim1
424# endif
425            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation
426            ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values
427            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
428               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
429         END DO
430      END DO
431      CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change
432
433      !                            !* Slope at u points
434# if defined key_vectopt_loop
435      DO jj = 1, 1
436         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
437# else
438      DO jj = 2, jpjm1
439         DO ji = 2, jpim1
440# endif
441            ! horizontal and vertical density gradient at u-points
442            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
443            ik = MIN( ik, jpkm1 )
444            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
445            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
446            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
447            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
448            zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) )
449            ! uslpml
450            uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
451         END DO
452      END DO
453      CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change)
454
455      !                            !* vertical density gradient for v-slope (from N^2)
456# if defined key_vectopt_loop
457      DO jj = 1, 1
458         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
459# else
460      DO jj = 1, jpjm1
461         DO ji = 1, jpim1
462# endif
463            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
464            ik = MIN( ik, jpkm1 )
465            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
466               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
467         END DO
468      END DO
469      CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change
470
471      !                            !* Slope at v points
472# if defined key_vectopt_loop
473      DO jj = 1, 1
474         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
475# else
476      DO jj = 2, jpjm1
477         DO ji = 2, jpim1
478# endif
479            ! horizontal and vertical density gradient at v-points
480            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
481            ik = MIN( ik,jpkm1 )
482            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
483            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
484            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
485            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
486            zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) )
487            ! vslpml
488            vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
489         END DO
490      END DO
491      CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change)
492
493
494      !                            !* vertical density gradient for w-slope (from N^2)
495# if defined key_vectopt_loop
496      DO jj = 1, 1
497         DO ji = 1, jpij   ! vector opt. (forced unrolling)
498# else
499      DO jj = 1, jpj
500         DO ji = 1, jpi
501# endif
502            ik = nmln(ji,jj) + 1
503            ik = MIN( ik, jpk )
504            ikm1 = MAX ( 1, ik-1)
505            zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
506               &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
507         END DO
508      END DO
509
510      !                            !* Slopes at w points
511# if defined key_vectopt_loop
512      DO jj = 1, 1
513         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
514# else
515      DO jj = 2, jpjm1
516         DO ji = 2, jpim1
517# endif
518            ik = nmln(ji,jj) + 1
519            ik = MIN( ik, jpk )
520            ikm1 = MAX ( 1, ik-1 )
521            ! horizontal density i-gradient at w-points
522            zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   &
523               &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
524            zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
525            zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   &
526               &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik)
527            ! horizontal density j-gradient at w-points
528            zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    &
529               &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) )
530            zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
531            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
532               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
533            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
534            !                   static instability:
535            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
536            zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
537            zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) )
538            ! wslpiml and wslpjml
539            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
540            wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
541         END DO
542      END DO
543      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions
544      !
545   END SUBROUTINE ldf_slp_mxl
546
547
548   SUBROUTINE ldf_slp_init
549      !!----------------------------------------------------------------------
550      !!                  ***  ROUTINE ldf_slp_init  ***
551      !!
552      !! ** Purpose :   Initialization for the isopycnal slopes computation
553      !!
554      !! ** Method  :   read the nammbf namelist and check the parameter
555      !!      values called by tra_dmp at the first timestep (nit000)
556      !!----------------------------------------------------------------------
557      INTEGER ::   ji, jj, jk   ! dummy loop indices
558      !!----------------------------------------------------------------------
559     
560      IF(lwp) THEN   
561         WRITE(numout,*)
562         WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
563         WRITE(numout,*) '~~~~~~~'
564      ENDIF
565
566      ! Direction of lateral diffusion (tracers and/or momentum)
567      ! ------------------------------
568      ! set the slope to zero (even in s-coordinates)
569
570      uslp (:,:,:) = 0.e0
571      vslp (:,:,:) = 0.e0
572      wslpi(:,:,:) = 0.e0
573      wslpj(:,:,:) = 0.e0
574
575      uslpml (:,:) = 0.e0
576      vslpml (:,:) = 0.e0
577      wslpiml(:,:) = 0.e0
578      wslpjml(:,:) = 0.e0
579
580      IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
581         IF(lwp) THEN
582            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
583         ENDIF
584
585         ! geopotential diffusion in s-coordinates on tracers and/or momentum
586         ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
587         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
588
589         ! set the slope of diffusion to the slope of s-surfaces
590         !      ( c a u t i o n : minus sign as fsdep has positive value )
591         DO jk = 1, jpk
592            DO jj = 2, jpjm1
593               DO ji = fs_2, fs_jpim1   ! vector opt.
594                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
595                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
596                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
597                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
598               END DO
599            END DO
600         END DO
601         ! Lateral boundary conditions on the slopes
602         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
603         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
604      ENDIF
605      !
606   END SUBROUTINE ldf_slp_init
607
608#else
609   !!------------------------------------------------------------------------
610   !!   Dummy module :                 NO Rotation of lateral mixing tensor
611   !!------------------------------------------------------------------------
612   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
613CONTAINS
614   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
615      INTEGER, INTENT(in) :: kt 
616      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
617      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
618   END SUBROUTINE ldf_slp
619#endif
620
621   !!======================================================================
622END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.