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/UKMO/NEMO_4.0_mirror/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/OCE/DYN/dynkeg.F90

Last change on this file was 11081, checked in by davestorkey, 5 years ago

UKMO/NEMO_4.0_mirror : update to be a copy of rev 11079 of release-4.0.
File size: 11.1 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   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983
34   !
35   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
36   
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dyn_keg( kt, kscheme )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dyn_keg  ***
49      !!
50      !! ** Purpose :   Compute the now momentum trend due to the horizontal
51      !!      gradient of the horizontal kinetic energy and add it to the
52      !!      general momentum trend.
53      !!
54      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
55      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
56      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
57      !!              * kscheme = nkeg_HW : Hollingsworth correction following
58      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
59      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
60      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
61      !!     
62      !!      Take its horizontal gradient and add it to the general momentum
63      !!      trend (ua,va).
64      !!         ua = ua - 1/e1u di[ zhke ]
65      !!         va = va - 1/e2v dj[ zhke ]
66      !!
67      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
68      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
69      !!
70      !! ** References : Arakawa, A., International Geophysics 2001.
71      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
72      !!----------------------------------------------------------------------
73      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
74      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
75      !
76      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices
77      INTEGER  ::   ifu, ifv, igrd, ib_bdy   ! local integers
78      REAL(wp) ::   zu, zv                   ! local scalars
79      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
80      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
81      REAL(wp)  :: zweightu, zweightv
82      !!----------------------------------------------------------------------
83      !
84      IF( ln_timing )   CALL timing_start('dyn_keg')
85      !
86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
89         IF(lwp) WRITE(numout,*) '~~~~~~~'
90      ENDIF
91
92      IF( l_trddyn ) THEN           ! Save the input trends
93         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
94         ztrdu(:,:,:) = ua(:,:,:) 
95         ztrdv(:,:,:) = va(:,:,:) 
96      ENDIF
97     
98      zhke(:,:,jpk) = 0._wp
99
100      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
101      !
102      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
103         DO jk = 1, jpkm1
104            DO jj = 2, jpj
105               DO ji = fs_2, jpi   ! vector opt.
106                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
107                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
108                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
109                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
110                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
111               END DO 
112            END DO
113         END DO
114         !
115         IF (ln_bdy) THEN
116            ! Maria Luneva & Fred Wobus: July-2016
117            ! compensate for lack of turbulent kinetic energy on liquid bdy points
118            DO ib_bdy = 1, nb_bdy
119               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
120                  igrd = 1           ! compensating null velocity on the bdy
121                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
122                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
123                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
124                     DO jk = 1, jpkm1
125                        zhke(ji,jj,jk) = 0._wp
126                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
127                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
128                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
129                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
130                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
131                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
132                     END DO
133                  END DO
134               END IF
135               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies
136            END DO
137         END IF
138         !
139      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
140         DO jk = 1, jpkm1
141            DO jj = 2, jpjm1       
142               DO ji = fs_2, jpim1   ! vector opt.
143                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
144                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
145                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
146                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
147                     !
148                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
149                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
150                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
151                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
152                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
153               END DO 
154            END DO
155         END DO
156         IF (ln_bdy) THEN
157            ! Maria Luneva & Fred Wobus: July-2016
158            ! compensate for lack of turbulent kinetic energy on liquid bdy points
159            DO ib_bdy = 1, nb_bdy
160               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
161                  igrd = 1           ! compensation null velocity on land at the bdy
162                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
163                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
164                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
165                     DO jk = 1, jpkm1
166                        zhke(ji,jj,jk) = 0._wp
167                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
168                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
169                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
170                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
171                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
172                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
173                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
174                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
175                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
176                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
177                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
178                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
179                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
180                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
181                     END DO
182                  END DO
183               END IF
184            END DO
185         END IF
186         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
187         !
188      END SELECT 
189      !
190      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
191         DO jj = 2, jpjm1
192            DO ji = fs_2, fs_jpim1   ! vector opt.
193               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
194               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
195            END DO
196         END DO
197      END DO
198      !
199      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
200         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
201         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
202         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
203         DEALLOCATE( ztrdu , ztrdv )
204      ENDIF
205      !
206      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
207         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
208      !
209      IF( ln_timing )   CALL timing_stop('dyn_keg')
210      !
211   END SUBROUTINE dyn_keg
212
213   !!======================================================================
214END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.