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_RK3.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_RK3.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: 40.0 KB
Line 
1MODULE domvvl_RK3
2   !!======================================================================
3   !!                       ***  MODULE domvvl_RK3   ***
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)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
158      e3f_n(:,:,jpk) = e3f_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      !    e3uw_0                                !: uw-vert. scale factor [m]
216      !    e3vw_0                                !: vw-vert. scale factor [m]
217      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
218      !
219      ! !  ref.   ! before  !   now   !
220      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
221      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
222      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
223      !
224      ! !  ref. ! before  !   now   !  after  !
225      !   ht_0            ,    ht_n             !: t-depth              [m]
226      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
227      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
228      !   hf_0                                  !: v-depth              [m]
229      ! r1_ht_0                                 !: inverse of u-depth [1/m]
230      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
231      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
232      ! r1_hf_0                                 !: inverse of v-depth [1/m]
233      !
234      !!----------------------------------------------------------------------
235      INTEGER, INTENT( in )           ::   kt      ! time step
236      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
237      !
238      INTEGER ::   ji, jj, jk   ! dummy loop indices
239      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
240      !!----------------------------------------------------------------------
241      !
242      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
243      !
244      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
245      !
246      IF( kt == nit000 ) THEN
247         IF(lwp) WRITE(numout,*)
248         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
249         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
250      ENDIF
251
252      !                             ! ------------------!
253      !                             ! z_star coordinate !     (and barotropic z-tilde part)
254      !                             ! ------------------!
255      !
256      !                                   !==  after ssh  ==!  (u- and v-points)
257      DO jj = 2, jpjm1
258         DO ji = 2, jpim1
259            zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
260            zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
261         END DO
262      END DO     
263      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
264      !
265      !                                   !==  after depths and its inverse  ==!
266         hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:)
267         hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:)
268      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
269      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
270      !
271      !                                   !==  after scale factors  ==!  (e3t , e3u , e3v)
272      zssht_h(:,:) = ssh    (:,:,Naa) * r1_ht_0(:,:)     ! t-point
273      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
274      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
275      DO jk = 1, jpkm1
276         e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
277         e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *      umask(:,:,jk)                     )
278         e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) *      vmask(:,:,jk)                     )
279         e3w_a(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
280      END DO
281      !
282      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
283      !
284   END SUBROUTINE dom_vvl_sf_nxt
285
286
287   SUBROUTINE dom_vvl_sf_swp( kt )
288      !!----------------------------------------------------------------------
289      !!                ***  ROUTINE dom_vvl_sf_swp  ***
290      !!                   
291      !! ** Purpose :  compute time filter and swap of scale factors
292      !!               compute all depths and related variables for next time step
293      !!               write outputs and restart file
294      !!
295      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
296      !!               - reconstruct scale factor at other grid points (interpolate)
297      !!               - recompute depths and water height fields
298      !!
299      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step
300      !!               - Recompute:
301      !!                    e3(u/v)_b       
302      !!                    e3w_n           
303      !!                    e3(u/v)w_b     
304      !!                    e3(u/v)w_n     
305      !!                    gdept_n, gdepw_n  and gde3w_n
306      !!                    h(u/v) and h(u/v)r
307      !!
308      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
309      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
310      !!----------------------------------------------------------------------
311      INTEGER, INTENT( in ) ::   kt   ! time step
312      !
313      INTEGER  ::   ji, jj, jk    ! dummy loop indices
314      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
315      !!----------------------------------------------------------------------
316      !
317      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
318      !
319      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
320      !
321      IF( kt == nit000 )   THEN
322         IF(lwp) WRITE(numout,*)
323         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
324         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
325      ENDIF
326      !
327      ! Swap and Compute all missing vertical scale factor and depths
328      ! =============================================================
329      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
330      ! - JC - hu_b, hv_b, hur_b, hvr_b also
331      !
332      ! - GM - to be updated :   e3f_n   ,  e3w_n  , e3w_b
333      !                          gdept_n , gdepw_n , gde3w_n
334      !                          ht_n
335      !        to be swap    :   hu_a , hv_a  , r1_hu_a , r1_hv_a
336      !
337      ! Local depth and Inverse of the local depth of the water
338      ! -------------------------------------------------------
339      !                       ! swap of depth and scale factors
340      !                       ! ===============================
341      DO jk = 1, jpkm1                          ! swap n--> a 
342         gdept_b(:,:,jk) = gdept_n(:,:,jk)         ! depth at t and w
343         gdepw_b(:,:,jk) = gdepw_n(:,:,jk)
344         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v, e3w
345         e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
346         e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
347         e3w_n  (:,:,jk) = e3w_a  (:,:,jk)
348      END DO
349      ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
350      !
351      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
352      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
353      !
354      !                    !==  before  ==!
355      !                                            !*  e3w_b
356      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
357      !
358      e3w_b(:,:,1) =  e3w_0(:,:,1) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,1) )
359      DO jk = 2, jpk
360         e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) ,  tmask(:,:,jk) ) )
361      END DO
362      !
363      zssht_h(:,:) = 1._wp + zssht_h(:,:)          !* gdept , gdepw
364      !
365      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
366         DO jk = 1, jpkm1
367            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
368            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
369         END DO
370      ELSE                    ! no ISF cavities
371         DO jk = 1, jpkm1
372            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
373            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
374         END DO
375      ENDIF
376      !     
377      !                    !==   now    ==!
378      !                                            !* ssh at f-points
379      DO jj = 1, jpjm1
380         DO ji = 1, jpim1            ! start from 1 for f-point
381            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
382               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
383         END DO
384      END DO     
385      CALL lbc_lnk( zsshf_h(:,:),'F', 1._wp )     
386      !
387      !                                            !* e3f_n
388      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)     ! f-point
389      !
390      DO jk = 1, jpkm1
391          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
392      END DO     
393      !
394      !                                            !* gdept_n , gdepw_n , gde3w_n
395      zssht_h(:,:) = 1._wp + ssh(:,:,Nnn) * r1_ht_0(:,:) 
396      !
397      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
398         DO jk = 1, jpkm1
399            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
400            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
401            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
402         END DO
403      ELSE                    ! no ISF cavities
404         DO jk = 1, jpkm1
405            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
406            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
407            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
408         END DO
409      ENDIF
410      !
411      ! write restart file
412      ! ==================
413      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
414      !
415      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
416      !
417   END SUBROUTINE dom_vvl_sf_swp
418
419
420   SUBROUTINE dom_vvl_rst( kt, cdrw )
421      !!---------------------------------------------------------------------
422      !!                   ***  ROUTINE dom_vvl_rst  ***
423      !!                     
424      !! ** Purpose :   Read or write VVL file in restart file
425      !!
426      !! ** Method  :   use of IOM library
427      !!                if the restart does not contain vertical scale factors,
428      !!                they are set to the _0 values
429      !!                if the restart does not contain vertical scale factors increments (z_tilde),
430      !!                they are set to 0.
431      !!----------------------------------------------------------------------
432      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
433      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
434      !
435      INTEGER ::   ji, jj, jk
436      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
437      !!----------------------------------------------------------------------
438      !
439      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
440         !                                   ! ===============
441         IF( ln_rstart ) THEN                   !* Read the restart file
442            CALL rst_read_open                  !  open the restart file if necessary
443            !
444            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. )
445            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. )
446            !
447            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
448               IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files'
449               !
450!!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read)
451!!              so that it will work with land processor suppression
452!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
453!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
454!!gm
455               CALL iom_get( numror, jpdom_data, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
456               CALL iom_get( numror, jpdom_data, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
457!!gm end
458               IF( l_1st_euler )   ssh(:,:,Nbb) = ssh(:,:,Nnn)
459            ELSE IF( id1 > 0 ) THEN
460               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
461               IF(lwp) write(numout,*) '   set ssh(Nnn) = ssh(Nbb)  and force l_1st_euler = true'
462!!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
463               CALL iom_get( numror, jpdom_data, 'sshb', ssh(:,:,Nbb), ldxios = lrxios )
464               ssh(:,:,Nnn) = ssh(:,:,Nbb)
465               l_1st_euler = .TRUE.
466            ELSE IF( id2 > 0 ) THEN
467               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
468               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn)  and force l_1st_euler = true'
469               CALL iom_get( numror, jpdom_data, 'sshn', ssh(:,:,Nnn), ldxios = lrxios )
470               ssh(:,:,Nbb) = ssh(:,:,Nnn)
471               l_1st_euler = .TRUE.
472            ELSE
473               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
474               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn) = 0  and force l_1st_euler = true'
475               ssh(:,:,Nbb) = 0._wp
476               ssh(:,:,Nnn) = 0._wp
477               l_1st_euler = .TRUE.
478            ENDIF
479         ELSE                                   !* Initialize at "rest"
480            !
481            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
482               !
483               IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case
484                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
485                  tsn(:,:,:,:) = tsb(:,:,:,:)            ! set now values from to before ones
486                  ssh(:,:,Nnn) = ssh(:,:,Nbb)
487                  un (:,:,:)   = ub (:,:,:)
488                  vn (:,:,:)   = vb (:,:,:)
489               ELSE                                   ! Not the test case
490                  ssh(:,:,Nnn) = -ssh_ref
491                  ssh(:,:,Nbb) = -ssh_ref
492                  !
493                  DO jj = 1, jpj
494                     DO ji = 1, jpi
495                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth
496                           ssh(ji,jj,Nbb) = rn_wdmin1 - (ht_0(ji,jj) )
497                           ssh(ji,jj,Nnn) = rn_wdmin1 - (ht_0(ji,jj) )
498                           ssh(ji,jj,Naa) = rn_wdmin1 - (ht_0(ji,jj) )
499                        ENDIF
500                     END DO
501                  END DO
502               ENDIF     ! If test case else
503               !
504               DO ji = 1, jpi
505                  DO jj = 1, jpj
506                     IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN
507                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
508                     ENDIF
509                  END DO
510               END DO 
511               !
512            ELSE
513               !
514               ! Just to read set ssh in fact, called latter once vertical grid is set up:
515!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
516               ssh(:,:,Nnn) = 0._wp
517               ssh(:,:,Nbb) = 0._wp
518               !
519            END IF
520            !
521         ENDIF
522         !
523      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
524         !                                   ! ===================
525
526!!gm      DO NOTHING,   ssh(Nbb) and ssh(Nnn)  are written in restart.F90
527
528      ENDIF
529      !
530   END SUBROUTINE dom_vvl_rst
531
532
533   SUBROUTINE dom_vvl_ctl
534      !!---------------------------------------------------------------------
535      !!                  ***  ROUTINE dom_vvl_ctl  ***
536      !!               
537      !! ** Purpose :   Control the consistency between namelist options
538      !!                for vertical coordinate
539      !!----------------------------------------------------------------------
540      INTEGER ::   ioptio, ios
541      !!
542      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
543         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
544         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
545      !!----------------------------------------------------------------------
546      !
547      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
548      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
549901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
550      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
551      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
552902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
553      IF(lwm) WRITE ( numond, nam_vvl )
554      !
555      IF(lwp) THEN                    ! Namelist print
556         WRITE(numout,*)
557         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
558         WRITE(numout,*) '~~~~~~~~~~~'
559         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
560         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
561         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
562         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
563         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
564         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
565         WRITE(numout,*) '      !'
566         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
567         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
568         IF( ln_vvl_ztilde_as_zstar ) THEN
569            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
570            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
571            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
572            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
573            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
574            WRITE(numout,*) '                         rn_lf_cutoff   = 1/rn_Dt'
575         ELSE
576            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
577            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
578         ENDIF
579         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
580      ENDIF
581      !
582
583!!gm
584      IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar )   CALL ctl_stop( 'z-tilde NOT available in this branch' )
585!!gm
586
587      !
588      ioptio = 0                      ! Parameter control
589      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
590      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
591      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
592      IF( ln_vvl_layer           )   ioptio = ioptio + 1
593      !
594      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
595      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
596      !
597      IF(lwp) THEN                   ! Print the choice
598         WRITE(numout,*)
599         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
600         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
601         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
602         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
603      ENDIF
604      !
605#if defined key_agrif
606      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
607#endif
608      !
609   END SUBROUTINE dom_vvl_ctl
610
611   
612   SUBROUTINE ssh2e3_now
613      !!----------------------------------------------------------------------
614      !!                  ***  ROUTINE ssh2e3_now  ***
615      ! !  ref.   ! before  !   now   ! after  !
616      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
617      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
618      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
619      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
620      !    e3uw_0                                !: uw-vert. scale factor [m]
621      !    e3vw_0                                !: vw-vert. scale factor [m]
622      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
623      !
624      ! !  ref.   ! before  !   now   !
625      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
626      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
627      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
628      !
629      ! !  ref. ! before  !   now   !  after  !
630      !   ht_0            ,    ht_n             !: t-depth              [m]
631      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
632      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
633      !   hf_0                                  !: v-depth              [m]
634      ! r1_ht_0                                 !: inverse of u-depth [1/m]
635      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
636      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
637      ! r1_hf_0                                 !: inverse of v-depth [1/m]
638      !
639      !!----------------------------------------------------------------------
640      INTEGER ::   ji, jj, jk
641      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
642      !!----------------------------------------------------------------------
643      !
644      !                             !==  ssh at u- and v-points  ==!
645      DO jj = 1, jpjm1
646         DO ji = 1, jpim1                 ! start from 1 due to f-point
647            zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
648            zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
649            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
650               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
651         END DO
652      END DO     
653      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )
654      !
655      !                             !==  ht, hu and hv  == !   (and their inverse)
656      ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
657      hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
658      hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:)
659      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
660      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
661      !     
662      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
663      zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
664      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
665      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
666      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
667      !
668      !                             !==  e3t  ,  e3u  ,  e3v  ,  e3f  ,  e3w  ==!
669      DO jk = 1, jpkm1
670          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
671          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
672          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
673          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
674      END DO
675      e3w_n(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
676      DO jk = 2, jpk
677         e3w_n(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
678      END DO
679      !
680      !                             !== depth of t- and w-points  ==!
681      !
682      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
683      !
684      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
685         DO jk = 1, jpk
686            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
687            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
688            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
689         END DO
690      ELSE                    ! no ISF cavities
691         DO jk = 1, jpk
692            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
693            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
694            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
695         END DO
696      ENDIF
697      !
698   END SUBROUTINE ssh2e3_now
699   
700   
701   SUBROUTINE ssh2e3_before
702      !!----------------------------------------------------------------------
703      !!                  ***  ROUTINE ssh2e3_before  ***
704      ! !  ref.   ! before  !   now   ! after  !
705      !     e3t_0 ,   e3t_b ,   e3t_n ,  e3t_a   !: t- vert. scale factor [m]
706      !     e3u_0 ,   e3u_b ,   e3u_n ,  e3u_a   !: u- vert. scale factor [m]
707      !     e3v_0 ,   e3v_b ,   e3v_n ,  e3v_a   !: v- vert. scale factor [m]
708      !     e3w_0 ,   e3w_b ,   e3w_n ,  e3w_a   !: w- vert. scale factor [m]
709      !    e3uw_0                                !: uw-vert. scale factor [m]
710      !    e3vw_0                                !: vw-vert. scale factor [m]
711      !     e3f_0           ,   e3f_n            !: f- vert. scale factor [m]
712      !
713      ! !  ref.   ! before  !   now   !
714      !   gdept_0 , gdept_b , gdept_n   !: t- depth              [m]
715      !   gdepw_0 , gdepw_b , gdepw_n   !: w- depth              [m]
716      !   gde3w_0           , gde3w_n   !: w- depth (sum of e3w) [m]
717      !
718      ! !  ref. ! before  !   now   !  after  !
719      !   ht_0            ,    ht_n             !: t-depth              [m]
720      !   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m]
721      !   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m]
722      !   hf_0                                  !: v-depth              [m]
723      ! r1_ht_0                                 !: inverse of u-depth [1/m]
724      ! r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m]
725      ! r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m]
726      ! r1_hf_0                                 !: inverse of v-depth [1/m]
727      !
728      !!----------------------------------------------------------------------
729      INTEGER ::   ji, jj, jk
730      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
731      !!----------------------------------------------------------------------
732      !
733      !                             !==  ssh at u- and v-points  ==!
734      DO jj = 2, jpjm1
735         DO ji = 2, jpim1
736            zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
737            zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
738         END DO
739      END DO     
740      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
741      !
742      !                             !==  ht, hu and hv  == !   (and their inverse)
743         hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:)
744         hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:)
745      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
746      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
747      !
748      !     
749      !                             !==  ssh / h  factor at t-, u- ,v-points  ==!
750      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
751      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
752      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
753      !
754      !                             !==  e3t  ,  e3u  ,  e3v ,  e3w  ==!
755      DO jk = 1, jpkm1
756          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
757          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
758          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
759      END DO
760      e3w_b(:,:,1) = e3w_0(:,:,1) * (  1._wp + zssht_h(:,:) *  tmask(:,:,1)  )
761      DO jk = 2, jpk
762         e3w_b(:,:,jk) = e3w_0(:,:,jk) * (  1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk-1) , tmask(:,:,jk) )  )
763      END DO
764      !   
765      !                             !== depth of t- and w-points  ==!
766      !
767      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
768      !
769      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
770         DO jk = 1, jpk
771            gdept_b(:,:,jk) = ( gdept_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
772            gdepw_b(:,:,jk) = ( gdepw_0(:,:,jk) - ht_isf(:,:) ) * zssht_h(:,:) + ht_isf(:,:)
773         END DO
774      ELSE                    ! no ISF cavities
775         DO jk = 1, jpk
776            gdept_b(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
777            gdepw_b(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
778         END DO
779      ENDIF
780      !
781   END SUBROUTINE ssh2e3_before
782   
783   !!======================================================================
784END MODULE domvvl_RK3
Note: See TracBrowser for help on using the repository browser.