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

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

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