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

Last change on this file since 13970 was 13970, checked in by andmirek, 3 years ago

Ticket #2462 into the trunk
File size: 68.6 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( 1, 1, 1, 1, 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( 1, 1, 1, 1 )
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 + nn_hls - 1   ;   ii1 = 111 + nn_hls - 1     
279                  ij0 = 128 + nn_hls       ;   ij1 = 135 + nn_hls
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( 1, 1, 1, 1, 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(:,:,:), ALLOCATABLE ::   ze3t
418      LOGICAL , DIMENSION(:,:,:), ALLOCATABLE ::   llmsk
419      !!----------------------------------------------------------------------
420      !
421      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
422      !
423      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
424      !
425      IF( kt == nit000 ) THEN
426         IF(lwp) WRITE(numout,*)
427         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
428         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
429      ENDIF
430
431      ll_do_bclinic = .TRUE.
432      IF( PRESENT(kcall) ) THEN
433         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
434      ENDIF
435
436      ! ******************************* !
437      ! After acale factors at t-points !
438      ! ******************************* !
439      !                                                ! --------------------------------------------- !
440      !                                                ! z_star coordinate and barotropic z-tilde part !
441      !                                                ! --------------------------------------------- !
442      !
443      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
444      DO jk = 1, jpkm1
445         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
446         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
447      END DO
448      !
449      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
450         !                                                               ! ------baroclinic part------ !
451         ! I - initialization
452         ! ==================
453
454         ! 1 - barotropic divergence
455         ! -------------------------
456         zhdiv(:,:) = 0._wp
457         zht(:,:)   = 0._wp
458         DO jk = 1, jpkm1
459            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
460            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
461         END DO
462         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
463
464         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
465         ! --------------------------------------------------
466         IF( ln_vvl_ztilde ) THEN
467            IF( kt > nit000 ) THEN
468               DO jk = 1, jpkm1
469                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
470                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
471               END DO
472            ENDIF
473         ENDIF
474
475         ! II - after z_tilde increments of vertical scale factors
476         ! =======================================================
477         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
478
479         ! 1 - High frequency divergence term
480         ! ----------------------------------
481         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
482            DO jk = 1, jpkm1
483               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
484            END DO
485         ELSE                         ! layer case
486            DO jk = 1, jpkm1
487               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
488            END DO
489         ENDIF
490
491         ! 2 - Restoring term (z-tilde case only)
492         ! ------------------
493         IF( ln_vvl_ztilde ) THEN
494            DO jk = 1, jpk
495               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
496            END DO
497         ENDIF
498
499         ! 3 - Thickness diffusion term
500         ! ----------------------------
501         zwu(:,:) = 0._wp
502         zwv(:,:) = 0._wp
503         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
504            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
505               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
506            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
507               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
508            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
509            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
510         END_3D
511         DO_2D( 1, 1, 1, 1 )
512            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
513            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
514         END_2D
515         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
516            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
517               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
518               &                                            ) * r1_e1e2t(ji,jj)
519         END_3D
520         !                       ! d - thickness diffusion transport: boundary conditions
521         !                             (stored for tracer advction and continuity equation)
522         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
523
524         ! 4 - Time stepping of baroclinic scale factors
525         ! ---------------------------------------------
526         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
527         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
528
529         ! Maximum deformation control
530         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
531         ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) )
532         ze3t(:,:,jpk) = 0._wp
533         DO jk = 1, jpkm1
534            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
535         END DO
536         !
537         llmsk(   1:Nis1,:,:) = .FALSE.   ! exclude halos from the checked region
538         llmsk(Nie1: jpi,:,:) = .FALSE.
539         llmsk(:,   1:Njs1,:) = .FALSE.
540         llmsk(:,Nje1: jpj,:) = .FALSE.
541         !
542         llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp                  ! define only the inner domain
543         z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_max( 'domvvl', z_tmax )   ! max over the global domain
544         z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_min( 'domvvl', z_tmin )   ! min over the global domain
545         ! - ML - test: for the moment, stop simulation for too large e3_t variations
546         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
547            CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max )
548            CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min )
549            IF (lwp) THEN
550               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
551               WRITE(numout, *) 'at i, j, k=', ijk_max
552               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
553               WRITE(numout, *) 'at i, j, k=', ijk_min           
554               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
555            ENDIF
556         ENDIF
557         DEALLOCATE( ze3t, llmsk )
558         ! - ML - end test
559         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
560         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
561         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
562
563         !
564         ! "tilda" change in the after scale factor
565         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
566         DO jk = 1, jpkm1
567            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
568         END DO
569         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
570         ! ==================================================================================
571         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
572         ! - ML - baroclinicity error should be better treated in the future
573         !        i.e. locally and not spread over the water column.
574         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
575         zht(:,:) = 0.
576         DO jk = 1, jpkm1
577            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
578         END DO
579         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
580         DO jk = 1, jpkm1
581            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
582         END DO
583
584      ENDIF
585
586      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
587      !                                           ! ---baroclinic part--------- !
588         DO jk = 1, jpkm1
589            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
590         END DO
591      ENDIF
592
593      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
594         !
595         IF( lwp ) WRITE(numout, *) 'kt =', kt
596         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
597            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
598            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
599            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
600         END IF
601         !
602         zht(:,:) = 0.0_wp
603         DO jk = 1, jpkm1
604            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
605         END DO
606         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
607         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
608         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
609         !
610         zht(:,:) = 0.0_wp
611         DO jk = 1, jpkm1
612            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
613         END DO
614         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
615         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
616         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
617         !
618         zht(:,:) = 0.0_wp
619         DO jk = 1, jpkm1
620            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
621         END DO
622         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
623         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
624         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
625         !
626         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
627         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
628         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
629         !
630         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
631         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
632         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
633         !
634         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
635         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
636         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
637      END IF
638
639      ! *********************************** !
640      ! After scale factors at u- v- points !
641      ! *********************************** !
642
643      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
644      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
645
646      ! *********************************** !
647      ! After depths at u- v points         !
648      ! *********************************** !
649
650      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
651      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
652      DO jk = 2, jpkm1
653         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
654         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
655      END DO
656      !                                        ! Inverse of the local depth
657!!gm BUG ?  don't understand the use of umask_i here .....
658      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
659      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
660      !
661      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
662      !
663   END SUBROUTINE dom_vvl_sf_nxt
664
665
666
667   SUBROUTINE dom_vvl_sf_nxt_st( kt, Kbb, Kmm, Kaa, kcall ) 
668      !!----------------------------------------------------------------------
669      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
670      !!                   
671      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
672      !!                 tranxt and dynspg routines
673      !!
674      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
675      !!               - z_tilde_case: after scale factor increment =
676      !!                                    high frequency part of horizontal divergence
677      !!                                  + retsoring towards the background grid
678      !!                                  + thickness difusion
679      !!                               Then repartition of ssh INCREMENT proportionnaly
680      !!                               to the "baroclinic" level thickness.
681      !!
682      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
683      !!               - tilde_e3t_a: after increment of vertical scale factor
684      !!                              in z_tilde case
685      !!               - e3(t/u/v)_a
686      !!
687      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
688      !!----------------------------------------------------------------------
689      INTEGER, INTENT( in )           ::   kt             ! time step
690      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
691      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
692      !
693      INTEGER                ::   ji, jj, jk            ! dummy loop indices
694      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
695      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
696      LOGICAL                ::   ll_do_bclinic         ! local logical
697      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
698      !!----------------------------------------------------------------------
699      !
700      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
701      !
702      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
703      !
704      IF( kt == nit000 ) THEN
705         IF(lwp) WRITE(numout,*)
706         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
707         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
708      ENDIF
709
710      ll_do_bclinic = .TRUE.
711      IF( PRESENT(kcall) ) THEN
712         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
713      ENDIF
714
715      ! ******************************* !
716      ! After acale factors at t-points !
717      ! ******************************* !
718      !
719      DO jk = 1, jpkm1
720         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) )
721         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) )
722         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) )
723      END DO
724      !
725      ! *********************************** !
726      ! After scale factors at u- v- points !
727      ! *********************************** !
728     
729      !!st CALL dom_vvl_interpol_st( r3u(:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
730      !!st CALL dom_vvl_interpol_st( r3v(:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
731
732      ! *********************************** !
733      ! After depths at u- v points         !
734      ! *********************************** !
735
736      !!st hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
737      !!st hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
738      !!st DO jk = 2, jpkm1
739      !!st    hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
740      !!st    hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
741      !!st   
742      !!st END DO
743      hu(:,:,Kaa) = (hu_0(:,:)*(1._wp+r3u(:,:,Kaa)))
744      hv(:,:,Kaa) = (hv_0(:,:)*(1._wp+r3v(:,:,Kaa)))
745      !                                        ! Inverse of the local depth
746!!gm BUG ?  don't understand the use of umask_i here .....
747      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
748      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
749      !
750      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
751      !
752   END SUBROUTINE dom_vvl_sf_nxt_st
753   
754
755
756   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
757      !!----------------------------------------------------------------------
758      !!                ***  ROUTINE dom_vvl_sf_update  ***
759      !!                   
760      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
761      !!               compute all depths and related variables for next time step
762      !!               write outputs and restart file
763      !!
764      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
765      !!               - reconstruct scale factor at other grid points (interpolate)
766      !!               - recompute depths and water height fields
767      !!
768      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
769      !!               - Recompute:
770      !!                    e3(u/v)_b       
771      !!                    e3w(:,:,:,Kmm)           
772      !!                    e3(u/v)w_b     
773      !!                    e3(u/v)w_n     
774      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
775      !!                    h(u/v) and h(u/v)r
776      !!
777      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
778      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
779      !!----------------------------------------------------------------------
780      INTEGER, INTENT( in ) ::   kt              ! time step
781      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
782      !
783      INTEGER  ::   ji, jj, jk   ! dummy loop indices
784      REAL(wp) ::   zcoef        ! local scalar
785      !!----------------------------------------------------------------------
786      !
787      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
788      !
789      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
790      !
791      IF( kt == nit000 )   THEN
792         IF(lwp) WRITE(numout,*)
793         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
794         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
795      ENDIF
796      !
797      ! Time filter and swap of scale factors
798      ! =====================================
799      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
800      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
801         IF( l_1st_euler ) THEN
802            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
803         ELSE
804            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
805            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
806         ENDIF
807         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
808      ENDIF
809
810      ! Compute all missing vertical scale factor and depths
811      ! ====================================================
812      ! Horizontal scale factor interpolations
813      ! --------------------------------------
814      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
815      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
816     
817      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F'  )
818     
819      ! Vertical scale factor interpolations
820      CALL dom_vvl_interpol( ssh(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
821      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
822      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
823      CALL dom_vvl_interpol( ssh(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
824      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
825      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
826
827      ! t- and w- points depth (set the isf depth as it is in the initial step)
828      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
829      gdepw(:,:,1,Kmm) = 0.0_wp
830      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
831      DO_3D( 1, 1, 1, 1, 2, jpk )
832        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
833                                                           ! 1 for jk = mikt
834         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
835         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
836         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
837             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
838         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
839      END_3D
840
841      ! Local depth and Inverse of the local depth of the water
842      ! -------------------------------------------------------
843      !
844      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
845      DO jk = 2, jpkm1
846         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
847      END DO
848
849      ! write restart file
850      ! ==================
851      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
852      !
853      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
854      !
855   END SUBROUTINE dom_vvl_sf_update
856   
857
858   SUBROUTINE dom_vvl_sf_update_st( kt, Kbb, Kmm, Kaa )
859      !!----------------------------------------------------------------------
860      !!                ***  ROUTINE dom_vvl_sf_update  ***
861      !!                   
862      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
863      !!               compute all depths and related variables for next time step
864      !!               write outputs and restart file
865      !!
866      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
867      !!               - reconstruct scale factor at other grid points (interpolate)
868      !!               - recompute depths and water height fields
869      !!
870      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
871      !!               - Recompute:
872      !!                    e3(u/v)_b       
873      !!                    e3w(:,:,:,Kmm)           
874      !!                    e3(u/v)w_b     
875      !!                    e3(u/v)w_n     
876      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
877      !!                    h(u/v) and h(u/v)r
878      !!
879      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
880      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
881      !!----------------------------------------------------------------------
882      INTEGER, INTENT( in ) ::   kt              ! time step
883      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
884      !
885      INTEGER  ::   ji, jj, jk   ! dummy loop indices
886      REAL(wp) ::   zcoef        ! local scalar
887      !!----------------------------------------------------------------------
888      !
889      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
890      !
891      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
892      !
893      IF( kt == nit000 )   THEN
894         IF(lwp) WRITE(numout,*)
895         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
896         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
897      ENDIF
898      !
899
900      ! Compute all missing vertical scale factor and depths
901      ! ====================================================
902      ! Horizontal scale factor interpolations
903      ! --------------------------------------
904      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
905      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
906     
907      CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F'  )
908
909      ! Vertical scale factor interpolations
910      CALL dom_vvl_interpol_st( r3t(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
911      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
912      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
913      CALL dom_vvl_interpol_st( r3t(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
914      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
915      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
916
917      ! t- and w- points depth (set the isf depth as it is in the initial step)
918      DO_3D( 1, 1, 1, 1, 1, jpk )
919         gdepw(ji,jj,jk,Kmm) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm))
920         gdept(ji,jj,jk,Kmm) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 
921         gde3w(ji,jj,jk    ) =  gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
922      END_3D
923
924      ! Local depth and Inverse of the local depth of the water
925      ! -------------------------------------------------------
926      !
927      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
928
929      ! write restart file
930      ! ==================
931      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
932      !
933      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
934      !
935   END SUBROUTINE dom_vvl_sf_update_st
936   
937
938   
939   SUBROUTINE dom_vvl_interpol_st( rc3, pe3, cdp )
940      !!---------------------------------------------------------------------
941      !!                  ***  ROUTINE dom_vvl__interpol  ***
942      !!
943      !! ** Purpose :   interpolate scale factors from one grid point to another
944      !!
945      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
946      !!                - horizontal interpolation: grid cell surface averaging
947      !!                - vertical interpolation: simple averaging
948      !!----------------------------------------------------------------------
949      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   rc3     ! input e3   NOT used here (ssh is used instead)
950      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m]
951      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' )
952      !
953      INTEGER ::   ji, jj, jk                 ! dummy loop indices
954      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace
955      !!----------------------------------------------------------------------
956      !
957      SELECT CASE ( cdp )     !==  type of interpolation  ==!
958         !
959      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
960         DO jk = 1, jpkm1
961            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
962         END DO
963         !
964      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
965         DO jk = 1, jpkm1
966            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
967         END DO
968         !
969      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
970         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh
971            e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + rc3(:,:) )
972         END DO
973         !
974      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
975         DO jk = 1, jpk
976            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
977         END DO
978         !
979      CASE( 'UW' )                  !* from U- to UW-point
980         DO jk = 1, jpk
981            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
982         END DO
983      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean
984         DO jk = 1, jpk
985            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
986         END DO
987         !
988      END SELECT
989      !
990   END SUBROUTINE dom_vvl_interpol_st
991   
992
993   SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp )
994      !!---------------------------------------------------------------------
995      !!                  ***  ROUTINE dom_vvl__interpol  ***
996      !!
997      !! ** Purpose :   interpolate scale factors from one grid point to another
998      !!
999      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
1000      !!                - horizontal interpolation: grid cell surface averaging
1001      !!                - vertical interpolation: simple averaging
1002      !!----------------------------------------------------------------------
1003      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pssh    ! input e3   NOT used here (ssh is used instead)
1004      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m]
1005      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' )
1006      !
1007      INTEGER ::   ji, jj, jk                 ! dummy loop indices
1008      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace
1009      !!----------------------------------------------------------------------
1010      !
1011      SELECT CASE ( cdp )     !==  type of interpolation  ==!
1012         !
1013      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
1014         DO_2D( 0, 0, 0, 0 )
1015            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
1016               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
1017         END_2D
1018         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp )
1019         !
1020         DO jk = 1, jpkm1
1021            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1022         END DO
1023         !
1024      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
1025         DO_2D( 0, 0, 0, 0 )
1026            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
1027               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
1028         END_2D
1029         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp )
1030         !
1031         DO jk = 1, jpkm1
1032            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1033         END DO
1034         !
1035      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
1036         DO_2D( 1, 0, 1, 0 )
1037            zc3(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
1038               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
1039               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
1040               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
1041         END_2D
1042         CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp )
1043         !
1044         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh
1045            e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + zc3(:,:) )
1046         END DO
1047         !
1048      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
1049         zc3(:,:) = pssh(:,:) * r1_ht_0(:,:)
1050         !
1051         DO jk = 1, jpk
1052            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1053         END DO
1054         !
1055      CASE( 'UW' )                  !* from U- to UW-point
1056         !
1057         DO_2D( 0, 0, 0, 0 )
1058            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
1059               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
1060         END_2D
1061         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp )
1062         !
1063         DO jk = 1, jpk
1064            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1065         END DO
1066      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean
1067         !
1068         DO_2D( 0, 0, 0, 0 )
1069            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
1070               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
1071         END_2D
1072         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp )
1073          !
1074         DO jk = 1, jpk
1075            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1076         END DO
1077         !
1078      END SELECT
1079      !
1080   END SUBROUTINE dom_vvl_interpol
1081   
1082
1083   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
1084      !!---------------------------------------------------------------------
1085      !!                   ***  ROUTINE dom_vvl_rst  ***
1086      !!                     
1087      !! ** Purpose :   Read or write VVL file in restart file
1088      !!
1089      !! ** Method  :   use of IOM library
1090      !!                if the restart does not contain vertical scale factors,
1091      !!                they are set to the _0 values
1092      !!                if the restart does not contain vertical scale factors increments (z_tilde),
1093      !!                they are set to 0.
1094      !!----------------------------------------------------------------------
1095      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
1096      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
1097      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
1098      !
1099      INTEGER ::   ji, jj, jk
1100      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
1101      !!----------------------------------------------------------------------
1102      !
1103      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1104         !                                   ! ===============
1105         IF( ln_rstart ) THEN                   !* Read the restart file
1106            CALL rst_read_open                  !  open the restart file if necessary
1107            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm)    )
1108            !
1109            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
1110            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
1111            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1112            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1113            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
1114            !
1115            !                             ! --------- !
1116            !                             ! all cases !
1117            !                             ! --------- !
1118            !
1119            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
1120               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
1121               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
1122               ! needed to restart if land processor not computed
1123               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
1124               WHERE ( tmask(:,:,:) == 0.0_wp ) 
1125                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
1126                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
1127               END WHERE
1128               IF( l_1st_euler ) THEN
1129                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1130               ENDIF
1131            ELSE IF( id1 > 0 ) THEN
1132               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
1133               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
1134               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1135               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
1136               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
1137               l_1st_euler = .true.
1138            ELSE IF( id2 > 0 ) THEN
1139               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
1140               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
1141               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1142               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
1143               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1144               l_1st_euler = .true.
1145            ELSE
1146               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
1147               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1148               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1149               DO jk = 1, jpk
1150                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
1151                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1152                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
1153               END DO
1154               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1155               l_1st_euler = .true.
1156            ENDIF
1157            !                             ! ----------- !
1158            IF( ln_vvl_zstar ) THEN       ! z_star case !
1159               !                          ! ----------- !
1160               IF( MIN( id3, id4 ) > 0 ) THEN
1161                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1162               ENDIF
1163               !                          ! ----------------------- !
1164            ELSE                          ! z_tilde and layer cases !
1165               !                          ! ----------------------- !
1166               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
1167                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
1168                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
1169               ELSE                            ! one at least array is missing
1170                  tilde_e3t_b(:,:,:) = 0.0_wp
1171                  tilde_e3t_n(:,:,:) = 0.0_wp
1172               ENDIF
1173               !                          ! ------------ !
1174               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
1175                  !                       ! ------------ !
1176                  IF( id5 > 0 ) THEN  ! required array exists
1177                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) )
1178                  ELSE                ! array is missing
1179                     hdiv_lf(:,:,:) = 0.0_wp
1180                  ENDIF
1181               ENDIF
1182            ENDIF
1183            !
1184         ELSE                                   !* Initialize at "rest"
1185            !
1186
1187            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
1188               !
1189               IF( cn_cfg == 'wad' ) THEN
1190                  ! Wetting and drying test case
1191                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
1192!!an                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
1193                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
1194                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
1195                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
1196               ELSE
1197                  ! if not test case
1198                  ssh(:,:,Kmm) = -ssh_ref
1199                  ssh(:,:,Kbb) = -ssh_ref
1200
1201                  DO_2D( 1, 1, 1, 1 )
1202                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
1203                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
1204                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
1205                     ENDIF
1206                  END_2D
1207               ENDIF !If test case else
1208
1209               ! Adjust vertical metrics for all wad
1210               DO jk = 1, jpk
1211                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
1212                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1213                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1214               END DO
1215               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1216
1217               DO ji = 1, jpi
1218                  DO jj = 1, jpj
1219                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
1220                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
1221                     ENDIF
1222                  END DO
1223               END DO 
1224               !
1225            ELSE
1226               !
1227               ! Just to read set ssh in fact, called latter once vertical grid
1228               ! is set up:
1229!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
1230!               !
1231!               DO jk=1,jpk
1232!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
1233!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
1234!               END DO
1235!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
1236                ssh(:,:,Kmm)=0._wp
1237                ssh(:,:,Kbb)=0._wp
1238                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
1239                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
1240               !
1241            END IF           ! end of ll_wd edits
1242
1243            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1244               tilde_e3t_b(:,:,:) = 0._wp
1245               tilde_e3t_n(:,:,:) = 0._wp
1246               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
1247            END IF
1248         ENDIF
1249         !
1250      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1251         !                                   ! ===================
1252         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1253         !                                           ! --------- !
1254         !                                           ! all cases !
1255         !                                           ! --------- !
1256         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb) )
1257         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm) )
1258         !                                           ! ----------------------- !
1259         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
1260            !                                        ! ----------------------- !
1261            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:))
1262            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:))
1263         END IF
1264         !                                           ! -------------!   
1265         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
1266            !                                        ! ------------ !
1267            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:))
1268         ENDIF
1269         !
1270      ENDIF
1271      !
1272   END SUBROUTINE dom_vvl_rst
1273
1274
1275   SUBROUTINE dom_vvl_ctl
1276      !!---------------------------------------------------------------------
1277      !!                  ***  ROUTINE dom_vvl_ctl  ***
1278      !!               
1279      !! ** Purpose :   Control the consistency between namelist options
1280      !!                for vertical coordinate
1281      !!----------------------------------------------------------------------
1282      INTEGER ::   ioptio, ios
1283      !!
1284      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1285         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
1286         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
1287      !!----------------------------------------------------------------------
1288      !
1289      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1290901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
1291      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1292902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
1293      IF(lwm) WRITE ( numond, nam_vvl )
1294      !
1295      IF(lwp) THEN                    ! Namelist print
1296         WRITE(numout,*)
1297         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1298         WRITE(numout,*) '~~~~~~~~~~~'
1299         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1300         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1301         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1302         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1303         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1304         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1305         WRITE(numout,*) '      !'
1306         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1307         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1308         IF( ln_vvl_ztilde_as_zstar ) THEN
1309            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1310            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1311            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1312            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1313            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1314            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
1315         ELSE
1316            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1317            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1318         ENDIF
1319         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1320      ENDIF
1321      !
1322      ioptio = 0                      ! Parameter control
1323      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1324      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1325      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1326      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1327      !
1328      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1329      !
1330      IF(lwp) THEN                   ! Print the choice
1331         WRITE(numout,*)
1332         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1333         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1334         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1335         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1336      ENDIF
1337      !
1338#if defined key_agrif
1339      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1340#endif
1341      !
1342   END SUBROUTINE dom_vvl_ctl
1343
1344#endif
1345!!stoops
1346
1347   !!======================================================================
1348END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.