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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 11107

Last change on this file since 11107 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 7.5 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
[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   !!----------------------------------------------------------------------
[5321]40   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
[6486]41   !! $Id$
[2715]42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[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      !
[503]76      INTEGER  ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp) ::   zu, zv       ! temporary scalars
[5328]78      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
79      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
[3]80      !!----------------------------------------------------------------------
[3294]81      !
[5321]82      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
[3294]83      !
[5328]84      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
[3294]85      !
[3]86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
[5321]88         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
[3]89         IF(lwp) WRITE(numout,*) '~~~~~~~'
[11101]90         IF(lflush) CALL flush(numout)
[3]91      ENDIF
[216]92
[503]93      IF( l_trddyn ) THEN           ! Save ua and va trends
[5321]94         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[5328]95         ztrdu(:,:,:) = ua(:,:,:) 
96         ztrdv(:,:,:) = va(:,:,:) 
[216]97      ENDIF
[5328]98     
99      zhke(:,:,jpk) = 0._wp
100     
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         !
[5328]116      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
117         DO jk = 1, jpkm1
118            DO jj = 2, jpjm1       
[5321]119               DO ji = fs_2, jpim1   ! vector opt.
120                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
121                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
122                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
123                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
124                     !
125                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
126                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
127                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
128                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
[5328]129                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
130               END DO 
[5321]131            END DO
[5328]132         END DO
133         CALL lbc_lnk( zhke, 'T', 1. )
[5321]134         !
[5328]135      END SELECT
136      !
137      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
[5321]138         DO jj = 2, jpjm1
[3]139            DO ji = fs_2, fs_jpim1   ! vector opt.
[5328]140               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
141               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
142            END DO
[3]143         END DO
[5321]144      END DO
145      !
146      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
[503]147         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
148         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]149         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
[5321]150         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[216]151      ENDIF
[503]152      !
153      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
154         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
155      !
[5328]156      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
[2715]157      !
[5321]158      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
[3294]159      !
[3]160   END SUBROUTINE dyn_keg
161
162   !!======================================================================
163END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.