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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90 @ 9250

Last change on this file since 9250 was 9250, checked in by gm, 7 years ago

dev_merge_2017 : HPC09_ZDF: drag finalization: move dyn_bfr routine in zdfdrg (renamed zdf_drg_exp) and delete dynbfr.F90 module - result unchanged
  • Property svn:keywords set to Id
File size: 20.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    ,ONLY: nldf, np_lap_i   ! dynamics: type of lateral mixing
22   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
23   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
24   USE trd_oce        ! trends: ocean variables
25   USE trddyn         ! trend manager: dynamics
26   !
27   USE in_out_manager ! I/O manager
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/OPA 4.0 , NEMO Consortium (2017)
43   !! $Id$
44   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
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), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
75      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
76      !!---------------------------------------------------------------------
77      !
78      IF( ln_timing )   CALL timing_start('dyn_zdf')
79      !
80      IF( kt == nit000 ) THEN       !* initialization
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
83         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
84         !
85         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
86         ELSE                   ;    r_vvl = 1._wp
87         ENDIF
88      ENDIF
89      !                             !* set time step
90      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
91      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
92      ENDIF
93      !
94      !                             !* explicit top/bottom drag case
95      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
96      !
97      !
98      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
99         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
100         ztrdu(:,:,:) = ua(:,:,:)
101         ztrdv(:,:,:) = va(:,:,:)
102      ENDIF
103      !
104      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
105      !
106      !                    ! time stepping except vertical diffusion
107      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
108         DO jk = 1, jpkm1
109            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
110            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
111         END DO
112      ELSE                                      ! applied on thickness weighted velocity
113         DO jk = 1, jpkm1
114            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
115               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
116            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
117               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
118         END DO
119      ENDIF
120      !                    ! add top/bottom friction
121      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
122      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
123      !                not lead to the effective stress seen over the whole barotropic loop.
124      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
125      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
126         DO jk = 1, jpkm1        ! remove barotropic velocities
127            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
128            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
129         END DO
130         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
131            DO ji = fs_2, fs_jpim1   ! vector opt.
132               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
133               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
134               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
135               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
136               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
137               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
138            END DO
139         END DO
140         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
141            DO jj = 2, jpjm1       
142               DO ji = fs_2, fs_jpim1   ! vector opt.
143                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
144                  ikv = mikv(ji,jj)         ! (first wet 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_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
148                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
149               END DO
150            END DO
151         END IF
152      ENDIF
153      !
154      !              !==  Vertical diffusion on u  ==!
155      !
156      !                    !* Matrix construction
157      zdt = r2dt * 0.5
158      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu)
159         DO jk = 1, jpkm1
160            DO jj = 2, jpjm1 
161               DO ji = fs_2, fs_jpim1   ! vector opt.
162                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
163                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
164                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
165                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
166                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
167                  zwi(ji,jj,jk) = zzwi
168                  zws(ji,jj,jk) = zzws
169                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
170               END DO
171            END DO
172         END DO
173      ELSE                          ! standard case
174         DO jk = 1, jpkm1
175            DO jj = 2, jpjm1 
176               DO ji = fs_2, fs_jpim1   ! vector opt.
177                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
178                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
179                  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)
180                  zwi(ji,jj,jk) = zzwi
181                  zws(ji,jj,jk) = zzws
182                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
183               END DO
184            END DO
185         END DO
186      ENDIF
187      !
188      DO jj = 2, jpjm1     !* Surface boundary conditions
189         DO ji = fs_2, fs_jpim1   ! vector opt.
190            zwi(ji,jj,1) = 0._wp
191            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
192         END DO
193      END DO
194      !
195      !              !==  Apply semi-implicit bottom friction  ==!
196      !
197      !     Only needed for semi-implicit bottom friction setup. The explicit
198      !     bottom friction has been included in "u(v)a" which act as the R.H.S
199      !     column vector of the tri-diagonal matrix equation
200      !
201      IF ( ln_drgimp ) THEN      ! implicit bottom friction
202         DO jj = 2, jpjm1
203            DO ji = 2, jpim1
204               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
205               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
206               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
207            END DO
208         END DO
209         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
210            DO jj = 2, jpjm1
211               DO ji = 2, jpim1
212                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
213                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
214                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
215                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
216               END DO
217            END DO
218         END IF
219      ENDIF
220      !
221      ! Matrix inversion starting from the first level
222      !-----------------------------------------------------------------------
223      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
224      !
225      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
226      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
227      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
228      !        (        ...               )( ...  ) ( ...  )
229      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
230      !
231      !   m is decomposed in the product of an upper and a lower triangular matrix
232      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
233      !   The solution (the after velocity) is in ua
234      !-----------------------------------------------------------------------
235      !
236      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
237         DO jj = 2, jpjm1   
238            DO ji = fs_2, fs_jpim1   ! vector opt.
239               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
240            END DO
241         END DO
242      END DO
243      !
244      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
245         DO ji = fs_2, fs_jpim1   ! vector opt.
246            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
247            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
248               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
249         END DO
250      END DO
251      DO jk = 2, jpkm1
252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1
254               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
255            END DO
256         END DO
257      END DO
258      !
259      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
260         DO ji = fs_2, fs_jpim1   ! vector opt.
261            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
262         END DO
263      END DO
264      DO jk = jpk-2, 1, -1
265         DO jj = 2, jpjm1
266            DO ji = fs_2, fs_jpim1
267               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
268            END DO
269         END DO
270      END DO
271      !
272      !              !==  Vertical diffusion on v  ==!
273      !
274      !                       !* Matrix construction
275      zdt = r2dt * 0.5
276      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu)
277         DO jk = 1, jpkm1
278            DO jj = 2, jpjm1   
279               DO ji = fs_2, fs_jpim1   ! vector opt.
280                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
281                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
282                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
283                  zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
284                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
285                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
286                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
287                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
288               END DO
289            END DO
290         END DO
291      ELSE                          ! standard case
292         DO jk = 1, jpkm1
293            DO jj = 2, jpjm1   
294               DO ji = fs_2, fs_jpim1   ! vector opt.
295                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
296                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
297                  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)
298                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
299                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
300                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
301               END DO
302            END DO
303         END DO
304      ENDIF
305      !
306      DO jj = 2, jpjm1        !* Surface boundary conditions
307         DO ji = fs_2, fs_jpim1   ! vector opt.
308            zwi(ji,jj,1) = 0._wp
309            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
310         END DO
311      END DO
312      !              !==  Apply semi-implicit top/bottom friction  ==!
313      !
314      !     Only needed for semi-implicit bottom friction setup. The explicit
315      !     bottom friction has been included in "u(v)a" which act as the R.H.S
316      !     column vector of the tri-diagonal matrix equation
317      !
318      IF( ln_drgimp ) THEN
319         DO jj = 2, jpjm1
320            DO ji = 2, jpim1
321               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
322               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
323               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
324            END DO
325         END DO
326         IF ( ln_isfcav ) THEN
327            DO jj = 2, jpjm1
328               DO ji = 2, jpim1
329                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
330                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
331                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
332               END DO
333            END DO
334         ENDIF
335      ENDIF
336
337      ! Matrix inversion
338      !-----------------------------------------------------------------------
339      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
340      !
341      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
342      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
343      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
344      !        (        ...               )( ...  ) ( ...  )
345      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
346      !
347      !   m is decomposed in the product of an upper and lower triangular matrix
348      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
349      !   The solution (after velocity) is in 2d array va
350      !-----------------------------------------------------------------------
351      !
352      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
353         DO jj = 2, jpjm1   
354            DO ji = fs_2, fs_jpim1   ! vector opt.
355               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
356            END DO
357         END DO
358      END DO
359      !
360      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
361         DO ji = fs_2, fs_jpim1   ! vector opt.         
362            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
363            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
364               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
365         END DO
366      END DO
367      DO jk = 2, jpkm1
368         DO jj = 2, jpjm1
369            DO ji = fs_2, fs_jpim1   ! vector opt.
370               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
371            END DO
372         END DO
373      END DO
374      !
375      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
376         DO ji = fs_2, fs_jpim1   ! vector opt.
377            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
378         END DO
379      END DO
380      DO jk = jpk-2, 1, -1
381         DO jj = 2, jpjm1
382            DO ji = fs_2, fs_jpim1
383               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
384            END DO
385         END DO
386      END DO
387      !
388      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
389         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
390         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
391         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
392         DEALLOCATE( ztrdu, ztrdv ) 
393      ENDIF
394      !                                          ! print mean trends (used for debugging)
395      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
396         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
397         !
398      IF( ln_timing )   CALL timing_stop('dyn_zdf')
399      !
400   END SUBROUTINE dyn_zdf
401
402   !!==============================================================================
403END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.