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

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

Added timing calls
  • 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   !!----------------------------------------------------------------------
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      USE timing,   ONLY: timing_start, timing_stop
66      !! DCSE_NEMO: need additional directives for renamed module variables
67!FTRANS zwd :I :I :z
68!FTRANS zws :I :I :z
69!FTRANS zwi :I :I :z
70      !!
71      INTEGER , INTENT(in) ::   kt    ! ocean time-step index
72      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
73      !!
74      INTEGER  ::   ji, jj, jk   ! dummy loop indices
75      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs, zzwibd   ! local scalars
76      !!----------------------------------------------------------------------
77
78      CALL timing_start('dyn_zdf_imp')
79
80      IF( wrk_in_use(3, 3) ) THEN
81         CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN
82      END IF
83
84      IF( kt == nit000 ) THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
87         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
88      ENDIF
89
90      ! 0. Local constant initialization
91      ! --------------------------------
92      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
93
94
95      ! 1. Vertical diffusion on u
96      ! ---------------------------
97      ! Matrix and second member construction
98      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
99      ! non zero value at the ocean bottom depending on the bottom friction
100      ! used but the bottom velocities have already been updated with the bottom
101      ! friction velocity in dyn_bfr using values from the previous timestep. There
102      ! is no need to include these in the implicit calculation.
103      !
104#if defined key_z_first
105      DO jj = 2, jpjm1 
106         DO ji = 2, jpim1
107            DO jk = 1, jpkfm1
108               zcoef = - p2dt / fse3u(ji,jj,jk)
109               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
110               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
111               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
112               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
113               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
114            END DO
115            ! Surface boundary conditions
116            zwi(ji,jj,1) = 0._wp
117            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
118         END DO
119      END DO
120#else
121      DO jk = 1, jpkfm1        ! Matrix
122         DO jj = 2, jpjm1 
123            DO ji = fs_2, fs_jpim1   ! vector opt.
124               zcoef = - p2dt / fse3u(ji,jj,jk)
125               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
126               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
127               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
128               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
129               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
130            END DO
131         END DO
132      END DO
133      DO jj = 2, jpjm1        ! Surface boudary conditions
134         DO ji = fs_2, fs_jpim1   ! vector opt.
135            zwi(ji,jj,1) = 0._wp
136            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
137         END DO
138      END DO
139#endif
140
141      ! Matrix inversion starting from the first level
142      !-----------------------------------------------------------------------
143      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
144      !
145      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
146      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
147      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
148      !        (        ...               )( ...  ) ( ...  )
149      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
150      !
151      !   m is decomposed in the product of an upper and a lower triangular matrix
152      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
153      !   The solution (the after velocity) is in ua
154      !-----------------------------------------------------------------------
155      !
156#if defined key_z_first
157      DO jj = 2, jpjm1 
158         DO ji = 2, jpim1
159            !== Do first and second recurrences in the same loop
160            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
161               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
162            DO jk = 2, jpkfm1
163               zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
164               !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
165               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
166               !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
167               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
168               ua(ji,jj,jk) = zrhs - zzwibd * ua(ji,jj,jk-1)
169            END DO
170            !==  third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
171            ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
172            DO jk = jpkf-2, 1, -1
173               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
174            END DO
175            ! Normalization to obtain the general momentum trend ua
176            DO jk = 1, jpkfm1
177               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
178            END DO
179         END DO
180      END DO
181#else
182      DO jk = 2, jpkfm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
183         DO jj = 2, jpjm1   
184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
186            END DO
187         END DO
188      END DO
189      !
190      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
191         DO ji = fs_2, fs_jpim1   ! vector opt.
192            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
193               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
194         END DO
195      END DO
196      DO jk = 2, jpkfm1
197         DO jj = 2, jpjm1   
198            DO ji = fs_2, fs_jpim1   ! vector opt.
199               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
200               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
201            END DO
202         END DO
203      END DO
204      !
205      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
206         DO ji = fs_2, fs_jpim1   ! vector opt.
207            ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
208         END DO
209      END DO
210      DO jk = jpkf-2, 1, -1
211         DO jj = 2, jpjm1   
212            DO ji = fs_2, fs_jpim1   ! vector opt.
213               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
214            END DO
215         END DO
216      END DO
217      ! Normalization to obtain the general momentum trend ua
218      DO jk = 1, jpkfm1
219         DO jj = 2, jpjm1   
220            DO ji = fs_2, fs_jpim1   ! vector opt.
221               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
222            END DO
223         END DO
224      END DO
225#endif
226
227      ! 2. Vertical diffusion on v
228      ! ---------------------------
229      ! Matrix and second member construction
230      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
231      ! non zero value at the ocean bottom depending on the bottom friction
232      ! used but the bottom velocities have already been updated with the bottom
233      ! friction velocity in dyn_bfr using values from the previous timestep. There
234      ! is no need to include these in the implicit calculation.
235      !
236#if defined key_z_first
237      DO jj = 2, jpjm1   
238         DO ji = 2, jpim1
239            DO jk = 1, jpkfm1        ! Matrix
240               zcoef = -p2dt / fse3v(ji,jj,jk)
241               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
242               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
243               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
244               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
245               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
246            END DO
247            ! Surface boundary conditions
248            zwi(ji,jj,1) = 0._wp
249            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
250         END DO
251      END DO
252#else
253      DO jk = 1, jpkfm1        ! Matrix
254         DO jj = 2, jpjm1   
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               zcoef = -p2dt / fse3v(ji,jj,jk)
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 boudary 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#endif
272
273      ! Matrix inversion
274      !-----------------------------------------------------------------------
275      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
276      !
277      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
278      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
279      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
280      !        (        ...               )( ...  ) ( ...  )
281      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
282      !
283      !   m is decomposed in the product of an upper and lower triangular matrix
284      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
285      !   The solution (after velocity) is in 2d array va
286      !-----------------------------------------------------------------------
287      !
288#if defined key_z_first
289      DO jj = 2, jpjm1   
290         DO ji = 2, jpim1
291            !== Do first and second recurrences in the same loop
292            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
293               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
294            DO jk = 2, jpkfm1
295               zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
296               !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
297               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
298               !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
299               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
300               va(ji,jj,jk) = zrhs - zzwibd * va(ji,jj,jk-1)
301            END DO
302            !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
303            va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
304            DO jk = jpkf-2, 1, -1
305               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
306            END DO
307            ! Normalization to obtain the general momentum trend va
308            DO jk = 1, jpkfm1
309               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
310            END DO
311         END DO
312      END DO
313#else
314      DO jk = 2, jpkfm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
315         DO jj = 2, jpjm1   
316            DO ji = fs_2, fs_jpim1   ! vector opt.
317               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
318            END DO
319         END DO
320      END DO
321      !
322      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
323         DO ji = fs_2, fs_jpim1   ! vector opt.
324            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
325               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
326         END DO
327      END DO
328      DO jk = 2, jpkfm1
329         DO jj = 2, jpjm1
330            DO ji = fs_2, fs_jpim1   ! vector opt.
331               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
332               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
333            END DO
334         END DO
335      END DO
336      !
337      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
338         DO ji = fs_2, fs_jpim1   ! vector opt.
339            va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
340         END DO
341      END DO
342      DO jk = jpkf-2, 1, -1
343         DO jj = 2, jpjm1   
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
346            END DO
347         END DO
348      END DO
349
350      ! Normalization to obtain the general momentum trend va
351      DO jk = 1, jpkfm1
352         DO jj = 2, jpjm1   
353            DO ji = fs_2, fs_jpim1   ! vector opt.
354               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
355            END DO
356         END DO
357      END DO
358#endif
359      !
360      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
361      !
362      CALL timing_stop('dyn_zdf_imp','section')
363      !
364   END SUBROUTINE dyn_zdf_imp
365
366   !!==============================================================================
367END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.