New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
ldfslp.F90 in trunk/NEMO/OFF_SRC/LDF – NEMO

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

Last change on this file since 1156 was 1156, checked in by rblod, 16 years ago

Update Id and licence information, see ticket #210

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