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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4747

Last change on this file since 4747 was 4747, checked in by mathiot, 10 years ago

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