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

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

source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9252

Last change on this file since 9252 was 9252, checked in by andmirek, 6 years ago

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