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/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk
  • Property svn:keywords set to Id
File size: 12.9 KB
Line 
1MODULE dynzdf_imp
2   !!======================================================================
3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!======================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
8   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
10   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffusion using a implicit time-stepping
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE zdf_oce         ! ocean vertical physics
20   USE phycst          ! physical constants
21   USE in_out_manager  ! I/O manager
22   USE lib_mpp         ! MPP library
23   USE zdfbfr          ! Bottom friction setup
24   USE wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dyn_zdf_imp   ! called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE dyn_zdf_imp( kt, p2dt )
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 (averaged over kt-1/2 & kt+1/2)
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 mixing trend.
59      !!---------------------------------------------------------------------
60      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
61      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
62      !!
63      INTEGER  ::   ji, jj, jk   ! dummy loop indices
64      INTEGER  ::   ikbu, ikbv   ! local integers
65      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
66      !!----------------------------------------------------------------------
67
68      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
69      REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv
70      !!----------------------------------------------------------------------
71      !
72      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
73      !
74      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
75      CALL wrk_alloc( jpi,jpj, zavmu, zavmv ) 
76      !
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
81      ENDIF
82
83      ! 0. Local constant initialization
84      ! --------------------------------
85      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
86
87      ! 1. Apply semi-implicit bottom friction
88      ! --------------------------------------
89      ! Only needed for semi-implicit bottom friction setup. The explicit
90      ! bottom friction has been included in "u(v)a" which act as the R.H.S
91      ! column vector of the tri-diagonal matrix equation
92      !
93
94      IF( ln_bfrimp ) THEN
95# if defined key_vectopt_loop
96      DO jj = 1, 1
97         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
98# else
99      DO jj = 2, jpjm1
100         DO ji = 2, jpim1
101# endif
102            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
103            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
104            zavmu(ji,jj) = avmu(ji,jj,ikbu+1)
105            zavmv(ji,jj) = avmv(ji,jj,ikbv+1)
106            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
107            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
108         END DO
109      END DO
110      ENDIF
111
112      ! 2. Vertical diffusion on u
113      ! ---------------------------
114      ! Matrix and second member construction
115      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
116      ! non zero value at the ocean bottom depending on the bottom friction used.
117      !
118      DO jk = 1, jpkm1        ! Matrix
119         DO jj = 2, jpjm1 
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               zcoef = - p2dt / fse3u(ji,jj,jk)
122               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
123               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
124               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
125               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
126               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
127            END DO
128         END DO
129      END DO
130      DO jj = 2, jpjm1        ! Surface boudary conditions
131         DO ji = fs_2, fs_jpim1   ! vector opt.
132            zwi(ji,jj,1) = 0._wp
133            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
134         END DO
135      END DO
136
137      ! Matrix inversion starting from the first level
138      !-----------------------------------------------------------------------
139      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
140      !
141      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
142      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
143      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
144      !        (        ...               )( ...  ) ( ...  )
145      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
146      !
147      !   m is decomposed in the product of an upper and a lower triangular matrix
148      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
149      !   The solution (the after velocity) is in ua
150      !-----------------------------------------------------------------------
151      !
152      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
153         DO jj = 2, jpjm1   
154            DO ji = fs_2, fs_jpim1   ! vector opt.
155               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
156            END DO
157         END DO
158      END DO
159      !
160      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
161         DO ji = fs_2, fs_jpim1   ! vector opt.
162            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
163               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
164         END DO
165      END DO
166      DO jk = 2, jpkm1
167         DO jj = 2, jpjm1   
168            DO ji = fs_2, fs_jpim1   ! vector opt.
169               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
170               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
171            END DO
172         END DO
173      END DO
174      !
175      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
176         DO ji = fs_2, fs_jpim1   ! vector opt.
177            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
178         END DO
179      END DO
180      DO jk = jpk-2, 1, -1
181         DO jj = 2, jpjm1   
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
184            END DO
185         END DO
186      END DO
187
188      ! Normalization to obtain the general momentum trend ua
189      DO jk = 1, jpkm1
190         DO jj = 2, jpjm1   
191            DO ji = fs_2, fs_jpim1   ! vector opt.
192               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
193            END DO
194         END DO
195      END DO
196
197
198      ! 3. Vertical diffusion on v
199      ! ---------------------------
200      ! Matrix and second member construction
201      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
202      ! non zero value at the ocean bottom depending on the bottom friction used
203      !
204      DO jk = 1, jpkm1        ! Matrix
205         DO jj = 2, jpjm1   
206            DO ji = fs_2, fs_jpim1   ! vector opt.
207               zcoef = -p2dt / fse3v(ji,jj,jk)
208               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
209               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
210               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
211               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
212               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
213            END DO
214         END DO
215      END DO
216      DO jj = 2, jpjm1        ! Surface boudary conditions
217         DO ji = fs_2, fs_jpim1   ! vector opt.
218            zwi(ji,jj,1) = 0._wp
219            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
220         END DO
221      END DO
222
223      ! Matrix inversion
224      !-----------------------------------------------------------------------
225      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
226      !
227      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
228      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
229      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
230      !        (        ...               )( ...  ) ( ...  )
231      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
232      !
233      !   m is decomposed in the product of an upper and lower triangular matrix
234      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
235      !   The solution (after velocity) is in 2d array va
236      !-----------------------------------------------------------------------
237      !
238      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
239         DO jj = 2, jpjm1   
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
242            END DO
243         END DO
244      END DO
245      !
246      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
247         DO ji = fs_2, fs_jpim1   ! vector opt.
248            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
249               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
250         END DO
251      END DO
252      DO jk = 2, jpkm1
253         DO jj = 2, jpjm1
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
256               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
257            END DO
258         END DO
259      END DO
260      !
261      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
262         DO ji = fs_2, fs_jpim1   ! vector opt.
263            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
264         END DO
265      END DO
266      DO jk = jpk-2, 1, -1
267         DO jj = 2, jpjm1   
268            DO ji = fs_2, fs_jpim1   ! vector opt.
269               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
270            END DO
271         END DO
272      END DO
273
274      ! Normalization to obtain the general momentum trend va
275      DO jk = 1, jpkm1
276         DO jj = 2, jpjm1   
277            DO ji = fs_2, fs_jpim1   ! vector opt.
278               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
279            END DO
280         END DO
281      END DO
282
283      !! restore bottom layer avmu(v)
284      IF( ln_bfrimp ) THEN
285# if defined key_vectopt_loop
286      DO jj = 1, 1
287         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
288# else
289      DO jj = 2, jpjm1
290         DO ji = 2, jpim1
291# endif
292            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
293            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
294            avmu(ji,jj,ikbu+1) = zavmu(ji,jj)
295            avmv(ji,jj,ikbv+1) = zavmv(ji,jj)
296         END DO
297      END DO
298      ENDIF
299      !
300      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
301      CALL wrk_dealloc( jpi,jpj, zavmu, zavmv) 
302      !
303      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
304      !
305   END SUBROUTINE dyn_zdf_imp
306
307   !!==============================================================================
308END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.