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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4354

Last change on this file since 4354 was 4354, checked in by jchanut, 10 years ago

Restore AGRIF and BDY compatibility, see ticket #1133
  • Property svn:keywords set to Id
File size: 16.0 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   USE dynspg_ts
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_zdf_imp   ! called by step.F90
35
36   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE dyn_zdf_imp( kt, p2dt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dyn_zdf_imp  ***
51      !!                   
52      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
53      !!      and the surface forcing, and add it to the general trend of
54      !!      the momentum equations.
55      !!
56      !! ** Method  :   The vertical momentum mixing trend is given by :
57      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
58      !!      backward time stepping
59      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
60      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
61      !!      Add this trend to the general trend ua :
62      !!         ua = ua + dz( avmu dz(u) )
63      !!
64      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
65      !!---------------------------------------------------------------------
66      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
67      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
68      !!
69      INTEGER  ::   ji, jj, jk   ! dummy loop indices
70      INTEGER  ::   ikbu, ikbv   ! local integers
71      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
72      REAL(wp) ::   ze3ua, ze3va
73      !!----------------------------------------------------------------------
74
75      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
76      !!----------------------------------------------------------------------
77      !
78      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
79      !
80      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
81      !
82      IF( kt == nit000 ) THEN
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
86         !
87         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
88         ELSE                ;    r_vvl = 0._wp       
89         ENDIF
90      ENDIF
91
92      ! 0. Local constant initialization
93      ! --------------------------------
94      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
95
96      ! 1. Apply semi-implicit bottom friction
97      ! --------------------------------------
98      ! Only needed for semi-implicit bottom friction setup. The explicit
99      ! bottom friction has been included in "u(v)a" which act as the R.H.S
100      ! column vector of the tri-diagonal matrix equation
101      !
102
103      IF( ln_bfrimp ) THEN
104# if defined key_vectopt_loop
105         DO jj = 1, 1
106            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
107# else
108         DO jj = 2, jpjm1
109            DO ji = 2, jpim1
110# endif
111               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
112               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
113               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
114               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
115            END DO
116         END DO
117      ENDIF
118
119#if defined key_dynspg_ts
120      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
121         DO jk = 1, jpkm1
122            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
123            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
124         END DO
125      ELSE                                            ! applied on thickness weighted velocity
126         DO jk = 1, jpkm1
127            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
128               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
129               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
130            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
131               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
132               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
133         END DO
134      ENDIF
135
136      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
137         ! remove barotropic velocities:
138         DO jk = 1, jpkm1
139            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
140            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
141         ENDDO
142         ! Add bottom stress due to barotropic component only:
143         DO jj = 2, jpjm1       
144            DO ji = fs_2, fs_jpim1   ! vector opt.
145               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
146               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
147               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
148               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
149               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
150               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
151            END DO
152         END DO
153      ENDIF
154#endif
155
156      ! 2. Vertical diffusion on u
157      ! ---------------------------
158      ! Matrix and second member construction
159      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
160      ! non zero value at the ocean bottom depending on the bottom friction used.
161      !
162      DO jk = 1, jpkm1        ! Matrix
163         DO jj = 2, jpjm1 
164            DO ji = fs_2, fs_jpim1   ! vector opt.
165               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
166               zcoef = - p2dt / ze3ua     
167               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
168               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
169               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
170               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
171               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
172            END DO
173         END DO
174      END DO
175      DO jj = 2, jpjm1        ! Surface boundary conditions
176         DO ji = fs_2, fs_jpim1   ! vector opt.
177            zwi(ji,jj,1) = 0._wp
178            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
179         END DO
180      END DO
181
182      ! Matrix inversion starting from the first level
183      !-----------------------------------------------------------------------
184      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
185      !
186      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
187      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
188      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
189      !        (        ...               )( ...  ) ( ...  )
190      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
191      !
192      !   m is decomposed in the product of an upper and a lower triangular matrix
193      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
194      !   The solution (the after velocity) is in ua
195      !-----------------------------------------------------------------------
196      !
197      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
198         DO jj = 2, jpjm1   
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
201            END DO
202         END DO
203      END DO
204      !
205      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
206         DO ji = fs_2, fs_jpim1   ! vector opt.
207            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
208#if defined key_dynspg_ts
209            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
210               &                                      / ( ze3ua * rau0 ) 
211#else
212            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
213               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
214#endif
215         END DO
216      END DO
217      DO jk = 2, jpkm1
218         DO jj = 2, jpjm1   
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220#if defined key_dynspg_ts
221               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
222#else
223               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
224#endif
225               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
226            END DO
227         END DO
228      END DO
229      !
230      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
231         DO ji = fs_2, fs_jpim1   ! vector opt.
232            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
233         END DO
234      END DO
235      DO jk = jpk-2, 1, -1
236         DO jj = 2, jpjm1   
237            DO ji = fs_2, fs_jpim1   ! vector opt.
238               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
239            END DO
240         END DO
241      END DO
242
243#if ! defined key_dynspg_ts
244      ! Normalization to obtain the general momentum trend ua
245      DO jk = 1, jpkm1
246         DO jj = 2, jpjm1   
247            DO ji = fs_2, fs_jpim1   ! vector opt.
248               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
249            END DO
250         END DO
251      END DO
252#endif
253
254      ! 3. Vertical diffusion on v
255      ! ---------------------------
256      ! Matrix and second member construction
257      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
258      ! non zero value at the ocean bottom depending on the bottom friction used
259      !
260      DO jk = 1, jpkm1        ! Matrix
261         DO jj = 2, jpjm1   
262            DO ji = fs_2, fs_jpim1   ! vector opt.
263               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
264               zcoef = - p2dt / ze3va
265               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
266               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
267               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
268               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
269               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
270            END DO
271         END DO
272      END DO
273      DO jj = 2, jpjm1        ! Surface boundary conditions
274         DO ji = fs_2, fs_jpim1   ! vector opt.
275            zwi(ji,jj,1) = 0._wp
276            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
277         END DO
278      END DO
279
280      ! Matrix inversion
281      !-----------------------------------------------------------------------
282      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
283      !
284      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
285      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
286      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
287      !        (        ...               )( ...  ) ( ...  )
288      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
289      !
290      !   m is decomposed in the product of an upper and lower triangular matrix
291      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
292      !   The solution (after velocity) is in 2d array va
293      !-----------------------------------------------------------------------
294      !
295      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
296         DO jj = 2, jpjm1   
297            DO ji = fs_2, fs_jpim1   ! vector opt.
298               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
299            END DO
300         END DO
301      END DO
302      !
303      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
304         DO ji = fs_2, fs_jpim1   ! vector opt.
305            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
306#if defined key_dynspg_ts           
307            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
308               &                                      / ( ze3va * rau0 ) 
309#else
310            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
311               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
312#endif
313         END DO
314      END DO
315      DO jk = 2, jpkm1
316         DO jj = 2, jpjm1
317            DO ji = fs_2, fs_jpim1   ! vector opt.
318#if defined key_dynspg_ts
319               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
320#else
321               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
322#endif
323               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
324            END DO
325         END DO
326      END DO
327      !
328      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
329         DO ji = fs_2, fs_jpim1   ! vector opt.
330            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
331         END DO
332      END DO
333      DO jk = jpk-2, 1, -1
334         DO jj = 2, jpjm1   
335            DO ji = fs_2, fs_jpim1   ! vector opt.
336               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
337            END DO
338         END DO
339      END DO
340
341      ! Normalization to obtain the general momentum trend va
342#if ! defined key_dynspg_ts
343      DO jk = 1, jpkm1
344         DO jj = 2, jpjm1   
345            DO ji = fs_2, fs_jpim1   ! vector opt.
346               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
347            END DO
348         END DO
349      END DO
350#endif
351
352      ! J. Chanut: Lines below are useless ?
353      !! restore bottom layer avmu(v)
354      IF( ln_bfrimp ) THEN
355# if defined key_vectopt_loop
356      DO jj = 1, 1
357         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
358# else
359      DO jj = 2, jpjm1
360         DO ji = 2, jpim1
361# endif
362            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
363            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
364            avmu(ji,jj,ikbu+1) = 0.e0
365            avmv(ji,jj,ikbv+1) = 0.e0
366         END DO
367      END DO
368      ENDIF
369      !
370      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
371      !
372      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
373      !
374   END SUBROUTINE dyn_zdf_imp
375
376   !!==============================================================================
377END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.