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_imp.F90 in branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5866

Last change on this file since 5866 was 5866, checked in by gm, 8 years ago

#1613: vvl by default: add ln_linssh and remove key_vvl
  • Property svn:keywords set to Id
File size: 17.1 KB
Line 
1MODULE dynzdf_imp
2   !!======================================================================
3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
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  : update the momentum trend with the vertical diffusion using a implicit time-stepping
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE domvvl          ! variable volume
19   USE sbc_oce         ! surface boundary condition: ocean
20   USE zdf_oce         ! ocean vertical physics
21   USE phycst          ! physical constants
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! MPP library
24   USE zdfbfr          ! Bottom friction setup
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27   USE dynadv          ! dynamics: vector invariant versus flux form
28   USE dynspg_oce, ONLY: lk_dynspg_ts
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   dyn_zdf_imp   ! called by step.F90
34
35   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if ln_linssh=F, =0 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      !!      and the surface forcing, and add it to the general trend of
52      !!      the momentum equations.
53      !!
54      !! ** Method  :   The vertical momentum mixing trend is given by :
55      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
56      !!      backward time stepping
57      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
58      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
59      !!      Add this trend to the general trend ua :
60      !!         ua = ua + dz( avmu dz(u) )
61      !!
62      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
63      !!---------------------------------------------------------------------
64      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
65      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
66      !!
67      INTEGER  ::   ji, jj, jk   ! dummy loop indices
68      INTEGER  ::   ikbu, ikbv   ! local integers
69      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
70      REAL(wp) ::   ze3ua, ze3va
71      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
72      !!----------------------------------------------------------------------
73      !
74      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
75      !
76      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
77      !
78      IF( kt == nit000 ) THEN
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
81         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
82         !
83         IF( .NOT.ln_linssh ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
84         ELSE                        ;    r_vvl = 0._wp       
85         ENDIF
86      ENDIF
87
88      ! 0. Local constant initialization
89      ! --------------------------------
90      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
91
92      ! 1. Apply semi-implicit bottom friction
93      ! --------------------------------------
94      ! Only needed for semi-implicit bottom friction setup. The explicit
95      ! bottom friction has been included in "u(v)a" which act as the R.H.S
96      ! column vector of the tri-diagonal matrix equation
97      !
98
99      IF( ln_bfrimp ) THEN
100         DO jj = 2, jpjm1
101            DO ji = 2, jpim1
102               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
103               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
104               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1)
105               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1)
106            END DO
107         END DO
108         IF ( ln_isfcav ) THEN
109            DO jj = 2, jpjm1
110               DO ji = 2, jpim1
111                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
112                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
113                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu)
114                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv)
115               END DO
116            END DO
117         END IF
118      ENDIF
119
120#if defined key_dynspg_ts
121      IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity
122         DO jk = 1, jpkm1
123            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
124            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
125         END DO
126      ELSE                                         ! applied on thickness weighted velocity
127         DO jk = 1, jpkm1
128            ua(:,:,jk) = (          ub(:,:,jk) * e3u_b(:,:,jk)    &
129               &           + p2dt * ua(:,:,jk) * e3u_n(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
130            va(:,:,jk) = (          vb(:,:,jk) * e3v_b(:,:,jk)    &
131               &           + p2dt * va(:,:,jk) * e3v_n(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
132         END DO
133      ENDIF
134
135      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
136         ! remove barotropic velocities:
137         DO jk = 1, jpkm1
138            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
139            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
140         END DO
141         ! Add bottom/top stress due to barotropic component only:
142         DO jj = 2, jpjm1       
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
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#endif
166
167      ! 2. Vertical diffusion on u
168      ! ---------------------------
169      ! Matrix and second member construction
170      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
171      ! non zero value at the ocean bottom depending on the bottom friction used.
172      !
173      DO jk = 1, jpkm1        ! Matrix
174         DO jj = 2, jpjm1 
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl   * e3u_a(ji,jj,jk)   ! after scale factor at T-point
177               zcoef = - p2dt / ze3ua     
178               zzwi  = zcoef * avmu(ji,jj,jk  ) / e3uw_n(ji,jj,jk  )
179               zzws  = zcoef * avmu(ji,jj,jk+1) / e3uw_n(ji,jj,jk+1) 
180               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
181               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
182               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
183            END DO
184         END DO
185      END DO
186      DO jj = 2, jpjm1        ! Surface boundary conditions
187         DO ji = fs_2, fs_jpim1   ! vector opt.
188            zwi(ji,jj,1) = 0._wp
189            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
190         END DO
191      END DO
192
193      ! Matrix inversion starting from the first level
194      !-----------------------------------------------------------------------
195      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
196      !
197      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
198      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
199      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
200      !        (        ...               )( ...  ) ( ...  )
201      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
202      !
203      !   m is decomposed in the product of an upper and a lower triangular matrix
204      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
205      !   The solution (the after velocity) is in ua
206      !-----------------------------------------------------------------------
207      !
208      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
209         DO jj = 2, jpjm1   
210            DO ji = fs_2, fs_jpim1   ! vector opt.
211               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
212            END DO
213         END DO
214      END DO
215      !
216      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
217         DO ji = fs_2, fs_jpim1   ! vector opt.
218#if defined key_dynspg_ts
219            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl   * e3u_a(ji,jj,1) 
220            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
221               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
222#else
223            ua(ji,jj,1) = ub(ji,jj,1)  &
224               &        + p2dt *( ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
225               &                                        / ( e3u_n(ji,jj,1) * rau0 ) * umask(ji,jj,1) ) 
226#endif
227         END DO
228      END DO
229      DO jk = 2, jpkm1
230         DO jj = 2, jpjm1
231            DO ji = fs_2, fs_jpim1
232#if defined key_dynspg_ts
233               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
234#else
235               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
236#endif
237               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
238            END DO
239         END DO
240      END DO
241      !
242      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
243         DO ji = fs_2, fs_jpim1   ! vector opt.
244            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
245         END DO
246      END DO
247      DO jk = jpk-2, 1, -1
248         DO jj = 2, jpjm1
249            DO ji = fs_2, fs_jpim1
250               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
251            END DO
252         END DO
253      END DO
254
255#if ! defined key_dynspg_ts
256      ! Normalization to obtain the general momentum trend ua
257      DO jk = 1, jpkm1
258         DO jj = 2, jpjm1   
259            DO ji = fs_2, fs_jpim1   ! vector opt.
260               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
261            END DO
262         END DO
263      END DO
264#endif
265
266      ! 3. Vertical diffusion on v
267      ! ---------------------------
268      ! Matrix and second member construction
269      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
270      ! non zero value at the ocean bottom depending on the bottom friction used
271      !
272      DO jk = 1, jpkm1        ! Matrix
273         DO jj = 2, jpjm1   
274            DO ji = fs_2, fs_jpim1   ! vector opt.
275               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk)  + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
276               zcoef = - p2dt / ze3va
277               zzwi          = zcoef * avmv (ji,jj,jk  ) / e3vw_n(ji,jj,jk  )
278               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
279               zzws          = zcoef * avmv (ji,jj,jk+1) / e3vw_n(ji,jj,jk+1)
280               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
281               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
282            END DO
283         END DO
284      END DO
285      DO jj = 2, jpjm1        ! Surface boundary conditions
286         DO ji = fs_2, fs_jpim1   ! vector opt.
287            zwi(ji,jj,1) = 0._wp
288            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
289         END DO
290      END DO
291
292      ! Matrix inversion
293      !-----------------------------------------------------------------------
294      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
295      !
296      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
297      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
298      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
299      !        (        ...               )( ...  ) ( ...  )
300      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
301      !
302      !   m is decomposed in the product of an upper and lower triangular matrix
303      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
304      !   The solution (after velocity) is in 2d array va
305      !-----------------------------------------------------------------------
306      !
307      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
308         DO jj = 2, jpjm1   
309            DO ji = fs_2, fs_jpim1   ! vector opt.
310               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
311            END DO
312         END DO
313      END DO
314      !
315      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
316         DO ji = fs_2, fs_jpim1   ! vector opt.
317#if defined key_dynspg_ts           
318            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl   * e3v_a(ji,jj,1) 
319            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
320               &                                      / ( ze3va * rau0 ) 
321#else
322            va(ji,jj,1) = vb(ji,jj,1)   &
323               &        + p2dt *( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
324               &                                       / ( e3v_n(ji,jj,1) * rau0 )         )
325#endif
326         END DO
327      END DO
328      DO jk = 2, jpkm1
329         DO jj = 2, jpjm1
330            DO ji = fs_2, fs_jpim1   ! vector opt.
331#if defined key_dynspg_ts
332               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
333#else
334               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
335#endif
336               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
337            END DO
338         END DO
339      END DO
340      !
341      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
342         DO ji = fs_2, fs_jpim1   ! vector opt.
343            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
344         END DO
345      END DO
346      DO jk = jpk-2, 1, -1
347         DO jj = 2, jpjm1
348            DO ji = fs_2, fs_jpim1
349               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
350            END DO
351         END DO
352      END DO
353
354      ! Normalization to obtain the general momentum trend va
355#if ! defined key_dynspg_ts
356      DO jk = 1, jpkm1
357         DO jj = 2, jpjm1   
358            DO ji = fs_2, fs_jpim1   ! vector opt.
359               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
360            END DO
361         END DO
362      END DO
363#endif
364
365      ! J. Chanut: Lines below are useless ?
366      !! restore bottom layer avmu(v)
367      IF( ln_bfrimp ) THEN
368        DO jj = 2, jpjm1
369           DO ji = 2, jpim1
370              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
371              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
372              avmu(ji,jj,ikbu+1) = 0.e0
373              avmv(ji,jj,ikbv+1) = 0.e0
374           END DO
375        END DO
376        IF (ln_isfcav) THEN
377           DO jj = 2, jpjm1
378              DO ji = 2, jpim1
379                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
380                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
381                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
382                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
383              END DO
384           END DO
385        END IF
386      ENDIF
387      !
388      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
389      !
390      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
391      !
392   END SUBROUTINE dyn_zdf_imp
393
394   !!==============================================================================
395END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.