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

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

#2385 restart with qco bug fix for euler 1st and RK3
File size: 14.4 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 steplf.F90
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
66   SUBROUTINE dom_qco_init( Kbb, Kmm, Kaa )
67      !!----------------------------------------------------------------------
68      !!                ***  ROUTINE dom_qco_init  ***
69      !!
70      !! ** Purpose :  Initialization of all ssh. to h._0 ratio
71      !!
72      !! ** Method  :  - use restart file and/or initialize
73      !!               - compute ssh. to h._0 ratio
74      !!
75      !! ** Action  : - r3(t/u/v)_b
76      !!              - r3(t/u/v/f)_n
77      !!
78      !!----------------------------------------------------------------------
79      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices
80      !!----------------------------------------------------------------------
81      !
82      IF(lwp) WRITE(numout,*)
83      IF(lwp) WRITE(numout,*) 'dom_qco_init : Variable volume activated'
84      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
85      !
86      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer)
87      !
88      CALL dom_qco_zgr( Kbb, Kmm )            ! interpolation scale factor, depth and water column
89      !
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
97      !
98   END SUBROUTINE dom_qco_init
99
100
101   SUBROUTINE dom_qco_zgr( Kbb, Kmm )
102      !!----------------------------------------------------------------------
103      !!                ***  ROUTINE dom_qco_init  ***
104      !!
105      !! ** Purpose :  Initialization of all r3. = ssh./h._0 ratios
106      !!
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
109      !!
110      !! ** Action  : - r3(t/u/v)(Kbb)
111      !!              - r3(t/u/v/f)(Kmm)
112      !!----------------------------------------------------------------------
113      INTEGER, INTENT(in) ::   Kbb, Kmm   ! time level indices
114      !!----------------------------------------------------------------------
115      !
116      !                    !== Set of all other vertical scale factors  ==!  (now and before)
117      !                                ! Horizontal interpolation of e3t
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
122      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
123      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
124#endif
125      !
126   END SUBROUTINE dom_qco_zgr
127
128
129   SUBROUTINE dom_qco_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
130      !!---------------------------------------------------------------------
131      !!                   ***  ROUTINE r3c  ***
132      !!
133      !! ** Purpose :   compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
134      !!
135      !! ** Method  : - compute the ssh at u- and v-points (f-point optional)
136      !!                   Vector Form : surface weighted averaging
137      !!                   Flux   Form : simple           averaging
138      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
139      !!----------------------------------------------------------------------
140      REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m]
141      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-]
142      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-]
143      !
144      INTEGER ::   ji, jj   ! dummy loop indices
145      !!----------------------------------------------------------------------
146      !
147      !
148      pr3t(:,:) = pssh(:,:) * r1_ht_0(:,:)   !==  ratio at t-point  ==!
149      !
150      !
151      !                                      !==  ratio at u-,v-point  ==!
152      !
153!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging)
154#if ! defined key_qcoTest_FluxForm
155      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
156         DO_2D( 0, 0, 0, 0 )
157            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
158               &                    + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
159            pr3v(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
160               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
161         END_2D
162!!st      ELSE                                         !- Flux Form   (simple averaging)
163#else
164         DO_2D( 0, 0, 0, 0 )
165            pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj)
166            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj)
167         END_2D
168!!st      ENDIF
169#endif         
170      !
171      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only
172         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
173         !
174         !
175      ELSE                                   !==  ratio at f-point  ==!
176         !
177!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging)
178#if ! defined key_qcoTest_FluxForm
179         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
180
181            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
182               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
183                  &                     + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
184                  &                     + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
185                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
186            END_2D
187!!st         ELSE                                      !- Flux Form   (simple averaging)
188#else
189            DO_2D( 0, 0, 0, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
190               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  &
191                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj)
192            END_2D
193!!st         ENDIF
194#endif
195         !                                                 ! lbc on ratio at u-,v-,f-points
196         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
197         !
198      ENDIF
199      !
200   END SUBROUTINE dom_qco_r3c
201
202
203   SUBROUTINE qco_ctl
204      !!---------------------------------------------------------------------
205      !!                  ***  ROUTINE qco_ctl  ***
206      !!
207      !! ** Purpose :   Control the consistency between namelist options
208      !!                for vertical coordinate
209      !!----------------------------------------------------------------------
210      INTEGER ::   ioptio, ios
211      !!
212      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
213         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
214         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
215      !!----------------------------------------------------------------------
216      !
217      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
218901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
219      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
220902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
221      IF(lwm) WRITE ( numond, nam_vvl )
222      !
223      IF(lwp) THEN                    ! Namelist print
224         WRITE(numout,*)
225         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate'
226         WRITE(numout,*) '~~~~~~~~'
227         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
228         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
229         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
230         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
231         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
232         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
233         WRITE(numout,*) '      !'
234         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
235         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
236         IF( ln_vvl_ztilde_as_zstar ) THEN
237            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
238            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
239            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
240            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
241            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
242            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
243         ELSE
244            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
245            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
246         ENDIF
247         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
248      ENDIF
249      !
250      ioptio = 0                      ! Parameter control
251      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
252      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
253      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
254      IF( ln_vvl_layer           )   ioptio = ioptio + 1
255      !
256      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
257      !
258      IF(lwp) THEN                   ! Print the choice
259         WRITE(numout,*)
260         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
261         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
262         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
263         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
264      ENDIF
265      !
266#if defined key_agrif
267      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
268#endif
269      !
270   END SUBROUTINE qco_ctl
271
272   !!======================================================================
273END MODULE domqco
Note: See TracBrowser for help on using the repository browser.