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.
domqco.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domqco.F90 @ 15627

Last change on this file since 15627 was 15627, checked in by techene, 3 years ago

#2605 qco r3. ratios management optimised (RK3) and some cleanning (MLF)
File size: 17.9 KB
Line 
1MODULE domqco
2   !!======================================================================
3   !!                       ***  MODULE domqco   ***
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) add time level indices for prognostic variables
11   !!             -   !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*)
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dom_qco_init  : define initial vertical scale factors, depths and column thickness
16   !!   dom_qco_zgr   : Set ssh/h_0 ratio at t
17   !!   dom_qco_r3c   : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points
18   !!       qco_ctl   : Check the vvl options
19   !!----------------------------------------------------------------------
20   USE oce            ! ocean dynamics and tracers
21   USE phycst         ! physical constant
22   USE dom_oce        ! ocean space and time domain
23   USE dynadv  , ONLY : ln_dynadv_vec
24   USE isf_oce        ! iceshelf cavities
25   USE sbc_oce        ! ocean surface boundary condition
26   USE wet_dry        ! wetting and drying
27   USE usrdef_istate  ! user defined initial state (wad only)
28   USE restart        ! ocean restart
29   !
30   USE in_out_manager ! I/O manager
31   USE iom            ! I/O manager library
32   USE lib_mpp        ! distributed memory computing library
33   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
34   USE timing         ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC  dom_qco_init       ! called by domain.F90
40   PUBLIC  dom_qco_zgr        ! called by isfcpl.F90
41   PUBLIC  dom_qco_r3c        ! called by stpmlf.F90
42   PUBLIC  dom_qco_r3c_RK3    ! called by stprk3_stg.F90
43
44   !                                                      !!* Namelist nam_vvl
45   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
46   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
51   !                                                       ! conservation: not used yet
52   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
53   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
54   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
55   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
56   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
57
58   !! * Substitutions
59#  include "do_loop_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id: domvvl.F90 12377 2020-02-12 14:39:06Z acc $
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE dom_qco_init( Kbb, Kmm, Kaa )
68      !!----------------------------------------------------------------------
69      !!                ***  ROUTINE dom_qco_init  ***
70      !!
71      !! ** Purpose :  Initialization of all ssh. to h._0 ratio
72      !!
73      !! ** Method  :  - use restart file and/or initialize
74      !!               - compute ssh. to h._0 ratio
75      !!
76      !! ** Action  : - r3(t/u/v)_b
77      !!              - r3(t/u/v/f)_n
78      !!
79      !!----------------------------------------------------------------------
80      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
81      !!----------------------------------------------------------------------
82      !
83      IF(lwp) WRITE(numout,*)
84      IF(lwp) WRITE(numout,*) 'dom_qco_init : Variable volume activated'
85      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
86      !
87      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer)
88      !
89      CALL dom_qco_zgr( Kbb, Kmm )            ! interpolation scale factor, depth and water column
90      !
91#if defined key_agrif
92      ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a
93      ! problem for the restartability...)
94      r3t(:,:,Kaa) = r3t(:,:,Kmm)
95      r3u(:,:,Kaa) = r3u(:,:,Kmm)
96      r3v(:,:,Kaa) = r3v(:,:,Kmm)
97#endif
98      !
99   END SUBROUTINE dom_qco_init
100
101
102   SUBROUTINE dom_qco_zgr( Kbb, Kmm )
103      !!----------------------------------------------------------------------
104      !!                ***  ROUTINE dom_qco_init  ***
105      !!
106      !! ** Purpose :  Initialization of all r3. = ssh./h._0 ratios
107      !!
108      !! ** Method  :  Call domqco using Kbb and Kmm
109      !!               NB: dom_qco_zgr is called by dom_qco_init it uses ssh from ssh_init
110      !!
111      !! ** Action  : - r3(t/u/v)(Kbb)
112      !!              - r3(t/u/v/f)(Kmm)
113      !!----------------------------------------------------------------------
114      INTEGER, INTENT(in) ::   Kbb, Kmm   ! time level indices
115      !!----------------------------------------------------------------------
116      !
117      !                    !== Set of all other vertical scale factors  ==!  (now and before)
118      !                                ! Horizontal interpolation of e3t
119#if defined key_RK3
120      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
121      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )  !! needed for agrif ???
122#else
123      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
124      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
125#endif
126      ! dom_qco_r3c defines over [nn_hls, nn_hls-1, nn_hls, nn_hls-1]
127      IF( nn_hls == 2 ) CALL lbc_lnk( 'dom_qco_zgr', r3u(:,:,Kbb), 'U', 1._wp, r3v(:,:,Kbb), 'V', 1._wp, &
128         &                                           r3u(:,:,Kmm), 'U', 1._wp, r3v(:,:,Kmm), 'V', 1._wp, r3f(:,:), 'F', 1._wp )
129      !                                                                                                ! r3f is needed for agrif
130   END SUBROUTINE dom_qco_zgr
131
132
133   SUBROUTINE dom_qco_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
134      !!---------------------------------------------------------------------
135      !!                   ***  ROUTINE r3c  ***
136      !!
137      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
138      !!
139      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
140      !!                   Vector Form : surface weighted averaging
141      !!                   Flux   Form : simple           averaging
142      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
143      !!----------------------------------------------------------------------
144      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
145      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
146      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
147      !
148      INTEGER ::   ji, jj   ! dummy loop indices
149      !!----------------------------------------------------------------------
150      !
151      !
152      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
153         pr3t(ji,jj) = pssh(ji,jj) * r1_ht_0(ji,jj)   !==  ratio at t-point  ==!
154      END_2D
155      !
156      !
157      !                                      !==  ratio at u-,v-point  ==!
158      !
159      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
160         pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
161            &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
162         pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
163            &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
164      END_2D     
165      !
166      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
167         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
168         !
169         !
170      ELSE                                   !==  ratio at f-point  ==!
171         !
172         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
173            ! round brackets added to fix the order of floating point operations
174            ! needed to ensure halo 1 - halo 2 compatibility
175            pr3f(ji,jj) = 0.25_wp * (   (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )      &
176               &                         + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  )   & ! bracket for halo 1 - halo 2 compatibility
177               &                     +  (  e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)      &
178               &                         + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
179               &                    ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
180         END_2D
181         !                                                 ! lbc on ratio at u-,v-,f-points
182         IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
183         !
184      ENDIF
185      !
186   END SUBROUTINE dom_qco_r3c
187
188   
189   SUBROUTINE dom_qco_r3c_RK3( pssh, pr3t, pr3u, pr3v, pr3f )
190      !!---------------------------------------------------------------------
191      !!                   ***  ROUTINE r3c  ***
192      !!
193      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
194      !!
195      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
196      !!                   Vector Form : surface weighted averaging
197      !!                   Flux   Form : simple           averaging
198      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
199      !!----------------------------------------------------------------------
200      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
201      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
202      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
203      !
204      INTEGER ::   ji, jj   ! dummy loop indices
205      !!----------------------------------------------------------------------
206      !
207      !
208      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
209         pr3t(ji,jj) = pssh(ji,jj) * r1_ht_0(ji,jj)   !==  ratio at t-point  ==!
210      END_2D
211      !
212      !
213      !                                      !==  ratio at u-,v-point  ==!
214      !
215!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
216#if ! defined key_qcoTest_FluxForm
217      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
218      DO_2D( 0, 0, 0, 0 )
219         pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
220            &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
221         pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
222            &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
223      END_2D
224!!st      ELSE                                         !- Flux Form   (simple averaging)
225#else
226      DO_2D( 0, 0, 0, 0 )
227         pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
228         pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
229      END_2D
230!!st      ENDIF
231#endif         
232      !
233      IF( PRESENT( pr3f ) ) THEN             !==  ratio at f-point  ==!
234         !
235!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
236#if ! defined key_qcoTest_FluxForm
237         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
238
239      DO_2D( 0, 0, 0, 0 )
240         ! round brackets added to fix the order of floating point operations
241         ! needed to ensure halo 1 - halo 2 compatibility
242         pr3f(ji,jj) = 0.25_wp * (   (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )      &
243            &                         + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  )   & ! bracket for halo 1 - halo 2 compatibility
244            &                     +  (  e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)      &
245            &                         + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
246            &                    ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
247      END_2D
248!!st         ELSE                                      !- Flux Form   (simple averaging)
249#else
250      DO_2D( 0, 0, 0, 0 )
251         ! round brackets added to fix the order of floating point operations
252         ! needed to ensure halo 1 - halo 2 compatibility
253         pr3f(ji,jj) = 0.25_wp * (   (  pssh(ji,jj  ) + pssh(ji+1,jj  )  )   &
254            &                     +  (  pssh(ji,jj+1) + pssh(ji+1,jj+1)  )   & ! bracket for halo 1 - halo 2 compatibility
255            &                    ) * r1_hf_0(ji,jj)
256      END_2D
257!!st         ENDIF
258#endif
259         !
260      ENDIF
261      !
262   END SUBROUTINE dom_qco_r3c_RK3
263
264   
265   SUBROUTINE qco_ctl
266      !!---------------------------------------------------------------------
267      !!                  ***  ROUTINE qco_ctl  ***
268      !!
269      !! ** Purpose :   Control the consistency between namelist options
270      !!                for vertical coordinate
271      !!----------------------------------------------------------------------
272      INTEGER ::   ioptio, ios
273      !!
274      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
275         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
276         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
277      !!----------------------------------------------------------------------
278      !
279      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
280901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
281      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
282902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
283      IF(lwm) WRITE ( numond, nam_vvl )
284      !
285      IF(lwp) THEN                    ! Namelist print
286         WRITE(numout,*)
287         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate'
288         WRITE(numout,*) '~~~~~~~~'
289         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
290         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
291         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
292         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
293         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
294         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
295         WRITE(numout,*) '      !'
296         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
297         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
298         IF( ln_vvl_ztilde_as_zstar ) THEN
299            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
300            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
301            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
302            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
303            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
304            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
305         ELSE
306            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
307            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
308         ENDIF
309         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
310      ENDIF
311      !
312      ioptio = 0                      ! Parameter control
313      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
314      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
315      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
316      IF( ln_vvl_layer           )   ioptio = ioptio + 1
317      !
318      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
319      !
320      IF(lwp) THEN                   ! Print the choice
321         WRITE(numout,*)
322         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
323         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
324         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
325         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
326      ENDIF
327      !
328#if defined key_agrif
329      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
330#endif
331      !
332   END SUBROUTINE qco_ctl
333
334   !!======================================================================
335END MODULE domqco
Note: See TracBrowser for help on using the repository browser.