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 @ 12644

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

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