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