source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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