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

Last change on this file since 9023 was 9023, checked in by timgraham, 6 years ago

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