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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

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