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 @ 10825

Last change on this file since 10825 was 10825, checked in by davestorkey, 5 years ago

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