source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4296

Last change on this file since 4296 was 4296, checked in by cetlod, 7 years ago

dev_MERGE_2013 : minor change in domvvl.F90

  • Property svn:keywords set to Id
File size: 69.7 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:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
15   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
16   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
17   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
18   !!   dom_vvl_rst      : read/write restart file
19   !!   dom_vvl_ctl      : Check the vvl options
20   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
21   !!                    : to account for manual changes to e[1,2][u,v] in some Straits
22   !!----------------------------------------------------------------------
23   !! * Modules used
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! ocean surface boundary condition
27   USE in_out_manager  ! I/O manager
28   USE iom             ! I/O manager library
29   USE restart         ! ocean restart
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC  dom_vvl_init       ! called by domain.F90
40   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
41   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
42   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
43   PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol
44
45   !!* Namelist nam_vvl
46   LOGICAL , PUBLIC                                      :: ln_vvl_zstar              ! zstar  vertical coordinate
47   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde             ! ztilde vertical coordinate
48   LOGICAL , PUBLIC                                      :: ln_vvl_layer              ! level  vertical coordinate
49   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor     ! ztilde vertical coordinate
51   LOGICAL , PUBLIC                                      :: ln_vvl_kepe               ! kinetic/potential energy transfer
52   !                                                                                           ! conservation: not used yet
53   REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient
54   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
55   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
56   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
57   LOGICAL , PUBLIC                                      :: ln_vvl_dbg                ! debug control prints
58
59   !! * Module variables
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
62   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a                        ! baroclinic scale factors
64   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
65   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
66
67   !! * Substitutions
68#  include "domzgr_substitute.h90"
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75
76CONTAINS
77
78   INTEGER FUNCTION dom_vvl_alloc()
79      !!----------------------------------------------------------------------
80      !!                ***  FUNCTION dom_vvl_alloc  ***
81      !!----------------------------------------------------------------------
82      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
83      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
84         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
85            &      un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     , STAT = dom_vvl_alloc        )
86         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
87         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
88         un_td = 0.0_wp
89         vn_td = 0.0_wp
90      ENDIF
91      IF( ln_vvl_ztilde ) THEN
92         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
93         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
94         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
95      ENDIF
96
97   END FUNCTION dom_vvl_alloc
98
99
100   SUBROUTINE dom_vvl_init
101      !!----------------------------------------------------------------------
102      !!                ***  ROUTINE dom_vvl_init  ***
103      !!                   
104      !! ** Purpose :  Initialization of all scale factors, depths
105      !!               and water column heights
106      !!
107      !! ** Method  :  - use restart file and/or initialize
108      !!               - interpolate scale factors
109      !!
110      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
111      !!              - Regrid: fse3(u/v)_n
112      !!                        fse3(u/v)_b       
113      !!                        fse3w_n           
114      !!                        fse3(u/v)w_b     
115      !!                        fse3(u/v)w_n     
116      !!                        fsdept_n, fsdepw_n and fsde3w_n
117      !!              - h(t/u/v)_0
118      !!              - frq_rst_e3t and frq_rst_hdv
119      !!
120      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
121      !!----------------------------------------------------------------------
122      USE phycst,  ONLY : rpi, rsmall, rad
123      !! * Local declarations
124      INTEGER ::   ji,jj,jk
125      INTEGER ::   ii0, ii1, ij0, ij1
126      !!----------------------------------------------------------------------
127      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
128
129      IF(lwp) WRITE(numout,*)
130      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
131      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
132
133      ! choose vertical coordinate (z_star, z_tilde or layer)
134      ! ==========================
135      CALL dom_vvl_ctl
136
137      ! Allocate module arrays
138      ! ======================
139      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
140
141      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
142      ! ============================================================================
143      CALL dom_vvl_rst( nit000, 'READ' )
144      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
145
146      ! Reconstruction of all vertical scale factors at now and before time steps
147      ! =============================================================================
148      ! Horizontal scale factor interpolations
149      ! --------------------------------------
150      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
151      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
152      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
153      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
154      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
155      ! Vertical scale factor interpolations
156      ! ------------------------------------
157      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
158      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
159      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
160      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
161      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
162      ! t- and w- points depth
163      ! ----------------------
164      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
165      fsdepw_n(:,:,1) = 0.0_wp
166      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
167      DO jk = 2, jpk
168         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
169         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
170         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
171      END DO
172      ! Reference water column height at t-, u- and v- point
173      ! ----------------------------------------------------
174      ht_0(:,:) = 0.0_wp
175      hu_0(:,:) = 0.0_wp
176      hv_0(:,:) = 0.0_wp
177      DO jk = 1, jpk
178         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
179         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
180         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
181      END DO
182
183      ! Restoring frequencies for z_tilde coordinate
184      ! ============================================
185      IF( ln_vvl_ztilde ) THEN
186         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
187         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
188         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
189         IF( ln_vvl_ztilde_as_zstar ) THEN
190            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
191            frq_rst_e3t(:,:) = 0.0_wp 
192            frq_rst_hdv(:,:) = 1.0_wp / rdt
193         ENDIF
194         IF ( ln_vvl_zstar_at_eqtor ) THEN
195            DO jj = 1, jpj
196               DO ji = 1, jpi
197                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
198                     ! values outside the equatorial band and transition zone (ztilde)
199                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
200                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
201                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
202                     ! values inside the equatorial band (ztilde as zstar)
203                     frq_rst_e3t(ji,jj) =  0.0_wp
204                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
205                  ELSE
206                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
207                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
208                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
209                        &                                          * 180._wp / 3.5_wp ) )
210                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
211                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
212                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
213                        &                                          * 180._wp / 3.5_wp ) )
214                  ENDIF
215               END DO
216            END DO
217            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
218               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
219               ij0 = 128   ;   ij1 = 135   ;   
220               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
221               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
222            ENDIF
223         ENDIF
224      ENDIF
225
226      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
227
228   END SUBROUTINE dom_vvl_init
229
230
231   SUBROUTINE dom_vvl_sf_nxt( kt ) 
232      !!----------------------------------------------------------------------
233      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
234      !!                   
235      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
236      !!                 tranxt and dynspg routines
237      !!
238      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
239      !!               - z_tilde_case: after scale factor increment =
240      !!                                    high frequency part of horizontal divergence
241      !!                                  + retsoring towards the background grid
242      !!                                  + thickness difusion
243      !!                               Then repartition of ssh INCREMENT proportionnaly
244      !!                               to the "baroclinic" level thickness.
245      !!
246      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
247      !!               - tilde_e3t_a: after increment of vertical scale factor
248      !!                              in z_tilde case
249      !!               - fse3(t/u/v)_a
250      !!
251      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
252      !!----------------------------------------------------------------------
253      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
254      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
255      !! * Arguments
256      INTEGER, INTENT( in )                  :: kt                    ! time step
257      !! * Local declarations
258      INTEGER                                :: ji, jj, jk            ! dummy loop indices
259      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
260      REAL(wp)                               :: z2dt                  ! temporary scalars
261      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
262      !!----------------------------------------------------------------------
263      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
264      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
265      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
266
267      IF(kt == nit000)   THEN
268         IF(lwp) WRITE(numout,*)
269         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
270         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
271      ENDIF
272
273      ! ******************************* !
274      ! After acale factors at t-points !
275      ! ******************************* !
276
277      !                                                ! ----------------- !
278      IF( ln_vvl_zstar ) THEN                          ! z_star coordinate !
279         !                                             ! ----------------- !
280
281         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
282         DO jk = 1, jpkm1
283            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
284         END DO
285
286      !                                                ! --------------------------- !
287      ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate !
288         !                                             ! --------------------------- !
289
290         ! I - initialization
291         ! ==================
292
293         ! 1 - barotropic divergence
294         ! -------------------------
295         zhdiv(:,:) = 0.
296         zht(:,:)   = 0.
297         DO jk = 1, jpkm1
298            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
299            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
300         END DO
301         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
302
303         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
304         ! --------------------------------------------------
305         IF( ln_vvl_ztilde ) THEN
306            IF( kt .GT. nit000 ) THEN
307               DO jk = 1, jpkm1
308                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
309                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
310               END DO
311            ENDIF
312         END IF
313
314         ! II - after z_tilde increments of vertical scale factors
315         ! =======================================================
316         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
317
318         ! 1 - High frequency divergence term
319         ! ----------------------------------
320         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
321            DO jk = 1, jpkm1
322               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
323            END DO
324         ELSE                         ! layer case
325            DO jk = 1, jpkm1
326               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
327            END DO
328         END IF
329
330         ! 2 - Restoring term (z-tilde case only)
331         ! ------------------
332         IF( ln_vvl_ztilde ) THEN
333            DO jk = 1, jpk
334               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
335            END DO
336         END IF
337
338         ! 3 - Thickness diffusion term
339         ! ----------------------------
340         zwu(:,:) = 0.0_wp
341         zwv(:,:) = 0.0_wp
342         ! a - first derivative: diffusive fluxes
343         DO jk = 1, jpkm1
344            DO jj = 1, jpjm1
345               DO ji = 1, fs_jpim1   ! vector opt.
346                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
347                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
348                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
349                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
350                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
351                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
352               END DO
353            END DO
354         END DO
355         ! b - correction for last oceanic u-v points
356         DO jj = 1, jpj
357            DO ji = 1, jpi
358               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
359               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
360            END DO
361         END DO
362         ! c - second derivative: divergence of diffusive fluxes
363         DO jk = 1, jpkm1
364            DO jj = 2, jpjm1
365               DO ji = fs_2, fs_jpim1   ! vector opt.
366                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
367                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
368                     &                                            ) * r1_e12t(ji,jj)
369               END DO
370            END DO
371         END DO
372         ! d - thickness diffusion transport: boundary conditions
373         !     (stored for tracer advction and continuity equation)
374         CALL lbc_lnk( un_td , 'U' , -1.)
375         CALL lbc_lnk( vn_td , 'V' , -1.)
376
377         ! 4 - Time stepping of baroclinic scale factors
378         ! ---------------------------------------------
379         ! Leapfrog time stepping
380         ! ~~~~~~~~~~~~~~~~~~~~~~
381         IF( neuler == 0 .AND. kt == nit000 ) THEN
382            z2dt =  rdt
383         ELSE
384            z2dt = 2.0_wp * rdt
385         ENDIF
386         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
387         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
388
389         ! Maximum deformation control
390         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
391         ze3t(:,:,jpk) = 0.0_wp
392         DO jk = 1, jpkm1
393            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
394         END DO
395         z_tmax = MAXVAL( ze3t(:,:,:) )
396         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
397         z_tmin = MINVAL( ze3t(:,:,:) )
398         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
399         ! - ML - test: for the moment, stop simulation for too large e3_t variations
400         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
401            IF( lk_mpp ) THEN
402               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
403               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
404            ELSE
405               ijk_max = MAXLOC( ze3t(:,:,:) )
406               ijk_max(1) = ijk_max(1) + nimpp - 1
407               ijk_max(2) = ijk_max(2) + njmpp - 1
408               ijk_min = MINLOC( ze3t(:,:,:) )
409               ijk_min(1) = ijk_min(1) + nimpp - 1
410               ijk_min(2) = ijk_min(2) + njmpp - 1
411            ENDIF
412            IF (lwp) THEN
413               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
414               WRITE(numout, *) 'at i, j, k=', ijk_max
415               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
416               WRITE(numout, *) 'at i, j, k=', ijk_min           
417               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
418            ENDIF
419         ENDIF
420         ! - ML - end test
421         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
422         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
423         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
424
425         ! Add "tilda" part to the after scale factor
426         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427         fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:)
428
429         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
430         ! ==================================================================================
431         ! add e3t(n-1) "star" Asselin-filtered
432         DO jk = 1, jpkm1
433            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk)
434         END DO
435         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
436         ! - ML - baroclinicity error should be better treated in the future
437         !        i.e. locally and not spread over the water column.
438         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
439         zht(:,:) = 0.
440         DO jk = 1, jpkm1
441            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
442         END DO
443         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
444         DO jk = 1, jpkm1
445            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
446         END DO
447
448      ENDIF
449
450      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
451         !
452         IF( lwp ) WRITE(numout, *) 'kt =', kt
453         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
454            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
455            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
456            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
457         END IF
458         !
459         zht(:,:) = 0.0_wp
460         DO jk = 1, jpkm1
461            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
462         END DO
463         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
464         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
465         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
466         !
467         zht(:,:) = 0.0_wp
468         DO jk = 1, jpkm1
469            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
470         END DO
471         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
472         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
473         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
474         !
475         zht(:,:) = 0.0_wp
476         DO jk = 1, jpkm1
477            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
478         END DO
479         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
480         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
481         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
482         !
483         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
484         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
485         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
486         !
487         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
488         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
489         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
490         !
491         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
492         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
493         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
494      END IF
495
496      ! *********************************** !
497      ! After scale factors at u- v- points !
498      ! *********************************** !
499
500      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
501      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
502
503      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
504      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
505
506      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
507
508   END SUBROUTINE dom_vvl_sf_nxt
509
510
511   SUBROUTINE dom_vvl_sf_swp( kt )
512      !!----------------------------------------------------------------------
513      !!                ***  ROUTINE dom_vvl_sf_swp  ***
514      !!                   
515      !! ** Purpose :  compute time filter and swap of scale factors
516      !!               compute all depths and related variables for next time step
517      !!               write outputs and restart file
518      !!
519      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
520      !!               - reconstruct scale factor at other grid points (interpolate)
521      !!               - recompute depths and water height fields
522      !!
523      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
524      !!               - Recompute:
525      !!                    fse3(u/v)_b       
526      !!                    fse3w_n           
527      !!                    fse3(u/v)w_b     
528      !!                    fse3(u/v)w_n     
529      !!                    fsdept_n, fsdepw_n  and fsde3w_n
530      !!                    h(u/v) and h(u/v)r
531      !!
532      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
533      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
534      !!----------------------------------------------------------------------
535      !! * Arguments
536      INTEGER, INTENT( in )               :: kt       ! time step
537      !! * Local declarations
538      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
539      INTEGER                             :: jk       ! dummy loop indices
540      !!----------------------------------------------------------------------
541
542      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
543      !
544      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
545      !
546      IF( kt == nit000 )   THEN
547         IF(lwp) WRITE(numout,*)
548         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
549         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
550      ENDIF
551      !
552      ! Time filter and swap of scale factors
553      ! =====================================
554      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
555      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
556         IF( neuler == 0 .AND. kt == nit000 ) THEN
557            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
558         ELSE
559            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
560            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
561         ENDIF
562         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
563      ENDIF
564      fse3t_n(:,:,:) = fse3t_a(:,:,:)
565      fse3u_n(:,:,:) = fse3u_a(:,:,:)
566      fse3v_n(:,:,:) = fse3v_a(:,:,:)
567
568      ! Compute all missing vertical scale factor and depths
569      ! ====================================================
570      ! Horizontal scale factor interpolations
571      ! --------------------------------------
572      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
573      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
574      ! Vertical scale factor interpolations
575      ! ------------------------------------
576      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
577      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
578      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
579      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
580      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
581      ! t- and w- points depth
582      ! ----------------------
583      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
584      fsdepw_n(:,:,1) = 0.0_wp
585      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
586      DO jk = 2, jpk
587         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
588         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
589         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
590      END DO
591      ! Local depth and Inverse of the local depth of the water column at u- and v- points
592      ! ----------------------------------------------------------------------------------
593      hu(:,:) = 0.
594      hv(:,:) = 0.
595      DO jk = 1, jpk
596         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
597         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
598      END DO
599      ! Inverse of the local depth
600      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
601      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) )
602
603      ! Write outputs
604      ! =============
605      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
606      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
607      CALL iom_put( "dept_n" , fsde3w_n (:,:,:) )
608      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
609
610      ! write restart file
611      ! ==================
612      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
613      !
614      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
615      !
616      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
617
618   END SUBROUTINE dom_vvl_sf_swp
619
620
621   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
622      !!---------------------------------------------------------------------
623      !!                  ***  ROUTINE dom_vvl__interpol  ***
624      !!
625      !! ** Purpose :   interpolate scale factors from one grid point to another
626      !!
627      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
628      !!                - horizontal interpolation: grid cell surface averaging
629      !!                - vertical interpolation: simple averaging
630      !!----------------------------------------------------------------------
631      !! * Arguments
632      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
633      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
634      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
635      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
636      !! * Local declarations
637      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
638      LOGICAL ::   l_is_orca                                           ! local logical
639      !!----------------------------------------------------------------------
640      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
641         !
642      l_is_orca = .FALSE.
643      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
644
645      SELECT CASE ( pout )
646         !               ! ------------------------------------- !
647      CASE( 'U' )        ! interpolation from T-point to U-point !
648         !               ! ------------------------------------- !
649         ! horizontal surface weighted interpolation
650         DO jk = 1, jpk
651            DO jj = 1, jpjm1
652               DO ji = 1, fs_jpim1   ! vector opt.
653                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
654                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
655                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
656               END DO
657            END DO
658         END DO
659         !
660         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
661         ! boundary conditions
662         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
663         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
664         !               ! ------------------------------------- !
665      CASE( 'V' )        ! interpolation from T-point to V-point !
666         !               ! ------------------------------------- !
667         ! horizontal surface weighted interpolation
668         DO jk = 1, jpk
669            DO jj = 1, jpjm1
670               DO ji = 1, fs_jpim1   ! vector opt.
671                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
672                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
673                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
674               END DO
675            END DO
676         END DO
677         !
678         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
679         ! boundary conditions
680         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
681         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
682         !               ! ------------------------------------- !
683      CASE( 'F' )        ! interpolation from U-point to F-point !
684         !               ! ------------------------------------- !
685         ! horizontal surface weighted interpolation
686         DO jk = 1, jpk
687            DO jj = 1, jpjm1
688               DO ji = 1, fs_jpim1   ! vector opt.
689                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
690                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
691                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
692               END DO
693            END DO
694         END DO
695         !
696         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
697         ! boundary conditions
698         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
699         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
700         !               ! ------------------------------------- !
701      CASE( 'W' )        ! interpolation from T-point to W-point !
702         !               ! ------------------------------------- !
703         ! vertical simple interpolation
704         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
705         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
706         DO jk = 2, jpk
707            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
708               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
709         END DO
710         !               ! -------------------------------------- !
711      CASE( 'UW' )       ! interpolation from U-point to UW-point !
712         !               ! -------------------------------------- !
713         ! vertical simple interpolation
714         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
715         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
716         DO jk = 2, jpk
717            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
718               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
719         END DO
720         !               ! -------------------------------------- !
721      CASE( 'VW' )       ! interpolation from V-point to VW-point !
722         !               ! -------------------------------------- !
723         ! vertical simple interpolation
724         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
725         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
726         DO jk = 2, jpk
727            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
728               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
729         END DO
730      END SELECT
731      !
732
733      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
734
735   END SUBROUTINE dom_vvl_interpol
736
737   SUBROUTINE dom_vvl_rst( kt, cdrw )
738      !!---------------------------------------------------------------------
739      !!                   ***  ROUTINE dom_vvl_rst  ***
740      !!                     
741      !! ** Purpose :   Read or write VVL file in restart file
742      !!
743      !! ** Method  :   use of IOM library
744      !!                if the restart does not contain vertical scale factors,
745      !!                they are set to the _0 values
746      !!                if the restart does not contain vertical scale factors increments (z_tilde),
747      !!                they are set to 0.
748      !!----------------------------------------------------------------------
749      !! * Arguments
750      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
751      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
752      !! * Local declarations
753      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
754      !!----------------------------------------------------------------------
755      !
756      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
757      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
758         !                                   ! ===============
759         IF( ln_rstart ) THEN                   !* Read the restart file
760            CALL rst_read_open                  !  open the restart file if necessary
761            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
762            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
763            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
764            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
765            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
766            !                             ! --------- !
767            !                             ! all cases !
768            !                             ! --------- !
769            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
770               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
771               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
772               IF( neuler == 0 ) THEN
773                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
774               ENDIF
775            ELSE IF( id1 > 0 ) THEN
776               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
777               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
778               fse3t_b(:,:,:) = fse3t_n(:,:,:)
779            ELSE                                 ! one at least array is missing
780               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
781            ENDIF
782            !                             ! ----------- !
783            IF( ln_vvl_zstar ) THEN       ! z_star case !
784               !                          ! ----------- !
785               IF( MIN( id3, id4 ) > 0 ) THEN
786                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
787               ENDIF
788               !                          ! ----------------------- !
789            ELSE                          ! z_tilde and layer cases !
790               !                          ! ----------------------- !
791               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
792                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
793                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
794               ELSE                            ! one at least array is missing
795                  tilde_e3t_b(:,:,:) = 0.0_wp
796                  tilde_e3t_n(:,:,:) = 0.0_wp
797               ENDIF
798               !                          ! ------------ !
799               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
800                  !                       ! ------------ !
801                  IF( id5 > 0 ) THEN  ! required array exists
802                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
803                  ELSE                ! array is missing
804                     hdiv_lf(:,:,:) = 0.0_wp
805                  ENDIF
806               ENDIF
807            ENDIF
808            !
809         ELSE                                   !* Initialize at "rest"
810            fse3t_b(:,:,:) = e3t_0(:,:,:)
811            fse3t_n(:,:,:) = e3t_0(:,:,:)
812            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
813               tilde_e3t_b(:,:,:) = 0.0_wp
814               tilde_e3t_n(:,:,:) = 0.0_wp
815               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
816            END IF
817         ENDIF
818
819      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
820         !                                   ! ===================
821         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
822         !                                           ! --------- !
823         !                                           ! all cases !
824         !                                           ! --------- !
825         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
826         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
827         !                                           ! ----------------------- !
828         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
829            !                                        ! ----------------------- !
830            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
831            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
832         END IF
833         !                                           ! -------------!   
834         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
835            !                                        ! ------------ !
836            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
837         ENDIF
838
839      ENDIF
840      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
841
842   END SUBROUTINE dom_vvl_rst
843
844
845   SUBROUTINE dom_vvl_ctl
846      !!---------------------------------------------------------------------
847      !!                  ***  ROUTINE dom_vvl_ctl  ***
848      !!               
849      !! ** Purpose :   Control the consistency between namelist options
850      !!                for vertical coordinate
851      !!----------------------------------------------------------------------
852      INTEGER ::   ioptio
853      INTEGER ::   ios
854
855      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
856                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
857                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
858      !!----------------------------------------------------------------------
859
860      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
861      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
862901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
863
864      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
865      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
866902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
867      WRITE ( numond, nam_vvl )
868
869      IF(lwp) THEN                    ! Namelist print
870         WRITE(numout,*)
871         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
872         WRITE(numout,*) '~~~~~~~~~~~'
873         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
874         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
875         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
876         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
877         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
878         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
879         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
880         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
881         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
882         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
883         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
884         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
885         IF( ln_vvl_ztilde_as_zstar ) THEN
886            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
887            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
888            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
889            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
890            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
891            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
892         ELSE
893            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
894            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
895            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
896            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
897         ENDIF
898         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
899         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
900      ENDIF
901
902      ioptio = 0                      ! Parameter control
903      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
904      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
905      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
906      IF( ln_vvl_layer           )        ioptio = ioptio + 1
907
908      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
909
910      IF(lwp) THEN                   ! Print the choice
911         WRITE(numout,*)
912         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
913         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
914         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
915         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
916         ! - ML - Option not developed yet
917         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
918         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
919      ENDIF
920
921   END SUBROUTINE dom_vvl_ctl
922
923   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
924      !!---------------------------------------------------------------------
925      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
926      !!                     
927      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
928      !!                scale factors at locations that have been individually
929      !!                modified in domhgr. Such modifications break the
930      !!                relationship between e12t and e1u*e2u etc.
931      !!                Recompute some scale factors ignoring the modified metric.
932      !!----------------------------------------------------------------------
933      !! * Arguments
934      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
935      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
936      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
937      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
938      !! * Local declarations
939      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
940      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
941      !! acc
942      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
943      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
944      !!
945      !                                                ! =====================
946      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
947         !                                             ! =====================
948      !! acc
949         IF( nn_cla == 0 ) THEN
950            !
951            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
952            ij0 = 102   ;   ij1 = 102
953            DO jk = 1, jpkm1
954               DO jj = mj0(ij0), mj1(ij1)
955                  DO ji = mi0(ii0), mi1(ii1)
956                     SELECT CASE ( pout )
957                     CASE( 'U' )
958                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
959                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
960                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
961                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
962                     CASE( 'F' )
963                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
964                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
965                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
966                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
967                     END SELECT
968                  END DO
969               END DO
970            END DO
971            !
972            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
973            ij0 =  88   ;   ij1 =  88
974            DO jk = 1, jpkm1
975               DO jj = mj0(ij0), mj1(ij1)
976                  DO ji = mi0(ii0), mi1(ii1)
977                     SELECT CASE ( pout )
978                     CASE( 'U' )
979                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
980                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
981                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
982                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
983                     CASE( 'V' )
984                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
985                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
986                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
987                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
988                     CASE( 'F' )
989                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
990                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
991                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
992                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
993                     END SELECT
994                  END DO
995               END DO
996            END DO
997         ENDIF
998
999         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1000         ij0 = 116   ;   ij1 = 116
1001         DO jk = 1, jpkm1
1002            DO jj = mj0(ij0), mj1(ij1)
1003               DO ji = mi0(ii0), mi1(ii1)
1004                  SELECT CASE ( pout )
1005                  CASE( 'U' )
1006                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1007                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1008                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1009                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1010                  CASE( 'F' )
1011                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1012                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1013                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1014                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1015                  END SELECT
1016               END DO
1017            END DO
1018         END DO
1019      ENDIF
1020      !
1021         !                                             ! =====================
1022      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1023         !                                             ! =====================
1024         !
1025         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1026         ij0 = 200   ;   ij1 = 200
1027         DO jk = 1, jpkm1
1028            DO jj = mj0(ij0), mj1(ij1)
1029               DO ji = mi0(ii0), mi1(ii1)
1030                  SELECT CASE ( pout )
1031                  CASE( 'U' )
1032                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1033                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1034                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1035                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1036                  CASE( 'F' )
1037                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1038                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1039                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1040                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1041                  END SELECT
1042               END DO
1043            END DO
1044         END DO
1045         !
1046         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1047         ij0 = 208   ;   ij1 = 208
1048         DO jk = 1, jpkm1
1049            DO jj = mj0(ij0), mj1(ij1)
1050               DO ji = mi0(ii0), mi1(ii1)
1051                  SELECT CASE ( pout )
1052                  CASE( 'U' )
1053                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1054                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1055                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1056                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1057                  CASE( 'F' )
1058                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1059                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1060                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1061                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1062                  END SELECT
1063               END DO
1064            END DO
1065         END DO
1066         !
1067         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1068         ij0 = 124   ;   ij1 = 125
1069         DO jk = 1, jpkm1
1070            DO jj = mj0(ij0), mj1(ij1)
1071               DO ji = mi0(ii0), mi1(ii1)
1072                  SELECT CASE ( pout )
1073                  CASE( 'V' )
1074                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1075                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1076                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1077                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1078                  END SELECT
1079               END DO
1080            END DO
1081         END DO
1082         !
1083         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1084         ij0 = 124   ;   ij1 = 125
1085         DO jk = 1, jpkm1
1086            DO jj = mj0(ij0), mj1(ij1)
1087               DO ji = mi0(ii0), mi1(ii1)
1088                  SELECT CASE ( pout )
1089                  CASE( 'V' )
1090                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1091                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1092                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1093                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1094                  END SELECT
1095               END DO
1096            END DO
1097         END DO
1098         !
1099         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1100         ij0 = 124   ;   ij1 = 125
1101         DO jk = 1, jpkm1
1102            DO jj = mj0(ij0), mj1(ij1)
1103               DO ji = mi0(ii0), mi1(ii1)
1104                  SELECT CASE ( pout )
1105                  CASE( 'V' )
1106                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1107                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1108                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1109                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1110                  END SELECT
1111               END DO
1112            END DO
1113         END DO
1114         !
1115         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1116         ij0 = 124   ;   ij1 = 125
1117         DO jk = 1, jpkm1
1118            DO jj = mj0(ij0), mj1(ij1)
1119               DO ji = mi0(ii0), mi1(ii1)
1120                  SELECT CASE ( pout )
1121                  CASE( 'V' )
1122                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1123                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1124                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1125                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1126                  END SELECT
1127               END DO
1128            END DO
1129         END DO
1130         !
1131         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1132         ij0 = 141   ;   ij1 = 142
1133         DO jk = 1, jpkm1
1134            DO jj = mj0(ij0), mj1(ij1)
1135               DO ji = mi0(ii0), mi1(ii1)
1136                  SELECT CASE ( pout )
1137                  CASE( 'V' )
1138                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1139                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1140                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1141                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1142                  END SELECT
1143               END DO
1144            END DO
1145         END DO
1146         !
1147         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1148         ij0 = 141   ;   ij1 = 142
1149         DO jk = 1, jpkm1
1150            DO jj = mj0(ij0), mj1(ij1)
1151               DO ji = mi0(ii0), mi1(ii1)
1152                  SELECT CASE ( pout )
1153                  CASE( 'V' )
1154                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1155                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1156                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1157                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1158                  END SELECT
1159               END DO
1160            END DO
1161         END DO
1162      ENDIF
1163         !                                             ! =====================
1164      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1165         !                                             ! =====================
1166         !
1167         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1168         ij0 = 327   ;   ij1 = 327
1169         DO jk = 1, jpkm1
1170            DO jj = mj0(ij0), mj1(ij1)
1171               DO ji = mi0(ii0), mi1(ii1)
1172                  SELECT CASE ( pout )
1173                  CASE( 'U' )
1174                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1175                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1176                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1177                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1178                  CASE( 'F' )
1179                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1180                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1181                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1182                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1183                  END SELECT
1184               END DO
1185            END DO
1186         END DO
1187         !
1188         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1189         ij0 = 343   ;   ij1 = 343
1190         DO jk = 1, jpkm1
1191            DO jj = mj0(ij0), mj1(ij1)
1192               DO ji = mi0(ii0), mi1(ii1)
1193                  SELECT CASE ( pout )
1194                  CASE( 'U' )
1195                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1196                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1197                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1198                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1199                  CASE( 'F' )
1200                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1201                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1202                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1203                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1204                  END SELECT
1205               END DO
1206            END DO
1207         END DO
1208         !
1209         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1210         ij0 = 232   ;   ij1 = 232
1211         DO jk = 1, jpkm1
1212            DO jj = mj0(ij0), mj1(ij1)
1213               DO ji = mi0(ii0), mi1(ii1)
1214                  SELECT CASE ( pout )
1215                  CASE( 'U' )
1216                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1217                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1218                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1219                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1220                  CASE( 'F' )
1221                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1222                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1223                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1224                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1225                  END SELECT
1226               END DO
1227            END DO
1228         END DO
1229         !
1230         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1231         ij0 = 232   ;   ij1 = 232
1232         DO jk = 1, jpkm1
1233            DO jj = mj0(ij0), mj1(ij1)
1234               DO ji = mi0(ii0), mi1(ii1)
1235                  SELECT CASE ( pout )
1236                  CASE( 'U' )
1237                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1238                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1239                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1240                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1241                  CASE( 'F' )
1242                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1243                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1244                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1245                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1246                  END SELECT
1247               END DO
1248            END DO
1249         END DO
1250         !
1251         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1252         ij0 = 270   ;   ij1 = 270
1253         DO jk = 1, jpkm1
1254            DO jj = mj0(ij0), mj1(ij1)
1255               DO ji = mi0(ii0), mi1(ii1)
1256                  SELECT CASE ( pout )
1257                  CASE( 'U' )
1258                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1259                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1260                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1261                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1262                  CASE( 'F' )
1263                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1264                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1265                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1266                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1267                  END SELECT
1268               END DO
1269            END DO
1270         END DO
1271         !
1272         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1273         ij0 = 232   ;   ij1 = 233
1274         DO jk = 1, jpkm1
1275            DO jj = mj0(ij0), mj1(ij1)
1276               DO ji = mi0(ii0), mi1(ii1)
1277                  SELECT CASE ( pout )
1278                  CASE( 'V' )
1279                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1280                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1281                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1282                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1283                  END SELECT
1284               END DO
1285            END DO
1286         END DO
1287         !
1288         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1289         ij0 = 276   ;   ij1 = 276
1290         DO jk = 1, jpkm1
1291            DO jj = mj0(ij0), mj1(ij1)
1292               DO ji = mi0(ii0), mi1(ii1)
1293                  SELECT CASE ( pout )
1294                  CASE( 'V' )
1295                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1296                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1297                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1298                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1299                  END SELECT
1300               END DO
1301            END DO
1302         END DO
1303      ENDIF
1304   END SUBROUTINE dom_vvl_orca_fix
1305
1306   !!======================================================================
1307END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.