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/trunk/src/OCE/DOM – NEMO

source: NEMO/trunk/src/OCE/DOM/domqco.F90 @ 15717

Last change on this file since 15717 was 15073, checked in by clem, 3 years ago

nn_hls=2 and debug mode: make agrif work and be repro

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