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.
dynkeg.F90 in NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90 @ 13005

Last change on this file since 13005 was 13005, checked in by gm, 4 years ago

ADE and more options to AM98 config
File size: 14.2 KB
Line 
1MODULE dynkeg
2   !!======================================================================
3   !!                       ***  MODULE  dynkeg  ***
4   !! Ocean dynamics:  kinetic energy gradient trend
5   !!======================================================================
6   !! History :  1.0  !  1987-09  (P. Andrich, M.-A. Foujols)  Original code
7   !!            7.0  !  1997-05  (G. Madec)  Split dynber into dynkeg and dynhpg
8   !!  NEMO      1.0  !  2002-07  (G. Madec)  F90: Free form and module
9   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option
10   !!----------------------------------------------------------------------
11   
12   !!----------------------------------------------------------------------
13   !!   dyn_keg      : update the momentum trend with the horizontal tke
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE trd_oce         ! trends: ocean variables
18   USE trddyn          ! trend manager: dynamics
19   !
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE lib_mpp         ! MPP library
23   USE prtctl          ! Print control
24   USE timing          ! Timing
25   USE bdy_oce         ! ocean open boundary conditions
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dyn_keg    ! routine called by step module
31   PUBLIC   dyn_kegAD   ! routine called by step module
32   
33   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme)
34!!an45
35   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2_wpo = 2   !: 2nd order centered scheme (wet point only : coastline at 45 degrees)
36   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW     = 1   !: Hollingsworth et al., QJRMS, 1983
37   !
38   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
39   
40   !! * Substitutions
41#  include "do_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id: dynkeg.F90 12377 2020-02-12 14:39:06Z acc $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE dyn_keg  ***
52      !!
53      !! ** Purpose :   Compute the now momentum trend due to the horizontal
54      !!      gradient of the horizontal kinetic energy and add it to the
55      !!      general momentum trend.
56      !!
57      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
58      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
59      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
60      !!              * kscheme = nkeg_HW : Hollingsworth correction following
61      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
62      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  )
63      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
64      !!     
65      !!      Take its horizontal gradient and add it to the general momentum
66      !!      trend.
67      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
68      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
69      !!
70      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
71      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
72      !!
73      !! ** References : Arakawa, A., International Geophysics 2001.
74      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
75      !!----------------------------------------------------------------------
76      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
77      INTEGER                             , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme
78      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices
79      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
80      !
81      INTEGER  ::   ji, jj, jk             ! dummy loop indices
82      REAL(wp) ::   zu, zv                   ! local scalars
83      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
84      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
85      !!----------------------------------------------------------------------
86      !
87      IF( ln_timing )   CALL timing_start('dyn_keg')
88      !
89      IF( kt == nit000 ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
92         IF(lwp) WRITE(numout,*) '~~~~~~~'
93      ENDIF
94
95      IF( l_trddyn ) THEN           ! Save the input trends
96         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
97         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
98         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
99      ENDIF
100     
101      zhke(:,:,jpk) = 0._wp
102
103      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
104!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)
105      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --!
106         DO_3D_01_01( 1, jpkm1 )
107            zu =  (   puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
108               &    + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )
109            zv =  (   pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
110               &    + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
111            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
112         END_3D
113!!an45         
114      !
115      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
116         DO_3D_01_01( 1, jpkm1 )
117            zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
118               &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)
119            zv =    pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
120               &  + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)
121            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
122         END_3D
123      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
124         DO_3D_00_00( 1, jpkm1 )
125            zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    &
126               &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  &
127               &   +     ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) )   &
128               &   +     ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
129               !
130            zv = 8._wp * ( pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)    &
131               &         + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) )  &
132               &  +      ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) )   &
133               &  +      ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) )
134            zhke(ji,jj,jk) = r1_48 * ( zv + zu )
135         END_3D
136         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
137         !
138      END SELECT 
139      !
140      DO_3D_00_00( 1, jpkm1 )
141         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
142         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
143      END_3D
144      !
145      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
146         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
147         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
148         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm )
149         DEALLOCATE( ztrdu , ztrdv )
150      ENDIF
151      !
152      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg  - Ua: ', mask1=umask,   &
153         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
154      !
155      IF( ln_timing )   CALL timing_stop('dyn_keg')
156      !
157   END SUBROUTINE dyn_keg
158   
159   
160   SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs )
161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE dyn_kegAD  ***
163      !!
164      !! ** Purpose :   Compute the now momentum trend due to the horizontal
165      !!      gradient of the horizontal kinetic energy and add it to the
166      !!      general momentum trend.
167      !!
168      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
169      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
170      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
171      !!              * kscheme = nkeg_HW : Hollingsworth correction following
172      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
173      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  )
174      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
175      !!     
176      !!      Take its horizontal gradient and add it to the general momentum
177      !!      trend.
178      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
179      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
180      !!
181      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
182      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
183      !!
184      !! ** References : Arakawa, A., International Geophysics 2001.
185      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
186      !!----------------------------------------------------------------------
187      !
188      INTEGER                                  , INTENT( in )  ::  kt               ! ocean time-step index
189      INTEGER                                  , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme
190      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::  puu, pvv         ! ocean velocities at Kmm
191      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::  pu_rhs, pv_rhs   ! RHS
192      !
193      INTEGER  ::   ji, jj, jk             ! dummy loop indices
194      REAL(wp) ::   zu, zv                   ! local scalars
195      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
196      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
197      !!----------------------------------------------------------------------
198      !
199      IF( ln_timing )   CALL timing_start('dyn_keg')
200      !
201      IF( kt == nit000 ) THEN
202         IF(lwp) WRITE(numout,*)
203         IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme
204         IF(lwp) WRITE(numout,*) '~~~~~~~~~'
205      ENDIF
206     
207      zhke(:,:,jpk) = 0._wp
208
209      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
210!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)
211      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --!
212         DO_3D_01_01( 1, jpkm1 )
213            zu =  (   puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   &
214               &    + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )
215            zv =  (   pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   &
216               &    + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
217            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
218         END_3D
219!!an45         
220      !
221      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
222         DO_3D_01_01( 1, jpkm1 )
223            zu =    puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   &
224               &  + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)
225            zv =    pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   &
226               &  + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)
227            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
228         END_3D
229!!an 00_00 ?
230      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
231         DO_3D_00_00( 1, jpkm1 )
232            zu = 8._wp * ( puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)    &
233               &         + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) )  &
234               &   +     ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) )   &
235               &   +     ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) * ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) )
236               !
237            zv = 8._wp * ( pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)    &
238               &         + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) )  &
239               &  +      ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) )   &
240               &  +      ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) * ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) )
241            zhke(ji,jj,jk) = r1_48 * ( zv + zu )
242         END_3D
243         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
244         !
245      END SELECT 
246      !
247         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***!
248            !
249            DO_3D_00_00( 1, jpkm1 )
250               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
251               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
252            END_3D
253            !
254         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***!
255            !
256            DO_3D_00_00( 1, jpkm1 )
257               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
258            END_3D
259            !
260         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***!
261            !
262            DO_3D_00_00( 1, jpkm1 )
263               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
264            END_3D
265            !
266         ENDIF
267      IF( ln_timing )   CALL timing_stop('dyn_kegAD')
268      !
269   END SUBROUTINE dyn_kegAD
270   !!======================================================================
271END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.