New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90 @ 14071

Last change on this file since 14071 was 13416, checked in by gm, 4 years ago

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