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_r11943_MERGE_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynzdf.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 23.3 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#  include "do_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, Kbb, Kmm, Krhs, puu, pvv, Kaa )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dyn_zdf  ***
51      !!
52      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
53      !!              together with the Leap-Frog time stepping using an
54      !!              implicit scheme.
55      !!
56      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
57      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
58      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise
59      !!               - update the after velocity with the implicit vertical mixing.
60      !!      This requires to solver the following system:
61      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ]
62      !!      with the following surface/top/bottom boundary condition:
63      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
64      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
65      !!
66      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity
67      !!---------------------------------------------------------------------
68      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
69      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
70      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
71      !
72      INTEGER  ::   ji, jj, jk         ! dummy loop indices
73      INTEGER  ::   iku, ikv           ! local integers
74      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
75      REAL(wp) ::   zzws, ze3va        !   -      -
76      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
77      REAL(wp) ::   zWu , zWv          !   -      -
78      REAL(wp) ::   zWui, zWvi         !   -      -
79      REAL(wp) ::   zWus, zWvs         !   -      -
80      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
81      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
82      !!---------------------------------------------------------------------
83      !
84      IF( ln_timing )   CALL timing_start('dyn_zdf')
85      !
86      IF( kt == nit000 ) THEN       !* initialization
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
89         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
90         !
91         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
92         ELSE                   ;    r_vvl = 1._wp
93         ENDIF
94      ENDIF
95      !                             !* set time step
96      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
97      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
98      ENDIF
99      !
100      !                             !* explicit top/bottom drag case
101      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
102      !
103      !
104      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
105         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
106         ztrdu(:,:,:) = puu(:,:,:,Krhs)
107         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
108      ENDIF
109      !
110      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
111      !
112      !                    ! time stepping except vertical diffusion
113      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
114         DO jk = 1, jpkm1
115            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)
116            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)
117         END DO
118      ELSE                                      ! applied on thickness weighted velocity
119         DO jk = 1, jpkm1
120            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
121               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
122            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
123               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * 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(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
131      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
132         DO jk = 1, jpkm1        ! remove barotropic velocities
133            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)
134            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)
135         END DO
136         DO_2D_00_00
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_2D
144         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
145            DO_2D_00_00
146               iku = miku(ji,jj)         ! top ocean level at u- and v-points
147               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
148               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
149               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
150               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
151               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
152            END_2D
153         END IF
154      ENDIF
155      !
156      !              !==  Vertical diffusion on u  ==!
157      !
158      !                    !* Matrix construction
159      zdt = r2dt * 0.5
160      IF( ln_zad_Aimp ) THEN   !!
161         SELECT CASE( nldf_dyn )
162         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
163            DO_3D_00_00( 1, jpkm1 )
164               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
165               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
166                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
167               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
168                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
169               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
170               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
171               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
172               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
173               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
174            END_3D
175         CASE DEFAULT               ! iso-level lateral mixing
176            DO_3D_00_00( 1, jpkm1 )
177               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
178               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
179               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)
180               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
181               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
182               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
183               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
184               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
185            END_3D
186         END SELECT
187         DO_2D_00_00
188            zwi(ji,jj,1) = 0._wp
189            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
190            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
191            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
192            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
193            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
194         END_2D
195      ELSE
196         SELECT CASE( nldf_dyn )
197         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
198            DO_3D_00_00( 1, jpkm1 )
199               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
200               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
201                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
202               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
203                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
204               zwi(ji,jj,jk) = zzwi
205               zws(ji,jj,jk) = zzws
206               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
207            END_3D
208         CASE DEFAULT               ! iso-level lateral mixing
209            DO_3D_00_00( 1, jpkm1 )
210               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
211               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
212               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)
213               zwi(ji,jj,jk) = zzwi
214               zws(ji,jj,jk) = zzws
215               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
216            END_3D
217         END SELECT
218         DO_2D_00_00
219            zwi(ji,jj,1) = 0._wp
220            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
221         END_2D
222      ENDIF
223      !
224      !
225      !              !==  Apply semi-implicit bottom friction  ==!
226      !
227      !     Only needed for semi-implicit bottom friction setup. The explicit
228      !     bottom friction has been included in "u(v)a" which act as the R.H.S
229      !     column vector of the tri-diagonal matrix equation
230      !
231      IF ( ln_drgimp ) THEN      ! implicit bottom friction
232         DO_2D_00_00
233            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
234            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
235            zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
236         END_2D
237         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
238            DO_2D_00_00
239               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
240               iku = miku(ji,jj)       ! ocean top level at u- and v-points
241               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
242               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
243            END_2D
244         END IF
245      ENDIF
246      !
247      ! Matrix inversion starting from the first level
248      !-----------------------------------------------------------------------
249      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
250      !
251      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
252      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
253      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
254      !        (        ...               )( ...  ) ( ...  )
255      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
256      !
257      !   m is decomposed in the product of an upper and a lower triangular matrix
258      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
259      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
260      !-----------------------------------------------------------------------
261      !
262      DO_3D_00_00( 2, jpkm1 )
263         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
264      END_3D
265      !
266      DO_2D_00_00
267         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
268         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
269            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
270      END_2D
271      DO_3D_00_00( 2, jpkm1 )
272         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)
273      END_3D
274      !
275      DO_2D_00_00
276         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
277      END_2D
278      DO_3DS_00_00( jpk-2, 1, -1 )
279         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
280      END_3D
281      !
282      !              !==  Vertical diffusion on v  ==!
283      !
284      !                       !* Matrix construction
285      zdt = r2dt * 0.5
286      IF( ln_zad_Aimp ) THEN   !!
287         SELECT CASE( nldf_dyn )
288         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
289            DO_3D_00_00( 1, jpkm1 )
290               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
291               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
292                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
293               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
294                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
295               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
296               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
297               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
298               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
299               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
300            END_3D
301         CASE DEFAULT               ! iso-level lateral mixing
302            DO_3D_00_00( 1, jpkm1 )
303               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
304               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
305               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)
306               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
307               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
308               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
309               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
310               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
311            END_3D
312         END SELECT
313         DO_2D_00_00
314            zwi(ji,jj,1) = 0._wp
315            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
316            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
317            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
318            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
319            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
320         END_2D
321      ELSE
322         SELECT CASE( nldf_dyn )
323         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
324            DO_3D_00_00( 1, jpkm1 )
325               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
326               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
327                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
328               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
329                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
330               zwi(ji,jj,jk) = zzwi
331               zws(ji,jj,jk) = zzws
332               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
333            END_3D
334         CASE DEFAULT               ! iso-level lateral mixing
335            DO_3D_00_00( 1, jpkm1 )
336               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
337               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
338               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)
339               zwi(ji,jj,jk) = zzwi
340               zws(ji,jj,jk) = zzws
341               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
342            END_3D
343         END SELECT
344         DO_2D_00_00
345            zwi(ji,jj,1) = 0._wp
346            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
347         END_2D
348      ENDIF
349      !
350      !              !==  Apply semi-implicit top/bottom friction  ==!
351      !
352      !     Only needed for semi-implicit bottom friction setup. The explicit
353      !     bottom friction has been included in "u(v)a" which act as the R.H.S
354      !     column vector of the tri-diagonal matrix equation
355      !
356      IF( ln_drgimp ) THEN
357         DO_2D_00_00
358            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
359            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
360            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
361         END_2D
362         IF ( ln_isfcav ) THEN
363            DO_2D_00_00
364               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
365               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
366               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
367            END_2D
368         ENDIF
369      ENDIF
370
371      ! Matrix inversion
372      !-----------------------------------------------------------------------
373      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
374      !
375      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
376      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
377      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
378      !        (        ...               )( ...  ) ( ...  )
379      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
380      !
381      !   m is decomposed in the product of an upper and lower triangular matrix
382      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
383      !   The solution (after velocity) is in 2d array va
384      !-----------------------------------------------------------------------
385      !
386      DO_3D_00_00( 2, jpkm1 )
387         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
388      END_3D
389      !
390      DO_2D_00_00
391         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
392         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
393            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
394      END_2D
395      DO_3D_00_00( 2, jpkm1 )
396         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)
397      END_3D
398      !
399      DO_2D_00_00
400         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
401      END_2D
402      DO_3DS_00_00( jpk-2, 1, -1 )
403         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
404      END_3D
405      !
406      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
407         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
408         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
409         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
410         DEALLOCATE( ztrdu, ztrdv ) 
411      ENDIF
412      !                                          ! print mean trends (used for debugging)
413      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
414         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
415         !
416      IF( ln_timing )   CALL timing_stop('dyn_zdf')
417      !
418   END SUBROUTINE dyn_zdf
419
420   !!==============================================================================
421END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.