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

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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

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