New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in branches/2017/dev_r8600_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2017/dev_r8600_xios_write/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 8644

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

ticket #1962 xios write functionality works

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