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

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 8168

Last change on this file since 8168 was 8168, checked in by glong, 7 years ago

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