source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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