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/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynzdf.F90 @ 12616

Last change on this file since 12616 was 12616, checked in by techene, 4 years ago

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