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

source: branches/DEV_r1784_mid_year_merge_2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 1954

Last change on this file since 1954 was 1662, checked in by rblod, 15 years ago

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