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 branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6986

Last change on this file since 6986 was 6986, checked in by acc, 8 years ago

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