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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4724

Last change on this file since 4724 was 4724, checked in by mathiot, 10 years ago

ISF branch: add comments, fix mpp and restar issues, add test to stop if incompatible options and fix mask issue in sbcice and sbcblk.

  • Property svn:keywords set to Id
File size: 74.2 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
15   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
16   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
17   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
18   !!   dom_vvl_rst      : read/write restart file
19   !!   dom_vvl_ctl      : Check the vvl options
20   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
21   !!                    : to account for manual changes to e[1,2][u,v] in some Straits
22   !!----------------------------------------------------------------------
23   !! * Modules used
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! ocean surface boundary condition
27   USE in_out_manager  ! I/O manager
28   USE iom             ! I/O manager library
29   USE restart         ! ocean restart
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC  dom_vvl_init       ! called by domain.F90
40   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
41   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
42   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
43   PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol
44
45   !!* Namelist nam_vvl
46   LOGICAL , PUBLIC                                      :: ln_vvl_zstar              ! zstar  vertical coordinate
47   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde             ! ztilde vertical coordinate
48   LOGICAL , PUBLIC                                      :: ln_vvl_layer              ! level  vertical coordinate
49   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor     ! ztilde vertical coordinate
51   LOGICAL , PUBLIC                                      :: ln_vvl_kepe               ! kinetic/potential energy transfer
52   !                                                                                           ! conservation: not used yet
53   REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient
54   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
55   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
56   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
57   LOGICAL , PUBLIC                                      :: ln_vvl_dbg                ! debug control prints
58
59   !! * Module variables
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
62   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
64   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
65   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
66
67   !! * Substitutions
68#  include "domzgr_substitute.h90"
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75
76CONTAINS
77
78   INTEGER FUNCTION dom_vvl_alloc()
79      !!----------------------------------------------------------------------
80      !!                ***  FUNCTION dom_vvl_alloc  ***
81      !!----------------------------------------------------------------------
82      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
83      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
84         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
85            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
86            &      STAT = dom_vvl_alloc        )
87         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
88         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
89         un_td = 0.0_wp
90         vn_td = 0.0_wp
91      ENDIF
92      IF( ln_vvl_ztilde ) THEN
93         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
94         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
95         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
96      ENDIF
97
98   END FUNCTION dom_vvl_alloc
99
100
101   SUBROUTINE dom_vvl_init
102      !!----------------------------------------------------------------------
103      !!                ***  ROUTINE dom_vvl_init  ***
104      !!                   
105      !! ** Purpose :  Initialization of all scale factors, depths
106      !!               and water column heights
107      !!
108      !! ** Method  :  - use restart file and/or initialize
109      !!               - interpolate scale factors
110      !!
111      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
112      !!              - Regrid: fse3(u/v)_n
113      !!                        fse3(u/v)_b       
114      !!                        fse3w_n           
115      !!                        fse3(u/v)w_b     
116      !!                        fse3(u/v)w_n     
117      !!                        fsdept_n, fsdepw_n and fsde3w_n
118      !!              - h(t/u/v)_0
119      !!              - frq_rst_e3t and frq_rst_hdv
120      !!
121      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
122      !!----------------------------------------------------------------------
123      USE phycst,  ONLY : rpi, rsmall, rad
124      !! * Local declarations
125      INTEGER ::   ji,jj,jk
126      INTEGER ::   ii0, ii1, ij0, ij1
127      !!----------------------------------------------------------------------
128      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
129
130      IF(lwp) WRITE(numout,*)
131      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
132      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
133
134      ! choose vertical coordinate (z_star, z_tilde or layer)
135      ! ==========================
136      CALL dom_vvl_ctl
137
138      ! Allocate module arrays
139      ! ======================
140      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
141
142      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
143      ! ============================================================================
144      CALL dom_vvl_rst( nit000, 'READ' )
145      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
146
147      ! Reconstruction of all vertical scale factors at now and before time steps
148      ! =============================================================================
149      ! Horizontal scale factor interpolations
150      ! --------------------------------------
151      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
152      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
153      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
154      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
155      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
156      ! Vertical scale factor interpolations
157      ! ------------------------------------
158      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
159      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
160      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
161      CALL dom_vvl_interpol( 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(:,:) ) * lmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - lmask(:,:) )
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.)
408         CALL lbc_lnk( vn_td , 'V' , -1.)
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. )
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. - lmask )
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. )
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. )
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. )
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, 'hdif_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               IF( neuler == 0 ) THEN
852                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
853               ENDIF
854            ELSE IF( id1 > 0 ) THEN
855               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
856               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
857               IF(lwp) write(numout,*) 'neuler is forced to 0'
858               fse3t_b(:,:,:) = fse3t_n(:,:,:)
859               neuler = 0
860            ELSE
861               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
862               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
863               IF(lwp) write(numout,*) 'neuler is forced to 0'
864               DO jk=1,jpk
865                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
866                      &                            / ( ht_0(:,:) + 1._wp - tmask_i(:,:) ) * tmask(:,:,jk) &
867                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
868               END DO
869               fse3t_b(:,:,:) = fse3t_n(:,:,:)
870               neuler = 0
871            ENDIF
872            !                             ! ----------- !
873            IF( ln_vvl_zstar ) THEN       ! z_star case !
874               !                          ! ----------- !
875               IF( MIN( id3, id4 ) > 0 ) THEN
876                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
877               ENDIF
878               !                          ! ----------------------- !
879            ELSE                          ! z_tilde and layer cases !
880               !                          ! ----------------------- !
881               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
882                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
883                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
884               ELSE                            ! one at least array is missing
885                  tilde_e3t_b(:,:,:) = 0.0_wp
886                  tilde_e3t_n(:,:,:) = 0.0_wp
887               ENDIF
888               !                          ! ------------ !
889               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
890                  !                       ! ------------ !
891                  IF( id5 > 0 ) THEN  ! required array exists
892                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
893                  ELSE                ! array is missing
894                     hdiv_lf(:,:,:) = 0.0_wp
895                  ENDIF
896               ENDIF
897            ENDIF
898            !
899         ELSE                                   !* Initialize at "rest"
900            fse3t_b(:,:,:) = e3t_0(:,:,:)
901            fse3t_n(:,:,:) = e3t_0(:,:,:)
902            sshn(:,:) = 0.0_wp
903            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
904               tilde_e3t_b(:,:,:) = 0.0_wp
905               tilde_e3t_n(:,:,:) = 0.0_wp
906               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
907            END IF
908         ENDIF
909
910      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
911         !                                   ! ===================
912         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
913         !                                           ! --------- !
914         !                                           ! all cases !
915         !                                           ! --------- !
916         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
917         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
918         !                                           ! ----------------------- !
919         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
920            !                                        ! ----------------------- !
921            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
922            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
923         END IF
924         !                                           ! -------------!   
925         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
926            !                                        ! ------------ !
927            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
928         ENDIF
929
930      ENDIF
931      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
932
933   END SUBROUTINE dom_vvl_rst
934
935
936   SUBROUTINE dom_vvl_ctl
937      !!---------------------------------------------------------------------
938      !!                  ***  ROUTINE dom_vvl_ctl  ***
939      !!               
940      !! ** Purpose :   Control the consistency between namelist options
941      !!                for vertical coordinate
942      !!----------------------------------------------------------------------
943      INTEGER ::   ioptio
944      INTEGER ::   ios
945
946      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
947                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
948                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
949      !!----------------------------------------------------------------------
950
951      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
952      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
953901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
954
955      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
956      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
957902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
958      IF(lwm) WRITE ( numond, nam_vvl )
959
960      IF(lwp) THEN                    ! Namelist print
961         WRITE(numout,*)
962         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
963         WRITE(numout,*) '~~~~~~~~~~~'
964         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
965         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
966         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
967         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
968         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
969         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
970         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
971         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
972         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
973         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
974         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
975         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
976         IF( ln_vvl_ztilde_as_zstar ) THEN
977            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
978            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
979            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
980            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
981            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
982            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
983         ELSE
984            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
985            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
986            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
987            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
988         ENDIF
989         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
990         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
991      ENDIF
992
993      ioptio = 0                      ! Parameter control
994      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
995      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
996      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
997      IF( ln_vvl_layer           )        ioptio = ioptio + 1
998
999      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1000      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0)   CALL ctl_stop( 'vvl_ztilde, vvl_layer, vvl_ztilde_as_zstar, vvl_zstar_at_eqtor not tested with ice shelf cavity (only vvl_zstar was tested)' )
1001
1002      IF(lwp) THEN                   ! Print the choice
1003         WRITE(numout,*)
1004         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1005         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1006         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1007         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1008         ! - ML - Option not developed yet
1009         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1010         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1011      ENDIF
1012
1013#if defined key_agrif
1014      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1015#endif
1016
1017   END SUBROUTINE dom_vvl_ctl
1018
1019   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1020      !!---------------------------------------------------------------------
1021      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1022      !!                     
1023      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1024      !!                scale factors at locations that have been individually
1025      !!                modified in domhgr. Such modifications break the
1026      !!                relationship between e12t and e1u*e2u etc.
1027      !!                Recompute some scale factors ignoring the modified metric.
1028      !!----------------------------------------------------------------------
1029      !! * Arguments
1030      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1031      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1032      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1033      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1034      !! * Local declarations
1035      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1036      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1037      !! acc
1038      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1039      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1040      !!
1041      !                                                ! =====================
1042      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1043         !                                             ! =====================
1044      !! acc
1045         IF( nn_cla == 0 ) THEN
1046            !
1047            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1048            ij0 = 102   ;   ij1 = 102
1049            DO jk = 1, jpkm1
1050               DO jj = mj0(ij0), mj1(ij1)
1051                  DO ji = mi0(ii0), mi1(ii1)
1052                     SELECT CASE ( pout )
1053                     CASE( 'U' )
1054                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1055                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1056                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1057                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1058                     CASE( 'F' )
1059                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1060                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1061                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1062                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1063                     END SELECT
1064                  END DO
1065               END DO
1066            END DO
1067            !
1068            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1069            ij0 =  88   ;   ij1 =  88
1070            DO jk = 1, jpkm1
1071               DO jj = mj0(ij0), mj1(ij1)
1072                  DO ji = mi0(ii0), mi1(ii1)
1073                     SELECT CASE ( pout )
1074                     CASE( 'U' )
1075                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1076                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1077                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1078                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1079                     CASE( 'V' )
1080                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1081                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1082                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1083                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1084                     CASE( 'F' )
1085                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1086                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1087                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1088                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1089                     END SELECT
1090                  END DO
1091               END DO
1092            END DO
1093         ENDIF
1094
1095         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1096         ij0 = 116   ;   ij1 = 116
1097         DO jk = 1, jpkm1
1098            DO jj = mj0(ij0), mj1(ij1)
1099               DO ji = mi0(ii0), mi1(ii1)
1100                  SELECT CASE ( pout )
1101                  CASE( 'U' )
1102                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1103                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1104                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1105                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1106                  CASE( 'F' )
1107                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1108                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1109                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1110                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1111                  END SELECT
1112               END DO
1113            END DO
1114         END DO
1115      ENDIF
1116      !
1117         !                                             ! =====================
1118      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1119         !                                             ! =====================
1120         !
1121         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1122         ij0 = 200   ;   ij1 = 200
1123         DO jk = 1, jpkm1
1124            DO jj = mj0(ij0), mj1(ij1)
1125               DO ji = mi0(ii0), mi1(ii1)
1126                  SELECT CASE ( pout )
1127                  CASE( 'U' )
1128                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1129                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1130                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1131                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1132                  CASE( 'F' )
1133                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1134                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1135                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1136                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1137                  END SELECT
1138               END DO
1139            END DO
1140         END DO
1141         !
1142         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1143         ij0 = 208   ;   ij1 = 208
1144         DO jk = 1, jpkm1
1145            DO jj = mj0(ij0), mj1(ij1)
1146               DO ji = mi0(ii0), mi1(ii1)
1147                  SELECT CASE ( pout )
1148                  CASE( 'U' )
1149                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1150                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1151                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1152                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1153                  CASE( 'F' )
1154                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1155                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1156                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1157                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1158                  END SELECT
1159               END DO
1160            END DO
1161         END DO
1162         !
1163         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1164         ij0 = 124   ;   ij1 = 125
1165         DO jk = 1, jpkm1
1166            DO jj = mj0(ij0), mj1(ij1)
1167               DO ji = mi0(ii0), mi1(ii1)
1168                  SELECT CASE ( pout )
1169                  CASE( 'V' )
1170                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1171                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1172                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1173                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1174                  END SELECT
1175               END DO
1176            END DO
1177         END DO
1178         !
1179         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1180         ij0 = 124   ;   ij1 = 125
1181         DO jk = 1, jpkm1
1182            DO jj = mj0(ij0), mj1(ij1)
1183               DO ji = mi0(ii0), mi1(ii1)
1184                  SELECT CASE ( pout )
1185                  CASE( 'V' )
1186                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1187                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1188                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1189                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1190                  END SELECT
1191               END DO
1192            END DO
1193         END DO
1194         !
1195         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1196         ij0 = 124   ;   ij1 = 125
1197         DO jk = 1, jpkm1
1198            DO jj = mj0(ij0), mj1(ij1)
1199               DO ji = mi0(ii0), mi1(ii1)
1200                  SELECT CASE ( pout )
1201                  CASE( 'V' )
1202                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1203                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1204                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1205                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1206                  END SELECT
1207               END DO
1208            END DO
1209         END DO
1210         !
1211         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1212         ij0 = 124   ;   ij1 = 125
1213         DO jk = 1, jpkm1
1214            DO jj = mj0(ij0), mj1(ij1)
1215               DO ji = mi0(ii0), mi1(ii1)
1216                  SELECT CASE ( pout )
1217                  CASE( 'V' )
1218                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1219                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1220                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1221                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1222                  END SELECT
1223               END DO
1224            END DO
1225         END DO
1226         !
1227         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1228         ij0 = 141   ;   ij1 = 142
1229         DO jk = 1, jpkm1
1230            DO jj = mj0(ij0), mj1(ij1)
1231               DO ji = mi0(ii0), mi1(ii1)
1232                  SELECT CASE ( pout )
1233                  CASE( 'V' )
1234                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1235                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1236                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1237                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1238                  END SELECT
1239               END DO
1240            END DO
1241         END DO
1242         !
1243         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1244         ij0 = 141   ;   ij1 = 142
1245         DO jk = 1, jpkm1
1246            DO jj = mj0(ij0), mj1(ij1)
1247               DO ji = mi0(ii0), mi1(ii1)
1248                  SELECT CASE ( pout )
1249                  CASE( 'V' )
1250                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1251                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1252                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1253                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1254                  END SELECT
1255               END DO
1256            END DO
1257         END DO
1258      ENDIF
1259         !                                             ! =====================
1260      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1261         !                                             ! =====================
1262         !
1263         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1264         ij0 = 327   ;   ij1 = 327
1265         DO jk = 1, jpkm1
1266            DO jj = mj0(ij0), mj1(ij1)
1267               DO ji = mi0(ii0), mi1(ii1)
1268                  SELECT CASE ( pout )
1269                  CASE( 'U' )
1270                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1271                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1272                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1273                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1274                  CASE( 'F' )
1275                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1276                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1277                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1278                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1279                  END SELECT
1280               END DO
1281            END DO
1282         END DO
1283         !
1284         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1285         ij0 = 343   ;   ij1 = 343
1286         DO jk = 1, jpkm1
1287            DO jj = mj0(ij0), mj1(ij1)
1288               DO ji = mi0(ii0), mi1(ii1)
1289                  SELECT CASE ( pout )
1290                  CASE( 'U' )
1291                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1292                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1293                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1294                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1295                  CASE( 'F' )
1296                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1297                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1298                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1299                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1300                  END SELECT
1301               END DO
1302            END DO
1303         END DO
1304         !
1305         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1306         ij0 = 232   ;   ij1 = 232
1307         DO jk = 1, jpkm1
1308            DO jj = mj0(ij0), mj1(ij1)
1309               DO ji = mi0(ii0), mi1(ii1)
1310                  SELECT CASE ( pout )
1311                  CASE( 'U' )
1312                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1313                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1314                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1315                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1316                  CASE( 'F' )
1317                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1318                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1319                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1320                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1321                  END SELECT
1322               END DO
1323            END DO
1324         END DO
1325         !
1326         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1327         ij0 = 232   ;   ij1 = 232
1328         DO jk = 1, jpkm1
1329            DO jj = mj0(ij0), mj1(ij1)
1330               DO ji = mi0(ii0), mi1(ii1)
1331                  SELECT CASE ( pout )
1332                  CASE( 'U' )
1333                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1334                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1335                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1336                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1337                  CASE( 'F' )
1338                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1339                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1340                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1341                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1342                  END SELECT
1343               END DO
1344            END DO
1345         END DO
1346         !
1347         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1348         ij0 = 270   ;   ij1 = 270
1349         DO jk = 1, jpkm1
1350            DO jj = mj0(ij0), mj1(ij1)
1351               DO ji = mi0(ii0), mi1(ii1)
1352                  SELECT CASE ( pout )
1353                  CASE( 'U' )
1354                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1355                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1356                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1357                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1358                  CASE( 'F' )
1359                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1360                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1361                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1362                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1363                  END SELECT
1364               END DO
1365            END DO
1366         END DO
1367         !
1368         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1369         ij0 = 232   ;   ij1 = 233
1370         DO jk = 1, jpkm1
1371            DO jj = mj0(ij0), mj1(ij1)
1372               DO ji = mi0(ii0), mi1(ii1)
1373                  SELECT CASE ( pout )
1374                  CASE( 'V' )
1375                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1376                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1377                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1378                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1379                  END SELECT
1380               END DO
1381            END DO
1382         END DO
1383         !
1384         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1385         ij0 = 276   ;   ij1 = 276
1386         DO jk = 1, jpkm1
1387            DO jj = mj0(ij0), mj1(ij1)
1388               DO ji = mi0(ii0), mi1(ii1)
1389                  SELECT CASE ( pout )
1390                  CASE( 'V' )
1391                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1392                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1393                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1394                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1395                  END SELECT
1396               END DO
1397            END DO
1398         END DO
1399      ENDIF
1400   END SUBROUTINE dom_vvl_orca_fix
1401
1402   !!======================================================================
1403END MODULE domvvl
1404
1405
1406
Note: See TracBrowser for help on using the repository browser.