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

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

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 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 trdmod          ! ocean dynamics trends
21   USE trdmod_oce      ! ocean variables trends
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC dyn_zdf_imp    ! called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !!   OPA 9.0 , LODYC-IPSL  (2003)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38
39   SUBROUTINE dyn_zdf_imp( kt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dyn_zdf_imp  ***
42      !!                   
43      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
44      !!      and the surface forcing, and add it to the general trend of
45      !!      the momentum equations.
46      !!
47      !! ** Method  :   The vertical momentum mixing trend is given by :
48      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
49      !!      backward time stepping
50      !!      Surface boundary conditions: wind stress input
51      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
52      !!      Add this trend to the general trend ua :
53      !!         ua = ua + dz( avmu dz(u) )
54      !!
55      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
56      !!               mixing trend.
57      !!             - Save the trends in (ztdua,ztdva) ('l_trddyn')
58      !!
59      !! History :
60      !!        !  90-10  (B. Blanke)  Original code
61      !!        !  97-05  (G. Madec)  vertical component of isopycnal
62      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
63      !!   9.0  !  04-08  (C. Talandier)  New trends organization
64      !!---------------------------------------------------------------------
65      !! * Modules used
66      USE oce, ONLY :  zwd   => ta,   &  ! use ta as workspace
67                       zws   => sa       ! use sa as workspace
68
69      !! * Arguments
70      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
71
72      !! * Local declarations
73      INTEGER ::   &
74         ji, jj, jk,                  &  ! dummy loop indices
75         ikbu, ikbum1, ikbv, ikbvm1      ! temporary integers
76      REAL(wp) ::   &
77         zrau0r, z2dt, zua, zva,      &  ! temporary scalars
78         z2dtf, zcoef, zzws, zrhs        !    "         "
79      REAL(wp), DIMENSION(jpi,jpj) ::   &
80         ztsx, ztsy, ztbx, ztby          ! temporary workspace arrays
81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
82         zwi, ztdua, ztdva               ! temporary workspace arrays
83      !!----------------------------------------------------------------------
84
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
89      ENDIF
90
91      ! 0. Local constant initialization
92      ! --------------------------------
93      zrau0r = 1. / rau0      ! inverse of the reference density
94      z2dt   = 2. * rdt       ! Leap-frog environnement
95      ztsx(:,:) = 0.e0
96      ztsy(:,:) = 0.e0 
97      ztbx(:,:) = 0.e0
98      ztby(:,:) = 0.e0
99      ! Euler time stepping when starting from rest
100      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
101
102      ! Save previous ua and va trends
103      IF( l_trddyn )   THEN
104         ztdua(:,:,:) = ua(:,:,:) 
105         ztdva(:,:,:) = va(:,:,:) 
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( l_trddyn )  THEN 
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               ztsx(ji,jj) = 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               ztbx(ji,jj) = - 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               ztdua(ji,jj,1     ) = ztdua(ji,jj,1     ) - ztsx(ji,jj)
206               ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj)
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               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / z2dt
216            END DO
217         END DO
218      END DO
219
220
221      ! 2. Vertical diffusion on v
222      ! ---------------------------
223      ! Matrix and second member construction
224      ! bottom boundary condition: only zws must be masked as avmv can take
225      ! non zero value at the ocean bottom depending on the bottom friction
226      ! used (see zdfmix.F)
227      DO jk = 1, jpkm1
228         DO jj = 2, jpjm1   
229            DO ji = fs_2, fs_jpim1   ! vector opt.
230               zcoef = -z2dt / fse3v(ji,jj,jk)
231               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
232               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
233               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
234               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
235            END DO
236         END DO
237      END DO
238
239      ! Surface boudary conditions
240      DO jj = 2, jpjm1   
241         DO ji = fs_2, fs_jpim1   ! vector opt.
242            zwi(ji,jj,1) = 0.e0
243            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
244         END DO
245      END DO
246
247      ! Matrix inversion
248      !-----------------------------------------------------------------------
249      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
250      !
251      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
252      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
253      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
254      !        (        ...               )( ...  ) ( ...  )
255      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
256      !
257      !   m is decomposed in the product of an upper and lower triangular
258      !   matrix
259      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
260      !   The solution (after velocity) is in 2d array va
261      !-----------------------------------------------------------------------
262
263      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
264      DO jk = 2, jpkm1
265         DO jj = 2, jpjm1   
266            DO ji = fs_2, fs_jpim1   ! vector opt.
267               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
268            END DO
269         END DO
270      END DO
271
272      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
273      DO jj = 2, jpjm1
274         DO ji = fs_2, fs_jpim1   ! vector opt.
275!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
276!!!         va(ji,jj,1) = vb(ji,jj,1)  &
277!!!                      + z2dt * ( va(ji,jj,1) + tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
278            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
279            va(ji,jj,1) = vb(ji,jj,1)  &
280                         + z2dt * va(ji,jj,1) + z2dtf * tauy(ji,jj)
281         END DO
282      END DO
283      DO jk = 2, jpkm1
284         DO jj = 2, jpjm1
285            DO ji = fs_2, fs_jpim1   ! vector opt.
286               zrhs = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)   ! zrhs=right hand side
287               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
288            END DO
289         END DO
290      END DO
291
292      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
293      DO jj = 2, jpjm1   
294         DO ji = fs_2, fs_jpim1   ! vector opt.
295            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
296         END DO
297      END DO
298      DO jk = jpk-2, 1, -1
299         DO jj = 2, jpjm1   
300            DO ji = fs_2, fs_jpim1   ! vector opt.
301               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
302            END DO
303         END DO
304      END DO
305
306      IF( l_trddyn )  THEN 
307         ! diagnose surface and bottom momentum fluxes
308         DO jj = 2, jpjm1   
309            DO ji = fs_2, fs_jpim1   ! vector opt.
310               ! save the surface forcing momentum fluxes
311               ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
312               ! save bottom friction momentum fluxes
313               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
314               ikbvm1 = MAX( ikbv-1, 1 )
315               ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   &
316                  / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
317               ! subtract surface forcing and bottom friction trend from vertical
318               ! diffusive momentum trend
319               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj)
320               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)
321            END DO
322         END DO
323      ENDIF
324
325      ! Normalization to obtain the general momentum trend va
326      DO jk = 1, jpkm1
327         DO jj = 2, jpjm1   
328            DO ji = fs_2, fs_jpim1   ! vector opt.
329               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / z2dt
330            END DO
331         END DO
332      END DO
333
334      ! save the vertical diffusion trends for diagnostic
335      ! momentum trends
336      IF( l_trddyn )  THEN
337         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
338         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
339
340         CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt)
341         ztdua(:,:,:) = 0.e0
342         ztdva(:,:,:) = 0.e0
343         ztdua(:,:,1) = ztsx(:,:)
344         ztdva(:,:,1) = ztsy(:,:)
345         CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt)
346         ztdua(:,:,:) = 0.e0
347         ztdva(:,:,:) = 0.e0
348         ztdua(:,:,1) = ztbx(:,:)
349         ztdva(:,:,1) = ztby(:,:)
350         CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt)
351      ENDIF
352
353      IF(l_ctl) THEN         ! print sum trends (used for debugging)
354         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
355         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
356         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
357         u_ctl = zua   ;   v_ctl = zva
358      ENDIF
359
360   END SUBROUTINE dyn_zdf_imp
361
362   !!==============================================================================
363END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.