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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 27.1 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 "vectopt_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 jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
136            DO ji = fs_2, fs_jpim1   ! vector opt.
137               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
138               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
139               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
140               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
141               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
142               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
143            END DO
144         END DO
145         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
146            DO jj = 2, jpjm1       
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
149                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
150                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
151                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
152                  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
153                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
154               END DO
155            END DO
156         END IF
157      ENDIF
158      !
159      !              !==  Vertical diffusion on u  ==!
160      !
161      !                    !* Matrix construction
162      zdt = r2dt * 0.5
163      IF( ln_zad_Aimp ) THEN   !!
164         SELECT CASE( nldf_dyn )
165         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
166            DO jk = 1, jpkm1
167               DO jj = 2, jpjm1 
168                  DO ji = fs_2, fs_jpim1   ! vector opt.
169                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
170                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
171                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
172                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
173                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
174                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
175                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
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 DO
180               END DO
181            END DO
182         CASE DEFAULT               ! iso-level lateral mixing
183            DO jk = 1, jpkm1
184               DO jj = 2, jpjm1 
185                  DO ji = fs_2, fs_jpim1   ! vector opt.
186                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
187                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
188                     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)
189                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
190                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
191                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
192                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
193                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
194                  END DO
195               END DO
196            END DO
197         END SELECT
198         DO jj = 2, jpjm1     !* Surface boundary conditions
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               zwi(ji,jj,1) = 0._wp
201               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
202               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
203               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) )
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 DO
207         END DO
208      ELSE
209         SELECT CASE( nldf_dyn )
210         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
211            DO jk = 1, jpkm1
212               DO jj = 2, jpjm1 
213                  DO ji = fs_2, fs_jpim1   ! vector opt.
214                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
215                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
216                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
217                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
218                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
219                     zwi(ji,jj,jk) = zzwi
220                     zws(ji,jj,jk) = zzws
221                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
222                  END DO
223               END DO
224            END DO
225         CASE DEFAULT               ! iso-level lateral mixing
226            DO jk = 1, jpkm1
227               DO jj = 2, jpjm1 
228                  DO ji = fs_2, fs_jpim1   ! vector opt.
229                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
230                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
231                     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)
232                     zwi(ji,jj,jk) = zzwi
233                     zws(ji,jj,jk) = zzws
234                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
235                  END DO
236               END DO
237            END DO
238         END SELECT
239         DO jj = 2, jpjm1     !* Surface boundary conditions
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zwi(ji,jj,1) = 0._wp
242               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
243            END DO
244         END DO
245      ENDIF
246      !
247      !
248      !              !==  Apply semi-implicit bottom friction  ==!
249      !
250      !     Only needed for semi-implicit bottom friction setup. The explicit
251      !     bottom friction has been included in "u(v)a" which act as the R.H.S
252      !     column vector of the tri-diagonal matrix equation
253      !
254      IF ( ln_drgimp ) THEN      ! implicit bottom friction
255         DO jj = 2, jpjm1
256            DO ji = 2, jpim1
257               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
258               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
259               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
260            END DO
261         END DO
262         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
263            DO jj = 2, jpjm1
264               DO ji = 2, jpim1
265                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
266                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
267                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
268                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
269               END DO
270            END DO
271         END IF
272      ENDIF
273      !
274      ! Matrix inversion starting from the first level
275      !-----------------------------------------------------------------------
276      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
277      !
278      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
279      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
280      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
281      !        (        ...               )( ...  ) ( ...  )
282      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
283      !
284      !   m is decomposed in the product of an upper and a lower triangular matrix
285      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
286      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
287      !-----------------------------------------------------------------------
288      !
289      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
290         DO jj = 2, jpjm1   
291            DO ji = fs_2, fs_jpim1   ! vector opt.
292               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
293            END DO
294         END DO
295      END DO
296      !
297      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
298         DO ji = fs_2, fs_jpim1   ! vector opt.
299            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
300            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
301               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
302         END DO
303      END DO
304      DO jk = 2, jpkm1
305         DO jj = 2, jpjm1
306            DO ji = fs_2, fs_jpim1
307               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)
308            END DO
309         END DO
310      END DO
311      !
312      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
313         DO ji = fs_2, fs_jpim1   ! vector opt.
314            puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
315         END DO
316      END DO
317      DO jk = jpk-2, 1, -1
318         DO jj = 2, jpjm1
319            DO ji = fs_2, fs_jpim1
320               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
321            END DO
322         END DO
323      END DO
324      !
325      !              !==  Vertical diffusion on v  ==!
326      !
327      !                       !* Matrix construction
328      zdt = r2dt * 0.5
329      IF( ln_zad_Aimp ) THEN   !!
330         SELECT CASE( nldf_dyn )
331         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
332            DO jk = 1, jpkm1
333               DO jj = 2, jpjm1 
334                  DO ji = fs_2, fs_jpim1   ! vector opt.
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  ) + akzv(ji,jj,jk  ) )   &
337                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
338                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
339                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
340                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
341                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
342                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
343                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
344                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
345                  END DO
346               END DO
347            END DO
348         CASE DEFAULT               ! iso-level lateral mixing
349            DO jk = 1, jpkm1
350               DO jj = 2, jpjm1 
351                  DO ji = fs_2, fs_jpim1   ! vector opt.
352                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
353                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
354                     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)
355                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
356                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
357                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
358                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
359                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
360                  END DO
361               END DO
362            END DO
363         END SELECT
364         DO jj = 2, jpjm1     !* Surface boundary conditions
365            DO ji = fs_2, fs_jpim1   ! vector opt.
366               zwi(ji,jj,1) = 0._wp
367               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
368               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
369               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) )
370               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
371               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
372            END DO
373         END DO
374      ELSE
375         SELECT CASE( nldf_dyn )
376         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
377            DO jk = 1, jpkm1
378               DO jj = 2, jpjm1   
379                  DO ji = fs_2, fs_jpim1   ! vector opt.
380                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
381                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
382                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
383                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
384                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
385                     zwi(ji,jj,jk) = zzwi
386                     zws(ji,jj,jk) = zzws
387                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
388                  END DO
389               END DO
390            END DO
391         CASE DEFAULT               ! iso-level lateral mixing
392            DO jk = 1, jpkm1
393               DO jj = 2, jpjm1   
394                  DO ji = fs_2, fs_jpim1   ! vector opt.
395                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
396                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
397                     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)
398                     zwi(ji,jj,jk) = zzwi
399                     zws(ji,jj,jk) = zzws
400                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
401                  END DO
402               END DO
403            END DO
404         END SELECT
405         DO jj = 2, jpjm1        !* Surface boundary conditions
406            DO ji = fs_2, fs_jpim1   ! vector opt.
407               zwi(ji,jj,1) = 0._wp
408               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
409            END DO
410         END DO
411      ENDIF
412      !
413      !              !==  Apply semi-implicit top/bottom friction  ==!
414      !
415      !     Only needed for semi-implicit bottom friction setup. The explicit
416      !     bottom friction has been included in "u(v)a" which act as the R.H.S
417      !     column vector of the tri-diagonal matrix equation
418      !
419      IF( ln_drgimp ) THEN
420         DO jj = 2, jpjm1
421            DO ji = 2, jpim1
422               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
423               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
424               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
425            END DO
426         END DO
427         IF ( ln_isfcav ) THEN
428            DO jj = 2, jpjm1
429               DO ji = 2, jpim1
430                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
431                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
432                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
433               END DO
434            END DO
435         ENDIF
436      ENDIF
437
438      ! Matrix inversion
439      !-----------------------------------------------------------------------
440      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
441      !
442      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
443      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
444      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
445      !        (        ...               )( ...  ) ( ...  )
446      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
447      !
448      !   m is decomposed in the product of an upper and lower triangular matrix
449      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
450      !   The solution (after velocity) is in 2d array pvv(:,:,:,Kaa)
451      !-----------------------------------------------------------------------
452      !
453      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
454         DO jj = 2, jpjm1   
455            DO ji = fs_2, fs_jpim1   ! vector opt.
456               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
457            END DO
458         END DO
459      END DO
460      !
461      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
462         DO ji = fs_2, fs_jpim1   ! vector opt.         
463            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
464            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
465               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
466         END DO
467      END DO
468      DO jk = 2, jpkm1
469         DO jj = 2, jpjm1
470            DO ji = fs_2, fs_jpim1   ! vector opt.
471               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)
472            END DO
473         END DO
474      END DO
475      !
476      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
477         DO ji = fs_2, fs_jpim1   ! vector opt.
478            pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
479         END DO
480      END DO
481      DO jk = jpk-2, 1, -1
482         DO jj = 2, jpjm1
483            DO ji = fs_2, fs_jpim1
484               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
485            END DO
486         END DO
487      END DO
488      !
489      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
490         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
491         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
492         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
493         DEALLOCATE( ztrdu, ztrdv ) 
494      ENDIF
495      !                                          ! print mean trends (used for debugging)
496      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
497         &                       tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
498         !
499      IF( ln_timing )   CALL timing_stop('dyn_zdf')
500      !
501   END SUBROUTINE dyn_zdf
502
503   !!==============================================================================
504END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.