New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90 @ 12482

Last change on this file since 12482 was 12482, checked in by techene, 4 years ago

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