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

source: branches/2016/dev_v3_6_STABLE_OMP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 6712

Last change on this file since 6712 was 6712, checked in by mocavero, 8 years ago

update some OMP directives
  • Property svn:keywords set to Id
File size: 18.4 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!$OMP PARALLEL DO schedule(static) private(jj, ji, ikbu, ikbv)
102         DO jj = 2, jpjm1
103            DO ji = 2, jpim1
104               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
105               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
106               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
107               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
108            END DO
109         END DO
110         IF ( ln_isfcav ) THEN
111!$OMP PARALLEL DO schedule(static) private(jj, ji, ikbu, ikbv)
112            DO jj = 2, jpjm1
113               DO ji = 2, jpim1
114                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
115                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
116                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
117                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
118               END DO
119            END DO
120         END IF
121      ENDIF
122
123#if defined key_dynspg_ts
124      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
125         DO jk = 1, jpkm1
126            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
127            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
128         END DO
129      ELSE                                            ! applied on thickness weighted velocity
130         DO jk = 1, jpkm1
131            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
132               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
133               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
134            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
135               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
136               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
137         END DO
138      ENDIF
139
140      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
141         ! remove barotropic velocities:
142         DO jk = 1, jpkm1
143            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
144            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
145         END DO
146         ! Add bottom/top stress due to barotropic component only:
147         DO jj = 2, jpjm1       
148            DO ji = fs_2, fs_jpim1   ! vector opt.
149               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
150               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
151               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
152               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
153               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
154               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
155            END DO
156         END DO
157         IF ( ln_isfcav ) THEN
158            DO jj = 2, jpjm1       
159               DO ji = fs_2, fs_jpim1   ! vector opt.
160                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
161                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
162                  ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
163                  ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
164                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
165                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
166               END DO
167            END DO
168         END IF
169      ENDIF
170#endif
171
172      ! 2. Vertical diffusion on u
173      ! ---------------------------
174      ! Matrix and second member construction
175      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
176      ! non zero value at the ocean bottom depending on the bottom friction used.
177      !
178!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ze3ua, zcoef, zzwi, zzws)
179      DO jk = 1, jpkm1        ! Matrix
180         DO jj = 2, jpjm1 
181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
183               zcoef = - p2dt / ze3ua     
184               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
185               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
186               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
187               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
188               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
189            END DO
190         END DO
191      END DO
192      DO jj = 2, jpjm1        ! Surface boundary conditions
193         DO ji = fs_2, fs_jpim1   ! vector opt.
194            zwi(ji,jj,1) = 0._wp
195            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
196         END DO
197      END DO
198
199      ! Matrix inversion starting from the first level
200      !-----------------------------------------------------------------------
201      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
202      !
203      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
204      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
205      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
206      !        (        ...               )( ...  ) ( ...  )
207      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
208      !
209      !   m is decomposed in the product of an upper and a lower triangular matrix
210      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
211      !   The solution (the after velocity) is in ua
212      !-----------------------------------------------------------------------
213      !
214      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
215      DO jk = 2, jpkm1
216!$OMP PARALLEL DO schedule(static) private(jj, ji)
217         DO jj = 2, jpjm1   
218            DO ji = fs_2, fs_jpim1   ! vector opt.
219               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
220            END DO
221         END DO
222      END DO
223      !
224!$OMP PARALLEL DO schedule(static) private(jj, ji, ze3ua)     
225      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
226         DO ji = fs_2, fs_jpim1   ! vector opt.
227#if defined key_dynspg_ts
228            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
229            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
230               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
231#else
232            ua(ji,jj,1) = ub(ji,jj,1) &
233               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
234               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) ) 
235#endif
236         END DO
237      END DO
238      DO jk = 2, jpkm1
239!$OMP PARALLEL DO schedule(static) private(jj, ji, zrhs)
240         DO jj = 2, jpjm1
241            DO ji = fs_2, fs_jpim1
242#if defined key_dynspg_ts
243               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
244#else
245               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
246#endif
247               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
248            END DO
249         END DO
250      END DO
251      !
252!$OMP PARALLEL DO schedule(static) private(jj, ji)
253      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
254         DO ji = fs_2, fs_jpim1   ! vector opt.
255            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
256         END DO
257      END DO
258      DO jk = jpk-2, 1, -1
259!$OMP PARALLEL DO schedule(static) private(jj, ji)
260         DO jj = 2, jpjm1
261            DO ji = fs_2, fs_jpim1
262               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
263            END DO
264         END DO
265      END DO
266
267#if ! defined key_dynspg_ts
268      ! Normalization to obtain the general momentum trend ua
269!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
270      DO jk = 1, jpkm1
271         DO jj = 2, jpjm1   
272            DO ji = fs_2, fs_jpim1   ! vector opt.
273               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
274            END DO
275         END DO
276      END DO
277#endif
278
279      ! 3. Vertical diffusion on v
280      ! ---------------------------
281      ! Matrix and second member construction
282      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
283      ! non zero value at the ocean bottom depending on the bottom friction used
284      !
285!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ze3va, zcoef, zzwi, zzws)
286      DO jk = 1, jpkm1        ! Matrix
287         DO jj = 2, jpjm1   
288            DO ji = fs_2, fs_jpim1   ! vector opt.
289               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
290               zcoef = - p2dt / ze3va
291               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
292               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
293               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
294               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
295               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
296            END DO
297         END DO
298      END DO
299!$OMP PARALLEL DO schedule(static) private(jj, ji)
300      DO jj = 2, jpjm1        ! Surface boundary conditions
301         DO ji = fs_2, fs_jpim1   ! vector opt.
302            zwi(ji,jj,1) = 0._wp
303            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
304         END DO
305      END DO
306
307      ! Matrix inversion
308      !-----------------------------------------------------------------------
309      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
310      !
311      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
312      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
313      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
314      !        (        ...               )( ...  ) ( ...  )
315      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
316      !
317      !   m is decomposed in the product of an upper and lower triangular matrix
318      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
319      !   The solution (after velocity) is in 2d array va
320      !-----------------------------------------------------------------------
321      !
322      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
323      DO jk = 2, jpkm1       
324!$OMP PARALLEL DO schedule(static) private(jj, ji)
325         DO jj = 2, jpjm1   
326            DO ji = fs_2, fs_jpim1   ! vector opt.
327               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
328            END DO
329         END DO
330      END DO
331      !
332!$OMP PARALLEL DO schedule(static) private(jj, ji, ze3va)
333      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
334         DO ji = fs_2, fs_jpim1   ! vector opt.
335#if defined key_dynspg_ts           
336            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
337            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
338               &                                      / ( ze3va * rau0 ) 
339#else
340            va(ji,jj,1) = vb(ji,jj,1) &
341               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
342               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
343#endif
344         END DO
345      END DO
346      DO jk = 2, jpkm1
347!$OMP PARALLEL DO schedule(static) private(jj, ji, zrhs)
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350#if defined key_dynspg_ts
351               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
352#else
353               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
354#endif
355               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
356            END DO
357         END DO
358      END DO
359      !
360!$OMP PARALLEL DO schedule(static) private(jj, ji)
361      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
362         DO ji = fs_2, fs_jpim1   ! vector opt.
363            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
364         END DO
365      END DO
366      DO jk = jpk-2, 1, -1
367!$OMP PARALLEL DO schedule(static) private(jj, ji)
368         DO jj = 2, jpjm1
369            DO ji = fs_2, fs_jpim1
370               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
371            END DO
372         END DO
373      END DO
374
375      ! Normalization to obtain the general momentum trend va
376#if ! defined key_dynspg_ts
377!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
378      DO jk = 1, jpkm1
379         DO jj = 2, jpjm1   
380            DO ji = fs_2, fs_jpim1   ! vector opt.
381               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
382            END DO
383         END DO
384      END DO
385#endif
386
387      ! J. Chanut: Lines below are useless ?
388      !! restore bottom layer avmu(v)
389      IF( ln_bfrimp ) THEN
390!$OMP PARALLEL DO schedule(static) private(jj, ji, ikbu, ikbv)
391        DO jj = 2, jpjm1
392           DO ji = 2, jpim1
393              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
394              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
395              avmu(ji,jj,ikbu+1) = 0.e0
396              avmv(ji,jj,ikbv+1) = 0.e0
397           END DO
398        END DO
399        IF (ln_isfcav) THEN
400!$OMP PARALLEL DO schedule(static) private(jj, ji, ikbu, ikbv)
401           DO jj = 2, jpjm1
402              DO ji = 2, jpim1
403                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
404                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
405                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
406                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
407              END DO
408           END DO
409        END IF
410      ENDIF
411      !
412      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
413      !
414      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
415      !
416   END SUBROUTINE dyn_zdf_imp
417
418   !!==============================================================================
419END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.