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_jki.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynzdf_imp_jki.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.5 KB
Line 
1MODULE dynzdf_imp_jki
2   !!==============================================================================
3   !!                    ***  MODULE  dynzdf_imp_jki  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!==============================================================================
6   !! History : 8.5  !  02-08  (G. Madec)  auto-tasking option
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dyn_zdf_imp_jki : update the momentum trend with the vertical
11   !!                     diffusion using an implicit time-stepping and
12   !!                     j-k-i loops.
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE zdf_oce         ! ocean vertical physics
19   USE in_out_manager  ! I/O manager
20   USE taumod          ! surface ocean stress
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC dyn_zdf_imp_jki     ! called by step.F90
27
28   !! * Substitutions
29#  include "domzgr_substitute.h90"
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !!   OPA 9.0 , LOCEAN-IPSL (2005)
33   !! $Header$
34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE dyn_zdf_imp_jki( kt, p2dt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dyn_zdf_imp_jki  ***
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      !!---------------------------------------------------------------------
58      !! * Arguments
59      INTEGER , INTENT( in ) ::   kt                          ! ocean time-step index
60      REAL(wp), INTENT( in ) ::   p2dt                        ! ocean time-step index
61
62      !! * Local declarations
63      INTEGER ::   ji, jj, jk                                 ! dummy loop indices
64      INTEGER ::   ikst, ikenm2, ikstp1                       ! temporary integers
65      REAL(wp) ::  zrau0r, z2dtf, zcoef, zzws                 ! temporary scalars
66      REAL(wp), DIMENSION(jpi,jpk) ::  zwx, zwy, zwz,       & ! workspace
67         &                             zwd, zws, zwi, zwt
68      !!----------------------------------------------------------------------
69
70      IF( kt == nit000 ) THEN
71         IF(lwp) WRITE(numout,*)
72         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp_jki : vertical momentum diffusion implicit operator'
73         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~   auto-task case (j-k-i loop)'
74      ENDIF
75
76
77      ! 0. Local constant initialization
78      ! --------------------------------
79      zrau0r = 1. / rau0      ! inverse of the reference density
80
81      !                                                ! ===============
82      DO jj = 2, jpjm1                                 !  Vertical slab
83         !                                             ! ===============
84         ! 1. Vertical diffusion on u
85         ! ---------------------------
86
87         ! Matrix and second member construction
88         ! bottom boundary condition: only zws must be masked as avmu can take
89         ! non zero value at the ocean bottom depending on the bottom friction
90         ! used (see zdfmix.F)
91         DO jk = 1, jpkm1
92            DO ji = 2, jpim1
93               zcoef = - p2dt / fse3u(ji,jj,jk)
94               zwi(ji,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
95               zzws       = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
96               zws(ji,jk) = zzws * umask(ji,jj,jk+1)
97               zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
98               zwy(ji,jk) = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
99            END DO
100         END DO
101
102         ! Surface boudary conditions
103         DO ji = 2, jpim1
104            z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 )
105            zwi(ji,1) = 0.
106            zwd(ji,1) = 1. - zws(ji,1)
107            zwy(ji,1) = zwy(ji,1) + z2dtf * taux(ji,jj)
108         END DO
109
110         ! Matrix inversion starting from the first level
111         ikst = 1
112!!----------------------------------------------------------------------
113!!         ZDF.MATRIXSOLVER
114!!       ********************
115!!----------------------------------------------------------------------
116!! Matrix inversion
117!   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
118!
119!        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
120!        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
121!        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
122!        (        ...               )( ...  ) ( ...  )
123!        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
124!
125!   m is decomposed in the product of an upper and lower triangular
126!   matrix
127!   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
128!   The second member is in 2d array zwy
129!   The solution is in 2d array zwx
130!   The 2d arry zwt and zwz are work space arrays
131!
132!   N.B. the starting vertical index (ikst) is equal to 1 except for
133!   the resolution of tke matrix where surface tke value is prescribed
134!   so that ikstrt=2.
135!!----------------------------------------------------------------------
136
137         ikstp1 = ikst + 1
138         ikenm2 = jpk - 2
139         DO ji = 2, jpim1
140            zwt(ji,ikst) = zwd(ji,ikst)
141         END DO
142         DO jk = ikstp1, jpkm1
143            DO ji = 2, jpim1
144               zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
145            END DO
146         END DO
147         DO ji = 2, jpim1
148            zwz(ji,ikst) = zwy(ji,ikst)
149         END DO
150         DO jk = ikstp1, jpkm1
151            DO ji = 2, jpim1
152               zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
153            END DO
154         END DO
155         DO ji = 2, jpim1
156            zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
157         END DO
158         DO jk = ikenm2, ikst, -1
159            DO ji = 2, jpim1
160            zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
161            END DO
162         END DO
163     
164         ! Normalization to obtain the general momentum trend ua
165         DO jk = 1, jpkm1
166            DO ji = 2, jpim1
167               ua(ji,jj,jk) = ( zwx(ji,jk) - ub(ji,jj,jk) ) / p2dt
168            END DO
169         END DO
170
171
172         ! 2. Vertical diffusion on v
173         ! ---------------------------
174         ! Matrix and second member construction
175         ! bottom boundary condition: only zws must be masked as avmv can take
176         ! non zero value at the ocean bottom depending on the bottom friction
177         ! used (see zdfmix.F)
178         DO jk = 1, jpkm1
179            DO ji = 2, jpim1
180               zcoef = -p2dt/fse3v(ji,jj,jk)
181               zwi(ji,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
182               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
183               zws(ji,jk) =  zzws * vmask(ji,jj,jk+1)
184               zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
185               zwy(ji,jk) = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
186            END DO
187         END DO
188
189         ! Surface boudary conditions
190         DO ji = 2, jpim1
191            z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 )
192            zwi(ji,1) = 0.e0
193            zwd(ji,1) = 1. - zws(ji,1)
194            zwy(ji,1) = zwy(ji,1) + z2dtf * tauy(ji,jj)
195         END DO
196
197         ! Matrix inversion starting from the first level
198         ikst = 1
199!!----------------------------------------------------------------------
200!!         ZDF.MATRIXSOLVER
201!!       ********************
202!!----------------------------------------------------------------------
203!! Matrix inversion
204!   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
205!
206!        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
207!        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
208!        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
209!        (        ...               )( ...  ) ( ...  )
210!        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
211!
212!   m is decomposed in the product of an upper and lower triangular
213!   matrix
214!   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
215!   The second member is in 2d array zwy
216!   The solution is in 2d array zwx
217!   The 2d arry zwt and zwz are work space arrays
218!
219!   N.B. the starting vertical index (ikst) is equal to 1 except for
220!   the resolution of tke matrix where surface tke value is prescribed
221!   so that ikstrt=2.
222!!----------------------------------------------------------------------
223
224         ikstp1 = ikst + 1
225         ikenm2 = jpk - 2
226         DO ji = 2, jpim1
227            zwt(ji,ikst) = zwd(ji,ikst)
228         END DO
229         DO jk = ikstp1, jpkm1
230            DO ji = 2, jpim1
231               zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
232            END DO
233         END DO
234         DO ji = 2, jpim1
235            zwz(ji,ikst) = zwy(ji,ikst)
236         END DO
237         DO jk = ikstp1, jpkm1
238            DO ji = 2, jpim1
239               zwz(ji,jk) = zwy(ji,jk) - zwi(ji,jk) / zwt(ji,jk-1) * zwz(ji,jk-1)
240            END DO
241         END DO
242         DO ji = 2, jpim1
243            zwx(ji,jpkm1) = zwz(ji,jpkm1) / zwt(ji,jpkm1)
244         END DO
245         DO jk = ikenm2, ikst, -1
246            DO ji = 2, jpim1
247               zwx(ji,jk) =( zwz(ji,jk) - zws(ji,jk) * zwx(ji,jk+1) ) / zwt(ji,jk)
248            END DO
249         END DO
250     
251         ! Normalization to obtain the general momentum trend va
252         DO jk = 1, jpkm1
253            DO ji = 2, jpim1
254               va(ji,jj,jk) = ( zwx(ji,jk) - vb(ji,jj,jk) ) / p2dt
255            END DO
256         END DO
257         !                                             ! ===============
258      END DO                                           !   End of slab
259      !                                                ! ===============
260
261   END SUBROUTINE dyn_zdf_imp_jki
262
263   !!==============================================================================
264END MODULE dynzdf_imp_jki
Note: See TracBrowser for help on using the repository browser.