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

source: branches/2014/dev_CNRS_CICE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4933

Last change on this file since 4933 was 4933, checked in by cetlod, 10 years ago

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