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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 8093

Last change on this file since 8093 was 8093, checked in by gm, 7 years ago

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