source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 16.0 KB
Line 
1MODULE dynzdf_imp
2   !!======================================================================
3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend, implicit scheme
5   !!======================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
8   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
10   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit scheme
15   !!                   together with the Leap-Frog time integration.
16   !!----------------------------------------------------------------------
17   USE oce            ! ocean dynamics and tracers
18   USE phycst         ! physical constants
19   USE dom_oce        ! ocean space and time domain
20   USE domvvl         ! variable volume
21   USE sbc_oce        ! surface boundary condition: ocean
22   USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form
23   USE zdf_oce        ! ocean vertical physics
24   USE zdfbfr         ! Bottom friction setup
25   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   dyn_zdf_imp   ! called by step.F90
34
35   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dyn_zdf_imp( kt, p2dt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dyn_zdf_imp  ***
49      !!                   
50      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
51      !!              together with the Leap-Frog time stepping using an
52      !!              implicit scheme.
53      !!
54      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
55      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
56      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
57      !!               - update the after velocity with the implicit vertical mixing.
58      !!      This requires to solver the following system:
59      !!         ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ]
60      !!      with the following surface/top/bottom boundary condition:
61      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
62      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F)
63      !!
64      !! ** Action :   (ua,va) after velocity
65      !!---------------------------------------------------------------------
66      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
67      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
68      !
69      INTEGER  ::   ji, jj, jk    ! dummy loop indices
70      INTEGER  ::   ikbu, ikbv    ! local integers
71      REAL(wp) ::   zzwi, ze3ua   ! local scalars
72      REAL(wp) ::   zzws, ze3va   !   -      -
73      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwd, zws
74      !!----------------------------------------------------------------------
75      !
76      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
77      !
78      !
79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
82         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
83         !
84         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
85         ELSE                   ;    r_vvl = 1._wp
86         ENDIF
87      ENDIF
88      !
89      !              !==  Time step dynamics  ==!
90      !
91      IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity
92         DO jk = 1, jpkm1
93            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
94            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
95         END DO
96      ELSE                                         ! applied on thickness weighted velocity
97         DO jk = 1, jpkm1
98            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
99               &          + p2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
100            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
101               &          + p2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
102         END DO
103      ENDIF
104      !
105      !              !==  Apply semi-implicit bottom friction  ==!
106      !
107      ! Only needed for semi-implicit bottom friction setup. The explicit
108      ! bottom friction has been included in "u(v)a" which act as the R.H.S
109      ! column vector of the tri-diagonal matrix equation
110      !
111      IF( ln_bfrimp ) THEN
112         DO jj = 2, jpjm1
113            DO ji = 2, jpim1
114               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
115               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
116               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1)
117               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1)
118            END DO
119         END DO
120         IF ( ln_isfcav ) THEN
121            DO jj = 2, jpjm1
122               DO ji = 2, jpim1
123                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
124                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
125                  IF( ikbu >= 2 )   avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu)
126                  IF( ikbv >= 2 )   avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv)
127               END DO
128            END DO
129         END IF
130      ENDIF
131      !
132      ! With split-explicit free surface, barotropic stress is treated explicitly
133      ! Update velocities at the bottom.
134      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
135      !            not lead to the effective stress seen over the whole barotropic loop.
136      ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
137      IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN
138         DO jk = 1, jpkm1        ! remove barotropic velocities
139            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
140            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
141         END DO
142         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
143            DO ji = fs_2, fs_jpim1   ! vector opt.
144               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
145               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
146               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
147               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
148               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
149               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
150            END DO
151         END DO
152         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
153            DO jj = 2, jpjm1       
154               DO ji = fs_2, fs_jpim1   ! vector opt.
155                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
156                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
157                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
158                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
159                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
160                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
161               END DO
162            END DO
163         END IF
164      ENDIF
165      !
166      !              !==  Vertical diffusion on u  ==!
167      !
168      ! Matrix and second member construction
169      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
170      ! non zero value at the ocean bottom depending on the bottom friction used.
171      !
172      DO jk = 1, jpkm1        ! Matrix
173         DO jj = 2, jpjm1 
174            DO ji = fs_2, fs_jpim1   ! vector opt.
175               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
176               zzwi = - p2dt * avmu(ji,jj,jk  ) / ( ze3ua * e3uw_n(ji,jj,jk  ) )
177               zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) )
178               zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk  )
179               zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1)
180               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
181            END DO
182         END DO
183      END DO
184      DO jj = 2, jpjm1        ! Surface boundary conditions
185         DO ji = fs_2, fs_jpim1   ! vector opt.
186            zwi(ji,jj,1) = 0._wp
187            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
188         END DO
189      END DO
190
191      ! Matrix inversion starting from the first level
192      !-----------------------------------------------------------------------
193      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
194      !
195      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
196      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
197      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
198      !        (        ...               )( ...  ) ( ...  )
199      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
200      !
201      !   m is decomposed in the product of an upper and a lower triangular matrix
202      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
203      !   The solution (the after velocity) is in ua
204      !-----------------------------------------------------------------------
205      !
206      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
207         DO jj = 2, jpjm1   
208            DO ji = fs_2, fs_jpim1   ! vector opt.
209               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
210            END DO
211         END DO
212      END DO
213      !
214      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
215         DO ji = fs_2, fs_jpim1   ! vector opt.
216            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
217            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
218               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
219         END DO
220      END DO
221      DO jk = 2, jpkm1
222         DO jj = 2, jpjm1
223            DO ji = fs_2, fs_jpim1
224               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
225            END DO
226         END DO
227      END DO
228      !
229      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
230         DO ji = fs_2, fs_jpim1   ! vector opt.
231            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
232         END DO
233      END DO
234      DO jk = jpk-2, 1, -1
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1
237               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
238            END DO
239         END DO
240      END DO
241      !
242      !              !==  Vertical diffusion on v  ==!
243      !
244      ! Matrix and second member construction
245      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
246      ! non zero value at the ocean bottom depending on the bottom friction used
247      !
248      DO jk = 1, jpkm1        ! Matrix
249         DO jj = 2, jpjm1   
250            DO ji = fs_2, fs_jpim1   ! vector opt.
251               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
252               zzwi = - p2dt * avmv (ji,jj,jk  ) / ( ze3va * e3vw_n(ji,jj,jk  ) )
253               zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) )
254               zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
255               zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
256               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
257            END DO
258         END DO
259      END DO
260      DO jj = 2, jpjm1        ! Surface boundary conditions
261         DO ji = fs_2, fs_jpim1   ! vector opt.
262            zwi(ji,jj,1) = 0._wp
263            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
264         END DO
265      END DO
266
267      ! Matrix inversion
268      !-----------------------------------------------------------------------
269      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
270      !
271      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
272      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
273      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
274      !        (        ...               )( ...  ) ( ...  )
275      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
276      !
277      !   m is decomposed in the product of an upper and lower triangular matrix
278      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
279      !   The solution (after velocity) is in 2d array va
280      !-----------------------------------------------------------------------
281      !
282      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
283         DO jj = 2, jpjm1   
284            DO ji = fs_2, fs_jpim1   ! vector opt.
285               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
286            END DO
287         END DO
288      END DO
289      !
290      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
291         DO ji = fs_2, fs_jpim1   ! vector opt.         
292            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
293            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
294               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
295         END DO
296      END DO
297      DO jk = 2, jpkm1
298         DO jj = 2, jpjm1
299            DO ji = fs_2, fs_jpim1   ! vector opt.
300               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
301            END DO
302         END DO
303      END DO
304      !
305      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
306         DO ji = fs_2, fs_jpim1   ! vector opt.
307            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
308         END DO
309      END DO
310      DO jk = jpk-2, 1, -1
311         DO jj = 2, jpjm1
312            DO ji = fs_2, fs_jpim1
313               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
314            END DO
315         END DO
316      END DO
317     
318      ! J. Chanut: Lines below are useless ?
319      !! restore bottom layer avmu(v)
320      !!gm  I almost sure it is !!!!
321      IF( ln_bfrimp ) THEN
322        DO jj = 2, jpjm1
323           DO ji = 2, jpim1
324              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
325              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
326              avmu(ji,jj,ikbu+1) = 0._wp
327              avmv(ji,jj,ikbv+1) = 0._wp
328           END DO
329        END DO
330        IF (ln_isfcav) THEN
331           DO jj = 2, jpjm1
332              DO ji = 2, jpim1
333                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
334                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
335                 IF( ikbu > 1 )   avmu(ji,jj,ikbu) = 0._wp
336                 IF( ikbv > 1 )   avmv(ji,jj,ikbv) = 0._wp
337              END DO
338           END DO
339        ENDIF
340      ENDIF
341      !
342      !
343      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp')
344      !
345   END SUBROUTINE dyn_zdf_imp
346
347   !!==============================================================================
348END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.