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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2450

Last change on this file since 2450 was 2450, checked in by gm, 14 years ago

v3.3beta: #766 share the deepest ocean level indices continuaton
  • Property svn:keywords set to Id
File size: 12.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     1.0  !  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 diffu-
14   !!                  sion using a implicit time-stepping.
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE zdf_oce         ! ocean vertical physics
20   USE phycst          ! physical constants
21   USE in_out_manager  ! I/O manager
22#if defined key_zdfgls
23   USE zdfbfr, ONLY : bfrua, bfrva, wbotu, wbotv ! bottom stresses
24   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
25#endif
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dyn_zdf_imp   ! called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dyn_zdf_imp( kt, p2dt )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dyn_zdf_imp  ***
46      !!                   
47      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
48      !!      and the surface forcing, and add it to the general trend of
49      !!      the momentum equations.
50      !!
51      !! ** Method  :   The vertical momentum mixing trend is given by :
52      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
53      !!      backward time stepping
54      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
55      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
56      !!      Add this trend to the general trend ua :
57      !!         ua = ua + dz( avmu dz(u) )
58      !!
59      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
60      !!---------------------------------------------------------------------
61      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
62      USE oce, ONLY :  zws   => sa      ! use sa as workspace
63      !!
64      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
65      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
66      !!
67      INTEGER  ::   ji, jj, jk             ! dummy loop indices
68      REAL(wp) ::   zrau0r, zcoef         ! temporary scalars
69      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
70      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
71      !!----------------------------------------------------------------------
72
73      IF( kt == nit000 ) THEN
74         IF(lwp) WRITE(numout,*)
75         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
76         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
77      ENDIF
78
79      ! 0. Local constant initialization
80      ! --------------------------------
81      zrau0r = 1. / rau0      ! inverse of the reference density
82
83      ! 1. Vertical diffusion on u
84      ! ---------------------------
85      ! Matrix and second member construction
86      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
87      ! non zero value at the ocean bottom depending on the bottom friction
88      ! used but the bottom velocities have already been updated with the bottom
89      ! friction velocity in dyn_bfr using values from the previous timestep. There
90      ! is no need to include these in the implicit calculation.
91      !
92      DO jk = 1, jpkm1        ! Matrix
93         DO jj = 2, jpjm1 
94            DO ji = fs_2, fs_jpim1   ! vector opt.
95               zcoef = - p2dt / fse3u(ji,jj,jk)
96               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
97               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
98               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
99               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
100               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
101            END DO
102         END DO
103      END DO
104      DO jj = 2, jpjm1        ! Surface boudary conditions
105         DO ji = fs_2, fs_jpim1   ! vector opt.
106            zwi(ji,jj,1) = 0.
107            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
108         END DO
109      END DO
110
111      ! Matrix inversion starting from the first level
112      !-----------------------------------------------------------------------
113      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
114      !
115      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
116      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
117      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
118      !        (        ...               )( ...  ) ( ...  )
119      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
120      !
121      !   m is decomposed in the product of an upper and a lower triangular matrix
122      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
123      !   The solution (the after velocity) is in ua
124      !-----------------------------------------------------------------------
125      !
126      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
127         DO jj = 2, jpjm1   
128            DO ji = fs_2, fs_jpim1   ! vector opt.
129               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
130            END DO
131         END DO
132      END DO
133      !
134      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
135         DO ji = fs_2, fs_jpim1   ! vector opt.
136            ua(ji,jj,1) = ub(ji,jj,1)   &
137               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
138         END DO
139      END DO
140      DO jk = 2, jpkm1
141         DO jj = 2, jpjm1   
142            DO ji = fs_2, fs_jpim1   ! vector opt.
143               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
144               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
145            END DO
146         END DO
147      END DO
148      !
149      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
150         DO ji = fs_2, fs_jpim1   ! vector opt.
151            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
152         END DO
153      END DO
154      DO jk = jpk-2, 1, -1
155         DO jj = 2, jpjm1   
156            DO ji = fs_2, fs_jpim1   ! vector opt.
157               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
158            END DO
159         END DO
160      END DO
161
162#if defined key_zdfgls
163      ! Save bottom stress for next time step
164      DO jj = 2, jpjm1
165         DO ji = fs_2, fs_jpim1   ! vector opt.
166            wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,mbku(ji,jj)) * umask(ji,jj,mbku(ji,jj))
167         END DO
168      END DO
169      CALL lbc_lnk( wbotu(:,:), 'U', -1. )
170#endif
171
172      ! Normalization to obtain the general momentum trend ua
173      DO jk = 1, jpkm1
174         DO jj = 2, jpjm1   
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
177            END DO
178         END DO
179      END DO
180
181
182      ! 2. Vertical diffusion on v
183      ! ---------------------------
184      ! Matrix and second member construction
185      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
186      ! non zero value at the ocean bottom depending on the bottom friction
187      ! used but the bottom velocities have already been updated with the bottom
188      ! friction velocity in dyn_bfr using values from the previous timestep. There
189      ! is no need to include these in the implicit calculation.
190      !
191      DO jk = 1, jpkm1        ! Matrix
192         DO jj = 2, jpjm1   
193            DO ji = fs_2, fs_jpim1   ! vector opt.
194               zcoef = -p2dt / fse3v(ji,jj,jk)
195               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
196               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
197               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
198               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
199               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
200            END DO
201         END DO
202      END DO
203      DO jj = 2, jpjm1        ! Surface boudary conditions
204         DO ji = fs_2, fs_jpim1   ! vector opt.
205            zwi(ji,jj,1) = 0.e0
206            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
207         END DO
208      END DO
209
210      ! Matrix inversion
211      !-----------------------------------------------------------------------
212      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
213      !
214      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
215      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
216      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
217      !        (        ...               )( ...  ) ( ...  )
218      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
219      !
220      !   m is decomposed in the product of an upper and lower triangular matrix
221      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
222      !   The solution (after velocity) is in 2d array va
223      !-----------------------------------------------------------------------
224      !
225      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
226         DO jj = 2, jpjm1   
227            DO ji = fs_2, fs_jpim1   ! vector opt.
228               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
229            END DO
230         END DO
231      END DO
232      !
233      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
234         DO ji = fs_2, fs_jpim1   ! vector opt.
235            va(ji,jj,1) = vb(ji,jj,1)   &
236               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
237         END DO
238      END DO
239      DO jk = 2, jpkm1
240         DO jj = 2, jpjm1
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
243               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
244            END DO
245         END DO
246      END DO
247      !
248      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
249         DO ji = fs_2, fs_jpim1   ! vector opt.
250            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
251         END DO
252      END DO
253      DO jk = jpk-2, 1, -1
254         DO jj = 2, jpjm1   
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
257            END DO
258         END DO
259      END DO
260
261#if defined key_zdfgls
262      ! Save bottom stress for next time step
263      DO jj = 2, jpjm1
264         DO ji = fs_2, fs_jpim1   ! vector opt.
265            wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,mbku(ji,jj))
266         END DO
267      END DO
268      CALL lbc_lnk( wbotv(:,:), 'V', -1. )
269#endif
270
271      ! Normalization to obtain the general momentum trend va
272      DO jk = 1, jpkm1
273         DO jj = 2, jpjm1   
274            DO ji = fs_2, fs_jpim1   ! vector opt.
275               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
276            END DO
277         END DO
278      END DO
279      !
280   END SUBROUTINE dyn_zdf_imp
281
282   !!==============================================================================
283END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.