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

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

update licence of all NEMO files...
  • Property svn:keywords set to Id
File size: 12.3 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#if defined key_zdfgls
72      INTEGER :: ikbu, ikbv, ikbum1, ikbvm1
73      REAL(wp) :: zcbcu, zcbcv
74#endif
75      !!----------------------------------------------------------------------
76
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
81      ENDIF
82
83      ! 0. Local constant initialization
84      ! --------------------------------
85      zrau0r = 1. / rau0      ! inverse of the reference density
86
87      ! 1. Vertical diffusion on u
88      ! ---------------------------
89      ! Matrix and second member construction
90      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
91      ! non zero value at the ocean bottom depending on the bottom friction
92      ! used but the bottom velocities have already been updated with the bottom
93      ! friction velocity in dyn_bfr using values from the previous timestep. There
94      ! is no need to include these in the implicit calculation.
95      !
96      DO jk = 1, jpkm1        ! Matrix
97         DO jj = 2, jpjm1 
98            DO ji = fs_2, fs_jpim1   ! vector opt.
99               zcoef = - p2dt / fse3u(ji,jj,jk)
100               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
101               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
102               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
103               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
104               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
105            END DO
106         END DO
107      END DO
108      DO jj = 2, jpjm1        ! Surface boudary conditions
109         DO ji = fs_2, fs_jpim1   ! vector opt.
110            zwi(ji,jj,1) = 0.
111            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
112         END DO
113      END DO
114
115      ! Matrix inversion starting from the first level
116      !-----------------------------------------------------------------------
117      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
118      !
119      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
120      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
121      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
122      !        (        ...               )( ...  ) ( ...  )
123      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
124      !
125      !   m is decomposed in the product of an upper and a lower triangular matrix
126      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
127      !   The solution (the after velocity) is in ua
128      !-----------------------------------------------------------------------
129      !
130      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
131         DO jj = 2, jpjm1   
132            DO ji = fs_2, fs_jpim1   ! vector opt.
133               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
134            END DO
135         END DO
136      END DO
137      !
138      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
139         DO ji = fs_2, fs_jpim1   ! vector opt.
140            ua(ji,jj,1) = ub(ji,jj,1)   &
141               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
142         END DO
143      END DO
144      DO jk = 2, jpkm1
145         DO jj = 2, jpjm1   
146            DO ji = fs_2, fs_jpim1   ! vector opt.
147               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
148               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
149            END DO
150         END DO
151      END DO
152      !
153      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
154         DO ji = fs_2, fs_jpim1   ! vector opt.
155            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
156         END DO
157      END DO
158      DO jk = jpk-2, 1, -1
159         DO jj = 2, jpjm1   
160            DO ji = fs_2, fs_jpim1   ! vector opt.
161               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
162            END DO
163         END DO
164      END DO
165
166#if defined key_zdfgls
167      ! Save bottom stress for next time step
168      DO jj = 2, jpjm1
169         DO ji = fs_2, fs_jpim1   ! vector opt.
170            ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
171            ikbum1 = MAX( ikbu-1, 1 )
172            wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1)
173         END DO
174      END DO
175      CALL lbc_lnk( wbotu(:,:), 'U', -1. )
176#endif
177
178      ! Normalization to obtain the general momentum trend ua
179      DO jk = 1, jpkm1
180         DO jj = 2, jpjm1   
181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
183            END DO
184         END DO
185      END DO
186
187
188      ! 2. Vertical diffusion on v
189      ! ---------------------------
190      ! Matrix and second member construction
191      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
192      ! non zero value at the ocean bottom depending on the bottom friction
193      ! used but the bottom velocities have already been updated with the bottom
194      ! friction velocity in dyn_bfr using values from the previous timestep. There
195      ! is no need to include these in the implicit calculation.
196      !
197      DO jk = 1, jpkm1        ! Matrix
198         DO jj = 2, jpjm1   
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               zcoef = -p2dt / fse3v(ji,jj,jk)
201               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
202               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
203               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
204               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
205               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
206            END DO
207         END DO
208      END DO
209      DO jj = 2, jpjm1        ! Surface boudary conditions
210         DO ji = fs_2, fs_jpim1   ! vector opt.
211            zwi(ji,jj,1) = 0.e0
212            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
213         END DO
214      END DO
215
216      ! Matrix inversion
217      !-----------------------------------------------------------------------
218      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
219      !
220      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
221      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
222      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
223      !        (        ...               )( ...  ) ( ...  )
224      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
225      !
226      !   m is decomposed in the product of an upper and lower triangular matrix
227      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
228      !   The solution (after velocity) is in 2d array va
229      !-----------------------------------------------------------------------
230      !
231      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
232         DO jj = 2, jpjm1   
233            DO ji = fs_2, fs_jpim1   ! vector opt.
234               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
235            END DO
236         END DO
237      END DO
238      !
239      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
240         DO ji = fs_2, fs_jpim1   ! vector opt.
241            va(ji,jj,1) = vb(ji,jj,1)   &
242               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
243         END DO
244      END DO
245      DO jk = 2, jpkm1
246         DO jj = 2, jpjm1
247            DO ji = fs_2, fs_jpim1   ! vector opt.
248               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
249               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
250            END DO
251         END DO
252      END DO
253      !
254      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
255         DO ji = fs_2, fs_jpim1   ! vector opt.
256            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
257         END DO
258      END DO
259      DO jk = jpk-2, 1, -1
260         DO jj = 2, jpjm1   
261            DO ji = fs_2, fs_jpim1   ! vector opt.
262               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
263            END DO
264         END DO
265      END DO
266
267#if defined key_zdfgls
268      ! Save bottom stress for next time step
269      DO jj = 2, jpjm1
270         DO ji = fs_2, fs_jpim1   ! vector opt.
271            ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
272            ikbvm1 = MAX( ikbv-1, 1 )
273            wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1)
274         END DO
275      END DO
276      CALL lbc_lnk( wbotv(:,:), 'V', -1. )
277#endif
278
279      ! Normalization to obtain the general momentum trend va
280      DO jk = 1, jpkm1
281         DO jj = 2, jpjm1   
282            DO ji = fs_2, fs_jpim1   ! vector opt.
283               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
284            END DO
285         END DO
286      END DO
287      !
288   END SUBROUTINE dyn_zdf_imp
289
290   !!==============================================================================
291END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.