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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

Ediag branche: #927 split TRA/DYN trd computation
  • 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 trd_oce         ! trends: ocean variables
17   USE trddyn          ! trend manager: dynamics
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      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
55      !!----------------------------------------------------------------------
56      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
57      !!
58      INTEGER  ::   ji, jj, jk   ! dummy loop indices
59      REAL(wp) ::   zu, zv       ! temporary scalars
60      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
61      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
62      !!----------------------------------------------------------------------
63      !
64      IF( nn_timing == 1 )  CALL timing_start('dyn_keg')
65      !
66      CALL wrk_alloc( jpi, jpj, jpk, zhke )
67      !
68      IF( kt == nit000 ) THEN
69         IF(lwp) WRITE(numout,*)
70         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
71         IF(lwp) WRITE(numout,*) '~~~~~~~'
72      ENDIF
73
74      IF( l_trddyn ) THEN           ! Save ua and va trends
75         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
76         ztrdu(:,:,:) = ua(:,:,:) 
77         ztrdv(:,:,:) = va(:,:,:) 
78      ENDIF
79     
80      !                                                ! ===============
81      DO jk = 1, jpkm1                                 ! Horizontal slab
82         !                                             ! ===============
83         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
84            DO ji = fs_2, jpi   ! vector opt.
85               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
86                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
87               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
88                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
89               zhke(ji,jj,jk) = zv + zu
90!!gm simplier coding  ==>>   ~ faster
91!    don't forget to suppress local zu zv scalars
92!               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
93!                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
94!                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
95!                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
96!!gm end <<==
97            END DO 
98         END DO 
99         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
100            DO ji = fs_2, fs_jpim1   ! vector opt.
101               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
102               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
103            END DO
104         END DO
105!!gm idea to be tested  ==>>   is it faster on scalar computers ?
106!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
107!            DO ji = fs_2, fs_jpim1   ! vector opt.
108!               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   &
109!                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   &
110!                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   &
111!                  !
112!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
113!                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
114!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj)
115!                  !
116!               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   &
117!                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   &
118!                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   &
119!                  !
120!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
121!                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
122!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj)
123!            END DO
124!         END DO
125!!gm en idea            <<==
126         !                                             ! ===============
127      END DO                                           !   End of slab
128      !                                                ! ===============
129
130      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
131         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
132         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
133         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
134         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
135      ENDIF
136      !
137      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
138         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
139      !
140      CALL wrk_dealloc( jpi, jpj, jpk, zhke )
141      !
142      IF( nn_timing == 1 )  CALL timing_stop('dyn_keg')
143      !
144   END SUBROUTINE dyn_keg
145
146   !!======================================================================
147END MODULE dynkeg
Note: See TracBrowser for help on using the repository browser.