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

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

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