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 @ 4370

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

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