source: trunk/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 13 years ago

get back to the nemo_v2_3 version for trunk

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