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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 9 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS
  • Property svn:keywords set to Id
File size: 7.4 KB
RevLine 
[3]1MODULE dynkeg
2   !!======================================================================
3   !!                       ***  MODULE  dynkeg  ***
4   !! Ocean dynamics:  kinetic energy gradient trend
5   !!======================================================================
[5600]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
[503]10   !!----------------------------------------------------------------------
[3]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
[5034]17   USE trd_oce         ! trends: ocean variables
18   USE trddyn          ! trend manager: dynamics
19   !
[2715]20   USE in_out_manager  ! I/O manager
[5600]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
[3]31   
[5600]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   
[3]37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
[503]39   !!----------------------------------------------------------------------
[5600]40   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
[1152]41   !! $Id$
[2715]42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[503]43   !!----------------------------------------------------------------------
[3]44CONTAINS
45
[5600]46   SUBROUTINE dyn_keg( kt, kscheme )
[3]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      !!
[5600]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 ) ]
[5600]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      !!     
[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
[5034]68      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
[5600]69      !!
70      !! ** References : Arakawa, A., International Geophysics 2001.
71      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
[503]72      !!----------------------------------------------------------------------
[5600]73      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
74      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
[5034]75      !
[503]76      INTEGER  ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp) ::   zu, zv       ! temporary scalars
[3294]78      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
79      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
[3]80      !!----------------------------------------------------------------------
[3294]81      !
[5600]82      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
[3294]83      !
[5600]84      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
[3294]85      !
[3]86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
[5600]88         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
[3]89         IF(lwp) WRITE(numout,*) '~~~~~~~'
90      ENDIF
[216]91
[503]92      IF( l_trddyn ) THEN           ! Save ua and va trends
[5600]93         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[503]94         ztrdu(:,:,:) = ua(:,:,:) 
95         ztrdv(:,:,:) = va(:,:,:) 
[216]96      ENDIF
[3]97     
[5600]98      zhke(:,:,jpk) = 0._wp
99     
100      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
101      !
102      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
103         DO jk = 1, jpkm1
104            DO jj = 2, jpj
105               DO ji = fs_2, jpi   ! vector opt.
106                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
107                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
108                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
109                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
110                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
111               END DO 
112            END DO
113         END DO
114         !
115      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
116         DO jk = 1, jpkm1
117            DO jj = 2, jpjm1       
118               DO ji = fs_2, jpim1   ! vector opt.
119                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
120                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
121                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
122                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
123                     !
124                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
125                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
126                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
127                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
128                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
129               END DO 
130            END DO
131         END DO
132         CALL lbc_lnk( zhke, 'T', 1. )
133         !
134      END SELECT
135      !
136      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
137         DO jj = 2, jpjm1
[3]138            DO ji = fs_2, fs_jpim1   ! vector opt.
[503]139               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
140               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
[3]141            END DO
142         END DO
[5600]143      END DO
144      !
145      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
[503]146         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
147         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[5034]148         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
[5600]149         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
[216]150      ENDIF
[503]151      !
152      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
153         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
154      !
[5600]155      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
[2715]156      !
[5600]157      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
[3294]158      !
[3]159   END SUBROUTINE dyn_keg
160
161   !!======================================================================
162END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.