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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5836

Last change on this file since 5836 was 5836, checked in by cetlod, 9 years ago

merge the simplification branch onto the trunk, see ticket #1612
  • Property svn:keywords set to Id
File size: 17.3 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 lk_vvl=T, =0 otherwise
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dyn_zdf_imp( kt, p2dt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_zdf_imp  ***
50      !!                   
51      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
52      !!      and the surface forcing, and add it to the general trend of
53      !!      the momentum equations.
54      !!
55      !! ** Method  :   The vertical momentum mixing trend is given by :
56      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
57      !!      backward time stepping
58      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
59      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
60      !!      Add this trend to the general trend ua :
61      !!         ua = ua + dz( avmu dz(u) )
62      !!
63      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
64      !!---------------------------------------------------------------------
65      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
66      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
67      !!
68      INTEGER  ::   ji, jj, jk   ! dummy loop indices
69      INTEGER  ::   ikbu, ikbv   ! local integers
70      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
71      REAL(wp) ::   ze3ua, ze3va
72      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
76      !
77      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
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( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
85         ELSE                ;    r_vvl = 0._wp       
86         ENDIF
87      ENDIF
88
89      ! 0. Local constant initialization
90      ! --------------------------------
91      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
92
93      ! 1. Apply semi-implicit bottom friction
94      ! --------------------------------------
95      ! Only needed for semi-implicit bottom friction setup. The explicit
96      ! bottom friction has been included in "u(v)a" which act as the R.H.S
97      ! column vector of the tri-diagonal matrix equation
98      !
99
100      IF( ln_bfrimp ) THEN
101         DO jj = 2, jpjm1
102            DO ji = 2, jpim1
103               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
104               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
105               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
106               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
107            END DO
108         END DO
109         IF ( ln_isfcav ) THEN
110            DO jj = 2, jpjm1
111               DO ji = 2, jpim1
112                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
113                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
114                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
115                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
116               END DO
117            END DO
118         END IF
119      ENDIF
120
121#if defined key_dynspg_ts
122      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
123         DO jk = 1, jpkm1
124            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
125            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
126         END DO
127      ELSE                                            ! applied on thickness weighted velocity
128         DO jk = 1, jpkm1
129            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
130               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
131               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
132            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
133               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
134               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
135         END DO
136      ENDIF
137
138      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
139         ! remove barotropic velocities:
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
142            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
143         END DO
144         ! Add bottom/top stress due to barotropic component only:
145         DO jj = 2, jpjm1       
146            DO ji = fs_2, fs_jpim1   ! vector opt.
147               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
148               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
149               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
150               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
151               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
152               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
153            END DO
154         END DO
155         IF ( ln_isfcav ) THEN
156            DO jj = 2, jpjm1       
157               DO ji = fs_2, fs_jpim1   ! vector opt.
158                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
159                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
160                  ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
161                  ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
162                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
163                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
164               END DO
165            END DO
166         END IF
167      ENDIF
168#endif
169
170      ! 2. Vertical diffusion on u
171      ! ---------------------------
172      ! Matrix and second member construction
173      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
174      ! non zero value at the ocean bottom depending on the bottom friction used.
175      !
176      DO jk = 1, jpkm1        ! Matrix
177         DO jj = 2, jpjm1 
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
180               zcoef = - p2dt / ze3ua     
181               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
182               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
183               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
184               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
185               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
186            END DO
187         END DO
188      END DO
189      DO jj = 2, jpjm1        ! Surface boundary conditions
190         DO ji = fs_2, fs_jpim1   ! vector opt.
191            zwi(ji,jj,1) = 0._wp
192            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
193         END DO
194      END DO
195
196      ! Matrix inversion starting from the first level
197      !-----------------------------------------------------------------------
198      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
199      !
200      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
201      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
202      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
203      !        (        ...               )( ...  ) ( ...  )
204      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
205      !
206      !   m is decomposed in the product of an upper and a lower triangular matrix
207      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
208      !   The solution (the after velocity) is in ua
209      !-----------------------------------------------------------------------
210      !
211      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
212         DO jj = 2, jpjm1   
213            DO ji = fs_2, fs_jpim1   ! vector opt.
214               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
215            END DO
216         END DO
217      END DO
218      !
219      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
220         DO ji = fs_2, fs_jpim1   ! vector opt.
221#if defined key_dynspg_ts
222            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
223            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
224               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
225#else
226            ua(ji,jj,1) = ub(ji,jj,1) &
227               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
228               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) ) 
229#endif
230         END DO
231      END DO
232      DO jk = 2, jpkm1
233         DO jj = 2, jpjm1
234            DO ji = fs_2, fs_jpim1
235#if defined key_dynspg_ts
236               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
237#else
238               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
239#endif
240               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
241            END DO
242         END DO
243      END DO
244      !
245      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
246         DO ji = fs_2, fs_jpim1   ! vector opt.
247            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
248         END DO
249      END DO
250      DO jk = jpk-2, 1, -1
251         DO jj = 2, jpjm1
252            DO ji = fs_2, fs_jpim1
253               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
254            END DO
255         END DO
256      END DO
257
258#if ! defined key_dynspg_ts
259      ! Normalization to obtain the general momentum trend ua
260      DO jk = 1, jpkm1
261         DO jj = 2, jpjm1   
262            DO ji = fs_2, fs_jpim1   ! vector opt.
263               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
264            END DO
265         END DO
266      END DO
267#endif
268
269      ! 3. Vertical diffusion on v
270      ! ---------------------------
271      ! Matrix and second member construction
272      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
273      ! non zero value at the ocean bottom depending on the bottom friction used
274      !
275      DO jk = 1, jpkm1        ! Matrix
276         DO jj = 2, jpjm1   
277            DO ji = fs_2, fs_jpim1   ! vector opt.
278               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
279               zcoef = - p2dt / ze3va
280               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
281               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
282               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
283               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
284               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
285            END DO
286         END DO
287      END DO
288      DO jj = 2, jpjm1        ! Surface boundary conditions
289         DO ji = fs_2, fs_jpim1   ! vector opt.
290            zwi(ji,jj,1) = 0._wp
291            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
292         END DO
293      END DO
294
295      ! Matrix inversion
296      !-----------------------------------------------------------------------
297      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
298      !
299      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
300      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
301      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
302      !        (        ...               )( ...  ) ( ...  )
303      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
304      !
305      !   m is decomposed in the product of an upper and lower triangular matrix
306      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
307      !   The solution (after velocity) is in 2d array va
308      !-----------------------------------------------------------------------
309      !
310      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
311         DO jj = 2, jpjm1   
312            DO ji = fs_2, fs_jpim1   ! vector opt.
313               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
314            END DO
315         END DO
316      END DO
317      !
318      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
319         DO ji = fs_2, fs_jpim1   ! vector opt.
320#if defined key_dynspg_ts           
321            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
322            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
323               &                                      / ( ze3va * rau0 ) 
324#else
325            va(ji,jj,1) = vb(ji,jj,1) &
326               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
327               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
328#endif
329         END DO
330      END DO
331      DO jk = 2, jpkm1
332         DO jj = 2, jpjm1
333            DO ji = fs_2, fs_jpim1   ! vector opt.
334#if defined key_dynspg_ts
335               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
336#else
337               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
338#endif
339               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
340            END DO
341         END DO
342      END DO
343      !
344      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
345         DO ji = fs_2, fs_jpim1   ! vector opt.
346            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
347         END DO
348      END DO
349      DO jk = jpk-2, 1, -1
350         DO jj = 2, jpjm1
351            DO ji = fs_2, fs_jpim1
352               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
353            END DO
354         END DO
355      END DO
356
357      ! Normalization to obtain the general momentum trend va
358#if ! defined key_dynspg_ts
359      DO jk = 1, jpkm1
360         DO jj = 2, jpjm1   
361            DO ji = fs_2, fs_jpim1   ! vector opt.
362               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
363            END DO
364         END DO
365      END DO
366#endif
367
368      ! J. Chanut: Lines below are useless ?
369      !! restore bottom layer avmu(v)
370      IF( ln_bfrimp ) THEN
371        DO jj = 2, jpjm1
372           DO ji = 2, jpim1
373              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
374              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
375              avmu(ji,jj,ikbu+1) = 0.e0
376              avmv(ji,jj,ikbv+1) = 0.e0
377           END DO
378        END DO
379        IF (ln_isfcav) THEN
380           DO jj = 2, jpjm1
381              DO ji = 2, jpim1
382                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
383                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
384                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
385                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
386              END DO
387           END DO
388        END IF
389      ENDIF
390      !
391      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
392      !
393      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
394      !
395   END SUBROUTINE dyn_zdf_imp
396
397   !!==============================================================================
398END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.