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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information

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