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 NEMO/branches/UKMO/BPC_miniapp/OpenMP – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/OpenMP/ldfslp.F90 @ 10838

Last change on this file since 10838 was 10838, checked in by wayne_gaudin, 5 years ago

Ticket #2197 - extracted versions added
File size: 28.3 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     1.0  ! 2002-10  (G. Madec)  Free form, F90
10   !!             -   ! 2005-10  (A. Beckmann)  correction for s-coordinates
11   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator
12   !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML
13   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  add limiter on triad slopes
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   ldf_slp       : calculates the slopes of neutral surface   (Madec operator)
18   !!   ldf_slp_mxl   : calculates the slopes at the base of the mixed layer (Madec operator)
19   !!----------------------------------------------------------------------
20   USE len_oce        ! ocean sizes
21   USE phycst         ! physical constants
22   USE in_out_manager ! I/O manager
23! mjb   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   ldf_slp         ! routine called by step.F90
29   PUBLIC   ldf_slp_mxl     
30
31   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp      !: slope limit                                     (nam_traldf namelist)
32   
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (./LICENSE)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE ldf_slp( ln_zps, ln_isfcav, prd, pn2, tmask, umask, vmask, wmask, gru, grv,              &
43   &                   e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n, mbku, mbkv,    &
44   &                   mikt, miku, mikv, nmln, hmlp, hmlpt, risfdep,                                   &
45   &                   omlmask, uslpml, vslpml, wslpiml, wslpjml, uslp, vslp, wslpi, wslpj )
46
47      !!----------------------------------------------------------------------
48      !!                 ***  ROUTINE ldf_slp  ***
49      !!
50      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
51      !!              surfaces referenced locally) (ln_traldf_iso=T).
52      !!
53      !! ** Method  :   The slope in the i-direction is computed at U- and
54      !!      W-points (uslp, wslpi) and the slope in the j-direction is
55      !!      computed at V- and W-points (vslp, wslpj).
56      !!      They are bounded by 1/100 over the whole ocean, and within the
57      !!      surface layer they are bounded by the distance to the surface
58      !!      ( slope<= depth/l  where l is the length scale of horizontal
59      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
60      !!      of 10cm/s)
61      !!        A horizontal shapiro filter is applied to the slopes
62      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
63      !!      the slope of the model level surface.
64      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
65      !!      [slopes already set to zero at level 1, and to zero or the ocean
66      !!      bottom slope (ln_sco=T) at level jpk in inildf]
67      !!
68      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
69      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
70      !!----------------------------------------------------------------------
71      LOGICAL , INTENT(in)                   ::   ln_zps, ln_isfcav
72      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   prd   ! in situ density
73      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
74      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   tmask, umask, vmask, wmask
75      REAL(wp), INTENT(in), DIMENSION(:,:)   ::   gru, grv, e1t, e2t, r1_e1u, r1_e2v
76      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   e3u_n, e3v_n, e3w_n, gdept_n, gdepw_n
77      INTEGER , INTENT(in), DIMENSION(:,:)   ::   mbku, mbkv, mikt, miku, mikv, nmln
78      REAL(wp), INTENT(in), DIMENSION(:,:)   ::   hmlp, hmlpt, risfdep
79      REAL(wp), INTENT(out),DIMENSION(:,:,:) ::   omlmask
80      REAL(wp), INTENT(out),DIMENSION(:,:)   ::   uslpml, vslpml, wslpiml, wslpjml
81      REAL(wp), INTENT(out),DIMENSION(:,:,:) ::   uslp, vslp, wslpi, wslpj
82
83      !!
84      INTEGER  ::   ji , jj , jk    ! dummy loop indices
85      INTEGER  ::   ii0, ii1        ! temporary integer
86      INTEGER  ::   ij0, ij1        ! temporary integer
87      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars
88      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
89      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
90      REAL(wp) ::   zck, zfk,      zbw             !   -      -
91      REAL(wp) ::   zdepu, zdepv                   !   -      -
92      REAL(wp),ALLOCATABLE, DIMENSION(:,:)     ::  zslpml_hmlpu, zslpml_hmlpv !  Make global scratch at some point
93      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)   ::  zgru, zwz, zdzr
94      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)   ::  zgrv, zww
95      !!----------------------------------------------------------------------
96      !
97      ALLOCATE(zslpml_hmlpu(jpi,jpj))
98      ALLOCATE(zslpml_hmlpv(jpi,jpj))
99      ALLOCATE(zgru(jpi,jpj,jpk))
100      ALLOCATE(zwz(jpi,jpj,jpk))
101      ALLOCATE(zdzr(jpi,jpj,jpk))
102      ALLOCATE(zgrv(jpi,jpj,jpk))
103      ALLOCATE(zww(jpi,jpj,jpk))
104! mjb      IF( ln_timing )   CALL timing_start('ldf_slp')
105      !
106      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
107      z1_16  =  1.0_wp / 16._wp
108      zm1_g  = -1.0_wp / grav
109      zm1_2g = -0.5_wp / grav
110      z1_slpmax = 1._wp / rn_slpmax
111      !
112      zww(:,:,:) = 0._wp
113      zwz(:,:,:) = 0._wp
114      !
115      DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
116         DO jj = 1, jpjm1
117            DO ji = 1, fs_jpim1   ! vector opt.
118               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )
119               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )
120            END DO
121         END DO
122      END DO
123      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
124         DO jj = 1, jpjm1
125            DO ji = 1, jpim1
126               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
127               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
128            END DO
129         END DO
130      ENDIF
131
132! MJB next lines commented out for simplicity
133!      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level
134!         DO jj = 1, jpjm1
135!            DO ji = 1, jpim1
136!               IF( miku(ji,jj) > 1 )   zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)
137!               IF( mikv(ji,jj) > 1 )   zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
138!            END DO
139!         END DO
140!      ENDIF
141      !
142      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
143      DO jk = 2, jpkm1
144         !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
145         !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0
146         !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1
147         !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
148         !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
149         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
150            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
151      END DO
152      !
153      !                          !==   Slopes just below the mixed layer   ==!
154      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, nmln, e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n,  & 
155      &                 tmask, umask, vmask, omlmask, uslpml, vslpml, wslpiml, wslpjml )
156
157
158      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
159      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
160      !
161
162      IF ( ln_isfcav ) THEN
163         DO jj = 2, jpjm1
164            DO ji = fs_2, fs_jpim1   ! vector opt.
165               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt  (ji,jj), hmlpt  (ji+1,jj  ), 5._wp) &
166                  &                                  - MAX(risfdep(ji,jj), risfdep(ji+1,jj  )       ) ) 
167               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt  (ji,jj), hmlpt  (ji  ,jj+1), 5._wp) &
168                  &                                  - MAX(risfdep(ji,jj), risfdep(ji  ,jj+1)       ) )
169            END DO
170         END DO
171      ELSE
172         DO jj = 2, jpjm1
173            DO ji = fs_2, fs_jpim1   ! vector opt.
174               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)
175               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)
176            END DO
177         END DO
178      END IF
179
180      DO jk = 2, jpkm1                            !* Slopes at u and v points
181         DO jj = 2, jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               !                                      ! horizontal and vertical density gradient at u- and v-points
184               zau = zgru(ji,jj,jk) * r1_e1u(ji,jj)
185               zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj)
186               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
187               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
188               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
189               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
190               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau )  )
191               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav )  )
192               !                                      ! uslp and vslp output in zwz and zww, resp.
193               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
194               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
195               ! thickness of water column between surface and level k at u/v point
196               zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) )                            &
197                                - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj))   )
198               zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) )                            &
199                                - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj))   )
200               !
201               zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps )                                     &
202                  &                      + zfi  * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk)
203               zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps )                                     &
204                  &                      + zfj  * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk)
205!!gm  modif to suppress omlmask.... (as in Griffies case)
206!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
207!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
208!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
209!               zci = 0.5 * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
210!               zcj = 0.5 * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
211!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
212!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
213!!gm end modif
214            END DO
215         END DO
216      END DO
217!MJB      CALL lbc_lnk_multi( zwz, 'U', -1.,  zww, 'V', -1. )      ! lateral boundary conditions
218      !
219      !                                            !* horizontal Shapiro filter
220      DO jk = 2, jpkm1
221         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
222            DO ji = 2, jpim1
223               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
224                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
225                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
226                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
227                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
228               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
229                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
230                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
231                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
232                  &                       + 4.*  zww(ji,jj    ,jk)                       )
233            END DO
234         END DO
235         DO jj = 3, jpj-2                               ! other rows
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
238                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
239                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
240                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
241                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
242               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
243                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
244                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
245                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
246                  &                       + 4.*  zww(ji,jj    ,jk)                       )
247            END DO
248         END DO
249         !                                        !* decrease along coastal boundaries
250         DO jj = 2, jpjm1
251            DO ji = fs_2, fs_jpim1   ! vector opt.
252               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   &
253                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp
254               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   &
255                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp
256            END DO
257         END DO
258      END DO
259
260
261      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
262      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
263      !
264      DO jk = 2, jpkm1
265         DO jj = 2, jpjm1
266            DO ji = fs_2, fs_jpim1   ! vector opt.
267               !                                  !* Local vertical density gradient evaluated from N^2
268               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
269               !                                  !* Slopes at w point
270               !                                        ! i- & j-gradient of density at w-points
271               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           &
272                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj)
273               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           &
274                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) * e2t(ji,jj)
275               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
276                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * wmask (ji,jj,jk)
277               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
278                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * wmask (ji,jj,jk)
279               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
280               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
281               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai )  )
282               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj )  )
283               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
284               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
285               zck = ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - gdepw_n(ji,jj,mikt(ji,jj)), 10._wp )
286               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
287               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk)
288
289!!gm  modif to suppress omlmask....  (as in Griffies operator)
290!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
291!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
292!               zck = gdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
293!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
294!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
295!!gm end modif
296            END DO
297         END DO
298      END DO
299! MJB      CALL lbc_lnk_multi( zwz, 'T', -1.,  zww, 'T', -1. )      ! lateral boundary conditions
300      !
301      !                                           !* horizontal Shapiro filter
302      DO jk = 2, jpkm1
303         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
304            DO ji = 2, jpim1
305               zcofw = wmask(ji,jj,jk) * z1_16
306               wslpi(ji,jj,jk) = (         zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
307                    &               +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
308                    &               + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
309                    &               +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
310                    &               + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
311
312               wslpj(ji,jj,jk) = (         zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
313                    &               +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
314                    &               + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
315                    &               +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
316                    &               + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
317            END DO
318         END DO
319         DO jj = 3, jpj-2                               ! other rows
320            DO ji = fs_2, fs_jpim1   ! vector opt.
321               zcofw = wmask(ji,jj,jk) * z1_16
322               wslpi(ji,jj,jk) = (         zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
323                    &               +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
324                    &               + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
325                    &               +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
326                    &               + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
327
328               wslpj(ji,jj,jk) = (         zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
329                    &               +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
330                    &               + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
331                    &               +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
332                    &               + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
333            END DO
334         END DO
335         !                                        !* decrease in vicinity of topography
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1   ! vector opt.
338               zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   &
339                  &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
340               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck
341               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
342            END DO
343         END DO
344      END DO
345
346
347      ! IV. Lateral boundary conditions
348      ! ===============================
349! mjb       CALL lbc_lnk_multi( uslp , 'U', -1. , vslp , 'V', -1. , wslpi, 'W', -1., wslpj, 'W', -1. )
350
351! mjb      IF( ln_timing )   CALL timing_stop('ldf_slp')
352      !
353   END SUBROUTINE ldf_slp
354
355   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, nmln, e1t, e2t, r1_e1u, r1_e2v, e3u_n, e3v_n, e3w_n, &
356  &                        tmask, umask, vmask, omlmask, uslpml, vslpml, wslpiml, wslpjml )
357      !!----------------------------------------------------------------------
358      !!                  ***  ROUTINE ldf_slp_mxl  ***
359      !!
360      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
361      !!              the mixed layer.
362      !!
363      !! ** Method  :   The slope in the i-direction is computed at u- & w-points
364      !!              (uslpml, wslpiml) and the slope in the j-direction is computed
365      !!              at v- and w-points (vslpml, wslpjml) with the same bounds as
366      !!              in ldf_slp.
367      !!
368      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces
369      !!                vslpml, wslpjml    just below the mixed layer
370      !!                omlmask         :  mixed layer mask
371      !!----------------------------------------------------------------------
372      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   prd            ! in situ density
373      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.)
374      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
375      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
376      INTEGER,  DIMENSION(:,:)  , INTENT(in) ::   nmln
377      REAL(wp), DIMENSION(:,:),   INTENT(in) ::   e1t, e2t, r1_e1u, r1_e2v
378      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   e3u_n, e3v_n, e3w_n, tmask, umask, vmask
379      REAL(wp), DIMENSION(:,:,:), INTENT(out)::   omlmask
380      REAL(wp), DIMENSION(:,:),   INTENT(out)::   uslpml, vslpml, wslpiml, wslpjml
381
382      !!
383      INTEGER  ::   ji , jj , jk                   ! dummy loop indices
384      INTEGER  ::   iku, ikv, ik, ikm1             ! local integers
385      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars
386      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
387      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
388      REAL(wp) ::   zck, zfk,      zbw             !   -      -
389      !!----------------------------------------------------------------------
390      !
391      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
392      zm1_g  = -1.0_wp / grav
393      zm1_2g = -0.5_wp / grav
394      z1_slpmax = 1._wp / rn_slpmax
395      !
396      uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp
397      vslpml (1,:) = 0._wp      ;      vslpml (jpi,:) = 0._wp
398      wslpiml(1,:) = 0._wp      ;      wslpiml(jpi,:) = 0._wp
399      wslpjml(1,:) = 0._wp      ;      wslpjml(jpi,:) = 0._wp
400      !
401      !                                            !==   surface mixed layer mask   !
402      DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise
403         DO jj = 1, jpj
404            DO ji = 1, jpi
405               ik = nmln(ji,jj) - 1
406               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
407               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
408               ENDIF
409            END DO
410         END DO
411      END DO
412
413
414      ! Slopes of isopycnal surfaces just before bottom of mixed layer
415      ! --------------------------------------------------------------
416      ! The slope are computed as in the 3D case.
417      ! A key point here is the definition of the mixed layer at u- and v-points.
418      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
419      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
420      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
421      ! induce velocity field near the base of the mixed layer.
422      !-----------------------------------------------------------------------
423      !
424      DO jj = 2, jpjm1
425         DO ji = 2, jpim1
426            !                        !==   Slope at u- & v-points just below the Mixed Layer   ==!
427            !
428            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
429            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
430            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !
431            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
432            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
433            !                        !- horizontal density gradient at u- & v-points
434            zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj)
435            zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj)
436            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
437            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
438            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau )  )
439            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav )  )
440            !                        !- Slope at u- & v-points (uslpml, vslpml)
441            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
442            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
443            !
444            !                        !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
445            !
446            ik   = MIN( nmln(ji,jj) + 1, jpk )
447            ikm1 = MAX( 1, ik-1 )
448            !                        !- vertical density gradient for w-slope (from N^2)
449            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
450            !                        !- horizontal density i- & j-gradient at w-points
451            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           &
452               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj)
453            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           &
454               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj)
455            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)             &
456               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik)
457            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
458               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask(ji,jj,ik)
459            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
460            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
461            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai )  )
462            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj )  )
463            !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
464            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
465            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
466         END DO
467      END DO
468      !!gm this lbc_lnk should be useless....
469! MJB      CALL lbc_lnk_multi( uslpml , 'U', -1. , vslpml , 'V', -1. , wslpiml, 'W', -1. , wslpjml, 'W', -1. )
470      !
471   END SUBROUTINE ldf_slp_mxl
472
473
474
475   !!======================================================================
476END MODULE ldfslp
Note: See TracBrowser for help on using the repository browser.