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_crs.F90 in branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp_crs.F90 @ 5601

Last change on this file since 5601 was 5601, checked in by cbricaud, 9 years ago

commit changes/bugfix/... for crs ; ok with time-splitting/fixed volume

  • Property svn:executable set to *
File size: 34.7 KB
Line 
1MODULE ldfslp_crs
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   !!----------------------------------------------------------------------
14#if   defined key_ldfslp   ||   defined key_esopa
15   !!----------------------------------------------------------------------
16   !!   'key_ldfslp'                      Rotation of lateral mixing tensor
17   !!----------------------------------------------------------------------
18   !!   ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)
19   !!   ldf_slp      : calculates the slopes of neutral surface   (Madec operator)
20   !!   ldf_slp_mxl  : calculates the slopes at the base of the mixed layer (Madec operator)
21   !!   ldf_slp_init : initialization of the slopes computation
22   !!----------------------------------------------------------------------
23   !USE oce             ! ocean dynamics and tracers
24   !USE dom_oce         ! ocean space and time domain
25   USE ldftra_oce      ! lateral diffusion: traceur
26   USE ldfdyn_oce      ! lateral diffusion: dynamics
27   USE phycst          ! physical constants
28   USE zdfmxl_crs      ! mixed layer depth
29   USE eosbn2_crs      ! equation of states
30   USE crslbclnk       ! ocean lateral boundary conditions (or mpp link)
31   USE in_out_manager  ! I/O manager
32   USE prtctl          ! Print control
33   USE wrk_nemo        ! work arrays
34   USE timing          ! Timing
35   USE crs
36   USE iom
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   ldf_slp_crs        ! routine called by step.F90
42   PUBLIC   ldf_slp_grif_crs   ! routine called by step.F90
43   PUBLIC   ldf_slp_init_crs   ! routine called by opa.F90
44
45 !  LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp_crs = .TRUE.     !: slopes flag
46   !                                                                             !! Madec operator
47   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
48   !                                                                !! Griffies operator
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wslp2                !: wslp**2 from Griffies quarter cells
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate
52
53   !                                                              !! Madec operator
54   !  Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
55   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   omlmask           ! mask of the surface mixed layer at T-pt
56   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
58
59   REAL(wp) ::   repsln = 1.e-25_wp       ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "ldftra_substitute.h90"
64#  include "ldfeiv_substitute.h90"
65#  include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
68   !! $Id: ldfslp.F90 3294 2012-01-28 16:44:18Z rblod $
69   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE ldf_slp_crs( kt, prd, pn2 )
74      !!----------------------------------------------------------------------
75      !!                 ***  ROUTINE ldf_slp  ***
76      !!
77      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
78      !!              surfaces referenced locally) (ln_traldf_iso=T).
79      !!
80      !! ** Method  :   The slope in the i-direction is computed at U- and
81      !!      W-points (uslp, wslpi) and the slope in the j-direction is
82      !!      computed at V- and W-points (vslp, wslpj).
83      !!      They are bounded by 1/100 over the whole ocean, and within the
84      !!      surface layer they are bounded by the distance to the surface
85      !!      ( slope<= depth/l  where l is the length scale of horizontal
86      !!      diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
87      !!      of 10cm/s)
88      !!        A horizontal shapiro filter is applied to the slopes
89      !!        ln_sco=T, s-coordinate, add to the previously computed slopes
90      !!      the slope of the model level surface.
91      !!        macro-tasked on horizontal slab (jk-loop)  (2, jpk-1)
92      !!      [slopes already set to zero at level 1, and to zero or the ocean
93      !!      bottom slope (ln_sco=T) at level jpk in inildf]
94      !!
95      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes
96      !!               of now neutral surfaces at u-, w- and v- w-points, resp.
97      !!----------------------------------------------------------------------
98      INTEGER , INTENT(in)                   ::   kt    ! ocean time-step index
99      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   prd   ! in situ density
100      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   pn2   ! Brunt-Vaisala frequency (locally ref.)
101      !!
102      INTEGER  ::   ji , jj , jk    ! dummy loop indices
103      INTEGER  ::   ii0, ii1, iku   ! temporary integer
104      INTEGER  ::   ij0, ij1, ikv   ! temporary integer
105      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars
106      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
107      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
108      REAL(wp) ::   zck, zfk,      zbw             !   -      -
109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr
111      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zgru, zgrv
112      !!----------------------------------------------------------------------
113      !
114      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_crs')
115      !
116      CALL wrk_alloc( jpi_crs,jpj_crs,jpk, zwz, zww, zdzr, zgru, zgrv )
117
118      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
119      z1_16  =  1.0_wp / 16._wp
120      zm1_g  = -1.0_wp / grav
121      zm1_2g = -0.5_wp / grav
122      !
123      zww(:,:,:) = 0._wp
124      zwz(:,:,:) = 0._wp
125      !
126      DO jk = 1, jpk             !==   i- & j-gradient of density   ==!
127         DO jj = 1, jpj_crsm1
128            DO ji = 1, jpi_crsm1   ! vector opt.
129               zgru(ji,jj,jk) = umask_crs(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )
130               zgrv(ji,jj,jk) = vmask_crs(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )
131            END DO
132         END DO
133      END DO
134      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level
135         DO jj = 1, jpj_crsm1
136            DO ji = 1, jpi_crsm1
137               zgru(ji,jj,mbku_crs(ji,jj)) = gru_crs(ji,jj)
138               zgrv(ji,jj,mbkv_crs(ji,jj)) = grv_crs(ji,jj)
139            END DO
140         END DO
141      ENDIF
142      !
143      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2)
144      DO jk = 2, jpkm1
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         !                               ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
148         !                               !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2
149         !                               ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster
150         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              &
151            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask_crs(:,:,jk+1) )
152      END DO
153      !
154      !                          !==   Slopes just below the mixed layer   ==!
155      CALL ldf_slp_mxl_crs( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml
156
157      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd )
158      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd )
159      !
160      DO jk = 2, jpkm1                            !* Slopes at u and v points
161         DO jj = 2, jpj_crsm1
162            DO ji = 2, jpi_crsm1   ! vector opt.
163               !                                      ! horizontal and vertical density gradient at u- and v-points
164               zau = zgru(ji,jj,jk) / e1u_crs(ji,jj)
165               zav = zgrv(ji,jj,jk) / e2v_crs(ji,jj)
166               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) )
167               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) )
168               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
169               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
170               zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/e3u_max_crs(ji,jj,jk)* ABS( zau )  )
171               zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/e3v_max_crs(ji,jj,jk)* ABS( zav )  )
172               !cc zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/e3u_crs(ji,jj,jk)* ABS( zau )  )
173               !cc zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/e3v_crs(ji,jj,jk)* ABS( zav )  )
174               !                                      ! uslp and vslp output in zwz and zww, resp.
175               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
176               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
177               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
178                  &                   + zfi  * uslpml(ji,jj)                                                     &
179                  &                          * 0.5_wp * ( gdept_crs(ji+1,jj,jk)+gdept_crs(ji,jj,jk) - e3u_max_crs(ji,jj,1) )   &
180                  &                          / MAX( hmlpt_crs(ji,jj), hmlpt_crs(ji+1,jj), 5._wp ) ) * umask_crs(ji,jj,jk)
181               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
182                  &                   + zfj  * vslpml(ji,jj)                                                     &
183                  &                          * 0.5_wp * ( gdept_crs(ji,jj+1,jk)+ gdept_crs(ji,jj,jk)-e3v_max_crs(ji,jj,1) )   &
184                  &                          / MAX( hmlpt_crs(ji,jj), hmlpt_crs(ji,jj+1), 5. ) ) * vmask_crs(ji,jj,jk)
185!!gm  modif to suppress omlmask.... (as in Griffies case)
186!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0.
187!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
188!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
189!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
190!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
191!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
192!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
193!!gm end modif
194            END DO
195         END DO
196      END DO
197      CALL crs_lbc_lnk( zwz, 'U', -1. )   ;   CALL crs_lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions
198      CALL iom_put("zwz_crs",zwz)
199      CALL iom_put("zww_crs",zww)
200      !
201      !                                            !* horizontal Shapiro filter
202      DO jk = 2, jpkm1
203         DO jj = 2, jpj_crsm1, MAX(1, jpj_crs-3)                        ! rows jj=2 and =jpjm1 only
204            DO ji = 2, jpi_crsm1
205               uslp_crs(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
206                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
207                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
208                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
209                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
210               vslp_crs(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
211                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
212                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
213                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
214                  &                       + 4.*  zww(ji,jj    ,jk)                       )
215            END DO
216         END DO
217         DO jj = 3, jpj_crs-2                               ! other rows
218            DO ji = 2, jpi_crsm1   ! vector opt.
219               uslp_crs(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      &
220                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      &
221                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      &
222                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    &
223                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       )
224               vslp_crs(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      &
225                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      &
226                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      &
227                  &                       +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )    &
228                  &                       + 4.*  zww(ji,jj    ,jk)                       )
229            END DO
230         END DO
231         !                                        !* decrease along coastal boundaries
232         DO jj = 2, jpj_crsm1
233            DO ji = 2, jpi_crsm1   ! vector opt.
234               uslp_crs(ji,jj,jk) = uslp_crs(ji,jj,jk) * ( umask_crs(ji,jj+1,jk) + umask_crs(ji,jj-1,jk  ) ) * 0.5_wp   &
235                  &                            * ( umask_crs(ji,jj  ,jk) + umask_crs(ji,jj  ,jk+1) ) * 0.5_wp
236               vslp_crs(ji,jj,jk) = vslp_crs(ji,jj,jk) * ( vmask_crs(ji+1,jj,jk) + vmask_crs(ji-1,jj,jk  ) ) * 0.5_wp   &
237                  &                            * ( vmask_crs(ji  ,jj,jk) + vmask_crs(ji  ,jj,jk+1) ) * 0.5_wp
238            END DO
239         END DO
240      END DO
241
242
243      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd )
244      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd )
245      !
246      DO jk = 2, jpkm1
247         DO jj = 2, jpj_crsm1
248            DO ji = 2, jpi_crsm1   ! vector opt.
249               !                                  !* Local vertical density gradient evaluated from N^2
250               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
251               !                                  !* Slopes at w point
252               !                                        ! i- & j-gradient of density at w-points
253               zci = MAX(  umask_crs(ji-1,jj,jk  ) + umask_crs(ji,jj,jk  )           &
254                  &      + umask_crs(ji-1,jj,jk-1) + umask_crs(ji,jj,jk-1) , zeps  ) * e1t_crs(ji,jj)
255               zcj = MAX(  vmask_crs(ji,jj-1,jk  ) + vmask_crs(ji,jj,jk-1)           &
256                  &      + vmask_crs(ji,jj-1,jk-1) + vmask_crs(ji,jj,jk  ) , zeps  ) *  e2t_crs(ji,jj)
257               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           &
258                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask_crs (ji,jj,jk)
259               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           &
260                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask_crs (ji,jj,jk)
261
262               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
263               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
264               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_max_crs(ji,jj,jk)* ABS( zai )  )
265               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_max_crs(ji,jj,jk)* ABS( zaj )  )
266               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
267               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
268               zck = gdepw_crs(ji,jj,jk) / MAX( hmlp_crs(ji,jj), 10._wp )
269               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask_crs(ji,jj,jk)
270               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask_crs(ji,jj,jk)
271!!gm  modif to suppress omlmask....  (as in Griffies operator)
272!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
273!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
274!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
275!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
276!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
277!!gm end modif
278            END DO
279         END DO
280      END DO
281      CALL crs_lbc_lnk( zwz, 'T', -1. )   ;    CALL crs_lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions
282      !
283      !                                           !* horizontal Shapiro filter
284      DO jk = 2, jpkm1
285         DO jj = 2, jpj_crsm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only
286            DO ji = 2, jpi_crsm1
287               zcofw = tmask_crs(ji,jj,jk) * z1_16
288               wslpi_crs(ji,jj,jk) = (      zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
289                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
290                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
291                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
292                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
293
294               wslpj_crs(ji,jj,jk) = (      zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
295                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
296                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
297                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
298                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
299            END DO
300         END DO
301         DO jj = 3, jpj_crs-2                               ! other rows
302            DO ji = 2, jpi_crsm1   ! vector opt.
303               zcofw = tmask_crs(ji,jj,jk) * z1_16
304               wslpi_crs(ji,jj,jk) = (      zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     &
305                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     &
306                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     &
307                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   &
308                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw
309
310               wslpj_crs(ji,jj,jk) = (      zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     &
311                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     &
312                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     &
313                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   &
314                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw
315            END DO
316         END DO
317         !                                        !* decrease along coastal boundaries
318         DO jj = 2, jpj_crsm1
319            DO ji = 2, jpi_crsm1   ! vector opt.
320               zck =   ( umask_crs(ji,jj,jk) + umask_crs(ji-1,jj,jk) )   &
321                  &  * ( vmask_crs(ji,jj,jk) + vmask_crs(ji,jj-1,jk) ) * 0.25
322               wslpi_crs(ji,jj,jk) = wslpi_crs(ji,jj,jk) * zck
323               wslpj_crs(ji,jj,jk) = wslpj_crs(ji,jj,jk) * zck
324            END DO
325         END DO
326      END DO
327
328      ! IV. Lateral boundary conditions
329      ! ===============================
330      CALL crs_lbc_lnk( uslp_crs , 'U', -1. )     
331      CALL crs_lbc_lnk( vslp_crs , 'V', -1. )
332      CALL crs_lbc_lnk( wslpi_crs, 'W', -1. )      ;      CALL crs_lbc_lnk( wslpj_crs, 'W', -1. )
333      !
334      CALL iom_swap( "nemo_crs" )    ! swap on the coarse grid
335      CALL iom_put("zgru_crs",zgru)
336      CALL iom_put("zgrv_crs",zgrv)
337      CALL iom_put("zdzr_crs",zdzr)
338      CALL iom_put("zwz_crs",zwz)
339      CALL iom_put("zww_crs",zww)
340      CALL iom_put("uslp_crs",uslp_crs)
341      CALL iom_put("vslp_crs",vslp_crs)
342      CALL iom_swap( "nemo" )    ! swap on the coarse grid
343      !
344      CALL wrk_dealloc( jpi_crs,jpj_crs,jpk, zwz, zww, zdzr, zgru, zgrv )
345      !
346      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_crs')
347      !
348   END SUBROUTINE ldf_slp_crs
349
350   SUBROUTINE ldf_slp_mxl_crs( prd, pn2, p_gru, p_grv, p_dzr )
351      !!----------------------------------------------------------------------
352      !!                  ***  ROUTINE ldf_slp_mxl  ***
353      !!
354      !! ** Purpose :   Compute the slopes of iso-neutral surface just below
355      !!              the mixed layer.
356      !!
357      !! ** Method  :   The slope in the i-direction is computed at u- & w-points
358      !!              (uslpml, wslpiml) and the slope in the j-direction is computed
359      !!              at v- and w-points (vslpml, wslpjml) with the same bounds as
360      !!              in ldf_slp.
361      !!
362      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces
363      !!                vslpml, wslpjml    just below the mixed layer
364      !!                omlmask         :  mixed layer mask
365      !!----------------------------------------------------------------------
366      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   prd            ! in situ density
367      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.)
368      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts)
369      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point)
370      !!
371      INTEGER  ::   ji , jj , jk                   ! dummy loop indices
372      INTEGER  ::   iku, ikv, ik, ikm1             ! local integers
373      REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars
374      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      -
375      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      -
376      REAL(wp) ::   zck, zfk,      zbw             !   -      -
377      !!----------------------------------------------------------------------
378      !
379      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_mxl')
380      !
381      zeps   =  1.e-20_wp        !==   Local constant initialization   ==!
382      zm1_g  = -1.0_wp / grav
383      zm1_2g = -0.5_wp / grav
384      !
385      uslpml (1,:) = 0._wp      ;      uslpml (jpi_crs,:) = 0._wp
386      vslpml (1,:) = 0._wp      ;      vslpml (jpi_crs,:) = 0._wp
387      wslpiml(1,:) = 0._wp      ;      wslpiml(jpi_crs,:) = 0._wp
388      wslpjml(1,:) = 0._wp      ;      wslpjml(jpi_crs,:) = 0._wp
389      !
390      !                                            !==   surface mixed layer mask   !
391      DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise
392         DO jj = 1, jpj_crs
393            DO ji = 1, jpi_crs
394               ik = nmln_crs(ji,jj) - 1
395               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp
396               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp
397               ENDIF
398            END DO
399         END DO
400      END DO
401
402
403      ! Slopes of isopycnal surfaces just before bottom of mixed layer
404      ! --------------------------------------------------------------
405      ! The slope are computed as in the 3D case.
406      ! A key point here is the definition of the mixed layer at u- and v-points.
407      ! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
408      ! Otherwise, a n2 value inside the mixed layer can be involved in the computation
409      ! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
410      ! induce velocity field near the base of the mixed layer.
411      !-----------------------------------------------------------------------
412      !
413      DO jj = 2, jpj_crsm1
414         DO ji = 2, jpi_crsm1
415            !                        !==   Slope at u- & v-points just below the Mixed Layer   ==!
416            !
417            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point)
418            iku = MIN(  MAX( 1, nmln_crs(ji,jj) , nmln_crs(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1)
419            ikv = MIN(  MAX( 1, nmln_crs(ji,jj) , nmln_crs(ji,jj+1) ) , jpkm1  )   !
420            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) )
421            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) )
422            !                        !- horizontal density gradient at u- & v-points
423            zau = p_gru(ji,jj,iku) / e1u_crs(ji,jj)
424            zav = p_grv(ji,jj,ikv) / e2v_crs(ji,jj)
425            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
426            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
427            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/e3u_max_crs(ji,jj,iku)* ABS( zau )  )
428            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/e3v_max_crs(ji,jj,ikv)* ABS( zav )  )
429            !                        !- Slope at u- & v-points (uslpml, vslpml)
430            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask_crs(ji,jj,iku)
431            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask_crs(ji,jj,ikv)
432            !
433            !                        !==   i- & j-slopes at w-points just below the Mixed Layer   ==!
434            !
435            ik   = MIN( nmln_crs(ji,jj) + 1, jpk )
436            ikm1 = MAX( 1, ik-1 )
437            !                        !- vertical density gradient for w-slope (from N^2)
438            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
439            !                        !- horizontal density i- & j-gradient at w-points
440            zci = MAX(   umask_crs(ji-1,jj,ik  ) + umask_crs(ji,jj,ik  )           &
441               &       + umask_crs(ji-1,jj,ikm1) + umask_crs(ji,jj,ikm1) , zeps  ) * e1t_crs(ji,jj)
442            zcj = MAX(   vmask_crs(ji,jj-1,ik  ) + vmask_crs(ji,jj,ik  )           &
443               &       + vmask_crs(ji,jj-1,ikm1) + vmask_crs(ji,jj,ikm1) , zeps  ) * e2t_crs(ji,jj)
444            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           &
445               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask_crs(ji,jj,ik)
446            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           &
447               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask_crs(ji,jj,ik)
448            !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
449            !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
450            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_max_crs(ji,jj,ik)* ABS( zai )  )
451            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_max_crs(ji,jj,ik)* ABS( zaj )  )
452            !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
453            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask_crs (ji,jj,ik)
454            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask_crs (ji,jj,ik)
455         END DO
456      END DO
457      !!gm this lbc_lnk should be useless....
458      CALL crs_lbc_lnk( uslpml , 'U', -1. )   ;   CALL crs_lbc_lnk( vslpml , 'V', -1. )   ! lateral boundary cond. (sign change)
459      CALL crs_lbc_lnk( wslpiml, 'W', -1. )   ;   CALL crs_lbc_lnk( wslpjml, 'W', -1. )   ! lateral boundary conditions
460      !
461      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_mxl')
462      !
463   END SUBROUTINE ldf_slp_mxl_crs
464
465
466   SUBROUTINE ldf_slp_init_crs
467      !!----------------------------------------------------------------------
468      !!                  ***  ROUTINE ldf_slp_init  ***
469      !!
470      !! ** Purpose :   Initialization for the isopycnal slopes computation
471      !!
472      !! ** Method  :   read the nammbf namelist and check the parameter
473      !!      values called by tra_dmp at the first timestep (nit000)
474      !!----------------------------------------------------------------------
475      INTEGER ::   ji, jj, jk   ! dummy loop indices
476      INTEGER ::   ierr         ! local integer
477      !!----------------------------------------------------------------------
478      !
479      IF( nn_timing == 1 )  CALL timing_start('ldf_slp_init')
480      !
481      IF(lwp) THEN
482         WRITE(numout,*)
483         WRITE(numout,*) 'ldf_slp_init_crs : direction of lateral mixing'
484         WRITE(numout,*) '~~~~~~~~~~~~'
485      ENDIF
486
487      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes
488         ALLOCATE( triadi_g(jpi_crs,jpj_crs,jpk,0:1,0:1) , triadj_g(jpi_crs,jpj_crs,jpk,0:1,0:1) , wslp2(jpi_crs,jpj_crs,jpk) , STAT=ierr )
489         ALLOCATE( triadi  (jpi_crs,jpj_crs,jpk,0:1,0:1) , triadj  (jpi_crs,jpj_crs,jpk,0:1,0:1)                      , STAT=ierr )
490         IF( ierr > 0             )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
491         !
492         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
493         !
494      ELSE                             ! Madec operator : slopes at u-, v-, and w-points
495         ALLOCATE( uslp_crs(jpi_crs,jpj_crs,jpk) , vslp_crs(jpi_crs,jpj_crs,jpk) , wslpi_crs(jpi_crs,jpj_crs,jpk) , wslpj_crs(jpi_crs,jpj_crs,jpk) ,  &
496            &   omlmask(jpi_crs,jpj_crs,jpk) , uslpml(jpi_crs,jpj_crs)   , vslpml(jpi_crs,jpj_crs)    , wslpiml(jpi_crs,jpj_crs)   , wslpjml(jpi_crs,jpj_crs) , STAT=ierr )
497         IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
498
499         ! Direction of lateral diffusion (tracers and/or momentum)
500         ! ------------------------------
501         uslp_crs (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp      ! set the slope to zero (even in s-coordinates)
502         vslp_crs (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp
503         wslpi_crs(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp
504         wslpj_crs(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
505
506         !!gm I no longer understand this.....
507         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
508            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
509
510            ! geopotential diffusion in s-coordinates on tracers and/or momentum
511            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
512            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
513
514            ! set the slope of diffusion to the slope of s-surfaces
515            !      ( c a u t i o n : minus sign as fsdep has positive value )
516            DO jk = 1, jpk
517               DO jj = 2, jpj_crsm1
518                  DO ji = 2, jpi_crsm1   ! vector opt.
519                  !cbr uslp_crs (ji,jj,jk) = -1./e1u_crs(ji,jj) * ( gdept_crs(ji+1,jj,jk) - gdept_crs(ji ,jj ,jk) ) * umask_crs(ji,jj,jk)
520                  !vslp_crs (ji,jj,jk) = -1./e2v_crs(ji,jj) * ( gdept_crs(ji,jj+1,jk) - gdept_crs(ji ,jj ,jk) ) * vmask_crs(ji,jj,jk)
521                  !wslpi_crs(ji,jj,jk) = -1./e1t_crs(ji,jj) * ( gdepw_crs(ji+1,jj,jk) - gdepw_crs(ji-1,jj,jk) ) * tmask_crs(ji,jj,jk) * 0.5
522                  !wslpj_crs(ji,jj,jk) = -1./e2t_crs(ji,jj) * ( gdepw_crs(ji,jj+1,jk) - gdepw_crs(ji,jj-1,jk) ) * tmask_crs(ji,jj,jk) * 0.5
523                  uslp_crs (ji,jj,jk) = -1. * ( gdept_crs(ji+1,jj,jk) - gdept_crs(ji ,jj ,jk) ) * umask_crs(ji,jj,jk)
524                  IF( e1u_crs(ji,jj) .NE. 0._wp ) uslp_crs (ji,jj,jk) = uslp_crs (ji,jj,jk) / e1u_crs(ji,jj)
525                  vslp_crs (ji,jj,jk) = -1. * ( gdept_crs(ji,jj+1,jk) - gdept_crs(ji ,jj ,jk) ) * vmask_crs(ji,jj,jk)
526                  IF( e2v_crs(ji,jj) .NE. 0._wp ) vslp_crs (ji,jj,jk) = vslp_crs (ji,jj,jk) / e2v_crs(ji,jj)
527                  wslpi_crs(ji,jj,jk) = -1. * ( gdepw_crs(ji+1,jj,jk) - gdepw_crs(ji-1,jj,jk) ) * tmask_crs(ji,jj,jk) * 0.5
528                  IF( e1t_crs(ji,jj) .NE. 0._wp ) wslpi_crs(ji,jj,jk) =  wslpi_crs(ji,jj,jk) / e1t_crs(ji,jj)
529                  wslpj_crs(ji,jj,jk) = -1. * ( gdepw_crs(ji,jj+1,jk) - gdepw_crs(ji,jj-1,jk) ) * tmask_crs(ji,jj,jk) * 0.5
530                  IF( e2t_crs(ji,jj) .NE. 0._wp ) wslpj_crs(ji,jj,jk) = wslpj_crs(ji,jj,jk) / e2t_crs(ji,jj)
531                  END DO
532               END DO
533            END DO
534            CALL crs_lbc_lnk( uslp_crs , 'U', -1. )   ;   CALL crs_lbc_lnk( vslp_crs , 'V', -1. )      ! Lateral boundary conditions
535            CALL crs_lbc_lnk( wslpi_crs, 'W', -1. )   ;   CALL crs_lbc_lnk( wslpj_crs, 'W', -1. )
536         ENDIF
537      ENDIF
538      !
539      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp_init')
540      !
541   END SUBROUTINE ldf_slp_init_crs
542
543#else
544   !!------------------------------------------------------------------------
545   !!   Dummy module :                 NO Rotation of lateral mixing tensor
546   !!------------------------------------------------------------------------
547   LOGICAL, PUBLIC, PARAMETER ::   lk_ldfslp_crs = .FALSE.    !: slopes flag
548CONTAINS
549   SUBROUTINE ldf_slp_crs( kt, prd, pn2 )   ! Dummy routine
550      INTEGER, INTENT(in) :: kt
551      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
552      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
553   END SUBROUTINE ldf_slp_crs
554   SUBROUTINE ldf_slp_init_crs              ! Dummy routine
555   END SUBROUTINE ldf_slp_init_crs
556#endif
557
558   SUBROUTINE ldf_slp_grif_crs( kt )        ! Dummy routine
559      INTEGER, INTENT(in) :: kt
560      WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt
561   END SUBROUTINE ldf_slp_grif_crs
562
563   !!======================================================================
564END MODULE ldfslp_crs
Note: See TracBrowser for help on using the repository browser.