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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90 @ 4400

Last change on this file since 4400 was 3432, checked in by trackstand2, 12 years ago

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