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 @ 15414

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

Transfer of a reusable section of subroutine dom_vvl_init into the new module subroutine dom_vvl_zgr (backport from trunk, ticket #2487)

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