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

Last change on this file since 8627 was 8627, checked in by gm, 3 years ago

#1963 : Correct implementation of CMIP6 KE dissipation + move EKE⇔PE diag in traadv_eiv (in both cases, no more need of ln_trd_KE=T)

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