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 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

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