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

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

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