source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4292

Last change on this file since 4292 was 4292, checked in by cetlod, 8 years ago

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

  • Property svn:keywords set to Id
File size: 69.4 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) &
347                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
348                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
349                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
350                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
351                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
352               END DO
353            END DO
354         END DO
355         ! b - correction for last oceanic u-v points
356         DO jj = 1, jpj
357            DO ji = 1, jpi
358               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
359               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
360            END DO
361         END DO
362         ! c - second derivative: divergence of diffusive fluxes
363         DO jk = 1, jpkm1
364            DO jj = 2, jpjm1
365               DO ji = fs_2, fs_jpim1   ! vector opt.
366                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
367                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
368                     &                                            ) * r1_e12t(ji,jj)
369               END DO
370            END DO
371         END DO
372         ! d - thickness diffusion transport: boundary conditions
373         !     (stored for tracer advction and continuity equation)
374         CALL lbc_lnk( un_td , 'U' , -1.)
375         CALL lbc_lnk( vn_td , 'V' , -1.)
376
377         ! 4 - Time stepping of baroclinic scale factors
378         ! ---------------------------------------------
379         ! Leapfrog time stepping
380         ! ~~~~~~~~~~~~~~~~~~~~~~
381         IF( neuler == 0 .AND. kt == nit000 ) THEN
382            z2dt =  rdt
383         ELSE
384            z2dt = 2.0_wp * rdt
385         ENDIF
386         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
387         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
388
389         ! Maximum deformation control
390         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
391         ze3t(:,:,jpk) = 0.0_wp
392         DO jk = 1, jpkm1
393            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
394         END DO
395         z_tmax = MAXVAL( ze3t(:,:,:) )
396         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
397         z_tmin = MINVAL( ze3t(:,:,:) )
398         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
399         ! - ML - test: for the moment, stop simulation for too large e3_t variations
400         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
401            IF( lk_mpp ) THEN
402               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
403               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
404            ELSE
405               ijk_max = MAXLOC( ze3t(:,:,:) )
406               ijk_max(1) = ijk_max(1) + nimpp - 1
407               ijk_max(2) = ijk_max(2) + njmpp - 1
408               ijk_min = MINLOC( ze3t(:,:,:) )
409               ijk_min(1) = ijk_min(1) + nimpp - 1
410               ijk_min(2) = ijk_min(2) + njmpp - 1
411            ENDIF
412            IF (lwp) THEN
413               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
414               WRITE(numout, *) 'at i, j, k=', ijk_max
415               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
416               WRITE(numout, *) 'at i, j, k=', ijk_min           
417               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
418            ENDIF
419         ENDIF
420         ! - ML - end test
421         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
422         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
423         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
424
425         ! Add "tilda" part to the after scale factor
426         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427         fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:)
428
429         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
430         ! ==================================================================================
431         ! add e3t(n-1) "star" Asselin-filtered
432         DO jk = 1, jpkm1
433            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk)
434         END DO
435         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
436         ! - ML - baroclinicity error should be better treated in the future
437         !        i.e. locally and not spread over the water column.
438         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
439         zht(:,:) = 0.
440         DO jk = 1, jpkm1
441            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
442         END DO
443         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
444         DO jk = 1, jpkm1
445            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
446         END DO
447
448      ENDIF
449
450      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
451         !
452         IF( lwp ) WRITE(numout, *) 'kt =', kt
453         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
454            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
455            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
456            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
457         END IF
458         !
459         zht(:,:) = 0.0_wp
460         DO jk = 1, jpkm1
461            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
462         END DO
463         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
464         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
465         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
466         !
467         zht(:,:) = 0.0_wp
468         DO jk = 1, jpkm1
469            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
470         END DO
471         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
472         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
473         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
474         !
475         zht(:,:) = 0.0_wp
476         DO jk = 1, jpkm1
477            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
478         END DO
479         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
480         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
481         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
482         !
483         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
484         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
485         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
486         !
487         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
488         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
489         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
490         !
491         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
492         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
493         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
494      END IF
495
496      ! *********************************** !
497      ! After scale factors at u- v- points !
498      ! *********************************** !
499
500      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
501      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
502
503      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
504      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
505
506      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
507
508   END SUBROUTINE dom_vvl_sf_nxt
509
510
511   SUBROUTINE dom_vvl_sf_swp( kt )
512      !!----------------------------------------------------------------------
513      !!                ***  ROUTINE dom_vvl_sf_swp  ***
514      !!                   
515      !! ** Purpose :  compute time filter and swap of scale factors
516      !!               compute all depths and related variables for next time step
517      !!               write outputs and restart file
518      !!
519      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
520      !!               - reconstruct scale factor at other grid points (interpolate)
521      !!               - recompute depths and water height fields
522      !!
523      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
524      !!               - Recompute:
525      !!                    fse3(u/v)_b       
526      !!                    fse3w_n           
527      !!                    fse3(u/v)w_b     
528      !!                    fse3(u/v)w_n     
529      !!                    fsdept_n, fsdepw_n  and fsde3w_n
530      !!                    h(u/v) and h(u/v)r
531      !!
532      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
533      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
534      !!----------------------------------------------------------------------
535      !! * Arguments
536      INTEGER, INTENT( in )               :: kt       ! time step
537      !! * Local declarations
538      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
539      INTEGER                             :: jk       ! dummy loop indices
540      !!----------------------------------------------------------------------
541
542      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
543      !
544      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
545      !
546      IF( kt == nit000 )   THEN
547         IF(lwp) WRITE(numout,*)
548         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
549         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
550      ENDIF
551      !
552      ! Time filter and swap of scale factors
553      ! =====================================
554      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
555      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
556         IF( neuler == 0 .AND. kt == nit000 ) THEN
557            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
558         ELSE
559            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
560            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
561         ENDIF
562         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
563      ENDIF
564      fse3t_n(:,:,:) = fse3t_a(:,:,:)
565      fse3u_n(:,:,:) = fse3u_a(:,:,:)
566      fse3v_n(:,:,:) = fse3v_a(:,:,:)
567
568      ! Compute all missing vertical scale factor and depths
569      ! ====================================================
570      ! Horizontal scale factor interpolations
571      ! --------------------------------------
572      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
573      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
574      ! Vertical scale factor interpolations
575      ! ------------------------------------
576      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
577      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
578      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
579      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
580      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
581      ! t- and w- points depth
582      ! ----------------------
583      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
584      fsdepw_n(:,:,1) = 0.0_wp
585      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
586      DO jk = 2, jpk
587         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
588         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
589         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
590      END DO
591      ! Local depth and Inverse of the local depth of the water column at u- and v- points
592      ! ----------------------------------------------------------------------------------
593      hu(:,:) = 0.
594      hv(:,:) = 0.
595      DO jk = 1, jpk
596         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
597         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
598      END DO
599      ! Inverse of the local depth
600      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
601      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) )
602
603      ! Write outputs
604      ! =============
605      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
606      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
607      CALL iom_put( "dept_n" , fsde3w_n (:,:,:) )
608      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
609
610      ! write restart file
611      ! ==================
612      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
613      !
614      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
615      !
616      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
617
618   END SUBROUTINE dom_vvl_sf_swp
619
620
621   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
622      !!---------------------------------------------------------------------
623      !!                  ***  ROUTINE dom_vvl__interpol  ***
624      !!
625      !! ** Purpose :   interpolate scale factors from one grid point to another
626      !!
627      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
628      !!                - horizontal interpolation: grid cell surface averaging
629      !!                - vertical interpolation: simple averaging
630      !!----------------------------------------------------------------------
631      !! * Arguments
632      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
633      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
634      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
635      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
636      !! * Local declarations
637      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
638      LOGICAL ::   l_is_orca                                           ! local logical
639      !!----------------------------------------------------------------------
640      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
641         !
642      l_is_orca = .FALSE.
643      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
644
645      SELECT CASE ( pout )
646         !               ! ------------------------------------- !
647      CASE( 'U' )        ! interpolation from T-point to U-point !
648         !               ! ------------------------------------- !
649         ! horizontal surface weighted interpolation
650         DO jk = 1, jpk
651            DO jj = 1, jpjm1
652               DO ji = 1, fs_jpim1   ! vector opt.
653                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
654                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
655                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
656               END DO
657            END DO
658         END DO
659         !
660         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
661         ! boundary conditions
662         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
663         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
664         !               ! ------------------------------------- !
665      CASE( 'V' )        ! interpolation from T-point to V-point !
666         !               ! ------------------------------------- !
667         ! horizontal surface weighted interpolation
668         DO jk = 1, jpk
669            DO jj = 1, jpjm1
670               DO ji = 1, fs_jpim1   ! vector opt.
671                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
672                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
673                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
674               END DO
675            END DO
676         END DO
677         !
678         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
679         ! boundary conditions
680         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
681         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
682         !               ! ------------------------------------- !
683      CASE( 'F' )        ! interpolation from U-point to F-point !
684         !               ! ------------------------------------- !
685         ! horizontal surface weighted interpolation
686         DO jk = 1, jpk
687            DO jj = 1, jpjm1
688               DO ji = 1, fs_jpim1   ! vector opt.
689                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
690                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
691                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
692               END DO
693            END DO
694         END DO
695         !
696         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
697         ! boundary conditions
698         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
699         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
700         !               ! ------------------------------------- !
701      CASE( 'W' )        ! interpolation from T-point to W-point !
702         !               ! ------------------------------------- !
703         ! vertical simple interpolation
704         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
705         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
706         DO jk = 2, jpk
707            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
708               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
709         END DO
710         !               ! -------------------------------------- !
711      CASE( 'UW' )       ! interpolation from U-point to UW-point !
712         !               ! -------------------------------------- !
713         ! vertical simple interpolation
714         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
715         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
716         DO jk = 2, jpk
717            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
718               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
719         END DO
720         !               ! -------------------------------------- !
721      CASE( 'VW' )       ! interpolation from V-point to VW-point !
722         !               ! -------------------------------------- !
723         ! vertical simple interpolation
724         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
725         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
726         DO jk = 2, jpk
727            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
728               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
729         END DO
730      END SELECT
731      !
732
733      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
734
735   END SUBROUTINE dom_vvl_interpol
736
737   SUBROUTINE dom_vvl_rst( kt, cdrw )
738      !!---------------------------------------------------------------------
739      !!                   ***  ROUTINE dom_vvl_rst  ***
740      !!                     
741      !! ** Purpose :   Read or write VVL file in restart file
742      !!
743      !! ** Method  :   use of IOM library
744      !!                if the restart does not contain vertical scale factors,
745      !!                they are set to the _0 values
746      !!                if the restart does not contain vertical scale factors increments (z_tilde),
747      !!                they are set to 0.
748      !!----------------------------------------------------------------------
749      !! * Arguments
750      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
751      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
752      !! * Local declarations
753      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
754      !!----------------------------------------------------------------------
755      !
756      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
757      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
758         !                                   ! ===============
759         IF( ln_rstart ) THEN                   !* Read the restart file
760            CALL rst_read_open                  !  open the restart file if necessary
761            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
762            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
763            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
764            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
765            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
766            !                             ! --------- !
767            !                             ! all cases !
768            !                             ! --------- !
769            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
770               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
771               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
772               IF( neuler == 0 ) THEN
773                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
774               ENDIF
775            ELSE IF( id1 > 0 ) THEN
776               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
777               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
778               fse3t_b(:,:,:) = fse3t_n(:,:,:)
779            ELSE                                 ! one at least array is missing
780               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
781            ENDIF
782            !                             ! ----------- !
783            IF( ln_vvl_zstar ) THEN       ! z_star case !
784               !                          ! ----------- !
785               IF( MIN( id3, id4 ) > 0 ) THEN
786                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
787               ENDIF
788               !                          ! ----------------------- !
789            ELSE                          ! z_tilde and layer cases !
790               !                          ! ----------------------- !
791               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
792                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
793                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
794               ELSE                            ! one at least array is missing
795                  tilde_e3t_b(:,:,:) = 0.0_wp
796                  tilde_e3t_n(:,:,:) = 0.0_wp
797               ENDIF
798               !                          ! ------------ !
799               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
800                  !                       ! ------------ !
801                  IF( id5 > 0 ) THEN  ! required array exists
802                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
803                  ELSE                ! array is missing
804                     hdiv_lf(:,:,:) = 0.0_wp
805                  ENDIF
806               ENDIF
807            ENDIF
808            !
809         ELSE                                   !* Initialize at "rest"
810            fse3t_b(:,:,:) = e3t_0(:,:,:)
811            fse3t_n(:,:,:) = e3t_0(:,:,:)
812            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
813               tilde_e3t_b(:,:,:) = 0.0_wp
814               tilde_e3t_n(:,:,:) = 0.0_wp
815               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
816            END IF
817         ENDIF
818
819      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
820         !                                   ! ===================
821         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
822         !                                           ! --------- !
823         !                                           ! all cases !
824         !                                           ! --------- !
825         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
826         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
827         !                                           ! ----------------------- !
828         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
829            !                                        ! ----------------------- !
830            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
831            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
832         END IF
833         !                                           ! -------------!   
834         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
835            !                                        ! ------------ !
836            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
837         ENDIF
838
839      ENDIF
840      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
841
842   END SUBROUTINE dom_vvl_rst
843
844
845   SUBROUTINE dom_vvl_ctl
846      !!---------------------------------------------------------------------
847      !!                  ***  ROUTINE dom_vvl_ctl  ***
848      !!               
849      !! ** Purpose :   Control the consistency between namelist options
850      !!                for vertical coordinate
851      !!----------------------------------------------------------------------
852      INTEGER ::   ioptio
853
854      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
855                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
856                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
857      !!----------------------------------------------------------------------
858
859      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
860      READ   ( numnam, nam_vvl )
861
862      IF(lwp) THEN                    ! Namelist print
863         WRITE(numout,*)
864         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
865         WRITE(numout,*) '~~~~~~~~~~~'
866         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
867         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
868         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
869         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
870         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
871         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
872         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
873         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
874         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
875         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
876         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
877         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
878         IF( ln_vvl_ztilde_as_zstar ) THEN
879            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
880            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
881            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
882            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
883            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
884            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
885         ELSE
886            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
887            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
888            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
889            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
890         ENDIF
891         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
892         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
893      ENDIF
894
895      ioptio = 0                      ! Parameter control
896      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
897      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
898      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
899      IF( ln_vvl_layer           )        ioptio = ioptio + 1
900
901      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
902
903      IF(lwp) THEN                   ! Print the choice
904         WRITE(numout,*)
905         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
906         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
907         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
908         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
909         ! - ML - Option not developed yet
910         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
911         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
912      ENDIF
913
914   END SUBROUTINE dom_vvl_ctl
915
916   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
917      !!---------------------------------------------------------------------
918      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
919      !!                     
920      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
921      !!                scale factors at locations that have been individually
922      !!                modified in domhgr. Such modifications break the
923      !!                relationship between e12t and e1u*e2u etc.
924      !!                Recompute some scale factors ignoring the modified metric.
925      !!----------------------------------------------------------------------
926      !! * Arguments
927      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
928      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
929      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
930      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
931      !! * Local declarations
932      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
933      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
934      !! acc
935      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
936      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
937      !!
938      !                                                ! =====================
939      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
940         !                                             ! =====================
941      !! acc
942         IF( nn_cla == 0 ) THEN
943            !
944            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
945            ij0 = 102   ;   ij1 = 102
946            DO jk = 1, jpkm1
947               DO jj = mj0(ij0), mj1(ij1)
948                  DO ji = mi0(ii0), mi1(ii1)
949                     SELECT CASE ( pout )
950                     CASE( 'U' )
951                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
952                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
953                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
954                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
955                     CASE( 'F' )
956                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
957                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
958                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
959                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
960                     END SELECT
961                  END DO
962               END DO
963            END DO
964            !
965            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
966            ij0 =  88   ;   ij1 =  88
967            DO jk = 1, jpkm1
968               DO jj = mj0(ij0), mj1(ij1)
969                  DO ji = mi0(ii0), mi1(ii1)
970                     SELECT CASE ( pout )
971                     CASE( 'U' )
972                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
973                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
974                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
975                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
976                     CASE( 'V' )
977                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
978                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
979                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
980                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
981                     CASE( 'F' )
982                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
983                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
984                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
985                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
986                     END SELECT
987                  END DO
988               END DO
989            END DO
990         ENDIF
991
992         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
993         ij0 = 116   ;   ij1 = 116
994         DO jk = 1, jpkm1
995            DO jj = mj0(ij0), mj1(ij1)
996               DO ji = mi0(ii0), mi1(ii1)
997                  SELECT CASE ( pout )
998                  CASE( 'U' )
999                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1000                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1001                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1002                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1003                  CASE( 'F' )
1004                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1005                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1006                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1007                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1008                  END SELECT
1009               END DO
1010            END DO
1011         END DO
1012      ENDIF
1013      !
1014         !                                             ! =====================
1015      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1016         !                                             ! =====================
1017         !
1018         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1019         ij0 = 200   ;   ij1 = 200
1020         DO jk = 1, jpkm1
1021            DO jj = mj0(ij0), mj1(ij1)
1022               DO ji = mi0(ii0), mi1(ii1)
1023                  SELECT CASE ( pout )
1024                  CASE( 'U' )
1025                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1026                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1027                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1028                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1029                  CASE( 'F' )
1030                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1031                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1032                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1033                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1034                  END SELECT
1035               END DO
1036            END DO
1037         END DO
1038         !
1039         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1040         ij0 = 208   ;   ij1 = 208
1041         DO jk = 1, jpkm1
1042            DO jj = mj0(ij0), mj1(ij1)
1043               DO ji = mi0(ii0), mi1(ii1)
1044                  SELECT CASE ( pout )
1045                  CASE( 'U' )
1046                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1047                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1048                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1049                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1050                  CASE( 'F' )
1051                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1052                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1053                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1054                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1055                  END SELECT
1056               END DO
1057            END DO
1058         END DO
1059         !
1060         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1061         ij0 = 124   ;   ij1 = 125
1062         DO jk = 1, jpkm1
1063            DO jj = mj0(ij0), mj1(ij1)
1064               DO ji = mi0(ii0), mi1(ii1)
1065                  SELECT CASE ( pout )
1066                  CASE( 'V' )
1067                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1068                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1069                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1070                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1071                  END SELECT
1072               END DO
1073            END DO
1074         END DO
1075         !
1076         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1077         ij0 = 124   ;   ij1 = 125
1078         DO jk = 1, jpkm1
1079            DO jj = mj0(ij0), mj1(ij1)
1080               DO ji = mi0(ii0), mi1(ii1)
1081                  SELECT CASE ( pout )
1082                  CASE( 'V' )
1083                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1084                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1085                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1086                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1087                  END SELECT
1088               END DO
1089            END DO
1090         END DO
1091         !
1092         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1093         ij0 = 124   ;   ij1 = 125
1094         DO jk = 1, jpkm1
1095            DO jj = mj0(ij0), mj1(ij1)
1096               DO ji = mi0(ii0), mi1(ii1)
1097                  SELECT CASE ( pout )
1098                  CASE( 'V' )
1099                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1100                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1101                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1102                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1103                  END SELECT
1104               END DO
1105            END DO
1106         END DO
1107         !
1108         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1109         ij0 = 124   ;   ij1 = 125
1110         DO jk = 1, jpkm1
1111            DO jj = mj0(ij0), mj1(ij1)
1112               DO ji = mi0(ii0), mi1(ii1)
1113                  SELECT CASE ( pout )
1114                  CASE( 'V' )
1115                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1116                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1117                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1118                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1119                  END SELECT
1120               END DO
1121            END DO
1122         END DO
1123         !
1124         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1125         ij0 = 141   ;   ij1 = 142
1126         DO jk = 1, jpkm1
1127            DO jj = mj0(ij0), mj1(ij1)
1128               DO ji = mi0(ii0), mi1(ii1)
1129                  SELECT CASE ( pout )
1130                  CASE( 'V' )
1131                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1132                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1133                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1134                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1135                  END SELECT
1136               END DO
1137            END DO
1138         END DO
1139         !
1140         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1141         ij0 = 141   ;   ij1 = 142
1142         DO jk = 1, jpkm1
1143            DO jj = mj0(ij0), mj1(ij1)
1144               DO ji = mi0(ii0), mi1(ii1)
1145                  SELECT CASE ( pout )
1146                  CASE( 'V' )
1147                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1148                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1149                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1150                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1151                  END SELECT
1152               END DO
1153            END DO
1154         END DO
1155      ENDIF
1156         !                                             ! =====================
1157      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1158         !                                             ! =====================
1159         !
1160         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1161         ij0 = 327   ;   ij1 = 327
1162         DO jk = 1, jpkm1
1163            DO jj = mj0(ij0), mj1(ij1)
1164               DO ji = mi0(ii0), mi1(ii1)
1165                  SELECT CASE ( pout )
1166                  CASE( 'U' )
1167                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1168                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1169                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1170                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1171                  CASE( 'F' )
1172                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1173                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1174                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1175                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1176                  END SELECT
1177               END DO
1178            END DO
1179         END DO
1180         !
1181         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1182         ij0 = 343   ;   ij1 = 343
1183         DO jk = 1, jpkm1
1184            DO jj = mj0(ij0), mj1(ij1)
1185               DO ji = mi0(ii0), mi1(ii1)
1186                  SELECT CASE ( pout )
1187                  CASE( 'U' )
1188                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1189                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1190                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1191                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1192                  CASE( 'F' )
1193                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1194                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1195                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1196                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1197                  END SELECT
1198               END DO
1199            END DO
1200         END DO
1201         !
1202         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1203         ij0 = 232   ;   ij1 = 232
1204         DO jk = 1, jpkm1
1205            DO jj = mj0(ij0), mj1(ij1)
1206               DO ji = mi0(ii0), mi1(ii1)
1207                  SELECT CASE ( pout )
1208                  CASE( 'U' )
1209                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1210                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1211                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1212                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1213                  CASE( 'F' )
1214                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1215                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1216                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1217                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1218                  END SELECT
1219               END DO
1220            END DO
1221         END DO
1222         !
1223         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1224         ij0 = 232   ;   ij1 = 232
1225         DO jk = 1, jpkm1
1226            DO jj = mj0(ij0), mj1(ij1)
1227               DO ji = mi0(ii0), mi1(ii1)
1228                  SELECT CASE ( pout )
1229                  CASE( 'U' )
1230                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1231                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1232                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1233                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1234                  CASE( 'F' )
1235                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1236                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1237                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1238                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1239                  END SELECT
1240               END DO
1241            END DO
1242         END DO
1243         !
1244         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1245         ij0 = 270   ;   ij1 = 270
1246         DO jk = 1, jpkm1
1247            DO jj = mj0(ij0), mj1(ij1)
1248               DO ji = mi0(ii0), mi1(ii1)
1249                  SELECT CASE ( pout )
1250                  CASE( 'U' )
1251                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1252                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1253                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1254                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1255                  CASE( 'F' )
1256                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1257                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1258                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1259                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1260                  END SELECT
1261               END DO
1262            END DO
1263         END DO
1264         !
1265         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1266         ij0 = 232   ;   ij1 = 233
1267         DO jk = 1, jpkm1
1268            DO jj = mj0(ij0), mj1(ij1)
1269               DO ji = mi0(ii0), mi1(ii1)
1270                  SELECT CASE ( pout )
1271                  CASE( 'V' )
1272                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1273                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1274                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1275                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1276                  END SELECT
1277               END DO
1278            END DO
1279         END DO
1280         !
1281         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1282         ij0 = 276   ;   ij1 = 276
1283         DO jk = 1, jpkm1
1284            DO jj = mj0(ij0), mj1(ij1)
1285               DO ji = mi0(ii0), mi1(ii1)
1286                  SELECT CASE ( pout )
1287                  CASE( 'V' )
1288                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1289                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1290                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1291                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1292                  END SELECT
1293               END DO
1294            END DO
1295         END DO
1296      ENDIF
1297   END SUBROUTINE dom_vvl_orca_fix
1298
1299   !!======================================================================
1300END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.