source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 8300

Last change on this file since 8300 was 8300, checked in by glong, 4 years ago

Changed diagnostics to calculate the contributions by taking after-before before going into the specific numerical scheme for the model.

File size: 7.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   USE divcur          ! for dyn_vrt_dia_3d
20   !
21   USE in_out_manager  ! I/O manager
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp         ! MPP library
24   USE prtctl          ! Print control
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   dyn_keg    ! routine called by step module
32   
33   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme)
34   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983
35   !
36   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
37   
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dyn_keg( kt, kscheme )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_keg  ***
50      !!
51      !! ** Purpose :   Compute the now momentum trend due to the horizontal
52      !!      gradient of the horizontal kinetic energy and add it to the
53      !!      general momentum trend.
54      !!
55      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
56      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
57      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
58      !!              * kscheme = nkeg_HW : Hollingsworth correction following
59      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
60      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
61      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
62      !!     
63      !!      Take its horizontal gradient and add it to the general momentum
64      !!      trend (ua,va).
65      !!         ua = ua - 1/e1u di[ zhke ]
66      !!         va = va - 1/e2v dj[ zhke ]
67      !!
68      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
69      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
70      !!
71      !! ** References : Arakawa, A., International Geophysics 2001.
72      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
73      !!----------------------------------------------------------------------
74      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
75      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
76      !
77      CHARACTER(len=3) :: id_vrt_dia_keg = "keg" ! TODO remove once flags set properly
78      INTEGER  ::   ji, jj, jk             ! dummy loop indices
79      REAL(wp) ::   zu, zv                 ! temporary scalars
80      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
81      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
85      !
86      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
87      !
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
91         IF(lwp) WRITE(numout,*) '~~~~~~~'
92      ENDIF
93
94      IF( l_trddyn .or. ( id_vrt_dia_keg == "keg" ) ) THEN           ! Save ua and va trends
95         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
96         ztrdu(:,:,:) = ua(:,:,:) 
97         ztrdv(:,:,:) = va(:,:,:) 
98      ENDIF
99     
100      zhke(:,:,jpk) = 0._wp
101     
102      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
103      !
104      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
105         DO jk = 1, jpkm1
106            DO jj = 2, jpj
107               DO ji = fs_2, jpi   ! vector opt.
108                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
109                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
110                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
111                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
112                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
113               END DO 
114            END DO
115         END DO
116         !
117      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
118         DO jk = 1, jpkm1
119            DO jj = 2, jpjm1       
120               DO ji = fs_2, jpim1   ! vector opt.
121                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
122                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
123                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
124                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
125                     !
126                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
127                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
128                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
129                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
130                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
131               END DO 
132            END DO
133         END DO
134         CALL lbc_lnk( zhke, 'T', 1. )
135         !
136      END SELECT
137      !
138      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
139         DO jj = 2, jpjm1
140            DO ji = fs_2, fs_jpim1   ! vector opt.
141               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
142               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
143            END DO
144         END DO
145      END DO
146      !
147      IF( l_trddyn .or. ( id_vrt_dia_keg == "keg" ) ) THEN      ! save the Kinetic Energy trends for diagnostic
148         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
149         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
150         !
151         IF( id_vrt_dia_keg == "keg" ) THEN
152             CALL dyn_vrt_dia_3d( ztrdu, ztrdv, id_vrt_dia_keg )
153         END IF
154         !
155         IF( l_trddyn ) THEN
156             CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
157         END IF
158         !
159         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
160      ENDIF
161      !
162      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
163         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
164      !
165      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
166      !
167      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
168      !
169   END SUBROUTINE dyn_keg
170
171   !!======================================================================
172END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.