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 @ 1514

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