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

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

CT : UPDATE001 : First major NEMO update

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