source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 54.5 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 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/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( cn_cfg == "orca" .AND. nn_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), DIMENSION(jpi,jpj,jpk) ::   ze3t
278      REAL(wp), DIMENSION(jpi,jpj) ::   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
286      IF( kt == nit000 ) THEN
287         IF(lwp) WRITE(numout,*)
288         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
289         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
290      ENDIF
291
292      ll_do_bclinic = .TRUE.
293      IF( PRESENT(kcall) ) THEN
294         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
295      ENDIF
296
297      ! ******************************* !
298      ! After acale factors at t-points !
299      ! ******************************* !
300      !                                                ! --------------------------------------------- !
301      !                                                ! z_star coordinate and barotropic z-tilde part !
302      !                                                ! --------------------------------------------- !
303      !
304      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
305      DO jk = 1, jpkm1
306         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
307         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
308      END DO
309      !
310      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
311         !                                                            ! ------baroclinic part------ !
312         ! I - initialization
313         ! ==================
314
315         ! 1 - barotropic divergence
316         ! -------------------------
317         zhdiv(:,:) = 0._wp
318         zht(:,:)   = 0._wp
319         DO jk = 1, jpkm1
320            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
321            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
322         END DO
323         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
324
325         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
326         ! --------------------------------------------------
327         IF( ln_vvl_ztilde ) THEN
328            IF( kt > nit000 ) THEN
329               DO jk = 1, jpkm1
330                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
331                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
332               END DO
333            ENDIF
334         ENDIF
335
336         ! II - after z_tilde increments of vertical scale factors
337         ! =======================================================
338         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
339
340         ! 1 - High frequency divergence term
341         ! ----------------------------------
342         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
343            DO jk = 1, jpkm1
344               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
345            END DO
346         ELSE                         ! layer case
347            DO jk = 1, jpkm1
348               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
349            END DO
350         ENDIF
351
352         ! 2 - Restoring term (z-tilde case only)
353         ! ------------------
354         IF( ln_vvl_ztilde ) THEN
355            DO jk = 1, jpk
356               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
357            END DO
358         ENDIF
359
360         ! 3 - Thickness diffusion term
361         ! ----------------------------
362         zwu(:,:) = 0._wp
363         zwv(:,:) = 0._wp
364         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
365            DO jj = 1, jpjm1
366               DO ji = 1, fs_jpim1   ! vector opt.
367                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
368                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
369                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
370                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
371                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
372                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
373               END DO
374            END DO
375         END DO
376         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
377            DO ji = 1, jpi
378               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
379               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
380            END DO
381         END DO
382         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
383            DO jj = 2, jpjm1
384               DO ji = fs_2, fs_jpim1   ! vector opt.
385                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
386                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
387                     &                                            ) * r1_e1e2t(ji,jj)
388               END DO
389            END DO
390         END DO
391         !                       ! d - thickness diffusion transport: boundary conditions
392         !                             (stored for tracer advction and continuity equation)
393         CALL lbc_lnk( un_td , 'U' , -1._wp)
394         CALL lbc_lnk( vn_td , 'V' , -1._wp)
395
396         ! 4 - Time stepping of baroclinic scale factors
397         ! ---------------------------------------------
398         ! Leapfrog time stepping
399         ! ~~~~~~~~~~~~~~~~~~~~~~
400         IF( neuler == 0 .AND. kt == nit000 ) THEN
401            z2dt =  rdt
402         ELSE
403            z2dt = 2.0_wp * rdt
404         ENDIF
405         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
406         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
407
408         ! Maximum deformation control
409         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
410         ze3t(:,:,jpk) = 0._wp
411         DO jk = 1, jpkm1
412            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
413         END DO
414         z_tmax = MAXVAL( ze3t(:,:,:) )
415         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
416         z_tmin = MINVAL( ze3t(:,:,:) )
417         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
418         ! - ML - test: for the moment, stop simulation for too large e3_t variations
419         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
420            IF( lk_mpp ) THEN
421               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
422               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
423            ELSE
424               ijk_max = MAXLOC( ze3t(:,:,:) )
425               ijk_max(1) = ijk_max(1) + nimpp - 1
426               ijk_max(2) = ijk_max(2) + njmpp - 1
427               ijk_min = MINLOC( ze3t(:,:,:) )
428               ijk_min(1) = ijk_min(1) + nimpp - 1
429               ijk_min(2) = ijk_min(2) + njmpp - 1
430            ENDIF
431            IF (lwp) THEN
432               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
433               WRITE(numout, *) 'at i, j, k=', ijk_max
434               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
435               WRITE(numout, *) 'at i, j, k=', ijk_min           
436               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
437            ENDIF
438         ENDIF
439         ! - ML - end test
440         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
441         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
442         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
443
444         !
445         ! "tilda" change in the after scale factor
446         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
447         DO jk = 1, jpkm1
448            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
449         END DO
450         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
451         ! ==================================================================================
452         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
453         ! - ML - baroclinicity error should be better treated in the future
454         !        i.e. locally and not spread over the water column.
455         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
456         zht(:,:) = 0.
457         DO jk = 1, jpkm1
458            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
459         END DO
460         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
461         DO jk = 1, jpkm1
462            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
463         END DO
464
465      ENDIF
466
467      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
468      !                                           ! ---baroclinic part--------- !
469         DO jk = 1, jpkm1
470            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
471         END DO
472      ENDIF
473
474      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
475         !
476         IF( lwp ) WRITE(numout, *) 'kt =', kt
477         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
478            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
479            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
480            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
481         END IF
482         !
483         zht(:,:) = 0.0_wp
484         DO jk = 1, jpkm1
485            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
486         END DO
487         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
488         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
489         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
490         !
491         zht(:,:) = 0.0_wp
492         DO jk = 1, jpkm1
493            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
494         END DO
495         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
496         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
497         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
498         !
499         zht(:,:) = 0.0_wp
500         DO jk = 1, jpkm1
501            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
502         END DO
503         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
504         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
505         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
506         !
507         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
508         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
509         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
510         !
511         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
512         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
513         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
514         !
515         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
516         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
517         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
518      END IF
519
520      ! *********************************** !
521      ! After scale factors at u- v- points !
522      ! *********************************** !
523
524      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
525      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
526
527      ! *********************************** !
528      ! After depths at u- v points         !
529      ! *********************************** !
530
531      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
532      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
533      DO jk = 2, jpkm1
534         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
535         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
536      END DO
537      !                                        ! Inverse of the local depth
538!!gm BUG ?  don't understand the use of umask_i here .....
539      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
540      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
541      !
542      !
543      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
544      !
545   END SUBROUTINE dom_vvl_sf_nxt
546
547
548   SUBROUTINE dom_vvl_sf_swp( kt )
549      !!----------------------------------------------------------------------
550      !!                ***  ROUTINE dom_vvl_sf_swp  ***
551      !!                   
552      !! ** Purpose :  compute time filter and swap of scale factors
553      !!               compute all depths and related variables for next time step
554      !!               write outputs and restart file
555      !!
556      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
557      !!               - reconstruct scale factor at other grid points (interpolate)
558      !!               - recompute depths and water height fields
559      !!
560      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
561      !!               - Recompute:
562      !!                    e3(u/v)_b       
563      !!                    e3w_n           
564      !!                    e3(u/v)w_b     
565      !!                    e3(u/v)w_n     
566      !!                    gdept_n, gdepw_n  and gde3w_n
567      !!                    h(u/v) and h(u/v)r
568      !!
569      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
570      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
571      !!----------------------------------------------------------------------
572      INTEGER, INTENT( in ) ::   kt   ! time step
573      !
574      INTEGER  ::   ji, jj, jk   ! dummy loop indices
575      REAL(wp) ::   zcoef        ! local scalar
576      !!----------------------------------------------------------------------
577      !
578      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
579      !
580      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
581      !
582      IF( kt == nit000 )   THEN
583         IF(lwp) WRITE(numout,*)
584         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
585         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
586      ENDIF
587      !
588      ! Time filter and swap of scale factors
589      ! =====================================
590      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
591      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
592         IF( neuler == 0 .AND. kt == nit000 ) THEN
593            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
594         ELSE
595            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
596            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
597         ENDIF
598         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
599      ENDIF
600      gdept_b(:,:,:) = gdept_n(:,:,:)
601      gdepw_b(:,:,:) = gdepw_n(:,:,:)
602
603      e3t_n(:,:,:) = e3t_a(:,:,:)
604      e3u_n(:,:,:) = e3u_a(:,:,:)
605      e3v_n(:,:,:) = e3v_a(:,:,:)
606
607      ! Compute all missing vertical scale factor and depths
608      ! ====================================================
609      ! Horizontal scale factor interpolations
610      ! --------------------------------------
611      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
612      ! - JC - hu_b, hv_b, hur_b, hvr_b also
613     
614      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
615     
616      ! Vertical scale factor interpolations
617      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
618      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
619      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
620      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
621      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
622      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
623
624      ! t- and w- points depth (set the isf depth as it is in the initial step)
625      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
626      gdepw_n(:,:,1) = 0.0_wp
627      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
628      DO jk = 2, jpk
629         DO jj = 1,jpj
630            DO ji = 1,jpi
631              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
632                                                                 ! 1 for jk = mikt
633               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
634               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
635               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
636                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
637               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
638            END DO
639         END DO
640      END DO
641
642      ! Local depth and Inverse of the local depth of the water
643      ! -------------------------------------------------------
644      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
645      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
646      !
647      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
648      DO jk = 2, jpkm1
649         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
650      END DO
651
652      ! write restart file
653      ! ==================
654      IF( lrst_oce )   CALL dom_vvl_rst( kt, 'WRITE' )
655      !
656      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_sf_swp')
657      !
658   END SUBROUTINE dom_vvl_sf_swp
659
660
661   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
662      !!---------------------------------------------------------------------
663      !!                  ***  ROUTINE dom_vvl__interpol  ***
664      !!
665      !! ** Purpose :   interpolate scale factors from one grid point to another
666      !!
667      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
668      !!                - horizontal interpolation: grid cell surface averaging
669      !!                - vertical interpolation: simple averaging
670      !!----------------------------------------------------------------------
671      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
672      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
673      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
674      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
675      !
676      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
677      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd = T/F
678      !!----------------------------------------------------------------------
679      !
680      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol')
681      !
682      IF(ln_wd) THEN
683        zlnwd = 1.0_wp
684      ELSE
685        zlnwd = 0.0_wp
686      END IF
687      !
688      SELECT CASE ( pout )    !==  type of interpolation  ==!
689         !
690      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
691         DO jk = 1, jpk
692            DO jj = 1, jpjm1
693               DO ji = 1, fs_jpim1   ! vector opt.
694                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
695                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
696                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
697               END DO
698            END DO
699         END DO
700         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
701         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
702         !
703      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
704         DO jk = 1, jpk
705            DO jj = 1, jpjm1
706               DO ji = 1, fs_jpim1   ! vector opt.
707                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
708                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
709                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
710               END DO
711            END DO
712         END DO
713         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
714         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
715         !
716      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
717         DO jk = 1, jpk
718            DO jj = 1, jpjm1
719               DO ji = 1, fs_jpim1   ! vector opt.
720                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
721                     &                       *    r1_e1e2f(ji,jj)                                                  &
722                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
723                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
724               END DO
725            END DO
726         END DO
727         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
728         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
729         !
730      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
731         !
732         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
733         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
734!!gm BUG? use here wmask in case of ISF ?  to be checked
735         DO jk = 2, jpk
736            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
737               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
738               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
739               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
740         END DO
741         !
742      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
743         !
744         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
745         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
746!!gm BUG? use here wumask in case of ISF ?  to be checked
747         DO jk = 2, jpk
748            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
749               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
750               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
751               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
752         END DO
753         !
754      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
755         !
756         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
757         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
758!!gm BUG? use here wvmask in case of ISF ?  to be checked
759         DO jk = 2, jpk
760            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
761               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
762               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
763               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
764         END DO
765      END SELECT
766      !
767      IF( nn_timing == 1 )   CALL timing_stop('dom_vvl_interpol')
768      !
769   END SUBROUTINE dom_vvl_interpol
770
771
772   SUBROUTINE dom_vvl_rst( kt, cdrw )
773      !!---------------------------------------------------------------------
774      !!                   ***  ROUTINE dom_vvl_rst  ***
775      !!                     
776      !! ** Purpose :   Read or write VVL file in restart file
777      !!
778      !! ** Method  :   use of IOM library
779      !!                if the restart does not contain vertical scale factors,
780      !!                they are set to the _0 values
781      !!                if the restart does not contain vertical scale factors increments (z_tilde),
782      !!                they are set to 0.
783      !!----------------------------------------------------------------------
784      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
785      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
786      !
787      INTEGER ::   ji, jj, jk
788      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
789      !!----------------------------------------------------------------------
790      !
791      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
792      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
793         !                                   ! ===============
794         IF( ln_rstart ) THEN                   !* Read the restart file
795            CALL rst_read_open                  !  open the restart file if necessary
796            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
797            !
798            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
799            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
800            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
801            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
802            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
803            !                             ! --------- !
804            !                             ! all cases !
805            !                             ! --------- !
806            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
807               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
808               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
809               ! needed to restart if land processor not computed
810               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
811               WHERE ( tmask(:,:,:) == 0.0_wp ) 
812                  e3t_n(:,:,:) = e3t_0(:,:,:)
813                  e3t_b(:,:,:) = e3t_0(:,:,:)
814               END WHERE
815               IF( neuler == 0 ) THEN
816                  e3t_b(:,:,:) = e3t_n(:,:,:)
817               ENDIF
818            ELSE IF( id1 > 0 ) THEN
819               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
820               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
821               IF(lwp) write(numout,*) 'neuler is forced to 0'
822               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:) )
823               e3t_n(:,:,:) = e3t_b(:,:,:)
824               neuler = 0
825            ELSE IF( id2 > 0 ) THEN
826               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
827               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
828               IF(lwp) write(numout,*) 'neuler is forced to 0'
829               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:) )
830               e3t_b(:,:,:) = e3t_n(:,:,:)
831               neuler = 0
832            ELSE
833               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
834               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
835               IF(lwp) write(numout,*) 'neuler is forced to 0'
836               DO jk = 1, jpk
837                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
838                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
839                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
840               END DO
841               e3t_b(:,:,:) = e3t_n(:,:,:)
842               neuler = 0
843            ENDIF
844            !                             ! ----------- !
845            IF( ln_vvl_zstar ) THEN       ! z_star case !
846               !                          ! ----------- !
847               IF( MIN( id3, id4 ) > 0 ) THEN
848                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
849               ENDIF
850               !                          ! ----------------------- !
851            ELSE                          ! z_tilde and layer cases !
852               !                          ! ----------------------- !
853               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
854                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
855                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
856               ELSE                            ! one at least array is missing
857                  tilde_e3t_b(:,:,:) = 0.0_wp
858                  tilde_e3t_n(:,:,:) = 0.0_wp
859               ENDIF
860               !                          ! ------------ !
861               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
862                  !                       ! ------------ !
863                  IF( id5 > 0 ) THEN  ! required array exists
864                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
865                  ELSE                ! array is missing
866                     hdiv_lf(:,:,:) = 0.0_wp
867                  ENDIF
868               ENDIF
869            ENDIF
870            !
871         ELSE                                   !* Initialize at "rest"
872            !
873            IF( ln_wd .AND. ( cn_cfg == 'wad' ) ) THEN
874              ! Wetting and drying test case
875              CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
876                       tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
877                       sshn (:,:)     = sshb(:,:)
878                       un   (:,:,:)   = ub  (:,:,:)
879                       vn   (:,:,:)   = vb  (:,:,:)
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                     ! potential bug
896                     ! Warning this assumes 2 layers only over wetting locations. needs investigating
897                     e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1
898                     e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1
899                     e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1
900                     sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)           !!gm I don't understand that !
901                     sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
902                     ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
903                  ENDIF
904                ENDDO
905              ENDDO
906               !
907            ELSE
908               !
909               e3t_b(:,:,:) = e3t_0(:,:,:)
910               e3t_n(:,:,:) = e3t_0(:,:,:)
911               sshn(:,:) = 0.0_wp
912               !
913            END IF
914
915            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
916               tilde_e3t_b(:,:,:) = 0._wp
917               tilde_e3t_n(:,:,:) = 0._wp
918               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
919            END IF
920         ENDIF
921         !
922      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
923         !                                   ! ===================
924         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
925         !                                           ! --------- !
926         !                                           ! all cases !
927         !                                           ! --------- !
928         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:) )
929         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:) )
930         !                                           ! ----------------------- !
931         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
932            !                                        ! ----------------------- !
933            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
934            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
935         END IF
936         !                                           ! -------------!   
937         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
938            !                                        ! ------------ !
939            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
940         ENDIF
941         !
942      ENDIF
943      !
944      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
945      !
946   END SUBROUTINE dom_vvl_rst
947
948
949   SUBROUTINE dom_vvl_ctl
950      !!---------------------------------------------------------------------
951      !!                  ***  ROUTINE dom_vvl_ctl  ***
952      !!               
953      !! ** Purpose :   Control the consistency between namelist options
954      !!                for vertical coordinate
955      !!----------------------------------------------------------------------
956      INTEGER ::   ioptio, ios
957      !!
958      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
959         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
960         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
961      !!----------------------------------------------------------------------
962      !
963      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
964      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
965901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
966      !
967      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
968      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
969902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
970      IF(lwm) WRITE ( numond, nam_vvl )
971      !
972      IF(lwp) THEN                    ! Namelist print
973         WRITE(numout,*)
974         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
975         WRITE(numout,*) '~~~~~~~~~~~'
976         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
977         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
978         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
979         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
980         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
981         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
982         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
983         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
984         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
985         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
986         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
987         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
988         IF( ln_vvl_ztilde_as_zstar ) THEN
989            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
990            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
991            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
992            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
993            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
994            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
995         ELSE
996            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
997            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
998            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
999            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
1000         ENDIF
1001         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1002         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1003      ENDIF
1004      !
1005      ioptio = 0                      ! Parameter control
1006      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1007      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1008      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1009      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1010      !
1011      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1012      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1013      !
1014      IF(lwp) THEN                   ! Print the choice
1015         WRITE(numout,*)
1016         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1017         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1018         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1019         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1020         ! - ML - Option not developed yet
1021         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1022         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1023      ENDIF
1024      !
1025#if defined key_agrif
1026      IF(.NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )
1027#endif
1028      !
1029   END SUBROUTINE dom_vvl_ctl
1030
1031   !!======================================================================
1032END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.