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 @ 11071

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

dev_r10984_HPC-13 : step 2, remove unneeded communications, see #2285
  • Property svn:keywords set to Id
File size: 11.9 KB
RevLine 
[3]1MODULE dynkeg
2   !!======================================================================
3   !!                       ***  MODULE  dynkeg  ***
4   !! Ocean dynamics:  kinetic energy gradient trend
5   !!======================================================================
[5321]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
[5328]9   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option
[503]10   !!----------------------------------------------------------------------
[5328]11   
[3]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
[4990]17   USE trd_oce         ! trends: ocean variables
18   USE trddyn          ! trend manager: dynamics
19   !
[2715]20   USE in_out_manager  ! I/O manager
[5321]21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[2715]22   USE lib_mpp         ! MPP library
[258]23   USE prtctl          ! Print control
[3294]24   USE timing          ! Timing
[7646]25   USE bdy_oce         ! ocean open boundary conditions
[3]26
27   IMPLICIT NONE
28   PRIVATE
29
[503]30   PUBLIC   dyn_keg    ! routine called by step module
[5328]31   
[5321]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)
[5328]36   
[3]37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
[503]39   !!----------------------------------------------------------------------
[9598]40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5328]41   !! $Id$
[10068]42   !! Software governed by the CeCILL license (see ./LICENSE)
[503]43   !!----------------------------------------------------------------------
[3]44CONTAINS
45
[5321]46   SUBROUTINE dyn_keg( kt, kscheme )
[3]47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dyn_keg  ***
49      !!
50      !! ** Purpose :   Compute the now momentum trend due to the horizontal
[5328]51      !!      gradient of the horizontal kinetic energy and add it to the
[3]52      !!      general momentum trend.
53      !!
[5328]54      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
55      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
[3]56      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
[5321]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  ) ]
[5328]61      !!     
[3]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
[4990]68      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
[5321]69      !!
70      !! ** References : Arakawa, A., International Geophysics 2001.
71      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
[503]72      !!----------------------------------------------------------------------
[5321]73      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
[5328]74      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
[4990]75      !
[10996]76      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices
[11024]77      INTEGER  ::   igrd, ib_bdy             ! local integers
[10996]78      REAL(wp) ::   zu, zv                   ! local scalars
[9019]79      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
80      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
[10996]81      REAL(wp)  :: zweightu, zweightv
[11071]82      LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how bdy communications are to be carried out
[3]83      !!----------------------------------------------------------------------
[3294]84      !
[9019]85      IF( ln_timing )   CALL timing_start('dyn_keg')
[3294]86      !
[3]87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
[5321]89         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
[3]90         IF(lwp) WRITE(numout,*) '~~~~~~~'
91      ENDIF
[216]92
[9019]93      IF( l_trddyn ) THEN           ! Save the input trends
94         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
[7753]95         ztrdu(:,:,:) = ua(:,:,:) 
96         ztrdv(:,:,:) = va(:,:,:) 
[216]97      ENDIF
[5328]98     
[7753]99      zhke(:,:,jpk) = 0._wp
[7646]100
[5328]101      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
102      !
103      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
104         DO jk = 1, jpkm1
[5321]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)
[5328]111                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
112               END DO 
[5321]113            END DO
[5328]114         END DO
[5321]115         !
[10996]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)
[11048]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
[11049]125                     IF( ji == 1 .OR. jj == 1 )  CYCLE
[10996]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
[11071]138            ! send jpi-1, jpj-1 and receive 1 used in the computation of the speed tendencies
139            llsend1(:) = .false.
140            llrecv1(:) = .false.
[11067]141            DO ib_bdy = 1, nb_bdy
[11071]142               llsend1(2) = llsend1(2) .OR. lsend_bdy(ib_bdy,igrd,2)   ! send east
143               llsend1(4) = llsend1(4) .OR. lsend_bdy(ib_bdy,igrd,4)   ! send north
144               llrecv1(1) = llrecv1(1) .OR. lrecv_bdy(ib_bdy,igrd,1)   ! receive west
145               llrecv1(3) = llrecv1(3) .OR. lrecv_bdy(ib_bdy,igrd,3)   ! receive south
[11067]146            END DO
[11071]147   
148            IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
149               CALL lbc_bdy_lnk( 'bdydyn2d', llsend1, llrecv1, zhke, 'T',  1. )
[11067]150            END IF
[10996]151         END IF
152         !
[5328]153      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
154         DO jk = 1, jpkm1
155            DO jj = 2, jpjm1       
[5321]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) )
[5328]166                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
167               END DO 
[5321]168            END DO
[5328]169         END DO
[10996]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)
[11048]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
[11049]179                     IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE
[10996]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
[10425]201         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
[5321]202         !
[10996]203      END SELECT 
[5328]204      !
205      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
[5321]206         DO jj = 2, jpjm1
[3]207            DO ji = fs_2, fs_jpim1   ! vector opt.
[5328]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
[3]211         END DO
[5321]212      END DO
213      !
214      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
[7753]215         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
216         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]217         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
[9019]218         DEALLOCATE( ztrdu , ztrdv )
[216]219      ENDIF
[503]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      !
[9019]224      IF( ln_timing )   CALL timing_stop('dyn_keg')
[2715]225      !
[3]226   END SUBROUTINE dyn_keg
227
228   !!======================================================================
229END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.