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

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

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