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_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynzdf_imp_tam.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007
  • Property svn:executable set to *
File size: 24.5 KB
Line 
1MODULE dynzdf_imp_tam
2#if defined key_tam
3   !!==============================================================================
4   !!                     ***  MODULE  dynzdf_imp_tam  ***
5   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
6   !!                 Tangent and Adjoint Module
7   !!==============================================================================
8   !! History of the direct module:
9   !!                !  90-10  (B. Blanke)  Original code
10   !!                !  97-05  (G. Madec)  vertical component of isopycnal
11   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
12   !! History of the TAM module:
13   !!           9.0  !  09-01  (A. Vidard) TAM of the 02-08 version
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffu-
18   !!                  sion using an implicit time-stepping scheme.
19   !!----------------------------------------------------------------------
20   !! * Modules used
21   USE par_oce
22   USE oce_tam
23   USE zdf_oce
24   USE dom_oce
25   USE phycst
26   USE in_out_manager
27   USE lib_mpp         ! MPP library
28   USE zdfbfr          ! Bottom friction setup
29   USE wrk_nemo        ! Memory Allocation
30   USE timing          ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   !! * Routine accessibility
36   PUBLIC dyn_zdf_imp_tan     ! called by dynzdf_tam.F90
37   PUBLIC dyn_zdf_imp_adj     ! called by dynzdf_tam.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46   SUBROUTINE dyn_zdf_imp_tan( kt, p2dt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dyn_zdf_imp_tan  ***
49      !!
50      !! ** Purpose of the direct routine:
51      !!      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 of the direct routine:
56      !!      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
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) )E
63      !!
64      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
65      !!               mixing trend.
66      !!---------------------------------------------------------------------
67      !! * Modules used
68      !! * Arguments
69      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
70      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step
71
72      !! * Local declarations
73      INTEGER ::   ji, jj, jk                          ! dummy loop indices
74      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzws, zrhstl ! temporary scalars
75      REAL(wp), POINTER, DIMENSION(:,:,:):: zwi, zws, zwd ! temporary workspace arrays
76      !!----------------------------------------------------------------------
77      !
78      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp_tan')
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_tan : vertical momentum diffusion explicit operator'
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ '
86      ENDIF
87      ! 0. Local constant initialization
88      ! --------------------------------
89      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
90
91      ! 1. Apply semi-implicit bottom friction
92      ! --------------------------------------
93      ! Only needed for semi-implicit bottom friction setup. The explicit
94      ! bottom friction has been included in "u(v)a" which act as the R.H.S
95      ! column vector of the tri-diagonal matrix equation
96      !
97
98      IF( ln_bfrimp ) THEN
99         !!!!!!!!!!!!!!!!!!!!!!!!!!!
100         ! avm* are unactivated for the current TAM
101         !!!!!!!!!!!!!!!!!!!!!!!!!!!
102!# if defined key_vectopt_loop
103      !DO jj = 1, 1
104         !DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
105!# else
106      !DO jj = 2, jpjm1
107         !DO ji = 2, jpim1
108!# endif
109            !ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
110            !ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
111            !zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
112            !zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
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
118      ENDIF
119
120      ! 2. Vertical diffusion on u
121      ! ---------------------------
122      ! Matrix and second member construction
123      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
124      ! non zero value at the ocean bottom depending on the bottom friction used.
125      DO jk = 1, jpkm1
126         DO jj = 2, jpjm1
127            DO ji = fs_2, fs_jpim1   ! vector opt.
128               zcoef = - p2dt / fse3u(ji,jj,jk)
129               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) * umask(ji,jj,jk)
130               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
131               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
132               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
133            END DO
134         END DO
135      END DO
136
137      ! Surface boudary conditions
138      DO jj = 2, jpjm1
139         DO ji = fs_2, fs_jpim1   ! vector opt.
140            zwi(ji,jj,1) = 0.0_wp
141            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
142         END DO
143      END DO
144
145      ! Matrix inversion starting from the first level
146      !-----------------------------------------------------------------------
147      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
148      !
149      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
150      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
151      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
152      !        (        ...               )( ...  ) ( ...  )
153      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
154      !
155      !   m is decomposed in the product of an upper and a lower triangular matrix
156      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
157      !   The solution (the after velocity) is in ua
158      !-----------------------------------------------------------------------
159
160      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
161      DO jk = 2, jpkm1
162         DO jj = 2, jpjm1
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
165            END DO
166         END DO
167      END DO
168
169      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
170      DO jj = 2, jpjm1
171         DO ji = fs_2, fs_jpim1   ! vector opt.
172            ua_tl(ji,jj,1) = ub_tl(ji,jj,1)  &
173                           + p2dt * ua_tl(ji,jj,1)
174         END DO
175      END DO
176      DO jk = 2, jpkm1
177         DO jj = 2, jpjm1
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               zrhstl          = ub_tl(ji,jj,jk) + p2dt * ua_tl(ji,jj,jk)   ! zrhs=right hand side
180               ua_tl(ji,jj,jk) = zrhstl - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua_tl(ji,jj,jk-1)
181            END DO
182         END DO
183      END DO
184
185      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
186      DO jj = 2, jpjm1
187         DO ji = fs_2, fs_jpim1   ! vector opt.
188            ua_tl(ji,jj,jpkm1) = ua_tl(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
189         END DO
190      END DO
191      DO jk = jpk-2, 1, -1
192         DO jj = 2, jpjm1
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               ua_tl(ji,jj,jk) = ( ua_tl(ji,jj,jk) - zws(ji,jj,jk) * ua_tl(ji,jj,jk+1) ) / zwd(ji,jj,jk)
195            END DO
196         END DO
197      END DO
198
199      ! Normalization to obtain the general momentum trend ua
200      DO jk = 1, jpkm1
201         DO jj = 2, jpjm1
202            DO ji = fs_2, fs_jpim1   ! vector opt.
203               ua_tl(ji,jj,jk) = ( ua_tl(ji,jj,jk) - ub_tl(ji,jj,jk) ) * z1_p2dt
204            END DO
205         END DO
206      END DO
207
208
209      ! 2. Vertical diffusion on v
210      ! ---------------------------
211      ! Matrix and second member construction
212      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
213      ! non zero value at the ocean bottom depending on the bottom friction
214      ! used but the bottom velocities have already been updated with the bottom
215      ! friction velocity in dyn_bfr using values from the previous timestep. There
216      ! is no need to include these in the implicit calculation.
217      DO jk = 1, jpkm1
218         DO jj = 2, jpjm1
219            DO ji = fs_2, fs_jpim1   ! vector opt.
220               zcoef          = -p2dt / fse3v(ji,jj,jk)
221               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) * vmask(ji,jj,jk)
222               zzws          = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
223               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
224               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
225            END DO
226         END DO
227      END DO
228
229      ! Surface boudary conditions
230      DO jj = 2, jpjm1
231         DO ji = fs_2, fs_jpim1   ! vector opt.
232            zwi(ji,jj,1) = 0._wp
233            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
234         END DO
235      END DO
236
237      ! Matrix inversion
238      !-----------------------------------------------------------------------
239      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
240      !
241      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
242      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
243      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
244      !        (        ...               )( ...  ) ( ...  )
245      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
246      !
247      !   m is decomposed in the product of an upper and lower triangular
248      !   matrix
249      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
250      !   The solution (after velocity) is in 2d array va
251      !-----------------------------------------------------------------------
252
253      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
254      DO jk = 2, jpkm1
255         DO jj = 2, jpjm1
256            DO ji = fs_2, fs_jpim1   ! vector opt.
257               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
258            END DO
259         END DO
260      END DO
261
262      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
263      DO jj = 2, jpjm1
264         DO ji = fs_2, fs_jpim1   ! vector opt.
265            va_tl(ji,jj,1) = vb_tl(ji,jj,1)  &
266                           + p2dt * va_tl(ji,jj,1)
267         END DO
268      END DO
269      DO jk = 2, jpkm1
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1   ! vector opt.
272               zrhstl          = vb_tl(ji,jj,jk) + p2dt * va_tl(ji,jj,jk)   ! zrhs=right hand side
273               va_tl(ji,jj,jk) = zrhstl - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va_tl(ji,jj,jk-1)
274            END DO
275         END DO
276      END DO
277
278      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
279      DO jj = 2, jpjm1
280         DO ji = fs_2, fs_jpim1   ! vector opt.
281            va_tl(ji,jj,jpkm1) = va_tl(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
282         END DO
283      END DO
284      DO jk = jpk-2, 1, -1
285         DO jj = 2, jpjm1
286            DO ji = fs_2, fs_jpim1   ! vector opt.
287               va_tl(ji,jj,jk) = ( va_tl(ji,jj,jk) - zws(ji,jj,jk) * va_tl(ji,jj,jk+1) ) / zwd(ji,jj,jk)
288            END DO
289         END DO
290      END DO
291
292      ! Normalization to obtain the general momentum trend va
293      DO jk = 1, jpkm1
294         DO jj = 2, jpjm1
295            DO ji = fs_2, fs_jpim1   ! vector opt.
296               va_tl(ji,jj,jk) = ( va_tl(ji,jj,jk) - vb_tl(ji,jj,jk) ) * z1_p2dt
297            END DO
298         END DO
299      END DO
300
301      !! restore bottom layer avmu(v)
302      IF( ln_bfrimp ) THEN
303         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304         ! avm* are unactivated in the current TAM
305         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
306!# if defined key_vectopt_loop
307      !DO jj = 1, 1
308         !DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
309!# else
310      !DO jj = 2, jpjm1
311         !DO ji = 2, jpim1
312!# endif
313            !ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
314            !ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
315            !avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
316            !avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
317         !END DO
318      !END DO
319      ENDIF
320      !
321      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)
322      !
323      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp_tan')
324      !
325   END SUBROUTINE dyn_zdf_imp_tan
326   SUBROUTINE dyn_zdf_imp_adj( kt, p2dt )
327      !!----------------------------------------------------------------------
328      !!                  ***  ROUTINE dyn_zdf_imp_adj  ***
329      !!
330      !! ** Purpose of the direct routine:
331      !!      Compute the trend due to the vert. momentum diffusion
332      !!      and the surface forcing, and add it to the general trend of
333      !!      the momentum equations.
334      !!
335      !! ** Method of the direct routine:
336      !!      The vertical momentum mixing trend is given by :
337      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
338      !!      backward time stepping
339      !!      Surface boundary conditions: wind stress input
340      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
341      !!      Add this trend to the general trend ua :
342      !!         ua = ua + dz( avmu dz(u) )E
343      !!
344      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
345      !!               mixing trend.
346      !!---------------------------------------------------------------------
347      !! * Modules used
348      !! * Arguments
349      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
350      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step
351
352      !! * Local declarations
353      !! * Local declarations
354      INTEGER ::   ji, jj, jk                          ! dummy loop indices
355      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzws, zrhsad ! temporary scalars
356      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zws, zwd! temporary workspace arrays
357      !!----------------------------------------------------------------------
358      !
359      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp_adj')
360      !
361      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )
362      !
363      IF( kt == nitend ) THEN
364         IF(lwp) WRITE(numout,*)
365         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp_adj : vertical momentum diffusion explicit operator'
366         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ '
367      ENDIF
368      ! 0. Local constant initialization
369      ! --------------------------------
370      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
371      zrhsad = 0.0_wp
372      !
373      !! restore bottom layer avmu(v)
374      IF( ln_bfrimp ) THEN
375!# if defined key_vectopt_loop
376      !DO jj = 1, 1
377         !DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
378!# else
379      !DO jj = 2, jpjm1
380         !DO ji = 2, jpim1
381!# endif
382            !ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
383            !ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
384            !avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
385            !avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
386         !END DO
387      !END DO
388      ENDIF
389      !
390      ! 2. Vertical diffusion on v
391      ! ---------------------------
392      ! Matrix and second member construction
393      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
394      ! non zero value at the ocean bottom depending on the bottom friction
395      ! used but the bottom velocities have already been updated with the bottom
396      ! friction velocity in dyn_bfr using values from the previous timestep. There
397      ! is no need to include these in the implicit calculation.
398      DO jk = 1, jpkm1
399         DO jj = 2, jpjm1
400            DO ji = fs_2, fs_jpim1   ! vector opt.
401               zcoef         = -p2dt / fse3v(ji,jj,jk)
402               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) * vmask(ji,jj,jk)
403               zzws          = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
404               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
405               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
406            END DO
407         END DO
408      END DO
409
410      ! Surface boudary conditions
411      DO jj = 2, jpjm1
412         DO ji = fs_2, fs_jpim1   ! vector opt.
413            zwi(ji,jj,1) = 0._wp
414            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
415         END DO
416      END DO
417
418      ! Matrix inversion
419      !-----------------------------------------------------------------------
420      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
421      !
422      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
423      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
424      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
425      !        (        ...               )( ...  ) ( ...  )
426      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
427      !
428      !   m is decomposed in the product of an upper and lower triangular
429      !   matrix
430      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
431      !   The solution (after velocity) is in 2d array va
432      !-----------------------------------------------------------------------
433
434      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
435      DO jk = 2, jpkm1
436         DO jj = 2, jpjm1
437            DO ji = fs_2, fs_jpim1   ! vector opt.
438               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
439            END DO
440         END DO
441      END DO
442
443      ! Normalization to obtain the general momentum trend va
444      DO jk = jpkm1, 1, -1
445         DO jj = jpjm1, 2, -1
446            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
447               vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) - va_ad(ji,jj,jk) / p2dt
448               va_ad(ji,jj,jk) = va_ad(ji,jj,jk) / p2dt
449            END DO
450         END DO
451      END DO
452      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
453      DO jk = 1, jpk-2
454         DO jj = jpjm1, 2, -1
455            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
456               va_ad(ji,jj,jk+1) = va_ad(ji,jj,jk+1) - zws(ji,jj,jk) * va_ad(ji,jj,jk) / zwd(ji,jj,jk)
457               va_ad(ji,jj,jk  ) = va_ad(ji,jj,jk  ) / zwd(ji,jj,jk)
458            END DO
459         END DO
460      END DO
461      DO jj = jpjm1, 2, -1
462         DO ji = fs_jpim1, fs_2, -1   ! vector opt.
463            va_ad(ji,jj,jpkm1) = va_ad(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
464         END DO
465      END DO
466      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
467      DO jk = jpkm1, 2, -1
468         DO jj = jpjm1, 2, -1
469            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
470               zrhsad            = zrhsad + va_ad(ji,jj,jk)
471               va_ad(ji,jj,jk-1) = va_ad(ji,jj,jk-1) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va_ad(ji,jj,jk)
472               va_ad(ji,jj,jk  ) = 0.0_wp
473               vb_ad(ji,jj,jk)   = vb_ad(ji,jj,jk) + zrhsad
474               va_ad(ji,jj,jk)   = va_ad(ji,jj,jk) + p2dt * zrhsad
475               zrhsad            = 0.0_wp
476            END DO
477         END DO
478      END DO
479      DO jj = jpjm1, 2, -1
480         DO ji = fs_jpim1, fs_2, -1   ! vector opt.
481            vb_ad(ji,jj,1) = vb_ad(ji,jj,1) + va_ad(ji,jj,1)
482            va_ad(ji,jj,1) = va_ad(ji,jj,1) * p2dt
483         END DO
484      END DO
485
486      ! 1. Vertical diffusion on u
487      ! ---------------------------
488      ! Matrix and second member construction
489      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
490      ! non zero value at the ocean bottom depending on the bottom friction
491      ! used but the bottom velocities have already been updated with the bottom
492      ! friction velocity in dyn_bfr using values from the previous timestep. There
493      ! is no need to include these in the implicit calculation.
494      DO jk = 1, jpkm1
495         DO jj = 2, jpjm1
496            DO ji = fs_2, fs_jpim1   ! vector opt.
497               zcoef = - p2dt / fse3u(ji,jj,jk)
498               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) * umask(ji,jj,jk)
499               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
500               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
501               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
502            END DO
503         END DO
504      END DO
505
506      ! Surface boudary conditions
507      DO jj = 2, jpjm1
508         DO ji = fs_2, fs_jpim1   ! vector opt.
509            zwi(ji,jj,1) = 0._wp
510            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
511         END DO
512      END DO
513
514      ! Matrix inversion starting from the first level
515      !-----------------------------------------------------------------------
516      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
517      !
518      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
519      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
520      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
521      !        (        ...               )( ...  ) ( ...  )
522      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
523      !
524      !   m is decomposed in the product of an upper and a lower triangular matrix
525      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
526      !   The solution (the after velocity) is in ua
527      !-----------------------------------------------------------------------
528
529     ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
530      DO jk = 2, jpkm1
531         DO jj = 2, jpjm1
532            DO ji = fs_2, fs_jpim1   ! vector opt.
533               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
534            END DO
535         END DO
536      END DO
537      ! Normalization to obtain the general momentum trend ua
538      DO jk = jpkm1, 1, -1
539         DO jj = jpjm1, 2, -1
540            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
541               ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) - ua_ad(ji,jj,jk) / p2dt
542               ua_ad(ji,jj,jk) = ua_ad(ji,jj,jk) / p2dt
543            END DO
544         END DO
545      END DO
546      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
547      DO jk = 1, jpk-2
548         DO jj = jpjm1, 2, -1
549            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
550               ua_ad(ji,jj,jk+1) = ua_ad(ji,jj,jk+1) - ua_ad(ji,jj,jk) * zws(ji,jj,jk) / zwd(ji,jj,jk)
551               ua_ad(ji,jj,jk)   = ua_ad(ji,jj,jk) / zwd(ji,jj,jk)
552            END DO
553         END DO
554      END DO
555      DO jj = jpjm1, 2, -1
556         DO ji = fs_jpim1, fs_2, -1
557            ua_ad(ji,jj,jpkm1) = ua_ad(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
558         END DO
559      END DO
560      DO jk = jpkm1, 2, -1
561         DO jj = jpjm1, 2, -1
562            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
563               zrhsad            = zrhsad + ua_ad(ji,jj,jk)
564               ua_ad(ji,jj,jk-1) = ua_ad(ji,jj,jk-1) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua_ad(ji,jj,jk)
565               ua_ad(ji,jj,jk)   = 0.0_wp
566               ub_ad(ji,jj,jk)   = ub_ad(ji,jj,jk) + zrhsad
567               ua_ad(ji,jj,jk)   = ua_ad(ji,jj,jk) + zrhsad * p2dt
568               zrhsad            = 0.0_wp
569            END DO
570         END DO
571      END DO
572      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
573      DO jj = 2, jpjm1
574         DO ji = fs_2, fs_jpim1   ! vector opt.
575            ub_ad(ji,jj,1) = ub_ad(ji,jj,1) + ua_ad(ji,jj,1)
576            ua_ad(ji,jj,1) = p2dt * ua_ad(ji,jj,1)
577         END DO
578      END DO
579      !
580      ! 1. Apply semi-implicit bottom friction
581      ! --------------------------------------
582      ! Only needed for semi-implicit bottom friction setup. The explicit
583      ! bottom friction has been included in "u(v)a" which act as the R.H.S
584      ! column vector of the tri-diagonal matrix equation
585      !
586      IF( ln_bfrimp ) THEN
587!# if defined key_vectopt_loop
588      !DO jj = 1, 1
589         !DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
590!# else
591      !DO jj = 2, jpjm1
592         !DO ji = 2, jpim1
593!# endif
594            !ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
595            !ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
596            !zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
597            !zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
598            !avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
599            !avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
600         !END DO
601      !END DO
602      ENDIF
603      !
604      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)
605      !
606      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp_adj')
607      !
608   END SUBROUTINE dyn_zdf_imp_adj
609#endif
610   !!==============================================================================
611END MODULE dynzdf_imp_tam
Note: See TracBrowser for help on using the repository browser.