source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 6079

Last change on this file since 6079 was 6079, checked in by jamesharle, 6 years ago

merge to trunk@5936

  • Property svn:keywords set to Id
File size: 16.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, ONLY: ln_dynadv_vec ! Momentum advection form
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dyn_zdf_imp   ! called by step.F90
33
34   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
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( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
84         ELSE                ;    r_vvl = 0._wp       
85         ENDIF
86      ENDIF
87
88      ! 1. Time step dynamics
89      ! ---------------------
90
91      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) 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) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
99               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
100               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
101            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
102               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
103               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
104         END DO
105      ENDIF
106
107      ! 2. Apply semi-implicit bottom friction
108      ! --------------------------------------
109      ! Only needed for semi-implicit bottom friction setup. The explicit
110      ! bottom friction has been included in "u(v)a" which act as the R.H.S
111      ! column vector of the tri-diagonal matrix equation
112      !
113      IF( ln_bfrimp ) THEN
114         DO jj = 2, jpjm1
115            DO ji = 2, jpim1
116               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
117               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
118               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
119               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
120            END DO
121         END DO
122         IF ( ln_isfcav ) THEN
123            DO jj = 2, jpjm1
124               DO ji = 2, jpim1
125                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
126                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
127                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
128                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
129               END DO
130            END DO
131         END IF
132      ENDIF
133
134      ! With split-explicit free surface, barotropic stress is treated explicitly
135      ! Update velocities at the bottom.
136      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
137      !            not lead to the effective stress seen over the whole barotropic loop.
138      IF ( ln_bfrimp .AND.ln_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
169      ! 3. Vertical diffusion on u
170      ! ---------------------------
171      ! Matrix and second member construction
172      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
173      ! non zero value at the ocean bottom depending on the bottom friction used.
174      !
175      DO jk = 1, jpkm1        ! Matrix
176         DO jj = 2, jpjm1 
177            DO ji = fs_2, fs_jpim1   ! vector opt.
178               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
179               zcoef = - p2dt / ze3ua     
180               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
181               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
182               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
183               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
184               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
185            END DO
186         END DO
187      END DO
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      ! Matrix inversion starting from the first level
196      !-----------------------------------------------------------------------
197      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
198      !
199      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
200      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
201      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
202      !        (        ...               )( ...  ) ( ...  )
203      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
204      !
205      !   m is decomposed in the product of an upper and a lower triangular matrix
206      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
207      !   The solution (the after velocity) is in ua
208      !-----------------------------------------------------------------------
209      !
210      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
211         DO jj = 2, jpjm1   
212            DO ji = fs_2, fs_jpim1   ! vector opt.
213               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
214            END DO
215         END DO
216      END DO
217      !
218      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
219         DO ji = fs_2, fs_jpim1   ! vector opt.
220            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
221            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
222               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
223         END DO
224      END DO
225      DO jk = 2, jpkm1
226         DO jj = 2, jpjm1
227            DO ji = fs_2, fs_jpim1
228               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
229               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
230            END DO
231         END DO
232      END DO
233      !
234      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
235         DO ji = fs_2, fs_jpim1   ! vector opt.
236            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
237         END DO
238      END DO
239      DO jk = jpk-2, 1, -1
240         DO jj = 2, jpjm1
241            DO ji = fs_2, fs_jpim1
242               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
243            END DO
244         END DO
245      END DO
246
247      ! 4. Vertical diffusion on v
248      ! ---------------------------
249      ! Matrix and second member construction
250      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
251      ! non zero value at the ocean bottom depending on the bottom friction used
252      !
253      DO jk = 1, jpkm1        ! Matrix
254         DO jj = 2, jpjm1   
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
257               zcoef = - p2dt / ze3va
258               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
259               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
260               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
261               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
262               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
263            END DO
264         END DO
265      END DO
266      DO jj = 2, jpjm1        ! Surface boundary conditions
267         DO ji = fs_2, fs_jpim1   ! vector opt.
268            zwi(ji,jj,1) = 0._wp
269            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
270         END DO
271      END DO
272
273      ! Matrix inversion
274      !-----------------------------------------------------------------------
275      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
276      !
277      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
278      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
279      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
280      !        (        ...               )( ...  ) ( ...  )
281      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
282      !
283      !   m is decomposed in the product of an upper and lower triangular matrix
284      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
285      !   The solution (after velocity) is in 2d array va
286      !-----------------------------------------------------------------------
287      !
288      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
289         DO jj = 2, jpjm1   
290            DO ji = fs_2, fs_jpim1   ! vector opt.
291               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
292            END DO
293         END DO
294      END DO
295      !
296      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
297         DO ji = fs_2, fs_jpim1   ! vector opt.         
298            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
299            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
300               &                                      / ( ze3va * rau0 ) 
301         END DO
302      END DO
303      DO jk = 2, jpkm1
304         DO jj = 2, jpjm1
305            DO ji = fs_2, fs_jpim1   ! vector opt.
306               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
307               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
308            END DO
309         END DO
310      END DO
311      !
312      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
313         DO ji = fs_2, fs_jpim1   ! vector opt.
314            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
315         END DO
316      END DO
317      DO jk = jpk-2, 1, -1
318         DO jj = 2, jpjm1
319            DO ji = fs_2, fs_jpim1
320               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
321            END DO
322         END DO
323      END DO
324
325      ! J. Chanut: Lines below are useless ?
326      !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces)
327      IF( ln_bfrimp ) THEN
328        DO jj = 2, jpjm1
329           DO ji = 2, jpim1
330              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
331              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
332              avmu(ji,jj,ikbu+1) = 0.e0
333              avmv(ji,jj,ikbv+1) = 0.e0
334           END DO
335        END DO
336        IF (ln_isfcav) THEN
337           DO jj = 2, jpjm1
338              DO ji = 2, jpim1
339                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
340                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
341                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
342                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
343              END DO
344           END DO
345        END IF
346      ENDIF
347      !
348      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
349      !
350      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
351      !
352   END SUBROUTINE dyn_zdf_imp
353
354   !!==============================================================================
355END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.