source: NEMO/branches/2019/fix_vvl_ticket1791/tests/VORTEX/MY_SRC/domvvl.F90 @ 11012

Last change on this file since 11012 was 11012, checked in by jchanut, 17 months ago

#1791, report domvvl change in VORTEX test case. Modification of domvvl there should be removed.

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