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

source: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.7 KB
Line 
1MODULE dynzdf_imp
2   !!==============================================================================
3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffu-
9   !!                  sion using a implicit time-stepping.
10   !!----------------------------------------------------------------------
11   !!   OPA 9.0 , LOCEAN-IPSL (2005)
12   !! $Header$
13   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE zdf_oce         ! ocean vertical physics
20   USE in_out_manager  ! I/O manager
21   USE taumod          ! surface ocean stress
22   USE trdmod          ! ocean dynamics trends
23   USE trdmod_oce      ! ocean variables trends
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC dyn_zdf_imp    ! called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LOCEAN-IPSL (2005)
36   !! $Header$
37   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42
43   SUBROUTINE dyn_zdf_imp( kt )
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
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
60      !!               mixing trend.
61      !!             - Save the trends in (ztdua,ztdva) ('l_trddyn')
62      !!
63      !! History :
64      !!        !  90-10  (B. Blanke)  Original code
65      !!        !  97-05  (G. Madec)  vertical component of isopycnal
66      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
67      !!   9.0  !  04-08  (C. Talandier)  New trends organization
68      !!---------------------------------------------------------------------
69      !! * Modules used
70      USE oce, ONLY :  zwd   => ta,   &  ! use ta as workspace
71                       zws   => sa       ! use sa as workspace
72
73      !! * Arguments
74      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
75
76      !! * Local declarations
77      INTEGER ::   &
78         ji, jj, jk,                  &  ! dummy loop indices
79         ikbu, ikbum1, ikbv, ikbvm1      ! temporary integers
80      REAL(wp) ::   &
81         zrau0r, z2dt, zua, zva,      &  ! temporary scalars
82         z2dtf, zcoef, zzws, zrhs        !    "         "
83      REAL(wp), DIMENSION(jpi,jpj) ::   &
84         ztsx, ztsy, ztbx, ztby          ! temporary workspace arrays
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
86         zwi, ztdua, ztdva               ! temporary workspace arrays
87      !!----------------------------------------------------------------------
88
89      IF( kt == nit000 ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
93      ENDIF
94
95      ! 0. Local constant initialization
96      ! --------------------------------
97      zrau0r = 1. / rau0      ! inverse of the reference density
98      z2dt   = 2. * rdt       ! Leap-frog environnement
99      ztsx(:,:) = 0.e0
100      ztsy(:,:) = 0.e0 
101      ztbx(:,:) = 0.e0
102      ztby(:,:) = 0.e0
103      ! Euler time stepping when starting from rest
104      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
105
106      ! Save previous ua and va trends
107      IF( l_trddyn )   THEN
108         ztdua(:,:,:) = ua(:,:,:) 
109         ztdva(:,:,:) = va(:,:,:) 
110      ENDIF
111
112      ! 1. Vertical diffusion on u
113      ! ---------------------------
114      ! Matrix and second member construction
115      ! bottom boundary condition: only zws must be masked as avmu can take
116      ! non zero value at the ocean bottom depending on the bottom friction
117      ! used (see zdfmix.F)
118      DO jk = 1, jpkm1
119         DO jj = 2, jpjm1 
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               zcoef = - z2dt / fse3u(ji,jj,jk)
122               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(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. - zwi(ji,jj,jk) - zzws
126            END DO
127         END DO
128      END DO
129
130      ! Surface boudary conditions
131      DO jj = 2, jpjm1 
132         DO ji = fs_2, fs_jpim1   ! vector opt.
133            zwi(ji,jj,1) = 0.
134            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
135         END DO
136      END DO
137
138      ! Matrix inversion starting from the first level
139      !-----------------------------------------------------------------------
140      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
141      !
142      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
143      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
144      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
145      !        (        ...               )( ...  ) ( ...  )
146      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
147      !
148      !   m is decomposed in the product of an upper and a lower triangular matrix
149      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
150      !   The solution (the after velocity) is in ua
151      !-----------------------------------------------------------------------
152
153      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
154      DO jk = 2, jpkm1
155         DO jj = 2, jpjm1   
156            DO ji = fs_2, fs_jpim1   ! vector opt.
157               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
158            END DO
159         END DO
160      END DO
161
162      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
163      DO jj = 2, jpjm1   
164         DO ji = fs_2, fs_jpim1   ! vector opt.
165!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
166!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
167!!!                      + z2dt * ( ua(ji,jj,1) + taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
168            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
169            ua(ji,jj,1) = ub(ji,jj,1)  &
170                         + z2dt *  ua(ji,jj,1) + z2dtf * taux(ji,jj)
171         END DO
172      END DO
173      DO jk = 2, jpkm1
174         DO jj = 2, jpjm1   
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               zrhs = ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)   ! zrhs=right hand side
177               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
178            END DO
179         END DO
180      END DO
181
182      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
183      DO jj = 2, jpjm1   
184         DO ji = fs_2, fs_jpim1   ! vector opt.
185            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
186         END DO
187      END DO
188      DO jk = jpk-2, 1, -1
189         DO jj = 2, jpjm1   
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
192            END DO
193         END DO
194      END DO
195
196      IF( l_trddyn )  THEN 
197         ! diagnose surface and bottom momentum fluxes
198         DO jj = 2, jpjm1   
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               ! save the surface forcing momentum fluxes
201               ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
202               ! save bottom friction momentum fluxes
203               ikbu   = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
204               ikbum1 = MAX( ikbu-1, 1 )
205               ztbx(ji,jj) = - avmu(ji,jj,ikbu) * ua(ji,jj,ikbum1)   &
206                  / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
207               ! subtract surface forcing and bottom friction trend from vertical
208               ! diffusive momentum trend
209               ztdua(ji,jj,1     ) = ztdua(ji,jj,1     ) - ztsx(ji,jj)
210               ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj)
211            END DO
212         END DO
213      ENDIF
214
215      ! Normalization to obtain the general momentum trend ua
216      DO jk = 1, jpkm1
217         DO jj = 2, jpjm1   
218            DO ji = fs_2, fs_jpim1   ! vector opt.
219               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / z2dt
220            END DO
221         END DO
222      END DO
223
224
225      ! 2. Vertical diffusion on v
226      ! ---------------------------
227      ! Matrix and second member construction
228      ! bottom boundary condition: only zws must be masked as avmv can take
229      ! non zero value at the ocean bottom depending on the bottom friction
230      ! used (see zdfmix.F)
231      DO jk = 1, jpkm1
232         DO jj = 2, jpjm1   
233            DO ji = fs_2, fs_jpim1   ! vector opt.
234               zcoef = -z2dt / fse3v(ji,jj,jk)
235               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
236               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
237               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
238               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
239            END DO
240         END DO
241      END DO
242
243      ! Surface boudary conditions
244      DO jj = 2, jpjm1   
245         DO ji = fs_2, fs_jpim1   ! vector opt.
246            zwi(ji,jj,1) = 0.e0
247            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
248         END DO
249      END DO
250
251      ! Matrix inversion
252      !-----------------------------------------------------------------------
253      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
254      !
255      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
256      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
257      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
258      !        (        ...               )( ...  ) ( ...  )
259      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
260      !
261      !   m is decomposed in the product of an upper and lower triangular
262      !   matrix
263      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
264      !   The solution (after velocity) is in 2d array va
265      !-----------------------------------------------------------------------
266
267      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
268      DO jk = 2, jpkm1
269         DO jj = 2, jpjm1   
270            DO ji = fs_2, fs_jpim1   ! vector opt.
271               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
272            END DO
273         END DO
274      END DO
275
276      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
277      DO jj = 2, jpjm1
278         DO ji = fs_2, fs_jpim1   ! vector opt.
279!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
280!!!         va(ji,jj,1) = vb(ji,jj,1)  &
281!!!                      + z2dt * ( va(ji,jj,1) + tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
282            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
283            va(ji,jj,1) = vb(ji,jj,1)  &
284                         + z2dt * va(ji,jj,1) + z2dtf * tauy(ji,jj)
285         END DO
286      END DO
287      DO jk = 2, jpkm1
288         DO jj = 2, jpjm1
289            DO ji = fs_2, fs_jpim1   ! vector opt.
290               zrhs = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)   ! zrhs=right hand side
291               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
292            END DO
293         END DO
294      END DO
295
296      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
297      DO jj = 2, jpjm1   
298         DO ji = fs_2, fs_jpim1   ! vector opt.
299            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
300         END DO
301      END DO
302      DO jk = jpk-2, 1, -1
303         DO jj = 2, jpjm1   
304            DO ji = fs_2, fs_jpim1   ! vector opt.
305               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
306            END DO
307         END DO
308      END DO
309
310      IF( l_trddyn )  THEN 
311         ! diagnose surface and bottom momentum fluxes
312         DO jj = 2, jpjm1   
313            DO ji = fs_2, fs_jpim1   ! vector opt.
314               ! save the surface forcing momentum fluxes
315               ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
316               ! save bottom friction momentum fluxes
317               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
318               ikbvm1 = MAX( ikbv-1, 1 )
319               ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   &
320                  / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
321               ! subtract surface forcing and bottom friction trend from vertical
322               ! diffusive momentum trend
323               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj)
324               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)
325            END DO
326         END DO
327      ENDIF
328
329      ! Normalization to obtain the general momentum trend va
330      DO jk = 1, jpkm1
331         DO jj = 2, jpjm1   
332            DO ji = fs_2, fs_jpim1   ! vector opt.
333               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / z2dt
334            END DO
335         END DO
336      END DO
337
338      ! save the vertical diffusion trends for diagnostic
339      ! momentum trends
340      IF( l_trddyn )  THEN
341         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
342         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
343
344         CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt)
345         ztdua(:,:,:) = 0.e0
346         ztdva(:,:,:) = 0.e0
347         ztdua(:,:,1) = ztsx(:,:)
348         ztdva(:,:,1) = ztsy(:,:)
349         CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt)
350         ztdua(:,:,:) = 0.e0
351         ztdva(:,:,:) = 0.e0
352         ztdua(:,:,1) = ztbx(:,:)
353         ztdva(:,:,1) = ztby(:,:)
354         CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt)
355      ENDIF
356
357      IF(l_ctl) THEN         ! print sum trends (used for debugging)
358         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
359         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
360         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
361         u_ctl = zua   ;   v_ctl = zva
362      ENDIF
363
364   END SUBROUTINE dyn_zdf_imp
365
366   !!==============================================================================
367END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.