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 NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM – NEMO

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90 @ 10060

Last change on this file since 10060 was 10060, checked in by jchanut, 6 years ago

ztilde update (1):
Add alternative vertical scale factor horizontal interpolations (dom_vvl_interpol). For simplicity, consider interpolations from T-points thicknesses, whatever the target point is. The default scheme is hardcoded and set to the actual interpolation scheme.
  • Property svn:keywords set to Id
File size: 64.8 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: vvl option includes z_star and z_tilde coordinates
9   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
14   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
15   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
16   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
17   !!   dom_vvl_rst      : read/write restart file
18   !!   dom_vvl_ctl      : Check the vvl options
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE phycst          ! physical constant
22   USE dom_oce         ! ocean space and time domain
23   USE sbc_oce         ! ocean surface boundary condition
24   USE wet_dry         ! wetting and drying
25   USE usrdef_istate   ! user defined initial state (wad only)
26   USE restart         ! ocean restart
27   !
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O manager library
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE timing          ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC  dom_vvl_init       ! called by domain.F90
38   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
39   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
40   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
41
42   !                                                      !!* Namelist nam_vvl
43   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
44   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
45   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
46   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
49   !                                                       ! conservation: not used yet
50   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
51   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
52   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
53   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
54   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
55
56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
57   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence
58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors
60   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
62
63   !! * Substitutions
64#  include "vectopt_loop_substitute.h90"
65   !!----------------------------------------------------------------------
66   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
67   !! $Id$
68   !! Software governed by the CeCILL licence     (./LICENSE)
69   !!----------------------------------------------------------------------
70CONTAINS
71
72   INTEGER FUNCTION dom_vvl_alloc()
73      !!----------------------------------------------------------------------
74      !!                ***  FUNCTION dom_vvl_alloc  ***
75      !!----------------------------------------------------------------------
76      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
77      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
78         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
79            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
80            &      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         un_td = 0._wp
84         vn_td = 0._wp
85      ENDIF
86      IF( ln_vvl_ztilde ) THEN
87         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
88         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
89         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
90      ENDIF
91      !
92   END FUNCTION dom_vvl_alloc
93
94
95   SUBROUTINE dom_vvl_init
96      !!----------------------------------------------------------------------
97      !!                ***  ROUTINE dom_vvl_init  ***
98      !!                   
99      !! ** Purpose :  Initialization of all scale factors, depths
100      !!               and water column heights
101      !!
102      !! ** Method  :  - use restart file and/or initialize
103      !!               - interpolate scale factors
104      !!
105      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
106      !!              - Regrid: e3(u/v)_n
107      !!                        e3(u/v)_b       
108      !!                        e3w_n           
109      !!                        e3(u/v)w_b     
110      !!                        e3(u/v)w_n     
111      !!                        gdept_n, gdepw_n and gde3w_n
112      !!              - h(t/u/v)_0
113      !!              - frq_rst_e3t and frq_rst_hdv
114      !!
115      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
116      !!----------------------------------------------------------------------
117      INTEGER ::   ji, jj, jk
118      INTEGER ::   ii0, ii1, ij0, ij1
119      REAL(wp)::   zcoef
120      !!----------------------------------------------------------------------
121      !
122      IF(lwp) WRITE(numout,*)
123      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
124      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
125      !
126      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
127      !
128      !                    ! Allocate module arrays
129      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
130      !
131      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
132      CALL dom_vvl_rst( nit000, 'READ' )
133      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
134      !
135      !                    !== Set of all other vertical scale factors  ==!  (now and before)
136      !                                ! Horizontal interpolation of e3t
137      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
138      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
139      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V
140      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
141      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
142      !                                ! Vertical interpolation of e3t,u,v
143      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
144      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
145      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
146      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
147      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
148      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
149
150      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
151      e3t_a(:,:,:) = e3t_n(:,:,:)
152      e3u_a(:,:,:) = e3u_n(:,:,:)
153      e3v_a(:,:,:) = e3v_n(:,:,:)
154      !
155      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
156      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
157      gdepw_n(:,:,1) = 0.0_wp
158      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
159      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
160      gdepw_b(:,:,1) = 0.0_wp
161      DO jk = 2, jpk                               ! vertical sum
162         DO jj = 1,jpj
163            DO ji = 1,jpi
164               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
165               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
166               !                             ! 0.5 where jk = mikt     
167!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
168               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
169               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
170               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
171                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
172               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
173               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
174               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
175                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
176            END DO
177         END DO
178      END DO
179      !
180      !                    !==  thickness of the water column  !!   (ocean portion only)
181      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
182      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
183      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
184      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
185      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
186      DO jk = 2, jpkm1
187         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
188         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
189         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
190         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
191         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
192      END DO
193      !
194      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
195      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
196      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
197      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
198      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
199
200      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
201      IF( ln_vvl_ztilde ) THEN
202!!gm : idea: add here a READ in a file of custumized restoring frequency
203         !                                   ! Values in days provided via the namelist
204         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
205         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
206         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
207         !
208         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
209            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
210            frq_rst_hdv(:,:) = 1._wp / rdt
211         ENDIF
212         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
213            DO jj = 1, jpj
214               DO ji = 1, jpi
215!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
216                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
217                     ! values outside the equatorial band and transition zone (ztilde)
218                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
219                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
220                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
221                     ! values inside the equatorial band (ztilde as zstar)
222                     frq_rst_e3t(ji,jj) =  0.0_wp
223                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
224                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
225                     !                                      ! (linearly transition from z-tilde to z-star)
226                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
227                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
228                        &                                          * 180._wp / 3.5_wp ) )
229                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
230                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
231                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
232                        &                                          * 180._wp / 3.5_wp ) )
233                  ENDIF
234               END DO
235            END DO
236            IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
237               ii0 = 103   ;   ii1 = 111       
238               ij0 = 128   ;   ij1 = 135   ;   
239               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
240               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
241            ENDIF
242         ENDIF
243      ENDIF
244      !
245      IF(lwxios) THEN
246! define variables in restart file when writing with XIOS
247         CALL iom_set_rstw_var_active('e3t_b')
248         CALL iom_set_rstw_var_active('e3t_n')
249         !                                           ! ----------------------- !
250         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
251            !                                        ! ----------------------- !
252            CALL iom_set_rstw_var_active('tilde_e3t_b')
253            CALL iom_set_rstw_var_active('tilde_e3t_n')
254         END IF
255         !                                           ! -------------!   
256         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
257            !                                        ! ------------ !
258            CALL iom_set_rstw_var_active('hdiv_lf')
259         ENDIF
260         !
261      ENDIF
262      !
263   END SUBROUTINE dom_vvl_init
264
265
266   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
267      !!----------------------------------------------------------------------
268      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
269      !!                   
270      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
271      !!                 tranxt and dynspg routines
272      !!
273      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
274      !!               - z_tilde_case: after scale factor increment =
275      !!                                    high frequency part of horizontal divergence
276      !!                                  + retsoring towards the background grid
277      !!                                  + thickness difusion
278      !!                               Then repartition of ssh INCREMENT proportionnaly
279      !!                               to the "baroclinic" level thickness.
280      !!
281      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
282      !!               - tilde_e3t_a: after increment of vertical scale factor
283      !!                              in z_tilde case
284      !!               - e3(t/u/v)_a
285      !!
286      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
287      !!----------------------------------------------------------------------
288      INTEGER, INTENT( in )           ::   kt      ! time step
289      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
290      !
291      INTEGER                ::   ji, jj, jk            ! dummy loop indices
292      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
293      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
294      LOGICAL                ::   ll_do_bclinic         ! local logical
295      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
296      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
297      !!----------------------------------------------------------------------
298      !
299      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
300      !
301      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
302      !
303      IF( kt == nit000 ) THEN
304         IF(lwp) WRITE(numout,*)
305         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
306         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
307      ENDIF
308
309      ll_do_bclinic = .TRUE.
310      IF( PRESENT(kcall) ) THEN
311         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
312      ENDIF
313
314      ! ******************************* !
315      ! After acale factors at t-points !
316      ! ******************************* !
317      !                                                ! --------------------------------------------- !
318      !                                                ! z_star coordinate and barotropic z-tilde part !
319      !                                                ! --------------------------------------------- !
320      !
321      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
322      DO jk = 1, jpkm1
323         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
324         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
325      END DO
326      !
327      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
328         !                                                            ! ------baroclinic part------ !
329         ! I - initialization
330         ! ==================
331
332         ! 1 - barotropic divergence
333         ! -------------------------
334         zhdiv(:,:) = 0._wp
335         zht(:,:)   = 0._wp
336         DO jk = 1, jpkm1
337            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
338            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
339         END DO
340         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
341
342         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
343         ! --------------------------------------------------
344         IF( ln_vvl_ztilde ) THEN
345            IF( kt > nit000 ) THEN
346               DO jk = 1, jpkm1
347                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
348                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
349               END DO
350            ENDIF
351         ENDIF
352
353         ! II - after z_tilde increments of vertical scale factors
354         ! =======================================================
355         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
356
357         ! 1 - High frequency divergence term
358         ! ----------------------------------
359         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
360            DO jk = 1, jpkm1
361               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
362            END DO
363         ELSE                         ! layer case
364            DO jk = 1, jpkm1
365               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
366            END DO
367         ENDIF
368
369         ! 2 - Restoring term (z-tilde case only)
370         ! ------------------
371         IF( ln_vvl_ztilde ) THEN
372            DO jk = 1, jpk
373               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
374            END DO
375         ENDIF
376
377         ! 3 - Thickness diffusion term
378         ! ----------------------------
379         zwu(:,:) = 0._wp
380         zwv(:,:) = 0._wp
381         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
382            DO jj = 1, jpjm1
383               DO ji = 1, fs_jpim1   ! vector opt.
384                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
385                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
386                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
387                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
388                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
389                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
390               END DO
391            END DO
392         END DO
393         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
394            DO ji = 1, jpi
395               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
396               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
397            END DO
398         END DO
399         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
400            DO jj = 2, jpjm1
401               DO ji = fs_2, fs_jpim1   ! vector opt.
402                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
403                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
404                     &                                            ) * r1_e1e2t(ji,jj)
405               END DO
406            END DO
407         END DO
408         !                       ! d - thickness diffusion transport: boundary conditions
409         !                             (stored for tracer advction and continuity equation)
410         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
411
412         ! 4 - Time stepping of baroclinic scale factors
413         ! ---------------------------------------------
414         ! Leapfrog time stepping
415         ! ~~~~~~~~~~~~~~~~~~~~~~
416         IF( neuler == 0 .AND. kt == nit000 ) THEN
417            z2dt =  rdt
418         ELSE
419            z2dt = 2.0_wp * rdt
420         ENDIF
421         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
422         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
423
424         ! Maximum deformation control
425         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
426         ze3t(:,:,jpk) = 0._wp
427         DO jk = 1, jpkm1
428            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
429         END DO
430         z_tmax = MAXVAL( ze3t(:,:,:) )
431         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
432         z_tmin = MINVAL( ze3t(:,:,:) )
433         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
434         ! - ML - test: for the moment, stop simulation for too large e3_t variations
435         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
436            IF( lk_mpp ) THEN
437               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
438               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
439            ELSE
440               ijk_max = MAXLOC( ze3t(:,:,:) )
441               ijk_max(1) = ijk_max(1) + nimpp - 1
442               ijk_max(2) = ijk_max(2) + njmpp - 1
443               ijk_min = MINLOC( ze3t(:,:,:) )
444               ijk_min(1) = ijk_min(1) + nimpp - 1
445               ijk_min(2) = ijk_min(2) + njmpp - 1
446            ENDIF
447            IF (lwp) THEN
448               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
449               WRITE(numout, *) 'at i, j, k=', ijk_max
450               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
451               WRITE(numout, *) 'at i, j, k=', ijk_min           
452               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
453            ENDIF
454         ENDIF
455         ! - ML - end test
456         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
457         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
458         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
459
460         !
461         ! "tilda" change in the after scale factor
462         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
463         DO jk = 1, jpkm1
464            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
465         END DO
466         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
467         ! ==================================================================================
468         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
469         ! - ML - baroclinicity error should be better treated in the future
470         !        i.e. locally and not spread over the water column.
471         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
472         zht(:,:) = 0.
473         DO jk = 1, jpkm1
474            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
475         END DO
476         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
477         DO jk = 1, jpkm1
478            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
479         END DO
480
481      ENDIF
482
483      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
484      !                                           ! ---baroclinic part--------- !
485         DO jk = 1, jpkm1
486            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
487         END DO
488      ENDIF
489
490      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
491         !
492         IF( lwp ) WRITE(numout, *) 'kt =', kt
493         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
494            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
495            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
496            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
497         END IF
498         !
499         zht(:,:) = 0.0_wp
500         DO jk = 1, jpkm1
501            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
502         END DO
503         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
504         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
505         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
506         !
507         zht(:,:) = 0.0_wp
508         DO jk = 1, jpkm1
509            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
510         END DO
511         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
512         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
513         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
514         !
515         zht(:,:) = 0.0_wp
516         DO jk = 1, jpkm1
517            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
518         END DO
519         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
520         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
521         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
522         !
523         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
524         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
525         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
526         !
527         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
528         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
529         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
530         !
531         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
532         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
533         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
534      END IF
535
536      ! *********************************** !
537      ! After scale factors at u- v- points !
538      ! *********************************** !
539
540      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
541      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
542
543      ! *********************************** !
544      ! After depths at u- v points         !
545      ! *********************************** !
546
547      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
548      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
549      DO jk = 2, jpkm1
550         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
551         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
552      END DO
553      !                                        ! Inverse of the local depth
554!!gm BUG ?  don't understand the use of umask_i here .....
555      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
556      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
557      !
558      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
559      !
560   END SUBROUTINE dom_vvl_sf_nxt
561
562
563   SUBROUTINE dom_vvl_sf_swp( kt )
564      !!----------------------------------------------------------------------
565      !!                ***  ROUTINE dom_vvl_sf_swp  ***
566      !!                   
567      !! ** Purpose :  compute time filter and swap of scale factors
568      !!               compute all depths and related variables for next time step
569      !!               write outputs and restart file
570      !!
571      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
572      !!               - reconstruct scale factor at other grid points (interpolate)
573      !!               - recompute depths and water height fields
574      !!
575      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
576      !!               - Recompute:
577      !!                    e3(u/v)_b       
578      !!                    e3w_n           
579      !!                    e3(u/v)w_b     
580      !!                    e3(u/v)w_n     
581      !!                    gdept_n, gdepw_n  and gde3w_n
582      !!                    h(u/v) and h(u/v)r
583      !!
584      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
585      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
586      !!----------------------------------------------------------------------
587      INTEGER, INTENT( in ) ::   kt   ! time step
588      !
589      INTEGER  ::   ji, jj, jk   ! dummy loop indices
590      REAL(wp) ::   zcoef        ! local scalar
591      !!----------------------------------------------------------------------
592      !
593      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
594      !
595      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
596      !
597      IF( kt == nit000 )   THEN
598         IF(lwp) WRITE(numout,*)
599         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
600         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
601      ENDIF
602      !
603      ! Time filter and swap of scale factors
604      ! =====================================
605      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
606      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
607         IF( neuler == 0 .AND. kt == nit000 ) THEN
608            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
609         ELSE
610            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
611            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
612         ENDIF
613         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
614      ENDIF
615      gdept_b(:,:,:) = gdept_n(:,:,:)
616      gdepw_b(:,:,:) = gdepw_n(:,:,:)
617
618      e3t_n(:,:,:) = e3t_a(:,:,:)
619      e3u_n(:,:,:) = e3u_a(:,:,:)
620      e3v_n(:,:,:) = e3v_a(:,:,:)
621
622      ! Compute all missing vertical scale factor and depths
623      ! ====================================================
624      ! Horizontal scale factor interpolations
625      ! --------------------------------------
626      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
627      ! - JC - hu_b, hv_b, hur_b, hvr_b also
628     
629      CALL dom_vvl_interpol( e3t_n(:,:,:), e3f_n(:,:,:), 'F'  )
630     
631      ! Vertical scale factor interpolations
632      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
633      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
634      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
635      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
636      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
637      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
638
639      ! t- and w- points depth (set the isf depth as it is in the initial step)
640      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
641      gdepw_n(:,:,1) = 0.0_wp
642      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
643      DO jk = 2, jpk
644         DO jj = 1,jpj
645            DO ji = 1,jpi
646              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
647                                                                 ! 1 for jk = mikt
648               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
649               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
650               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
651                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
652               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
653            END DO
654         END DO
655      END DO
656
657      ! Local depth and Inverse of the local depth of the water
658      ! -------------------------------------------------------
659      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
660      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
661      !
662      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
663      DO jk = 2, jpkm1
664         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
665      END DO
666
667      ! write restart file
668      ! ==================
669      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
670      !
671      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
672      !
673   END SUBROUTINE dom_vvl_sf_swp
674
675
676   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
677      !!---------------------------------------------------------------------
678      !!                  ***  ROUTINE dom_vvl__interpol  ***
679      !!
680      !! ** Purpose :   interpolate scale factors from one grid point to another
681      !!
682      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
683      !!                - horizontal interpolation: grid cell surface averaging
684      !!                - vertical interpolation: simple averaging
685      !!----------------------------------------------------------------------
686      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
687      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
688      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
689      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
690      !
691      INTEGER ::   ji, jj, jk, jkbot                                ! dummy loop indices
692      INTEGER ::   nmet                                             ! horizontal interpolation method
693      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
694      REAL(wp) ::  ztap, zsmall                                     ! Parameters defining minimum thicknesses UVF-points
695      REAL(wp) ::  zmin
696      REAL(wp) ::  zdo, zup                                         ! Lower and upper interfaces depths anomalies
697      REAL(wp), DIMENSION(jpi,jpj) :: zs                            ! Surface interface depth anomaly
698      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw                        ! Interface depth anomaly
699      !!----------------------------------------------------------------------
700      !
701      nmet = 0  ! Original method (Surely wrong)
702!     nmet = 1  ! Interface interpolation
703!     nmet = 2  ! Internal interfaces interpolation only, spread barotropic increment
704!     Note that we kept surface weighted interpolation for barotropic increment to be compliant
705!     with what is done in surface pressure module.
706      !
707      IF(ln_wd_il) THEN
708        zlnwd = 1.0_wp
709      ELSE
710        zlnwd = 0.0_wp
711      END IF
712      !
713      ztap   = 0.0_wp   ! Minimum fraction of T-point thickness at cell interfaces
714      zsmall = 1.e-3_wp ! Minimum thickness at U or V points (m)
715      !
716      IF ( (nmet==1).OR.(nmet==2) ) THEN
717         SELECT CASE ( pout )
718            !
719         CASE( 'U', 'V', 'F' )
720            ! Compute interface depth anomaly at T-points
721            !
722            zw(:,:,:) =  0._wp   
723            !
724            DO jk=2,jpk
725               zw(:,:,jk) = zw(:,:,jk-1) + pe3_in(:,:,jk-1)*tmask(:,:,jk-1)
726            END DO 
727            ! Interface depth anomalies:
728            DO jk=1,jpkm1
729               zw(:,:,jk) = zw(:,:,jk) - zw(:,:,jpk) + ht_0(:,:)
730            END DO
731            zw(:,:,jpk) = ht_0(:,:)
732            !
733            IF (nmet==2) THEN        ! Consider "internal" interfaces only
734               zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh)
735               !
736               DO jj = 1, jpj
737                  DO ji = 1, jpi
738                     DO jk=1,jpk
739                        zw(ji,jj,jk) = (zw(ji,jj,jk) + zs(ji,jj))                                          &
740                                     & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1))  &
741                                     & * tmask(ji,jj,jk)
742                     END DO
743                  END DO
744               END DO
745            ENDIF 
746            !
747         END SELECT
748      END IF
749      !
750      pe3_out(:,:,:) = 0.0_wp
751      !
752      SELECT CASE ( pout )    !==  type of interpolation  ==!
753         !
754      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
755         IF (nmet==0) THEN
756            DO jk = 1, jpk
757               DO jj = 1, jpjm1
758                  DO ji = 1, fs_jpim1   ! vector opt.
759                     pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
760                        &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
761                        &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
762                  END DO
763               END DO
764            END DO
765         ELSE
766            DO jj = 1, jpjm1
767               DO ji = 1, fs_jpim1   ! vector opt.
768                  ! Correction at last level:
769                  jkbot = mbku(ji,jj)
770                  zdo   = hu_0(ji,jj)
771                  DO jk=jkbot,1,-1
772                     zup = 0.5_wp * ( zw(ji  ,jj,jk) + zw(ji+1,jj,jk) ) 
773                     !
774                     ! If there is a step, taper bottom interface:
775                     IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(zup>zdo)) THEN
776                        IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN
777                           zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))
778                        ELSE
779                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
780                        ENDIF
781                        zup = MIN(zup, zdo-zmin)
782                     ENDIF
783                     zup = MIN(zup, zdo-zsmall)
784                     pe3_out(ji,jj,jk) = (zdo - zup - e3u_0(ji,jj,jk)) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
785                     zdo = zup
786                  END DO
787               END DO
788            END DO
789         END IF
790         !
791         IF (nmet==2) THEN           ! Spread sea level anomaly
792            DO jj = 1, jpjm1
793               DO ji = 1, fs_jpim1   ! vector opt.
794                  DO jk=1,jpk
795                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                   &
796                           &               + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) )               & 
797                           &               / ( hu_0(ji,jj) + 1._wp - umask(ji,jj,1) )              &
798                           &               * 0.5_wp * r1_e1e2u(ji,jj)                              &
799                           &               * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )        &
800                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji+1,jj)*zs(ji+1,jj) )       
801                  END DO
802               END DO
803            END DO
804            !
805         ENDIF
806         !
807         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
808         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
809         !
810      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
811         IF (nmet==0) THEN
812            DO jk = 1, jpk
813               DO jj = 1, jpjm1
814                  DO ji = 1, fs_jpim1   ! vector opt.
815                     pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
816                        &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
817                        &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
818                  END DO
819               END DO
820            END DO
821         ELSE
822            DO jj = 1, jpjm1
823               DO ji = 1, fs_jpim1   ! vector opt.
824                  ! Correction at last level:
825                  jkbot = mbkv(ji,jj)
826                  zdo   = hv_0(ji,jj)
827                  DO jk=jkbot,1,-1
828                     zup = 0.5_wp * ( zw(ji,jj  ,jk) + zw(ji,jj+1,jk) ) 
829                     !
830                     ! If there is a step, taper bottom interface:
831                     IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN
832                        IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN
833                           zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))
834                        ELSE
835                           zmin = ztap * (zw(ji  ,jj,jk+1)-zw(ji  ,jj,jk))
836                        ENDIF
837                        zup = MIN(zup, zdo-zmin)
838                     ENDIF
839                     zup = MIN(zup, zdo-zsmall)
840                     pe3_out(ji,jj,jk) = (zdo - zup - e3v_0(ji,jj,jk)) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )
841                     zdo = zup
842                  END DO
843               END DO
844            END DO
845         END IF
846         !
847         IF (nmet==2) THEN           ! Spread sea level anomaly
848            DO jj = 1, jpjm1
849               DO ji = 1, fs_jpim1   ! vector opt.
850                  DO jk=1,jpk
851                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                                       &
852                           &               + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) )                                   & 
853                           &               / ( hv_0(ji,jj) + 1._wp - vmask(ji,jj,1) )                                  &
854                           &               * 0.5_wp * r1_e1e2v(ji,jj)                                                  &
855                                           * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )                            &
856                           &               * ( e1e2t(ji,jj)*zs(ji,jj) + e1e2t(ji,jj+1)*zs(ji,jj+1) )       
857                  END DO
858               END DO
859            END DO
860            !
861         ENDIF
862         !
863         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
864         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
865         !
866      CASE( 'F' )                   !* from T-point to F-point : hor. surface weighted mean
867         IF (nmet==0) THEN
868            DO jk=1,jpk
869               DO jj = 1, jpjm1
870                  DO ji = 1, fs_jpim1   ! vector opt.
871                     pe3_out(ji,jj,jk) = 0.25_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
872                           &                     *    r1_e1e2f(ji,jj)                                                  &
873                           &                     * (  e1e2t(ji  ,jj  ) * ( pe3_in(ji  ,jj  ,jk)-e3t_0(ji  ,jj  ,jk) )  & 
874                           &                        + e1e2t(ji  ,jj+1) * ( pe3_in(ji  ,jj+1,jk)-e3t_0(ji  ,jj+1,jk) )  &
875                           &                        + e1e2t(ji+1,jj  ) * ( pe3_in(ji+1,jj  ,jk)-e3t_0(ji+1,jj  ,jk) )  & 
876                           &                        + e1e2t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) ) )                                                 
877                  END DO
878               END DO
879            END DO
880         ELSE
881            DO jj = 1, jpjm1
882               DO ji = 1, fs_jpim1   ! vector opt.
883                  ! bottom correction:
884                  jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 
885                  zdo = hf_0(ji,jj)
886                  DO jk=jkbot,1,-1
887                     zup =  0.25_wp * (  zw(ji  ,jj  ,jk)  & 
888                           &           + zw(ji+1,jj  ,jk)  & 
889                           &           + zw(ji  ,jj+1,jk)  & 
890                           &           + zw(ji+1,jj+1,jk)  )
891                     !
892                     ! If there is a step, taper bottom interface:
893                     IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj  ) + hu_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN
894                        IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN
895                           IF ( ht_0(ji+1,jj+1) < ht_0(ji  ,jj+1) ) THEN
896                              zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))
897                           ELSE
898                              zmin = ztap * (zw(ji  ,jj+1,jk+1)-zw(ji  ,jj+1,jk))
899                           ENDIF
900                        ELSE
901                           IF ( ht_0(ji+1,jj  ) < ht_0(ji  ,jj  ) ) THEN
902                              zmin = ztap * (zw(ji+1,jj  ,jk+1)-zw(ji+1,jj  ,jk))
903                           ELSE
904                              zmin = ztap * (zw(ji  ,jj  ,jk+1)-zw(ji  ,jj  ,jk))
905                           ENDIF
906                        ENDIF
907                        zup = MIN(zup, zdo-zmin)
908                     ENDIF
909                     zup = MIN(zup, zdo-zsmall)
910                     !
911                     pe3_out(ji,jj,jk) = ( zdo - zup - e3f_0(ji,jj,jk) ) & 
912                                      & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd )
913                     zdo = zup
914                  END DO
915               END DO
916            END DO
917         END IF
918         !
919         IF (nmet==2) THEN           ! Spread sea level anomaly
920            !
921            DO jj = 1, jpjm1
922               DO ji = 1, fs_jpim1   ! vector opt.
923                  DO jk=1,jpk
924                     pe3_out(ji,jj,jk) =       pe3_out(ji,jj,jk)                                           &
925                           &               + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) )                       & 
926                           &               / ( hf_0(ji,jj) + 1._wp - umask(ji,jj,1)*umask(ji,jj+1,1) )     &
927                           &               * 0.25_wp * r1_e1e2f(ji,jj)                                     & 
928                           &               * ( umask(ji,jj,jk)*umask(ji,jj+1,jk)*(1.0_wp - zlnwd) + zlnwd )&
929                           &               * ( e1e2t(ji  ,jj)*zs(ji  ,jj) + e1e2t(ji  ,jj+1)*zs(ji  ,jj+1) &
930                           &                  +e1e2t(ji+1,jj)*zs(ji+1,jj) + e1e2t(ji+1,jj+1)*zs(ji+1,jj+1) )               
931                  END DO
932               END DO
933            END DO
934         END IF
935         !
936         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
937         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
938         !
939      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
940         !
941         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
942         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
943!!gm BUG? use here wmask in case of ISF ?  to be checked
944         DO jk = 2, jpk
945            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
946               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
947               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
948               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
949         END DO
950         !
951      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
952         !
953         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
954         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
955!!gm BUG? use here wumask in case of ISF ?  to be checked
956         DO jk = 2, jpk
957            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
958               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
959               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
960               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
961         END DO
962         !
963      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
964         !
965         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
966         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
967!!gm BUG? use here wvmask in case of ISF ?  to be checked
968         DO jk = 2, jpk
969            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
970               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
971               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
972               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
973         END DO
974      END SELECT
975      !
976   END SUBROUTINE dom_vvl_interpol
977
978
979   SUBROUTINE dom_vvl_rst( kt, cdrw )
980      !!---------------------------------------------------------------------
981      !!                   ***  ROUTINE dom_vvl_rst  ***
982      !!                     
983      !! ** Purpose :   Read or write VVL file in restart file
984      !!
985      !! ** Method  :   use of IOM library
986      !!                if the restart does not contain vertical scale factors,
987      !!                they are set to the _0 values
988      !!                if the restart does not contain vertical scale factors increments (z_tilde),
989      !!                they are set to 0.
990      !!----------------------------------------------------------------------
991      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
992      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
993      !
994      INTEGER ::   ji, jj, jk
995      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
996      !!----------------------------------------------------------------------
997      !
998      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
999         !                                   ! ===============
1000         IF( ln_rstart ) THEN                   !* Read the restart file
1001            CALL rst_read_open                  !  open the restart file if necessary
1002            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
1003            !
1004            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
1005            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
1006            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1007            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1008            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
1009            !                             ! --------- !
1010            !                             ! all cases !
1011            !                             ! --------- !
1012            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
1013               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1014               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1015               ! needed to restart if land processor not computed
1016               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
1017               WHERE ( tmask(:,:,:) == 0.0_wp ) 
1018                  e3t_n(:,:,:) = e3t_0(:,:,:)
1019                  e3t_b(:,:,:) = e3t_0(:,:,:)
1020               END WHERE
1021               IF( neuler == 0 ) THEN
1022                  e3t_b(:,:,:) = e3t_n(:,:,:)
1023               ENDIF
1024            ELSE IF( id1 > 0 ) THEN
1025               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
1026               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
1027               IF(lwp) write(numout,*) 'neuler is forced to 0'
1028               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
1029               e3t_n(:,:,:) = e3t_b(:,:,:)
1030               neuler = 0
1031            ELSE IF( id2 > 0 ) THEN
1032               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
1033               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
1034               IF(lwp) write(numout,*) 'neuler is forced to 0'
1035               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
1036               e3t_b(:,:,:) = e3t_n(:,:,:)
1037               neuler = 0
1038            ELSE
1039               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
1040               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1041               IF(lwp) write(numout,*) 'neuler is forced to 0'
1042               DO jk = 1, jpk
1043                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
1044                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1045                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
1046               END DO
1047               e3t_b(:,:,:) = e3t_n(:,:,:)
1048               neuler = 0
1049            ENDIF
1050            !                             ! ----------- !
1051            IF( ln_vvl_zstar ) THEN       ! z_star case !
1052               !                          ! ----------- !
1053               IF( MIN( id3, id4 ) > 0 ) THEN
1054                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1055               ENDIF
1056               !                          ! ----------------------- !
1057            ELSE                          ! z_tilde and layer cases !
1058               !                          ! ----------------------- !
1059               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
1060                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
1061                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
1062               ELSE                            ! one at least array is missing
1063                  tilde_e3t_b(:,:,:) = 0.0_wp
1064                  tilde_e3t_n(:,:,:) = 0.0_wp
1065               ENDIF
1066               !                          ! ------------ !
1067               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
1068                  !                       ! ------------ !
1069                  IF( id5 > 0 ) THEN  ! required array exists
1070                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
1071                  ELSE                ! array is missing
1072                     hdiv_lf(:,:,:) = 0.0_wp
1073                  ENDIF
1074               ENDIF
1075            ENDIF
1076            !
1077         ELSE                                   !* Initialize at "rest"
1078            !
1079
1080            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
1081               !
1082               IF( cn_cfg == 'wad' ) THEN
1083                  ! Wetting and drying test case
1084                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
1085                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
1086                  sshn (:,:)     = sshb(:,:)
1087                  un   (:,:,:)   = ub  (:,:,:)
1088                  vn   (:,:,:)   = vb  (:,:,:)
1089               ELSE
1090                  ! if not test case
1091                  sshn(:,:) = -ssh_ref
1092                  sshb(:,:) = -ssh_ref
1093
1094                  DO jj = 1, jpj
1095                     DO ji = 1, jpi
1096                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
1097
1098                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1099                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1100                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
1101                        ENDIF
1102                     ENDDO
1103                  ENDDO
1104               ENDIF !If test case else
1105
1106               ! Adjust vertical metrics for all wad
1107               DO jk = 1, jpk
1108                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
1109                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1110                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1111               END DO
1112               e3t_b(:,:,:) = e3t_n(:,:,:)
1113
1114               DO ji = 1, jpi
1115                  DO jj = 1, jpj
1116                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
1117                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
1118                     ENDIF
1119                  END DO
1120               END DO 
1121               !
1122            ELSE
1123               !
1124               ! Just to read set ssh in fact, called latter once vertical grid
1125               ! is set up:
1126!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
1127!               !
1128!               DO jk=1,jpk
1129!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
1130!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
1131!               END DO
1132!               e3t_n(:,:,:) = e3t_b(:,:,:)
1133                sshn(:,:)=0._wp
1134                e3t_n(:,:,:)=e3t_0(:,:,:)
1135                e3t_b(:,:,:)=e3t_0(:,:,:)
1136               !
1137            END IF           ! end of ll_wd edits
1138
1139            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1140               tilde_e3t_b(:,:,:) = 0._wp
1141               tilde_e3t_n(:,:,:) = 0._wp
1142               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
1143            END IF
1144         ENDIF
1145         !
1146      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1147         !                                   ! ===================
1148         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1149         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1150         !                                           ! --------- !
1151         !                                           ! all cases !
1152         !                                           ! --------- !
1153         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
1154         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
1155         !                                           ! ----------------------- !
1156         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
1157            !                                        ! ----------------------- !
1158            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
1159            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
1160         END IF
1161         !                                           ! -------------!   
1162         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
1163            !                                        ! ------------ !
1164            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
1165         ENDIF
1166         !
1167         IF( lwxios ) CALL iom_swap(      cxios_context          )
1168      ENDIF
1169      !
1170   END SUBROUTINE dom_vvl_rst
1171
1172
1173   SUBROUTINE dom_vvl_ctl
1174      !!---------------------------------------------------------------------
1175      !!                  ***  ROUTINE dom_vvl_ctl  ***
1176      !!               
1177      !! ** Purpose :   Control the consistency between namelist options
1178      !!                for vertical coordinate
1179      !!----------------------------------------------------------------------
1180      INTEGER ::   ioptio, ios
1181      !!
1182      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1183         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
1184         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
1185      !!----------------------------------------------------------------------
1186      !
1187      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
1188      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1189901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
1190      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
1191      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1192902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1193      IF(lwm) WRITE ( numond, nam_vvl )
1194      !
1195      IF(lwp) THEN                    ! Namelist print
1196         WRITE(numout,*)
1197         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1198         WRITE(numout,*) '~~~~~~~~~~~'
1199         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1200         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1201         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1202         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1203         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1204         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1205         WRITE(numout,*) '      !'
1206         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1207         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1208         IF( ln_vvl_ztilde_as_zstar ) THEN
1209            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1210            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1211            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1212            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1213            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1214            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1215         ELSE
1216            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1217            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1218         ENDIF
1219         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1220      ENDIF
1221      !
1222      ioptio = 0                      ! Parameter control
1223      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1224      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1225      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1226      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1227      !
1228      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1229      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1230      !
1231      IF(lwp) THEN                   ! Print the choice
1232         WRITE(numout,*)
1233         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1234         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1235         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1236         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1237      ENDIF
1238      !
1239#if defined key_agrif
1240      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1241#endif
1242      !
1243   END SUBROUTINE dom_vvl_ctl
1244
1245   !!======================================================================
1246END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.