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/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN – NEMO

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

Last change on this file since 11195 was 11195, checked in by girrmann, 5 years ago

dev_r10984_HPC-13 : update of trc_bdy following [11191], merge of lbc_lnk and lbc_bdy_lnk in a single lbc_lnk routine, see #2285
  • Property svn:keywords set to Id
File size: 12.0 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      LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how bdy communications are to be carried out
83      !!----------------------------------------------------------------------
84      !
85      IF( ln_timing )   CALL timing_start('dyn_keg')
86      !
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
90         IF(lwp) WRITE(numout,*) '~~~~~~~'
91      ENDIF
92
93      IF( l_trddyn ) THEN           ! Save the input trends
94         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
95         ztrdu(:,:,:) = ua(:,:,:) 
96         ztrdv(:,:,:) = va(:,:,:) 
97      ENDIF
98     
99      zhke(:,:,jpk) = 0._wp
100
101      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
102      !
103      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
104         DO jk = 1, jpkm1
105            DO jj = 2, jpj
106               DO ji = fs_2, jpi   ! vector opt.
107                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
108                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
109                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
110                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
111                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
112               END DO 
113            END DO
114         END DO
115         !
116         IF (ln_bdy) THEN
117            ! Maria Luneva & Fred Wobus: July-2016
118            ! compensate for lack of turbulent kinetic energy on liquid bdy points
119            DO ib_bdy = 1, nb_bdy
120               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
121                  igrd = 1           ! compensating null velocity on the bdy
122                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
123                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi
124                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
125                     IF( ji == 1 .OR. jj == 1 )  CYCLE
126                     DO jk = 1, jpkm1
127                        zhke(ji,jj,jk) = 0._wp
128                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
129                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
130                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
131                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
132                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
133                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
134                     END DO
135                  END DO
136               END IF
137            END DO
138            ! send jpi-1, jpj-1 and receive 1 used in the computation of the speed tendencies
139            llsend1(:) = .false.
140            llrecv1(:) = .false.
141            DO ib_bdy = 1, nb_bdy
142               llsend1(2) = llsend1(2) .OR. ANY(lsend_bdy(ib_bdy,igrd,2,:))   ! send east
143               llsend1(4) = llsend1(4) .OR. ANY(lsend_bdy(ib_bdy,igrd,4,:))   ! send north
144               llrecv1(1) = llrecv1(1) .OR. ANY(lrecv_bdy(ib_bdy,igrd,1,:))   ! receive west
145               llrecv1(3) = llrecv1(3) .OR. ANY(lrecv_bdy(ib_bdy,igrd,3,:))   ! receive south
146            END DO
147   
148            IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
149               CALL lbc_lnk( 'bdydyn2d', zhke, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
150            END IF
151         END IF
152         !
153      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
154         DO jk = 1, jpkm1
155            DO jj = 2, jpjm1       
156               DO ji = fs_2, jpim1   ! vector opt.
157                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
158                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
159                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
160                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
161                     !
162                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
163                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
164                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
165                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
166                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
167               END DO 
168            END DO
169         END DO
170         IF (ln_bdy) THEN
171            ! Maria Luneva & Fred Wobus: July-2016
172            ! compensate for lack of turbulent kinetic energy on liquid bdy points
173            DO ib_bdy = 1, nb_bdy
174               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
175                  igrd = 1           ! compensation null velocity on land at the bdy
176                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
177                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi
178                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
179                     IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE
180                     DO jk = 1, jpkm1
181                        zhke(ji,jj,jk) = 0._wp
182                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
183                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
184                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
185                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
186                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
187                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
188                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
189                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
190                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
191                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
192                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
193                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
194                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
195                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
196                     END DO
197                  END DO
198               END IF
199            END DO
200         END IF
201         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
202         !
203      END SELECT 
204      !
205      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
206         DO jj = 2, jpjm1
207            DO ji = fs_2, fs_jpim1   ! vector opt.
208               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
209               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
210            END DO
211         END DO
212      END DO
213      !
214      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
215         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
216         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
217         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
218         DEALLOCATE( ztrdu , ztrdv )
219      ENDIF
220      !
221      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
222         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
223      !
224      IF( ln_timing )   CALL timing_stop('dyn_keg')
225      !
226   END SUBROUTINE dyn_keg
227
228   !!======================================================================
229END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.