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

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago

add wmask, wumask, wvmask and restore loop order and add flag to ignore specific isf code if no isf
  • Property svn:keywords set to Id
File size: 17.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          ! 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               ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
154               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
155               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
156               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
157               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
158               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
159            END DO
160         END DO
161      ENDIF
162#endif
163
164      ! 2. Vertical diffusion on u
165      ! ---------------------------
166      ! Matrix and second member construction
167      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
168      ! non zero value at the ocean bottom depending on the bottom friction used.
169      !
170      DO jk = 1, jpkm1        ! Matrix
171         DO jj = 2, jpjm1 
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
174               zcoef = - p2dt / ze3ua     
175               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
176               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk)
177               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
178               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
179               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
180            END DO
181         END DO
182      END DO
183      DO jj = 2, jpjm1        ! Surface boundary conditions
184         DO ji = fs_2, fs_jpim1   ! vector opt.
185            zwi(ji,jj,1) = 0._wp
186            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
187         END DO
188      END DO
189
190      ! Matrix inversion starting from the first level
191      !-----------------------------------------------------------------------
192      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
193      !
194      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
195      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
196      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
197      !        (        ...               )( ...  ) ( ...  )
198      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
199      !
200      !   m is decomposed in the product of an upper and a lower triangular matrix
201      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
202      !   The solution (the after velocity) is in ua
203      !-----------------------------------------------------------------------
204      !
205      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
206      DO jk = 2, jpkm1
207         DO jj = 2, jpjm1   
208            DO ji = fs_2, fs_jpim1   ! vector opt.
209               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
210            END DO
211         END DO
212      END DO
213      !
214      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
215         DO ji = fs_2, fs_jpim1   ! vector opt.
216#if defined key_dynspg_ts
217            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
218            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
219               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
220#else
221            ua(ji,jj,1) = ub(ji,jj,1) &
222               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
223               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) ) 
224#endif
225         END DO
226      END DO
227      DO jk = 2, jpkm1
228         DO jj = 2, jpjm1
229            DO ji = fs_2, fs_jpim1
230#if defined key_dynspg_ts
231               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
232#else
233               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
234#endif
235               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
236            END DO
237         END DO
238      END DO
239      !
240      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
241         DO ji = fs_2, fs_jpim1   ! vector opt.
242            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
243         END DO
244      END DO
245      DO jk = jpk-2, 1, -1
246         DO jj = 2, jpjm1
247            DO ji = fs_2, fs_jpim1
248               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
249            END DO
250         END DO
251      END DO
252
253#if ! defined key_dynspg_ts
254      ! Normalization to obtain the general momentum trend ua
255      DO jk = 1, jpkm1
256         DO jj = 2, jpjm1   
257            DO ji = fs_2, fs_jpim1   ! vector opt.
258               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
259            END DO
260         END DO
261      END DO
262#endif
263
264      ! 3. Vertical diffusion on v
265      ! ---------------------------
266      ! Matrix and second member construction
267      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
268      ! non zero value at the ocean bottom depending on the bottom friction used
269      !
270      DO jk = 1, jpkm1        ! Matrix
271         DO jj = 2, jpjm1   
272            DO ji = fs_2, fs_jpim1   ! vector opt.
273               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
274               zcoef = - p2dt / ze3va
275               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
276               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
277               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
278               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
279               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
280            END DO
281         END DO
282      END DO
283      DO jj = 2, jpjm1        ! Surface boundary conditions
284         DO ji = fs_2, fs_jpim1   ! vector opt.
285            zwi(ji,jj,1) = 0._wp
286            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
287         END DO
288      END DO
289
290      ! Matrix inversion
291      !-----------------------------------------------------------------------
292      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
293      !
294      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
295      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
296      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
297      !        (        ...               )( ...  ) ( ...  )
298      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
299      !
300      !   m is decomposed in the product of an upper and lower triangular matrix
301      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
302      !   The solution (after velocity) is in 2d array va
303      !-----------------------------------------------------------------------
304      !
305      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
306      DO jk = 2, jpkm1       
307         DO jj = 2, jpjm1   
308            DO ji = fs_2, fs_jpim1   ! vector opt.
309               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
310            END DO
311         END DO
312      END DO
313      !
314      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
315         DO ji = fs_2, fs_jpim1   ! vector opt.
316#if defined key_dynspg_ts           
317            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
318            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
319               &                                      / ( ze3va * rau0 ) 
320#else
321            va(ji,jj,1) = vb(ji,jj,1) &
322               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
323               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
324#endif
325         END DO
326      END DO
327      DO jk = 2, jpkm1
328         DO jj = 2, jpjm1
329            DO ji = fs_2, fs_jpim1   ! vector opt.
330#if defined key_dynspg_ts
331               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
332#else
333               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
334#endif
335               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
336            END DO
337         END DO
338      END DO
339      !
340      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
341         DO ji = fs_2, fs_jpim1   ! vector opt.
342            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
343         END DO
344      END DO
345      DO jk = jpk-2, 1, -1
346         DO jj = 2, jpjm1
347            DO ji = fs_2, fs_jpim1
348               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
349            END DO
350         END DO
351      END DO
352
353      ! Normalization to obtain the general momentum trend va
354#if ! defined key_dynspg_ts
355      DO jk = 1, jpkm1
356         DO jj = 2, jpjm1   
357            DO ji = fs_2, fs_jpim1   ! vector opt.
358               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
359            END DO
360         END DO
361      END DO
362#endif
363
364      ! J. Chanut: Lines below are useless ?
365      !! restore bottom layer avmu(v)
366      IF( ln_bfrimp ) THEN
367        DO jj = 2, jpjm1
368           DO ji = 2, jpim1
369              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
370              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
371              avmu(ji,jj,ikbu+1) = 0.e0
372              avmv(ji,jj,ikbv+1) = 0.e0
373           END DO
374        END DO
375        IF (ln_isfcav) THEN
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.e0
381                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
382              END DO
383           END DO
384        END IF
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.