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/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90 @ 13005

Last change on this file since 13005 was 13005, checked in by gm, 4 years ago

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