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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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