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

source: branches/2017/dev_r8600_xios_read/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 8612

Last change on this file since 8612 was 8612, checked in by andmirek, 6 years ago

#1953 read single file restart with XIOS

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