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 tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/LDF – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/LDF/ldfslp.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

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