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

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

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