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 branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 11186

Last change on this file since 11186 was 11186, checked in by kingr, 5 years ago

Merging changes from AMM15_v3_6_STABLE_package_collate r10728:11063
File size: 9.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 wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
26#if defined key_bdy 
27   USE bdy_oce         ! ocean open boundary conditions
28#endif
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   dyn_keg    ! routine called by step module
34   
35   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme)
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 "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE dyn_keg( kt, kscheme )
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 + ((un(j+1)+un(j-1))/2)^2  )
63      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
64      !!     
65      !!      Take its horizontal gradient and add it to the general momentum
66      !!      trend (ua,va).
67      !!         ua = ua - 1/e1u di[ zhke ]
68      !!         va = va - 1/e2v dj[ zhke ]
69      !!
70      !! ** Action : - Update the (ua, va) 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   type of KEG scheme
78      !
79      INTEGER  ::   ji, jj, jk   ! dummy loop indices
80      REAL(wp) ::   zu, zv       ! temporary scalars
81      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
82      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
83#if defined key_bdy 
84      INTEGER  ::   jb                 ! dummy loop indices
85      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
86      INTEGER  ::   fu, fv 
87#endif 
88      !!----------------------------------------------------------------------
89      !
90      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
91      !
92      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
93      !
94      IF( kt == nit000 ) THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
97         IF(lwp) WRITE(numout,*) '~~~~~~~'
98      ENDIF
99
100      IF( l_trddyn ) THEN           ! Save ua and va trends
101         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
102         ztrdu(:,:,:) = ua(:,:,:) 
103         ztrdv(:,:,:) = va(:,:,:) 
104      ENDIF
105     
106      zhke(:,:,jpk) = 0._wp
107
108#if defined key_bdy 
109      ! Maria Luneva & Fred Wobus: July-2016
110      ! compensate for lack of turbulent kinetic energy on liquid bdy points
111      DO ib_bdy = 1, nb_bdy 
112         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
113            igrd = 2           ! Copying normal velocity into points outside bdy
114            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
115               DO jk = 1, jpkm1 
116                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
117                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
118                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
119                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
120               END DO
121            END DO 
122            !
123            igrd = 3           ! Copying normal velocity into points outside bdy
124            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
125               DO jk = 1, jpkm1 
126                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
127                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
128                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
129                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
130               END DO
131            END DO
132         ENDIF
133      ENDDO 
134!CEOD ADD
135      CALL lbc_lnk( un,'U',-1. )
136      CALL lbc_lnk( vn,'V',-1. )
137!CEOD ADD
138#endif       
139
140      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
141      !
142      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
143         DO jk = 1, jpkm1
144            DO jj = 2, jpj
145               DO ji = fs_2, jpi   ! vector opt.
146                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
147                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
148                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
149                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
150                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
151               END DO 
152            END DO
153         END DO
154         !
155         CALL lbc_lnk( zhke, 'T', 1. )
156         !
157      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
158         DO jk = 1, jpkm1
159            DO jj = 2, jpjm1       
160               DO ji = fs_2, jpim1   ! vector opt.
161                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
162                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
163                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
164                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
165                     !
166                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
167                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
168                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
169                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
170                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
171               END DO 
172            END DO
173         END DO
174         CALL lbc_lnk( zhke, 'T', 1. )
175         !
176      END SELECT
177      !
178
179#if defined key_bdy 
180      ! restore velocity masks at points outside boundary
181      un(:,:,:) = un(:,:,:) * umask(:,:,:) 
182      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
183#endif 
184
185      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
186         DO jj = 2, jpjm1
187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
189               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
190            END DO
191         END DO
192      END DO
193      !
194      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
195         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
196         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
197         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
198         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
199      ENDIF
200      !
201      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
202         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
203      !
204      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
205      !
206      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
207      !
208   END SUBROUTINE dyn_keg
209
210   !!======================================================================
211END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.