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

Last change on this file since 4338 was 4338, checked in by acc, 7 years ago

Branch dev_MERGE_2013. #1209. Fix to compute after scale factors before wn computation in time-splitting case. dom_vvl_sf_nxt is now called twice but correctly computes the baroclinic contribution only once. Still need to sort out the asselin filtering of the separate components

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