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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4406

Last change on this file since 4406 was 4406, checked in by trackstand2, 10 years ago

Move from jpk to jpkf - trim sub-domains in z
  • Property svn:keywords set to Id
File size: 15.9 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   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffusion using a implicit time-stepping
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE zdf_oce         ! ocean vertical physics
19   USE phycst          ! physical constants
20   USE in_out_manager  ! I/O manager
21   USE lib_mpp         ! MPP library
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dyn_zdf_imp   ! called by step.F90
27
28   !! * Control permutation of array indices
29#  include "oce_ftrans.h90"
30#  include "dom_oce_ftrans.h90"
31#  include "sbc_oce_ftrans.h90"
32#  include "zdf_oce_ftrans.h90"
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dyn_zdf_imp( kt, p2dt )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE dyn_zdf_imp  ***
47      !!                   
48      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
49      !!      and the surface forcing, and add it to the general trend of
50      !!      the momentum equations.
51      !!
52      !! ** Method  :   The vertical momentum mixing trend is given by :
53      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
54      !!      backward time stepping
55      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
56      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
57      !!      Add this trend to the general trend ua :
58      !!         ua = ua + dz( avmu dz(u) )
59      !!
60      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
61      !!---------------------------------------------------------------------
62      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
63      USE oce     , ONLY:  zwd  => ta       , zws   => sa   ! (ta,sa) used as 3D workspace
64      USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 3D workspace
65      !! DCSE_NEMO: need additional directives for renamed module variables
66!FTRANS zwd :I :I :z
67!FTRANS zws :I :I :z
68!FTRANS zwi :I :I :z
69      !!
70      INTEGER , INTENT(in) ::   kt    ! ocean time-step index
71      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
72      !!
73      INTEGER  ::   ji, jj, jk   ! dummy loop indices
74      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs, zzwibd   ! local scalars
75      !!----------------------------------------------------------------------
76
77      IF( wrk_in_use(3, 3) ) THEN
78         CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN
79      END IF
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      ENDIF
86
87      ! 0. Local constant initialization
88      ! --------------------------------
89      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
90
91
92      ! 1. Vertical diffusion on u
93      ! ---------------------------
94      ! Matrix and second member construction
95      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
96      ! non zero value at the ocean bottom depending on the bottom friction
97      ! used but the bottom velocities have already been updated with the bottom
98      ! friction velocity in dyn_bfr using values from the previous timestep. There
99      ! is no need to include these in the implicit calculation.
100      !
101#if defined key_z_first
102      DO jj = 2, jpjm1 
103         DO ji = 2, jpim1
104            DO jk = 1, jpkfm1
105               zcoef = - p2dt / fse3u(ji,jj,jk)
106               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
107               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
108               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
109               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
110               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
111            END DO
112            ! Surface boundary conditions
113            zwi(ji,jj,1) = 0._wp
114            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
115         END DO
116      END DO
117#else
118      DO jk = 1, jpkfm1        ! Matrix
119         DO jj = 2, jpjm1 
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               zcoef = - p2dt / fse3u(ji,jj,jk)
122               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
123               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
124               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
125               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
126               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
127            END DO
128         END DO
129      END DO
130      DO jj = 2, jpjm1        ! Surface boudary conditions
131         DO ji = fs_2, fs_jpim1   ! vector opt.
132            zwi(ji,jj,1) = 0._wp
133            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
134         END DO
135      END DO
136#endif
137
138      ! Matrix inversion starting from the first level
139      !-----------------------------------------------------------------------
140      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
141      !
142      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
143      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
144      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
145      !        (        ...               )( ...  ) ( ...  )
146      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
147      !
148      !   m is decomposed in the product of an upper and a lower triangular matrix
149      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
150      !   The solution (the after velocity) is in ua
151      !-----------------------------------------------------------------------
152      !
153#if defined key_z_first
154      DO jj = 2, jpjm1 
155         DO ji = 2, jpim1
156            !== Do first and second recurrences in the same loop
157            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
158               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
159            DO jk = 2, jpkfm1
160               zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
161               !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
162               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
163               !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
164               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
165               ua(ji,jj,jk) = zrhs - zzwibd * ua(ji,jj,jk-1)
166            END DO
167            !==  third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
168            ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
169            DO jk = jpkf-2, 1, -1
170               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
171            END DO
172            ! Normalization to obtain the general momentum trend ua
173            DO jk = 1, jpkfm1
174               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
175            END DO
176         END DO
177      END DO
178#else
179      DO jk = 2, jpkfm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
180         DO jj = 2, jpjm1   
181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
183            END DO
184         END DO
185      END DO
186      !
187      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
188         DO ji = fs_2, fs_jpim1   ! vector opt.
189            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
190               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
191         END DO
192      END DO
193      DO jk = 2, jpkfm1
194         DO jj = 2, jpjm1   
195            DO ji = fs_2, fs_jpim1   ! vector opt.
196               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
197               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
198            END DO
199         END DO
200      END DO
201      !
202      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
203         DO ji = fs_2, fs_jpim1   ! vector opt.
204            ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
205         END DO
206      END DO
207      DO jk = jpkf-2, 1, -1
208         DO jj = 2, jpjm1   
209            DO ji = fs_2, fs_jpim1   ! vector opt.
210               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
211            END DO
212         END DO
213      END DO
214      ! Normalization to obtain the general momentum trend ua
215      DO jk = 1, jpkfm1
216         DO jj = 2, jpjm1   
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
219            END DO
220         END DO
221      END DO
222#endif
223
224      ! 2. Vertical diffusion on v
225      ! ---------------------------
226      ! Matrix and second member construction
227      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
228      ! non zero value at the ocean bottom depending on the bottom friction
229      ! used but the bottom velocities have already been updated with the bottom
230      ! friction velocity in dyn_bfr using values from the previous timestep. There
231      ! is no need to include these in the implicit calculation.
232      !
233#if defined key_z_first
234      DO jj = 2, jpjm1   
235         DO ji = 2, jpim1
236            DO jk = 1, jpkfm1        ! Matrix
237               zcoef = -p2dt / fse3v(ji,jj,jk)
238               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
239               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
240               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
241               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
242               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
243            END DO
244            ! Surface boundary conditions
245            zwi(ji,jj,1) = 0._wp
246            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
247         END DO
248      END DO
249#else
250      DO jk = 1, jpkfm1        ! Matrix
251         DO jj = 2, jpjm1   
252            DO ji = fs_2, fs_jpim1   ! vector opt.
253               zcoef = -p2dt / fse3v(ji,jj,jk)
254               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
255               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
256               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
257               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
258               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
259            END DO
260         END DO
261      END DO
262      DO jj = 2, jpjm1        ! Surface boudary conditions
263         DO ji = fs_2, fs_jpim1   ! vector opt.
264            zwi(ji,jj,1) = 0._wp
265            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
266         END DO
267      END DO
268#endif
269
270      ! Matrix inversion
271      !-----------------------------------------------------------------------
272      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
273      !
274      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
275      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
276      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
277      !        (        ...               )( ...  ) ( ...  )
278      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
279      !
280      !   m is decomposed in the product of an upper and lower triangular matrix
281      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
282      !   The solution (after velocity) is in 2d array va
283      !-----------------------------------------------------------------------
284      !
285#if defined key_z_first
286      DO jj = 2, jpjm1   
287         DO ji = 2, jpim1
288            !== Do first and second recurrences in the same loop
289            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
290               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
291            DO jk = 2, jpkfm1
292               zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
293               !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
294               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
295               !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
296               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
297               va(ji,jj,jk) = zrhs - zzwibd * va(ji,jj,jk-1)
298            END DO
299            !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
300            va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
301            DO jk = jpkf-2, 1, -1
302               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
303            END DO
304            ! Normalization to obtain the general momentum trend va
305            DO jk = 1, jpkfm1
306               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
307            END DO
308         END DO
309      END DO
310#else
311      DO jk = 2, jpkfm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
312         DO jj = 2, jpjm1   
313            DO ji = fs_2, fs_jpim1   ! vector opt.
314               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
315            END DO
316         END DO
317      END DO
318      !
319      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
320         DO ji = fs_2, fs_jpim1   ! vector opt.
321            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
322               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
323         END DO
324      END DO
325      DO jk = 2, jpkfm1
326         DO jj = 2, jpjm1
327            DO ji = fs_2, fs_jpim1   ! vector opt.
328               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
329               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
330            END DO
331         END DO
332      END DO
333      !
334      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
335         DO ji = fs_2, fs_jpim1   ! vector opt.
336            va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
337         END DO
338      END DO
339      DO jk = jpkf-2, 1, -1
340         DO jj = 2, jpjm1   
341            DO ji = fs_2, fs_jpim1   ! vector opt.
342               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
343            END DO
344         END DO
345      END DO
346
347      ! Normalization to obtain the general momentum trend va
348      DO jk = 1, jpkfm1
349         DO jj = 2, jpjm1   
350            DO ji = fs_2, fs_jpim1   ! vector opt.
351               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
352            END DO
353         END DO
354      END DO
355#endif
356      !
357      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
358      !
359   END SUBROUTINE dyn_zdf_imp
360
361   !!==============================================================================
362END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.