source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90 @ 10001

Last change on this file since 10001 was 10001, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step I.1 and I.2 (see wiki page)

  • Property svn:keywords set to Id
File size: 26.9 KB
Line 
1MODULE dynzdf
2   !!==============================================================================
3   !!                 ***  MODULE  dynzdf  ***
4   !! Ocean dynamics :  vertical component of the momentum mixing trend
5   !!==============================================================================
6   !! History :  1.0  !  2005-11  (G. Madec)  Original code
7   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
8   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
13   !!       zdf_trd   : diagnose the zdf velocity trends and the KE dissipation trend
14!!gm                        ==>>> zdf_trd currently not used
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE phycst         ! physical constants
18   USE dom_oce        ! ocean space and time domain variables
19   USE sbc_oce        ! surface boundary condition: ocean
20   USE zdf_oce        ! ocean vertical physics variables
21   USE zdfdrg         ! vertical physics: top/bottom drag coef.
22   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
23   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
24   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
25   USE trd_oce        ! trends: ocean variables
26   USE trddyn         ! trend manager: dynamics
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE iom             ! IOM library
31   USE prtctl         ! Print control
32   USE timing         ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   dyn_zdf   !  routine called by step.F90
38
39   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
40
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL licence (./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49   
50   SUBROUTINE dyn_zdf( kt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE dyn_zdf  ***
53      !!
54      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
55      !!              together with the Leap-Frog time stepping using an
56      !!              implicit scheme.
57      !!
58      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
59      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
60      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
61      !!               - update the after velocity with the implicit vertical mixing.
62      !!      This requires to solver the following system:
63      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
64      !!      with the following surface/top/bottom boundary condition:
65      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
66      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
67      !!
68      !! ** Action :   (ua,va)   after velocity
69      !!---------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
71      !
72      INTEGER  ::   ji, jj, jk            ! dummy loop indices
73      INTEGER  ::   iku, ikv              ! local integers
74      REAL(wp) ::   zzwi, ze3ua, z2dt_2   ! local scalars
75      REAL(wp) ::   zzws, ze3va           !   -      -
76      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace
77      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      -
78      !!---------------------------------------------------------------------
79      !
80      IF( ln_timing )   CALL timing_start('dyn_zdf')
81      !
82      IF( kt == nit000 ) THEN       !* initialization
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
86         !
87         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
88         ELSE                   ;    r_vvl = 1._wp
89         ENDIF
90      ENDIF
91      !
92      z2dt_2 = rDt * 0.5_wp        !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2
93      !
94      !
95      !                             !* explicit top/bottom drag case
96      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
97      !
98      !
99      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
100         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
101         ztrdu(:,:,:) = ua(:,:,:)
102         ztrdv(:,:,:) = va(:,:,:)
103      ENDIF
104      !
105      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
106      !
107      !                    ! time stepping except vertical diffusion
108      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
109         DO jk = 1, jpkm1
110            ua(:,:,jk) = ( ub(:,:,jk) + rDt * ua(:,:,jk) ) * umask(:,:,jk)
111            va(:,:,jk) = ( vb(:,:,jk) + rDt * va(:,:,jk) ) * vmask(:,:,jk)
112         END DO
113      ELSE                                      ! applied on thickness weighted velocity
114         DO jk = 1, jpkm1
115            ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + rDt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
116            va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + rDt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
117         END DO
118      ENDIF
119      !                    ! add top/bottom friction
120      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
121      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
122      !                not lead to the effective stress seen over the whole barotropic loop.
123      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
124      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
125         DO jk = 1, jpkm1        ! remove barotropic velocities
126            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
127            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
128         END DO
129         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
130            DO ji = fs_2, fs_jpim1   ! vector opt.
131               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
132               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
133               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
134               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
135               ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
136               va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
137            END DO
138         END DO
139         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
140            DO jj = 2, jpjm1       
141               DO ji = fs_2, fs_jpim1   ! vector opt.
142                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
143                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
144                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
145                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
146                  ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
147                  va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
148               END DO
149            END DO
150         END IF
151      ENDIF
152      !
153      !              !==  Vertical diffusion on u  ==!
154      !
155      SELECT CASE( nldf_dyn )    !* Matrix construction
156      !
157      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
158         DO jk = 1, jpkm1
159            DO jj = 2, jpjm1 
160               DO ji = fs_2, fs_jpim1   ! vector opt.
161                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
162                  zzwi = - rDt * ( 0.5 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) + akzu(ji,jj,jk  ) )   &
163                     &            / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
164                  zzws = - rDt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) )   &
165                     &            / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
166                  zwi(ji,jj,jk) = zzwi
167                  zws(ji,jj,jk) = zzws
168                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
169               END DO
170            END DO
171         END DO
172      CASE DEFAULT                  ! iso-level lateral mixing
173         DO jk = 1, jpkm1
174            DO jj = 2, jpjm1 
175               DO ji = fs_2, fs_jpim1   ! vector opt.
176                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
177                  zzwi = - z2dt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
178                  zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
179                  zwi(ji,jj,jk) = zzwi
180                  zws(ji,jj,jk) = zzws
181                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
182               END DO
183            END DO
184         END DO
185      END SELECT
186      !
187      DO jj = 2, jpjm1           !* Surface boundary conditions
188         DO ji = fs_2, fs_jpim1     ! vector opt.
189            zwi(ji,jj,1) = 0._wp
190            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
191         END DO
192      END DO
193      !
194      !              !==  Apply semi-implicit bottom friction  ==!
195      !
196      !     Only needed for semi-implicit bottom friction setup. The explicit
197      !     bottom friction has been included in "u(v)a" which act as the R.H.S
198      !     column vector of the tri-diagonal matrix equation
199      !
200      IF ( ln_drgimp ) THEN      ! implicit bottom friction
201         DO jj = 2, jpjm1
202            DO ji = 2, jpim1
203               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
204               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
205               zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
206            END DO
207         END DO
208         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
209            DO jj = 2, jpjm1
210               DO ji = 2, jpim1
211                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
212                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
213                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
214                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
215               END DO
216            END DO
217         END IF
218      ENDIF
219      !
220      ! Matrix inversion starting from the first level
221      !-----------------------------------------------------------------------
222      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
223      !
224      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
225      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
226      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
227      !        (        ...               )( ...  ) ( ...  )
228      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
229      !
230      !   m is decomposed in the product of an upper and a lower triangular matrix
231      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
232      !   The solution (the after velocity) is in ua
233      !-----------------------------------------------------------------------
234      !
235      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
236         DO jj = 2, jpjm1   
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
239            END DO
240         END DO
241      END DO
242      !
243      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
244         DO ji = fs_2, fs_jpim1   ! vector opt.
245            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
246            ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rho0 ) * umask(ji,jj,1)
247         END DO
248      END DO
249      DO jk = 2, jpkm1
250         DO jj = 2, jpjm1
251            DO ji = fs_2, fs_jpim1
252               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
253            END DO
254         END DO
255      END DO
256      !
257      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
258         DO ji = fs_2, fs_jpim1   ! vector opt.
259            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
260         END DO
261      END DO
262      DO jk = jpk-2, 1, -1
263         DO jj = 2, jpjm1
264            DO ji = fs_2, fs_jpim1
265               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
266            END DO
267         END DO
268      END DO
269      !
270      !              !==  Vertical diffusion on v  ==!
271      !
272      !                     
273      SELECT CASE( nldf_dyn )    !* Matrix construction
274      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
275         DO jk = 1, jpkm1
276            DO jj = 2, jpjm1   
277               DO ji = fs_2, fs_jpim1   ! vector opt.
278                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
279                  zzwi = - rDt * ( 0.5 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) + akzv(ji,jj,jk  ) )   &
280                     &            / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
281                  zzws = - rDt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) )   &
282                     &            / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
283                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
284                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
285                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
286               END DO
287            END DO
288         END DO
289      CASE DEFAULT                  ! iso-level lateral mixing
290         DO jk = 1, jpkm1
291            DO jj = 2, jpjm1   
292               DO ji = fs_2, fs_jpim1   ! vector opt.
293                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
294                  zzwi = - z2dt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
295                  zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
296                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
297                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
298                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
299               END DO
300            END DO
301         END DO
302      END SELECT
303      !
304      DO jj = 2, jpjm1           !* Surface boundary conditions
305         DO ji = fs_2, fs_jpim1     ! vector opt.
306            zwi(ji,jj,1) = 0._wp
307            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
308         END DO
309      END DO
310      !              !==  Apply semi-implicit top/bottom friction  ==!
311      !
312      !     Only needed for semi-implicit bottom friction setup. The explicit
313      !     bottom friction has been included in "u(v)a" which act as the R.H.S
314      !     column vector of the tri-diagonal matrix equation
315      !
316      IF( ln_drgimp ) THEN
317         DO jj = 2, jpjm1
318            DO ji = 2, jpim1
319               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
320               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
321               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
322            END DO
323         END DO
324         IF ( ln_isfcav ) THEN
325            DO jj = 2, jpjm1
326               DO ji = 2, jpim1
327                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
328                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
329                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
330               END DO
331            END DO
332         ENDIF
333      ENDIF
334
335      ! Matrix inversion
336      !-----------------------------------------------------------------------
337      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
338      !
339      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
340      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
341      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
342      !        (        ...               )( ...  ) ( ...  )
343      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
344      !
345      !   m is decomposed in the product of an upper and lower triangular matrix
346      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
347      !   The solution (after velocity) is in 2d array va
348      !-----------------------------------------------------------------------
349      !
350      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
351         DO jj = 2, jpjm1   
352            DO ji = fs_2, fs_jpim1   ! vector opt.
353               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
354            END DO
355         END DO
356      END DO
357      !
358      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
359         DO ji = fs_2, fs_jpim1   ! vector opt.         
360            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
361            va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rho0 ) * vmask(ji,jj,1)
362         END DO
363      END DO
364      DO jk = 2, jpkm1
365         DO jj = 2, jpjm1
366            DO ji = fs_2, fs_jpim1   ! vector opt.
367               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
368            END DO
369         END DO
370      END DO
371      !
372      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
373         DO ji = fs_2, fs_jpim1   ! vector opt.
374            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
375         END DO
376      END DO
377      DO jk = jpk-2, 1, -1
378         DO jj = 2, jpjm1
379            DO ji = fs_2, fs_jpim1
380               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
381            END DO
382         END DO
383      END DO
384      !
385      IF( l_trddyn ) THEN                        ! save the vertical diffusive trends for further diagnostics
386         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
387            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ztrdu(:,:,:)
388            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ztrdv(:,:,:)
389         ELSE                                      ! applied on thickness weighted velocity
390            ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ztrdu(:,:,:)
391            ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ztrdv(:,:,:)
392         ENDIF
393         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
394         DEALLOCATE( ztrdu, ztrdv ) 
395      ENDIF
396      !                                          ! print mean trends (used for debugging)
397      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
398         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
399         !
400      IF( ln_timing )   CALL timing_stop('dyn_zdf')
401      !
402   END SUBROUTINE dyn_zdf
403
404!!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing
405
406   SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt )
407      !!----------------------------------------------------------------------
408      !!                  ***  ROUTINE zdf_trd  ***
409      !!
410      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
411      !!              together with the Leap-Frog time stepping using an
412      !!              implicit scheme.
413      !!
414      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
415      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
416      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
417      !!               - update the after velocity with the implicit vertical mixing.
418      !!      This requires to solver the following system:
419      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
420      !!      with the following surface/top/bottom boundary condition:
421      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
422      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
423      !!
424      !! ** Action :   (ua,va)   after velocity
425      !!---------------------------------------------------------------------
426      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptrdu, ptrdv   ! 3D work arrays use for zdf trends diag
427      INTEGER , INTENT(in   )                         ::   kt             ! ocean time-step index
428      !
429      INTEGER  ::   ji, jj, jk       ! dummy loop indices
430      REAL(wp) ::   zzz              ! local scalar
431      REAL(wp) ::   zavmu, zavmum1   !   -      -
432      REAL(wp) ::   zavmv, zavmvm1   !   -      -
433      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2d    !  -      -
434      !!---------------------------------------------------------------------
435      !
436      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )   ! apply lateral boundary condition on (ua,va)
437      !
438      !
439      !                 !==  momentum trend due to vertical diffusion  ==!
440      !
441      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
442         ptrdu(:,:,:) = (              ua(:,:,:) -              ub(:,:,:) )                * r1_Dt - ptrdu(:,:,:)
443         ptrdv(:,:,:) = (              va(:,:,:) -              vb(:,:,:) )                * r1_Dt - ptrdv(:,:,:)
444      ELSE                                      ! applied on thickness weighted velocity
445         ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ptrdu(:,:,:)
446         ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ptrdv(:,:,:)
447      ENDIF
448      CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt )
449      !
450      !
451      !                 !==  KE dissipation trend due to vertical diffusion  ==!
452      !
453      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area
454         !                                ! due to v friction (v=vertical)
455         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points)
456         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of
457         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured).
458         !                                ! Note also that now e3 scale factors are used as after one are not computed !
459         !
460         ALLOCATE( z2d(jpi,jpj) )
461         z2d(:,:) = 0._wp
462         DO jk = 1, jpkm1
463            DO jj = 2, jpjm1
464               DO ji = 2, jpim1
465                  zavmu   = 0.5 * ( avm(ji+1,jj,jk) + avm(ji  ,jj,jk) )
466                  zavmum1 = 0.5 * ( avm(ji  ,jj,jk) + avm(ji-1,jj,jk) )
467                  zavmv   = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj  ,jk) )
468                  zavmvm1 = 0.5 * ( avm(ji,jj  ,jk) + avm(ji,jj-1,jk) )
469               
470                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
471                     &   zavmu   * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / e3uw_n(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
472                     & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
473                     & + zavmv   * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / e3vw_n(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
474                     & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
475                     &                        )
476!!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step  ===>>> use it ?
477!!                                                                                     No since in zdfshé only kz tke (or gls) is used
478!!
479!!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of
480!!gm     now by before shears, i.e. the source term of TKE (local positivity is not ensured).
481!!       CAUTION: requires to compute e3uw_a and e3vw_a !!!
482!                  z2d(ji,jj) = z2d(ji,jj)  + (                                                                                 &
483!                     &    avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) / e3uw_n(ji  ,jj,jk)                        &
484!                     &                     * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) ) / e3uw_a(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
485!                     &  + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk)                        &
486!                     &                       ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
487!                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) / e3vw_n(ji,jj  ,jk)                        &
488!                     &                       ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) ) / e3vw_a(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
489!                     &  + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk)                        &
490!                     &                       ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
491!                     &                       )
492!!gm end
493               END DO
494            END DO
495         END DO
496         zzz= - 0.5_wp* rho0           ! caution sign minus here
497         z2d(:,:) = zzz * z2d(:,:) 
498         CALL lbc_lnk( z2d,'T', 1. )
499         CALL iom_put( 'dispkevfo', z2d )
500         DEALLOCATE( z2d )
501      ENDIF
502      !
503   END SUBROUTINE zdf_trd
504   
505!!gm end not used
506
507   !!==============================================================================
508END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.