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

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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