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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/tests/VORTEX/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/tests/VORTEX/MY_SRC/domvvl.F90 @ 12202

Last change on this file since 12202 was 12202, checked in by cetlod, 4 years ago

dev_merge_option2 : merge in dev_r11613_ENHANCE-04_namelists_as_internalfiles

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