New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 6 years ago

Remove svn keywords

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