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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

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