source: branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5066

Last change on this file since 5066 was 5066, checked in by rfurner, 6 years ago

added current state of wetting and drying code to test…note it does not work

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