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

source: branches/UKMO/dev_r5107_test_branch/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5274

Last change on this file since 5274 was 5274, checked in by davestorkey, 9 years ago

UKMO test branch: Remove keyword updating.

File size: 16.7 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 domvvl          ! variable volume
19   USE sbc_oce         ! surface boundary condition: ocean
20   USE zdf_oce         ! ocean vertical physics
21   USE phycst          ! physical constants
22   USE in_out_manager  ! I/O manager
23   USE lib_mpp         ! MPP library
24   USE zdfbfr          ! Bottom friction setup
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27   USE dynadv          ! dynamics: vector invariant versus flux form
28   USE dynspg_oce, ONLY: lk_dynspg_ts
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   dyn_zdf_imp   ! called by step.F90
34
35   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id: dynzdf_imp.F90 4990 2014-12-15 16:42:49Z timgraham $
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dyn_zdf_imp( kt, p2dt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_zdf_imp  ***
50      !!                   
51      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
52      !!      and the surface forcing, and add it to the general trend of
53      !!      the momentum equations.
54      !!
55      !! ** Method  :   The vertical momentum mixing trend is given by :
56      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
57      !!      backward time stepping
58      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
59      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
60      !!      Add this trend to the general trend ua :
61      !!         ua = ua + dz( avmu dz(u) )
62      !!
63      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
64      !!---------------------------------------------------------------------
65      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
66      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
67      !!
68      INTEGER  ::   ji, jj, jk   ! dummy loop indices
69      INTEGER  ::   ikbu, ikbv   ! local integers
70      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
71      REAL(wp) ::   ze3ua, ze3va
72      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
76      !
77      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
78      !
79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
82         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
83         !
84         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
85         ELSE                ;    r_vvl = 0._wp       
86         ENDIF
87      ENDIF
88
89      ! 0. Local constant initialization
90      ! --------------------------------
91      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
92
93      ! 1. Apply semi-implicit bottom friction
94      ! --------------------------------------
95      ! Only needed for semi-implicit bottom friction setup. The explicit
96      ! bottom friction has been included in "u(v)a" which act as the R.H.S
97      ! column vector of the tri-diagonal matrix equation
98      !
99
100      IF( ln_bfrimp ) THEN
101         DO jj = 2, jpjm1
102            DO ji = 2, jpim1
103               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
104               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
105               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
106               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
107               ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
108               ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
109               IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
110               IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
111            END DO
112         END DO
113      ENDIF
114
115#if defined key_dynspg_ts
116      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
117         DO jk = 1, jpkm1
118            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
119            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
120         END DO
121      ELSE                                            ! applied on thickness weighted velocity
122         DO jk = 1, jpkm1
123            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
124               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
125               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
126            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
127               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
128               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
129         END DO
130      ENDIF
131
132      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
133         ! remove barotropic velocities:
134         DO jk = 1, jpkm1
135            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
136            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
137         END DO
138         ! Add bottom/top stress due to barotropic component only:
139         DO jj = 2, jpjm1       
140            DO ji = fs_2, fs_jpim1   ! vector opt.
141               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
142               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
143               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
144               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
145               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
146               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
147               ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
148               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
149               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
150               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
151               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
152               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
153            END DO
154         END DO
155      ENDIF
156#endif
157
158      ! 2. Vertical diffusion on u
159      ! ---------------------------
160      ! Matrix and second member construction
161      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
162      ! non zero value at the ocean bottom depending on the bottom friction used.
163      !
164      DO jk = 1, jpkm1        ! Matrix
165         DO jj = 2, jpjm1 
166            DO ji = fs_2, fs_jpim1   ! vector opt.
167               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
168               zcoef = - p2dt / ze3ua     
169               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
170               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
171               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
172               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
173               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
174            END DO
175         END DO
176      END DO
177      DO jj = 2, jpjm1        ! Surface boundary conditions
178         DO ji = fs_2, fs_jpim1   ! vector opt.
179            zwi(ji,jj,1) = 0._wp
180            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
181         END DO
182      END DO
183
184      ! Matrix inversion starting from the first level
185      !-----------------------------------------------------------------------
186      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
187      !
188      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
189      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
190      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
191      !        (        ...               )( ...  ) ( ...  )
192      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
193      !
194      !   m is decomposed in the product of an upper and a lower triangular matrix
195      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
196      !   The solution (the after velocity) is in ua
197      !-----------------------------------------------------------------------
198      !
199      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
200      DO jj = 2, jpjm1   
201         DO ji = fs_2, fs_jpim1   ! vector opt.
202            DO jk = miku(ji,jj)+1, jpkm1
203               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
204            END DO
205         END DO
206      END DO
207      !
208      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
209         DO ji = fs_2, fs_jpim1   ! vector opt.
210            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj)) 
211#if defined key_dynspg_ts
212            ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
213               &                                      / ( ze3ua * rau0 ) 
214#else
215            ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) &
216               &                   + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
217               &                                  / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) ) 
218#endif
219            DO jk = miku(ji,jj)+1, jpkm1
220#if defined key_dynspg_ts
221               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
222#else
223               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
224#endif
225               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
226            END DO
227         END DO
228      END DO
229      !
230      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
231         DO ji = fs_2, fs_jpim1   ! vector opt.
232            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
233            DO jk = jpk-2, miku(ji,jj), -1
234               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
235            END DO
236         END DO
237      END DO
238
239#if ! defined key_dynspg_ts
240      ! Normalization to obtain the general momentum trend ua
241      DO jk = 1, jpkm1
242         DO jj = 2, jpjm1   
243            DO ji = fs_2, fs_jpim1   ! vector opt.
244               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
245            END DO
246         END DO
247      END DO
248#endif
249
250      ! 3. Vertical diffusion on v
251      ! ---------------------------
252      ! Matrix and second member construction
253      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
254      ! non zero value at the ocean bottom depending on the bottom friction used
255      !
256      DO jk = 1, jpkm1        ! Matrix
257         DO jj = 2, jpjm1   
258            DO ji = fs_2, fs_jpim1   ! vector opt.
259               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
260               zcoef = - p2dt / ze3va
261               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
262               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
263               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
264               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
265               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
266            END DO
267         END DO
268      END DO
269      DO jj = 2, jpjm1        ! Surface boundary conditions
270         DO ji = fs_2, fs_jpim1   ! vector opt.
271            zwi(ji,jj,1) = 0._wp
272            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
273         END DO
274      END DO
275
276      ! Matrix inversion
277      !-----------------------------------------------------------------------
278      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
279      !
280      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
281      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
282      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
283      !        (        ...               )( ...  ) ( ...  )
284      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
285      !
286      !   m is decomposed in the product of an upper and lower triangular matrix
287      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
288      !   The solution (after velocity) is in 2d array va
289      !-----------------------------------------------------------------------
290      !
291      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
292      DO jj = 2, jpjm1   
293         DO ji = fs_2, fs_jpim1   ! vector opt.
294            DO jk = mikv(ji,jj)+1, jpkm1       
295               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
296            END DO
297         END DO
298      END DO
299      !
300      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
301         DO ji = fs_2, fs_jpim1   ! vector opt.
302            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj)) 
303#if defined key_dynspg_ts           
304            va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
305               &                                      / ( ze3va * rau0 ) 
306#else
307            va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) &
308               &                   + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
309               &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  )
310#endif
311            DO jk = mikv(ji,jj)+1, jpkm1
312#if defined key_dynspg_ts
313               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
314#else
315               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
316#endif
317               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
318            END DO
319         END DO
320      END DO
321      !
322      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
323         DO ji = fs_2, fs_jpim1   ! vector opt.
324            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
325            DO jk = jpk-2, mikv(ji,jj), -1
326               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
327            END DO
328         END DO
329      END DO
330
331      ! Normalization to obtain the general momentum trend va
332#if ! defined key_dynspg_ts
333      DO jk = 1, jpkm1
334         DO jj = 2, jpjm1   
335            DO ji = fs_2, fs_jpim1   ! vector opt.
336               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
337            END DO
338         END DO
339      END DO
340#endif
341
342      ! J. Chanut: Lines below are useless ?
343      !! restore bottom layer avmu(v)
344      IF( ln_bfrimp ) THEN
345        DO jj = 2, jpjm1
346           DO ji = 2, jpim1
347              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
348              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
349              avmu(ji,jj,ikbu+1) = 0.e0
350              avmv(ji,jj,ikbv+1) = 0.e0
351              ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
352              ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
353              IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
354              IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
355           END DO
356        END DO
357      ENDIF
358      !
359      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
362      !
363   END SUBROUTINE dyn_zdf_imp
364
365   !!==============================================================================
366END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.