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

source: trunk/NEMO/OFF_SRC/LDF/ldfslp.F90 @ 343

Last change on this file since 343 was 343, checked in by opalod, 18 years ago

nemo_v1_update_O29:RB: add header for OFFLINE component

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