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

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

nemo_v1_update_048:RB: reorganization of dynamics part, in addition change atsk to jki, suppress dynhpg_atsk.F90 dynzdf_imp_atsk.F90 dynzdf_iso.F90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.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 , LOCEAN-IPSL (2005)
12   !! $Header$
13   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE phycst          ! physical constants
19   USE zdf_oce         ! ocean vertical physics
20   USE in_out_manager  ! I/O manager
21   USE taumod          ! surface ocean stress
22   USE prtctl          ! Print control
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC dyn_zdf_imp    ! called by step.F90
29
30   !! * Substitutions
31#  include "domzgr_substitute.h90"
32#  include "vectopt_loop_substitute.h90"
33   !!----------------------------------------------------------------------
34   !!   OPA 9.0 , LOCEAN-IPSL (2005)
35   !! $Header$
36   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41
42   SUBROUTINE dyn_zdf_imp( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dyn_zdf_imp  ***
45      !!                   
46      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
47      !!      and the surface forcing, and add it to the general trend of
48      !!      the momentum equations.
49      !!
50      !! ** Method  :   The vertical momentum mixing trend is given by :
51      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
52      !!      backward time stepping
53      !!      Surface boundary conditions: wind stress input
54      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
55      !!      Add this trend to the general trend ua :
56      !!         ua = ua + dz( avmu dz(u) )
57      !!
58      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive
59      !!               mixing trend.
60      !!
61      !! History :
62      !!        !  90-10  (B. Blanke)  Original code
63      !!        !  97-05  (G. Madec)  vertical component of isopycnal
64      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
65      !!   9.0  !  04-08  (C. Talandier)  New trends organization
66      !!---------------------------------------------------------------------
67      !! * Modules used
68      USE oce, ONLY :  zwd   => ta,   &  ! use ta as workspace
69                       zws   => sa       ! use sa as workspace
70
71      !! * Arguments
72      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
73
74      !! * Local declarations
75      INTEGER ::   ji, jj, jk            ! dummy loop indices
76      REAL(wp) ::   &
77         zrau0r, z2dt,                &  ! temporary scalars
78         z2dtf, zcoef, zzws, zrhs        !    "         "
79      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
80         zwi                             ! temporary workspace arrays
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
94      ! Euler time stepping when starting from rest
95      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
96
97      ! 1. Vertical diffusion on u
98      ! ---------------------------
99      ! Matrix and second member construction
100      ! bottom boundary condition: only zws must be masked as avmu can take
101      ! non zero value at the ocean bottom depending on the bottom friction
102      ! used (see zdfmix.F)
103      DO jk = 1, jpkm1
104         DO jj = 2, jpjm1 
105            DO ji = fs_2, fs_jpim1   ! vector opt.
106               zcoef = - z2dt / fse3u(ji,jj,jk)
107               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
108               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
109               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
110               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
111            END DO
112         END DO
113      END DO
114
115      ! Surface boudary conditions
116      DO jj = 2, jpjm1 
117         DO ji = fs_2, fs_jpim1   ! vector opt.
118            zwi(ji,jj,1) = 0.
119            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
120         END DO
121      END DO
122
123      ! Matrix inversion starting from the first level
124      !-----------------------------------------------------------------------
125      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
126      !
127      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
128      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
129      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
130      !        (        ...               )( ...  ) ( ...  )
131      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
132      !
133      !   m is decomposed in the product of an upper and a lower triangular matrix
134      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
135      !   The solution (the after velocity) is in ua
136      !-----------------------------------------------------------------------
137
138      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
139      DO jk = 2, jpkm1
140         DO jj = 2, jpjm1   
141            DO ji = fs_2, fs_jpim1   ! vector opt.
142               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
143            END DO
144         END DO
145      END DO
146
147      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
148      DO jj = 2, jpjm1   
149         DO ji = fs_2, fs_jpim1   ! vector opt.
150!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
151!!!         ua(ji,jj,1) = ub(ji,jj,1)  &
152!!!                      + z2dt * ( ua(ji,jj,1) + taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) )
153            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
154            ua(ji,jj,1) = ub(ji,jj,1)  &
155                         + z2dt *  ua(ji,jj,1) + z2dtf * taux(ji,jj)
156         END DO
157      END DO
158      DO jk = 2, jpkm1
159         DO jj = 2, jpjm1   
160            DO ji = fs_2, fs_jpim1   ! vector opt.
161               zrhs = ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)   ! zrhs=right hand side
162               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
163            END DO
164         END DO
165      END DO
166
167      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
168      DO jj = 2, jpjm1   
169         DO ji = fs_2, fs_jpim1   ! vector opt.
170            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
171         END DO
172      END DO
173      DO jk = jpk-2, 1, -1
174         DO jj = 2, jpjm1   
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180
181      ! Normalization to obtain the general momentum trend ua
182      DO jk = 1, jpkm1
183         DO jj = 2, jpjm1   
184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / z2dt
186            END DO
187         END DO
188      END DO
189
190
191      ! 2. Vertical diffusion on v
192      ! ---------------------------
193      ! Matrix and second member construction
194      ! bottom boundary condition: only zws must be masked as avmv can take
195      ! non zero value at the ocean bottom depending on the bottom friction
196      ! used (see zdfmix.F)
197      DO jk = 1, jpkm1
198         DO jj = 2, jpjm1   
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               zcoef = -z2dt / fse3v(ji,jj,jk)
201               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
202               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
203               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
204               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
205            END DO
206         END DO
207      END DO
208
209      ! Surface boudary conditions
210      DO jj = 2, jpjm1   
211         DO ji = fs_2, fs_jpim1   ! vector opt.
212            zwi(ji,jj,1) = 0.e0
213            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
214         END DO
215      END DO
216
217      ! Matrix inversion
218      !-----------------------------------------------------------------------
219      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
220      !
221      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
222      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
223      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
224      !        (        ...               )( ...  ) ( ...  )
225      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
226      !
227      !   m is decomposed in the product of an upper and lower triangular
228      !   matrix
229      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
230      !   The solution (after velocity) is in 2d array va
231      !-----------------------------------------------------------------------
232
233      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
234      DO jk = 2, jpkm1
235         DO jj = 2, jpjm1   
236            DO ji = fs_2, fs_jpim1   ! vector opt.
237               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
238            END DO
239         END DO
240      END DO
241
242      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
243      DO jj = 2, jpjm1
244         DO ji = fs_2, fs_jpim1   ! vector opt.
245!!! change les resultats (derniers digit, pas significativement + rapide 1* de moins)
246!!!         va(ji,jj,1) = vb(ji,jj,1)  &
247!!!                      + z2dt * ( va(ji,jj,1) + tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) )
248            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
249            va(ji,jj,1) = vb(ji,jj,1)  &
250                         + z2dt * va(ji,jj,1) + z2dtf * tauy(ji,jj)
251         END DO
252      END DO
253      DO jk = 2, jpkm1
254         DO jj = 2, jpjm1
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               zrhs = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)   ! zrhs=right hand side
257               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
258            END DO
259         END DO
260      END DO
261
262      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
263      DO jj = 2, jpjm1   
264         DO ji = fs_2, fs_jpim1   ! vector opt.
265            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
266         END DO
267      END DO
268      DO jk = jpk-2, 1, -1
269         DO jj = 2, jpjm1   
270            DO ji = fs_2, fs_jpim1   ! vector opt.
271               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
272            END DO
273         END DO
274      END DO
275
276! flux de surface doit etre calcule dans trdmod et boootom stress
277! deduit par integration verticale dans trmod pour jpdtdzdf
278!RB      IF( l_trddyn )  THEN
279!         ! diagnose surface and bottom momentum fluxes
280!         DO jj = 2, jpjm1   
281!            DO ji = fs_2, fs_jpim1   ! vector opt.
282!               ! save the surface forcing momentum fluxes
283!               ztsy(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
284!               ztsx(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
285!               ! save bottom friction momentum fluxes
286!               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
287!               ikbvm1 = MAX( ikbv-1, 1 )
288!               ztby(ji,jj) = - avmv(ji,jj,ikbv) * va(ji,jj,ikbvm1)   &
289!                  / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
290!               ! subtract surface forcing and bottom friction trend from vertical
291!               ! diffusive momentum trend
292!               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj)
293!               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)
294!            END DO
295!         END DO
296!      ENDIF
297
298      ! Normalization to obtain the general momentum trend va
299      DO jk = 1, jpkm1
300         DO jj = 2, jpjm1   
301            DO ji = fs_2, fs_jpim1   ! vector opt.
302               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / z2dt
303            END DO
304         END DO
305      END DO
306
307   END SUBROUTINE dyn_zdf_imp
308
309   !!==============================================================================
310END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.