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/releases/r4.0/r4.0-HEAD/src/OCE/DOM – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DOM/domvvl.F90 @ 15563

Last change on this file since 15563 was 15563, checked in by jchanut, 3 years ago

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