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 NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90 @ 9863

Last change on this file since 9863 was 9863, checked in by gm, 7 years ago

#1911 (ENHANCE-04): simplified implementation of the Euler stepping at nit000

  • 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)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors
59   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_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/OCE 4.0 , NEMO Consortium (2018)
67   !! $Id$
68   !! Software governed by the CeCILL licence     (./LICENSE)
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( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   &
79            &      dte3t_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 te3t_(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), te3t_(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      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
151      e3t_a(:,:,:) = e3t_n(:,:,:)
152      e3u_a(:,:,:) = e3u_n(:,:,:)
153      e3v_a(:,:,:) = e3v_n(:,:,:)
154      !
155      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
156      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
157      gdepw_n(:,:,1) = 0.0_wp
158      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
159      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
160      gdepw_b(:,:,1) = 0.0_wp
161      DO jk = 2, jpk                               ! vertical sum
162         DO jj = 1,jpj
163            DO ji = 1,jpi
164               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
165               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
166               !                             ! 0.5 where jk = mikt     
167!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
168               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
169               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
170               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
171                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
172               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
173               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
174               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
175                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
176            END DO
177         END DO
178      END DO
179      !
180      !                    !==  thickness of the water column  !!   (ocean portion only)
181      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
182      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
183      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
184      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
185      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
186      DO jk = 2, jpkm1
187         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
188         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
189         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
190         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
191         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
192      END DO
193      !
194      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
195      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
196      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
197      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
198      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
199
200      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
201      IF( ln_vvl_ztilde ) THEN
202!!gm : idea: add here a READ in a file of custumized restoring frequency
203         !                                   ! Values in days provided via the namelist
204         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
205         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
206         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
207         !
208         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
209            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
210            frq_rst_hdv(:,:) = 1._wp / rdt
211         ENDIF
212         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
213            DO jj = 1, jpj
214               DO ji = 1, jpi
215!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
216                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
217                     ! values outside the equatorial band and transition zone (ztilde)
218                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
219                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
220                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
221                     ! values inside the equatorial band (ztilde as zstar)
222                     frq_rst_e3t(ji,jj) =  0.0_wp
223                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
224                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
225                     !                                      ! (linearly transition from z-tilde to z-star)
226                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
227                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
228                        &                                          * 180._wp / 3.5_wp ) )
229                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
230                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
231                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
232                        &                                          * 180._wp / 3.5_wp ) )
233                  ENDIF
234               END DO
235            END DO
236            IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
237               ii0 = 103   ;   ii1 = 111       
238               ij0 = 128   ;   ij1 = 135   ;   
239               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
240               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
241            ENDIF
242         ENDIF
243      ENDIF
244      !
245      IF(lwxios) THEN
246! define variables in restart file when writing with XIOS
247         CALL iom_set_rstw_var_active('e3t_b')
248         CALL iom_set_rstw_var_active('e3t_n')
249         !                                           ! ----------------------- !
250         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
251            !                                        ! ----------------------- !
252            CALL iom_set_rstw_var_active('tilde_e3t_b')
253            CALL iom_set_rstw_var_active('tilde_e3t_n')
254         END IF
255         !                                           ! -------------!   
256         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
257            !                                        ! ------------ !
258            CALL iom_set_rstw_var_active('hdiv_lf')
259         ENDIF
260         !
261      ENDIF
262      !
263   END SUBROUTINE dom_vvl_init
264
265
266   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
267      !!----------------------------------------------------------------------
268      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
269      !!                   
270      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
271      !!                 tranxt and dynspg routines
272      !!
273      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
274      !!               - z_tilde_case: after scale factor increment =
275      !!                                    high frequency part of horizontal divergence
276      !!                                  + retsoring towards the background grid
277      !!                                  + thickness difusion
278      !!                               Then repartition of ssh INCREMENT proportionnaly
279      !!                               to the "baroclinic" level thickness.
280      !!
281      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
282      !!               - te3t_a: after increment of vertical scale factor
283      !!                              in z_tilde case
284      !!               - e3(t/u/v)_a
285      !!
286      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
287      !!----------------------------------------------------------------------
288      INTEGER, INTENT( in )           ::   kt      ! time step
289      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
290      !
291      INTEGER                ::   ji, jj, jk            ! dummy loop indices
292      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
293      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
294      LOGICAL                ::   ll_do_bclinic         ! local logical
295      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
296      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
297      !!----------------------------------------------------------------------
298      !
299      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
300      !
301      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
302      !
303      IF( kt == nit000 ) THEN
304         IF(lwp) WRITE(numout,*)
305         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
306         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
307      ENDIF
308
309      ll_do_bclinic = .TRUE.
310      IF( PRESENT(kcall) ) THEN
311         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
312      ENDIF
313
314      ! ******************************* !
315      ! After acale factors at t-points !
316      ! ******************************* !
317      !                                                ! --------------------------------------------- !
318      !                                                ! z_star coordinate and barotropic z-tilde part !
319      !                                                ! --------------------------------------------- !
320      !
321      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
322      DO jk = 1, jpkm1
323         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
324         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
325      END DO
326      !
327      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
328         !                                                            ! ------baroclinic part------ !
329         ! I - initialization
330         ! ==================
331
332         ! 1 - barotropic divergence
333         ! -------------------------
334         zhdiv(:,:) = 0._wp
335         zht(:,:)   = 0._wp
336         DO jk = 1, jpkm1
337            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
338            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
339         END DO
340         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
341
342         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
343         ! --------------------------------------------------
344         IF( ln_vvl_ztilde ) THEN
345            IF( kt > nit000 ) THEN
346               DO jk = 1, jpkm1
347                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
348                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
349               END DO
350            ENDIF
351         ENDIF
352
353         ! II - after z_tilde increments of vertical scale factors
354         ! =======================================================
355         te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms
356
357         ! 1 - High frequency divergence term
358         ! ----------------------------------
359         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
360            DO jk = 1, jpkm1
361               te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
362            END DO
363         ELSE                         ! layer case
364            DO jk = 1, jpkm1
365               te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
366            END DO
367         ENDIF
368
369         ! 2 - Restoring term (z-tilde case only)
370         ! ------------------
371         IF( ln_vvl_ztilde ) THEN
372            DO jk = 1, jpk
373               te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk)
374            END DO
375         ENDIF
376
377         ! 3 - Thickness diffusion term
378         ! ----------------------------
379         zwu(:,:) = 0._wp
380         zwv(:,:) = 0._wp
381         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
382            DO jj = 1, jpjm1
383               DO ji = 1, fs_jpim1   ! vector opt.
384                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
385                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) )
386                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
387                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji  ,jj+1,jk) )
388                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
389                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
390               END DO
391            END DO
392         END DO
393         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
394            DO ji = 1, jpi
395               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
396               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
397            END DO
398         END DO
399         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
400            DO jj = 2, jpjm1
401               DO ji = fs_2, fs_jpim1   ! vector opt.
402                  te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
403                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
404                     &                                            ) * r1_e1e2t(ji,jj)
405               END DO
406            END DO
407         END DO
408         !                       ! d - thickness diffusion transport: boundary conditions
409         !                             (stored for tracer advction and continuity equation)
410         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
411
412         ! 4 - Time stepping of baroclinic scale factors
413         ! ---------------------------------------------
414         ! Leapfrog time stepping
415         ! ~~~~~~~~~~~~~~~~~~~~~~
416         CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp )
417         te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:)
418
419         ! Maximum deformation control
420         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
421         ze3t(:,:,jpk) = 0._wp
422         DO jk = 1, jpkm1
423            ze3t(:,:,jk) = te3t_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( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
444               WRITE(numout, *) 'at i, j, k=', ijk_max
445               WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
446               WRITE(numout, *) 'at i, j, k=', ijk_min           
447               CALL ctl_warn('MAX( ABS( te3t_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         te3t_a(:,:,:) = MIN( te3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
453         te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
454
455         !
456         ! "tilda" change in the after scale factor
457         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
458         DO jk = 1, jpkm1
459            dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_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._wp
468         DO jk = 1, jpkm1
469            zht(:,:)  = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk)
470         END DO
471         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
472         DO jk = 1, jpkm1
473            dte3t_a(:,:,jk) = dte3t_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) + dte3t_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(te3t_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), te3t_(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, ze3f   ! 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( l_1st_euler ) THEN
603            te3t_n(:,:,:) = te3t_a(:,:,:)
604         ELSE
605            DO jk = 1, jpk
606               DO jj = 1, jpj
607                  DO ji = 1, jpi
608                     ze3f = te3t_n(ji,jj,jk)   &
609                        & + atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) )
610                     te3t_b(ji,jj,jk) = ze3f
611                     te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk)
612                  END DO
613               END DO
614            END DO
615         ENDIF
616      ENDIF
617      gdept_b(:,:,:) = gdept_n(:,:,:)
618      gdepw_b(:,:,:) = gdepw_n(:,:,:)
619
620      e3t_n(:,:,:) = e3t_a(:,:,:)
621      e3u_n(:,:,:) = e3u_a(:,:,:)
622      e3v_n(:,:,:) = e3v_a(:,:,:)
623
624      ! Compute all missing vertical scale factor and depths
625      ! ====================================================
626      ! Horizontal scale factor interpolations
627      ! --------------------------------------
628      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
629      ! - JC - hu_b, hv_b, hur_b, hvr_b also
630     
631      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
632     
633      ! Vertical scale factor interpolations
634      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
635      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
636      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
637      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
638      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
639      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
640
641      ! t- and w- points depth (set the isf depth as it is in the initial step)
642      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
643      gdepw_n(:,:,1) = 0.0_wp
644      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
645      DO jk = 2, jpk
646         DO jj = 1,jpj
647            DO ji = 1,jpi
648              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
649                                                                 ! 1 for jk = mikt
650               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
651               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
652               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
653                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
654               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
655            END DO
656         END DO
657      END DO
658
659      ! Local depth and Inverse of the local depth of the water
660      ! -------------------------------------------------------
661      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
662      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
663      !
664      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
665      DO jk = 2, jpkm1
666         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
667      END DO
668
669      ! write restart file
670      ! ==================
671      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
672      !
673      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
674      !
675   END SUBROUTINE dom_vvl_sf_swp
676
677
678   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
679      !!---------------------------------------------------------------------
680      !!                  ***  ROUTINE dom_vvl__interpol  ***
681      !!
682      !! ** Purpose :   interpolate scale factors from one grid point to another
683      !!
684      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
685      !!                - horizontal interpolation: grid cell surface averaging
686      !!                - vertical interpolation: simple averaging
687      !!----------------------------------------------------------------------
688      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
689      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
690      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
691      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
692      !
693      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
694      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
695      !!----------------------------------------------------------------------
696      !
697      IF(ln_wd_il) THEN
698        zlnwd = 1.0_wp
699      ELSE
700        zlnwd = 0.0_wp
701      END IF
702      !
703      SELECT CASE ( pout )    !==  type of interpolation  ==!
704         !
705      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
706         DO jk = 1, jpk
707            DO jj = 1, jpjm1
708               DO ji = 1, fs_jpim1   ! vector opt.
709                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
710                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
711                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
712               END DO
713            END DO
714         END DO
715         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
716         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
717         !
718      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
719         DO jk = 1, jpk
720            DO jj = 1, jpjm1
721               DO ji = 1, fs_jpim1   ! vector opt.
722                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
723                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
724                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
725               END DO
726            END DO
727         END DO
728         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
729         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
730         !
731      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
732         DO jk = 1, jpk
733            DO jj = 1, jpjm1
734               DO ji = 1, fs_jpim1   ! vector opt.
735                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
736                     &                       *    r1_e1e2f(ji,jj)                                                  &
737                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
738                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
739               END DO
740            END DO
741         END DO
742         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
743         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
744         !
745      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
746         !
747         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
748         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
749!!gm BUG? use here wmask in case of ISF ?  to be checked
750         DO jk = 2, jpk
751            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
752               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
753               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
754               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
755         END DO
756         !
757      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
758         !
759         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
760         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
761!!gm BUG? use here wumask in case of ISF ?  to be checked
762         DO jk = 2, jpk
763            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
764               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
765               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
766               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
767         END DO
768         !
769      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
770         !
771         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
772         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
773!!gm BUG? use here wvmask in case of ISF ?  to be checked
774         DO jk = 2, jpk
775            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
776               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
777               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
778               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
779         END DO
780      END SELECT
781      !
782   END SUBROUTINE dom_vvl_interpol
783
784
785   SUBROUTINE dom_vvl_rst( kt, cdrw )
786      !!---------------------------------------------------------------------
787      !!                   ***  ROUTINE dom_vvl_rst  ***
788      !!                     
789      !! ** Purpose :   Read or write VVL file in restart file
790      !!
791      !! ** Method  :   use of IOM library
792      !!                if the restart does not contain vertical scale factors,
793      !!                they are set to the _0 values
794      !!                if the restart does not contain vertical scale factors increments (z_tilde),
795      !!                they are set to 0.
796      !!----------------------------------------------------------------------
797      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
798      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
799      !
800      INTEGER ::   ji, jj, jk
801      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
802      !!----------------------------------------------------------------------
803      !
804      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
805         !                                   ! ===============
806         IF( ln_rstart ) THEN                   !* Read the restart file
807            CALL rst_read_open                  !  open the restart file if necessary
808            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
809            !
810            id1 = iom_varid( numror, 'e3t_b'      , ldstop = .FALSE. )
811            id2 = iom_varid( numror, 'e3t_n'      , ldstop = .FALSE. )
812            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
813            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
814            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. )
815            !                             ! --------- !
816            !                             ! all cases !
817            !                             ! --------- !
818            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
819               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
820               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
821               ! needed to restart if land processor not computed
822               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
823               WHERE ( tmask(:,:,:) == 0.0_wp ) 
824                  e3t_n(:,:,:) = e3t_0(:,:,:)
825                  e3t_b(:,:,:) = e3t_0(:,:,:)
826               END WHERE
827               IF( l_1st_euler ) THEN
828                  e3t_b(:,:,:) = e3t_n(:,:,:)
829               ENDIF
830            ELSE IF( id1 > 0 ) THEN
831               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
832               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
833               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
834               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
835               e3t_n(:,:,:) = e3t_b(:,:,:)
836               l_1st_euler = .TRUE.
837            ELSE IF( id2 > 0 ) THEN
838               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
839               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
840               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
841               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
842               e3t_b(:,:,:) = e3t_n(:,:,:)
843               l_1st_euler = .TRUE.
844            ELSE
845               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
846               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
847               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
848               DO jk = 1, jpk
849                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
850                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
851                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
852               END DO
853               e3t_b(:,:,:) = e3t_n(:,:,:)
854               l_1st_euler = .TRUE.
855            ENDIF
856            !                             ! ----------- !
857            IF( ln_vvl_zstar ) THEN       ! z_star case !
858               !                          ! ----------- !
859               IF( MIN( id3, id4 ) > 0 ) THEN
860                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
861               ENDIF
862               !                          ! ----------------------- !
863            ELSE                          ! z_tilde and layer cases !
864               !                          ! ----------------------- !
865               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
866                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios )
867                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios )
868               ELSE                            ! one at least array is missing
869                  te3t_b(:,:,:) = 0.0_wp
870                  te3t_n(:,:,:) = 0.0_wp
871               ENDIF
872               !                          ! ------------ !
873               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
874                  !                       ! ------------ !
875                  IF( id5 > 0 ) THEN  ! required array exists
876                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
877                  ELSE                ! array is missing
878                     hdiv_lf(:,:,:) = 0.0_wp
879                  ENDIF
880               ENDIF
881            ENDIF
882            !
883         ELSE                                   !* Initialize at "rest"
884            !
885
886            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
887               !
888               IF( cn_cfg == 'wad' ) THEN
889                  ! Wetting and drying test case
890                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
891                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
892                  sshn (:,:)     = sshb(:,:)
893                  un   (:,:,:)   = ub  (:,:,:)
894                  vn   (:,:,:)   = vb  (:,:,:)
895               ELSE
896                  ! if not test case
897                  sshn(:,:) = -ssh_ref
898                  sshb(:,:) = -ssh_ref
899
900                  DO jj = 1, jpj
901                     DO ji = 1, jpi
902                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
903
904                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
905                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
906                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
907                        ENDIF
908                     ENDDO
909                  ENDDO
910               ENDIF !If test case else
911
912               ! Adjust vertical metrics for all wad
913               DO jk = 1, jpk
914                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
915                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
916                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
917               END DO
918               e3t_b(:,:,:) = e3t_n(:,:,:)
919
920               DO ji = 1, jpi
921                  DO jj = 1, jpj
922                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
923                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
924                     ENDIF
925                  END DO
926               END DO 
927               !
928            ELSE
929               !
930               ! Just to read set ssh in fact, called latter once vertical grid
931               ! is set up:
932!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
933!               !
934!               DO jk=1,jpk
935!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
936!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
937!               END DO
938!               e3t_n(:,:,:) = e3t_b(:,:,:)
939                sshn(:,:)=0._wp
940                e3t_n(:,:,:)=e3t_0(:,:,:)
941                e3t_b(:,:,:)=e3t_0(:,:,:)
942               !
943            END IF           ! end of ll_wd edits
944
945            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
946               te3t_b(:,:,:) = 0._wp
947               te3t_n(:,:,:) = 0._wp
948               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
949            END IF
950         ENDIF
951         !
952      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
953         !                                   ! ===================
954         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
955         IF( lwxios ) CALL iom_swap(      cwxios_context          )
956         !                                           ! --------- !
957         !                                           ! all cases !
958         !                                           ! --------- !
959         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
960         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
961         !                                           ! ----------------------- !
962         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
963            !                                        ! ----------------------- !
964            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios)
965            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios)
966         END IF
967         !                                           ! -------------!   
968         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
969            !                                        ! ------------ !
970            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
971         ENDIF
972         !
973         IF( lwxios ) CALL iom_swap(      cxios_context          )
974      ENDIF
975      !
976   END SUBROUTINE dom_vvl_rst
977
978
979   SUBROUTINE dom_vvl_ctl
980      !!---------------------------------------------------------------------
981      !!                  ***  ROUTINE dom_vvl_ctl  ***
982      !!               
983      !! ** Purpose :   Control the consistency between namelist options
984      !!                for vertical coordinate
985      !!----------------------------------------------------------------------
986      INTEGER ::   ioptio, ios
987      !!
988      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
989         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
990         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
991      !!----------------------------------------------------------------------
992      !
993      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
994      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
995901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
996      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
997      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
998902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
999      IF(lwm) WRITE ( numond, nam_vvl )
1000      !
1001      IF(lwp) THEN                    ! Namelist print
1002         WRITE(numout,*)
1003         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1004         WRITE(numout,*) '~~~~~~~~~~~'
1005         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1006         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1007         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1008         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1009         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1010         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1011         WRITE(numout,*) '      !'
1012         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1013         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1014         IF( ln_vvl_ztilde_as_zstar ) THEN
1015            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1016            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1017            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1018            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1019            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1020            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1021         ELSE
1022            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1023            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1024         ENDIF
1025         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1026      ENDIF
1027      !
1028      ioptio = 0                      ! Parameter control
1029      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1030      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1031      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1032      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1033      !
1034      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1035      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1036      !
1037      IF(lwp) THEN                   ! Print the choice
1038         WRITE(numout,*)
1039         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1040         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1041         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1042         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1043      ENDIF
1044      !
1045#if defined key_agrif
1046      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1047#endif
1048      !
1049   END SUBROUTINE dom_vvl_ctl
1050
1051   !!======================================================================
1052END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.