source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90 @ 10874

Last change on this file since 10874 was 10874, checked in by davestorkey, 20 months ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Revert all changes so far in preparation for implementation of new design.

  • Property svn:keywords set to Id
File size: 8.8 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  ::   ii, ifu, ib_bdy   ! local integers
78      INTEGER  ::   ij, ifv, igrd     !   -       -
79      REAL(wp) ::   zu, zv            ! local scalars
80      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
81      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
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      IF (ln_bdy) THEN
101         ! Maria Luneva & Fred Wobus: July-2016
102         ! compensate for lack of turbulent kinetic energy on liquid bdy points
103         DO ib_bdy = 1, nb_bdy
104            IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
105               igrd = 2           ! Copying normal velocity into points outside bdy
106               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
107                  DO jk = 1, jpkm1
108                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
109                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
110                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
111                     un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
112                  END DO
113               END DO
114               !
115               igrd = 3           ! Copying normal velocity into points outside bdy
116               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
117                  DO jk = 1, jpkm1
118                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
119                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
120                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
121                     vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
122                  END DO
123               END DO
124            ENDIF
125         ENDDO 
126      ENDIF
127
128      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
129      !
130      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
131         DO jk = 1, jpkm1
132            DO jj = 2, jpj
133               DO ji = fs_2, jpi   ! vector opt.
134                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
135                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
136                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
137                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
138                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
139               END DO 
140            END DO
141         END DO
142         !
143      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
144         DO jk = 1, jpkm1
145            DO jj = 2, jpjm1       
146               DO ji = fs_2, jpim1   ! vector opt.
147                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
148                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
149                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
150                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
151                     !
152                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
153                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
154                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
155                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
156                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
157               END DO 
158            END DO
159         END DO
160         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
161         !
162      END SELECT
163
164      IF (ln_bdy) THEN
165         ! restore velocity masks at points outside boundary
166         un(:,:,:) = un(:,:,:) * umask(:,:,:)
167         vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
168      ENDIF     
169
170      !
171      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
172         DO jj = 2, jpjm1
173            DO ji = fs_2, fs_jpim1   ! vector opt.
174               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
175               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
176            END DO
177         END DO
178      END DO
179      !
180      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
181         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
182         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
183         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
184         DEALLOCATE( ztrdu , ztrdv )
185      ENDIF
186      !
187      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
188         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
189      !
190      IF( ln_timing )   CALL timing_stop('dyn_keg')
191      !
192   END SUBROUTINE dyn_keg
193
194   !!======================================================================
195END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.