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

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

dynatfQCO.F90, stepLF.F90 : fixed (remove pe3. from dyn_atf_qco input arguments), all : remove e3. tables and include gurvan's feedbacks

File size: 40.9 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   !!       qe_rst_read : 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!!st   PUBLIC  dom_qe_sf_nxt     ! called by steplf.F90
45!!st   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 qe_rst_read( nit000, Kbb, Kmm )
106!!st      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! !!st
151!       DO jk = 1, jpkm1                    ! Horizontal interpolation of e3t
152!          e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) * tmask(:,:,jk) )   ! Kbb time level
153!          e3u(:,:,jk,Kbb) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) * umask(:,:,jk) )
154!          e3v(:,:,jk,Kbb) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) * vmask(:,:,jk) )
155!          e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) * tmask(:,:,jk) )   ! Kmm time level
156!          e3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) * umask(:,:,jk) )
157!          e3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) * vmask(:,:,jk) )
158!          e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + r3f(:,:)     * fmask(:,:,jk) )
159!       END DO
160!       !
161!       DO jk = 1, jpk                      ! Vertical interpolation of e3t,u,v
162!          !                                   ! The ratio does not have to be masked at w-level
163!          e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
164!          e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
165!          e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
166!          e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
167!          e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
168!          e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
169!       END DO
170!       !
171!       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
172!       e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
173!       e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
174!       e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
175!!st end
176      !
177      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
178      IF( ln_isf ) THEN          !** IceShelF cavities
179      !                             ! to be created depending of the new names in isf
180      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
181      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
182!!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
183         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
184         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
185         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
186         DO jk = 2, jpk
187            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
188                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
189            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
190                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
191            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
192            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
193                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
194            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
195                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
196         END DO
197         !
198      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
199         !
200!!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
201         DO jk = 1, jpk
202            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
203            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
204            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
205            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
206            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
207         END DO
208         !
209      ENDIF
210      !
211      !                    !==  thickness of the water column  !!   (ocean portion only)
212      ht(:,:)     = ht_0(:,:)           + ssh(:,:,Kmm)
213      hu(:,:,Kbb) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kbb) )
214      hu(:,:,Kmm) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kmm) )
215      hv(:,:,Kbb) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kbb) )
216      hv(:,:,Kmm) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kmm) )
217      !
218      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
219      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
220      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
221      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
222      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
223      !
224   END SUBROUTINE dom_qe_zgr
225
226! !!st
227!    SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
228!       !!----------------------------------------------------------------------
229!       !!                ***  ROUTINE dom_qe_sf_nxt  ***
230!       !!
231!       !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
232!       !!                 tranxt and dynspg routines
233!       !!
234!       !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
235!       !!
236!       !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
237!       !!               - tilde_e3t_a: after increment of vertical scale factor
238!       !!                              in z_tilde case
239!       !!               - e3(t/u/v)_a
240!       !!
241!       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
242!       !!----------------------------------------------------------------------
243!       INTEGER, INTENT( in )           ::   kt             ! time step
244!       INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
245!       INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
246!       !
247!       INTEGER                ::   ji, jj, jk            ! dummy loop indices
248!       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
249!       REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
250!       LOGICAL                ::   ll_do_bclinic         ! local logical
251!       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
252!       !!----------------------------------------------------------------------
253!       !
254!       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
255!       !
256!       IF( ln_timing )   CALL timing_start('dom_qe_sf_nxt')
257!       !
258!       IF( kt == nit000 ) THEN
259!          IF(lwp) WRITE(numout,*)
260!          IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors'
261!          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
262!       ENDIF
263!
264!
265!       ! ******************************* !
266!       ! After acale factors at t-points !
267!       ! ******************************* !
268!       !                                                ! --------------------------------------------- !
269!       !                                                ! z_star coordinate and barotropic z-tilde part !
270!       !                                                ! --------------------------------------------- !
271!       !
272!       !
273!       ! *********************************** !
274!       ! After scale factors at u- v- points !
275!       ! *********************************** !
276!       !
277!       DO jk = 1, jpkm1
278!          e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) * tmask(:,:,jk) )
279!          e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) * umask(:,:,jk) )
280!          e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) * vmask(:,:,jk) )
281!       END DO
282!       !
283!       ! *********************************** !
284!       ! After depths at u- v points         !
285!       ! *********************************** !
286!       hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
287!       hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
288!       !                                        ! Inverse of the local depth
289!       r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
290!       r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
291!       !
292!       IF( ln_timing )   CALL timing_stop('dom_qe_sf_nxt')
293!       !
294!    END SUBROUTINE dom_qe_sf_nxt
295!!st end
296
297   SUBROUTINE dom_h_nxt( kt, Kbb, Kmm, Kaa, kcall )
298      !!----------------------------------------------------------------------
299      !!                ***  ROUTINE dom_qe_sf_nxt  ***
300      !!
301      !! ** Purpose :  - compute the after water heigh used in tra_zdf, dynnxt,
302      !!                 tranxt and dynspg routines
303      !!
304      !! ** Method  :  - z_star case:  Proportionnaly to the water column thickness.
305      !!
306      !! ** Action  :  - h(u/v) update wrt ssh/h(u/v)_0
307      !!
308      !!----------------------------------------------------------------------
309      INTEGER, INTENT( in )           ::   kt             ! time step
310      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
311      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
312      !
313      !!----------------------------------------------------------------------
314      !
315      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
316      !
317      IF( ln_timing )   CALL timing_start('dom_h_nxt')
318      !
319      IF( kt == nit000 ) THEN
320         IF(lwp) WRITE(numout,*)
321         IF(lwp) WRITE(numout,*) 'dom_h_nxt : compute after scale factors'
322         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
323      ENDIF
324      !
325      ! *********************************** !
326      ! After depths at u- v points         !
327      ! *********************************** !
328      hu(:,:,Kaa) = hu_0(:,:) * ( 1._wp + r3u(:,:,Kaa) )
329      hv(:,:,Kaa) = hv_0(:,:) * ( 1._wp + r3v(:,:,Kaa) )
330      !                                        ! Inverse of the local depth
331      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
332      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
333      !
334      IF( ln_timing )   CALL timing_stop('dom_h_nxt')
335      !
336   END SUBROUTINE dom_h_nxt
337
338! !!st
339!    SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
340!       !!----------------------------------------------------------------------
341!       !!                ***  ROUTINE dom_qe_sf_update  ***
342!       !!
343!       !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
344!       !!               compute all depths and related variables for next time step
345!       !!               write outputs and restart file
346!       !!
347!       !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
348!       !!               - recompute depths and water height fields
349!       !!
350!       !! ** Action  :  - Recompute:
351!       !!                    e3(u/v)_b
352!       !!                    e3w(:,:,:,Kmm)
353!       !!                    e3(u/v)w_b
354!       !!                    e3(u/v)w_n
355!       !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
356!       !!                    h(u/v) and h(u/v)r
357!       !!
358!       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
359!       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
360!       !!----------------------------------------------------------------------
361!       INTEGER, INTENT( in ) ::   kt              ! time step
362!       INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
363!       !
364!       INTEGER  ::   ji, jj, jk   ! dummy loop indices
365!       REAL(wp) ::   zcoef        ! local scalar
366!       !!----------------------------------------------------------------------
367!       !
368!       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
369!       !
370!       IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
371!       !
372!       IF( kt == nit000 )   THEN
373!          IF(lwp) WRITE(numout,*)
374!          IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
375!          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
376!       ENDIF
377!       !
378!       ! Compute all missing vertical scale factor and depths
379!       ! ====================================================
380!       ! Horizontal scale factor interpolations
381!       ! --------------------------------------
382!       ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
383!       ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
384!
385!
386!       ! Scale factor computation
387!       DO jk = 1, jpk             ! Horizontal interpolation
388!          e3f(:,:,jk)      =  e3f_0(:,:,jk) * ( 1._wp + r3f(:,:) * fmask(:,:,jk) )   ! Kmm time level
389!          !                       ! Vertical interpolation
390!        !                                   ! The ratio does not have to be masked at w-level
391!          e3w (:,:,jk,Kmm) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )   ! Kmm time level
392!          e3uw(:,:,jk,Kmm) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kmm) )
393!          e3vw(:,:,jk,Kmm) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kmm) )
394!          e3w (:,:,jk,Kbb) =  e3w_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )   ! Kbb time level
395!          e3uw(:,:,jk,Kbb) = e3uw_0(:,:,jk) * ( 1._wp + r3u(:,:,Kbb) )
396!          e3vw(:,:,jk,Kbb) = e3vw_0(:,:,jk) * ( 1._wp + r3v(:,:,Kbb) )
397!       END DO
398!
399!
400!       IF( ln_isf ) THEN          !** IceShelF cavities
401!       !                             ! to be created depending of the new names in isf
402!       !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
403!       !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
404! !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
405!          gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
406!          gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
407!          gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
408!          DO jk = 2, jpk
409!             gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
410!                               + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
411!             gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
412!                               + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
413!             gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
414!             gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
415!                               + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
416!             gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
417!                               + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
418!          END DO
419!          !
420!       ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
421!          !
422! !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
423!          DO jk = 1, jpk
424!             gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
425!             gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
426!             gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
427!             gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
428!             gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
429!          END DO
430!          !
431!       ENDIF
432!
433!       ! Local depth and Inverse of the local depth of the water
434!       ! -------------------------------------------------------
435!       !
436!       ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
437!
438!       ! write restart file
439!       ! ==================
440!       IF( lrst_oce  )   CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )
441!       !
442!       IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
443!       !
444!    END SUBROUTINE dom_qe_sf_update
445!!st end
446
447   SUBROUTINE dom_h_update( kt, Kbb, Kmm, Kaa )
448      !!----------------------------------------------------------------------
449      !!                ***  ROUTINE dom_qe_sf_update  ***
450      !!
451      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
452      !!               compute all depths and related variables for next time step
453      !!               write outputs and restart file
454      !!
455      !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
456      !!               - recompute depths and water height fields
457      !!
458      !! ** Action  :  - Recompute:
459      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
460      !!                    h(u/v) and h(u/v)r
461      !!
462      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
463      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
464      !!----------------------------------------------------------------------
465      INTEGER, INTENT( in ) ::   kt              ! time step
466      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
467      !
468      INTEGER  ::   ji, jj, jk   ! dummy loop indices
469      REAL(wp) ::   zcoef        ! local scalar
470      !!----------------------------------------------------------------------
471      !
472      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
473      !
474      IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
475      !
476      IF( kt == nit000 )   THEN
477         IF(lwp) WRITE(numout,*)
478         IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
479         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
480      ENDIF
481      !
482      ! Compute all missing vertical scale factor and depths
483      ! ====================================================
484      ! Horizontal scale factor interpolations
485      ! --------------------------------------
486      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
487
488      IF( ln_isf ) THEN          !** IceShelF cavities
489      !                             ! to be created depending of the new names in isf
490      !                             ! it should be something like that :  (with h_isf = thickness of iceshelf)
491      !                             !  in fact currently, h_isf(:,:) is called : risfdep(:,:)
492   !!gm - depth idea 0  : just realize that mask is not needed ===>>>> with ISF, rescale all grid point position below ISF : no mask !
493         gdept(:,:,1,Kmm) = gdept_0(:,:,1) * ( 1._wp + r3t(:,:,Kmm) )
494         gdepw(:,:,1,Kmm) = 0._wp                            ! Initialized to zero one for all
495         gde3w(:,:,1)     = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
496         DO jk = 2, jpk
497            gdept(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
498                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
499            gdepw(:,:,jk,Kmm) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
500                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kmm) )
501            gde3w(:,:,jk)     = gdept(:,:,jk,Kmm) - ssh(:,:,Kmm)
502            gdept(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdept_0(:,:,jk) )            &
503                              + MAX(   0._wp    , gdept_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
504            gdepw(:,:,jk,Kbb) = MIN( risfdep(:,:) , gdepw_0(:,:,jk) )    &
505                              + MAX(   0._wp    , gdepw_0(:,:,jk)-risfdep(:,:) ) * ( 1._wp + r3t(:,:,Kbb) )
506         END DO
507         !
508      ELSE                       !** No cavities   (all depth rescaled, even inside topography: no mask)
509         !
510   !!gm idea 0 :   just realize that mask is not needed ===>>>> without ISF, rescale all grid point position : no mask !
511         DO jk = 1, jpk
512            gdept(:,:,jk,Kmm) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
513            gdepw(:,:,jk,Kmm) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kmm) )
514            gde3w(:,:,jk)     = gdept  (:,:,jk,Kmm)       - ssh(:,:,Kmm)
515            gdept(:,:,jk,Kbb) = gdept_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
516            gdepw(:,:,jk,Kbb) = gdepw_0(:,:,jk) * ( 1._wp + r3t(:,:,Kbb) )
517         END DO
518         !
519      ENDIF
520
521      ! Local depth and Inverse of the local depth of the water
522      ! -------------------------------------------------------
523      !
524      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
525
526      ! write restart file
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 qe_rst_read( kt, Kbb, Kmm )
599      !!---------------------------------------------------------------------
600      !!                   ***  ROUTINE qe_rst_read  ***
601      !!
602      !! ** Purpose :   Read ssh in restart file
603      !!
604      !! ** Method  :   use of IOM library
605      !!                if the restart does not contain ssh,
606      !!                it is set to the _0 values.
607      !!----------------------------------------------------------------------
608      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
609      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
610      !
611      INTEGER ::   ji, jj, jk
612      INTEGER ::   id1, id2     ! local integers
613      !!----------------------------------------------------------------------
614      !
615         IF( ln_rstart ) THEN                   !* Read the restart file
616            CALL rst_read_open                  !  open the restart file if necessary
617            !
618            id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. )
619            id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. )
620            !
621            !                             ! --------- !
622            !                             ! all cases !
623            !                             ! --------- !
624            !
625            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
626               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    )
627               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
628               ! needed to restart if land processor not computed
629               IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files'
630               WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh when it was required on e3
631                  ssh(:,:,Kmm) = 0._wp
632                  ssh(:,:,Kbb) = 0._wp
633               END WHERE
634               IF( neuler == 0 ) THEN
635                  ssh(:,:,Kbb) = ssh(:,:,Kmm)
636               ENDIF
637            ELSE IF( id1 > 0 ) THEN
638               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files'
639               IF(lwp) write(numout,*) 'sshn set equal to sshb.'
640               IF(lwp) write(numout,*) 'neuler is forced to 0'
641               CALL iom_get( numror, jpdom_autoglo, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
642               ssh(:,:,Kmm) = ssh(:,:,Kbb)
643               neuler = 0
644            ELSE IF( id2 > 0 ) THEN
645               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files'
646               IF(lwp) write(numout,*) 'sshb set equal to sshn.'
647               IF(lwp) write(numout,*) 'neuler is forced to 0'
648               CALL iom_get( numror, jpdom_autoglo, 'sshn', ssh(:,:,Kmm), ldxios = lrxios )
649               ssh(:,:,Kbb) = ssh(:,:,Kmm)
650               neuler = 0
651            ELSE
652               IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file'
653               IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero'
654               IF(lwp) write(numout,*) 'neuler is forced to 0'
655               ssh(:,:,:) = 0._wp
656               neuler = 0
657            ENDIF
658            !
659         ELSE                                   !* Initialize at "rest"
660            !
661            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
662               !
663               IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case
664                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
665                  ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
666                  ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)
667                  uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb)
668                  vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb)
669               ELSE                                  ! if not test case
670                  ssh(:,:,Kmm) = -ssh_ref
671                  ssh(:,:,Kbb) = -ssh_ref
672                  !
673                  DO_2D_11_11
674                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
675                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
676                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
677                     ENDIF
678                  END_2D
679               ENDIF
680
681               DO ji = 1, jpi
682                  DO jj = 1, jpj
683                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
684                       CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' )
685                     ENDIF
686                  END DO
687               END DO
688               !
689            ELSE
690               !
691               ! Just to read set ssh in fact, called latter once vertical grid
692               ! is set up:
693!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
694!               !
695                ssh(:,:,:) = 0._wp
696               !
697            ENDIF           ! end of ll_wd edits
698            !
699         ENDIF
700      !
701   END SUBROUTINE qe_rst_read
702
703
704   SUBROUTINE dom_qe_ctl
705      !!---------------------------------------------------------------------
706      !!                  ***  ROUTINE dom_qe_ctl  ***
707      !!
708      !! ** Purpose :   Control the consistency between namelist options
709      !!                for vertical coordinate
710      !!----------------------------------------------------------------------
711      INTEGER ::   ioptio, ios
712      !!
713      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
714         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
715         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
716      !!----------------------------------------------------------------------
717      !
718      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
719901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
720      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
721902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
722      IF(lwm) WRITE ( numond, nam_vvl )
723      !
724      IF(lwp) THEN                    ! Namelist print
725         WRITE(numout,*)
726         WRITE(numout,*) 'dom_qe_ctl : choice/control of the variable vertical coordinate'
727         WRITE(numout,*) '~~~~~~~~~~~'
728         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
729         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
730         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
731         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
732         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
733         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
734         WRITE(numout,*) '      !'
735         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
736         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
737         IF( ln_vvl_ztilde_as_zstar ) THEN
738            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
739            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
740            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
741            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
742            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
743            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
744         ELSE
745            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
746            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
747         ENDIF
748         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
749      ENDIF
750      !
751      ioptio = 0                      ! Parameter control
752      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
753      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
754      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
755      IF( ln_vvl_layer           )   ioptio = ioptio + 1
756      !
757      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
758      !
759      IF(lwp) THEN                   ! Print the choice
760         WRITE(numout,*)
761         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
762         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
763         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
764         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
765      ENDIF
766      !
767#if defined key_agrif
768      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
769#endif
770      !
771   END SUBROUTINE dom_qe_ctl
772
773   !!======================================================================
774END MODULE domqe
Note: See TracBrowser for help on using the repository browser.