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

source: branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2132

Last change on this file since 2132 was 2132, checked in by cbricaud, 14 years ago

add change from DEV_r1784_GLS

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