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

Last change on this file since 9119 was 9090, checked in by flavoni, 6 years ago

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