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/trunk/src/SWE – NEMO

source: NEMO/trunk/src/SWE/domvvl.F90 @ 13237

Last change on this file since 13237 was 12983, checked in by techene, 4 years ago

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