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/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90 @ 10030

Last change on this file since 10030 was 10030, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.3 remove e3uw_$ e3vw_$, except e3.w_0 and use only after e3 in dyn/trazdf
  • Property svn:keywords set to Id
File size: 39.6 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  5.0  !  2018-07  (G. Madec)  Flux Form with Kinetic Energy conservation
7   !!                             ==>>>  here z* and s* only (no z-tilde)
8   
9   ! 1- remove   z-tilde          ==>>>  pure z-star (or s-star)
10   ! 2- remove   dom_vvl_interpol 
11   ! 3-
12   
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
17   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
18   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
19   !!   dom_vvl_rst      : read/write restart file
20   !!   dom_vvl_ctl      : Check the vvl options
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE phycst          ! physical constant
24   USE dom_oce         ! ocean space and time domain
25   USE sbc_oce         ! ocean surface boundary condition
26   USE wet_dry         ! wetting and drying
27   USE usrdef_istate   ! user defined initial state (wad only)
28   USE restart         ! ocean restart
29   !
30   USE in_out_manager  ! I/O manager
31   USE iom             ! I/O manager library
32   USE lib_mpp         ! distributed memory computing library
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE timing          ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC  dom_vvl_init       ! called by domain.F90
40   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
41   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
42   PUBLIC  ssh2e3_before      ! ...
43   PUBLIC  ssh2e3_now         ! ...
44
45   !                                                      !!* Namelist nam_vvl
46   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
51   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
52   !                                                       ! conservation: not used yet
53   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
54   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
55   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
56   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
57   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
58
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td      ! thickness diffusion transport
60   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf           ! low frequency part of hz divergence
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors
62   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a   ! baroclinic scale factors
63   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t       ! retoring period for scale factors
64   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv       ! retoring period for low freq. divergence
65
66
67!!gm add
68!!gm
69
70
71
72   !! * Substitutions
73#  include "vectopt_loop_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
76   !! $Id$
77   !! Software governed by the CeCILL licence     (./LICENSE)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   INTEGER FUNCTION dom_vvl_alloc()
82      !!----------------------------------------------------------------------
83      !!                ***  FUNCTION dom_vvl_alloc  ***
84      !!----------------------------------------------------------------------
85      !
86      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
87      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
88         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   &
89            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
90            &      STAT = dom_vvl_alloc        )
91         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
92         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
93         un_td = 0._wp
94         vn_td = 0._wp
95      ENDIF
96      IF( ln_vvl_ztilde ) THEN
97         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
98         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
99         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
100      ENDIF
101      !
102   END FUNCTION dom_vvl_alloc
103
104
105   SUBROUTINE dom_vvl_init
106      !!----------------------------------------------------------------------
107      !!                ***  ROUTINE dom_vvl_init  ***
108      !!                   
109      !! ** Purpose :  Initialization of all scale factors, depths
110      !!               and water column heights
111      !!
112      !! ** Method  :  - use restart file and/or initialize
113      !!               - interpolate scale factors
114      !!
115      !! ** Action  : - e3t_(n/b) and te3t_(n/b)
116      !!              - Regrid: e3(u/v)_n
117      !!                        e3(u/v)_b       
118      !!                        e3w_n           
119      !!                        e3(u/v)w_b     
120      !!                        e3(u/v)w_n     
121      !!                        gdept_n, gdepw_n and gde3w_n
122      !!              - h(t/u/v)_0
123      !!              - frq_rst_e3t and frq_rst_hdv
124      !!
125      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
126      !!----------------------------------------------------------------------
127      !!----------------------------------------------------------------------
128      !
129      IF(lwp) WRITE(numout,*)
130      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
131      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
132      !
133      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
134      !
135      !                    ! Allocate module arrays
136      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
137      !
138      !                    ! Read or initialize ssh(Nbb) & ssh(Nnn)
139      CALL dom_vvl_rst( nit000, 'READ' )
140      !
141      !                    !== Set of all other vertical mesh fields  ==!  (now and before)     
142      !
143      !                          !* BEFORE fields :
144      CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv
145      !                                !      e3t, e3u , e3v              (from 1 to jpkm1)
146      !                                !      e3w, gdept, gdepw           (from 1 to jpk  )
147      !
148      !                                ! set jpk level one to the e3._0 values
149      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3u_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
150      !
151      !                          !* NOW fields :
152      CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
153      !                                !      e3t, e3u , e3v, e3f         (from 1 to jpkm1)
154      !                                !      e3w, gdept, gdepw, gde3w    (from 1 to jpk  )
155      !
156      !                                ! set one for all last level to the e3._0 value
157      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3u_0(:,:,jpk)
158      e3f_n(:,:,jpk) = e3f_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
159      !
160      !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
161      e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)
162      e3w_a(:,:,:) = e3w_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
163     
164!!gm           
165!!gm        ===>>>>   below:  some issues to think about !!!
166!!gm 
167!!gm   fmask definition checked (0 or 1 nothing else)
168!!        in z-tilde or ALE    e3._0  should be the time varying fields step forward with an euler scheme
169!!gm   e3w_b  & gdept_b are not used   except its update in agrif 
170!!     and used to compute before slope of surface in dynldf_iso ==>>>  remove it !!!!
171!!         NB: in triads on TRA, gdept_n  is used !!!!   BUG?
172!!gm   e3f_n  almost not used  ===>>>>  verify whether it can be removed or not...
173!!       verify the use of wumask & wvmask   mau be replaced by a multiplication by umask(k)*umask(k+1)  ???
174!!
175!!gm  ISF case to be checked by Pierre Mathiot
176!!
177!!gm  setting of e3._a  for agrif....  TO BE CHECKED
178
179      !
180      IF(lwxios) THEN      ! define variables in restart file when writing with XIOS
181         CALL iom_set_rstw_var_active('e3t_b')
182         CALL iom_set_rstw_var_active('e3t_n')
183      ENDIF
184      !
185   END SUBROUTINE dom_vvl_init
186
187
188   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
189      !!----------------------------------------------------------------------
190      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
191      !!                   
192      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
193      !!                 tranxt and dynspg routines
194      !!
195      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
196      !!               - z_tilde_case: after scale factor increment =
197      !!                                    high frequency part of horizontal divergence
198      !!                                  + retsoring towards the background grid
199      !!                                  + thickness difusion
200      !!                               Then repartition of ssh INCREMENT proportionnaly
201      !!                               to the "baroclinic" level thickness.
202      !!
203      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
204      !!               - te3t_a: after increment of vertical scale factor
205      !!                              in z_tilde case
206      !!               - e3(t/u/v)_a
207      !!
208      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
209      !!
210      ! !  ref.   ! before  !   now   ! after  !
211      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
212      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
213      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
214      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
215      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
216      !
217      ! !  ref.   ! before  !   now   !
218      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
219      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
220      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
221      !
222      ! !  ref. ! before  !   now   !  after  !
223      !   ht_0            ,    ht_n             !: t-depth              [m]
224      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
225      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
226      !   hf_0                                  !: v-depth              [m]
227      ! r1_ht_0                                 !: inverse of u-depth [1/m]
228      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
229      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
230      ! r1_hf_0                                 !: inverse of v-depth [1/m]
231      !
232      !!----------------------------------------------------------------------
233      INTEGER, INTENT( in )           ::   kt      ! time step
234      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
235      !
236      INTEGER ::   ji, jj, jk   ! dummy loop indices
237      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
238      !!----------------------------------------------------------------------
239      !
240      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
241      !
242      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
243      !
244      IF( kt == nit000 ) THEN
245         IF(lwp) WRITE(numout,*)
246         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
247         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
248      ENDIF
249
250      !                             ! ------------------!
251      !                             ! z_star coordinate !     (and barotropic z-tilde part)
252      !                             ! ------------------!
253      !
254      !                                   !==  after ssh  ==!  (u- and v-points)
255      DO jj = 2, jpjm1
256         DO ji = 2, jpim1
257            zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
258            zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
259         END DO
260      END DO     
261      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
262      !
263      !                                   !==  after depths and its inverse  ==!
264         hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:)
265         hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:)
266      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
267      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
268      !
269      !                                   !==  after scale factors  ==!  (e3t , e3u , e3v)
270      zssht_h(:,:) = ssh    (:,:,Naa) * r1_ht_0(:,:)     ! t-point
271      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
272      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
273      DO jk = 1, jpkm1
274         e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
275         e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *      umask(:,:,jk)                     )
276         e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) *      vmask(:,:,jk)                     )
277         e3w_a(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
278      END DO
279      !
280      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
281      !
282   END SUBROUTINE dom_vvl_sf_nxt
283
284
285   SUBROUTINE dom_vvl_sf_swp( kt )
286      !!----------------------------------------------------------------------
287      !!                ***  ROUTINE dom_vvl_sf_swp  ***
288      !!                   
289      !! ** Purpose :  compute time filter and swap of scale factors
290      !!               compute all depths and related variables for next time step
291      !!               write outputs and restart file
292      !!
293      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
294      !!               - reconstruct scale factor at other grid points (interpolate)
295      !!               - recompute depths and water height fields
296      !!
297      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step
298      !!               - Recompute:
299      !!                    e3(u/v)_b       
300      !!                    e3w_n           
301      !!                    e3(u/v)w_b     
302      !!                    e3(u/v)w_n     
303      !!                    gdept_n, gdepw_n  and gde3w_n
304      !!                    h(u/v) and h(u/v)r
305      !!
306      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
307      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
308      !!----------------------------------------------------------------------
309      INTEGER, INTENT( in ) ::   kt   ! time step
310      !
311      INTEGER  ::   ji, jj, jk    ! dummy loop indices
312      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
313      !!----------------------------------------------------------------------
314      !
315      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
316      !
317      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
318      !
319      IF( kt == nit000 )   THEN
320         IF(lwp) WRITE(numout,*)
321         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
322         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
323      ENDIF
324      !
325      ! Swap and Compute all missing vertical scale factor and depths
326      ! =============================================================
327      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
328      ! - JC - hu_b, hv_b, hur_b, hvr_b also
329      !
330      ! - GM - to be updated :   e3f_n   ,  e3w_n  , e3w_b
331      !                          gdept_n , gdepw_n , gde3w_n
332      !                          ht_n
333      !        to be swap    :   hu_a , hv_a  , r1_hu_a , r1_hv_a
334      !
335      ! Local depth and Inverse of the local depth of the water
336      ! -------------------------------------------------------
337      !                       ! swap of depth and scale factors
338      !                       ! ===============================
339      DO jk = 1, jpkm1                          ! swap n--> a 
340         gdept_b(:,:,jk) = gdept_n(:,:,jk)         ! depth at t and w
341         gdepw_b(:,:,jk) = gdepw_n(:,:,jk)
342         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v, e3w
343         e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
344         e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
345         e3w_n  (:,:,jk) = e3w_a  (:,:,jk)
346      END DO
347      ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
348      !
349      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
350      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
351      !
352      !                    !==  before  ==!
353      !                                            !*  e3w_b
354      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
355      !
356      e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
357      DO jk = 2, jpk
358         e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
359      END DO
360      !
361      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw
362      !
363      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
364         DO jk = 1, jpkm1
365            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
366            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
367         END DO
368      ELSE                    ! no ISF cavities
369         DO jk = 1, jpkm1
370            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
371            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
372         END DO
373      ENDIF
374      !     
375      !                    !==   now    ==!
376      !                                            !* ssh at f-points
377      DO jj = 1, jpjm1
378         DO ji = 1, jpim1            ! start from 1 for f-point
379            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
380               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
381         END DO
382      END DO     
383      CALL lbc_lnk( zsshf_h(:,:),'F', 1._wp )     
384      !
385      !                                            !* e3f_n
386      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)     ! f-point
387      !
388      DO jk = 1, jpkm1
389          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
390      END DO     
391      !
392      !                                            !* gdept_n , gdepw_n , gde3w_n
393      zssht_h(:,:) = 1._wp + ssh(:,:,Nnn) * r1_ht_0(:,:) 
394      !
395      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
396         DO jk = 1, jpkm1
397            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
398            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
399            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
400         END DO
401      ELSE                    ! no ISF cavities
402         DO jk = 1, jpkm1
403            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
404            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
405            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
406         END DO
407      ENDIF
408      !
409      ! write restart file
410      ! ==================
411      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
412      !
413      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
414      !
415   END SUBROUTINE dom_vvl_sf_swp
416
417
418   SUBROUTINE dom_vvl_rst( kt, cdrw )
419      !!---------------------------------------------------------------------
420      !!                   ***  ROUTINE dom_vvl_rst  ***
421      !!                     
422      !! ** Purpose :   Read or write VVL file in restart file
423      !!
424      !! ** Method  :   use of IOM library
425      !!                if the restart does not contain vertical scale factors,
426      !!                they are set to the _0 values
427      !!                if the restart does not contain vertical scale factors increments (z_tilde),
428      !!                they are set to 0.
429      !!----------------------------------------------------------------------
430      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
431      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
432      !
433      INTEGER ::   ji, jj, jk
434      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
435      !!----------------------------------------------------------------------
436      !
437      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
438         !                                   ! ===============
439         IF( ln_rstart ) THEN                   !* Read the restart file
440            CALL rst_read_open                  !  open the restart file if necessary
441            !
442            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. )
443            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. )
444            !
445            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
446               IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files'
447               !
448!!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read)
449!!              so that it will work with land processor suppression
450!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
451!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
452!!gm
453               CALL iom_get( numror, jpdom_data, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
454               CALL iom_get( numror, jpdom_data, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
455!!gm end
456               IF( l_1st_euler )   ssh(:,:,Nbb) = ssh(:,:,Nnn)
457            ELSE IF( id1 > 0 ) THEN
458               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
459               IF(lwp) write(numout,*) '   set ssh(Nnn) = ssh(Nbb)  and force l_1st_euler = true'
460!!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
461               CALL iom_get( numror, jpdom_data, 'sshb', ssh(:,:,Nbb), ldxios = lrxios )
462               ssh(:,:,Nnn) = ssh(:,:,Nbb)
463               l_1st_euler = .TRUE.
464            ELSE IF( id2 > 0 ) THEN
465               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
466               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn)  and force l_1st_euler = true'
467               CALL iom_get( numror, jpdom_data, 'sshn', ssh(:,:,Nnn), ldxios = lrxios )
468               ssh(:,:,Nbb) = ssh(:,:,Nnn)
469               l_1st_euler = .TRUE.
470            ELSE
471               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
472               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn) = 0  and force l_1st_euler = true'
473               ssh(:,:,Nbb) = 0._wp
474               ssh(:,:,Nnn) = 0._wp
475               l_1st_euler = .TRUE.
476            ENDIF
477         ELSE                                   !* Initialize at "rest"
478            !
479            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
480               !
481               IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case
482                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
483                  tsn(:,:,:,:) = tsb(:,:,:,:)            ! set now values from to before ones
484                  ssh(:,:,Nnn) = ssh(:,:,Nbb)
485                  un (:,:,:)   = ub (:,:,:)
486                  vn (:,:,:)   = vb (:,:,:)
487               ELSE                                   ! Not the test case
488                  ssh(:,:,Nnn) = -ssh_ref
489                  ssh(:,:,Nbb) = -ssh_ref
490                  !
491                  DO jj = 1, jpj
492                     DO ji = 1, jpi
493                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth
494                           ssh(ji,jj,Nbb) = rn_wdmin1 - (ht_0(ji,jj) )
495                           ssh(ji,jj,Nnn) = rn_wdmin1 - (ht_0(ji,jj) )
496                           ssh(ji,jj,Naa) = rn_wdmin1 - (ht_0(ji,jj) )
497                        ENDIF
498                     END DO
499                  END DO
500               ENDIF     ! If test case else
501               !
502               DO ji = 1, jpi
503                  DO jj = 1, jpj
504                     IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN
505                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
506                     ENDIF
507                  END DO
508               END DO 
509               !
510            ELSE
511               !
512               ! Just to read set ssh in fact, called latter once vertical grid is set up:
513!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
514               ssh(:,:,Nnn) = 0._wp
515               ssh(:,:,Nbb) = 0._wp
516               !
517            END IF
518            !
519         ENDIF
520         !
521      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
522         !                                   ! ===================
523
524!!gm      DO NOTHING,   ssh(Nbb) and ssh(Nnn)  are written in restart.F90
525
526      ENDIF
527      !
528   END SUBROUTINE dom_vvl_rst
529
530
531   SUBROUTINE dom_vvl_ctl
532      !!---------------------------------------------------------------------
533      !!                  ***  ROUTINE dom_vvl_ctl  ***
534      !!               
535      !! ** Purpose :   Control the consistency between namelist options
536      !!                for vertical coordinate
537      !!----------------------------------------------------------------------
538      INTEGER ::   ioptio, ios
539      !!
540      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
541         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
542         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
543      !!----------------------------------------------------------------------
544      !
545      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
546      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
547901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
548      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
549      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
550902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
551      IF(lwm) WRITE ( numond, nam_vvl )
552      !
553      IF(lwp) THEN                    ! Namelist print
554         WRITE(numout,*)
555         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
556         WRITE(numout,*) '~~~~~~~~~~~'
557         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
558         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
559         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
560         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
561         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
562         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
563         WRITE(numout,*) '      !'
564         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
565         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
566         IF( ln_vvl_ztilde_as_zstar ) THEN
567            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
568            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
569            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
570            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
571            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
572            WRITE(numout,*) '                         rn_lf_cutoff   = 1/rn_Dt'
573         ELSE
574            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
575            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
576         ENDIF
577         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
578      ENDIF
579      !
580
581!!gm
582      IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar )   CALL ctl_stop( 'z-tilde NOT available in this branch' )
583!!gm
584
585      !
586      ioptio = 0                      ! Parameter control
587      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
588      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
589      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
590      IF( ln_vvl_layer           )   ioptio = ioptio + 1
591      !
592      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
593      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
594      !
595      IF(lwp) THEN                   ! Print the choice
596         WRITE(numout,*)
597         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
598         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
599         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
600         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
601      ENDIF
602      !
603#if defined key_agrif
604      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
605#endif
606      !
607   END SUBROUTINE dom_vvl_ctl
608
609   
610   SUBROUTINE ssh2e3_now
611      !!----------------------------------------------------------------------
612      !!                  ***  ROUTINE ssh2e3_now  ***
613      ! !  ref.   ! before  !   now   ! after  !
614      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
615      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
616      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
617      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
618      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
619      !
620      ! !  ref.   ! before  !   now   !
621      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
622      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
623      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
624      !
625      ! !  ref. ! before  !   now   !  after  !
626      !   ht_0            ,    ht_n             !: t-depth              [m]
627      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
628      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
629      !   hf_0                                  !: v-depth              [m]
630      ! r1_ht_0                                 !: inverse of u-depth [1/m]
631      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
632      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
633      ! r1_hf_0                                 !: inverse of v-depth [1/m]
634      !
635      !!----------------------------------------------------------------------
636      INTEGER ::   ji, jj, jk
637      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
638      !!----------------------------------------------------------------------
639      !
640      !                             !==  ssh at u- and v-points  ==!
641      DO jj = 1, jpjm1
642         DO ji = 1, jpim1                 ! start from 1 due to f-point
643            zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
644            zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
645            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
646               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
647         END DO
648      END DO     
649      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )
650      !
651      !                             !==  ht, hu and hv  == !   (and their inverse)
652      ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
653      hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
654      hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:)
655      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
656      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
657      !     
658      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
659      zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
660      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
661      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
662      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
663      !
664      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ,  e3w  ==!
665      DO jk = 1, jpkm1
666          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
667          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
668          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
669          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
670      END DO
671      e3w_n(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
672      DO jk = 2, jpk
673         e3w_n(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
674      END DO
675      !
676      !                             !== depth of t- and w-points  ==!
677      !
678      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
679      !
680      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
681         DO jk = 1, jpk
682            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
683            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
684            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
685         END DO
686      ELSE                    ! no ISF cavities
687         DO jk = 1, jpk
688            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
689            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
690            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
691         END DO
692      ENDIF
693      !
694   END SUBROUTINE ssh2e3_now
695   
696   
697   SUBROUTINE ssh2e3_before
698      !!----------------------------------------------------------------------
699      !!                  ***  ROUTINE ssh2e3_before  ***
700      ! !  ref.   ! before  !   now   ! after  !
701      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
702      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
703      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
704      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
705      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
706      !
707      ! !  ref.   ! before  !   now   !
708      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
709      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
710      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
711      !
712      ! !  ref. ! before  !   now   !  after  !
713      !   ht_0            ,    ht_n             !: t-depth              [m]
714      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
715      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
716      !   hf_0                                  !: v-depth              [m]
717      ! r1_ht_0                                 !: inverse of u-depth [1/m]
718      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
719      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
720      ! r1_hf_0                                 !: inverse of v-depth [1/m]
721      !
722      !!----------------------------------------------------------------------
723      INTEGER ::   ji, jj, jk
724      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
725      !!----------------------------------------------------------------------
726      !
727      !                             !==  ssh at u- and v-points  ==!
728      DO jj = 2, jpjm1
729         DO ji = 2, jpim1
730            zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
731            zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
732         END DO
733      END DO     
734      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
735      !
736      !                             !==  ht, hu and hv  == !   (and their inverse)
737         hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:)
738         hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:)
739      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
740      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
741      !
742      !     
743      !                             !==  ssh / h  factor at t-, u- ,v-points  ==!
744      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
745      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
746      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
747      !
748      !                             !==  e3t  ,  e3u  ,  e3v ,  e3w  ==!
749      DO jk = 1, jpkm1
750          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
751          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
752          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
753      END DO
754      e3w_b(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
755      DO jk = 2, jpk
756         e3w_b(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
757      END DO
758      !   
759      !                             !== depth of t- and w-points  ==!
760      !
761      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
762      !
763      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
764         DO jk = 1, jpk
765            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
766            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
767         END DO
768      ELSE                    ! no ISF cavities
769         DO jk = 1, jpk
770            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
771            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
772         END DO
773      ENDIF
774      !
775   END SUBROUTINE ssh2e3_before
776   
777   !!======================================================================
778END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.