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/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago

add wmask, wumask, wvmask and restore loop order and add flag to ignore specific isf code if no isf
  • Property svn:keywords set to Id
File size: 75.0 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      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
586      INTEGER                             :: ji,jj,jk       ! dummy loop indices
587      REAL(wp)                            :: zcoef
588      !!----------------------------------------------------------------------
589
590      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
591      !
592      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
593      !
594      IF( kt == nit000 )   THEN
595         IF(lwp) WRITE(numout,*)
596         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
597         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
598      ENDIF
599      !
600      ! Time filter and swap of scale factors
601      ! =====================================
602      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
603      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
604         IF( neuler == 0 .AND. kt == nit000 ) THEN
605            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
606         ELSE
607            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
608            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
609         ENDIF
610         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
611      ENDIF
612      fsdept_b(:,:,:) = fsdept_n(:,:,:)
613      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
614
615      fse3t_n(:,:,:) = fse3t_a(:,:,:)
616      fse3u_n(:,:,:) = fse3u_a(:,:,:)
617      fse3v_n(:,:,:) = fse3v_a(:,:,:)
618
619      ! Compute all missing vertical scale factor and depths
620      ! ====================================================
621      ! Horizontal scale factor interpolations
622      ! --------------------------------------
623      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
624      ! - JC - hu_b, hv_b, hur_b, hvr_b also
625      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
626      ! Vertical scale factor interpolations
627      ! ------------------------------------
628      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
629      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
630      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
631      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
632      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
633      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
634      ! t- and w- points depth
635      ! ----------------------
636      ! set the isf depth as it is in the initial step
637      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
638      fsdepw_n(:,:,1) = 0.0_wp
639      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
640
641      DO jk = 2, jpk
642         DO jj = 1,jpj
643            DO ji = 1,jpi
644              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
645                                                                 ! 1 for jk = mikt
646               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
647               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
648               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
649                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk)) 
650               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj)
651            END DO
652         END DO
653      END DO
654
655      ! Local depth and Inverse of the local depth of the water column at u- and v- points
656      ! ----------------------------------------------------------------------------------
657      hu (:,:) = hu_a (:,:)
658      hv (:,:) = hv_a (:,:)
659
660      ! Inverse of the local depth
661      hur(:,:) = hur_a(:,:)
662      hvr(:,:) = hvr_a(:,:)
663
664      ! Local depth of the water column at t- points
665      ! --------------------------------------------
666      ht(:,:) = 0.
667      DO jk = 1, jpkm1
668         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
669      END DO
670
671      ! Write outputs
672      ! =============
673      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
674      CALL iom_put( "cellthc" , fse3t_n  (:,:,:) )
675      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
676      CALL iom_put( "e3tdef"  , z_e3t_def(:,:,:) )
677
678      ! write restart file
679      ! ==================
680      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
681      !
682      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
683      !
684      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
685
686   END SUBROUTINE dom_vvl_sf_swp
687
688
689   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
690      !!---------------------------------------------------------------------
691      !!                  ***  ROUTINE dom_vvl__interpol  ***
692      !!
693      !! ** Purpose :   interpolate scale factors from one grid point to another
694      !!
695      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
696      !!                - horizontal interpolation: grid cell surface averaging
697      !!                - vertical interpolation: simple averaging
698      !!----------------------------------------------------------------------
699      !! * Arguments
700      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
701      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
702      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
703      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
704      !! * Local declarations
705      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
706      LOGICAL ::   l_is_orca                                           ! local logical
707      !!----------------------------------------------------------------------
708      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
709         !
710      l_is_orca = .FALSE.
711      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
712
713      SELECT CASE ( pout )
714         !               ! ------------------------------------- !
715      CASE( 'U' )        ! interpolation from T-point to U-point !
716         !               ! ------------------------------------- !
717         ! horizontal surface weighted interpolation
718         DO jk = 1, jpk
719            DO jj = 1, jpjm1
720               DO ji = 1, fs_jpim1   ! vector opt.
721                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
722                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
723                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
724               END DO
725            END DO
726         END DO
727         !
728         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
729         ! boundary conditions
730         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
731         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
732         !               ! ------------------------------------- !
733      CASE( 'V' )        ! interpolation from T-point to V-point !
734         !               ! ------------------------------------- !
735         ! horizontal surface weighted interpolation
736         DO jk = 1, jpk
737            DO jj = 1, jpjm1
738               DO ji = 1, fs_jpim1   ! vector opt.
739                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
740                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
741                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
742               END DO
743            END DO
744         END DO
745         !
746         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
747         ! boundary conditions
748         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
749         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
750         !               ! ------------------------------------- !
751      CASE( 'F' )        ! interpolation from U-point to F-point !
752         !               ! ------------------------------------- !
753         ! horizontal surface weighted interpolation
754         DO jk = 1, jpk
755            DO jj = 1, jpjm1
756               DO ji = 1, fs_jpim1   ! vector opt.
757                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
758                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
759                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
760               END DO
761            END DO
762         END DO
763         !
764         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
765         ! boundary conditions
766         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
767         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
768         !               ! ------------------------------------- !
769      CASE( 'W' )        ! interpolation from T-point to W-point !
770         !               ! ------------------------------------- !
771         ! vertical simple interpolation
772         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
773         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
774         DO jk = 2, jpk
775            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
776               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
777         END DO
778         !               ! -------------------------------------- !
779      CASE( 'UW' )       ! interpolation from U-point to UW-point !
780         !               ! -------------------------------------- !
781         ! vertical simple interpolation
782         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
783         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
784         DO jk = 2, jpk
785            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
786               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
787         END DO
788         !               ! -------------------------------------- !
789      CASE( 'VW' )       ! interpolation from V-point to VW-point !
790         !               ! -------------------------------------- !
791         ! vertical simple interpolation
792         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
793         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
794         DO jk = 2, jpk
795            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
796               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
797         END DO
798      END SELECT
799      !
800
801      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
802
803   END SUBROUTINE dom_vvl_interpol
804
805   SUBROUTINE dom_vvl_rst( kt, cdrw )
806      !!---------------------------------------------------------------------
807      !!                   ***  ROUTINE dom_vvl_rst  ***
808      !!                     
809      !! ** Purpose :   Read or write VVL file in restart file
810      !!
811      !! ** Method  :   use of IOM library
812      !!                if the restart does not contain vertical scale factors,
813      !!                they are set to the _0 values
814      !!                if the restart does not contain vertical scale factors increments (z_tilde),
815      !!                they are set to 0.
816      !!----------------------------------------------------------------------
817      !! * Arguments
818      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
819      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
820      !! * Local declarations
821      INTEGER ::   jk
822      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
823      !!----------------------------------------------------------------------
824      !
825      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
826      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
827         !                                   ! ===============
828         IF( ln_rstart ) THEN                   !* Read the restart file
829            CALL rst_read_open                  !  open the restart file if necessary
830            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
831            !
832            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
833            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
834            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
835            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
836            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
837            !                             ! --------- !
838            !                             ! all cases !
839            !                             ! --------- !
840            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
841               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
842               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
843               ! needed to restart if land processor not computed
844               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
845               WHERE ( tmask(:,:,:) == 0.0_wp ) 
846                  fse3t_n(:,:,:) = e3t_0(:,:,:)
847                  fse3t_b(:,:,:) = e3t_0(:,:,:)
848               END WHERE
849               IF( neuler == 0 ) THEN
850                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
851               ENDIF
852            ELSE IF( id1 > 0 ) THEN
853               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
854               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
855               IF(lwp) write(numout,*) 'neuler is forced to 0'
856               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
857               fse3t_n(:,:,:) = fse3t_b(:,:,:)
858               neuler = 0
859            ELSE IF( id2 > 0 ) THEN
860               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
861               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
862               IF(lwp) write(numout,*) 'neuler is forced to 0'
863               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
864               fse3t_b(:,:,:) = fse3t_n(:,:,:)
865               neuler = 0
866            ELSE
867               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
868               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
869               IF(lwp) write(numout,*) 'neuler is forced to 0'
870               DO jk=1,jpk
871                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
872                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
873                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
874               END DO
875               fse3t_b(:,:,:) = fse3t_n(:,:,:)
876               neuler = 0
877            ENDIF
878            !                             ! ----------- !
879            IF( ln_vvl_zstar ) THEN       ! z_star case !
880               !                          ! ----------- !
881               IF( MIN( id3, id4 ) > 0 ) THEN
882                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
883               ENDIF
884               !                          ! ----------------------- !
885            ELSE                          ! z_tilde and layer cases !
886               !                          ! ----------------------- !
887               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
888                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
889                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
890               ELSE                            ! one at least array is missing
891                  tilde_e3t_b(:,:,:) = 0.0_wp
892                  tilde_e3t_n(:,:,:) = 0.0_wp
893               ENDIF
894               !                          ! ------------ !
895               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
896                  !                       ! ------------ !
897                  IF( id5 > 0 ) THEN  ! required array exists
898                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
899                  ELSE                ! array is missing
900                     hdiv_lf(:,:,:) = 0.0_wp
901                  ENDIF
902               ENDIF
903            ENDIF
904            !
905         ELSE                                   !* Initialize at "rest"
906            fse3t_b(:,:,:) = e3t_0(:,:,:)
907            fse3t_n(:,:,:) = e3t_0(:,:,:)
908            sshn(:,:) = 0.0_wp
909            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
910               tilde_e3t_b(:,:,:) = 0.0_wp
911               tilde_e3t_n(:,:,:) = 0.0_wp
912               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
913            END IF
914         ENDIF
915
916      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
917         !                                   ! ===================
918         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
919         !                                           ! --------- !
920         !                                           ! all cases !
921         !                                           ! --------- !
922         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
923         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
924         !                                           ! ----------------------- !
925         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
926            !                                        ! ----------------------- !
927            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
928            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
929         END IF
930         !                                           ! -------------!   
931         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
932            !                                        ! ------------ !
933            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
934         ENDIF
935
936      ENDIF
937      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
938
939   END SUBROUTINE dom_vvl_rst
940
941
942   SUBROUTINE dom_vvl_ctl
943      !!---------------------------------------------------------------------
944      !!                  ***  ROUTINE dom_vvl_ctl  ***
945      !!               
946      !! ** Purpose :   Control the consistency between namelist options
947      !!                for vertical coordinate
948      !!----------------------------------------------------------------------
949      INTEGER ::   ioptio
950      INTEGER ::   ios
951
952      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
953                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
954                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
955      !!----------------------------------------------------------------------
956
957      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
958      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
959901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
960
961      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
962      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
963902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
964      IF(lwm) WRITE ( numond, nam_vvl )
965
966      IF(lwp) THEN                    ! Namelist print
967         WRITE(numout,*)
968         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
969         WRITE(numout,*) '~~~~~~~~~~~'
970         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
971         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
972         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
973         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
974         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
975         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
976         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
977         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
978         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
979         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
980         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
981         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
982         IF( ln_vvl_ztilde_as_zstar ) THEN
983            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
984            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
985            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
986            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
987            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
988            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
989         ELSE
990            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
991            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
992            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
993            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
994         ENDIF
995         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
996         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
997      ENDIF
998
999      ioptio = 0                      ! Parameter control
1000      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1001      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
1002      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
1003      IF( ln_vvl_layer           )        ioptio = ioptio + 1
1004
1005      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1006      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1007
1008      IF(lwp) THEN                   ! Print the choice
1009         WRITE(numout,*)
1010         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1011         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1012         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1013         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1014         ! - ML - Option not developed yet
1015         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1016         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1017      ENDIF
1018
1019#if defined key_agrif
1020      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1021#endif
1022
1023   END SUBROUTINE dom_vvl_ctl
1024
1025   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1026      !!---------------------------------------------------------------------
1027      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1028      !!                     
1029      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1030      !!                scale factors at locations that have been individually
1031      !!                modified in domhgr. Such modifications break the
1032      !!                relationship between e12t and e1u*e2u etc.
1033      !!                Recompute some scale factors ignoring the modified metric.
1034      !!----------------------------------------------------------------------
1035      !! * Arguments
1036      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1037      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1038      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1039      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1040      !! * Local declarations
1041      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1042      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1043      !! acc
1044      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1045      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1046      !!
1047      !                                                ! =====================
1048      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1049         !                                             ! =====================
1050      !! acc
1051         IF( nn_cla == 0 ) THEN
1052            !
1053            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1054            ij0 = 102   ;   ij1 = 102
1055            DO jk = 1, jpkm1
1056               DO jj = mj0(ij0), mj1(ij1)
1057                  DO ji = mi0(ii0), mi1(ii1)
1058                     SELECT CASE ( pout )
1059                     CASE( 'U' )
1060                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1061                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1062                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1063                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1064                     CASE( 'F' )
1065                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1066                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1067                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1068                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1069                     END SELECT
1070                  END DO
1071               END DO
1072            END DO
1073            !
1074            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1075            ij0 =  88   ;   ij1 =  88
1076            DO jk = 1, jpkm1
1077               DO jj = mj0(ij0), mj1(ij1)
1078                  DO ji = mi0(ii0), mi1(ii1)
1079                     SELECT CASE ( pout )
1080                     CASE( 'U' )
1081                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1082                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1083                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1084                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1085                     CASE( 'V' )
1086                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1087                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1088                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1089                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1090                     CASE( 'F' )
1091                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1092                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1093                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1094                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1095                     END SELECT
1096                  END DO
1097               END DO
1098            END DO
1099         ENDIF
1100
1101         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1102         ij0 = 116   ;   ij1 = 116
1103         DO jk = 1, jpkm1
1104            DO jj = mj0(ij0), mj1(ij1)
1105               DO ji = mi0(ii0), mi1(ii1)
1106                  SELECT CASE ( pout )
1107                  CASE( 'U' )
1108                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1109                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1110                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1111                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1112                  CASE( 'F' )
1113                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1114                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1115                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1116                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1117                  END SELECT
1118               END DO
1119            END DO
1120         END DO
1121      ENDIF
1122      !
1123         !                                             ! =====================
1124      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1125         !                                             ! =====================
1126         !
1127         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1128         ij0 = 200   ;   ij1 = 200
1129         DO jk = 1, jpkm1
1130            DO jj = mj0(ij0), mj1(ij1)
1131               DO ji = mi0(ii0), mi1(ii1)
1132                  SELECT CASE ( pout )
1133                  CASE( 'U' )
1134                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1135                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1136                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1137                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1138                  CASE( 'F' )
1139                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1140                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1141                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1142                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1143                  END SELECT
1144               END DO
1145            END DO
1146         END DO
1147         !
1148         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1149         ij0 = 208   ;   ij1 = 208
1150         DO jk = 1, jpkm1
1151            DO jj = mj0(ij0), mj1(ij1)
1152               DO ji = mi0(ii0), mi1(ii1)
1153                  SELECT CASE ( pout )
1154                  CASE( 'U' )
1155                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1156                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1157                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1158                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1159                  CASE( 'F' )
1160                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1161                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1162                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1163                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1164                  END SELECT
1165               END DO
1166            END DO
1167         END DO
1168         !
1169         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1170         ij0 = 124   ;   ij1 = 125
1171         DO jk = 1, jpkm1
1172            DO jj = mj0(ij0), mj1(ij1)
1173               DO ji = mi0(ii0), mi1(ii1)
1174                  SELECT CASE ( pout )
1175                  CASE( 'V' )
1176                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1177                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1178                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1179                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1180                  END SELECT
1181               END DO
1182            END DO
1183         END DO
1184         !
1185         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1186         ij0 = 124   ;   ij1 = 125
1187         DO jk = 1, jpkm1
1188            DO jj = mj0(ij0), mj1(ij1)
1189               DO ji = mi0(ii0), mi1(ii1)
1190                  SELECT CASE ( pout )
1191                  CASE( 'V' )
1192                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1193                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1194                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1195                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1196                  END SELECT
1197               END DO
1198            END DO
1199         END DO
1200         !
1201         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1202         ij0 = 124   ;   ij1 = 125
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( 'V' )
1208                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1209                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1210                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1211                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1212                  END SELECT
1213               END DO
1214            END DO
1215         END DO
1216         !
1217         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1218         ij0 = 124   ;   ij1 = 125
1219         DO jk = 1, jpkm1
1220            DO jj = mj0(ij0), mj1(ij1)
1221               DO ji = mi0(ii0), mi1(ii1)
1222                  SELECT CASE ( pout )
1223                  CASE( 'V' )
1224                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1225                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1226                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1227                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1228                  END SELECT
1229               END DO
1230            END DO
1231         END DO
1232         !
1233         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1234         ij0 = 141   ;   ij1 = 142
1235         DO jk = 1, jpkm1
1236            DO jj = mj0(ij0), mj1(ij1)
1237               DO ji = mi0(ii0), mi1(ii1)
1238                  SELECT CASE ( pout )
1239                  CASE( 'V' )
1240                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1241                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1242                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1243                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1244                  END SELECT
1245               END DO
1246            END DO
1247         END DO
1248         !
1249         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1250         ij0 = 141   ;   ij1 = 142
1251         DO jk = 1, jpkm1
1252            DO jj = mj0(ij0), mj1(ij1)
1253               DO ji = mi0(ii0), mi1(ii1)
1254                  SELECT CASE ( pout )
1255                  CASE( 'V' )
1256                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1257                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1258                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1259                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1260                  END SELECT
1261               END DO
1262            END DO
1263         END DO
1264      ENDIF
1265         !                                             ! =====================
1266      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1267         !                                             ! =====================
1268         !
1269         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1270         ij0 = 327   ;   ij1 = 327
1271         DO jk = 1, jpkm1
1272            DO jj = mj0(ij0), mj1(ij1)
1273               DO ji = mi0(ii0), mi1(ii1)
1274                  SELECT CASE ( pout )
1275                  CASE( 'U' )
1276                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1277                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1278                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1279                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1280                  CASE( 'F' )
1281                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1282                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1283                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1284                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1285                  END SELECT
1286               END DO
1287            END DO
1288         END DO
1289         !
1290         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1291         ij0 = 343   ;   ij1 = 343
1292         DO jk = 1, jpkm1
1293            DO jj = mj0(ij0), mj1(ij1)
1294               DO ji = mi0(ii0), mi1(ii1)
1295                  SELECT CASE ( pout )
1296                  CASE( 'U' )
1297                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1298                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1299                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1300                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1301                  CASE( 'F' )
1302                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1303                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1304                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1305                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1306                  END SELECT
1307               END DO
1308            END DO
1309         END DO
1310         !
1311         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1312         ij0 = 232   ;   ij1 = 232
1313         DO jk = 1, jpkm1
1314            DO jj = mj0(ij0), mj1(ij1)
1315               DO ji = mi0(ii0), mi1(ii1)
1316                  SELECT CASE ( pout )
1317                  CASE( 'U' )
1318                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1319                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1320                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1321                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1322                  CASE( 'F' )
1323                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1324                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1325                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1326                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1327                  END SELECT
1328               END DO
1329            END DO
1330         END DO
1331         !
1332         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1333         ij0 = 232   ;   ij1 = 232
1334         DO jk = 1, jpkm1
1335            DO jj = mj0(ij0), mj1(ij1)
1336               DO ji = mi0(ii0), mi1(ii1)
1337                  SELECT CASE ( pout )
1338                  CASE( 'U' )
1339                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1340                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1341                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1342                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1343                  CASE( 'F' )
1344                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1345                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1346                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1347                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1348                  END SELECT
1349               END DO
1350            END DO
1351         END DO
1352         !
1353         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1354         ij0 = 270   ;   ij1 = 270
1355         DO jk = 1, jpkm1
1356            DO jj = mj0(ij0), mj1(ij1)
1357               DO ji = mi0(ii0), mi1(ii1)
1358                  SELECT CASE ( pout )
1359                  CASE( 'U' )
1360                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1361                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1362                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1363                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1364                  CASE( 'F' )
1365                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1366                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1367                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1368                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1369                  END SELECT
1370               END DO
1371            END DO
1372         END DO
1373         !
1374         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1375         ij0 = 232   ;   ij1 = 233
1376         DO jk = 1, jpkm1
1377            DO jj = mj0(ij0), mj1(ij1)
1378               DO ji = mi0(ii0), mi1(ii1)
1379                  SELECT CASE ( pout )
1380                  CASE( 'V' )
1381                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1382                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1383                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1384                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1385                  END SELECT
1386               END DO
1387            END DO
1388         END DO
1389         !
1390         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1391         ij0 = 276   ;   ij1 = 276
1392         DO jk = 1, jpkm1
1393            DO jj = mj0(ij0), mj1(ij1)
1394               DO ji = mi0(ii0), mi1(ii1)
1395                  SELECT CASE ( pout )
1396                  CASE( 'V' )
1397                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1398                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1399                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1400                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1401                  END SELECT
1402               END DO
1403            END DO
1404         END DO
1405      ENDIF
1406   END SUBROUTINE dom_vvl_orca_fix
1407
1408   !!======================================================================
1409END MODULE domvvl
1410
1411
1412
Note: See TracBrowser for help on using the repository browser.