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

Last change on this file since 456 was 456, checked in by opalod, 18 years ago

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