source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5018

Last change on this file since 5018 was 5018, checked in by hliu, 6 years ago

a small change for Wetting/Drying? branch

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