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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 10.0 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 wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
[7646]26   USE bdy_oce         ! ocean open boundary conditions
[3]27
28   IMPLICIT NONE
29   PRIVATE
30
[503]31   PUBLIC   dyn_keg    ! routine called by step module
[5328]32   
[5321]33   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme)
34   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983
35   !
36   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
[5328]37   
[3]38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
[503]40   !!----------------------------------------------------------------------
[5321]41   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
[5328]42   !! $Id$
[2715]43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[503]44   !!----------------------------------------------------------------------
[3]45CONTAINS
46
[5321]47   SUBROUTINE dyn_keg( kt, kscheme )
[3]48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_keg  ***
50      !!
51      !! ** Purpose :   Compute the now momentum trend due to the horizontal
[5328]52      !!      gradient of the horizontal kinetic energy and add it to the
[3]53      !!      general momentum trend.
54      !!
[5328]55      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
56      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
[3]57      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
[5321]58      !!              * kscheme = nkeg_HW : Hollingsworth correction following
59      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
60      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
61      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
[5328]62      !!     
[3]63      !!      Take its horizontal gradient and add it to the general momentum
64      !!      trend (ua,va).
65      !!         ua = ua - 1/e1u di[ zhke ]
66      !!         va = va - 1/e2v dj[ zhke ]
67      !!
68      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
[4990]69      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
[5321]70      !!
71      !! ** References : Arakawa, A., International Geophysics 2001.
72      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
[503]73      !!----------------------------------------------------------------------
[5321]74      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
[5328]75      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
[4990]76      !
[503]77      INTEGER  ::   ji, jj, jk   ! dummy loop indices
78      REAL(wp) ::   zu, zv       ! temporary scalars
[5328]79      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
80      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
[7646]81      INTEGER  ::   jb                 ! dummy loop indices
82      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
83      INTEGER  ::   fu, fv
[3]84      !!----------------------------------------------------------------------
[3294]85      !
[5321]86      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
[3294]87      !
[5328]88      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
[3294]89      !
[3]90      IF( kt == nit000 ) THEN
91         IF(lwp) WRITE(numout,*)
[5321]92         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
[3]93         IF(lwp) WRITE(numout,*) '~~~~~~~'
94      ENDIF
[216]95
[503]96      IF( l_trddyn ) THEN           ! Save ua and va trends
[5321]97         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[7698]98!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
99         DO jk = 1, jpk
100            DO jj = 1, jpj
101               DO ji = 1, jpi
102                  ztrdu(ji,jj,jk) = ua(ji,jj,jk)
103                  ztrdv(ji,jj,jk) = va(ji,jj,jk)
104               END DO
105            END DO
106         END DO
[216]107      ENDIF
[7698]108!$OMP PARALLEL DO schedule(static) private(jj, ji)
109      DO jj = 1, jpj
110         DO ji = 1, jpi
111            zhke(ji,jj,jpk) = 0._wp
112         END DO
113      END DO
[5328]114     
[7646]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 = 2           ! Copying normal velocity into points outside bdy
121               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
122                  DO jk = 1, jpkm1
123                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
124                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
125                     fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
126                     un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
127                  END DO
128               END DO
129               !
130               igrd = 3           ! Copying normal velocity into points outside bdy
131               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
132                  DO jk = 1, jpkm1
133                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
134                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
135                     fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
136                     vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
137                  END DO
138               END DO
139            ENDIF
140         ENDDO 
141      ENDIF
142
[5328]143      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
144      !
145      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
[7698]146!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zu, zv)
[5328]147         DO jk = 1, jpkm1
[5321]148            DO jj = 2, jpj
149               DO ji = fs_2, jpi   ! vector opt.
150                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
151                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
152                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
153                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
[5328]154                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
155               END DO 
[5321]156            END DO
[5328]157         END DO
[5321]158         !
[5328]159      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
[7698]160!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zu, zv)
[5328]161         DO jk = 1, jpkm1
162            DO jj = 2, jpjm1       
[5321]163               DO ji = fs_2, jpim1   ! vector opt.
164                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
165                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
166                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
167                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
168                     !
169                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
170                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
171                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
172                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
[5328]173                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
174               END DO 
[5321]175            END DO
[5328]176         END DO
177         CALL lbc_lnk( zhke, 'T', 1. )
[5321]178         !
[5328]179      END SELECT
[7646]180
181      IF (ln_bdy) THEN
182         ! restore velocity masks at points outside boundary
[7698]183!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
184         DO jk = 1, jpk
185            DO jj = 1, jpj
186               DO ji = 1, jpi
187                  un(ji,jj,jk) = un(ji,jj,jk) * umask(ji,jj,jk)
188                  vn(ji,jj,jk) = vn(ji,jj,jk) * vmask(ji,jj,jk)
189               END DO
190            END DO
191         END DO
192      ENDIF
[7646]193
[5328]194      !
[7698]195!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
[5328]196      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
[5321]197         DO jj = 2, jpjm1
[3]198            DO ji = fs_2, fs_jpim1   ! vector opt.
[5328]199               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
200               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
201            END DO
[3]202         END DO
[5321]203      END DO
204      !
205      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
[7698]206!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
207           DO jk = 1, jpk
208              DO jj = 1, jpj
209                 DO ji = 1, jpi
210                    ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
211                    ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
212                 END DO
213              END DO
214           END DO
[4990]215         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
[5321]216         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[216]217      ENDIF
[503]218      !
219      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
220         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
221      !
[5328]222      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
[2715]223      !
[5321]224      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
[3294]225      !
[3]226   END SUBROUTINE dyn_keg
227
228   !!======================================================================
229END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.