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

source: NEMO/trunk/src/SWE/dynkeg.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 4 years ago

Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision
File size: 8.3 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   
32   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme)
33!!an45
34   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2_wpo = 2   !: 2nd order centered scheme (wet point only : coastline at 45 degrees)
35   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW     = 1   !: Hollingsworth et al., QJRMS, 1983
36   !
37   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
38   
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id: dynkeg.F90 12377 2020-02-12 14:39:06Z acc $
44   !! Software governed by the CeCILL license (see ./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dyn_keg  ***
51      !!
52      !! ** Purpose :   Compute the now momentum trend due to the horizontal
53      !!      gradient of the horizontal kinetic energy and add it to the
54      !!      general momentum trend.
55      !!
56      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
57      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
58      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
59      !!              * kscheme = nkeg_HW : Hollingsworth correction following
60      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
61      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  )
62      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
63      !!     
64      !!      Take its horizontal gradient and add it to the general momentum
65      !!      trend.
66      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
67      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
68      !!
69      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
70      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
71      !!
72      !! ** References : Arakawa, A., International Geophysics 2001.
73      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
74      !!----------------------------------------------------------------------
75      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
76      INTEGER                             , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme
77      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices
78      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
79      !
80      INTEGER  ::   ji, jj, jk             ! dummy loop indices
81      REAL(wp) ::   zu, zv                   ! local scalars
82      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
83      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
84      !!----------------------------------------------------------------------
85      !
86      IF( ln_timing )   CALL timing_start('dyn_keg')
87      !
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
91         IF(lwp) WRITE(numout,*) '~~~~~~~'
92      ENDIF
93
94      IF( l_trddyn ) THEN           ! Save the input trends
95         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
96         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
97         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
98      ENDIF
99     
100      zhke(:,:,jpk) = 0._wp
101
102      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
103!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)
104      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --!
105         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
106            zu =  (   puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
107               &    + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )
108            zv =  (   pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
109               &    + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
110            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
111         END_3D
112!!an45         
113      !
114      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
115         DO_3D( 0, 1, 0, 1, 1, jpkm1 )
116            zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
117               &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)
118            zv =    pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
119               &  + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)
120            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
121         END_3D
122      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
123         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
124            zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    &
125               &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  &
126               &   +     ( 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) )   &
127               &   +     ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
128               !
129            zv = 8._wp * ( pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)    &
130               &         + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) )  &
131               &  +      ( 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) )   &
132               &  +      ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) )
133            zhke(ji,jj,jk) = r1_48 * ( zv + zu )
134         END_3D
135         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
136         !
137      END SELECT 
138      !
139      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
140         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
141         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
142      END_3D
143      !
144      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
145         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
146         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
147         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm )
148         DEALLOCATE( ztrdu , ztrdv )
149      ENDIF
150      !
151      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg  - Ua: ', mask1=umask,   &
152         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
153      !
154      IF( ln_timing )   CALL timing_stop('dyn_keg')
155      !
156   END SUBROUTINE dyn_keg
157
158   !!======================================================================
159END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.