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/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation
  • Property svn:keywords set to Id
File size: 15.8 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.7 , NEMO Consortium (2014)
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      !!----------------------------------------------------------------------
73
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( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
87         ELSE                ;    r_vvl = 0._wp       
88         ENDIF
89      ENDIF
90
91      ! 0. Local constant initialization
92      ! --------------------------------
93      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
94
95      ! 1. Apply semi-implicit bottom friction
96      ! --------------------------------------
97      ! Only needed for semi-implicit bottom friction setup. The explicit
98      ! bottom friction has been included in "u(v)a" which act as the R.H.S
99      ! column vector of the tri-diagonal matrix equation
100      !
101
102      IF( ln_bfrimp ) THEN
103         DO jj = 2, jpjm1
104            DO ji = 2, jpim1
105               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
106               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
107               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
108               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
109            END DO
110         END DO
111      ENDIF
112
113#if defined key_dynspg_ts
114      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
115         DO jk = 1, jpkm1
116            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
117            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
118         END DO
119      ELSE                                            ! applied on thickness weighted velocity
120         DO jk = 1, jpkm1
121            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
122               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
123               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
124            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
125               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
126               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
127         END DO
128      ENDIF
129
130      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
131         ! remove barotropic velocities:
132         DO jk = 1, jpkm1
133            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
134            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
135         ENDDO
136         ! Add bottom stress due to barotropic component only:
137         DO jj = 2, jpjm1       
138            DO ji = fs_2, fs_jpim1   ! vector opt.
139               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
140               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
141               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
142               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
143               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
144               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
145            END DO
146         END DO
147      ENDIF
148#endif
149
150      ! 2. Vertical diffusion on u
151      ! ---------------------------
152      ! Matrix and second member construction
153      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
154      ! non zero value at the ocean bottom depending on the bottom friction used.
155      !
156      DO jk = 1, jpkm1        ! Matrix
157         DO jj = 2, jpjm1 
158            DO ji = fs_2, fs_jpim1   ! vector opt.
159               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
160               zcoef = - p2dt / ze3ua     
161               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
162               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
163               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
164               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
165               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
166            END DO
167         END DO
168      END DO
169      DO jj = 2, jpjm1        ! Surface boundary conditions
170         DO ji = fs_2, fs_jpim1   ! vector opt.
171            zwi(ji,jj,1) = 0._wp
172            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
173         END DO
174      END DO
175
176      ! Matrix inversion starting from the first level
177      !-----------------------------------------------------------------------
178      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
179      !
180      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
181      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
182      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
183      !        (        ...               )( ...  ) ( ...  )
184      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
185      !
186      !   m is decomposed in the product of an upper and a lower triangular matrix
187      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
188      !   The solution (the after velocity) is in ua
189      !-----------------------------------------------------------------------
190      !
191      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
192         DO jj = 2, jpjm1   
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
195            END DO
196         END DO
197      END DO
198      !
199      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
200         DO ji = fs_2, fs_jpim1   ! vector opt.
201            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
202#if defined key_dynspg_ts
203            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
204               &                                      / ( ze3ua * rau0 ) 
205#else
206            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
207               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
208#endif
209         END DO
210      END DO
211      DO jk = 2, jpkm1
212         DO jj = 2, jpjm1   
213            DO ji = fs_2, fs_jpim1   ! vector opt.
214#if defined key_dynspg_ts
215               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
216#else
217               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
218#endif
219               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
220            END DO
221         END DO
222      END DO
223      !
224      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
225         DO ji = fs_2, fs_jpim1   ! vector opt.
226            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
227         END DO
228      END DO
229      DO jk = jpk-2, 1, -1
230         DO jj = 2, jpjm1   
231            DO ji = fs_2, fs_jpim1   ! vector opt.
232               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
233            END DO
234         END DO
235      END DO
236
237#if ! defined key_dynspg_ts
238      ! Normalization to obtain the general momentum trend ua
239      DO jk = 1, jpkm1
240         DO jj = 2, jpjm1   
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
243            END DO
244         END DO
245      END DO
246#endif
247
248      ! 3. Vertical diffusion on v
249      ! ---------------------------
250      ! Matrix and second member construction
251      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
252      ! non zero value at the ocean bottom depending on the bottom friction used
253      !
254      DO jk = 1, jpkm1        ! Matrix
255         DO jj = 2, jpjm1   
256            DO ji = fs_2, fs_jpim1   ! vector opt.
257               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
258               zcoef = - p2dt / ze3va
259               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
260               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
261               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
262               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
263               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
264            END DO
265         END DO
266      END DO
267      DO jj = 2, jpjm1        ! Surface boundary conditions
268         DO ji = fs_2, fs_jpim1   ! vector opt.
269            zwi(ji,jj,1) = 0._wp
270            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
271         END DO
272      END DO
273
274      ! Matrix inversion
275      !-----------------------------------------------------------------------
276      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
277      !
278      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
279      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
280      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
281      !        (        ...               )( ...  ) ( ...  )
282      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
283      !
284      !   m is decomposed in the product of an upper and lower triangular matrix
285      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
286      !   The solution (after velocity) is in 2d array va
287      !-----------------------------------------------------------------------
288      !
289      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
290         DO jj = 2, jpjm1   
291            DO ji = fs_2, fs_jpim1   ! vector opt.
292               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
293            END DO
294         END DO
295      END DO
296      !
297      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
298         DO ji = fs_2, fs_jpim1   ! vector opt.
299            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
300#if defined key_dynspg_ts           
301            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
302               &                                      / ( ze3va * rau0 ) 
303#else
304            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
305               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
306#endif
307         END DO
308      END DO
309      DO jk = 2, jpkm1
310         DO jj = 2, jpjm1
311            DO ji = fs_2, fs_jpim1   ! vector opt.
312#if defined key_dynspg_ts
313               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
314#else
315               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
316#endif
317               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
318            END DO
319         END DO
320      END DO
321      !
322      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
323         DO ji = fs_2, fs_jpim1   ! vector opt.
324            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
325         END DO
326      END DO
327      DO jk = jpk-2, 1, -1
328         DO jj = 2, jpjm1   
329            DO ji = fs_2, fs_jpim1   ! vector opt.
330               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
331            END DO
332         END DO
333      END DO
334
335      ! Normalization to obtain the general momentum trend va
336#if ! defined key_dynspg_ts
337      DO jk = 1, jpkm1
338         DO jj = 2, jpjm1   
339            DO ji = fs_2, fs_jpim1   ! vector opt.
340               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
341            END DO
342         END DO
343      END DO
344#endif
345
346      ! J. Chanut: Lines below are useless ?
347      !! restore bottom layer avmu(v)
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      ENDIF
358      !
359      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
362      !
363   END SUBROUTINE dyn_zdf_imp
364
365   !!==============================================================================
366END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.