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/2017/dev_r8600_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2017/dev_r8600_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 8651

Last change on this file since 8651 was 8651, checked in by andmirek, 7 years ago

change names of the new subroutines

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