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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 7508

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

changes on code duplication and workshare construct
  • Property svn:keywords set to Id
File size: 8.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
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/OPA 3.6 , NEMO Consortium (2015)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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   ! dummy loop indices
77      REAL(wp) ::   zu, zv       ! temporary scalars
78      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
79      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
80      !!----------------------------------------------------------------------
81      !
82      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
83      !
84      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
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 ua and va trends
93         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
94!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
95        DO jk = 1, jpk
96           DO jj = 1, jpj
97              DO ji = 1, jpi
98                 ztrdu(ji,jj,jk) = ua(ji,jj,jk)
99                 ztrdv(ji,jj,jk) = va(ji,jj,jk)
100              END DO
101           END DO
102        END DO
103      ENDIF
104!$OMP PARALLEL DO schedule(static) private(jj, ji)
105      DO jj = 1, jpj
106         DO ji = 1, jpi
107            zhke(ji,jj,jpk) = 0._wp
108         END DO
109      END DO
110     
111      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
112      !
113      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
114!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zu, zv)
115         DO jk = 1, jpkm1
116            DO jj = 2, jpj
117               DO ji = fs_2, jpi   ! vector opt.
118                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
119                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
120                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
121                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
122                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
123               END DO 
124            END DO
125         END DO
126         !
127      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
128!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zu, zv)
129         DO jk = 1, jpkm1
130            DO jj = 2, jpjm1       
131               DO ji = fs_2, jpim1   ! vector opt.
132                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
133                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
134                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
135                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
136                     !
137                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
138                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
139                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
140                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
141                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
142               END DO 
143            END DO
144         END DO
145         CALL lbc_lnk( zhke, 'T', 1. )
146         !
147      END SELECT
148      !
149!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
150      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
151         DO jj = 2, jpjm1
152            DO ji = fs_2, fs_jpim1   ! vector opt.
153               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
154               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
155            END DO
156         END DO
157      END DO
158      !
159      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
160!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
161           DO jk = 1, jpk
162              DO jj = 1, jpj
163                 DO ji = 1, jpi
164                    ztrdu(ji,jj,jk) = ua(ji,jj,jk) - ztrdu(ji,jj,jk)
165                    ztrdv(ji,jj,jk) = va(ji,jj,jk) - ztrdv(ji,jj,jk)
166                 END DO
167              END DO
168           END DO
169         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
170         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
171      ENDIF
172      !
173      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
174         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
175      !
176      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
177      !
178      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
179      !
180   END SUBROUTINE dyn_keg
181
182   !!======================================================================
183END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.