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 branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/LDF – NEMO

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 2083

Last change on this file since 2083 was 2027, checked in by cetlod, 14 years ago

Reorganisation of the initialisation phase, see ticket:695
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 31.9 KB
Line 
1MODULE ldfslp
2   !!======================================================================
3   !!                       ***  MODULE  ldfslp  ***
4   !! Ocean physics: slopes of neutral surfaces
5   !!======================================================================
6   !! History :  OPA  ! 1994-12  (G. Madec, M. Imbard)  Original code
7   !!            8.0  ! 1997-06  (G. Madec)  optimization, lbc
8   !!            8.1  ! 1999-10  (A. Jouzeau)  NEW profile in the mixed layer
9   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
10   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
11   !!----------------------------------------------------------------------
12#if   defined key_ldfslp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
15   !!----------------------------------------------------------------------
16   !!   ldf_slp      : compute the slopes of neutral surface
17   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface
18   !!   ldf_slp_init : initialization of the slopes computation
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE ldftra_oce
23   USE ldfdyn_oce
24   USE phycst          ! physical constants
25   USE zdfmxl          ! mixed layer depth
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE in_out_manager  ! I/O manager
28   USE prtctl          ! Print control
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   ldf_slp         ! routine called by step.F90
34   PUBLIC   ldf_slp_init    ! routine called by opa.F90
35
36   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag
37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points
38   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points
39   
40   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt
41   REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
42   REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
49   !! $Id$
50   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
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 isopycnal
60      !!              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      USE oce , zgru  => ua   ! use ua as workspace
81      USE oce , zgrv  => va   ! use va as workspace
82      USE oce , zwy   => ta   ! use ta as workspace
83      USE oce , zwz   => sa   ! use sa as workspace
84      !!
85      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
86      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   prd   ! in situ density
87      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
88      !!
89      INTEGER  ::   ji , jj , jk    ! dummy loop indices
90      INTEGER  ::   ii0, ii1, iku   ! temporary integer
91      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
92      REAL(wp) ::   zeps, zmg, zm05g, zalpha        ! temporary scalars
93      REAL(wp) ::   zcoef1, zcoef2, zcoef3          !    -         -
94      REAL(wp) ::   zcofu , zcofv , zcofw           !    -         -
95      REAL(wp) ::   zau, zbu, zai, zbi, z1u, z1wu   !    -         -
96      REAL(wp) ::   zav, zbv, zaj, zbj, z1v, z1wv   !
97      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww     ! 3D workspace
98      !!----------------------------------------------------------------------
99     
100      zeps  =  1.e-20                             ! Local constant initialization
101      zmg   = -1.0 / grav
102      zm05g = -0.5 / grav
103      !
104      zww(:,:,:) = 0.e0
105      zwz(:,:,:) = 0.e0
106      !                                           ! horizontal density gradient computation
107      DO jk = 1, jpk
108         DO jj = 1, jpjm1
109            DO ji = 1, fs_jpim1   ! vector opt.
110               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
111               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
112            END DO
113         END DO
114      END DO
115      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
116# if defined key_vectopt_loop 
117         DO jj = 1, 1
118            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
119# else
120         DO jj = 1, jpjm1
121            DO ji = 1, jpim1
122# endif
123               iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1               ! last ocean level
124               ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
125               zgru(ji,jj,iku) = gru(ji,jj) 
126               zgrv(ji,jj,ikv) = grv(ji,jj)               
127            END DO
128         END DO
129      ENDIF
130     
131      CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer
132
133     
134      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
135      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
136      !               
137      !                                           !* Local vertical density gradient evaluated from N^2
138      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
139         DO jj = 1, jpj
140            DO ji = 1, jpi
141               zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   &
142                  &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. )
143            END DO
144         END DO
145      END DO
146      DO jk = 2, jpkm1                            !* Slopes at u and v points
147         DO jj = 2, jpjm1
148            DO ji = fs_2, fs_jpim1   ! vector opt.
149               ! horizontal and vertical density gradient at u- and v-points
150               zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk)
151               zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk)
152               zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj  ,jk) )
153               zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji  ,jj+1,jk) )
154               ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
155               !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
156               zbu = MIN( zbu, -100.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) )
157               zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) )
158               ! uslp and vslp output in zwz and zww, resp.
159               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
160               zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           &
161                  &                    + zalpha  * uslpml(ji,jj)                                                  &
162                  &                              * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
163                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk)
164               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
165               zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           &
166                  &                    + zalpha  * vslpml(ji,jj)                                                  &
167                  &                              * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
168                  &                              / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
169            END DO
170         END DO
171      END DO
172      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
173      !
174      zcofu = 1. / 16.                            !* horizontal Shapiro filter
175      zcofv = 1. / 16.
176      DO jk = 2, jpkm1
177         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only
178            DO ji = 2, jpim1 
179               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
180                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
181                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
182                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
183                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
184               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
185                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
186                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
187                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
188                  &                       + 4.*  zww(ji,jj    ,jk)                       )
189            END DO
190         END DO
191         DO jj = 3, jpj-2                               ! other rows
192            DO ji = fs_2, fs_jpim1   ! vector opt.
193               uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
194                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
195                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
196                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
197                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
198               vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
199                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
200                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
201                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
202                  &                       + 4.*  zww(ji,jj    ,jk)                       )
203            END DO
204         END DO
205         !                                        !* decrease along coastal boundaries
206         DO jj = 2, jpjm1
207            DO ji = fs_2, fs_jpim1   ! vector opt.
208               z1u  = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5
209               z1v  = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5
210               z1wu = ( umask(ji,jj,jk)   + umask(ji,jj,jk+1) )*.5
211               z1wv = ( vmask(ji,jj,jk)   + vmask(ji,jj,jk+1) )*.5
212               uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu
213               vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv
214            END DO
215         END DO
216      END DO
217
218
219      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
220      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
221      !               
222      !                                           !* Local vertical density gradient evaluated from N^2
223      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point
224         DO jj = 1, jpj
225            DO ji = 1, jpi
226               zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
227            END DO
228         END DO
229      END DO
230      DO jk = 2, jpkm1                            !* Slopes at w point
231         DO jj = 2, jpjm1
232            DO ji = fs_2, fs_jpim1   ! vector opt.
233               !                                        ! horizontal density i-gradient at w-points
234               zcoef1 = MAX( zeps, umask(ji-1,jj,jk  )+umask(ji,jj,jk  )    &
235                  &               +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) )
236               zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
237               zai = zcoef1 * (  zgru(ji  ,jj,jk  ) + zgru(ji  ,jj,jk-1)   &
238                  &            + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk  ) ) * tmask (ji,jj,jk)
239               !                                        ! horizontal density j-gradient at w-points
240               zcoef2 = MAX( zeps, vmask(ji,jj-1,jk  )+vmask(ji,jj,jk-1)   &
241                  &               +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk  ) )
242               zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
243               zaj = zcoef2 * (  zgrv(ji,jj  ,jk  ) + zgrv(ji,jj  ,jk-1)   &
244                  &            + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk  ) ) * tmask (ji,jj,jk)
245               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
246               !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
247               zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) )
248               zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) )
249               !                                        ! wslpi and wslpj output in zwz and zww, resp.
250               zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )
251               zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
252               zwz(ji,jj,jk) = (     zai / ( zbi - zeps)  * ( 1. - zalpha )   &
253                  &             + zcoef3 * wslpiml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
254               zww(ji,jj,jk) = (     zaj / ( zbj - zeps)  * ( 1. - zalpha )   &
255                  &             + zcoef3 * wslpjml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk)
256            END DO
257         END DO
258      END DO
259      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww
260      !
261      !                                           !* horizontal Shapiro filter
262      DO jk = 2, jpkm1
263         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1
264            DO ji = 2, jpim1
265               zcofw = tmask(ji,jj,jk) / 16.
266               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
267                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
268                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
269                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
270                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
271
272               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
273                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
274                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
275                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
276                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
277            END DO
278         END DO 
279         DO jj = 3, jpj-2                               ! other rows
280            DO ji = fs_2, fs_jpim1   ! vector opt.
281               zcofw = tmask(ji,jj,jk) / 16.
282               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
283                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
284                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
285                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
286                  &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw
287
288               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
289                  &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
290                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
291                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
292                  &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw
293            END DO
294         END DO
295         !                                        !* decrease along coastal boundaries
296         DO jj = 2, jpjm1
297            DO ji = fs_2, fs_jpim1   ! vector opt.
298               z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5
299               z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5
300               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v
301               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v
302            END DO
303         END DO
304      END DO
305         
306
307      ! III.  Specific grid points     
308      ! ===========================
309      !               
310      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area
311         !                                                    ! Gibraltar Strait
312         ij0 =  50   ;   ij1 =  53
313         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
314         ij0 =  51   ;   ij1 =  53
315         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
316         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
317         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
318         !
319         !                                                    ! Mediterrannean Sea
320         ij0 =  49   ;   ij1 =  56
321         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
322         ij0 =  50   ;   ij1 =  56
323         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
324         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
325         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0
326      ENDIF
327
328     
329      ! IV. Lateral boundary conditions
330      ! ===============================
331      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
332      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
333
334
335      IF(ln_ctl) THEN
336         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk)
337         CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
338      ENDIF
339      !
340   END SUBROUTINE ldf_slp
341
342
343   SUBROUTINE ldf_slp_mxl( prd, pn2 )
344      !!----------------------------------------------------------------------
345      !!                  ***  ROUTINE ldf_slp_mxl  ***
346      !!
347      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
348      !!              the mixed layer.
349      !!
350      !! ** Method :
351      !!      The slope in the i-direction is computed at u- and w-points
352      !!   (uslp, wslpi) and the slope in the j-direction is computed at
353      !!   v- and w-points (vslp, wslpj).
354      !!   They are bounded by 1/100 over the whole ocean, and within the
355      !!   surface layer they are bounded by the distance to the surface
356      !!   ( slope<= depth/l  where l is the length scale of horizontal
357      !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
358      !!   of 10cm/s)
359      !!
360      !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
361      !!   of now neutral surfaces at u-, w- and v- w-points, resp.
362      !!----------------------------------------------------------------------
363      USE oce , zgru  => ua   ! ua, va used as workspace and set to hor.
364      USE oce , zgrv  => va   ! density gradient in ldf_slp
365      !!
366      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density
367      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
368      !!
369      INTEGER  ::   ji, jj, jk           ! dummy loop indices
370      INTEGER  ::   ik, ikm1             ! temporary integers
371      REAL(wp) ::   zeps, zmg, zm05g     ! temporary scalars
372      REAL(wp) ::   zcoef1, zcoef2       !    -         -
373      REAL(wp) ::   zau, zbu, zai, zbi   !    -         -
374      REAL(wp) ::   zav, zbv, zaj, zbj   !    -         -
375      REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace
376      !!----------------------------------------------------------------------
377
378      zeps  =  1.e-20                    ! Local constant initialization
379      zmg   = -1.0 / grav
380      zm05g = -0.5 / grav
381      !
382      uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0
383      vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0
384      wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0
385      wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0
386
387      !                                  ! surface mixed layer mask
388      DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise
389# if defined key_vectopt_loop
390         DO jj = 1, 1
391            DO ji = 1, jpij   ! vector opt. (forced unrolling)
392# else
393         DO jj = 1, jpj
394            DO ji = 1, jpi
395# endif
396               ik = nmln(ji,jj) - 1
397               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0
398               ELSE                  ;   omlmask(ji,jj,jk) = 0.e0
399               ENDIF
400            END DO
401         END DO
402      END DO
403
404
405      ! Slopes of isopycnal surfaces just before bottom of mixed layer
406      ! --------------------------------------------------------------
407      ! The slope are computed as in the 3D case.
408      ! A key point here is the definition of the mixed layer at u- and v-points.
409      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
410      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
411      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
412      ! induce velocity field near the base of the mixed layer.
413      !-----------------------------------------------------------------------
414      !
415      zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2)
416      zwy(jpi,:) = 0.e0
417# if defined key_vectopt_loop
418      DO jj = 1, 1
419         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
420# else
421      DO jj = 1, jpjm1
422         DO ji = 1, jpim1
423# endif
424            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation
425            ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values
426            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
427               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
428         END DO
429      END DO
430      CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change
431
432      !                            !* Slope at u points
433# if defined key_vectopt_loop
434      DO jj = 1, 1
435         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
436# else
437      DO jj = 2, jpjm1
438         DO ji = 2, jpim1
439# endif
440            ! horizontal and vertical density gradient at u-points
441            ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )
442            ik = MIN( ik, jpkm1 )
443            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik)
444            zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) )
445            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
446            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
447            zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) )
448            ! uslpml
449            uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik)
450         END DO
451      END DO
452      CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change)
453
454      !                            !* vertical density gradient for v-slope (from N^2)
455# if defined key_vectopt_loop
456      DO jj = 1, 1
457         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
458# else
459      DO jj = 1, jpjm1
460         DO ji = 1, jpim1
461# endif
462            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
463            ik = MIN( ik, jpkm1 )
464            zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   &
465               &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. )
466         END DO
467      END DO
468      CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change
469
470      !                            !* Slope at v points
471# if defined key_vectopt_loop
472      DO jj = 1, 1
473         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
474# else
475      DO jj = 2, jpjm1
476         DO ji = 2, jpim1
477# endif
478            ! horizontal and vertical density gradient at v-points
479            ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) )
480            ik = MIN( ik,jpkm1 )
481            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik)
482            zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) )
483            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
484            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
485            zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) )
486            ! vslpml
487            vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik)
488         END DO
489      END DO
490      CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change)
491
492
493      !                            !* vertical density gradient for w-slope (from N^2)
494# if defined key_vectopt_loop
495      DO jj = 1, 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            ik = nmln(ji,jj) + 1
502            ik = MIN( ik, jpk )
503            ikm1 = MAX ( 1, ik-1)
504            zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     &
505               &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
506         END DO
507      END DO
508
509      !                            !* Slopes at w points
510# if defined key_vectopt_loop
511      DO jj = 1, 1
512         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
513# else
514      DO jj = 2, jpjm1
515         DO ji = 2, jpim1
516# endif
517            ik = nmln(ji,jj) + 1
518            ik = MIN( ik, jpk )
519            ikm1 = MAX ( 1, ik-1 )
520            ! horizontal density i-gradient at w-points
521            zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   &
522               &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) )
523            zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) )
524            zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   &
525               &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik)
526            ! horizontal density j-gradient at w-points
527            zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    &
528               &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) )
529            zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) )
530            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   &
531               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik)
532            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
533            !                   static instability:
534            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
535            zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) )
536            zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) )
537            ! wslpiml and wslpjml
538            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik)
539            wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik)
540         END DO
541      END DO
542      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions
543      !
544   END SUBROUTINE ldf_slp_mxl
545
546
547   SUBROUTINE ldf_slp_init
548      !!----------------------------------------------------------------------
549      !!                  ***  ROUTINE ldf_slp_init  ***
550      !!
551      !! ** Purpose :   Initialization for the isopycnal slopes computation
552      !!
553      !! ** Method  :   read the nammbf namelist and check the parameter
554      !!      values called by tra_dmp at the first timestep (nit000)
555      !!----------------------------------------------------------------------
556      INTEGER ::   ji, jj, jk   ! dummy loop indices
557      !!----------------------------------------------------------------------
558     
559      IF(lwp) THEN   
560         WRITE(numout,*)
561         WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
562         WRITE(numout,*) '~~~~~~~'
563      ENDIF
564
565      ! Direction of lateral diffusion (tracers and/or momentum)
566      ! ------------------------------
567      ! set the slope to zero (even in s-coordinates)
568
569      uslp (:,:,:) = 0.e0
570      vslp (:,:,:) = 0.e0
571      wslpi(:,:,:) = 0.e0
572      wslpj(:,:,:) = 0.e0
573
574      uslpml (:,:) = 0.e0
575      vslpml (:,:) = 0.e0
576      wslpiml(:,:) = 0.e0
577      wslpjml(:,:) = 0.e0
578
579      IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
580         IF(lwp) THEN
581            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
582         ENDIF
583
584         ! geopotential diffusion in s-coordinates on tracers and/or momentum
585         ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
586         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
587
588         ! set the slope of diffusion to the slope of s-surfaces
589         !      ( c a u t i o n : minus sign as fsdep has positive value )
590         DO jk = 1, jpk
591            DO jj = 2, jpjm1
592               DO ji = fs_2, fs_jpim1   ! vector opt.
593                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
594                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
595                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
596                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
597               END DO
598            END DO
599         END DO
600         ! Lateral boundary conditions on the slopes
601         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
602         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
603      ENDIF
604      !
605   END SUBROUTINE ldf_slp_init
606
607#else
608   !!------------------------------------------------------------------------
609   !!   Dummy module :                 NO Rotation of lateral mixing tensor
610   !!------------------------------------------------------------------------
611   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp = .FALSE.    !: slopes flag
612CONTAINS
613   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine
614      INTEGER, INTENT(in) :: kt 
615      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
616      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
617   END SUBROUTINE ldf_slp
618#endif
619
620   !!======================================================================
621END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.