source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90 @ 11048

Last change on this file since 11048 was 11048, checked in by girrmann, 19 months ago

dev_r10984_HPC-13 : Step 1, boundary is now detected all over the local domain, this does not change the result. Improve bdy treatment for bdy_rnf in bdytra.F90, this changes the result when keyword runoff is specified in namelist

  • Property svn:keywords set to Id
File size: 11.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   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  ::   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 1 to jpi
123                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
124                     IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE   ! to remove?
125                     DO jk = 1, jpkm1
126                        zhke(ji,jj,jk) = 0._wp
127                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
128                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
129                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
130                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
131                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
132                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
133                     END DO
134                  END DO
135               END IF
136               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies
137            END DO
138         END IF
139         !
140      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
141         DO jk = 1, jpkm1
142            DO jj = 2, jpjm1       
143               DO ji = fs_2, jpim1   ! vector opt.
144                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
145                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
146                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
147                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
148                     !
149                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
150                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
151                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
152                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
153                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
154               END DO 
155            END DO
156         END DO
157         IF (ln_bdy) THEN
158            ! Maria Luneva & Fred Wobus: July-2016
159            ! compensate for lack of turbulent kinetic energy on liquid bdy points
160            DO ib_bdy = 1, nb_bdy
161               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
162                  igrd = 1           ! compensation null velocity on land at the bdy
163                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
164                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi
165                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
166                     IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE   ! to remove
167                     DO jk = 1, jpkm1
168                        zhke(ji,jj,jk) = 0._wp
169                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
170                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
171                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
172                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
173                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
174                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
175                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
176                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
177                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
178                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
179                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
180                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
181                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
182                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
183                     END DO
184                  END DO
185               END IF
186            END DO
187         END IF
188         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
189         !
190      END SELECT 
191      !
192      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
193         DO jj = 2, jpjm1
194            DO ji = fs_2, fs_jpim1   ! vector opt.
195               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
196               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
197            END DO
198         END DO
199      END DO
200      !
201      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
202         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
203         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
204         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
205         DEALLOCATE( ztrdu , ztrdv )
206      ENDIF
207      !
208      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
209         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
210      !
211      IF( ln_timing )   CALL timing_stop('dyn_keg')
212      !
213   END SUBROUTINE dyn_keg
214
215   !!======================================================================
216END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.