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.
trdken.F90 in NEMO/trunk/src/OCE/TRD – NEMO

source: NEMO/trunk/src/OCE/TRD/trdken.F90

Last change on this file was 15104, checked in by clem, 3 years ago

nn_hls=2: repair some loops.

  • Property svn:keywords set to Id
File size: 12.6 KB
RevLine 
[4619]1MODULE trdken
2   !!======================================================================
3   !!                       ***  MODULE  trdken  ***
4   !! Ocean diagnostics:  compute and output 3D kinetic energy trends
5   !!=====================================================================
6   !! History :  3.5  !  2012-02  (G. Madec) original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   trd_ken       : compute and output 3D Kinetic energy trends using IOM
11   !!   trd_ken_init  : initialisation
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers variables
14   USE dom_oce        ! ocean space and time domain variables
[9019]15   USE phycst         ! physical constants
[5836]16   USE sbc_oce        ! surface boundary condition: ocean
[4619]17   USE zdf_oce        ! ocean vertical physics variables
[9256]18!!gm   USE zdfdrg         ! ocean vertical physics: bottom friction
[9019]19   USE ldftra         ! ocean active tracers lateral physics
[4619]20   USE trd_oce        ! trends: ocean variables
21   USE trdvor         ! ocean vorticity trends
22   USE trdglo         ! trends:global domain averaged
[5836]23   USE trdmxl         ! ocean active mixed layer tracers trends
24   !
[4619]25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager library
27   USE lib_mpp        ! MPP library
[7646]28   USE ldfslp         ! Isopycnal slopes
[4619]29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   trd_ken       ! called by trddyn module
34   PUBLIC   trd_ken_init  ! called by trdini module
35
36   INTEGER ::   nkstp       ! current time step
37
38   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   bu, bv   ! volume of u- and v-boxes
39   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   r1_bt    ! inverse of t-box volume
40
41   !! * Substitutions
[12377]42#  include "do_loop_substitute.h90"
[13237]43#  include "domzgr_substitute.h90"
[4619]44   !!----------------------------------------------------------------------
[9598]45   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]46   !! $Id$
[10068]47   !! Software governed by the CeCILL license (see ./LICENSE)
[4619]48   !!----------------------------------------------------------------------
49CONTAINS
50
51   INTEGER FUNCTION trd_ken_alloc()
52      !!---------------------------------------------------------------------
53      !!                  ***  FUNCTION trd_ken_alloc  ***
54      !!---------------------------------------------------------------------
55      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
56      !
[10425]57      CALL mpp_sum ( 'trdken', trd_ken_alloc )
58      IF( trd_ken_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trd_ken_alloc: failed to allocate arrays' )
[4619]59   END FUNCTION trd_ken_alloc
60
61
[12377]62   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt, Kmm )
[4619]63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE trd_ken  ***
65      !!
66      !! ** Purpose :   output 3D Kinetic Energy trends using IOM
67      !!
68      !! ** Method  : - apply lbc to the input masked velocity trends
69      !!              - compute the associated KE trend:
[12377]70      !!          zke = 0.5 * (  mi-1[ uu(Kmm) * putrd * bu ] + mj-1[ vv(Kmm) * pvtrd * bv]  ) / bt
[4619]71      !!      where bu, bv, bt are the volume of u-, v- and t-boxes.
72      !!              - vertical diffusion case (jpdyn_zdf):
73      !!          diagnose separately the KE trend associated with wind stress
74      !!              - bottom friction case (jpdyn_bfr):
[9019]75      !!          explicit case (ln_drgimp=F): bottom trend put in the 1st level
[4619]76      !!                                       of putrd, pvtrd
77      !
78      !
79      !!----------------------------------------------------------------------
80      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
81      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
82      INTEGER                   , INTENT(in   ) ::   kt             ! time step
[12377]83      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
[4619]84      !
85      INTEGER ::   ji, jj, jk       ! dummy loop indices
86      INTEGER ::   ikbu  , ikbv     ! local integers
87      INTEGER ::   ikbum1, ikbvm1   !   -       -
[9019]88      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2dx, z2dy, zke2d   ! 2D workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zke                 ! 3D workspace
[4619]90      !!----------------------------------------------------------------------
91      !
[14433]92      CALL lbc_lnk( 'trdken', putrd, 'U', -1.0_wp , pvtrd, 'V', -1.0_wp )      ! lateral boundary conditions
[4619]93      !
[6140]94      nkstp = kt
95      DO jk = 1, jpkm1
[12377]96         bu   (:,:,jk) =    e1e2u(:,:) * e3u(:,:,jk,Kmm)
97         bv   (:,:,jk) =    e1e2v(:,:) * e3v(:,:,jk,Kmm)
98         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
[6140]99      END DO
[4619]100      !
101      zke(:,:,jpk) = 0._wp
[15104]102      zke(1:nn_hls,:, : ) = 0._wp
103      zke(:,1:nn_hls, : ) = 0._wp
104      DO_3D( 0, nn_hls, 0, nn_hls, 1, jpkm1 )
[12489]105         zke(ji,jj,jk) = 0.5_wp * rho0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
[12377]106            &                           + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
107            &                           + vv(ji,jj  ,jk,Kmm) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
108            &                           + vv(ji,jj-1,jk,Kmm) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
109      END_3D
[4619]110      !
111      SELECT CASE( ktrd )
[6140]112         CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg"   , zke )    ! hydrostatic pressure gradient
113         CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg"   , zke )    ! surface pressure gradient
114         CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo"   , zke )    ! planetary vorticity
115         CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo"   , zke )    ! relative  vorticity     (or metric term)
116         CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg"   , zke )    ! Kinetic Energy gradient (or had)
117         CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad"   , zke )    ! vertical   advection
118         CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf"   , zke )    ! lateral diffusion
119         CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion
120         !                   !                                          ! wind stress trends
[9019]121                                 ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
[12377]122                           z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
123                           z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
[15104]124                           zke2d(1:nn_hls,:) = 0._wp   ;   zke2d(:,1:nn_hls) = 0._wp
125                           DO_2D( 0, nn_hls, 0, nn_hls )
[12489]126                              zke2d(ji,jj) = r1_rho0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
[12377]127                              &                                   + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
128                           END_2D
[6140]129                                 CALL iom_put( "ketrd_tau"   , zke2d )  !
[9019]130                                 DEALLOCATE( z2dx , z2dy , zke2d )
[6140]131         CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case)
[4619]132!!gm TO BE DONE properly
[9019]133!!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
134!         IF(.NOT. ln_drgimp) THEN
[4619]135!            DO jj = 1, jpj    !   
136!               DO ji = 1, jpi
137!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
138!                  ikbv = mbkv(ji,jj)   
[12377]139!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
140!                  z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
[4619]141!               END DO
142!            END DO
143!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
144!            DO jj = 2, jpj
145!               DO ji = 2, jpi
146!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
147!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
148!               END DO
149!            END DO
[6140]150!                                    CALL iom_put( "ketrd_bfr"  , zke2d )   ! bottom friction (explicit case)
[4619]151!         ENDIF
152!!gm end
[6140]153         CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf"   , zke )    ! asselin filter trends
[4619]154!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdyn_bfr avant le swap !!!
[12377]155!! reflechir a une possible sauvegarde du "vrai" uu(Kmm),vv(Kmm) pour le calcul de atf....
[4619]156!
[9019]157!         IF( ln_drgimp ) THEN                                          ! bottom friction (implicit case)
[4619]158!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
159!               DO ji = 1, jpi
160!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
161!                  ikbv = mbkv(ji,jj)
[12377]162!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm)
163!                  z2dy(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm)
[4619]164!               END DO
165!            END DO
166!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
167!            DO jj = 2, jpj
168!               DO ji = 2, jpi
169!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
170!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   )
171!               END DO
172!            END DO
173!                              CALL iom_put( "ketrd_bfri", zke2d )
174!         ENDIF
[7646]175        CASE( jpdyn_ken )   ;   ! kinetic energy
176                    ! called in dynnxt.F90 before asselin time filter
[12377]177                    ! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
[7646]178                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
179                    CALL iom_put( "KE", zke )
180                    !
[12377]181                    CALL ken_p2k( kt , zke, Kmm )
[7646]182                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
[4619]183         !
184      END SELECT
185      !
186   END SUBROUTINE trd_ken
187
188
[12377]189   SUBROUTINE ken_p2k( kt , pconv, Kmm )
[4619]190      !!---------------------------------------------------------------------
191      !!                 ***  ROUTINE ken_p2k  ***
192      !!                   
193      !! ** Purpose :   compute rate of conversion from potential to kinetic energy
194      !!
195      !! ** Method  : - compute conv defined as -rau*g*w on T-grid points
196      !!
197      !! ** Work only for full steps and partial steps (ln_hpg_zco or ln_hpg_zps)
198      !!----------------------------------------------------------------------
[9019]199      INTEGER                   , INTENT(in   ) ::   kt      ! ocean time-step index
[12377]200      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
[9019]201      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pconv   !
[4619]202      !
[9019]203      INTEGER  ::   ji, jj, jk   ! dummy loop indices
204      INTEGER  ::   iku, ikv     ! local integers
205      REAL(wp) ::   zcoef        ! local scalars
206      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zconv  ! 3D workspace
[4619]207      !!----------------------------------------------------------------------
208      !
209      ! Local constant initialization
[12489]210      zcoef = - rho0 * grav * 0.5_wp     
[4619]211     
212      !  Surface value (also valid in partial step case)
[12377]213      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
[4619]214
215      ! interior value (2=<jk=<jpkm1)
216      DO jk = 2, jpk
[12377]217         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
[4619]218      END DO
219
220      ! conv value on T-point
[15104]221      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
[12377]222         zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
223         pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
224      END_3D
[4619]225      !
226   END SUBROUTINE ken_p2k
227
228
229   SUBROUTINE trd_ken_init
230      !!---------------------------------------------------------------------
231      !!                  ***  ROUTINE trd_ken_init  ***
232      !!
233      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
234      !!----------------------------------------------------------------------
235      INTEGER  ::   ji, jj, jk   ! dummy loop indices
236      !!----------------------------------------------------------------------
237      !
238      IF(lwp) THEN
239         WRITE(numout,*)
240         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends'
241         WRITE(numout,*) '~~~~~~~~~~~~~'
242      ENDIF
243      !                           ! allocate box volume arrays
[5836]244      IF( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
[4619]245      !
246   END SUBROUTINE trd_ken_init
247
248   !!======================================================================
249END MODULE trdken
Note: See TracBrowser for help on using the repository browser.