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_12905_xios_restart/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/dev_12905_xios_restart/tests/CANAL/MY_SRC/domvvl.F90 @ 12950

Last change on this file since 12950 was 12950, checked in by andmirek, 4 years ago

Ticket #2462: new XIOS restart read/write interfaces
  • Property svn:keywords set to Id
File size: 54.9 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( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
169      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
170      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
171      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
172      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
173      !                                ! Vertical interpolation of e3t,u,v
174      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
175      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
176      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
177      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
178      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
179      CALL dom_vvl_interpol( e3v(:,:,:,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   END SUBROUTINE dom_vvl_zgr
273
274
275   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
276      !!----------------------------------------------------------------------
277      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
278      !!                   
279      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
280      !!                 tranxt and dynspg routines
281      !!
282      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
283      !!               - z_tilde_case: after scale factor increment =
284      !!                                    high frequency part of horizontal divergence
285      !!                                  + retsoring towards the background grid
286      !!                                  + thickness difusion
287      !!                               Then repartition of ssh INCREMENT proportionnaly
288      !!                               to the "baroclinic" level thickness.
289      !!
290      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
291      !!               - tilde_e3t_a: after increment of vertical scale factor
292      !!                              in z_tilde case
293      !!               - e3(t/u/v)_a
294      !!
295      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
296      !!----------------------------------------------------------------------
297      INTEGER, INTENT( in )           ::   kt             ! time step
298      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
299      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
300      !
301      INTEGER                ::   ji, jj, jk            ! dummy loop indices
302      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
303      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
304      LOGICAL                ::   ll_do_bclinic         ! local logical
305      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
306      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
307      !!----------------------------------------------------------------------
308      !
309      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
310      !
311      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
312      !
313      IF( kt == nit000 ) THEN
314         IF(lwp) WRITE(numout,*)
315         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
316         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
317      ENDIF
318
319      ll_do_bclinic = .TRUE.
320      IF( PRESENT(kcall) ) THEN
321         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
322      ENDIF
323
324      ! ******************************* !
325      ! After acale factors at t-points !
326      ! ******************************* !
327      !                                                ! --------------------------------------------- !
328      !                                                ! z_star coordinate and barotropic z-tilde part !
329      !                                                ! --------------------------------------------- !
330      !
331      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
332      DO jk = 1, jpkm1
333         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
334         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
335      END DO
336      !
337      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
338         !                                                               ! ------baroclinic part------ !
339         ! I - initialization
340         ! ==================
341
342         ! 1 - barotropic divergence
343         ! -------------------------
344         zhdiv(:,:) = 0._wp
345         zht(:,:)   = 0._wp
346         DO jk = 1, jpkm1
347            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
348            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
349         END DO
350         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
351
352         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
353         ! --------------------------------------------------
354         IF( ln_vvl_ztilde ) THEN
355            IF( kt > nit000 ) THEN
356               DO jk = 1, jpkm1
357                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
358                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
359               END DO
360            ENDIF
361         ENDIF
362
363         ! II - after z_tilde increments of vertical scale factors
364         ! =======================================================
365         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
366
367         ! 1 - High frequency divergence term
368         ! ----------------------------------
369         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
370            DO jk = 1, jpkm1
371               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
372            END DO
373         ELSE                         ! layer case
374            DO jk = 1, jpkm1
375               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
376            END DO
377         ENDIF
378
379         ! 2 - Restoring term (z-tilde case only)
380         ! ------------------
381         IF( ln_vvl_ztilde ) THEN
382            DO jk = 1, jpk
383               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
384            END DO
385         ENDIF
386
387         ! 3 - Thickness diffusion term
388         ! ----------------------------
389         zwu(:,:) = 0._wp
390         zwv(:,:) = 0._wp
391         DO_3D_10_10( 1, jpkm1 )
392            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
393               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
394            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
395               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
396            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
397            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
398         END_3D
399         DO_2D_11_11
400            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
401            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
402         END_2D
403         DO_3D_00_00( 1, jpkm1 )
404            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
405               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
406               &                                            ) * r1_e1e2t(ji,jj)
407         END_3D
408         !                       ! d - thickness diffusion transport: boundary conditions
409         !                             (stored for tracer advction and continuity equation)
410         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
411
412         ! 4 - Time stepping of baroclinic scale factors
413         ! ---------------------------------------------
414         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
415         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
416
417         ! Maximum deformation control
418         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
419         ze3t(:,:,jpk) = 0._wp
420         DO jk = 1, jpkm1
421            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
422         END DO
423         z_tmax = MAXVAL( ze3t(:,:,:) )
424         CALL mpp_max( 'domvvl', z_tmax )                 ! max over the global domain
425         z_tmin = MINVAL( ze3t(:,:,:) )
426         CALL mpp_min( 'domvvl', z_tmin )                 ! min over the global domain
427         ! - ML - test: for the moment, stop simulation for too large e3_t variations
428         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
429            IF( lk_mpp ) THEN
430               CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max )
431               CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min )
432            ELSE
433               ijk_max = MAXLOC( ze3t(:,:,:) )
434               ijk_max(1) = ijk_max(1) + nimpp - 1
435               ijk_max(2) = ijk_max(2) + njmpp - 1
436               ijk_min = MINLOC( ze3t(:,:,:) )
437               ijk_min(1) = ijk_min(1) + nimpp - 1
438               ijk_min(2) = ijk_min(2) + njmpp - 1
439            ENDIF
440            IF (lwp) THEN
441               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
442               WRITE(numout, *) 'at i, j, k=', ijk_max
443               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
444               WRITE(numout, *) 'at i, j, k=', ijk_min           
445               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
446            ENDIF
447         ENDIF
448         ! - ML - end test
449         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
450         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
451         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
452
453         !
454         ! "tilda" change in the after scale factor
455         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456         DO jk = 1, jpkm1
457            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
458         END DO
459         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
460         ! ==================================================================================
461         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
462         ! - ML - baroclinicity error should be better treated in the future
463         !        i.e. locally and not spread over the water column.
464         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
465         zht(:,:) = 0.
466         DO jk = 1, jpkm1
467            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
468         END DO
469         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
470         DO jk = 1, jpkm1
471            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
472         END DO
473
474      ENDIF
475
476      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
477      !                                           ! ---baroclinic part--------- !
478         DO jk = 1, jpkm1
479            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
480         END DO
481      ENDIF
482
483      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
484         !
485         IF( lwp ) WRITE(numout, *) 'kt =', kt
486         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
487            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
488            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
489            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
490         END IF
491         !
492         zht(:,:) = 0.0_wp
493         DO jk = 1, jpkm1
494            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
495         END DO
496         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
497         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
498         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
499         !
500         zht(:,:) = 0.0_wp
501         DO jk = 1, jpkm1
502            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
503         END DO
504         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
505         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
506         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
507         !
508         zht(:,:) = 0.0_wp
509         DO jk = 1, jpkm1
510            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
511         END DO
512         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
513         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
514         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
515         !
516         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
517         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
518         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
519         !
520         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
521         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
522         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
523         !
524         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
525         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
526         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
527      END IF
528
529      ! *********************************** !
530      ! After scale factors at u- v- points !
531      ! *********************************** !
532
533      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
534      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
535
536      ! *********************************** !
537      ! After depths at u- v points         !
538      ! *********************************** !
539
540      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
541      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
542      DO jk = 2, jpkm1
543         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
544         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
545      END DO
546      !                                        ! Inverse of the local depth
547!!gm BUG ?  don't understand the use of umask_i here .....
548      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
549      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
550      !
551      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
552      !
553   END SUBROUTINE dom_vvl_sf_nxt
554
555
556   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
557      !!----------------------------------------------------------------------
558      !!                ***  ROUTINE dom_vvl_sf_update  ***
559      !!                   
560      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
561      !!               compute all depths and related variables for next time step
562      !!               write outputs and restart file
563      !!
564      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
565      !!               - reconstruct scale factor at other grid points (interpolate)
566      !!               - recompute depths and water height fields
567      !!
568      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
569      !!               - Recompute:
570      !!                    e3(u/v)_b       
571      !!                    e3w(:,:,:,Kmm)           
572      !!                    e3(u/v)w_b     
573      !!                    e3(u/v)w_n     
574      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
575      !!                    h(u/v) and h(u/v)r
576      !!
577      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
578      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
579      !!----------------------------------------------------------------------
580      INTEGER, INTENT( in ) ::   kt              ! time step
581      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
582      !
583      INTEGER  ::   ji, jj, jk   ! dummy loop indices
584      REAL(wp) ::   zcoef        ! local scalar
585      !!----------------------------------------------------------------------
586      !
587      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
588      !
589      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
590      !
591      IF( kt == nit000 )   THEN
592         IF(lwp) WRITE(numout,*)
593         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
594         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
595      ENDIF
596      !
597      ! Time filter and swap of scale factors
598      ! =====================================
599      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
600      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
601         IF( l_1st_euler ) THEN
602            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
603         ELSE
604            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
605            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
606         ENDIF
607         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
608      ENDIF
609
610      ! Compute all missing vertical scale factor and depths
611      ! ====================================================
612      ! Horizontal scale factor interpolations
613      ! --------------------------------------
614      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
615      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
616     
617      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
618     
619      ! Vertical scale factor interpolations
620      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
621      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
622      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
623      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
624      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
625      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
626
627      ! t- and w- points depth (set the isf depth as it is in the initial step)
628      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
629      gdepw(:,:,1,Kmm) = 0.0_wp
630      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
631      DO_3D_11_11( 2, jpk )
632        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
633                                                           ! 1 for jk = mikt
634         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
635         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
636         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
637             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
638         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
639      END_3D
640
641      ! Local depth and Inverse of the local depth of the water
642      ! -------------------------------------------------------
643      !
644      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
645      DO jk = 2, jpkm1
646         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
647      END DO
648
649      ! write restart file
650      ! ==================
651      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
652      !
653      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
654      !
655   END SUBROUTINE dom_vvl_sf_update
656
657
658   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
659      !!---------------------------------------------------------------------
660      !!                  ***  ROUTINE dom_vvl__interpol  ***
661      !!
662      !! ** Purpose :   interpolate scale factors from one grid point to another
663      !!
664      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
665      !!                - horizontal interpolation: grid cell surface averaging
666      !!                - vertical interpolation: simple averaging
667      !!----------------------------------------------------------------------
668      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
669      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
670      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
671      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
672      !
673      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
674      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
675      !!----------------------------------------------------------------------
676      !
677      IF(ln_wd_il) THEN
678        zlnwd = 1.0_wp
679      ELSE
680        zlnwd = 0.0_wp
681      END IF
682      !
683      SELECT CASE ( pout )    !==  type of interpolation  ==!
684         !
685      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
686         DO_3D_10_10( 1, jpk )
687            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
688               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
689               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
690         END_3D
691         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
692         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
693         !
694      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
695         DO_3D_10_10( 1, jpk )
696            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
697               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
698               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
699         END_3D
700         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
701         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
702         !
703      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
704         DO_3D_10_10( 1, jpk )
705            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
706               &                       *    r1_e1e2f(ji,jj)                                                  &
707               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
708               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
709         END_3D
710         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
711         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
712         !
713      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
714         !
715         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
716         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
717!!gm BUG? use here wmask in case of ISF ?  to be checked
718         DO jk = 2, jpk
719            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
720               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
721               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
722               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
723         END DO
724         !
725      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
726         !
727         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
728         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
729!!gm BUG? use here wumask in case of ISF ?  to be checked
730         DO jk = 2, jpk
731            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
732               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
733               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
734               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
735         END DO
736         !
737      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
738         !
739         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
740         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
741!!gm BUG? use here wvmask in case of ISF ?  to be checked
742         DO jk = 2, jpk
743            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
744               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
745               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
746               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
747         END DO
748      END SELECT
749      !
750   END SUBROUTINE dom_vvl_interpol
751
752
753   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
754      !!---------------------------------------------------------------------
755      !!                   ***  ROUTINE dom_vvl_rst  ***
756      !!                     
757      !! ** Purpose :   Read or write VVL file in restart file
758      !!
759      !! ** Method  :   use of IOM library
760      !!                if the restart does not contain vertical scale factors,
761      !!                they are set to the _0 values
762      !!                if the restart does not contain vertical scale factors increments (z_tilde),
763      !!                they are set to 0.
764      !!----------------------------------------------------------------------
765      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
766      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
767      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
768      !
769      INTEGER ::   ji, jj, jk
770      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
771      !!----------------------------------------------------------------------
772      !
773      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
774         !                                   ! ===============
775         IF( ln_rstart ) THEN                   !* Read the restart file
776            CALL rst_read_open                  !  open the restart file if necessary
777            IF( lrxios ) CALL iom_swap( TRIM(crxios_context) )
778            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
779            !
780            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
781            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
782            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
783            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
784            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
785            !
786            !                             ! --------- !
787            !                             ! all cases !
788            !                             ! --------- !
789            !
790            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
791               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
792               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
793               ! needed to restart if land processor not computed
794               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
795               WHERE ( tmask(:,:,:) == 0.0_wp ) 
796                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
797                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
798               END WHERE
799               IF( l_1st_euler ) THEN
800                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
801               ENDIF
802            ELSE IF( id1 > 0 ) THEN
803               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
804               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
805               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
806               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
807               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
808               l_1st_euler = .true.
809            ELSE IF( id2 > 0 ) THEN
810               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
811               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
812               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
813               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
814               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
815               l_1st_euler = .true.
816            ELSE
817               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
818               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
819               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
820               DO jk = 1, jpk
821                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
822                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
823                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
824               END DO
825               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
826               l_1st_euler = .true.
827            ENDIF
828            !                             ! ----------- !
829            IF( ln_vvl_zstar ) THEN       ! z_star case !
830               !                          ! ----------- !
831               IF( MIN( id3, id4 ) > 0 ) THEN
832                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
833               ENDIF
834               !                          ! ----------------------- !
835            ELSE                          ! z_tilde and layer cases !
836               !                          ! ----------------------- !
837               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
838                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
839                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
840               ELSE                            ! one at least array is missing
841                  tilde_e3t_b(:,:,:) = 0.0_wp
842                  tilde_e3t_n(:,:,:) = 0.0_wp
843               ENDIF
844               !                          ! ------------ !
845               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
846                  !                       ! ------------ !
847                  IF( id5 > 0 ) THEN  ! required array exists
848                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
849                  ELSE                ! array is missing
850                     hdiv_lf(:,:,:) = 0.0_wp
851                  ENDIF
852               ENDIF
853            ENDIF
854            IF( lrxios ) CALL iom_swap( TRIM(cxios_context) )
855            !
856         ELSE                                   !* Initialize at "rest"
857            !
858
859            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
860               !
861               IF( cn_cfg == 'wad' ) THEN
862                  ! Wetting and drying test case
863                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
864                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
865                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
866                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
867                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
868               ELSE
869                  ! if not test case
870                  ssh(:,:,Kmm) = -ssh_ref
871                  ssh(:,:,Kbb) = -ssh_ref
872
873                  DO_2D_11_11
874                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
875                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
876                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
877                     ENDIF
878                  END_2D
879               ENDIF !If test case else
880
881               ! Adjust vertical metrics for all wad
882               DO jk = 1, jpk
883                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
884                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
885                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
886               END DO
887               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
888
889               DO_2D_11_11
890                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
891                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
892                  ENDIF
893               END_2D
894               !
895            ELSE
896               !
897               ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm)
898               !
899               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
900               !
901               ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v
902               !
903               DO jk=1,jpk
904                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
905                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
906                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points
907               END DO
908               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
909               ssh(:,:,Kmm) = ssh(:,:,Kbb)                                     ! needed later for gde3w
910               !
911            END IF           ! end of ll_wd edits
912
913            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
914               tilde_e3t_b(:,:,:) = 0._wp
915               tilde_e3t_n(:,:,:) = 0._wp
916               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
917            END IF
918         ENDIF
919         !
920      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
921         !                                   ! ===================
922         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
923         IF( lwxios ) CALL iom_swap(      cwxios_context          )
924         !                                           ! --------- !
925         !                                           ! all cases !
926         !                                           ! --------- !
927         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
928         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
929         !                                           ! ----------------------- !
930         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
931            !                                        ! ----------------------- !
932            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
933            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
934         END IF
935         !                                           ! -------------!   
936         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
937            !                                        ! ------------ !
938            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
939         ENDIF
940         !
941         IF( lwxios ) CALL iom_swap(      cxios_context          )
942      ENDIF
943      !
944   END SUBROUTINE dom_vvl_rst
945
946
947   SUBROUTINE dom_vvl_ctl
948      !!---------------------------------------------------------------------
949      !!                  ***  ROUTINE dom_vvl_ctl  ***
950      !!               
951      !! ** Purpose :   Control the consistency between namelist options
952      !!                for vertical coordinate
953      !!----------------------------------------------------------------------
954      INTEGER ::   ioptio, ios
955      !!
956      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
957         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
958         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
959      !!----------------------------------------------------------------------
960      !
961      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
962901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
963      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
964902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
965      IF(lwm) WRITE ( numond, nam_vvl )
966      !
967      IF(lwp) THEN                    ! Namelist print
968         WRITE(numout,*)
969         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
970         WRITE(numout,*) '~~~~~~~~~~~'
971         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
972         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
973         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
974         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
975         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
976         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
977         WRITE(numout,*) '      !'
978         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
979         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
980         IF( ln_vvl_ztilde_as_zstar ) THEN
981            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
982            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
983            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
984            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
985            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
986            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
987         ELSE
988            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
989            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
990         ENDIF
991         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
992      ENDIF
993      !
994      ioptio = 0                      ! Parameter control
995      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
996      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
997      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
998      IF( ln_vvl_layer           )   ioptio = ioptio + 1
999      !
1000      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1001      !
1002      IF(lwp) THEN                   ! Print the choice
1003         WRITE(numout,*)
1004         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1005         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1006         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1007         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1008      ENDIF
1009      !
1010#if defined key_agrif
1011      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1012#endif
1013      !
1014   END SUBROUTINE dom_vvl_ctl
1015
1016   !!======================================================================
1017END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.