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/UKMO/NEMO_4.0.4_momentum_trends/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_momentum_trends/src/OCE/DYN/dynzdf.F90 @ 15169

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

UKMO/NEMO_4.0.4_momentum_trends: Bug fix. Remove double counting of trends
due to top and bottom stress from barotropic currents in ZDF trend.
File size: 29.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 iom             ! IOM library
28   USE lib_mpp        ! MPP library
29   USE prtctl         ! Print control
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_zdf   !  routine called by step.F90
36
37   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
38
39   !! * Substitutions
40#  include "vectopt_loop_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 )
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      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
58      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
59      !!               - update the after velocity with the implicit vertical mixing.
60      !!      This requires to solver the following system:
61      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a 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 :   (ua,va)   after velocity
67      !!---------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
69      !
70      INTEGER  ::   ji, jj, jk         ! dummy loop indices
71      INTEGER  ::   iku, ikv           ! local integers
72      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
73      REAL(wp) ::   zzws, ze3va        !   -      -
74      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
75      REAL(wp) ::   zWu , zWv          !   -      -
76      REAL(wp) ::   zWui, zWvi         !   -      -
77      REAL(wp) ::   zWus, zWvs         !   -      -
78      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
79      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv         !  -      -
80      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu_fr, ztrdv_fr   !  -      -
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, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
101      !
102      !
103      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
104      !
105      !                    ! time stepping except vertical diffusion
106      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
107         DO jk = 1, jpkm1
108            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
109            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
110         END DO
111      ELSE                                      ! applied on thickness weighted velocity
112         DO jk = 1, jpkm1
113            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
114               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
115            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
116               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
117         END DO
118      ENDIF
119      !
120      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
121         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
122         ztrdu(:,:,:) = ua(:,:,:)
123         ztrdv(:,:,:) = va(:,:,:)
124      ENDIF
125      !
126      !                    ! add top/bottom friction
127      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
128      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
129      !                not lead to the effective stress seen over the whole barotropic loop.
130      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
131      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
132         DO jk = 1, jpkm1        ! remove barotropic velocities
133            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
134            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
135         END DO
136         IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
137            ALLOCATE( ztrdu_fr(jpi,jpj,jpk), ztrdv_fr(jpi,jpj,jpk) ) 
138            ztrdu_fr(:,:,:) = ua(:,:,:)
139            ztrdv_fr(:,:,:) = va(:,:,:)
140         ENDIF
141         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
142            DO ji = fs_2, fs_jpim1   ! vector opt.
143               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
144               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
145               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
146               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
147               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
148               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
149            END DO
150         END DO
151         IF( l_trddyn )   THEN    ! save baroclinic bottom friction trends due to barotropic currents 
152            ztrdu_fr(:,:,:) = ( ua(:,:,:) - ztrdu_fr(:,:,:) ) / r2dt 
153            ztrdv_fr(:,:,:) = ( va(:,:,:) - ztrdv_fr(:,:,:) ) / r2dt 
154            CALL trd_dyn( ztrdu_fr, ztrdv_fr, jpdyn_bfre_bt, kt )
155         ENDIF
156         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
157            IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
158               ztrdu_fr(:,:,:) = ua(:,:,:)
159               ztrdv_fr(:,:,:) = va(:,:,:)
160            ENDIF
161            DO jj = 2, jpjm1       
162               DO ji = fs_2, fs_jpim1   ! vector opt.
163                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
164                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
165                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
166                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
167                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
168                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
169               END DO
170            END DO
171            IF( l_trddyn )   THEN   ! save baroclinic top friction trends due to barotropic currents
172               ztrdu_fr(:,:,:) = ( ua(:,:,:) - ztrdu_fr(:,:,:) ) / r2dt
173               ztrdv_fr(:,:,:) = ( va(:,:,:) - ztrdv_fr(:,:,:) ) / r2dt 
174               CALL trd_dyn( ztrdu_fr, ztrdv_fr, jpdyn_tfre_bt, kt )
175            ENDIF
176         END IF
177     ENDIF
178      !
179      !              !==  Vertical diffusion on u  ==!
180      !
181      !                    !* Matrix construction
182      zdt = r2dt * 0.5
183      IF( ln_zad_Aimp ) THEN   !!
184         SELECT CASE( nldf_dyn )
185         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
186            DO jk = 1, jpkm1
187               DO jj = 2, jpjm1 
188                  DO ji = fs_2, fs_jpim1   ! vector opt.
189                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
190                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
191                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
192                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
193                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
194                     zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
195                     zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
196                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
197                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
198                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
199                  END DO
200               END DO
201            END DO
202         CASE DEFAULT               ! iso-level lateral mixing
203            DO jk = 1, jpkm1
204               DO jj = 2, jpjm1 
205                  DO ji = fs_2, fs_jpim1   ! vector opt.
206                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
207                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
208                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
209                     zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
210                     zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
211                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
212                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
213                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
214                  END DO
215               END DO
216            END DO
217         END SELECT
218         DO jj = 2, jpjm1     !* Surface boundary conditions
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220               zwi(ji,jj,1) = 0._wp
221               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
222               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
223               zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
224               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
225               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
226            END DO
227         END DO
228      ELSE
229         SELECT CASE( nldf_dyn )
230         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
231            DO jk = 1, jpkm1
232               DO jj = 2, jpjm1 
233                  DO ji = fs_2, fs_jpim1   ! vector opt.
234                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
235                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
236                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
237                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
238                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
239                     zwi(ji,jj,jk) = zzwi
240                     zws(ji,jj,jk) = zzws
241                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
242                  END DO
243               END DO
244            END DO
245         CASE DEFAULT               ! iso-level lateral mixing
246            DO jk = 1, jpkm1
247               DO jj = 2, jpjm1 
248                  DO ji = fs_2, fs_jpim1   ! vector opt.
249                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
250                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
251                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
252                     zwi(ji,jj,jk) = zzwi
253                     zws(ji,jj,jk) = zzws
254                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
255                  END DO
256               END DO
257            END DO
258         END SELECT
259         DO jj = 2, jpjm1     !* Surface boundary conditions
260            DO ji = fs_2, fs_jpim1   ! vector opt.
261               zwi(ji,jj,1) = 0._wp
262               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
263            END DO
264         END DO
265      ENDIF
266      !
267      !
268      !              !==  Apply semi-implicit bottom friction  ==!
269      !
270      !     Only needed for semi-implicit bottom friction setup. The explicit
271      !     bottom friction has been included in "u(v)a" which act as the R.H.S
272      !     column vector of the tri-diagonal matrix equation
273      !
274      IF ( ln_drgimp ) THEN      ! implicit bottom friction
275         DO jj = 2, jpjm1
276            DO ji = 2, jpim1
277               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
278               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
279               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
280            END DO
281         END DO
282         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
283            DO jj = 2, jpjm1
284               DO ji = 2, jpim1
285                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
286                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
287                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
288                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
289               END DO
290            END DO
291         END IF
292      ENDIF
293      !
294      ! Matrix inversion starting from the first level
295      !-----------------------------------------------------------------------
296      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
297      !
298      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
299      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
300      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
301      !        (        ...               )( ...  ) ( ...  )
302      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
303      !
304      !   m is decomposed in the product of an upper and a lower triangular matrix
305      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
306      !   The solution (the after velocity) is in ua
307      !-----------------------------------------------------------------------
308      !
309      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
310         DO jj = 2, jpjm1   
311            DO ji = fs_2, fs_jpim1   ! vector opt.
312               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
313            END DO
314         END DO
315      END DO
316      !
317      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
318         DO ji = fs_2, fs_jpim1   ! vector opt.
319            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
320            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
321               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
322         END DO
323      END DO
324      DO jk = 2, jpkm1
325         DO jj = 2, jpjm1
326            DO ji = fs_2, fs_jpim1
327               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
328            END DO
329         END DO
330      END DO
331      !
332      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
333         DO ji = fs_2, fs_jpim1   ! vector opt.
334            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
335         END DO
336      END DO
337      DO jk = jpk-2, 1, -1
338         DO jj = 2, jpjm1
339            DO ji = fs_2, fs_jpim1
340               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
341            END DO
342         END DO
343      END DO
344      !
345      !              !==  Vertical diffusion on v  ==!
346      !
347      !                       !* Matrix construction
348      zdt = r2dt * 0.5
349      IF( ln_zad_Aimp ) THEN   !!
350         SELECT CASE( nldf_dyn )
351         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
352            DO jk = 1, jpkm1
353               DO jj = 2, jpjm1 
354                  DO ji = fs_2, fs_jpim1   ! vector opt.
355                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
356                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
357                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
358                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
359                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
360                     zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
361                     zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
362                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
363                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
364                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
365                  END DO
366               END DO
367            END DO
368         CASE DEFAULT               ! iso-level lateral mixing
369            DO jk = 1, jpkm1
370               DO jj = 2, jpjm1 
371                  DO ji = fs_2, fs_jpim1   ! vector opt.
372                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
373                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
374                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
375                     zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
376                     zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
377                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
378                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
379                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
380                  END DO
381               END DO
382            END DO
383         END SELECT
384         DO jj = 2, jpjm1     !* Surface boundary conditions
385            DO ji = fs_2, fs_jpim1   ! vector opt.
386               zwi(ji,jj,1) = 0._wp
387               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
388               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
389               zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
390               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
391               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
392            END DO
393         END DO
394      ELSE
395         SELECT CASE( nldf_dyn )
396         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
397            DO jk = 1, jpkm1
398               DO jj = 2, jpjm1   
399                  DO ji = fs_2, fs_jpim1   ! vector opt.
400                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
401                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
402                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
403                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
404                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
405                     zwi(ji,jj,jk) = zzwi
406                     zws(ji,jj,jk) = zzws
407                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
408                  END DO
409               END DO
410            END DO
411         CASE DEFAULT               ! iso-level lateral mixing
412            DO jk = 1, jpkm1
413               DO jj = 2, jpjm1   
414                  DO ji = fs_2, fs_jpim1   ! vector opt.
415                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
416                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
417                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
418                     zwi(ji,jj,jk) = zzwi
419                     zws(ji,jj,jk) = zzws
420                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
421                  END DO
422               END DO
423            END DO
424         END SELECT
425         DO jj = 2, jpjm1        !* Surface boundary conditions
426            DO ji = fs_2, fs_jpim1   ! vector opt.
427               zwi(ji,jj,1) = 0._wp
428               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
429            END DO
430         END DO
431      ENDIF
432      !
433      !              !==  Apply semi-implicit top/bottom friction  ==!
434      !
435      !     Only needed for semi-implicit bottom friction setup. The explicit
436      !     bottom friction has been included in "u(v)a" which act as the R.H.S
437      !     column vector of the tri-diagonal matrix equation
438      !
439      IF( ln_drgimp ) THEN
440         DO jj = 2, jpjm1
441            DO ji = 2, jpim1
442               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
443               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
444               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
445            END DO
446         END DO
447         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
448            DO jj = 2, jpjm1
449               DO ji = 2, jpim1
450                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
451                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
452                  zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
453               END DO
454            END DO
455         ENDIF
456      ENDIF
457
458      ! Matrix inversion
459      !-----------------------------------------------------------------------
460      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
461      !
462      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
463      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
464      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
465      !        (        ...               )( ...  ) ( ...  )
466      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
467      !
468      !   m is decomposed in the product of an upper and lower triangular matrix
469      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
470      !   The solution (after velocity) is in 2d array va
471      !-----------------------------------------------------------------------
472      !
473      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
474         DO jj = 2, jpjm1   
475            DO ji = fs_2, fs_jpim1   ! vector opt.
476               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
477            END DO
478         END DO
479      END DO
480      !
481      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
482         DO ji = fs_2, fs_jpim1   ! vector opt.         
483            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
484            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
485               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
486         END DO
487      END DO
488      DO jk = 2, jpkm1
489         DO jj = 2, jpjm1
490            DO ji = fs_2, fs_jpim1   ! vector opt.
491               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
492            END DO
493         END DO
494      END DO
495      !
496      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
497         DO ji = fs_2, fs_jpim1   ! vector opt.
498            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
499         END DO
500      END DO
501      DO jk = jpk-2, 1, -1
502         DO jj = 2, jpjm1
503            DO ji = fs_2, fs_jpim1
504               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
505            END DO
506         END DO
507      END DO
508      !
509      IF( l_trddyn )   THEN     
510         ! Save the vertical diffusive trends for further diagnostics
511         ! then calculate implicit drag trends
512         ! (The order of operation is important here).
513         !
514         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt 
515         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt 
516         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
517         !
518         IF( ln_drgimp ) THEN
519            ztrdu_fr(:,:,:) = 0._wp    ;   ztrdv_fr(:,:,:) = 0._wp
520            DO jj = 2, jpjm1
521               DO ji = 2, jpim1
522                  iku = mbku(ji,jj)          ! deepest ocean u- & v-levels
523                  ikv = mbkv(ji,jj)
524                  ztrdu_fr(ji,jj,iku) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) )  & 
525     &                                 * un(ji,jj,iku) / e3u_n(ji,jj,iku)
526                  ztrdv_fr(ji,jj,ikv) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) )  &
527     &                                 * vn(ji,jj,ikv) / e3v_n(ji,jj,ikv)
528               END DO
529            END DO
530            CALL trd_dyn( ztrdu_fr, ztrdv_fr, jpdyn_bfri, kt )
531            IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF) or implicit ice-ocean drag
532               ztrdu_fr(:,:,:) = 0._wp    ;   ztrdv_fr(:,:,:) = 0._wp
533               DO jj = 2, jpjm1
534                  DO ji = 2, jpim1
535                     iku = miku(ji,jj)       ! ocean top level at u-points
536                     ikv = mikv(ji,jj)       ! ocean top level at v-points
537                     ztrdu_fr(ji,jj,iku) = 0.5 * ( rCdU_top(ji+1,jj) + rCdU_top(ji,jj) )  & 
538     &                                      * un(ji,jj,iku) / e3u_n(ji,jj,iku)
539                     ztrdv_fr(ji,jj,ikv) = 0.5 * ( rCdU_top(ji,jj+1) + rCdU_top(ji,jj) )  &
540     &                                      * vn(ji,jj,ikv) / e3v_n(ji,jj,ikv)
541                  END DO
542               END DO
543               CALL trd_dyn( ztrdu_fr, ztrdv_fr, jpdyn_tfri, kt )
544            ENDIF
545         ENDIF
546         !
547         IF( ALLOCATED(ztrdu_fr) ) DEALLOCATE( ztrdu_fr, ztrdv_fr ) 
548         DEALLOCATE( ztrdu, ztrdv ) 
549      ENDIF
550      !                                          ! print mean trends (used for debugging)
551      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
552         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
553         !
554      IF( ln_timing )   CALL timing_stop('dyn_zdf')
555      !
556   END SUBROUTINE dyn_zdf
557
558   !!==============================================================================
559END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.