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 trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4795

Last change on this file since 4795 was 4795, checked in by jchanut, 10 years ago

Read hdiv_lf with ztilde option, see ticket #1387

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