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

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 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), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
76      !!----------------------------------------------------------------------
77
78      IF( kt == nit000 ) THEN
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
81         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
82      ENDIF
83
84      ! 0. Local constant initialization
85      ! --------------------------------
86      zrau0r = 1. / rau0      ! inverse of the reference density
87
88      ! 1. Vertical diffusion on u
89      ! ---------------------------
90      ! Matrix and second member construction
91      ! bottom boundary condition: only zws must be masked as avmu can take
92      ! non zero value at the ocean bottom depending on the bottom friction
93      ! used (see zdfmix.F)
94      DO jk = 1, jpkm1
95         DO jj = 2, jpjm1 
96            DO ji = fs_2, fs_jpim1   ! vector opt.
97               zcoef = - p2dt / fse3u(ji,jj,jk)
98               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
99               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
100               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
101               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
102            END DO
103         END DO
104      END DO
105
106      ! Surface boudary conditions
107      DO jj = 2, jpjm1 
108         DO ji = fs_2, fs_jpim1   ! vector opt.
109            zwi(ji,jj,1) = 0.
110            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
111         END DO
112      END DO
113
114      ! Matrix inversion starting from the first level
115      !-----------------------------------------------------------------------
116      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
117      !
118      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
119      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
120      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
121      !        (        ...               )( ...  ) ( ...  )
122      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
123      !
124      !   m is decomposed in the product of an upper and a lower triangular matrix
125      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
126      !   The solution (the after velocity) is in ua
127      !-----------------------------------------------------------------------
128
129      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
130      DO jk = 2, jpkm1
131         DO jj = 2, jpjm1   
132            DO ji = fs_2, fs_jpim1   ! vector opt.
133               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
134            END DO
135         END DO
136      END DO
137
138      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
139      DO jj = 2, jpjm1   
140         DO ji = fs_2, fs_jpim1   ! vector opt.
141!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
142!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
143!!!                      + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
144            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
145            ua(ji,jj,1) = ub(ji,jj,1)  &
146                         + p2dt *  ua(ji,jj,1) + z2dtf * utau(ji,jj)
147         END DO
148      END DO
149      DO jk = 2, jpkm1
150         DO jj = 2, jpjm1   
151            DO ji = fs_2, fs_jpim1   ! vector opt.
152               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
153               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
154            END DO
155         END DO
156      END DO
157
158      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
159      DO jj = 2, jpjm1   
160         DO ji = fs_2, fs_jpim1   ! vector opt.
161            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
162         END DO
163      END DO
164      DO jk = jpk-2, 1, -1
165         DO jj = 2, jpjm1   
166            DO ji = fs_2, fs_jpim1   ! vector opt.
167               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
168            END DO
169         END DO
170      END DO
171
172      ! Normalization to obtain the general momentum trend ua
173      DO jk = 1, jpkm1
174         DO jj = 2, jpjm1   
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
177            END DO
178         END DO
179      END DO
180
181
182      ! 2. Vertical diffusion on v
183      ! ---------------------------
184      ! Matrix and second member construction
185      ! bottom boundary condition: only zws must be masked as avmv can take
186      ! non zero value at the ocean bottom depending on the bottom friction
187      ! used (see zdfmix.F)
188      DO jk = 1, jpkm1
189         DO jj = 2, jpjm1   
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               zcoef = -p2dt / fse3v(ji,jj,jk)
192               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
193               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
194               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
195               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
196            END DO
197         END DO
198      END DO
199
200      ! Surface boudary conditions
201      DO jj = 2, jpjm1   
202         DO ji = fs_2, fs_jpim1   ! vector opt.
203            zwi(ji,jj,1) = 0.e0
204            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
205         END DO
206      END DO
207
208      ! Matrix inversion
209      !-----------------------------------------------------------------------
210      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
211      !
212      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
213      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
214      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
215      !        (        ...               )( ...  ) ( ...  )
216      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
217      !
218      !   m is decomposed in the product of an upper and lower triangular
219      !   matrix
220      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
221      !   The solution (after velocity) is in 2d array va
222      !-----------------------------------------------------------------------
223
224      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
225      DO jk = 2, jpkm1
226         DO jj = 2, jpjm1   
227            DO ji = fs_2, fs_jpim1   ! vector opt.
228               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
229            END DO
230         END DO
231      END DO
232
233      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
234      DO jj = 2, jpjm1
235         DO ji = fs_2, fs_jpim1   ! vector opt.
236!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
237!!!         va(ji,jj,1) = vb(ji,jj,1)  &
238!!!                      + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
239            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
240            va(ji,jj,1) = vb(ji,jj,1)  &
241                         + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj)
242         END DO
243      END DO
244      DO jk = 2, jpkm1
245         DO jj = 2, jpjm1
246            DO ji = fs_2, fs_jpim1   ! vector opt.
247               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
248               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
249            END DO
250         END DO
251      END DO
252
253      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
254      DO jj = 2, jpjm1   
255         DO ji = fs_2, fs_jpim1   ! vector opt.
256            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
257         END DO
258      END DO
259      DO jk = jpk-2, 1, -1
260         DO jj = 2, jpjm1   
261            DO ji = fs_2, fs_jpim1   ! vector opt.
262               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
263            END DO
264         END DO
265      END DO
266
267      ! Normalization to obtain the general momentum trend va
268      DO jk = 1, jpkm1
269         DO jj = 2, jpjm1   
270            DO ji = fs_2, fs_jpim1   ! vector opt.
271               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
272            END DO
273         END DO
274      END DO
275
276   END SUBROUTINE dyn_zdf_imp
277
278   !!==============================================================================
279END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.