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/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/DYN/dynzdf.F90 @ 12397

Last change on this file since 12397 was 12397, checked in by davestorkey, 4 years ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2 : Consolidation of code to
handle initial Euler timestep in the context of leapfrog
timestepping. This version passes all SETTE tests but fails to bit
compare with the control for several tests (ORCA2_ICE_PISCES, AMM12,
ISOMIP, AGRIF_DEMO, SPITZ12).
  • Property svn:keywords set to Id
File size: 23.0 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   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46   
47   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_zdf  ***
50      !!
51      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
52      !!              together with the Leap-Frog time stepping using an
53      !!              implicit scheme.
54      !!
55      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
56      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
57      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise
58      !!               - update the after velocity with the implicit vertical mixing.
59      !!      This requires to solver the following system:
60      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ]
61      !!      with the following surface/top/bottom boundary condition:
62      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
63      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
64      !!
65      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity
66      !!---------------------------------------------------------------------
67      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
68      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
69      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
70      !
71      INTEGER  ::   ji, jj, jk         ! dummy loop indices
72      INTEGER  ::   iku, ikv           ! local integers
73      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
74      REAL(wp) ::   zzws, ze3va        !   -      -
75      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
76      REAL(wp) ::   zWu , zWv          !   -      -
77      REAL(wp) ::   zWui, zWvi         !   -      -
78      REAL(wp) ::   zWus, zWvs         !   -      -
79      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
80      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
81      !!---------------------------------------------------------------------
82      !
83      IF( ln_timing )   CALL timing_start('dyn_zdf')
84      !
85      IF( kt == nit000 ) THEN       !* initialization
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
89         !
90         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
91         ELSE                   ;    r_vvl = 1._wp
92         ENDIF
93      ENDIF
94      !                             !* explicit top/bottom drag case
95      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
96      !
97      !
98      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
99         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
100         ztrdu(:,:,:) = puu(:,:,:,Krhs)
101         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
102      ENDIF
103      !
104      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
105      !
106      !                    ! time stepping except vertical diffusion
107      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
108         DO jk = 1, jpkm1
109            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)
110            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)
111         END DO
112      ELSE                                      ! applied on thickness weighted velocity
113         DO jk = 1, jpkm1
114            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
115               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
116            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
117               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
118         END DO
119      ENDIF
120      !                    ! add top/bottom friction
121      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
122      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
123      !                not lead to the effective stress seen over the whole barotropic loop.
124      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
125      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
126         DO jk = 1, jpkm1        ! remove barotropic velocities
127            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)
128            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)
129         END DO
130         DO_2D_00_00
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(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
134            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
135            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
136            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
137         END_2D
138         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
139            DO_2D_00_00
140               iku = miku(ji,jj)         ! top ocean level at u- and v-points
141               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
142               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
143               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
144               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
145               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
146            END_2D
147         END IF
148      ENDIF
149      !
150      !              !==  Vertical diffusion on u  ==!
151      !
152      !                    !* Matrix construction
153      zdt = r2dt * 0.5
154      IF( ln_zad_Aimp ) THEN   !!
155         SELECT CASE( nldf_dyn )
156         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
157            DO_3D_00_00( 1, jpkm1 )
158               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
159               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
160                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
161               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
162                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
163               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
164               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
165               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
166               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
167               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
168            END_3D
169         CASE DEFAULT               ! iso-level lateral mixing
170            DO_3D_00_00( 1, jpkm1 )
171               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
172               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
173               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
174               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
175               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
176               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
177               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
178               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
179            END_3D
180         END SELECT
181         DO_2D_00_00
182            zwi(ji,jj,1) = 0._wp
183            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
184            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
185            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
186            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
187            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
188         END_2D
189      ELSE
190         SELECT CASE( nldf_dyn )
191         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
192            DO_3D_00_00( 1, jpkm1 )
193               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
194               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
195                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
196               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
197                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
198               zwi(ji,jj,jk) = zzwi
199               zws(ji,jj,jk) = zzws
200               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
201            END_3D
202         CASE DEFAULT               ! iso-level lateral mixing
203            DO_3D_00_00( 1, jpkm1 )
204               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
205               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
206               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
207               zwi(ji,jj,jk) = zzwi
208               zws(ji,jj,jk) = zzws
209               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
210            END_3D
211         END SELECT
212         DO_2D_00_00
213            zwi(ji,jj,1) = 0._wp
214            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
215         END_2D
216      ENDIF
217      !
218      !
219      !              !==  Apply semi-implicit bottom friction  ==!
220      !
221      !     Only needed for semi-implicit bottom friction setup. The explicit
222      !     bottom friction has been included in "u(v)a" which act as the R.H.S
223      !     column vector of the tri-diagonal matrix equation
224      !
225      IF ( ln_drgimp ) THEN      ! implicit bottom friction
226         DO_2D_00_00
227            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
228            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
229            zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
230         END_2D
231         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
232            DO_2D_00_00
233               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
234               iku = miku(ji,jj)       ! ocean top level at u- and v-points
235               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
236               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
237            END_2D
238         END IF
239      ENDIF
240      !
241      ! Matrix inversion starting from the first level
242      !-----------------------------------------------------------------------
243      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
244      !
245      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
246      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
247      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
248      !        (        ...               )( ...  ) ( ...  )
249      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
250      !
251      !   m is decomposed in the product of an upper and a lower triangular matrix
252      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
253      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
254      !-----------------------------------------------------------------------
255      !
256      DO_3D_00_00( 2, jpkm1 )
257         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
258      END_3D
259      !
260      DO_2D_00_00
261         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
262         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
263            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
264      END_2D
265      DO_3D_00_00( 2, jpkm1 )
266         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)
267      END_3D
268      !
269      DO_2D_00_00
270         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
271      END_2D
272      DO_3DS_00_00( jpk-2, 1, -1 )
273         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
274      END_3D
275      !
276      !              !==  Vertical diffusion on v  ==!
277      !
278      !                       !* Matrix construction
279      zdt = r2dt * 0.5
280      IF( ln_zad_Aimp ) THEN   !!
281         SELECT CASE( nldf_dyn )
282         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
283            DO_3D_00_00( 1, jpkm1 )
284               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
285               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
286                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
287               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
288                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
289               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
290               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
291               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
292               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
293               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
294            END_3D
295         CASE DEFAULT               ! iso-level lateral mixing
296            DO_3D_00_00( 1, jpkm1 )
297               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
298               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
299               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
300               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
301               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
302               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
303               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
304               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
305            END_3D
306         END SELECT
307         DO_2D_00_00
308            zwi(ji,jj,1) = 0._wp
309            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
310            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
311            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
312            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
313            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
314         END_2D
315      ELSE
316         SELECT CASE( nldf_dyn )
317         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
318            DO_3D_00_00( 1, jpkm1 )
319               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
320               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
321                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
322               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
323                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
324               zwi(ji,jj,jk) = zzwi
325               zws(ji,jj,jk) = zzws
326               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
327            END_3D
328         CASE DEFAULT               ! iso-level lateral mixing
329            DO_3D_00_00( 1, jpkm1 )
330               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
331               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
332               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
333               zwi(ji,jj,jk) = zzwi
334               zws(ji,jj,jk) = zzws
335               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
336            END_3D
337         END SELECT
338         DO_2D_00_00
339            zwi(ji,jj,1) = 0._wp
340            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
341         END_2D
342      ENDIF
343      !
344      !              !==  Apply semi-implicit top/bottom friction  ==!
345      !
346      !     Only needed for semi-implicit bottom friction setup. The explicit
347      !     bottom friction has been included in "u(v)a" which act as the R.H.S
348      !     column vector of the tri-diagonal matrix equation
349      !
350      IF( ln_drgimp ) THEN
351         DO_2D_00_00
352            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
353            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
354            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
355         END_2D
356         IF ( ln_isfcav ) THEN
357            DO_2D_00_00
358               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
359               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
360               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
361            END_2D
362         ENDIF
363      ENDIF
364
365      ! Matrix inversion
366      !-----------------------------------------------------------------------
367      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
368      !
369      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
370      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
371      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
372      !        (        ...               )( ...  ) ( ...  )
373      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
374      !
375      !   m is decomposed in the product of an upper and lower triangular matrix
376      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
377      !   The solution (after velocity) is in 2d array va
378      !-----------------------------------------------------------------------
379      !
380      DO_3D_00_00( 2, jpkm1 )
381         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
382      END_3D
383      !
384      DO_2D_00_00
385         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
386         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
387            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
388      END_2D
389      DO_3D_00_00( 2, jpkm1 )
390         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)
391      END_3D
392      !
393      DO_2D_00_00
394         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
395      END_2D
396      DO_3DS_00_00( jpk-2, 1, -1 )
397         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
398      END_3D
399      !
400      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
401         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
402         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
403         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
404         DEALLOCATE( ztrdu, ztrdv ) 
405      ENDIF
406      !                                          ! print mean trends (used for debugging)
407      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
408         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
409         !
410      IF( ln_timing )   CALL timing_stop('dyn_zdf')
411      !
412   END SUBROUTINE dyn_zdf
413
414   !!==============================================================================
415END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.