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

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

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