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 @ 3953

Last change on this file since 3953 was 3951, checked in by acc, 11 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Minor fix to dynspg_ts.F90 and introduction of an option to suppress ztilde coordinate in the equatorial band

  • Property svn:keywords set to Id
File size: 69.0 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!----------------------------------------------------------------------
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                                 ! one at least array is missing
773               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
774            ENDIF
775            !                             ! ----------- !
776            IF( ln_vvl_zstar ) THEN       ! z_star case !
777               !                          ! ----------- !
778               IF( MIN( id3, id4 ) > 0 ) THEN
779                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
780               ENDIF
781               !                          ! ----------------------- !
782            ELSE                          ! z_tilde and layer cases !
783               !                          ! ----------------------- !
784               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
785                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
786                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
787               ELSE                            ! one at least array is missing
788                  tilde_e3t_b(:,:,:) = 0.0_wp
789                  tilde_e3t_n(:,:,:) = 0.0_wp
790               ENDIF
791               !                          ! ------------ !
792               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
793                  !                       ! ------------ !
794                  IF( id5 > 0 ) THEN  ! required array exists
795                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
796                  ELSE                ! array is missing
797                     hdiv_lf(:,:,:) = 0.0_wp
798                  ENDIF
799               ENDIF
800            ENDIF
801            !
802         ELSE                                   !* Initialize at "rest"
803            fse3t_b(:,:,:) = e3t_0(:,:,:)
804            fse3t_n(:,:,:) = e3t_0(:,:,:)
805            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
806               tilde_e3t_b(:,:,:) = 0.0_wp
807               tilde_e3t_n(:,:,:) = 0.0_wp
808               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
809            END IF
810         ENDIF
811
812      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
813         !                                   ! ===================
814         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
815         !                                           ! --------- !
816         !                                           ! all cases !
817         !                                           ! --------- !
818         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
819         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
820         !                                           ! ----------------------- !
821         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
822            !                                        ! ----------------------- !
823            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
824            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
825         END IF
826         !                                           ! -------------!   
827         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
828            !                                        ! ------------ !
829            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
830         ENDIF
831
832      ENDIF
833      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
834
835   END SUBROUTINE dom_vvl_rst
836
837
838   SUBROUTINE dom_vvl_ctl
839      !!---------------------------------------------------------------------
840      !!                  ***  ROUTINE dom_vvl_ctl  ***
841      !!               
842      !! ** Purpose :   Control the consistency between namelist options
843      !!                for vertical coordinate
844      !!----------------------------------------------------------------------
845      INTEGER ::   ioptio
846
847      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
848                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
849                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
850      !!----------------------------------------------------------------------
851
852      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
853      READ   ( numnam, nam_vvl )
854
855      IF(lwp) THEN                    ! Namelist print
856         WRITE(numout,*)
857         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
858         WRITE(numout,*) '~~~~~~~~~~~'
859         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
860         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
861         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
862         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
863         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
864         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
865         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
866         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
867         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
868         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
869         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
870         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
871         IF( ln_vvl_ztilde_as_zstar ) THEN
872            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
873            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
874            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
875            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
876            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
877            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
878         ELSE
879            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
880            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
881            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
882            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
883         ENDIF
884         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
885         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
886      ENDIF
887
888      ioptio = 0                      ! Parameter control
889      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
890      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
891      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
892      IF( ln_vvl_layer           )        ioptio = ioptio + 1
893
894      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
895
896      IF(lwp) THEN                   ! Print the choice
897         WRITE(numout,*)
898         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
899         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
900         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
901         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
902         ! - ML - Option not developed yet
903         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
904         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
905      ENDIF
906
907   END SUBROUTINE dom_vvl_ctl
908
909   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
910      !!---------------------------------------------------------------------
911      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
912      !!                     
913      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
914      !!                scale factors at locations that have been individually
915      !!                modified in domhgr. Such modifications break the
916      !!                relationship between e12t and e1u*e2u etc.
917      !!                Recompute some scale factors ignoring the modified metric.
918      !!----------------------------------------------------------------------
919      !! * Arguments
920      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
921      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
922      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
923      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
924      !! * Local declarations
925      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
926      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
927      !! acc
928      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
929      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
930      !!
931      !                                                ! =====================
932      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
933         !                                             ! =====================
934      !! acc
935         IF( nn_cla == 0 ) THEN
936            !
937            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
938            ij0 = 102   ;   ij1 = 102
939            DO jk = 1, jpkm1
940               DO jj = mj0(ij0), mj1(ij1)
941                  DO ji = mi0(ii0), mi1(ii1)
942                     SELECT CASE ( pout )
943                     CASE( 'U' )
944                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
945                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
946                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
947                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
948                     CASE( 'F' )
949                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
950                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
951                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
952                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
953                     END SELECT
954                  END DO
955               END DO
956            END DO
957            !
958            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
959            ij0 =  88   ;   ij1 =  88
960            DO jk = 1, jpkm1
961               DO jj = mj0(ij0), mj1(ij1)
962                  DO ji = mi0(ii0), mi1(ii1)
963                     SELECT CASE ( pout )
964                     CASE( 'U' )
965                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
966                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
967                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
968                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
969                     CASE( 'V' )
970                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
971                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
972                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
973                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
974                     CASE( 'F' )
975                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
976                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
977                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
978                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
979                     END SELECT
980                  END DO
981               END DO
982            END DO
983         ENDIF
984
985         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
986         ij0 = 116   ;   ij1 = 116
987         DO jk = 1, jpkm1
988            DO jj = mj0(ij0), mj1(ij1)
989               DO ji = mi0(ii0), mi1(ii1)
990                  SELECT CASE ( pout )
991                  CASE( 'U' )
992                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
993                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
994                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
995                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
996                  CASE( 'F' )
997                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
998                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
999                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1000                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1001                  END SELECT
1002               END DO
1003            END DO
1004         END DO
1005      ENDIF
1006      !
1007         !                                             ! =====================
1008      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1009         !                                             ! =====================
1010         !
1011         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1012         ij0 = 200   ;   ij1 = 200
1013         DO jk = 1, jpkm1
1014            DO jj = mj0(ij0), mj1(ij1)
1015               DO ji = mi0(ii0), mi1(ii1)
1016                  SELECT CASE ( pout )
1017                  CASE( 'U' )
1018                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1019                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1020                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1021                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1022                  CASE( 'F' )
1023                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1024                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1025                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1026                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1027                  END SELECT
1028               END DO
1029            END DO
1030         END DO
1031         !
1032         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1033         ij0 = 208   ;   ij1 = 208
1034         DO jk = 1, jpkm1
1035            DO jj = mj0(ij0), mj1(ij1)
1036               DO ji = mi0(ii0), mi1(ii1)
1037                  SELECT CASE ( pout )
1038                  CASE( 'U' )
1039                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1040                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1041                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1042                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1043                  CASE( 'F' )
1044                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1045                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1046                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1047                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1048                  END SELECT
1049               END DO
1050            END DO
1051         END DO
1052         !
1053         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1054         ij0 = 124   ;   ij1 = 125
1055         DO jk = 1, jpkm1
1056            DO jj = mj0(ij0), mj1(ij1)
1057               DO ji = mi0(ii0), mi1(ii1)
1058                  SELECT CASE ( pout )
1059                  CASE( 'V' )
1060                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1061                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1062                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1063                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1064                  END SELECT
1065               END DO
1066            END DO
1067         END DO
1068         !
1069         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1070         ij0 = 124   ;   ij1 = 125
1071         DO jk = 1, jpkm1
1072            DO jj = mj0(ij0), mj1(ij1)
1073               DO ji = mi0(ii0), mi1(ii1)
1074                  SELECT CASE ( pout )
1075                  CASE( 'V' )
1076                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1077                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1078                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1079                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1080                  END SELECT
1081               END DO
1082            END DO
1083         END DO
1084         !
1085         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1086         ij0 = 124   ;   ij1 = 125
1087         DO jk = 1, jpkm1
1088            DO jj = mj0(ij0), mj1(ij1)
1089               DO ji = mi0(ii0), mi1(ii1)
1090                  SELECT CASE ( pout )
1091                  CASE( 'V' )
1092                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1093                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1094                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1095                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1096                  END SELECT
1097               END DO
1098            END DO
1099         END DO
1100         !
1101         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1102         ij0 = 124   ;   ij1 = 125
1103         DO jk = 1, jpkm1
1104            DO jj = mj0(ij0), mj1(ij1)
1105               DO ji = mi0(ii0), mi1(ii1)
1106                  SELECT CASE ( pout )
1107                  CASE( 'V' )
1108                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1109                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1110                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1111                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1112                  END SELECT
1113               END DO
1114            END DO
1115         END DO
1116         !
1117         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1118         ij0 = 141   ;   ij1 = 142
1119         DO jk = 1, jpkm1
1120            DO jj = mj0(ij0), mj1(ij1)
1121               DO ji = mi0(ii0), mi1(ii1)
1122                  SELECT CASE ( pout )
1123                  CASE( 'V' )
1124                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1125                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1126                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1127                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1128                  END SELECT
1129               END DO
1130            END DO
1131         END DO
1132         !
1133         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1134         ij0 = 141   ;   ij1 = 142
1135         DO jk = 1, jpkm1
1136            DO jj = mj0(ij0), mj1(ij1)
1137               DO ji = mi0(ii0), mi1(ii1)
1138                  SELECT CASE ( pout )
1139                  CASE( 'V' )
1140                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1141                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1142                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1143                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1144                  END SELECT
1145               END DO
1146            END DO
1147         END DO
1148      ENDIF
1149         !                                             ! =====================
1150      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1151         !                                             ! =====================
1152         !
1153         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1154         ij0 = 327   ;   ij1 = 327
1155         DO jk = 1, jpkm1
1156            DO jj = mj0(ij0), mj1(ij1)
1157               DO ji = mi0(ii0), mi1(ii1)
1158                  SELECT CASE ( pout )
1159                  CASE( 'U' )
1160                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1161                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1162                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1163                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1164                  CASE( 'F' )
1165                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1166                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1167                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1168                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1169                  END SELECT
1170               END DO
1171            END DO
1172         END DO
1173         !
1174         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1175         ij0 = 343   ;   ij1 = 343
1176         DO jk = 1, jpkm1
1177            DO jj = mj0(ij0), mj1(ij1)
1178               DO ji = mi0(ii0), mi1(ii1)
1179                  SELECT CASE ( pout )
1180                  CASE( 'U' )
1181                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1182                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1183                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1184                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1185                  CASE( 'F' )
1186                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1187                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1188                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1189                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1190                  END SELECT
1191               END DO
1192            END DO
1193         END DO
1194         !
1195         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1196         ij0 = 232   ;   ij1 = 232
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( 'U' )
1202                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1203                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1204                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1205                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1206                  CASE( 'F' )
1207                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1208                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1209                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1210                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1211                  END SELECT
1212               END DO
1213            END DO
1214         END DO
1215         !
1216         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1217         ij0 = 232   ;   ij1 = 232
1218         DO jk = 1, jpkm1
1219            DO jj = mj0(ij0), mj1(ij1)
1220               DO ji = mi0(ii0), mi1(ii1)
1221                  SELECT CASE ( pout )
1222                  CASE( 'U' )
1223                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1224                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1225                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1226                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1227                  CASE( 'F' )
1228                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1229                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1230                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1231                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1232                  END SELECT
1233               END DO
1234            END DO
1235         END DO
1236         !
1237         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1238         ij0 = 270   ;   ij1 = 270
1239         DO jk = 1, jpkm1
1240            DO jj = mj0(ij0), mj1(ij1)
1241               DO ji = mi0(ii0), mi1(ii1)
1242                  SELECT CASE ( pout )
1243                  CASE( 'U' )
1244                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1245                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1246                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1247                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1248                  CASE( 'F' )
1249                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1250                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1251                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1252                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1253                  END SELECT
1254               END DO
1255            END DO
1256         END DO
1257         !
1258         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1259         ij0 = 232   ;   ij1 = 233
1260         DO jk = 1, jpkm1
1261            DO jj = mj0(ij0), mj1(ij1)
1262               DO ji = mi0(ii0), mi1(ii1)
1263                  SELECT CASE ( pout )
1264                  CASE( 'V' )
1265                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1266                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1267                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1268                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1269                  END SELECT
1270               END DO
1271            END DO
1272         END DO
1273         !
1274         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1275         ij0 = 276   ;   ij1 = 276
1276         DO jk = 1, jpkm1
1277            DO jj = mj0(ij0), mj1(ij1)
1278               DO ji = mi0(ii0), mi1(ii1)
1279                  SELECT CASE ( pout )
1280                  CASE( 'V' )
1281                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1282                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1283                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1284                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1285                  END SELECT
1286               END DO
1287            END DO
1288         END DO
1289      ENDIF
1290   END SUBROUTINE dom_vvl_orca_fix
1291
1292   !!======================================================================
1293END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.