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

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago

Dev. branch CMCC4_simplification ticket #1456

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