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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4286

Last change on this file since 4286 was 4286, checked in by davestorkey, 10 years ago

Fix to allow VVL runs to read old VVL restarts.

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