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

Last change on this file since 3211 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

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