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.
dynzdf.F90 in NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN – NEMO

source: NEMO/branches/2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends/src/OCE/DYN/dynzdf.F90 @ 14725

Last change on this file since 14725 was 14725, checked in by davestorkey, 3 years ago

2021/ENHANCE-01_davestorkey_fix_3D_momentum_trends : bug fix
  • Property svn:keywords set to Id
File size: 26.4 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   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain variables
17   USE sbc_oce        ! surface boundary condition: ocean
18   USE zdf_oce        ! ocean vertical physics variables
19   USE zdfdrg         ! vertical physics: top/bottom drag coef.
20   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
21   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
22   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
23   USE trd_oce        ! trends: ocean variables
24   USE trddyn         ! trend manager: dynamics
25   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE prtctl         ! Print control
29   USE timing         ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_zdf   !  routine called by step.F90
35
36   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
37
38   !! * Substitutions
39#  include "do_loop_substitute.h90"
40#  include "domzgr_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47   
48   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dyn_zdf  ***
51      !!
52      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
53      !!              together with the Leap-Frog time stepping using an
54      !!              implicit scheme.
55      !!
56      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
57      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
58      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after   otherwise
59      !!               - update the after velocity with the implicit vertical mixing.
60      !!      This requires to solver the following system:
61      !!         u(after) = u(after) + 1/e3u_after  dk+1[ mi(avm) / e3uw_after dk[ua] ]
62      !!      with the following surface/top/bottom boundary condition:
63      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
64      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
65      !!
66      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity
67      !!---------------------------------------------------------------------
68      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
69      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
70      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
71      !
72      INTEGER  ::   ji, jj, jk         ! dummy loop indices
73      INTEGER  ::   iku, ikv           ! local integers
74      INTEGER  ::   ikbu, ikbv         ! local integers
75      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
76      REAL(wp) ::   zzws, ze3va        !   -      -
77      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
78      REAL(wp) ::   zWu , zWv          !   -      -
79      REAL(wp) ::   zWui, zWvi         !   -      -
80      REAL(wp) ::   zWus, zWvs         !   -      -
81      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
82      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
83      !!---------------------------------------------------------------------
84      !
85      IF( ln_timing )   CALL timing_start('dyn_zdf')
86      !
87      IF( kt == nit000 ) THEN       !* initialization
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
90         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
91         !
92         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
93         ELSE                   ;    r_vvl = 1._wp
94         ENDIF
95      ENDIF
96      !                             !* explicit top/bottom drag case
97      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, Kaa, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
98      !
99      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
100      !
101      !                    ! time stepping except vertical diffusion
102      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
103         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
104            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
105            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
106         END_3D
107      ELSE                                      ! applied on thickness weighted velocity
108         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
109            puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  &
110               &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) &
111               &                        / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
112            pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb )  &
113               &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) &
114               &                        / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
115         END_3D
116      ENDIF
117      !                    ! add top/bottom friction
118      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
119      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
120      !                not lead to the effective stress seen over the whole barotropic loop.
121      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
122      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
123         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      ! remove barotropic velocities
124            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
125            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
126         END_3D
127         IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
128            ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
129            ztrdu(:,:,:) = puu(:,:,:,Kaa)
130            ztrdv(:,:,:) = pvv(:,:,:,Kaa)
131         ENDIF
132         DO_2D( 0, 0, 0, 0 )      ! Add bottom/top stress due to barotropic component only
133            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
134            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
135            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
136               &             + r_vvl   * e3u(ji,jj,iku,Kaa)
137            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
138               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
139            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
140            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
141         END_2D
142         IF( l_trddyn )   THEN                      ! save explicit part of bottom friction trends
143            ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - ztrdu(:,:,:) ) / rDt 
144            ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - ztrdv(:,:,:) ) / rDt 
145            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfr, kt, Kmm, Kaa )
146         ENDIF
147         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
148            DO_2D( 0, 0, 0, 0 )
149               iku = miku(ji,jj)         ! top ocean level at u- and v-points
150               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
151               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
152                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)
153               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
154                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
155               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
156               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
157            END_2D
158         END IF
159      ENDIF
160      !
161      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
162         IF( .NOT. ALLOCATED(ztrdu) ) ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
163         ztrdu(:,:,:) = puu(:,:,:,Kaa)
164         ztrdv(:,:,:) = pvv(:,:,:,Kaa)
165      ENDIF
166      !
167      !              !==  Vertical diffusion on u  ==!
168      !
169      !                    !* Matrix construction
170      zdt = rDt * 0.5
171      IF( ln_zad_Aimp ) THEN   !!
172         SELECT CASE( nldf_dyn )
173         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
174            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
175               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
176                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
177               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
178                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
179               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
180                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
181               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
182               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
183               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
184               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
185               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
186            END_3D
187         CASE DEFAULT               ! iso-level lateral mixing
188            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
189               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point
190                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)
191               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
192                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
193               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
194                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
195               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
196               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
197               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
198               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
199               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
200            END_3D
201         END SELECT
202         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
203            zwi(ji,jj,1) = 0._wp
204            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
205               &             + r_vvl   * e3u(ji,jj,1,Kaa)
206            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
207               &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
208            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
209            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
210            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
211         END_2D
212      ELSE
213         SELECT CASE( nldf_dyn )
214         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
215            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
216               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
217                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
218               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
219                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
220               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
221                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
222               zwi(ji,jj,jk) = zzwi
223               zws(ji,jj,jk) = zzws
224               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
225            END_3D
226         CASE DEFAULT               ! iso-level lateral mixing
227            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
228               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
229                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
230               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
231                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
232               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
233                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
234               zwi(ji,jj,jk) = zzwi
235               zws(ji,jj,jk) = zzws
236               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
237            END_3D
238         END SELECT
239         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
240            zwi(ji,jj,1) = 0._wp
241            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
242         END_2D
243      ENDIF
244      !
245      !
246      !              !==  Apply semi-implicit bottom friction  ==!
247      !
248      !     Only needed for semi-implicit bottom friction setup. The explicit
249      !     bottom friction has been included in "u(v)a" which act as the R.H.S
250      !     column vector of the tri-diagonal matrix equation
251      !
252      IF ( ln_drgimp ) THEN      ! implicit bottom friction
253         DO_2D( 0, 0, 0, 0 )
254            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
255            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
256               &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
257            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
258         END_2D
259         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
260            DO_2D( 0, 0, 0, 0 )
261               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
262               iku = miku(ji,jj)       ! ocean top level at u- and v-points
263               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
264                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
265               zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
266            END_2D
267         END IF
268      ENDIF
269      !
270      ! Matrix inversion starting from the first level
271      !-----------------------------------------------------------------------
272      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
273      !
274      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
275      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
276      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
277      !        (        ...               )( ...  ) ( ...  )
278      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
279      !
280      !   m is decomposed in the product of an upper and a lower triangular matrix
281      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
282      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
283      !-----------------------------------------------------------------------
284      !
285      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
286         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
287      END_3D
288      !
289      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
290         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
291            &             + r_vvl   * e3u(ji,jj,1,Kaa) 
292         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
293            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
294      END_2D
295      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
296         puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
297      END_3D
298      !
299      DO_2D( 0, 0, 0, 0 )             !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
300         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
301      END_2D
302      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
303         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
304      END_3D
305      !
306      !              !==  Vertical diffusion on v  ==!
307      !
308      !                       !* Matrix construction
309      zdt = rDt * 0.5
310      IF( ln_zad_Aimp ) THEN   !!
311         SELECT CASE( nldf_dyn )
312         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
313            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
314               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
315                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
316               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
317                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
318               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
319                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
320               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
321               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
322               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
323               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
324               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
325            END_3D
326         CASE DEFAULT               ! iso-level lateral mixing
327            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
328               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
329                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
330               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
331                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
332               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
333                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
334               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
335               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
336               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
337               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
338               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
339            END_3D
340         END SELECT
341         DO_2D( 0, 0, 0, 0 )   !* Surface boundary conditions
342            zwi(ji,jj,1) = 0._wp
343            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
344               &             + r_vvl   * e3v(ji,jj,1,Kaa)
345            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
346               &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
347            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
348            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
349            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
350         END_2D
351      ELSE
352         SELECT CASE( nldf_dyn )
353         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
354            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
355               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
356                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
357               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
358                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
359               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
360                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
361               zwi(ji,jj,jk) = zzwi
362               zws(ji,jj,jk) = zzws
363               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
364            END_3D
365         CASE DEFAULT               ! iso-level lateral mixing
366            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
367               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
368                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
369               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
370                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
371               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
372                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
373               zwi(ji,jj,jk) = zzwi
374               zws(ji,jj,jk) = zzws
375               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
376            END_3D
377         END SELECT
378         DO_2D( 0, 0, 0, 0 )        !* Surface boundary conditions
379            zwi(ji,jj,1) = 0._wp
380            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
381         END_2D
382      ENDIF
383      !
384      !              !==  Apply semi-implicit top/bottom friction  ==!
385      !
386      !     Only needed for semi-implicit bottom friction setup. The explicit
387      !     bottom friction has been included in "u(v)a" which act as the R.H.S
388      !     column vector of the tri-diagonal matrix equation
389      !
390      IF( ln_drgimp ) THEN
391         DO_2D( 0, 0, 0, 0 )
392            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
393            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
394               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
395            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
396         END_2D
397         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
398            DO_2D( 0, 0, 0, 0 )
399               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
400               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
401                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
402               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
403            END_2D
404         ENDIF
405      ENDIF
406
407      ! Matrix inversion
408      !-----------------------------------------------------------------------
409      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
410      !
411      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
412      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
413      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
414      !        (        ...               )( ...  ) ( ...  )
415      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
416      !
417      !   m is decomposed in the product of an upper and lower triangular matrix
418      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
419      !   The solution (after velocity) is in 2d array va
420      !-----------------------------------------------------------------------
421      !
422      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
423         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
424      END_3D
425      !
426      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
427         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
428            &             + r_vvl   * e3v(ji,jj,1,Kaa) 
429         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
430            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
431      END_2D
432      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
433         pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
434      END_3D
435      !
436      DO_2D( 0, 0, 0, 0 )             !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
437         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
438      END_2D
439      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
440         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
441      END_3D
442      !
443      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
444         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - ztrdu(:,:,:) ) / rDt 
445         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - ztrdv(:,:,:) ) / rDt
446         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm, Kaa )
447         !
448         IF( ln_drgimp ) THEN
449            ztrdu(:,:,:) = 0._wp    ;   ztrdv(:,:,:) = 0._wp
450            DO_2D( 0, 0, 0, 0 )
451               ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
452               ikbv = mbkv(ji,jj)
453               ztrdu(ji,jj,ikbu) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) )  & 
454     &                         * puu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm)
455               ztrdv(ji,jj,ikbv) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) )  &
456     &                         * pvv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm)
457            END_2D
458            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfri, kt, Kmm )
459         ENDIF
460         !
461         DEALLOCATE( ztrdu, ztrdv ) 
462      ENDIF
463      !                                          ! print mean trends (used for debugging)
464      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
465         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
466         !
467      IF( ln_timing )   CALL timing_stop('dyn_zdf')
468      !
469   END SUBROUTINE dyn_zdf
470
471   !!==============================================================================
472END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.