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

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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4486

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

Finalize Time split and AGRIF (tickets #106 and #107) + ticket #1240

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