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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 9286

Last change on this file since 9286 was 9286, checked in by cetlod, 6 years ago

bugfix on ocean kinetic energy dissipation due to vertical friction diag., see ticket #2005
  • Property svn:keywords set to Id
File size: 19.6 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 dynadv          ! dynamics: vector invariant versus flux form
23   USE dynspg_oce, ONLY: lk_dynspg_ts
24   USE zdfbfr          ! Bottom friction setup
25   !
26   USE in_out_manager  ! I/O manager
27   USE iom             ! I/O library
28   USE lib_mpp         ! MPP library
29   USE wrk_nemo        ! Memory Allocation
30   USE timing          ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_zdf_imp   ! called by step.F90
36
37   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE dyn_zdf_imp( kt, p2dt )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE dyn_zdf_imp  ***
52      !!                   
53      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
54      !!      and the surface forcing, and add it to the general trend of
55      !!      the momentum equations.
56      !!
57      !! ** Method  :   The vertical momentum mixing trend is given by :
58      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
59      !!      backward time stepping
60      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
61      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
62      !!      Add this trend to the general trend ua :
63      !!         ua = ua + dz( avmu dz(u) )
64      !!
65      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
66      !!---------------------------------------------------------------------
67      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
68      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
69      !!
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, zzz
74      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d
75      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
76      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d
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      !
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      ! 3. Vertical diffusion on v
264      ! ---------------------------
265      ! Matrix and second member construction
266      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
267      ! non zero value at the ocean bottom depending on the bottom friction used
268      !
269      DO jk = 1, jpkm1        ! Matrix
270         DO jj = 2, jpjm1   
271            DO ji = fs_2, fs_jpim1   ! vector opt.
272               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
273               zcoef = - p2dt / ze3va
274               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
275               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
276               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
277               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
278               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
279            END DO
280         END DO
281      END DO
282      DO jj = 2, jpjm1        ! Surface boundary conditions
283         DO ji = fs_2, fs_jpim1   ! vector opt.
284            zwi(ji,jj,1) = 0._wp
285            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
286         END DO
287      END DO
288
289      ! Matrix inversion
290      !-----------------------------------------------------------------------
291      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
292      !
293      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
294      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
295      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
296      !        (        ...               )( ...  ) ( ...  )
297      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
298      !
299      !   m is decomposed in the product of an upper and lower triangular matrix
300      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
301      !   The solution (after velocity) is in 2d array va
302      !-----------------------------------------------------------------------
303      !
304      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
305      DO jk = 2, jpkm1       
306         DO jj = 2, jpjm1   
307            DO ji = fs_2, fs_jpim1   ! vector opt.
308               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
309            END DO
310         END DO
311      END DO
312      !
313      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
314         DO ji = fs_2, fs_jpim1   ! vector opt.
315#if defined key_dynspg_ts           
316            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
317            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
318               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)
319#else
320            va(ji,jj,1) = vb(ji,jj,1) &
321               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
322               &                                      / ( fse3v(ji,jj,1) * rau0     ) * vmask(ji,jj,1) )
323#endif
324         END DO
325      END DO
326      DO jk = 2, jpkm1
327         DO jj = 2, jpjm1
328            DO ji = fs_2, fs_jpim1   ! vector opt.
329#if defined key_dynspg_ts
330               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
331#else
332               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
333#endif
334               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
335            END DO
336         END DO
337      END DO
338      !
339      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
340         DO ji = fs_2, fs_jpim1   ! vector opt.
341            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
342         END DO
343      END DO
344      DO jk = jpk-2, 1, -1
345         DO jj = 2, jpjm1
346            DO ji = fs_2, fs_jpim1
347               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
348            END DO
349         END DO
350      END DO
351
352      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area
353         !                                ! due to v friction (v=vertical)
354         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points)
355         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of
356         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured).
357         CALL wrk_alloc(jpi, jpj,      z2d )
358         CALL wrk_alloc(jpi, jpj, jpk, z3d )
359         z2d(:,:) = 0._wp
360         z3d(:,:,:) = ua(:,:,:)     ;      CALL lbc_lnk( z3d,'U', -1. )
361         DO jk = 2, jpkm1
362            DO jj = 2, jpjm1
363               DO ji = fs_2, fs_jpim1   ! vector opt.
364                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
365                     &   avmu(ji  ,jj,jk) * ( z3d(ji  ,jj,jk-1) - z3d(ji  ,jj,jk) )**2 / fse3uw(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
366                     & + avmu(ji-1,jj,jk) * ( z3d(ji-1,jj,jk-1) - z3d(ji-1,jj,jk) )**2 / fse3uw(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
367                     &                        )
368               END DO
369            END DO
370         END DO
371         z3d(:,:,:) = va(:,:,:)     ;      CALL lbc_lnk( z3d,'V', -1. )
372         DO jk = 2, jpkm1
373            DO jj = 2, jpjm1
374               DO ji = fs_2, fs_jpim1   ! vector opt.
375                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
376                     &   avmv(ji,jj  ,jk) * ( z3d(ji,jj  ,jk-1) - z3d(ji,jj  ,jk) )**2 / fse3vw(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
377                     & + avmv(ji,jj-1,jk) * ( z3d(ji,jj-1,jk-1) - z3d(ji,jj-1,jk) )**2 / fse3vw(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
378                     &                        )
379               END DO
380            END DO
381         END DO
382         zzz= - 0.5_wp* rau0           ! caution sign minus here
383         z2d(:,:) = zzz * z2d(:,:) 
384         CALL lbc_lnk( z2d,'T', 1. )
385         CALL iom_put( 'dispkevfo', z2d )
386         !
387         CALL wrk_dealloc(jpi, jpj,     z2d  )
388         CALL wrk_dealloc(jpi, jpj, jpk, z3d )
389      ENDIF
390
391#if ! defined key_dynspg_ts
392!!gm this can be removed if tranxt is changed like in the trunk so that implicit outcome with
393!!gm the after velocity, not a trend
394      ! Normalization to obtain the general momentum trend ua
395      DO jk = 1, jpkm1
396         DO jj = 2, jpjm1   
397            DO ji = fs_2, fs_jpim1   ! vector opt.
398               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
399               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
400            END DO
401         END DO
402      END DO
403#endif
404
405      ! J. Chanut: Lines below are useless ?
406      !! restore bottom layer avmu(v)
407      IF( ln_bfrimp ) THEN
408        DO jj = 2, jpjm1
409           DO ji = 2, jpim1
410              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
411              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
412              avmu(ji,jj,ikbu+1) = 0.e0
413              avmv(ji,jj,ikbv+1) = 0.e0
414           END DO
415        END DO
416        IF (ln_isfcav) THEN
417           DO jj = 2, jpjm1
418              DO ji = 2, jpim1
419                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
420                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
421                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
422                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
423              END DO
424           END DO
425        END IF
426      ENDIF
427      !
428      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
429      !
430      !
431      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
432      !
433   END SUBROUTINE dyn_zdf_imp
434
435   !!==============================================================================
436END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.