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

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

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