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 NEMO/branches/2020/ticket2487/src/OCE/DOM – NEMO

source: NEMO/branches/2020/ticket2487/src/OCE/DOM/domvvl.F90 @ 15418

Last change on this file since 15418 was 15418, checked in by smueller, 12 months ago

Addition of vertical-scale-factor and water-column-height updates to complete initial SSH adjustments due to non-zero sea-ice/snow mass (ticket #2487)

  • Property svn:keywords set to Id
File size: 55.8 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
9   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
14   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
15   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
16   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
17   !!   dom_vvl_rst      : read/write restart file
18   !!   dom_vvl_ctl      : Check the vvl options
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers
21   USE phycst          ! physical constant
22   USE dom_oce         ! ocean space and time domain
23   USE sbc_oce         ! ocean surface boundary condition
24   USE wet_dry         ! wetting and drying
25   USE usrdef_istate   ! user defined initial state (wad only)
26   USE restart         ! ocean restart
27   !
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O manager library
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE timing          ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC  dom_vvl_init       ! called by domain.F90
38   PUBLIC  dom_vvl_zgr        ! called in modules iceistate and sbcice_cice
39   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
40   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
41   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
42
43   !                                                      !!* Namelist nam_vvl
44   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
45   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
46   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
50   !                                                       ! conservation: not used yet
51   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
52   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
53   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
54   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
55   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
56
57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
58   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
62   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
63
64   !! * Substitutions
65#  include "vectopt_loop_substitute.h90"
66   !!----------------------------------------------------------------------
67   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
68   !! $Id$
69   !! Software governed by the CeCILL license (see ./LICENSE)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   INTEGER FUNCTION dom_vvl_alloc()
74      !!----------------------------------------------------------------------
75      !!                ***  FUNCTION dom_vvl_alloc  ***
76      !!----------------------------------------------------------------------
77      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
78      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
79         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
80            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
81            &      STAT = dom_vvl_alloc        )
82         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
83         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
84         un_td = 0._wp
85         vn_td = 0._wp
86      ENDIF
87      IF( ln_vvl_ztilde ) THEN
88         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
89         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
90         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
91      ENDIF
92      !
93   END FUNCTION dom_vvl_alloc
94
95
96   SUBROUTINE dom_vvl_init
97      !!----------------------------------------------------------------------
98      !!                ***  ROUTINE dom_vvl_init  ***
99      !!                   
100      !! ** Purpose :  Initialization of all scale factors, depths
101      !!               and water column heights
102      !!
103      !! ** Method  :  - use restart file and/or initialize
104      !!               - interpolate scale factors
105      !!
106      !! ** Action  : - call dom_vvl_zgr
107      !!
108      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
109      !!----------------------------------------------------------------------
110      !
111      IF(lwp) WRITE(numout,*)
112      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
113      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
114      !
115      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
116      !
117      !                    ! Allocate module arrays
118      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
119      !
120      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
121      CALL dom_vvl_rst( nit000, 'READ' )
122      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
123      !
124      CALL dom_vvl_zgr
125      !
126      IF(lwxios) THEN
127! define variables in restart file when writing with XIOS
128         CALL iom_set_rstw_var_active('e3t_b')
129         CALL iom_set_rstw_var_active('e3t_n')
130         !                                           ! ----------------------- !
131         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
132            !                                        ! ----------------------- !
133            CALL iom_set_rstw_var_active('tilde_e3t_b')
134            CALL iom_set_rstw_var_active('tilde_e3t_n')
135         END IF
136         !                                           ! -------------!   
137         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
138            !                                        ! ------------ !
139            CALL iom_set_rstw_var_active('hdiv_lf')
140         ENDIF
141         !
142      ENDIF
143      !
144   END SUBROUTINE dom_vvl_init
145
146   SUBROUTINE dom_vvl_zgr
147      !!----------------------------------------------------------------------
148      !!                ***  ROUTINE dom_vvl_zgr  ***
149      !!                   
150      !! ** Purpose :  Interpolation of scale factors, update of depths
151      !!               and water column heights
152      !!
153      !! ** Method  :  - interpolate scale factors
154      !!
155      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
156      !!              - Regrid: e3(u/v)_n
157      !!                        e3(u/v)_b       
158      !!                        e3w_n           
159      !!                        e3(u/v)w_b     
160      !!                        e3(u/v)w_n     
161      !!                        gdept_n, gdepw_n and gde3w_n
162      !!              - h(t/u/v)_0
163      !!              - frq_rst_e3t and frq_rst_hdv
164      !!
165      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
166      !!----------------------------------------------------------------------
167      INTEGER ::   ji, jj, jk
168      INTEGER ::   ii0, ii1, ij0, ij1
169      REAL(wp)::   zcoef
170      !!----------------------------------------------------------------------
171      !
172      !                    !== Set vertical scale factors  ==!  (now and before)
173      !                                ! Horizontal interpolation of e3t
174      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
175      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
176      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V
177      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
178      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
179      !                                ! Vertical interpolation of e3t,u,v
180      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
181      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
182      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
183      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
184      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
185      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
186
187      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
188      e3t_a(:,:,:) = e3t_n(:,:,:)
189      e3u_a(:,:,:) = e3u_n(:,:,:)
190      e3v_a(:,:,:) = e3v_n(:,:,:)
191      !
192      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
193      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
194      gdepw_n(:,:,1) = 0.0_wp
195      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
196      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
197      gdepw_b(:,:,1) = 0.0_wp
198      DO jk = 2, jpk                               ! vertical sum
199         DO jj = 1,jpj
200            DO ji = 1,jpi
201               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
202               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
203               !                             ! 0.5 where jk = mikt     
204!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
205               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
206               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
207               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
208                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
209               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
210               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
211               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
212                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
213            END DO
214         END DO
215      END DO
216      !
217      !                    !==  thickness of the water column  !!   (ocean portion only)
218      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
219      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
220      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
221      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
222      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
223      DO jk = 2, jpkm1
224         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
225         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
226         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
227         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
228         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
229      END DO
230      !
231      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
232      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
233      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
234      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
235      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
236
237      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
238      IF( ln_vvl_ztilde ) THEN
239!!gm : idea: add here a READ in a file of custumized restoring frequency
240         !                                   ! Values in days provided via the namelist
241         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
242         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
243         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
244         !
245         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
246            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
247            frq_rst_hdv(:,:) = 1._wp / rdt
248         ENDIF
249         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
250            DO jj = 1, jpj
251               DO ji = 1, jpi
252!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
253                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
254                     ! values outside the equatorial band and transition zone (ztilde)
255                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
256                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
257                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
258                     ! values inside the equatorial band (ztilde as zstar)
259                     frq_rst_e3t(ji,jj) =  0.0_wp
260                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
261                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
262                     !                                      ! (linearly transition from z-tilde to z-star)
263                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
264                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
265                        &                                          * 180._wp / 3.5_wp ) )
266                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
267                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
268                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
269                        &                                          * 180._wp / 3.5_wp ) )
270                  ENDIF
271               END DO
272            END DO
273            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
274               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
275                  ii0 = 103   ;   ii1 = 111       
276                  ij0 = 128   ;   ij1 = 135   ;   
277                  frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
278                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
279               ENDIF
280            ENDIF
281         ENDIF
282      ENDIF
283      !
284   END SUBROUTINE dom_vvl_zgr
285
286   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
287      !!----------------------------------------------------------------------
288      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
289      !!                   
290      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
291      !!                 tranxt and dynspg routines
292      !!
293      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
294      !!               - z_tilde_case: after scale factor increment =
295      !!                                    high frequency part of horizontal divergence
296      !!                                  + retsoring towards the background grid
297      !!                                  + thickness difusion
298      !!                               Then repartition of ssh INCREMENT proportionnaly
299      !!                               to the "baroclinic" level thickness.
300      !!
301      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
302      !!               - tilde_e3t_a: after increment of vertical scale factor
303      !!                              in z_tilde case
304      !!               - e3(t/u/v)_a
305      !!
306      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
307      !!----------------------------------------------------------------------
308      INTEGER, INTENT( in )           ::   kt      ! time step
309      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
310      !
311      INTEGER                ::   ji, jj, jk            ! dummy loop indices
312      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
313      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
314      LOGICAL                ::   ll_do_bclinic         ! local logical
315      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
316      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
317      !!----------------------------------------------------------------------
318      !
319      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
320      !
321      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
322      !
323      IF( kt == nit000 ) THEN
324         IF(lwp) WRITE(numout,*)
325         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
326         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
327      ENDIF
328
329      ll_do_bclinic = .TRUE.
330      IF( PRESENT(kcall) ) THEN
331         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
332      ENDIF
333
334      ! ******************************* !
335      ! After acale factors at t-points !
336      ! ******************************* !
337      !                                                ! --------------------------------------------- !
338      !                                                ! z_star coordinate and barotropic z-tilde part !
339      !                                                ! --------------------------------------------- !
340      !
341      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
342      DO jk = 1, jpkm1
343         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
344         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
345      END DO
346      !
347      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
348         !                                                               ! ------baroclinic part------ !
349         ! I - initialization
350         ! ==================
351
352         ! 1 - barotropic divergence
353         ! -------------------------
354         zhdiv(:,:) = 0._wp
355         zht(:,:)   = 0._wp
356         DO jk = 1, jpkm1
357            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
358            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
359         END DO
360         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
361
362         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
363         ! --------------------------------------------------
364         IF( ln_vvl_ztilde ) THEN
365            IF( kt > nit000 ) THEN
366               DO jk = 1, jpkm1
367                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
368                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
369               END DO
370            ENDIF
371         ENDIF
372
373         ! II - after z_tilde increments of vertical scale factors
374         ! =======================================================
375         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
376
377         ! 1 - High frequency divergence term
378         ! ----------------------------------
379         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
380            DO jk = 1, jpkm1
381               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
382            END DO
383         ELSE                         ! layer case
384            DO jk = 1, jpkm1
385               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
386            END DO
387         ENDIF
388
389         ! 2 - Restoring term (z-tilde case only)
390         ! ------------------
391         IF( ln_vvl_ztilde ) THEN
392            DO jk = 1, jpk
393               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
394            END DO
395         ENDIF
396
397         ! 3 - Thickness diffusion term
398         ! ----------------------------
399         zwu(:,:) = 0._wp
400         zwv(:,:) = 0._wp
401         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
402            DO jj = 1, jpjm1
403               DO ji = 1, fs_jpim1   ! vector opt.
404                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
405                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
406                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
407                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
408                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
409                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
410               END DO
411            END DO
412         END DO
413         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
414            DO ji = 1, jpi
415               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
416               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
417            END DO
418         END DO
419         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
420            DO jj = 2, jpjm1
421               DO ji = fs_2, fs_jpim1   ! vector opt.
422                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
423                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
424                     &                                            ) * r1_e1e2t(ji,jj)
425               END DO
426            END DO
427         END DO
428         !                       ! d - thickness diffusion transport: boundary conditions
429         !                             (stored for tracer advction and continuity equation)
430         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
431
432         ! 4 - Time stepping of baroclinic scale factors
433         ! ---------------------------------------------
434         ! Leapfrog time stepping
435         ! ~~~~~~~~~~~~~~~~~~~~~~
436         IF( neuler == 0 .AND. kt == nit000 ) THEN
437            z2dt =  rdt
438         ELSE
439            z2dt = 2.0_wp * rdt
440         ENDIF
441         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
442         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
443
444         ! Maximum deformation control
445         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
446         ze3t(:,:,jpk) = 0._wp
447         DO jk = 1, jpkm1
448            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
449         END DO
450         z_tmax = MAXVAL( ze3t(:,:,:) )
451         CALL mpp_max( 'domvvl', z_tmax )                 ! max over the global domain
452         z_tmin = MINVAL( ze3t(:,:,:) )
453         CALL mpp_min( 'domvvl', z_tmin )                 ! min over the global domain
454         ! - ML - test: for the moment, stop simulation for too large e3_t variations
455         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
456            IF( lk_mpp ) THEN
457               CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max )
458               CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min )
459            ELSE
460               ijk_max = MAXLOC( ze3t(:,:,:) )
461               ijk_max(1) = ijk_max(1) + nimpp - 1
462               ijk_max(2) = ijk_max(2) + njmpp - 1
463               ijk_min = MINLOC( ze3t(:,:,:) )
464               ijk_min(1) = ijk_min(1) + nimpp - 1
465               ijk_min(2) = ijk_min(2) + njmpp - 1
466            ENDIF
467            IF (lwp) THEN
468               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
469               WRITE(numout, *) 'at i, j, k=', ijk_max
470               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
471               WRITE(numout, *) 'at i, j, k=', ijk_min           
472               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
473            ENDIF
474         ENDIF
475         ! - ML - end test
476         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
477         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
478         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
479
480         !
481         ! "tilda" change in the after scale factor
482         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
483         DO jk = 1, jpkm1
484            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
485         END DO
486         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
487         ! ==================================================================================
488         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
489         ! - ML - baroclinicity error should be better treated in the future
490         !        i.e. locally and not spread over the water column.
491         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
492         zht(:,:) = 0.
493         DO jk = 1, jpkm1
494            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
495         END DO
496         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
497         DO jk = 1, jpkm1
498            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
499         END DO
500
501      ENDIF
502
503      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
504      !                                           ! ---baroclinic part--------- !
505         DO jk = 1, jpkm1
506            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
507         END DO
508      ENDIF
509
510      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
511         !
512         IF( lwp ) WRITE(numout, *) 'kt =', kt
513         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
514            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
515            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
516            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
517         END IF
518         !
519         zht(:,:) = 0.0_wp
520         DO jk = 1, jpkm1
521            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
522         END DO
523         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
524         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
525         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
526         !
527         zht(:,:) = 0.0_wp
528         DO jk = 1, jpkm1
529            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
530         END DO
531         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
532         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
533         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
534         !
535         zht(:,:) = 0.0_wp
536         DO jk = 1, jpkm1
537            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
538         END DO
539         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
540         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
541         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
542         !
543         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
544         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
545         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
546         !
547         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
548         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
549         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
550         !
551         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
552         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
553         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
554      END IF
555
556      ! *********************************** !
557      ! After scale factors at u- v- points !
558      ! *********************************** !
559
560      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
561      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
562
563      ! *********************************** !
564      ! After depths at u- v points         !
565      ! *********************************** !
566
567      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
568      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
569      DO jk = 2, jpkm1
570         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
571         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
572      END DO
573      !                                        ! Inverse of the local depth
574!!gm BUG ?  don't understand the use of umask_i here .....
575      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
576      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
577      !
578      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
579      !
580   END SUBROUTINE dom_vvl_sf_nxt
581
582
583   SUBROUTINE dom_vvl_sf_swp( kt )
584      !!----------------------------------------------------------------------
585      !!                ***  ROUTINE dom_vvl_sf_swp  ***
586      !!                   
587      !! ** Purpose :  compute time filter and swap of scale factors
588      !!               compute all depths and related variables for next time step
589      !!               write outputs and restart file
590      !!
591      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
592      !!               - reconstruct scale factor at other grid points (interpolate)
593      !!               - recompute depths and water height fields
594      !!
595      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
596      !!               - Recompute:
597      !!                    e3(u/v)_b       
598      !!                    e3w_n           
599      !!                    e3(u/v)w_b     
600      !!                    e3(u/v)w_n     
601      !!                    gdept_n, gdepw_n  and gde3w_n
602      !!                    h(u/v) and h(u/v)r
603      !!
604      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
605      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
606      !!----------------------------------------------------------------------
607      INTEGER, INTENT( in ) ::   kt   ! time step
608      !
609      INTEGER  ::   ji, jj, jk   ! dummy loop indices
610      REAL(wp) ::   zcoef        ! local scalar
611      !!----------------------------------------------------------------------
612      !
613      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
614      !
615      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
616      !
617      IF( kt == nit000 )   THEN
618         IF(lwp) WRITE(numout,*)
619         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
620         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
621      ENDIF
622      !
623      ! Time filter and swap of scale factors
624      ! =====================================
625      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
626      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
627         IF( neuler == 0 .AND. kt == nit000 ) THEN
628            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
629         ELSE
630            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
631            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
632         ENDIF
633         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
634      ENDIF
635      gdept_b(:,:,:) = gdept_n(:,:,:)
636      gdepw_b(:,:,:) = gdepw_n(:,:,:)
637
638      e3t_n(:,:,:) = e3t_a(:,:,:)
639      e3u_n(:,:,:) = e3u_a(:,:,:)
640      e3v_n(:,:,:) = e3v_a(:,:,:)
641
642      ! Compute all missing vertical scale factor and depths
643      ! ====================================================
644      ! Horizontal scale factor interpolations
645      ! --------------------------------------
646      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
647      ! - JC - hu_b, hv_b, hur_b, hvr_b also
648     
649      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
650     
651      ! Vertical scale factor interpolations
652      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
653      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
654      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
655      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
656      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
657      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
658
659      ! t- and w- points depth (set the isf depth as it is in the initial step)
660      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
661      gdepw_n(:,:,1) = 0.0_wp
662      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
663      DO jk = 2, jpk
664         DO jj = 1,jpj
665            DO ji = 1,jpi
666              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
667                                                                 ! 1 for jk = mikt
668               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
669               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
670               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
671                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
672               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
673            END DO
674         END DO
675      END DO
676
677      ! Local depth and Inverse of the local depth of the water
678      ! -------------------------------------------------------
679      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
680      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
681      !
682      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
683      DO jk = 2, jpkm1
684         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
685      END DO
686
687      ! write restart file
688      ! ==================
689      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
690      !
691      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
692      !
693   END SUBROUTINE dom_vvl_sf_swp
694
695
696   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
697      !!---------------------------------------------------------------------
698      !!                  ***  ROUTINE dom_vvl__interpol  ***
699      !!
700      !! ** Purpose :   interpolate scale factors from one grid point to another
701      !!
702      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
703      !!                - horizontal interpolation: grid cell surface averaging
704      !!                - vertical interpolation: simple averaging
705      !!----------------------------------------------------------------------
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      !
711      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
712      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
713      !!----------------------------------------------------------------------
714      !
715      IF(ln_wd_il) THEN
716        zlnwd = 1.0_wp
717      ELSE
718        zlnwd = 0.0_wp
719      END IF
720      !
721      SELECT CASE ( pout )    !==  type of interpolation  ==!
722         !
723      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
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) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
728                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
729                     &                           + e1e2t(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         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
734         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
735         !
736      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
737         DO jk = 1, jpk
738            DO jj = 1, jpjm1
739               DO ji = 1, fs_jpim1   ! vector opt.
740                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
741                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
742                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
743               END DO
744            END DO
745         END DO
746         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
747         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
748         !
749      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
750         DO jk = 1, jpk
751            DO jj = 1, jpjm1
752               DO ji = 1, fs_jpim1   ! vector opt.
753                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
754                     &                       *    r1_e1e2f(ji,jj)                                                  &
755                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
756                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
757               END DO
758            END DO
759         END DO
760         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
761         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
762         !
763      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
764         !
765         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
766         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
767!!gm BUG? use here wmask in case of ISF ?  to be checked
768         DO jk = 2, jpk
769            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
770               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
771               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
772               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
773         END DO
774         !
775      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
776         !
777         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
778         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
779!!gm BUG? use here wumask in case of ISF ?  to be checked
780         DO jk = 2, jpk
781            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
782               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
783               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
784               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
785         END DO
786         !
787      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
788         !
789         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
790         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
791!!gm BUG? use here wvmask in case of ISF ?  to be checked
792         DO jk = 2, jpk
793            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
794               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
795               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
796               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
797         END DO
798      END SELECT
799      !
800   END SUBROUTINE dom_vvl_interpol
801
802
803   SUBROUTINE dom_vvl_rst( kt, cdrw )
804      !!---------------------------------------------------------------------
805      !!                   ***  ROUTINE dom_vvl_rst  ***
806      !!                     
807      !! ** Purpose :   Read or write VVL file in restart file
808      !!
809      !! ** Method  :   use of IOM library
810      !!                if the restart does not contain vertical scale factors,
811      !!                they are set to the _0 values
812      !!                if the restart does not contain vertical scale factors increments (z_tilde),
813      !!                they are set to 0.
814      !!----------------------------------------------------------------------
815      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
816      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
817      !
818      INTEGER ::   ji, jj, jk
819      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
820      !!----------------------------------------------------------------------
821      !
822      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
823         !                                   ! ===============
824         IF( ln_rstart ) THEN                   !* Read the restart file
825            CALL rst_read_open                  !  open the restart file if necessary
826            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
827            !
828            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
829            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
830            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
831            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
832            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
833            !                             ! --------- !
834            !                             ! all cases !
835            !                             ! --------- !
836            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
837               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
838               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
839               ! needed to restart if land processor not computed
840               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
841               WHERE ( tmask(:,:,:) == 0.0_wp ) 
842                  e3t_n(:,:,:) = e3t_0(:,:,:)
843                  e3t_b(:,:,:) = e3t_0(:,:,:)
844               END WHERE
845               IF( neuler == 0 ) THEN
846                  e3t_b(:,:,:) = e3t_n(:,:,:)
847               ENDIF
848            ELSE IF( id1 > 0 ) THEN
849               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
850               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
851               IF(lwp) write(numout,*) 'neuler is forced to 0'
852               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
853               e3t_n(:,:,:) = e3t_b(:,:,:)
854               neuler = 0
855            ELSE IF( id2 > 0 ) THEN
856               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
857               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
858               IF(lwp) write(numout,*) 'neuler is forced to 0'
859               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
860               e3t_b(:,:,:) = e3t_n(:,:,:)
861               neuler = 0
862            ELSE
863               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
864               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
865               IF(lwp) write(numout,*) 'neuler is forced to 0'
866               DO jk = 1, jpk
867                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
868                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
869                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
870               END DO
871               e3t_b(:,:,:) = e3t_n(:,:,:)
872               neuler = 0
873            ENDIF
874            !                             ! ----------- !
875            IF( ln_vvl_zstar ) THEN       ! z_star case !
876               !                          ! ----------- !
877               IF( MIN( id3, id4 ) > 0 ) THEN
878                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
879               ENDIF
880               !                          ! ----------------------- !
881            ELSE                          ! z_tilde and layer cases !
882               !                          ! ----------------------- !
883               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
884                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
885                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
886               ELSE                            ! one at least array is missing
887                  tilde_e3t_b(:,:,:) = 0.0_wp
888                  tilde_e3t_n(:,:,:) = 0.0_wp
889               ENDIF
890               !                          ! ------------ !
891               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
892                  !                       ! ------------ !
893                  IF( id5 > 0 ) THEN  ! required array exists
894                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
895                  ELSE                ! array is missing
896                     hdiv_lf(:,:,:) = 0.0_wp
897                  ENDIF
898               ENDIF
899            ENDIF
900            !
901         ELSE                                   !* Initialize at "rest"
902            !
903
904            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
905               !
906               IF( cn_cfg == 'wad' ) THEN
907                  ! Wetting and drying test case
908                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
909                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
910                  sshn (:,:)     = sshb(:,:)
911                  un   (:,:,:)   = ub  (:,:,:)
912                  vn   (:,:,:)   = vb  (:,:,:)
913               ELSE
914                  ! if not test case
915                  sshn(:,:) = -ssh_ref
916                  sshb(:,:) = -ssh_ref
917
918                  DO jj = 1, jpj
919                     DO ji = 1, jpi
920                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
921
922                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
923                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
924                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
925                        ENDIF
926                     ENDDO
927                  ENDDO
928               ENDIF !If test case else
929
930               ! Adjust vertical metrics for all wad
931               DO jk = 1, jpk
932                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
933                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
934                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
935               END DO
936               e3t_b(:,:,:) = e3t_n(:,:,:)
937
938               DO ji = 1, jpi
939                  DO jj = 1, jpj
940                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
941                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
942                     ENDIF
943                  END DO
944               END DO 
945               !
946            ELSE
947               !
948               ! Just to read set ssh in fact, called latter once vertical grid
949               ! is set up:
950!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
951!               !
952!               DO jk=1,jpk
953!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
954!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
955!               END DO
956!               e3t_n(:,:,:) = e3t_b(:,:,:)
957                sshn(:,:)=0._wp
958                e3t_n(:,:,:)=e3t_0(:,:,:)
959                e3t_b(:,:,:)=e3t_0(:,:,:)
960               !
961            END IF           ! end of ll_wd edits
962
963            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
964               tilde_e3t_b(:,:,:) = 0._wp
965               tilde_e3t_n(:,:,:) = 0._wp
966               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
967            END IF
968         ENDIF
969         !
970      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
971         !                                   ! ===================
972         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
973         IF( lwxios ) CALL iom_swap(      cwxios_context          )
974         !                                           ! --------- !
975         !                                           ! all cases !
976         !                                           ! --------- !
977         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
978         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
979         !                                           ! ----------------------- !
980         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
981            !                                        ! ----------------------- !
982            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
983            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
984         END IF
985         !                                           ! -------------!   
986         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
987            !                                        ! ------------ !
988            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
989         ENDIF
990         !
991         IF( lwxios ) CALL iom_swap(      cxios_context          )
992      ENDIF
993      !
994   END SUBROUTINE dom_vvl_rst
995
996
997   SUBROUTINE dom_vvl_ctl
998      !!---------------------------------------------------------------------
999      !!                  ***  ROUTINE dom_vvl_ctl  ***
1000      !!               
1001      !! ** Purpose :   Control the consistency between namelist options
1002      !!                for vertical coordinate
1003      !!----------------------------------------------------------------------
1004      INTEGER ::   ioptio, ios
1005      !!
1006      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1007         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
1008         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
1009      !!----------------------------------------------------------------------
1010      !
1011      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
1012      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1013901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
1014      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
1015      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1016902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
1017      IF(lwm) WRITE ( numond, nam_vvl )
1018      !
1019      IF(lwp) THEN                    ! Namelist print
1020         WRITE(numout,*)
1021         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1022         WRITE(numout,*) '~~~~~~~~~~~'
1023         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1024         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1025         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1026         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1027         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1028         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1029         WRITE(numout,*) '      !'
1030         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1031         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1032         IF( ln_vvl_ztilde_as_zstar ) THEN
1033            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1034            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1035            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1036            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1037            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1038            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1039         ELSE
1040            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1041            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1042         ENDIF
1043         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1044      ENDIF
1045      !
1046      ioptio = 0                      ! Parameter control
1047      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1048      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1049      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1050      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1051      !
1052      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1053      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1054      !
1055      IF(lwp) THEN                   ! Print the choice
1056         WRITE(numout,*)
1057         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1058         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1059         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1060         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1061      ENDIF
1062      !
1063#if defined key_agrif
1064      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1065#endif
1066      !
1067   END SUBROUTINE dom_vvl_ctl
1068
1069   !!======================================================================
1070END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.