source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/TRD/trdken.F90 @ 10170

Last change on this file since 10170 was 10170, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of lbc_lnk, see #2133

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