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 utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src – NEMO

source: utils/tools_AGRIF_CMEMS_2020/DOMAINcfg/src/domvvl.F90 @ 10727

Last change on this file since 10727 was 10727, checked in by rblod, 5 years ago

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