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/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90 @ 11423

Last change on this file since 11423 was 11423, checked in by mathiot, 5 years ago

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