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 trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.5 KB
Line 
1MODULE dynkeg
2   !!======================================================================
3   !!                       ***  MODULE  dynkeg  ***
4   !! Ocean dynamics:  kinetic energy gradient trend
5   !!======================================================================
6   
7   !!----------------------------------------------------------------------
8   !!   dyn_keg      : update the momentum trend with the horizontal tke
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers
12   USE dom_oce         ! ocean space and time domain
13   USE in_out_manager  ! I/O manager
14   USE trdmod          ! ocean dynamics trends
15   USE trdmod_oce      ! ocean variables trends
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Accessibility
21   PUBLIC dyn_keg                ! routine called by step.F90
22   
23   !! * Substitutions
24#  include "vectopt_loop_substitute.h90"
25   !!---------------------------------------------------------------------------------
26   !!   OPA 9.0 , LOCEAN-IPSL (2005)
27   !! $Header$
28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
29   !!---------------------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE dyn_keg( kt )
34      !!----------------------------------------------------------------------
35      !!                  ***  ROUTINE dyn_keg  ***
36      !!
37      !! ** Purpose :   Compute the now momentum trend due to the horizontal
38      !!      gradient of the horizontal kinetic energy and add it to the
39      !!      general momentum trend.
40      !!
41      !! ** Method  :   Compute the now horizontal kinetic energy:
42      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
43      !!      Take its horizontal gradient and add it to the general momentum
44      !!      trend (ua,va).
45      !!         ua = ua - 1/e1u di[ zhke ]
46      !!         va = va - 1/e2v dj[ zhke ]
47      !!
48      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
49      !!             - Save the trends in (utrd,vtrd) ('key_trddyn')
50      !!
51      !! History :
52      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
53      !!   7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
54      !!   9.0  !  02-07  (G. Madec)  F90: Free form and module
55      !!    "   !  04-08  (C. Talandier) New trends organization
56      !!----------------------------------------------------------------------
57      !! * Modules used     
58      USE oce, ONLY :    ztdua => ta,   & ! use ta as 3D workspace   
59                         ztdva => sa      ! use sa as 3D workspace   
60
61      !! * Arguments
62      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
63
64      !! * Local declarations
65      INTEGER  ::   ji, jj, jk          ! dummy loop indices
66      REAL(wp) ::   zua, zva, zu, zv    ! temporary scalars
67      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
68         zhke                           ! temporary workspace
69      !!----------------------------------------------------------------------
70
71      IF( kt == nit000 ) THEN
72         IF(lwp) WRITE(numout,*)
73         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
74         IF(lwp) WRITE(numout,*) '~~~~~~~'
75      ENDIF
76
77      ! Save ua and va trends
78      IF( l_trddyn )   THEN
79         ztdua(:,:,:) = ua(:,:,:) 
80         ztdva(:,:,:) = va(:,:,:) 
81      ENDIF
82     
83      !                                                ! ===============
84      DO jk = 1, jpkm1                                 ! Horizontal slab
85         !                                             ! ===============
86         ! Horizontal kinetic energy at T-point
87         DO jj = 2, jpj
88            DO ji = fs_2, jpi   ! vector opt.
89               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
90                            + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
91               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
92                            + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
93               zhke(ji,jj,jk) = zv + zu
94            END DO 
95         END DO 
96         
97         ! Horizontal gradient of Horizontal kinetic energy
98         DO jj = 2, jpjm1
99            DO ji = fs_2, fs_jpim1   ! vector opt.
100               ! gradient of kinetic energy
101               zua = -( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
102               zva = -( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
103               ! add to the general momentum trends
104               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
105               va(ji,jj,jk) = va(ji,jj,jk) + zva
106            END DO
107         END DO
108         !                                             ! ===============
109      END DO                                           !   End of slab
110      !                                                ! ===============
111
112      ! save the Kinetic Energy trends for diagnostic
113      ! momentum trends
114      IF( l_trddyn )   THEN
115         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
116         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
117
118         CALL trd_mod(ztdua, ztdva, jpdtdkeg, 'DYN', kt)
119      ENDIF
120
121      IF(l_ctl) THEN         ! print sum trends (used for debugging)
122         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
123         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
124         WRITE(numout,*) ' keg  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
125         u_ctl = zua   ;   v_ctl = zva
126      ENDIF
127
128   END SUBROUTINE dyn_keg
129
130   !!======================================================================
131END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.