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.
dyncor.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dyncor.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
  • Property svn:keywords set to Id
File size: 17.6 KB
Line 
1MODULE dyncor
2   !!======================================================================
3   !!                       ***  MODULE  dyncor  ***
4   !! Ocean dynamics: Update the momentum trend with the planetary vorticity trends
5   !!======================================================================
6   !! History :  5.0  ! 2018-07  (G. Madec)  Coriolis trend for Flux Form
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dyn_cor       : Update the momentum trend with the vorticity trend
11   !!       vor_ens   : enstrophy conserving scheme       (ln_dynvor_ens=T)
12   !!       vor_ene   : energy conserving scheme          (ln_dynvor_ene=T)
13   !!       vor_een   : energy and enstrophy conserving   (ln_dynvor_een=T)
14   !!   dyn_cor_init  : set and control of the different vorticity option
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE dommsk         ! ocean mask
19   USE dynadv         ! momentum advection
20   USE trd_oce        ! trends: ocean variables
21   USE trddyn         ! trend manager: dynamics
22   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
23   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
24   !
25   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
26   USE prtctl         ! Print control
27   USE in_out_manager ! I/O manager
28   USE lib_mpp        ! MPP library
29   USE timing         ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_cor        ! routine called by step.F90
35   PUBLIC   dyn_cor_init   ! routine called by nemogcm.F90
36
37   !                                   !!* Namelist namdyn_vor: vorticity term
38   LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS)
39   LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE)
40   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT)
41   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET)
42   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN)
43   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
44   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX)
45   LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
46
47   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme
48   !                                       ! associated indices:
49   INTEGER, PUBLIC, PARAMETER ::   np_ENS = 0   ! ENS scheme
50   INTEGER, PUBLIC, PARAMETER ::   np_ENE = 1   ! ENE scheme
51   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity)
52   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t)
53   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme
54   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme
55
56   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity
57   !                               ! associated indices:
58   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary)
59   INTEGER, PUBLIC, PARAMETER ::   np_RVO = 2         ! relative vorticity
60   INTEGER, PUBLIC, PARAMETER ::   np_MET = 3         ! metric term
61   INTEGER, PUBLIC, PARAMETER ::   np_CRV = 4         ! relative + planetary (total vorticity)
62   INTEGER, PUBLIC, PARAMETER ::   np_CME = 5         ! Coriolis + metric term
63
64   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2u_2        ! = di(e2u)/2          used in T-point metric term calculation
65   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1v_2        ! = dj(e1v)/2           -        -      -       -
66   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   di_e2v_2e1e2f   ! = di(e2u)/(2*e1e2f)  used in F-point metric term calculation
67   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   dj_e1u_2e1e2f   ! = dj(e1v)/(2*e1e2f)   -        -      -       -
68   
69   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4
70   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8
71   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12
72   
73   !! * Substitutions
74#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
77   !! $Id$
78   !! Software governed by the CeCILL licence     (./LICENSE)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   SUBROUTINE dyn_cor( kt )
83      !!----------------------------------------------------------------------
84      !!
85      !! ** Purpose :   compute the lateral ocean tracer physics.
86      !!
87      !! ** Action : - Update (ua,va) with the now vorticity term trend
88      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
89      !!               and planetary vorticity trends) and send them to trd_dyn
90      !!               for futher diagnostics (l_trddyn=T)
91      !!----------------------------------------------------------------------
92      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
93      !
94      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv
95      !!----------------------------------------------------------------------
96      !
97      IF( ln_timing )   CALL timing_start('dyn_cor')
98      !
99      IF( l_trddyn ) THEN     !==  trend diagnostics case : split the added trend in two parts  ==!
100         !
101         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
102         !
103         ztrdu(:,:,:) = ua(:,:,:)            !* planetary vorticity trend (including Stokes-Coriolis force)
104         ztrdv(:,:,:) = va(:,:,:)
105
106         CALL cor_ene( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (T-pts)
107
108         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
109         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
110         CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
111         !
112         IF( n_dynadv /= np_LIN_dyn ) THEN   !* relative vorticity or metric trend (only in non-linear case)
113            ztrdu(:,:,:) = ua(:,:,:)
114            ztrdv(:,:,:) = va(:,:,:)
115            CALL cor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts)
116            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
117            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
118            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
119         ENDIF
120         !
121         DEALLOCATE( ztrdu, ztrdv )
122         !
123      ELSE              !==  Coriolis (+metric) trend added to the general trend  ==!
124         !
125         !                    !* energy conserving scheme  (T-pts)
126         IF( ln_stcor )   CALL cor_ene( kt, ntot, usd, vsd , ua, va )   !  Stokes drift
127                          CALL cor_ene( kt, ncor, un , vn  , ua, va )   ! 
128         !
129      ENDIF
130      !
131      !                       ! print sum trends (used for debugging)
132      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
133         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
134      !
135      IF( ln_timing )   CALL timing_stop('dyn_cor')
136      !
137   END SUBROUTINE dyn_cor
138
139
140   SUBROUTINE cor_ene( kt, kvor, pu, pv, pu_rhs, pv_rhs )
141      !!----------------------------------------------------------------------
142      !!                  ***  ROUTINE vor_enT  ***
143      !!
144      !! ** Purpose :   Compute the now Coriolis (+ metric term) trend and 
145      !!      add it to the general trend of the momentum equation.
146      !!
147      !! ** Method  :   Trend evaluated using now fields (centered in time)
148      !!       and t-point evaluation of vorticity (planetary and relative).
149      !!       conserves the horizontal kinetic energy.
150      !!         The general trend of momentum is increased due to the vorticity
151      !!       term which is given by:
152      !!          voru = 1/bu  mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t  mj[vn] ]
153      !!          vorv = 1/bv  mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f  mj[un] ]
154      !!       where rvor is the relative vorticity at f-point
155      !!
156      !! ** Action : - Update (ua,va) with the now vorticity term trend
157      !!----------------------------------------------------------------------
158      INTEGER                         , INTENT(in   ) ::   kt               ! ocean time-step index
159      INTEGER                         , INTENT(in   ) ::   kvor             ! total, planetary, relative, or metric
160      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv           ! now velocities
161      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! total v-trend
162      !
163      INTEGER  ::   ji, jj, jk           ! dummy loop indices
164      REAL(wp) ::   zx1, zy1, zx2, zy2   ! local scalars
165      REAL(wp), DIMENSION(jpi,jpj) ::   zwx, zwy, zwz, zwt   ! 2D workspace
166      !!----------------------------------------------------------------------
167      !
168      IF( kt == nit000 ) THEN
169         IF(lwp) WRITE(numout,*)
170         IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
171         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
172      ENDIF
173      !
174      !                                                ! ===============
175      DO jk = 1, jpkm1                                 ! Horizontal slab
176         !                                             ! ===============
177         !
178         SELECT CASE( kvor )                 !==  volume weighted vorticity considered  ==!
179         CASE ( np_COR )                           !* Coriolis (planetary vorticity)
180            zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk)
181         CASE ( np_MET )                           !* metric term
182            DO jj = 2, jpj
183               DO ji = 2, jpi
184                  zwt(ji,jj) = (   ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)   &
185                     &           - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)   ) * e3t_n(ji,jj,jk)
186               END DO
187            END DO
188         CASE ( np_CME )                           !* Coriolis + metric
189            DO jj = 2, jpj
190               DO ji = 2, jpi   ! vector opt.
191                  zwt(ji,jj) = (  ff_t(ji,jj) * e1e2t(ji,jj)                           &
192                       &        + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj)  &
193                       &        - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj)  ) * e3t_n(ji,jj,jk)
194               END DO
195            END DO
196         CASE DEFAULT                                             ! error
197            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  )
198         END SELECT
199         !
200         !                                   !==  compute and add the vorticity term trend  =!
201         DO jj = 2, jpjm1
202            DO ji = 2, jpim1   ! vector opt.
203               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)                    &
204                  &                                * (  zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) )   &
205                  &                                   + zwt(ji  ,jj) * ( pv(ji  ,jj,jk) + pv(ji  ,jj-1,jk) )   )
206                  !
207               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)                    &
208                  &                                * (  zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) )   & 
209                  &                                   + zwt(ji,jj  ) * ( pu(ji,jj  ,jk) + pu(ji-1,jj  ,jk) )   ) 
210            END DO 
211         END DO 
212         !                                             ! ===============
213      END DO                                           !   End of slab
214      !                                                ! ===============
215   END SUBROUTINE cor_ene
216
217
218   SUBROUTINE dyn_cor_init
219      !!---------------------------------------------------------------------
220      !!                  ***  ROUTINE dyn_cor_init  ***
221      !!
222      !! ** Purpose :   Control the consistency between cpp options for
223      !!              tracer advection schemes
224      !!----------------------------------------------------------------------
225      INTEGER ::   ji, jj, jk    ! dummy loop indices
226      INTEGER ::   ioptio, ios   ! local integer
227      !!
228      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   &
229         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk
230      !!----------------------------------------------------------------------
231      !
232      IF(lwp) THEN
233         WRITE(numout,*)
234         WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
235         WRITE(numout,*) '~~~~~~~~~~~~'
236      ENDIF
237      !
238      REWIND( numnam_ref )              ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
239      READ  ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
240901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
241      REWIND( numnam_cfg )              ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
242      READ  ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
243902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
244      IF(lwm) WRITE ( numond, namdyn_vor )
245      !
246      IF(lwp) THEN                    ! Namelist print
247         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme'
248         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens
249         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene
250         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT
251         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT
252         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een
253         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f
254         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix
255         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk
256      ENDIF
257
258      IF( ln_dynvor_msk )   CALL ctl_stop( 'dyn_vor_init:   masked vorticity is not currently not available')
259
260!!gm  this should be removed when choosing a unique strategy for fmask at the coast
261      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
262      ! at angles with three ocean points and one land point
263      IF(lwp) WRITE(numout,*)
264      IF(lwp) WRITE(numout,*) '      change fmask value in the angles (T)           ln_vorlat = ', ln_vorlat
265      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
266         DO jk = 1, jpk
267            DO jj = 1, jpjm1
268               DO ji = 1, jpim1
269                  IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
270                     & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
271               END DO
272            END DO
273         END DO
274         !
275         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
276         !
277      ENDIF
278!!gm end
279
280      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme)
281      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF
282      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF
283      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF
284      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF
285      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF
286      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF
287      !
288      IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
289      !                     
290      IF(lwp) WRITE(numout,*)        ! type of calculated vorticity (set ncor, nrvm, ntot)
291      ncor = np_COR                       ! planetary vorticity
292      SELECT CASE( n_dynadv )
293      CASE( np_LIN_dyn )
294         IF(lwp) WRITE(numout,*) '   ==>>>   linear dynamics : only Coriolis, no metric term'
295         nrvm = np_COR        ! planetary vorticity
296         ntot = np_COR        !     -         -
297      CASE( np_VEC_c2  )
298         CALL ctl_stop( 'dyncor_init : cor_ene requires FLUX form dynamics, not VECTOR form' )
299      CASE( np_FLX_c2 , np_FLX_ubs  )
300         IF(lwp) WRITE(numout,*) '   ==>>>   flux form dynamics : total vorticity = Coriolis + metric term'
301         nrvm = np_MET        ! metric term
302         ntot = np_CME        ! Coriolis + metric term
303         !
304         !                    ! pre-computed gradients for the metric term:
305         !                         !* T-point metric term :   pre-compute di(e2u)/2 and dj(e1v)/2
306            ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
307            DO jj = 2, jpjm1
308               DO ji = 2, jpim1
309                  di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj  ) ) * 0.5_wp
310                  dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji  ,jj-1) ) * 0.5_wp
311               END DO
312            END DO
313            CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. )   ! Lateral boundary conditions
314            !
315         !
316      END SELECT
317      !     
318      IF(lwp) WRITE(numout,*)
319      IF(lwp) WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)'
320      !
321   END SUBROUTINE dyn_cor_init
322
323   !!==============================================================================
324END MODULE dyncor
Note: See TracBrowser for help on using the repository browser.