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

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

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