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

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

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