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.
domvvl.F90 in branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 2917

Last change on this file since 2917 was 2917, checked in by mlelod, 13 years ago

save memory and cpu in the layer case, see ticket/863?
  • Property svn:keywords set to Id
File size: 40.6 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10#if defined key_vvl
11   !!----------------------------------------------------------------------
12   !!   'key_vvl'                              variable volume
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
16   !!   dom_vvl_nxt      : Compute next vertical scale factors
17   !!   dom_vvl_swp      : Swap vertical scale factors and update the vertical grid
18   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
19   !!   dom_vvl_rst      : read/write restart file
20   !!   dom_vvl_ctl      : Check the vvl options
21   !!----------------------------------------------------------------------
22   !! * Modules used
23   USE oce             ! ocean dynamics and tracers
24   USE dom_oce         ! ocean space and time domain
25   USE sbc_oce         ! ocean surface boundary condition
26   USE in_out_manager  ! I/O manager
27   USE iom             ! I/O manager library
28   USE restart         ! ocean restart
29   USE lib_mpp         ! distributed memory computing library
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31
32   IMPLICIT NONE
33   PRIVATE
34
35   !! * Routine accessibility
36   PUBLIC dom_vvl_init       ! called by domain.F90
37   PUBLIC dom_vvl_sf_nxt     ! called by step.F90
38   PUBLIC dom_vvl_sf_swp     ! called by step.F90
39   PUBLIC dom_vvl_interpol   ! called by dynnxt.F90
40
41   !!* Namelist nam_vvl
42   LOGICAL , PUBLIC                                      :: ln_vvl_zstar  = .TRUE.    ! zstar  vertical coordinate
43   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.   ! ztilde vertical coordinate
44   LOGICAL , PUBLIC                                      :: ln_vvl_layer  = .FALSE.   ! level  vertical coordinate
45   LOGICAL , PUBLIC                                      :: ln_vvl_kepe   = .FALSE.   ! kinetic/potential energy transfer
46   !                                                                                  ! conservation: not used yet
47   REAL(wp), PUBLIC                                      :: ahe3          =  0.e0     ! thickness diffusion coefficient
48
49   !! * Module variables
50   INTEGER                                               :: nvvl                      ! choice of vertical coordinate
51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
52   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                   ! low frequency part of hz divergence
53   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors
54   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_e3t           ! retoring period for scale factors
55   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_restore_hdv           ! retoring period for low freq. divergence
56
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66
67CONTAINS
68
69   INTEGER FUNCTION dom_vvl_alloc()
70      !!----------------------------------------------------------------------
71      !!                ***  FUNCTION dom_vvl_alloc  ***
72      !!----------------------------------------------------------------------
73      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
74         ALLOCATE( e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) , e3t_t_b(jpi,jpj,jpk) ,   &
75            &      un_td  (jpi,jpj,jpk) , vn_td  (jpi,jpj,jpk) , STAT = dom_vvl_alloc        )
76         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
77         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
78      ENDIF
79      IF( ln_vvl_ztilde ) THEN
80         ALLOCATE( frq_restore_e3t(jpi,jpj) , frq_restore_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
81         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
82         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
83      ENDIF
84
85   END FUNCTION dom_vvl_alloc
86
87
88   SUBROUTINE dom_vvl_init
89      !!----------------------------------------------------------------------
90      !!                ***  ROUTINE dom_vvl_init  ***
91      !!                   
92      !! ** Purpose :  Initialization of all scale factors, depths
93      !!               and water column heights
94      !!
95      !! ** Method  :  - use restart file and/or initialize
96      !!               - interpolate scale factors
97      !!
98      !! ** Action  : - fse3t_(n/b) and e3t_t_(n/b)
99      !!              - Regrid: fse3(u/v)_n
100      !!                        fse3(u/v)_b       
101      !!                        fse3w_n           
102      !!                        fse3(u/v)w_b     
103      !!                        fse3(u/v)w_n     
104      !!                        fsdept_n, fsdepw_n and fsde3w_n
105      !!                        h(u/v) and h(u/v)r
106      !!              - ht_0 and ht1_0
107      !!
108      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
109      !!----------------------------------------------------------------------
110      !! * Local declarations
111      INTEGER ::   jk
112      !!----------------------------------------------------------------------
113
114      IF(lwp) WRITE(numout,*)
115      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
116      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
117
118#if defined key_zco
119      CALL ctl_stop( 'dom_vvl_init : options key_zco/key_dynspg_rl are incompatible with variable volume option key_vvl')
120#endif
121
122      ! choose vertical coordinate (z_star, z_tilde or layer)
123      ! ==========================
124      CALL dom_vvl_ctl
125
126      ! Allocate module arrays
127      ! ======================
128      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
129
130      ! read or initialize e3t_t_(b/n) and fse3t_(b/n)
131      ! ==============================================
132      CALL dom_vvl_rst( nit000, 'READ' )
133
134      ! Reconstruction of all vertical scale factors at now and before time steps
135      ! =========================================================================
136      ! Horizontal scale factor interpolations
137      ! --------------------------------------
138      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
139      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
140      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
141      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
142      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
143      ! Vertical scale factor interpolations
144      ! ------------------------------------
145      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
146      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
147      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
148      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
149      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
150      ! t- and w- points depth
151      ! ----------------------
152      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
153      fsdepw_n(:,:,1) = 0.e0
154      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
155      DO jk = 2, jpk
156         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
157         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
158         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
159      END DO
160      ! Reference water column height at t-, u- and v- point
161      ! ----------------------------------------------------
162      ht_0(:,:) = 0.e0
163      hu_0(:,:) = 0.e0
164      hv_0(:,:) = 0.e0
165      DO jk = 1, jpk
166         ht_0(:,:) = ht_0(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)
167         hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)
168         hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)
169      END DO
170
171      ! Restoring frequencies for z_tilde coordinate
172      ! ============================================
173      IF( ln_vvl_ztilde ) THEN
174         ! - ML - In the future, this should be tunable by the user
175         ! DO jj = 1, jpj
176         !    DO ji = 1, jpi
177         !       frq_restore_hdv(ji,jj) = 2.e0 * rpi / 86400.e0 / 5.e0   &
178         !          &                     * MAX( SIN( ABS( gphit(ji,jj) ) / .5e0, 1.e0 / 6.e0) )
179         !    END DO
180         ! END DO
181         ! frq_restore_e3t(:,:) = ( 1.e0 / 6.e0 ) * frq_restore_hdv(:,:)
182         frq_restore_e3t(:,:) = 2.e0 * rpi / ( 30.e0 * 86400.e0 )
183         frq_restore_hdv(:,:) = 2.e0 * rpi / (  5.e0 * 86400.e0 )
184         ! frq_restore_hdv(:,:) = 2.e0 * rpi / (  2.e0 * 86400.e0 )
185      ENDIF
186
187   END SUBROUTINE dom_vvl_init
188
189
190   SUBROUTINE dom_vvl_sf_nxt( kt ) 
191      !!----------------------------------------------------------------------
192      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
193      !!                   
194      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
195      !!                 tranxt and dynspg routines
196      !!
197      !! ** Method  :  - z_tilde_case: after scale factor increment computed with
198      !!                               high frequency part of horizontal divergence + retsoring to
199      !!                               towards the background grid + thickness difusion.
200      !!                               Then repartition of ssh INCREMENT proportionnaly
201      !!                               to the "baroclinic" level thickness.
202      !!               - z_star case:  Repartition of ssh proportionnaly to the level thickness.
203      !!
204      !! ** Action  :  - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case
205      !!               - e3t_t_a: after increment of vertical scale factor
206      !!                          in z_tilde case
207      !!               - fse3(t/u/v)_a
208      !!
209      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
210      !!----------------------------------------------------------------------
211      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
212      USE oce     , ONLY:   ze3t    => ta         ! ua used as workspace
213      USE wrk_nemo, ONLY:   zht     => wrk_2d_1   ! 2D REAL workspace
214      USE wrk_nemo, ONLY:   z_scale => wrk_2d_2   ! 2D REAL workspace
215      USE wrk_nemo, ONLY:   zwu     => wrk_2d_3   ! 2D REAL workspace
216      USE wrk_nemo, ONLY:   zwv     => wrk_2d_4   ! 2D REAL workspace
217      USE wrk_nemo, ONLY:   zhdiv   => wrk_2d_5   ! 2D REAL workspace
218      !! * Arguments
219      INTEGER, INTENT( in )                  :: kt                    ! time step
220      !! * Local declarations
221      INTEGER                                :: ji, jj, jk            ! dummy loop indices
222      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
223      REAL(wp)                               :: z2dt                  ! temporary scalars
224      REAL(wp)                               :: z_def_max             ! temporary scalar
225      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
226      LOGICAL                                :: ln_debug = .FALSE.    ! local logical for debug prints
227      !!----------------------------------------------------------------------
228      IF( wrk_in_use(2, 1, 2, 3, 4, 5) ) THEN
229         CALL ctl_stop('dom_vvl_sf_nxt: requested workspace arrays unavailable')   ;   RETURN
230      END IF
231
232      IF(kt == nit000)   THEN
233         IF(lwp) WRITE(numout,*)
234         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
235         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
236      ENDIF
237
238      ! ******************************* !
239      ! After acale factors at t-points !
240      ! ******************************* !
241
242      !                               ! ----------------- !
243      IF( ln_vvl_zstar ) THEN         ! z_star coordinate !
244         !                            ! ----------------- !
245         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
246         DO jk = 1, jpkm1
247            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:)
248         END DO
249
250      !                               ! --------------------------- !
251      ELSE                            ! z_tilde or layer coordinate !
252         !                            ! --------------------------- !
253
254         ! I - Low frequency horizontal divergence
255         ! =======================================
256
257         ! 1 - barotropic divergence
258         ! -------------------------
259         zhdiv(:,:) = 0.
260         zht(:,:)   = 0.
261         DO jk = 1, jpkm1
262            zhdiv(:,:) = zhdiv(:,:) + hdivn(:,:,jk) * fse3t_n(:,:,jk)
263            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
264         END DO
265         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
266
267         ! 2 - restoring equation  (z-tilde case only)
268         ! ----------------------
269         IF( ln_vvl_ztilde ) THEN
270            IF( kt .GT. nit000 ) THEN
271               DO jk = 1, jpkm1
272                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_restore_hdv(:,:)   &
273                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
274               END DO
275            ENDIF
276         END IF
277
278         ! II - after z_tilde increments of vertical scale factors
279         ! =======================================================
280         e3t_t_a(:,:,:) = 0.e0   ! e3t_t_a used to store tendency terms
281
282         ! 1 - High frequency divergence term
283         ! ----------------------------------
284         IF( ln_vvl_ztilde ) THEN   ! z_tilde case
285            DO jk = 1, jpkm1
286               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
287            END DO                  ! layer case
288         ELSE
289            DO jk = 1, jpkm1
290               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
291            END DO
292         END IF
293
294         ! 2 - Restoring term (z-tilde case only)
295         ! ------------------
296         IF( ln_vvl_ztilde ) THEN
297            DO jk = 1, jpk
298               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_restore_e3t(:,:) * e3t_t_b(:,:,jk)
299            END DO
300         END IF
301
302         ! 3 - Thickness diffusion term
303         ! ----------------------------
304         zwu(:,:) = 0.e0
305         zwv(:,:) = 0.e0
306         ! a - first derivative: diffusive fluxes
307         DO jk = 1, jpkm1
308            DO jj = 1, jpjm1
309               DO ji = 1, fs_jpim1   ! vector opt.
310                  un_td(ji,jj,jk) = ahe3 * umask(ji,jj,jk) * e1ur(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj  ,jk) )
311                  vn_td(ji,jj,jk) = ahe3 * vmask(ji,jj,jk) * e2vr(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji  ,jj+1,jk) )
312                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
313                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
314               END DO
315            END DO
316         END DO
317         ! b - correction for last oceanic u-v points
318         DO jj = 1, jpj
319            DO ji = 1, jpi
320               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
321               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
322            END DO
323         END DO
324         ! c - second derivative: divergence of diffusive fluxes
325         DO jk = 1, jpkm1
326            DO jj = 2, jpjm1
327               DO ji = fs_2, fs_jpim1   ! vector opt.
328                  e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
329                     &                                      + vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk) )  &
330                     &                                    * e12t_1(ji,jj)
331               END DO
332            END DO
333         END DO
334         ! d - thickness diffusion equivalent transport: boundary conditions
335         !     (stored for tracer advction and continuity equation)
336         CALL lbc_lnk( un_td , 'U' , -1.)
337         CALL lbc_lnk( vn_td , 'V' , -1.)
338
339         ! 4 - Time stepping of baroclinic scale factors
340         ! ---------------------------------------------
341         ! Leapfrog time stepping
342         ! ~~~~~~~~~~~~~~~~~~~~~~
343         IF( neuler == 0 .AND. kt == nit000 ) THEN
344            z2dt =  rdt
345         ELSE
346            z2dt = 2.e0 * rdt
347         ENDIF
348         CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )
349         e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:)
350         fse3t_a(:,:,:) = fse3t_0(:,:,:) + e3t_t_a(:,:,:)
351
352         ! Maximum deformation control
353         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
354         ! - ML - Should perhaps be put in the namelist
355         z_def_max = 0.9e0
356         ze3t(:,:,jpk) = 0.e0
357         DO jk = 1, jpkm1
358            ze3t(:,:,jk) = e3t_t_a(:,:,jk) / fse3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
359         END DO
360         z_tmax = MAXVAL( ze3t(:,:,:) )
361         z_tmin = MINVAL( ze3t(:,:,:) )
362         ! - ML - test: for the moment, stop simulation for too large e3_t variations
363         IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN
364            ijk_max = MAXLOC( ze3t(:,:,:) )
365            ijk_min = MINLOC( ze3t(:,:,:) )
366            WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmax
367            WRITE(numout, *) 'at i, j, k=', ijk_max
368            WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / fse3t_0(:,:,:) ) =', z_tmin
369            WRITE(numout, *) 'at i, j, k=', ijk_min           
370            CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / fse3t_0(:,:,:) ) too high')
371         ENDIF
372         ! - ML - end test
373         ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used
374         e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:), ( 1.e0 + z_def_max ) * fse3t_0(:,:,:) )
375         e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), ( 1.e0 - z_def_max ) * fse3t_0(:,:,:) )
376
377         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
378         ! ==================================================================================
379         ! add e3t(n-1) "star" Asselin-filtered
380         DO jk = 1, jpkm1
381            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - fse3t_0(:,:,jk) - e3t_t_b(:,:,jk)
382         END DO
383         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
384         ! - ML - baroclinicity error should be better treated in the future
385         !        i.e. locally and not spread over the water column.
386         zht(:,:) = 0.
387         DO jk = 1, jpkm1
388            zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk)
389         END DO
390         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
391         DO jk = 1, jpkm1
392            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
393         END DO
394
395      ENDIF
396
397      IF( ln_debug ) THEN   ! - ML - test: control prints for debuging
398         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
399            WRITE(numout, *) 'kt  =', kt
400            WRITE(numout, *) 'MAXVAL(abs(ht_0-SUM(e3t_t_a))) =',   &
401               &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) - zht(:,:) ) )
402         END IF
403         zht(:,:) = 0.e0
404         DO jk = 1, jpkm1
405            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
406         END DO
407         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =',   &
408            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
409         zht(:,:) = 0.e0
410         DO jk = 1, jpkm1
411            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
412         END DO
413         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_a))) =',   &
414            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
415      END IF
416
417      ! *********************************** !
418      ! After scale factors at u- v- points !
419      ! *********************************** !
420
421      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
422      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
423
424      IF( wrk_not_released(2, 1, 2, 3, 4, 5) ) THEN
425         CALL ctl_stop( 'dom_vvl_sf_nxt: failed to release workspace arrays' )
426      ENDIF
427
428   END SUBROUTINE dom_vvl_sf_nxt
429
430
431   SUBROUTINE dom_vvl_sf_swp( kt )
432      !!----------------------------------------------------------------------
433      !!                ***  ROUTINE dom_vvl_sf_swp  ***
434      !!                   
435      !! ** Purpose :  compute time filter and swap of scale factors
436      !!               compute all depths and related variables for next time step
437      !!               write outputs and restart file
438      !!
439      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
440      !!               - reconstruct scale factor at other grid points (interpolate)
441      !!               - recompute depths and water height fields
442      !!
443      !! ** Action  :  - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step
444      !!               - Recompute:
445      !!                    fse3(u/v)_b       
446      !!                    fse3w_n           
447      !!                    fse3(u/v)w_b     
448      !!                    fse3(u/v)w_n     
449      !!                    fsdept_n, fsdepw_n  and fsde3w_n
450      !!                    h(u/v) and h(u/v)r
451      !!
452      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
453      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
454      !!----------------------------------------------------------------------
455      USE oce, ONLY:   z_e3t_def => ta                ! square of baroclinic scale factors deformation (in %)
456      !! * Arguments
457      INTEGER, INTENT( in )               :: kt       ! time step
458      !! * Local declarations
459      INTEGER                             :: jk       ! dummy loop indices
460      !!----------------------------------------------------------------------
461
462      IF( kt == nit000 )   THEN
463         IF(lwp) WRITE(numout,*)
464         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
465         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
466      ENDIF
467
468      ! Time filter and swap of scale factors
469      ! =====================================
470      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
471      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
472         IF( neuler == 0 .AND. kt == nit000 ) THEN
473            e3t_t_b(:,:,:) = e3t_t_n(:,:,:)
474         ELSE
475            e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * (e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) )
476         ENDIF
477         e3t_t_n(:,:,:) = e3t_t_a(:,:,:)
478      ENDIF
479      fse3t_n(:,:,:) = fse3t_a(:,:,:)
480      fse3u_n(:,:,:) = fse3u_a(:,:,:)
481      fse3v_n(:,:,:) = fse3v_a(:,:,:)
482
483      ! Compute all missing vertical scale factor and depths
484      ! ====================================================
485      ! Horizontal scale factor interpolations
486      ! --------------------------------------
487      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
488      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
489      ! Vertical scale factor interpolations
490      ! ------------------------------------
491      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
492      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
493      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
494      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
495      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
496      ! t- and w- points depth
497      ! ----------------------
498      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
499      fsdepw_n(:,:,1) = 0.e0
500      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
501      DO jk = 2, jpk
502         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
503         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
504         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
505      END DO
506      ! Local depth and Inverse of the local depth of the water column at u- and v- points
507      ! ----------------------------------------------------------------------------------
508      hu(:,:) = 0.
509      hv(:,:) = 0.
510      DO jk = 1, jpk
511         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
512         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
513      END DO
514      ! Inverse of the local depth
515      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
516      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - umask(:,:,1) )
517
518      ! Write outputs
519      ! =============
520      ! - ML - add output variables in xml file for all configurations
521      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - fse3t_0(:,:,:) ) / fse3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
522      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
523      CALL iom_put( "dept"   , fsdept_n (:,:,:) )
524      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
525
526      ! write restart file
527      ! ==================
528      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
529
530   END SUBROUTINE dom_vvl_sf_swp
531
532
533   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
534      !!---------------------------------------------------------------------
535      !!                  ***  ROUTINE dom_vvl__interpol  ***
536      !!
537      !! ** Purpose :   interpolate scale factors from one grid point to another
538      !!
539      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
540      !!                - horizontal interpolation: grid cell surface averaging
541      !!                - vertical interpolation: simple averaging
542      !!----------------------------------------------------------------------
543      !! * Arguments
544      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
545      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
546      CHARACTER(len=1), INTENT( in )                    ::  pout       ! grid point of out scale factors
547      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
548      !! * Local declarations
549      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
550      !!----------------------------------------------------------------------
551      SELECT CASE ( pout )
552         !               ! ------------------------------------- !
553      CASE( 'U' )        ! interpolation from T-point to U-point !
554         !               ! ------------------------------------- !
555         ! horizontal surface weighted interpolation
556         DO jk = 1, jpkm1
557            DO jj = 2, jpjm1
558               DO ji = 1, fs_jpim1   ! vector opt.
559                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12u_1(ji,jj)                                          &
560                     &                  * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - fse3t_0(ji  ,jj,jk) )     &
561                     &                      + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - fse3t_0(ji+1,jj,jk) ) )
562               END DO
563            END DO
564         END DO
565         ! boundary conditions
566         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
567         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3u_0(:,:,:)
568         !               ! ------------------------------------- !
569      CASE( 'V' )        ! interpolation from T-point to V-point !
570         !               ! ------------------------------------- !
571         ! horizontal surface weighted interpolation
572         DO jk = 1, jpkm1
573            DO jj = 1, jpjm1
574               DO ji = fs_2, fs_jpim1   ! vector opt.
575                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * e12v_1(ji,jj)                                          &
576                     &                  * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3t_0(ji,jj  ,jk) )     &
577                     &                      + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3t_0(ji,jj+1,jk) ) )
578               END DO
579            END DO
580         END DO
581         ! boundary conditions
582         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
583         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3v_0(:,:,:)
584         !               ! ------------------------------------- !
585      CASE( 'F' )        ! interpolation from U-point to F-point !
586         !               ! ------------------------------------- !
587         ! horizontal surface weighted interpolation
588         DO jk = 1, jpkm1
589            DO jj = 1, jpjm1
590               DO ji = 1, fs_jpim1   ! vector opt.
591                  pe3_out(ji,jj,jk) = umask(ji,jj,jk) * umask(ji,jj+1,jk) * e12f_1(ji,jj)                      &
592                     &                  * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - fse3u_0(ji,jj  ,jk) )     &
593                     &                      + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - fse3u_0(ji,jj+1,jk) ) )
594               END DO
595            END DO
596         END DO
597         ! boundary conditions
598         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
599         pe3_out(:,:,:) = pe3_out(:,:,:) + fse3f_0(:,:,:)
600         !               ! ------------------------------------- !
601      CASE( 'W' )        ! interpolation from T-point to W-point !
602         !               ! ------------------------------------- !
603         ! vertical simple interpolation
604         pe3_out(:,:,1) = fse3w_0(:,:,1) + pe3_in(:,:,1) - fse3t_0(:,:,1)
605         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
606         DO jk = 2, jpk
607            pe3_out(:,:,jk) = fse3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3t_0(:,:,jk-1) )   &
608               &                              + (      0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3t_0(:,:,jk  ) )
609         END DO
610         !               ! -------------------------------------- !
611      CASE( 'UW' )       ! interpolation from U-point to UW-point !
612         !               ! -------------------------------------- !
613         ! vertical simple interpolation
614         pe3_out(:,:,1) = fse3uw_0(:,:,1) + pe3_in(:,:,1) - fse3u_0(:,:,1)
615         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
616         DO jk = 2, jpk
617            pe3_out(:,:,jk) = fse3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3u_0(:,:,jk-1) )   &
618               &                               + (      0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3u_0(:,:,jk  ) )
619         END DO
620         !               ! -------------------------------------- !
621      CASE( 'VW' )       ! interpolation from V-point to VW-point !
622         !               ! -------------------------------------- !
623         ! vertical simple interpolation
624         pe3_out(:,:,1) = fse3vw_0(:,:,1) + pe3_in(:,:,1) - fse3v_0(:,:,1)
625         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without undirect adressing
626         DO jk = 2, jpk
627            pe3_out(:,:,jk) = fse3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - fse3v_0(:,:,jk-1) )   &
628               &                               + (      0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk  ) - fse3v_0(:,:,jk  ) )
629         END DO
630      END SELECT
631
632   END SUBROUTINE dom_vvl_interpol
633
634
635   SUBROUTINE dom_vvl_rst( kt, cdrw )
636      !!---------------------------------------------------------------------
637      !!                   ***  ROUTINE dom_vvl_rst  ***
638      !!                     
639      !! ** Purpose :   Read or write VVL file in restart file
640      !!
641      !! ** Method  :   use of IOM library
642      !!                if the restart does not contain vertical scale factors,
643      !!                they are set to the _0 values
644      !!                if the restart does not contain vertical scale factors increments (z_tilde),
645      !!                they are set to 0.
646      !!----------------------------------------------------------------------
647      !! * Arguments
648      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
649      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
650      !! * Local declarations
651      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
652      !!----------------------------------------------------------------------
653      !
654      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
655         !                                   ! ===============
656         IF( ln_rstart ) THEN                   !* Read the restart file
657            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
658            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
659            id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. )
660            id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. )
661            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
662            !                             ! --------- !
663            !                             ! all cases !
664            !                             ! --------- !
665            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
666               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
667               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
668               IF( neuler == 0 ) THEN
669                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
670               ENDIF
671            ELSE                                 ! one at least array is missing
672               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
673            ENDIF
674            !                             ! ----------- !
675            IF( ln_vvl_zstar ) THEN       ! z_star case !
676               !                          ! ----------- !
677               IF( MIN( id3, id4 ) > 0 ) THEN
678                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
679               ENDIF
680               !                          ! ----------------------- !
681            ELSE                          ! z_tilde and layer cases !
682               !                          ! ----------------------- !
683               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
684                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) )
685                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) )
686               ELSE                            ! one at least array is missing
687                  e3t_t_b(:,:,:) = 0.e0
688                  e3t_t_n(:,:,:) = 0.e0
689               ENDIF
690               !                          ! ------------ !
691               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
692                  !                       ! ------------ !
693                  IF( id5 > 0 ) THEN  ! required array exists
694                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
695                  ELSE                ! array is missing
696                     hdiv_lf(:,:,:) = 0.e0
697                  ENDIF
698               ENDIF
699            ENDIF
700            !
701         ELSE                                   !* Initialize at "rest"
702            fse3t_b(:,:,:) = fse3t_0(:,:,:)
703            fse3t_n(:,:,:) = fse3t_0(:,:,:)
704            e3t_t_b(:,:,:) = 0.e0
705            e3t_t_n(:,:,:) = 0.e0
706            hdiv_lf(:,:,:) = 0.e0
707         ENDIF
708
709      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
710         !                                   ! ===================
711         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
712         !                                           ! --------- !
713         !                                           ! all cases !
714         !                                           ! --------- !
715         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
716         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
717         !                                           ! ----------------------- !
718         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
719            !                                        ! ----------------------- !
720            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) )
721            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) )
722         END IF
723         !                                           ! -------------!   
724         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
725            !                                        ! ------------ !
726            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
727         ENDIF
728
729      ENDIF
730
731   END SUBROUTINE dom_vvl_rst
732
733
734   SUBROUTINE dom_vvl_ctl
735      !!---------------------------------------------------------------------
736      !!                  ***  ROUTINE dom_vvl_ctl  ***
737      !!               
738      !! ** Purpose :   Control the consistency between namelist options for
739      !!                vertical coordinate and set nvvl
740      !!----------------------------------------------------------------------
741      INTEGER ::   ioptio
742
743      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ahe3! , ln_vvl_kepe
744      !!----------------------------------------------------------------------
745
746      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
747      READ   ( numnam, nam_vvl )
748
749      IF(lwp) THEN                    ! Namelist print
750         WRITE(numout,*)
751         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
752         WRITE(numout,*) '~~~~~~~~~~~'
753         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
754         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
755         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
756         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
757         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
758         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
759         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
760         WRITE(numout,*) '                                         ahe3           = ', ahe3
761      ENDIF
762
763      ioptio = 0                      ! Parameter control
764      IF( ln_vvl_zstar     )   ioptio = ioptio + 1
765      IF( ln_vvl_ztilde    )   ioptio = ioptio + 1
766      IF( ln_vvl_layer     )   ioptio = ioptio + 1
767
768      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
769
770      IF( ln_vvl_zstar  )   nvvl = 1
771      IF( ln_vvl_ztilde )   nvvl = 2
772      IF( ln_vvl_layer  )   nvvl = 3
773
774      IF(lwp) THEN                   ! Print the choice
775         WRITE(numout,*)
776         IF( nvvl ==  1        ) WRITE(numout,*) '              zstar vertical coordinate is used'
777         IF( nvvl ==  2        ) WRITE(numout,*) '              ztilde vertical coordinate is used'
778         IF( nvvl ==  3        ) WRITE(numout,*) '              layer vertical coordinate is used'
779         ! - ML - Option not developed yet
780         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
781         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
782      ENDIF
783
784   END SUBROUTINE dom_vvl_ctl
785
786
787#else
788
789
790   !!----------------------------------------------------------------------
791   !!   Default option :                                      Empty routine
792   !!----------------------------------------------------------------------
793   SUBROUTINE dom_vvl_init
794      WRITE(*,*) 'dom_vvl_init: You should not have seen this print! error?'
795   END SUBROUTINE dom_vvl_init
796   SUBROUTINE dom_vvl_sf_nxt( kt ) 
797      !! * Arguments
798      INTEGER, INTENT( in )  ::    kt
799      WRITE(*,*) 'dom_vvl_sf_nxt: You should not have seen this print! error?', kt
800   END SUBROUTINE dom_vvl_sf_nxt
801   SUBROUTINE dom_vvl_sf_swp( kt ) 
802      !! * Arguments
803      INTEGER, INTENT( in )  ::    kt
804      WRITE(*,*) 'dom_vvl_sf_swp: You should not have seen this print! error?', kt
805   END SUBROUTINE dom_vvl_sf_swp
806   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
807      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in
808      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out
809      CHARACTER(len=1), INTENT( in )                    ::  pout
810      WRITE(*,*) 'dom_vvl_interpol: You should not have seen this print! error?'
811   END SUBROUTINE dom_vvl_interpol
812
813
814#endif
815
816   !!======================================================================
817END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.