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

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

CT : UPDATE023 : Addition of new diagnostics controled with logical key l_ctl
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.7 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      IF(l_ctl) THEN
481         WRITE(numout,*) ' slp  - u : ', SUM( uslp (:,:,:) ), ' v : ', SUM( vslp (:,:,:) )
482         WRITE(numout,*) '      - wi: ', SUM( wslpi(:,:,:) ), ' wj: ', SUM( wslpj(:,:,:) )
483      ENDIF
484
485   END SUBROUTINE ldf_slp
486
487
488   SUBROUTINE ldf_slp_mxl( prd, pn2 )
489      !!----------------------------------------------------------------------
490      !!                  ***  ROUTINE ldf_slp_mxl  ***
491      !! ** Purpose :
492      !!     Compute the slopes of iso-neutral surface (slope of isopycnal
493      !!   surfaces referenced locally) just above the mixed layer.
494      !!
495      !! ** Method :
496      !!      The slope in the i-direction is computed at u- and w-points
497      !!   (uslp, wslpi) and the slope in the j-direction is computed at
498      !!   v- and w-points (vslp, wslpj).
499      !!   They are bounded by 1/100 over the whole ocean, and within the
500      !!   surface layer they are bounded by the distance to the surface
501      !!   ( slope<= depth/l  where l is the length scale of horizontal
502      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
503      !!   of 10cm/s)
504      !!
505      !! ** Action :
506      !!      Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
507      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
508      !!
509      !! History :
510      !!   8.1  !  99-10  (A. Jouzeau)  Original code
511      !!   8.5  !  99-10  (G. Madec)  Free form, F90
512      !!----------------------------------------------------------------------
513      !! * Modules used
514      USE oce           , zgru  => ua,  &  ! ua, va used as workspace and set to hor.
515                          zgrv  => va      ! density gradient in ldf_slp
516
517      !! * Arguments
518      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
519         prd,                           &  ! in situ density
520         pn2                              ! Brunt-Vaisala frequency (locally ref.)
521
522      !! * Local declarations
523      INTEGER  ::   ji, jj, jk             ! dummy loop indices
524      INTEGER  ::   ik, ikm1               ! temporary integers
525      REAL(wp), DIMENSION(jpi,jpj) ::   &
526         zwy                               ! temporary workspace
527      REAL(wp) ::   &
528         zeps, zmg, zm05g,              &  ! temporary scalars
529         zcoef1, zcoef2,                &  !    "         "
530         zau, zbu, zav, zbv,            &  !    "         "
531         zai, zbi, zaj, zbj                !    "         "
532      !!----------------------------------------------------------------------
533
534
535      ! 0. Local constant initialization
536      ! --------------------------------
537
538      zeps  =  1.e-20
539      zmg   = -1.0 / grav
540      zm05g = -0.5 / grav
541
542
543      uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
544      vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0
545      wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0
546      wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
547
548      ! surface mixed layer mask
549
550      ! mask for mixed layer
551      DO jk = 1, jpk
552# if defined key_vectopt_loop   &&   ! defined key_autotasking
553         jj = 1
554         DO ji = 1, jpij   ! vector opt. (forced unrolling)
555# else
556         DO jj = 1, jpj
557            DO ji = 1, jpi
558# endif
559               ! mixed layer interior (mask = 1) and exterior (mask = 0)
560               ik = nmln(ji,jj) - 1
561               IF( jk <= ik ) THEN
562                  omlmask(ji,jj,jk) = 1.e0
563               ELSE
564                  omlmask(ji,jj,jk) = 0.e0
565               ENDIF
566# if ! defined key_vectopt_loop   ||   defined key_autotasking
567            END DO
568# endif
569         END DO
570      END DO
571
572
573      ! Slopes of isopycnal surfaces just before bottom of mixed layer
574      ! --------------------------------------------------------------
575      ! uslpml = d/di( prd ) / d/dz( prd )
576      ! vslpml = d/dj( prd ) / d/dz( prd )
577
578      ! Local vertical density gradient evaluated from N^2
579      ! zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
580
581      !-----------------------------------------------------------------------
582      zwy(:,jpj) = 0.e0
583      zwy(jpi,:) = 0.e0
584# if defined key_vectopt_loop   &&   ! defined key_autotasking
585      jj = 1
586      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
587# else
588      DO jj = 1, jpjm1
589         DO ji = 1, jpim1
590# endif
591            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
592            ! if ik = jpk take jpkm1 values
593            ik = MIN( ik,jpkm1 )
594            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. )   &
595               &             * ( pn2(ji,jj,ik) + pn2(ji,jj,ik+1) )   &
596               &             / MAX( tmask(ji,jj,ik) + tmask (ji,jj,ik+1), 1. )
597# if ! defined key_vectopt_loop   ||   defined key_autotasking
598         END DO
599# endif
600      END DO
601      ! lateral boundary conditions on zwy
602      CALL lbc_lnk( zwy, 'U', -1. )
603
604      ! Slope at u points
605# if defined key_vectopt_loop   &&   ! defined key_autotasking
606      jj = 1
607      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
608# else
609      DO jj = 2, jpjm1
610         DO ji = 2, jpim1
611# endif
612            ! horizontal and vertical density gradient at u-points
613            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
614            ik = MIN( ik,jpkm1 )
615            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
616            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
617            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
618            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
619            zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) )
620            ! uslpml
621            uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
622# if ! defined key_vectopt_loop   ||   defined key_autotasking
623         END DO
624# endif
625      END DO
626
627      ! lateral boundary conditions on uslpml
628      CALL lbc_lnk( uslpml, 'U', -1. )
629
630      ! Local vertical density gradient evaluated from N^2
631      !     zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
632      zwy ( :, jpj) = 0.e0
633      zwy ( jpi, :) = 0.e0
634# if defined key_vectopt_loop   &&   ! defined key_autotasking
635      jj = 1
636      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
637# else
638      DO jj = 1, jpjm1
639         DO ji = 1, jpim1
640# endif
641            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
642            ik = MIN( ik,jpkm1 )
643            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. )   &
644               &             * ( pn2(ji,jj,ik) + pn2(ji,jj,ik+1) )   &
645               &             / MAX( tmask(ji,jj,ik) + tmask (ji,jj,ik+1), 1. )
646# if ! defined key_vectopt_loop   ||   defined key_autotasking
647         END DO
648# endif
649      END DO
650
651      ! lateral boundary conditions on zwy
652      CALL lbc_lnk( zwy, 'V', -1. )
653
654      ! Slope at v points
655# if defined key_vectopt_loop   &&   ! defined key_autotasking
656      jj = 1
657      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
658# else
659      DO jj = 2, jpjm1
660         DO ji = 2, jpim1
661# endif
662            ! horizontal and vertical density gradient at v-points
663            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
664            ik = MIN( ik,jpkm1 )
665            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
666            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
667            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
668            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
669            zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) )
670            ! vslpml
671            vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
672# if ! defined key_vectopt_loop   ||   defined key_autotasking
673         END DO
674# endif
675      END DO
676
677      ! lateral boundary conditions on vslpml
678      CALL lbc_lnk( vslpml, 'V', -1. )
679
680      ! wslpiml = mij( d/di( prd ) / d/dz( prd )
681      ! wslpjml = mij( d/dj( prd ) / d/dz( prd )
682
683
684      ! Local vertical density gradient evaluated from N^2
685      ! zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
686# if defined key_vectopt_loop   &&   ! defined key_autotasking
687      jj = 1
688      DO ji = 1, jpij   ! vector opt. (forced unrolling)
689# else
690      DO jj = 1, jpj
691         DO ji = 1, jpi
692# endif
693            ik = nmln(ji,jj)+1
694            ik = MIN( ik,jpk )
695            ikm1 = MAX ( 1, ik-1)
696            zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
697               &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
698# if ! defined key_vectopt_loop   ||   defined key_autotasking
699         END DO
700# endif
701      END DO
702
703      ! Slope at w point
704# if defined key_vectopt_loop   &&   ! defined key_autotasking
705      jj = 1
706      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
707# else
708      DO jj = 2, jpjm1
709         DO ji = 2, jpim1
710# endif
711            ik = nmln(ji,jj)+1
712            ik = MIN( ik,jpk )
713            ikm1 = MAX ( 1, ik-1)
714            ! horizontal density i-gradient at w-points
715            zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   &
716               &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
717            zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
718            zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   &
719               &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik)
720            ! horizontal density j-gradient at w-points
721            zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    &
722               &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) )
723            zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
724            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
725               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
726            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
727            !                   static instability:
728            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
729            zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
730            zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) )
731            ! wslpiml and wslpjml
732            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
733            wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
734# if ! defined key_vectopt_loop   ||   defined key_autotasking
735         END DO
736# endif
737      END DO
738
739      ! lateral boundary conditions on wslpiml and wslpjml
740      CALL lbc_lnk( wslpiml, 'W', -1. )
741      CALL lbc_lnk( wslpjml, 'W', -1. )
742
743   END SUBROUTINE ldf_slp_mxl
744
745
746   SUBROUTINE ldf_slp_init
747      !!----------------------------------------------------------------------
748      !!                  ***  ROUTINE ldf_slp_init  ***
749      !!
750      !! ** Purpose :   Initialization for the isopycnal slopes computation
751      !!
752      !! ** Method  :   read the nammbf namelist and check the parameter
753      !!      values called by tra_dmp at the first timestep (nit000)
754      !!
755      !! History :
756      !!  8.5  ! 02-06 (G. Madec) original code
757      !!----------------------------------------------------------------------
758      !! * local declarations
759      INTEGER ::   ji, jj, jk   ! dummy loop indices
760      !!----------------------------------------------------------------------
761     
762     
763      ! Parameter control and print
764      ! ---------------------------
765      IF(lwp) THEN
766         WRITE(numout,*)
767         WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
768         WRITE(numout,*) '~~~~~~~'
769      ENDIF
770
771      ! Direction of lateral diffusion (tracers and/or momentum)
772      ! ------------------------------
773      ! set the slope to zero (even in s-coordinates)
774
775      uslp (:,:,:) = 0.e0
776      vslp (:,:,:) = 0.e0
777      wslpi(:,:,:) = 0.e0
778      wslpj(:,:,:) = 0.e0
779
780      uslpml (:,:) = 0.e0
781      vslpml (:,:) = 0.e0
782      wslpiml(:,:) = 0.e0
783      wslpjml(:,:) = 0.e0
784
785      IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN
786
787         ! geopotential diffusion in s-coordinates on tracers and/or momentum
788         ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
789         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
790
791         ! set the slope of diffusion to the slope of s-surfaces
792         !      ( c a u t i o n : minus sign as fsdep has positive value )
793         DO jk = 1, jpk
794            DO jj = 2, jpjm1
795               DO ji = fs_2, fs_jpim1   ! vector opt.
796                  uslp (ji,jj,jk) = -1. / e1u(ji,jj) * umask(ji,jj,jk)   &
797                     &                               * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) )
798                  vslp (ji,jj,jk) = -1. / e2v(ji,jj) * vmask(ji,jj,jk)   &
799                     &                               * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) )
800                  wslpi(ji,jj,jk) = -1. / e1t(ji,jj) * tmask(ji,jj,jk)   &
801                     &                               * ( fsdepuw(ji+1,jj,jk) - fsdepuw(ji,jj,jk) )
802                  wslpj(ji,jj,jk) = -1. / e2t(ji,jj) * tmask(ji,jj,jk)   &
803                     &                               * ( fsdepvw(ji,jj+1,jk) - fsdepvw(ji,jj,jk) )
804               END DO
805            END DO
806         END DO
807
808         ! Lateral boundary conditions on the slopes
809         CALL lbc_lnk( uslp , 'U', -1. )
810         CALL lbc_lnk( vslp , 'V', -1. )
811         CALL lbc_lnk( wslpi, 'W', -1. )
812         CALL lbc_lnk( wslpj, 'W', -1. )
813      ENDIF
814
815   END SUBROUTINE ldf_slp_init
816
817#else
818   !!------------------------------------------------------------------------
819   !!   Dummy module :                 NO Rotation of lateral mixing tensor
820   !!------------------------------------------------------------------------
821   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
822CONTAINS
823   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
824      INTEGER, INTENT(in) :: kt 
825      REAL,DIMENSION(:,:,:), INTENT(in) :: prd, pn2
826      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
827   END SUBROUTINE ldf_slp
828#endif
829
830   !!======================================================================
831END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.