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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo
  • 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   USE wrk_nemo        ! Memory Allocation
22   USE timing          ! Timing
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   dyn_keg    ! routine called by step module
28   
29   !! * Substitutions
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE dyn_keg( kt )
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE dyn_keg  ***
41      !!
42      !! ** Purpose :   Compute the now momentum trend due to the horizontal
43      !!      gradient of the horizontal kinetic energy and add it to the
44      !!      general momentum trend.
45      !!
46      !! ** Method  :   Compute the now horizontal kinetic energy
47      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
48      !!      Take its horizontal gradient and add it to the general momentum
49      !!      trend (ua,va).
50      !!         ua = ua - 1/e1u di[ zhke ]
51      !!         va = va - 1/e2v dj[ zhke ]
52      !!
53      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
54      !!             - save this trends (l_trddyn=T) for post-processing
55      !!----------------------------------------------------------------------
56      USE oce     , ONLY:   tsa             ! tsa used as 2 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(:,:,:) :: zhke
63      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
64      !!----------------------------------------------------------------------
65      !
66      IF( nn_timing == 1 )  CALL timing_start('dyn_keg')
67      !
68      CALL wrk_alloc( jpi, jpj, jpk, zhke )
69      !
70      IF( kt == nit000 ) THEN
71         IF(lwp) WRITE(numout,*)
72         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
73         IF(lwp) WRITE(numout,*) '~~~~~~~'
74      ENDIF
75
76      IF( l_trddyn ) THEN           ! Save ua and va trends
77         ztrdu => tsa(:,:,:,1) 
78         ztrdv => tsa(:,:,:,2) 
79         !
80         ztrdu(:,:,:) = ua(:,:,:) 
81         ztrdv(:,:,:) = va(:,:,:) 
82      ENDIF
83     
84      !                                                ! ===============
85      DO jk = 1, jpkm1                                 ! Horizontal slab
86         !                                             ! ===============
87         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
88            DO ji = fs_2, jpi   ! vector opt.
89               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
90                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
91               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
92                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
93               zhke(ji,jj,jk) = zv + zu
94!!gm simplier coding  ==>>   ~ faster
95!    don't forget to suppress local zu zv scalars
96!               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
97!                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
98!                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
99!                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
100!!gm end <<==
101            END DO 
102         END DO 
103         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
104            DO ji = fs_2, fs_jpim1   ! vector opt.
105               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
106               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
107            END DO
108         END DO
109!!gm idea to be tested  ==>>   is it faster on scalar computers ?
110!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
111!            DO ji = fs_2, fs_jpim1   ! vector opt.
112!               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   &
113!                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   &
114!                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   &
115!                  !
116!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
117!                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
118!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj)
119!                  !
120!               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   &
121!                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   &
122!                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   &
123!                  !
124!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
125!                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
126!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj)
127!            END DO
128!         END DO
129!!gm en idea            <<==
130         !                                             ! ===============
131      END DO                                           !   End of slab
132      !                                                ! ===============
133
134      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
135         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
136         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
137         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )
138      ENDIF
139      !
140      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
141         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
142      !
143      CALL wrk_dealloc( jpi, jpj, jpk, zhke )
144      !
145      IF( nn_timing == 1 )  CALL timing_stop('dyn_keg')
146      !
147   END SUBROUTINE dyn_keg
148
149   !!======================================================================
150END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.