source: NEMO/trunk/src/OCE/DYN/dynzdf.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 9 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 23.3 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      !                             !* set time step
95      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
96      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
97      ENDIF
98      !
99      !                             !* explicit top/bottom drag case
100      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))
101      !
102      !
103      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
104         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
105         ztrdu(:,:,:) = puu(:,:,:,Krhs)
106         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
107      ENDIF
108      !
109      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
110      !
111      !                    ! time stepping except vertical diffusion
112      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
113         DO jk = 1, jpkm1
114            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)
115            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)
116         END DO
117      ELSE                                      ! applied on thickness weighted velocity
118         DO jk = 1, jpkm1
119            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
120               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
121            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
122               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
123         END DO
124      ENDIF
125      !                    ! add top/bottom friction
126      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
127      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
128      !                not lead to the effective stress seen over the whole barotropic loop.
129      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
130      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
131         DO jk = 1, jpkm1        ! remove barotropic velocities
132            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)
133            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)
134         END DO
135         DO_2D_00_00
136            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
137            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
138            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
139            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
140            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
141            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
142         END_2D
143         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
144            DO_2D_00_00
145               iku = miku(ji,jj)         ! top ocean level at u- and v-points
146               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
147               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
148               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
149               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
150               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
151            END_2D
152         END IF
153      ENDIF
154      !
155      !              !==  Vertical diffusion on u  ==!
156      !
157      !                    !* Matrix construction
158      zdt = r2dt * 0.5
159      IF( ln_zad_Aimp ) THEN   !!
160         SELECT CASE( nldf_dyn )
161         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
162            DO_3D_00_00( 1, jpkm1 )
163               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
164               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
165                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
166               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
167                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
168               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
169               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
170               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
171               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
172               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
173            END_3D
174         CASE DEFAULT               ! iso-level lateral mixing
175            DO_3D_00_00( 1, jpkm1 )
176               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + 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  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
178               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)
179               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
180               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
181               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
182               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
183               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
184            END_3D
185         END SELECT
186         DO_2D_00_00
187            zwi(ji,jj,1) = 0._wp
188            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
189            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
190            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
191            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
192            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
193         END_2D
194      ELSE
195         SELECT CASE( nldf_dyn )
196         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
197            DO_3D_00_00( 1, jpkm1 )
198               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
199               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
200                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
201               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
202                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
203               zwi(ji,jj,jk) = zzwi
204               zws(ji,jj,jk) = zzws
205               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
206            END_3D
207         CASE DEFAULT               ! iso-level lateral mixing
208            DO_3D_00_00( 1, jpkm1 )
209               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
210               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
211               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)
212               zwi(ji,jj,jk) = zzwi
213               zws(ji,jj,jk) = zzws
214               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
215            END_3D
216         END SELECT
217         DO_2D_00_00
218            zwi(ji,jj,1) = 0._wp
219            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
220         END_2D
221      ENDIF
222      !
223      !
224      !              !==  Apply semi-implicit bottom friction  ==!
225      !
226      !     Only needed for semi-implicit bottom friction setup. The explicit
227      !     bottom friction has been included in "u(v)a" which act as the R.H.S
228      !     column vector of the tri-diagonal matrix equation
229      !
230      IF ( ln_drgimp ) THEN      ! implicit bottom friction
231         DO_2D_00_00
232            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
233            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
234            zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
235         END_2D
236         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
237            DO_2D_00_00
238               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
239               iku = miku(ji,jj)       ! ocean top level at u- and v-points
240               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
241               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
242            END_2D
243         END IF
244      ENDIF
245      !
246      ! Matrix inversion starting from the first level
247      !-----------------------------------------------------------------------
248      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
249      !
250      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
251      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
252      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
253      !        (        ...               )( ...  ) ( ...  )
254      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
255      !
256      !   m is decomposed in the product of an upper and a lower triangular matrix
257      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
258      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
259      !-----------------------------------------------------------------------
260      !
261      DO_3D_00_00( 2, jpkm1 )
262         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
263      END_3D
264      !
265      DO_2D_00_00
266         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
267         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
268            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
269      END_2D
270      DO_3D_00_00( 2, jpkm1 )
271         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)
272      END_3D
273      !
274      DO_2D_00_00
275         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
276      END_2D
277      DO_3DS_00_00( jpk-2, 1, -1 )
278         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
279      END_3D
280      !
281      !              !==  Vertical diffusion on v  ==!
282      !
283      !                       !* Matrix construction
284      zdt = r2dt * 0.5
285      IF( ln_zad_Aimp ) THEN   !!
286         SELECT CASE( nldf_dyn )
287         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
288            DO_3D_00_00( 1, jpkm1 )
289               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
290               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
291                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
292               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
293                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
294               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
295               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
296               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
297               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
298               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
299            END_3D
300         CASE DEFAULT               ! iso-level lateral mixing
301            DO_3D_00_00( 1, jpkm1 )
302               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
303               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
304               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)
305               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
306               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
307               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
308               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
309               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
310            END_3D
311         END SELECT
312         DO_2D_00_00
313            zwi(ji,jj,1) = 0._wp
314            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
315            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
316            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
317            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
318            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
319         END_2D
320      ELSE
321         SELECT CASE( nldf_dyn )
322         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
323            DO_3D_00_00( 1, jpkm1 )
324               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
325               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
326                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
327               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
328                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
329               zwi(ji,jj,jk) = zzwi
330               zws(ji,jj,jk) = zzws
331               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
332            END_3D
333         CASE DEFAULT               ! iso-level lateral mixing
334            DO_3D_00_00( 1, jpkm1 )
335               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
336               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
337               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)
338               zwi(ji,jj,jk) = zzwi
339               zws(ji,jj,jk) = zzws
340               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
341            END_3D
342         END SELECT
343         DO_2D_00_00
344            zwi(ji,jj,1) = 0._wp
345            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
346         END_2D
347      ENDIF
348      !
349      !              !==  Apply semi-implicit top/bottom friction  ==!
350      !
351      !     Only needed for semi-implicit bottom friction setup. The explicit
352      !     bottom friction has been included in "u(v)a" which act as the R.H.S
353      !     column vector of the tri-diagonal matrix equation
354      !
355      IF( ln_drgimp ) THEN
356         DO_2D_00_00
357            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
358            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
359            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
360         END_2D
361         IF ( ln_isfcav ) THEN
362            DO_2D_00_00
363               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
364               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
365               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
366            END_2D
367         ENDIF
368      ENDIF
369
370      ! Matrix inversion
371      !-----------------------------------------------------------------------
372      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
373      !
374      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
375      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
376      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
377      !        (        ...               )( ...  ) ( ...  )
378      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
379      !
380      !   m is decomposed in the product of an upper and lower triangular matrix
381      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
382      !   The solution (after velocity) is in 2d array va
383      !-----------------------------------------------------------------------
384      !
385      DO_3D_00_00( 2, jpkm1 )
386         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
387      END_3D
388      !
389      DO_2D_00_00
390         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
391         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
392            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
393      END_2D
394      DO_3D_00_00( 2, jpkm1 )
395         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)
396      END_3D
397      !
398      DO_2D_00_00
399         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
400      END_2D
401      DO_3DS_00_00( jpk-2, 1, -1 )
402         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
403      END_3D
404      !
405      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
406         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
407         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
408         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
409         DEALLOCATE( ztrdu, ztrdv ) 
410      ENDIF
411      !                                          ! print mean trends (used for debugging)
412      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
413         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
414         !
415      IF( ln_timing )   CALL timing_stop('dyn_zdf')
416      !
417   END SUBROUTINE dyn_zdf
418
419   !!==============================================================================
420END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.