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

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

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