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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

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

Last change on this file since 10789 was 10789, checked in by davestorkey, 5 years ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: Convert first batch of DYN routines and "wn" -> "ww".
  • Property svn:keywords set to Id
File size: 9.4 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
[10789]46   SUBROUTINE dyn_keg( kt, ktlev, kscheme, pu_rhs, pv_rhs )
[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
[10789]56      !!         zhke = 1/2 [ mi-1( uu(:,:,:,ktlev)^2 ) + mj-1( vv(:,:,:,ktlev)^2 ) ]
[5321]57      !!              * kscheme = nkeg_HW : Hollingsworth correction following
58      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
[10789]59      !!         zhke = 1/6 [ mi-1(  2 * uu(:,:,:,ktlev)^2 + ((uu(j+1,ktlev)+uu(j-1,ktlev))/2)^2  )
60      !!                    + mj-1(  2 * vv(:,:,:,ktlev)^2 + ((vv(i+1,ktlev)+vv(i-1,ktlev))/2)^2  ) ]
[5328]61      !!     
[3]62      !!      Take its horizontal gradient and add it to the general momentum
[10789]63      !!      trend (pu_rhs,pv_rhs).
64      !!         pu_rhs = pu_rhs - 1/e1u di[ zhke ]
65      !!         pv_rhs = pv_rhs - 1/e2v dj[ zhke ]
[3]66      !!
[10789]67      !! ** Action : - Update the (pu_rhs, pv_rhs) 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
[10789]74      INTEGER, INTENT( in ) ::   ktlev     ! time level index for source terms
[5328]75      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
[10789]76      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
[4990]77      !
[9019]78      INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices
79      INTEGER  ::   ii, ifu, ib_bdy   ! local integers
80      INTEGER  ::   ij, ifv, igrd     !   -       -
81      REAL(wp) ::   zu, zv            ! local scalars
82      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
83      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
[3]84      !!----------------------------------------------------------------------
[3294]85      !
[9019]86      IF( ln_timing )   CALL timing_start('dyn_keg')
[3294]87      !
[3]88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
[5321]90         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
[3]91         IF(lwp) WRITE(numout,*) '~~~~~~~'
92      ENDIF
[216]93
[9019]94      IF( l_trddyn ) THEN           ! Save the input trends
95         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
[10789]96         ztrdu(:,:,:) = pu_rhs(:,:,:) 
97         ztrdv(:,:,:) = pv_rhs(:,:,:) 
[216]98      ENDIF
[5328]99     
[7753]100      zhke(:,:,jpk) = 0._wp
101     
[7646]102      IF (ln_bdy) THEN
103         ! Maria Luneva & Fred Wobus: July-2016
104         ! compensate for lack of turbulent kinetic energy on liquid bdy points
105         DO ib_bdy = 1, nb_bdy
106            IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
107               igrd = 2           ! Copying normal velocity into points outside bdy
108               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
109                  DO jk = 1, jpkm1
110                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
111                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
[9019]112                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
[10789]113                     uu(ii-ifu,ij,jk,ktlev) = uu(ii,ij,jk,ktlev) * umask(ii,ij,jk)
[7646]114                  END DO
115               END DO
116               !
117               igrd = 3           ! Copying normal velocity into points outside bdy
118               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
119                  DO jk = 1, jpkm1
120                     ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
121                     ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
[9019]122                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
[10789]123                     vv(ii,ij-ifv,jk,ktlev) = vv(ii,ij,jk,ktlev) * vmask(ii,ij,jk)
[7646]124                  END DO
125               END DO
126            ENDIF
127         ENDDO 
128      ENDIF
129
[5328]130      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
131      !
132      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
133         DO jk = 1, jpkm1
[5321]134            DO jj = 2, jpj
135               DO ji = fs_2, jpi   ! vector opt.
[10789]136                  zu =    uu(ji-1,jj  ,jk,ktlev) * uu(ji-1,jj  ,jk,ktlev)   &
137                     &  + uu(ji  ,jj  ,jk,ktlev) * uu(ji  ,jj  ,jk,ktlev)
138                  zv =    vv(ji  ,jj-1,jk,ktlev) * vv(ji  ,jj-1,jk,ktlev)   &
139                     &  + vv(ji  ,jj  ,jk,ktlev) * vv(ji  ,jj  ,jk,ktlev)
[5328]140                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
141               END DO 
[5321]142            END DO
[5328]143         END DO
[5321]144         !
[5328]145      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
146         DO jk = 1, jpkm1
147            DO jj = 2, jpjm1       
[5321]148               DO ji = fs_2, jpim1   ! vector opt.
[10789]149                  zu = 8._wp * ( uu(ji-1,jj  ,jk,ktlev) * uu(ji-1,jj  ,jk,ktlev)    &
150                     &         + uu(ji  ,jj  ,jk,ktlev) * uu(ji  ,jj  ,jk,ktlev) )  &
151                     &   +     ( uu(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) ) * ( uu(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) )   &
152                     &   +     ( uu(ji  ,jj-1,jk,ktlev) + uu(ji  ,jj+1,jk,ktlev) ) * ( uu(ji  ,jj-1,jk,ktlev) + uu(ji  ,jj+1,jk,ktlev) )
[5321]153                     !
[10789]154                  zv = 8._wp * ( vv(ji  ,jj-1,jk,ktlev) * vv(ji  ,jj-1,jk,ktlev)    &
155                     &         + vv(ji  ,jj  ,jk,ktlev) * vv(ji  ,jj  ,jk,ktlev) )  &
156                     &  +      ( vv(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) ) * ( vv(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) )   &
157                     &  +      ( vv(ji-1,jj  ,jk,ktlev) + vv(ji+1,jj  ,jk,ktlev) ) * ( vv(ji-1,jj  ,jk,ktlev) + vv(ji+1,jj  ,jk,ktlev) )
[5328]158                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
159               END DO 
[5321]160            END DO
[5328]161         END DO
[10425]162         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
[5321]163         !
[5328]164      END SELECT
[7646]165
166      IF (ln_bdy) THEN
167         ! restore velocity masks at points outside boundary
[10789]168         uu(:,:,:,ktlev) = uu(:,:,:,ktlev) * umask(:,:,:)
169         vv(:,:,:,ktlev) = vv(:,:,:,ktlev) * vmask(:,:,:)
[7753]170      ENDIF     
[7646]171
[5328]172      !
173      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
[5321]174         DO jj = 2, jpjm1
[3]175            DO ji = fs_2, fs_jpim1   ! vector opt.
[10789]176               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
177               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
[5328]178            END DO
[3]179         END DO
[5321]180      END DO
181      !
182      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
[10789]183         ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)
184         ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)
[4990]185         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
[9019]186         DEALLOCATE( ztrdu , ztrdv )
[216]187      ENDIF
[503]188      !
189      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
190         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
191      !
[9019]192      IF( ln_timing )   CALL timing_stop('dyn_keg')
[2715]193      !
[3]194   END SUBROUTINE dyn_keg
195
196   !!======================================================================
197END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.