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 @ 5866

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

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