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/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

Implementation of the merge of TRA/TRP : first guess, see ticket #842
  • Property svn:keywords set to Id
File size: 7.3 KB
Line 
1MODULE dynkeg
2   !!======================================================================
3   !!                       ***  MODULE  dynkeg  ***
4   !! Ocean dynamics:  kinetic energy gradient trend
5   !!======================================================================
6   !! History :  1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
7   !!            7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
8   !!            9.0  !  02-07  (G. Madec)  F90: Free form and module
9   !!----------------------------------------------------------------------
10   
11   !!----------------------------------------------------------------------
12   !!   dyn_keg      : update the momentum trend with the horizontal tke
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE trdmod          ! ocean dynamics trends
17   USE trdmod_oce      ! ocean variables trends
18   USE in_out_manager  ! I/O manager
19   USE lib_mpp         ! MPP library
20   USE prtctl          ! Print control
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   dyn_keg    ! routine called by step module
26   
27   !! * Substitutions
28#  include "vectopt_loop_substitute.h90"
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE dyn_keg( kt )
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE dyn_keg  ***
39      !!
40      !! ** Purpose :   Compute the now momentum trend due to the horizontal
41      !!      gradient of the horizontal kinetic energy and add it to the
42      !!      general momentum trend.
43      !!
44      !! ** Method  :   Compute the now horizontal kinetic energy
45      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
46      !!      Take its horizontal gradient and add it to the general momentum
47      !!      trend (ua,va).
48      !!         ua = ua - 1/e1u di[ zhke ]
49      !!         va = va - 1/e2v dj[ zhke ]
50      !!
51      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
52      !!             - save this trends (l_trddyn=T) for post-processing
53      !!----------------------------------------------------------------------
54      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
55      USE oce     , ONLY:   tsa             ! tsa used as 2 3D workspace
56      USE wrk_nemo, ONLY:   zhke  => wrk_3d_1                 ! 3D workspace
57      !!
58      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
59      !!
60      INTEGER  ::   ji, jj, jk   ! dummy loop indices
61      REAL(wp) ::   zu, zv       ! temporary scalars
62      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
63      !!----------------------------------------------------------------------
64
65      IF( wrk_in_use(3,1) ) THEN
66         CALL ctl_stop('dyn_keg: requested workspace array is unavailable')   ;   RETURN
67      ENDIF
68
69      IF( kt == nit000 ) THEN
70         IF(lwp) WRITE(numout,*)
71         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
72         IF(lwp) WRITE(numout,*) '~~~~~~~'
73      ENDIF
74
75      IF( l_trddyn ) THEN           ! Save ua and va trends
76         ztrdu => tsa(:,:,:,1) 
77         ztrdv => tsa(:,:,:,2) 
78         !
79         ztrdu(:,:,:) = ua(:,:,:) 
80         ztrdv(:,:,:) = va(:,:,:) 
81      ENDIF
82     
83      !                                                ! ===============
84      DO jk = 1, jpkm1                                 ! Horizontal slab
85         !                                             ! ===============
86         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
87            DO ji = fs_2, jpi   ! vector opt.
88               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
89                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
90               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
91                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
92               zhke(ji,jj,jk) = zv + zu
93!!gm simplier coding  ==>>   ~ faster
94!    don't forget to suppress local zu zv scalars
95!               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
96!                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
97!                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
98!                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
99!!gm end <<==
100            END DO 
101         END DO 
102         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
103            DO ji = fs_2, fs_jpim1   ! vector opt.
104               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
105               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
106            END DO
107         END DO
108!!gm idea to be tested  ==>>   is it faster on scalar computers ?
109!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
110!            DO ji = fs_2, fs_jpim1   ! vector opt.
111!               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   &
112!                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   &
113!                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   &
114!                  !
115!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
116!                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
117!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj)
118!                  !
119!               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   &
120!                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   &
121!                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   &
122!                  !
123!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
124!                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
125!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj)
126!            END DO
127!         END DO
128!!gm en idea            <<==
129         !                                             ! ===============
130      END DO                                           !   End of slab
131      !                                                ! ===============
132
133      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
134         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
135         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
136         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )
137      ENDIF
138      !
139      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
140         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
141      !
142      IF( wrk_not_released(3, 1) )   CALL ctl_stop('dyn_keg: failed to release workspace array')
143      !
144   END SUBROUTINE dyn_keg
145
146   !!======================================================================
147END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.