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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 3278

Last change on this file since 3278 was 3278, checked in by hliu, 12 years ago

keep avmu(v) intact in dynzdf_imp.F90, ref: ticket #918
  • Property svn:keywords set to Id
File size: 12.4 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, zavmu, zavmv
69      !!----------------------------------------------------------------------
70      !
71      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
72      !
73      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws, zavmu, zavmv ) 
74      !
75      IF( kt == nit000 ) THEN
76         IF(lwp) WRITE(numout,*)
77         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
78         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
79      ENDIF
80
81      ! 0. Local constant initialization
82      ! --------------------------------
83      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
84
85      ! 1. Apply semi-implicit bottom friction
86      ! --------------------------------------
87      ! Only needed for semi-implicit bottom friction setup. The explicit
88      ! bottom friction has been included in "u(v)a" which act as the R.H.S
89      ! column vector of the tri-diagonal matrix equation
90      !
91      IF( ln_bfrimp ) THEN
92      DO jk = 1, jpk
93         DO jj = 1, jpj
94            DO ji =1, jpi
95               zavmu(ji,jj,jk) = avmu(ji,jj,jk)
96               zavmv(ji,jj,jk) = avmv(ji,jj,jk)
97            END DO
98         END DO
99       END DO
100
101# if defined key_vectopt_loop
102      DO jj = 1, 1
103         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
104# else
105      DO jj = 2, jpjm1
106         DO ji = 2, jpim1
107# endif
108            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
109            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
110            zavmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
111            zavmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
112         END DO
113      END DO
114      ENDIF
115
116      ! 2. Vertical diffusion on u
117      ! ---------------------------
118      ! Matrix and second member construction
119      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
120      ! non zero value at the ocean bottom depending on the bottom friction used.
121      !
122      DO jk = 1, jpkm1        ! Matrix
123         DO jj = 2, jpjm1 
124            DO ji = fs_2, fs_jpim1   ! vector opt.
125               zcoef = - p2dt / fse3u(ji,jj,jk)
126               zzwi          = zcoef * zavmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
127               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
128               zzws          = zcoef * zavmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
129               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
130               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
131            END DO
132         END DO
133      END DO
134      DO jj = 2, jpjm1        ! Surface boudary conditions
135         DO ji = fs_2, fs_jpim1   ! vector opt.
136            zwi(ji,jj,1) = 0._wp
137            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
138         END DO
139      END DO
140
141      ! Matrix inversion starting from the first level
142      !-----------------------------------------------------------------------
143      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
144      !
145      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
146      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
147      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
148      !        (        ...               )( ...  ) ( ...  )
149      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
150      !
151      !   m is decomposed in the product of an upper and a lower triangular matrix
152      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
153      !   The solution (the after velocity) is in ua
154      !-----------------------------------------------------------------------
155      !
156      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
157         DO jj = 2, jpjm1   
158            DO ji = fs_2, fs_jpim1   ! vector opt.
159               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
160            END DO
161         END DO
162      END DO
163      !
164      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
165         DO ji = fs_2, fs_jpim1   ! vector opt.
166            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
167               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
168         END DO
169      END DO
170      DO jk = 2, jpkm1
171         DO jj = 2, jpjm1   
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
174               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
175            END DO
176         END DO
177      END DO
178      !
179      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
180         DO ji = fs_2, fs_jpim1   ! vector opt.
181            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
182         END DO
183      END DO
184      DO jk = jpk-2, 1, -1
185         DO jj = 2, jpjm1   
186            DO ji = fs_2, fs_jpim1   ! vector opt.
187               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
188            END DO
189         END DO
190      END DO
191
192      ! Normalization to obtain the general momentum trend ua
193      DO jk = 1, jpkm1
194         DO jj = 2, jpjm1   
195            DO ji = fs_2, fs_jpim1   ! vector opt.
196               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
197            END DO
198         END DO
199      END DO
200
201
202      ! 3. Vertical diffusion on v
203      ! ---------------------------
204      ! Matrix and second member construction
205      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
206      ! non zero value at the ocean bottom depending on the bottom friction used
207      !
208      DO jk = 1, jpkm1        ! Matrix
209         DO jj = 2, jpjm1   
210            DO ji = fs_2, fs_jpim1   ! vector opt.
211               zcoef = -p2dt / fse3v(ji,jj,jk)
212               zzwi          = zcoef * zavmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
213               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
214               zzws          = zcoef * zavmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
215               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
216               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
217            END DO
218         END DO
219      END DO
220      DO jj = 2, jpjm1        ! Surface boudary conditions
221         DO ji = fs_2, fs_jpim1   ! vector opt.
222            zwi(ji,jj,1) = 0._wp
223            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
224         END DO
225      END DO
226
227      ! Matrix inversion
228      !-----------------------------------------------------------------------
229      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
230      !
231      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
232      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
233      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
234      !        (        ...               )( ...  ) ( ...  )
235      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
236      !
237      !   m is decomposed in the product of an upper and lower triangular matrix
238      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
239      !   The solution (after velocity) is in 2d array va
240      !-----------------------------------------------------------------------
241      !
242      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
243         DO jj = 2, jpjm1   
244            DO ji = fs_2, fs_jpim1   ! vector opt.
245               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
246            END DO
247         END DO
248      END DO
249      !
250      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
251         DO ji = fs_2, fs_jpim1   ! vector opt.
252            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
253               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
254         END DO
255      END DO
256      DO jk = 2, jpkm1
257         DO jj = 2, jpjm1
258            DO ji = fs_2, fs_jpim1   ! vector opt.
259               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
260               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
261            END DO
262         END DO
263      END DO
264      !
265      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
266         DO ji = fs_2, fs_jpim1   ! vector opt.
267            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
268         END DO
269      END DO
270      DO jk = jpk-2, 1, -1
271         DO jj = 2, jpjm1   
272            DO ji = fs_2, fs_jpim1   ! vector opt.
273               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
274            END DO
275         END DO
276      END DO
277
278      ! Normalization to obtain the general momentum trend va
279      DO jk = 1, jpkm1
280         DO jj = 2, jpjm1   
281            DO ji = fs_2, fs_jpim1   ! vector opt.
282               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
283            END DO
284         END DO
285      END DO
286      !
287      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws ) 
288      !
289      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
290      !
291   END SUBROUTINE dyn_zdf_imp
292
293   !!==============================================================================
294END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.