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

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5235

Last change on this file since 5235 was 5235, checked in by davestorkey, 9 years ago

Clear svn keyword information in branch
2014/dev_r4650_UKMO11_restart_functionality

File size: 75.0 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 = .FALSE.              ! zstar  vertical coordinate
47   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde = .FALSE.             ! ztilde vertical coordinate
48   LOGICAL , PUBLIC                                      :: ln_vvl_layer = .FALSE.              ! level  vertical coordinate
49   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor = .FALSE.     ! ztilde vertical coordinate
51   LOGICAL , PUBLIC                                      :: ln_vvl_kepe = .FALSE.               ! 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 = .FALSE.      ! 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._wp)
408         CALL lbc_lnk( vn_td , 'V' , -1._wp)
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._wp )
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      INTEGER                             :: ji,jj,jk       ! dummy loop indices
591      !!----------------------------------------------------------------------
592
593      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
594      !
595      IF( kt == nit000 )   THEN
596         IF(lwp) WRITE(numout,*)
597         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
598         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
599      ENDIF
600      !
601      ! Time filter and swap of scale factors
602      ! =====================================
603      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
604      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
605         IF( neuler == 0 .AND. kt == nit000 ) THEN
606            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
607         ELSE
608            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
609            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
610         ENDIF
611         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
612      ENDIF
613      fsdept_b(:,:,:) = fsdept_n(:,:,:)
614      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
615
616      fse3t_n(:,:,:) = fse3t_a(:,:,:)
617      fse3u_n(:,:,:) = fse3u_a(:,:,:)
618      fse3v_n(:,:,:) = fse3v_a(:,:,:)
619
620      ! Compute all missing vertical scale factor and depths
621      ! ====================================================
622      ! Horizontal scale factor interpolations
623      ! --------------------------------------
624      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
625      ! - JC - hu_b, hv_b, hur_b, hvr_b also
626      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
627      ! Vertical scale factor interpolations
628      ! ------------------------------------
629      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
630      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
631      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
632      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
633      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
634      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
635      ! t- and w- points depth
636      ! ----------------------
637      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
638      fsdepw_n(:,:,1) = 0.0_wp
639      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
640      DO jj = 1,jpj
641         DO ji = 1,jpi
642            DO jk = 2,mikt(ji,jj)-1
643               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
644               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
645               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
646            END DO
647            IF (mikt(ji,jj) .GT. 1) THEN
648               jk = mikt(ji,jj)
649               fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
650               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
651               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
652            END IF
653            DO jk = mikt(ji,jj)+1, jpk
654               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
655               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
656               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
657            END DO
658         END DO
659      END DO
660      ! Local depth and Inverse of the local depth of the water column at u- and v- points
661      ! ----------------------------------------------------------------------------------
662      hu (:,:) = hu_a (:,:)
663      hv (:,:) = hv_a (:,:)
664
665      ! Inverse of the local depth
666      hur(:,:) = hur_a(:,:)
667      hvr(:,:) = hvr_a(:,:)
668
669      ! Local depth of the water column at t- points
670      ! --------------------------------------------
671      ht(:,:) = 0.
672      DO jk = 1, jpkm1
673         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
674      END DO
675
676      ! Write outputs
677      ! =============
678      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) )
679      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) )
680      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) )
681      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) )
682      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
683      IF( iom_use("e3tdef") )   &
684         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
685
686      ! write restart file
687      ! ==================
688      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
689      !
690      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
691
692   END SUBROUTINE dom_vvl_sf_swp
693
694
695   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
696      !!---------------------------------------------------------------------
697      !!                  ***  ROUTINE dom_vvl__interpol  ***
698      !!
699      !! ** Purpose :   interpolate scale factors from one grid point to another
700      !!
701      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
702      !!                - horizontal interpolation: grid cell surface averaging
703      !!                - vertical interpolation: simple averaging
704      !!----------------------------------------------------------------------
705      !! * Arguments
706      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
707      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
708      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
709      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
710      !! * Local declarations
711      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
712      LOGICAL ::   l_is_orca                                           ! local logical
713      !!----------------------------------------------------------------------
714      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
715         !
716      l_is_orca = .FALSE.
717      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
718
719      SELECT CASE ( pout )
720         !               ! ------------------------------------- !
721      CASE( 'U' )        ! interpolation from T-point to U-point !
722         !               ! ------------------------------------- !
723         ! horizontal surface weighted interpolation
724         DO jk = 1, jpk
725            DO jj = 1, jpjm1
726               DO ji = 1, fs_jpim1   ! vector opt.
727                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
728                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
729                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
730               END DO
731            END DO
732         END DO
733         !
734         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
735         ! boundary conditions
736         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
737         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
738         !               ! ------------------------------------- !
739      CASE( 'V' )        ! interpolation from T-point to V-point !
740         !               ! ------------------------------------- !
741         ! horizontal surface weighted interpolation
742         DO jk = 1, jpk
743            DO jj = 1, jpjm1
744               DO ji = 1, fs_jpim1   ! vector opt.
745                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
746                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
747                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
748               END DO
749            END DO
750         END DO
751         !
752         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
753         ! boundary conditions
754         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
755         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
756         !               ! ------------------------------------- !
757      CASE( 'F' )        ! interpolation from U-point to F-point !
758         !               ! ------------------------------------- !
759         ! horizontal surface weighted interpolation
760         DO jk = 1, jpk
761            DO jj = 1, jpjm1
762               DO ji = 1, fs_jpim1   ! vector opt.
763                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
764                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
765                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
766               END DO
767            END DO
768         END DO
769         !
770         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
771         ! boundary conditions
772         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
773         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
774         !               ! ------------------------------------- !
775      CASE( 'W' )        ! interpolation from T-point to W-point !
776         !               ! ------------------------------------- !
777         ! vertical simple interpolation
778         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
779         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
780         DO jk = 2, jpk
781            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
782               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
783         END DO
784         !               ! -------------------------------------- !
785      CASE( 'UW' )       ! interpolation from U-point to UW-point !
786         !               ! -------------------------------------- !
787         ! vertical simple interpolation
788         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
789         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
790         DO jk = 2, jpk
791            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
792               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
793         END DO
794         !               ! -------------------------------------- !
795      CASE( 'VW' )       ! interpolation from V-point to VW-point !
796         !               ! -------------------------------------- !
797         ! vertical simple interpolation
798         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
799         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
800         DO jk = 2, jpk
801            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
802               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
803         END DO
804      END SELECT
805      !
806
807      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
808
809   END SUBROUTINE dom_vvl_interpol
810
811   SUBROUTINE dom_vvl_rst( kt, cdrw )
812      !!---------------------------------------------------------------------
813      !!                   ***  ROUTINE dom_vvl_rst  ***
814      !!                     
815      !! ** Purpose :   Read or write VVL file in restart file
816      !!
817      !! ** Method  :   use of IOM library
818      !!                if the restart does not contain vertical scale factors,
819      !!                they are set to the _0 values
820      !!                if the restart does not contain vertical scale factors increments (z_tilde),
821      !!                they are set to 0.
822      !!----------------------------------------------------------------------
823      !! * Arguments
824      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
825      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
826      !! * Local declarations
827      INTEGER ::   jk
828      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
829      !!----------------------------------------------------------------------
830      !
831      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
832      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
833         !                                   ! ===============
834         IF( ln_rstart ) THEN                   !* Read the restart file
835            CALL rst_read_open                  !  open the restart file if necessary
836            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
837            !
838            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
839            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
840            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
841            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
842            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
843            !                             ! --------- !
844            !                             ! all cases !
845            !                             ! --------- !
846            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
847               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
848               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
849               ! needed to restart if land processor not computed
850               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
851               WHERE ( tmask(:,:,:) == 0.0_wp ) 
852                  fse3t_n(:,:,:) = e3t_0(:,:,:)
853                  fse3t_b(:,:,:) = e3t_0(:,:,:)
854               END WHERE
855               IF( neuler == 0 ) THEN
856                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
857               ENDIF
858            ELSE IF( id1 > 0 ) THEN
859               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
860               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
861               IF(lwp) write(numout,*) 'neuler is forced to 0'
862               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
863               fse3t_n(:,:,:) = fse3t_b(:,:,:)
864               neuler = 0
865            ELSE IF( id2 > 0 ) THEN
866               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
867               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
868               IF(lwp) write(numout,*) 'neuler is forced to 0'
869               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
870               fse3t_b(:,:,:) = fse3t_n(:,:,:)
871               neuler = 0
872            ELSE
873               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
874               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
875               IF(lwp) write(numout,*) 'neuler is forced to 0'
876               DO jk=1,jpk
877                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
878                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
879                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
880               END DO
881               fse3t_b(:,:,:) = fse3t_n(:,:,:)
882               neuler = 0
883            ENDIF
884            !                             ! ----------- !
885            IF( ln_vvl_zstar ) THEN       ! z_star case !
886               !                          ! ----------- !
887               IF( MIN( id3, id4 ) > 0 ) THEN
888                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
889               ENDIF
890               !                          ! ----------------------- !
891            ELSE                          ! z_tilde and layer cases !
892               !                          ! ----------------------- !
893               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
894                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
895                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
896               ELSE                            ! one at least array is missing
897                  tilde_e3t_b(:,:,:) = 0.0_wp
898                  tilde_e3t_n(:,:,:) = 0.0_wp
899               ENDIF
900               !                          ! ------------ !
901               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
902                  !                       ! ------------ !
903                  IF( id5 > 0 ) THEN  ! required array exists
904                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
905                  ELSE                ! array is missing
906                     hdiv_lf(:,:,:) = 0.0_wp
907                  ENDIF
908               ENDIF
909            ENDIF
910            !
911         ELSE                                   !* Initialize at "rest"
912            fse3t_b(:,:,:) = e3t_0(:,:,:)
913            fse3t_n(:,:,:) = e3t_0(:,:,:)
914            sshn(:,:) = 0.0_wp
915            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
916               tilde_e3t_b(:,:,:) = 0.0_wp
917               tilde_e3t_n(:,:,:) = 0.0_wp
918               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
919            END IF
920         ENDIF
921
922      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
923         !                                   ! ===================
924         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
925         !                                           ! --------- !
926         !                                           ! all cases !
927         !                                           ! --------- !
928         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
929         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
930         !                                           ! ----------------------- !
931         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
932            !                                        ! ----------------------- !
933            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
934            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
935         END IF
936         !                                           ! -------------!   
937         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
938            !                                        ! ------------ !
939            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
940         ENDIF
941
942      ENDIF
943      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
944
945   END SUBROUTINE dom_vvl_rst
946
947
948   SUBROUTINE dom_vvl_ctl
949      !!---------------------------------------------------------------------
950      !!                  ***  ROUTINE dom_vvl_ctl  ***
951      !!               
952      !! ** Purpose :   Control the consistency between namelist options
953      !!                for vertical coordinate
954      !!----------------------------------------------------------------------
955      INTEGER ::   ioptio
956      INTEGER ::   ios
957
958      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
959                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
960                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
961      !!----------------------------------------------------------------------
962
963      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
964      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
965901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
966
967      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
968      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
969902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
970      IF(lwm) WRITE ( numond, nam_vvl )
971
972      IF(lwp) THEN                    ! Namelist print
973         WRITE(numout,*)
974         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
975         WRITE(numout,*) '~~~~~~~~~~~'
976         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
977         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
978         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
979         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
980         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
981         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
982         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
983         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
984         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
985         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
986         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
987         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
988         IF( ln_vvl_ztilde_as_zstar ) THEN
989            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
990            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
991            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
992            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
993            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
994            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
995         ELSE
996            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
997            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
998            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
999            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
1000         ENDIF
1001         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1002         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1003      ENDIF
1004
1005      ioptio = 0                      ! Parameter control
1006      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1007      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
1008      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
1009      IF( ln_vvl_layer           )        ioptio = ioptio + 1
1010
1011      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1012      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1013
1014      IF(lwp) THEN                   ! Print the choice
1015         WRITE(numout,*)
1016         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1017         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1018         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1019         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1020         ! - ML - Option not developed yet
1021         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1022         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1023      ENDIF
1024
1025#if defined key_agrif
1026      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1027#endif
1028
1029   END SUBROUTINE dom_vvl_ctl
1030
1031   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1032      !!---------------------------------------------------------------------
1033      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1034      !!                     
1035      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1036      !!                scale factors at locations that have been individually
1037      !!                modified in domhgr. Such modifications break the
1038      !!                relationship between e12t and e1u*e2u etc.
1039      !!                Recompute some scale factors ignoring the modified metric.
1040      !!----------------------------------------------------------------------
1041      !! * Arguments
1042      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1043      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1044      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1045      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1046      !! * Local declarations
1047      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1048      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1049      !! acc
1050      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1051      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1052      !!
1053      !                                                ! =====================
1054      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1055         !                                             ! =====================
1056      !! acc
1057         IF( nn_cla == 0 ) THEN
1058            !
1059            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1060            ij0 = 102   ;   ij1 = 102
1061            DO jk = 1, jpkm1
1062               DO jj = mj0(ij0), mj1(ij1)
1063                  DO ji = mi0(ii0), mi1(ii1)
1064                     SELECT CASE ( pout )
1065                     CASE( 'U' )
1066                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1067                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1068                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1069                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1070                     CASE( 'F' )
1071                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1072                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1073                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1074                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1075                     END SELECT
1076                  END DO
1077               END DO
1078            END DO
1079            !
1080            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1081            ij0 =  88   ;   ij1 =  88
1082            DO jk = 1, jpkm1
1083               DO jj = mj0(ij0), mj1(ij1)
1084                  DO ji = mi0(ii0), mi1(ii1)
1085                     SELECT CASE ( pout )
1086                     CASE( 'U' )
1087                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1088                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1089                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1090                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1091                     CASE( 'V' )
1092                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1093                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1094                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1095                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1096                     CASE( 'F' )
1097                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1098                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1099                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1100                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1101                     END SELECT
1102                  END DO
1103               END DO
1104            END DO
1105         ENDIF
1106
1107         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1108         ij0 = 116   ;   ij1 = 116
1109         DO jk = 1, jpkm1
1110            DO jj = mj0(ij0), mj1(ij1)
1111               DO ji = mi0(ii0), mi1(ii1)
1112                  SELECT CASE ( pout )
1113                  CASE( 'U' )
1114                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1115                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1116                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1117                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1118                  CASE( 'F' )
1119                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1120                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1121                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1122                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1123                  END SELECT
1124               END DO
1125            END DO
1126         END DO
1127      ENDIF
1128      !
1129         !                                             ! =====================
1130      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1131         !                                             ! =====================
1132         !
1133         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1134         ij0 = 200   ;   ij1 = 200
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( 'U' )
1140                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1141                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1142                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1143                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1144                  CASE( 'F' )
1145                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1146                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1147                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1148                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1149                  END SELECT
1150               END DO
1151            END DO
1152         END DO
1153         !
1154         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1155         ij0 = 208   ;   ij1 = 208
1156         DO jk = 1, jpkm1
1157            DO jj = mj0(ij0), mj1(ij1)
1158               DO ji = mi0(ii0), mi1(ii1)
1159                  SELECT CASE ( pout )
1160                  CASE( 'U' )
1161                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1162                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1163                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1164                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1165                  CASE( 'F' )
1166                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1167                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1168                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1169                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1170                  END SELECT
1171               END DO
1172            END DO
1173         END DO
1174         !
1175         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1176         ij0 = 124   ;   ij1 = 125
1177         DO jk = 1, jpkm1
1178            DO jj = mj0(ij0), mj1(ij1)
1179               DO ji = mi0(ii0), mi1(ii1)
1180                  SELECT CASE ( pout )
1181                  CASE( 'V' )
1182                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1183                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1184                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1185                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1186                  END SELECT
1187               END DO
1188            END DO
1189         END DO
1190         !
1191         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1192         ij0 = 124   ;   ij1 = 125
1193         DO jk = 1, jpkm1
1194            DO jj = mj0(ij0), mj1(ij1)
1195               DO ji = mi0(ii0), mi1(ii1)
1196                  SELECT CASE ( pout )
1197                  CASE( 'V' )
1198                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1199                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1200                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1201                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1202                  END SELECT
1203               END DO
1204            END DO
1205         END DO
1206         !
1207         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1208         ij0 = 124   ;   ij1 = 125
1209         DO jk = 1, jpkm1
1210            DO jj = mj0(ij0), mj1(ij1)
1211               DO ji = mi0(ii0), mi1(ii1)
1212                  SELECT CASE ( pout )
1213                  CASE( 'V' )
1214                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1215                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1216                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1217                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1218                  END SELECT
1219               END DO
1220            END DO
1221         END DO
1222         !
1223         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1224         ij0 = 124   ;   ij1 = 125
1225         DO jk = 1, jpkm1
1226            DO jj = mj0(ij0), mj1(ij1)
1227               DO ji = mi0(ii0), mi1(ii1)
1228                  SELECT CASE ( pout )
1229                  CASE( 'V' )
1230                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1231                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1232                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1233                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1234                  END SELECT
1235               END DO
1236            END DO
1237         END DO
1238         !
1239         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1240         ij0 = 141   ;   ij1 = 142
1241         DO jk = 1, jpkm1
1242            DO jj = mj0(ij0), mj1(ij1)
1243               DO ji = mi0(ii0), mi1(ii1)
1244                  SELECT CASE ( pout )
1245                  CASE( 'V' )
1246                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1247                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1248                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1249                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1250                  END SELECT
1251               END DO
1252            END DO
1253         END DO
1254         !
1255         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1256         ij0 = 141   ;   ij1 = 142
1257         DO jk = 1, jpkm1
1258            DO jj = mj0(ij0), mj1(ij1)
1259               DO ji = mi0(ii0), mi1(ii1)
1260                  SELECT CASE ( pout )
1261                  CASE( 'V' )
1262                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1263                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1264                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1265                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1266                  END SELECT
1267               END DO
1268            END DO
1269         END DO
1270      ENDIF
1271         !                                             ! =====================
1272      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1273         !                                             ! =====================
1274         !
1275         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1276         ij0 = 327   ;   ij1 = 327
1277         DO jk = 1, jpkm1
1278            DO jj = mj0(ij0), mj1(ij1)
1279               DO ji = mi0(ii0), mi1(ii1)
1280                  SELECT CASE ( pout )
1281                  CASE( 'U' )
1282                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1283                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1284                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1285                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1286                  CASE( 'F' )
1287                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1288                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1289                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1290                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1291                  END SELECT
1292               END DO
1293            END DO
1294         END DO
1295         !
1296         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1297         ij0 = 343   ;   ij1 = 343
1298         DO jk = 1, jpkm1
1299            DO jj = mj0(ij0), mj1(ij1)
1300               DO ji = mi0(ii0), mi1(ii1)
1301                  SELECT CASE ( pout )
1302                  CASE( 'U' )
1303                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1304                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1305                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1306                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1307                  CASE( 'F' )
1308                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1309                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1310                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1311                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1312                  END SELECT
1313               END DO
1314            END DO
1315         END DO
1316         !
1317         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1318         ij0 = 232   ;   ij1 = 232
1319         DO jk = 1, jpkm1
1320            DO jj = mj0(ij0), mj1(ij1)
1321               DO ji = mi0(ii0), mi1(ii1)
1322                  SELECT CASE ( pout )
1323                  CASE( 'U' )
1324                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1325                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1326                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1327                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1328                  CASE( 'F' )
1329                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1330                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1331                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1332                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1333                  END SELECT
1334               END DO
1335            END DO
1336         END DO
1337         !
1338         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1339         ij0 = 232   ;   ij1 = 232
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( 'U' )
1345                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1346                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1347                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1348                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1349                  CASE( 'F' )
1350                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1351                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1352                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1353                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1354                  END SELECT
1355               END DO
1356            END DO
1357         END DO
1358         !
1359         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1360         ij0 = 270   ;   ij1 = 270
1361         DO jk = 1, jpkm1
1362            DO jj = mj0(ij0), mj1(ij1)
1363               DO ji = mi0(ii0), mi1(ii1)
1364                  SELECT CASE ( pout )
1365                  CASE( 'U' )
1366                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1367                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1368                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1369                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1370                  CASE( 'F' )
1371                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1372                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1373                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1374                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1375                  END SELECT
1376               END DO
1377            END DO
1378         END DO
1379         !
1380         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1381         ij0 = 232   ;   ij1 = 233
1382         DO jk = 1, jpkm1
1383            DO jj = mj0(ij0), mj1(ij1)
1384               DO ji = mi0(ii0), mi1(ii1)
1385                  SELECT CASE ( pout )
1386                  CASE( 'V' )
1387                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1388                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1389                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1390                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1391                  END SELECT
1392               END DO
1393            END DO
1394         END DO
1395         !
1396         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1397         ij0 = 276   ;   ij1 = 276
1398         DO jk = 1, jpkm1
1399            DO jj = mj0(ij0), mj1(ij1)
1400               DO ji = mi0(ii0), mi1(ii1)
1401                  SELECT CASE ( pout )
1402                  CASE( 'V' )
1403                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1404                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1405                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1406                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1407                  END SELECT
1408               END DO
1409            END DO
1410         END DO
1411      ENDIF
1412   END SUBROUTINE dom_vvl_orca_fix
1413
1414   !!======================================================================
1415END MODULE domvvl
1416
1417
1418
Note: See TracBrowser for help on using the repository browser.