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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9441

Last change on this file since 9441 was 9404, checked in by mathiot, 6 years ago

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