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.
domqe.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90 @ 12584

Last change on this file since 12584 was 12584, checked in by techene, 4 years ago

OCE/DOM/domqe.F90: remove dom_qe_r3c call from dom_qe_sf_nxt and make it public, OCE/steplf.F90: add dom_qe_r3c and update r3 factors (bug)
File size: 42.0 KB
Line 
1MODULE domqe
2   !!======================================================================
3   !!                       ***  MODULE domqe   ***
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   !!            4.x  ! 2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dom_qe_init     : define initial vertical scale factors, depths and column thickness
16   !!   dom_qe_sf_nxt   : Compute next vertical scale factors
17   !!   dom_qe_sf_update   : Swap vertical scale factors and update the vertical grid
18   !!   dom_qe_interpol : Interpolate vertical scale factors from one grid point to another
19   !!   dom_qe_r3c      : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points
20   !!   dom_qe_rst      : read/write restart file
21   !!   dom_qe_ctl      : Check the vvl options
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers
24   USE phycst          ! physical constant
25   USE dom_oce         ! ocean space and time domain
26   USE dynadv  , ONLY : ln_dynadv_vec
27   USE isf_oce         ! iceshelf cavities
28   USE sbc_oce         ! ocean surface boundary condition
29   USE wet_dry         ! wetting and drying
30   USE usrdef_istate   ! user defined initial state (wad only)
31   USE restart         ! ocean restart
32   !
33   USE in_out_manager  ! I/O manager
34   USE iom             ! I/O manager library
35   USE lib_mpp         ! distributed memory computing library
36   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
37   USE timing          ! Timing
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC  dom_qe_init       ! called by domain.F90
43   PUBLIC  dom_qe_zgr        ! called by isfcpl.F90
44   PUBLIC  dom_qe_sf_nxt     ! called by steplf.F90
45   PUBLIC  dom_qe_sf_update  ! called by steplf.F90
46   PUBLIC  dom_qe_interpol   ! called by dynnxt.F90
47   PUBLIC  dom_qe_r3c        ! called by steplf.F90
48
49   !                                                      !!* Namelist nam_vvl
50   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
51   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
52   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
53   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
54   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
55   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
56   !                                                       ! conservation: not used yet
57   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
58   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
59   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
60   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
61   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
62
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
64
65   !! * Substitutions
66#  include "do_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
69   !! $Id: domvvl.F90 12377 2020-02-12 14:39:06Z acc $
70   !! Software governed by the CeCILL license (see ./LICENSE)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE dom_qe_init( Kbb, Kmm, Kaa )
75      !!----------------------------------------------------------------------
76      !!                ***  ROUTINE dom_qe_init  ***
77      !!
78      !! ** Purpose :  Initialization of all scale factors, depths
79      !!               and water column heights
80      !!
81      !! ** Method  :  - use restart file and/or initialize
82      !!               - interpolate scale factors
83      !!
84      !! ** Action  : - e3t_(n/b)
85      !!              - Regrid: e3[u/v](:,:,:,Kmm)
86      !!                        e3[u/v](:,:,:,Kmm)
87      !!                        e3w(:,:,:,Kmm)
88      !!                        e3[u/v]w_b
89      !!                        e3[u/v]w_n
90      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
91      !!              - h(t/u/v)_0
92      !!
93      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
94      !!----------------------------------------------------------------------
95      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
96      !
97      IF(lwp) WRITE(numout,*)
98      IF(lwp) WRITE(numout,*) 'dom_qe_init : Variable volume activated'
99      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
100      !
101      CALL dom_qe_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
102      !
103      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
104      CALL dom_qe_rst( nit000, Kbb, Kmm, 'READ' )
105      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
106      !
107      CALL dom_qe_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
108      !
109      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
110         CALL iom_set_rstw_var_active('e3t_b')
111         CALL iom_set_rstw_var_active('e3t_n')
112      ENDIF
113      !
114   END SUBROUTINE dom_qe_init
115
116
117   SUBROUTINE dom_qe_zgr(Kbb, Kmm, Kaa)
118      !!----------------------------------------------------------------------
119      !!                ***  ROUTINE dom_qe_init  ***
120      !!
121      !! ** Purpose :  Interpolation of all scale factors,
122      !!               depths and water column heights
123      !!
124      !! ** Method  :  - interpolate scale factors
125      !!
126      !! ** Action  : - e3t_(n/b)
127      !!              - Regrid: e3(u/v)_n
128      !!                        e3(u/v)_b
129      !!                        e3w_n
130      !!                        e3(u/v)w_b
131      !!                        e3(u/v)w_n
132      !!                        gdept_n, gdepw_n and gde3w_n
133      !!              - h(t/u/v)_0
134      !!
135      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
136      !!----------------------------------------------------------------------
137      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
138      !!----------------------------------------------------------------------
139      INTEGER ::   ji, jj, jk
140      INTEGER ::   ii0, ii1, ij0, ij1
141      REAL(wp)::   zcoef
142      !!----------------------------------------------------------------------
143      !
144      !                    !== Set of all other vertical scale factors  ==!  (now and before)
145      !                                ! Horizontal interpolation of e3t
146      CALL dom_qe_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
147      CALL dom_qe_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
148      !
149      DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
150         e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) )   ! Kbb time level
151         e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) )
152         e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) )
153         e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )   ! Kmm time level
154         e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) )
155         e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) )
156         e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )
157      END DO
158      !
159      DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
160         !                                   ! The ratio does not have to be masked at w-level
161         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
162         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
163         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
164         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
165         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
166         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
167      END DO
168      !
169      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
170      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
171      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
172      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
173      !
174      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
175      IF( ln_isf ) THEN          !** IceShelF cavities
176      !                             ! to be created depending of the new names in isf
177      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
178      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
179!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
180         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
181         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
182         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
183         DO jk = 2, jpk
184            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
185                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
186            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
187                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
188            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
189            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
190                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
191            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
192                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
193         END DO
194         !
195      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
196         !
197!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
198         DO jk = 1, jpk
199            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
200            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
201            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
202            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
203            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
204         END DO
205         !
206      ENDIF
207      !
208      !                    !==  thickness of the water column  !!   (ocean portion only)
209      ht(:,:)     = ht_0(:,:)           + ssh(:,:,Kmm)
210      hu(:,:,Kbb) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kbb) )
211      hu(:,:,Kmm) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kmm) )
212      hv(:,:,Kbb) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kbb) )
213      hv(:,:,Kmm) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kmm) )
214      !
215      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
216      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
217      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
218      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
219      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
220      !
221   END SUBROUTINE dom_qe_zgr
222
223
224   SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
225      !!----------------------------------------------------------------------
226      !!                ***  ROUTINE dom_qe_sf_nxt  ***
227      !!
228      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
229      !!                 tranxt and dynspg routines
230      !!
231      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
232      !!
233      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
234      !!               - tilde_e3t_a: after increment of vertical scale factor
235      !!                              in z_tilde case
236      !!               - e3(t/u/v)_a
237      !!
238      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
239      !!----------------------------------------------------------------------
240      INTEGER, INTENT( in )           ::   kt             ! time step
241      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
242      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
243      !
244      INTEGER                ::   ji, jj, jk            ! dummy loop indices
245      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
246      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
247      LOGICAL                ::   ll_do_bclinic         ! local logical
248      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
249      !!----------------------------------------------------------------------
250      !
251      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
252      !
253      IF( ln_timing )   CALL timing_start('dom_qe_sf_nxt')
254      !
255      IF( kt == nit000 ) THEN
256         IF(lwp) WRITE(numout,*)
257         IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors'
258         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
259      ENDIF
260
261
262      ! ******************************* !
263      ! After acale factors at t-points !
264      ! ******************************* !
265      !                                                ! --------------------------------------------- !
266      !                                                ! z_star coordinate and barotropic z-tilde part !
267      !                                                ! --------------------------------------------- !
268      !
269      !
270      ! *********************************** !
271      ! After scale factors at u- v- points !
272      ! *********************************** !
273      !
274      DO jk = 1, jpkm1
275         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )
276         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )
277         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )
278      END DO
279      !
280      ! *********************************** !
281      ! After depths at u- v points         !
282      ! *********************************** !
283      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
284      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
285      !                                        ! Inverse of the local depth
286      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
287      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
288      !
289      IF( ln_timing )   CALL timing_stop('dom_qe_sf_nxt')
290      !
291   END SUBROUTINE dom_qe_sf_nxt
292
293
294   SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
295      !!----------------------------------------------------------------------
296      !!                ***  ROUTINE dom_qe_sf_update  ***
297      !!
298      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
299      !!               compute all depths and related variables for next time step
300      !!               write outputs and restart file
301      !!
302      !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
303      !!               - recompute depths and water height fields
304      !!
305      !! ** Action  :  - Recompute:
306      !!                    e3(u/v)_b
307      !!                    e3w(:,:,:,Kmm)
308      !!                    e3(u/v)w_b
309      !!                    e3(u/v)w_n
310      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
311      !!                    h(u/v) and h(u/v)r
312      !!
313      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
314      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
315      !!----------------------------------------------------------------------
316      INTEGER, INTENT( in ) ::   kt              ! time step
317      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
318      !
319      INTEGER  ::   ji, jj, jk   ! dummy loop indices
320      REAL(wp) ::   zcoef        ! local scalar
321      !!----------------------------------------------------------------------
322      !
323      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
324      !
325      IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
326      !
327      IF( kt == nit000 )   THEN
328         IF(lwp) WRITE(numout,*)
329         IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
330         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
331      ENDIF
332      !
333      ! Compute all missing vertical scale factor and depths
334      ! ====================================================
335      ! Horizontal scale factor interpolations
336      ! --------------------------------------
337      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
338      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
339
340
341      ! Scale factor computation
342      DO jk = 1, jpk             ! Horizontal interpolation
343         e3f(:,:,jk)      =  e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) )   ! Kmm time level
344         !                       ! Vertical interpolation
345         !                                   ! The ratio does not have to be masked at w-level
346         e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
347         e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
348         e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
349         e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
350         e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
351         e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
352      END DO
353
354
355      IF( ln_isf ) THEN          !** IceShelF cavities
356      !                             ! to be created depending of the new names in isf
357      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
358      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
359!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
360         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
361         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
362         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
363         DO jk = 2, jpk
364            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
365                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
366            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
367                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
368            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
369            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
370                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
371            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
372                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
373         END DO
374         !
375      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
376         !
377!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
378         DO jk = 1, jpk
379            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
380            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
381            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
382            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
383            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
384         END DO
385         !
386      ENDIF
387
388      ! Local depth and Inverse of the local depth of the water
389      ! -------------------------------------------------------
390      !
391      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
392
393      ! write restart file
394      ! ==================
395      IF( lrst_oce  )   CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )
396      !
397      IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
398      !
399   END SUBROUTINE dom_qe_sf_update
400
401
402   SUBROUTINE dom_qe_interpol( pe3_in, pe3_out, pout )
403      !!---------------------------------------------------------------------
404      !!                  ***  ROUTINE dom_qe_interpol  ***
405      !!
406      !! ** Purpose :   interpolate scale factors from one grid point to another
407      !!
408      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
409      !!                - horizontal interpolation: grid cell surface averaging
410      !!                - vertical interpolation: simple averaging
411      !!----------------------------------------------------------------------
412      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
413      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
414      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
415      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
416      !
417      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
418      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
419      !!----------------------------------------------------------------------
420      !
421      IF(ln_wd_il) THEN
422        zlnwd = 1.0_wp
423      ELSE
424        zlnwd = 0.0_wp
425      END IF
426      !
427      SELECT CASE ( pout )    !==  type of interpolation  ==!
428         !
429      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
430         DO_3D_10_10( 1, jpk )
431            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
432               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
433               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
434         END_3D
435         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'U', 1._wp )
436         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
437         !
438      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
439         DO_3D_10_10( 1, jpk )
440            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
441               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
442               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
443         END_3D
444         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'V', 1._wp )
445         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
446         !
447      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
448         DO_3D_10_10( 1, jpk )
449            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
450               &                       *    r1_e1e2f(ji,jj)                                                  &
451               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
452               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
453         END_3D
454         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'F', 1._wp )
455         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
456         !
457      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
458         !
459         !zlnwd = 1.0_wp
460         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
461         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
462!!gm BUG? use here wmask in case of ISF ?  to be checked
463         DO jk = 2, jpk
464            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
465               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
466               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
467               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
468         END DO
469         !
470      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
471         !
472         !zlnwd = 1.0_wp
473         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
474         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
475!!gm BUG? use here wumask in case of ISF ?  to be checked
476         DO jk = 2, jpk
477            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
478               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
479               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
480               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
481         END DO
482         !
483      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
484         !
485         !zlnwd = 1.0_wp
486         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
487         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
488!!gm BUG? use here wvmask in case of ISF ?  to be checked
489         DO jk = 2, jpk
490            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
491               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
492               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
493               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
494         END DO
495      END SELECT
496      !
497   END SUBROUTINE dom_qe_interpol
498
499
500   SUBROUTINE dom_qe_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
501      !!---------------------------------------------------------------------
502      !!                   ***  ROUTINE r3c  ***
503      !!
504      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
505      !!
506      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
507      !!                   Vector Form : surface weighted averaging
508      !!                   Flux   Form : simple           averaging
509      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
510      !!----------------------------------------------------------------------
511      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
512      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
513      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
514      !
515      INTEGER ::   ji, jj   ! dummy loop indices
516      !!----------------------------------------------------------------------
517      !
518      !
519      pr3t(:,:) = pssh(:,:) * r1_ht_0(:,:)   !==  ratio at t-point  ==!
520      !
521      !
522      !                                      !==  ratio at u-,v-point  ==!
523      !
524      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
525         DO_2D_11_11
526            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
527               &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
528            pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
529               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
530         END_2D
531      ELSE                                         !- Flux Form   (simple averaging)
532         DO_2D_11_11
533            pr3u(ji,jj) = 0.5_wp * (  pssh(ji  ,jj) + pssh(ji+1,jj)  ) * r1_hu_0(ji,jj)
534            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj  ) + pssh(ji,jj+1)  ) * r1_hv_0(ji,jj)
535         END_2D
536      ENDIF
537      !
538      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
539         CALL lbc_lnk_multi( 'dom_qe_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
540         !
541         !
542      ELSE                                   !==  ratio at f-point  ==!
543         !
544         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
545            DO_2D_01_01                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
546               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
547                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
548                  &                     + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
549                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
550            END_2D
551         ELSE                                      !- Flux Form   (simple averaging)
552            DO_2D_01_01                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
553               pr3f(ji,jj) = 0.25_wp * (  pssh(ji  ,jj  ) + pssh(ji+1,jj  )  &
554                  &                     + pssh(ji  ,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
555            END_2D
556         ENDIF
557         !                                                 ! lbc on ratio at u-,v-,f-points
558         CALL lbc_lnk_multi( 'dom_qe_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
559         !
560      ENDIF
561      !
562   END SUBROUTINE dom_qe_r3c
563
564
565   SUBROUTINE dom_qe_rst( kt, Kbb, Kmm, cdrw )
566      !!---------------------------------------------------------------------
567      !!                   ***  ROUTINE dom_qe_rst  ***
568      !!
569      !! ** Purpose :   Read or write VVL file in restart file
570      !!
571      !! ** Method  :   use of IOM library
572      !!                if the restart does not contain vertical scale factors,
573      !!                they are set to the _0 values
574      !!                if the restart does not contain vertical scale factors increments (z_tilde),
575      !!                they are set to 0.
576      !!----------------------------------------------------------------------
577      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
578      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
579      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
580      !
581      INTEGER ::   ji, jj, jk
582      INTEGER ::   id1, id2     ! local integers
583      !!----------------------------------------------------------------------
584      !
585      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
586         !                                   ! ===============
587         IF( ln_rstart ) THEN                   !* Read the restart file
588            CALL rst_read_open                  !  open the restart file if necessary
589            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
590            !
591            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
592            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
593            !
594            !                             ! --------- !
595            !                             ! all cases !
596            !                             ! --------- !
597            !
598            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
599               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
600               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
601               ! needed to restart if land processor not computed
602               IF(lwp) write(numout,*) 'dom_qe_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
603               WHERE ( tmask(:,:,:) == 0.0_wp )
604                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
605                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
606               END WHERE
607               IF( neuler == 0 ) THEN
608                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
609               ENDIF
610            ELSE IF( id1 > 0 ) THEN
611               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
612               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
613               IF(lwp) write(numout,*) 'neuler is forced to 0'
614               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
615               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
616               neuler = 0
617            ELSE IF( id2 > 0 ) THEN
618               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
619               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
620               IF(lwp) write(numout,*) 'neuler is forced to 0'
621               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
622               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
623               neuler = 0
624            ELSE
625               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
626               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
627               IF(lwp) write(numout,*) 'neuler is forced to 0'
628               DO jk = 1, jpk
629                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
630                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
631                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
632               END DO
633               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
634               neuler = 0
635            ENDIF
636            !
637         ELSE                                   !* Initialize at "rest"
638            !
639            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
640               !
641               IF( cn_cfg == 'wad' ) THEN
642                  ! Wetting and drying test case
643                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
644                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
645                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
646                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
647                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
648               ELSE
649                  ! if not test case
650                  ssh(:,:,Kmm) = -ssh_ref
651                  ssh(:,:,Kbb) = -ssh_ref
652
653                  DO_2D_11_11
654                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
655                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
656                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
657                     ENDIF
658                  END_2D
659               ENDIF !If test case else
660
661               ! Adjust vertical metrics for all wad
662               DO jk = 1, jpk
663                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )
664               END DO
665               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
666
667               DO ji = 1, jpi
668                  DO jj = 1, jpj
669                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
670                       CALL ctl_stop( 'dom_qe_rst: ht_0 must be positive at potentially wet points' )
671                     ENDIF
672                  END DO
673               END DO
674               !
675            ELSE
676               !
677               ! Just to read set ssh in fact, called latter once vertical grid
678               ! is set up:
679!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
680!               !
681!               DO jk=1,jpk
682!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
683!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
684!               END DO
685!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
686                ssh(:,:,Kmm)=0._wp
687                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
688                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
689               !
690            ENDIF           ! end of ll_wd edits
691            !
692         ENDIF
693         !
694      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
695         !                                   ! ===================
696         IF(lwp) WRITE(numout,*) '---- dom_qe_rst ----'
697         IF( lwxios ) CALL iom_swap(      cwxios_context          )
698         !                                           ! --------- !
699         !                                           ! all cases !
700         !                                           ! --------- !
701         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
702         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
703         !
704         IF( lwxios ) CALL iom_swap(      cxios_context          )
705      ENDIF
706      !
707   END SUBROUTINE dom_qe_rst
708
709
710   SUBROUTINE dom_qe_ctl
711      !!---------------------------------------------------------------------
712      !!                  ***  ROUTINE dom_qe_ctl  ***
713      !!
714      !! ** Purpose :   Control the consistency between namelist options
715      !!                for vertical coordinate
716      !!----------------------------------------------------------------------
717      INTEGER ::   ioptio, ios
718      !!
719      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
720         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
721         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
722      !!----------------------------------------------------------------------
723      !
724      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
725901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
726      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
727902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
728      IF(lwm) WRITE ( numond, nam_vvl )
729      !
730      IF(lwp) THEN                    ! Namelist print
731         WRITE(numout,*)
732         WRITE(numout,*) 'dom_qe_ctl : choice/control of the variable vertical coordinate'
733         WRITE(numout,*) '~~~~~~~~~~~'
734         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
735         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
736         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
737         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
738         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
739         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
740         WRITE(numout,*) '      !'
741         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
742         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
743         IF( ln_vvl_ztilde_as_zstar ) THEN
744            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
745            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
746            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
747            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
748            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
749            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
750         ELSE
751            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
752            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
753         ENDIF
754         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
755      ENDIF
756      !
757      ioptio = 0                      ! Parameter control
758      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
759      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
760      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
761      IF( ln_vvl_layer           )   ioptio = ioptio + 1
762      !
763      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
764      !
765      IF(lwp) THEN                   ! Print the choice
766         WRITE(numout,*)
767         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
768         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
769         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
770         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
771      ENDIF
772      !
773#if defined key_agrif
774      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
775#endif
776      !
777   END SUBROUTINE dom_qe_ctl
778
779   !!======================================================================
780END MODULE domqe
Note: See TracBrowser for help on using the repository browser.