source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90

Last change on this file was 11738, checked in by marc, 18 months ago

The Dr Hook changes from my perl code.

File size: 16.1 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 zdf_oce        ! ocean vertical physics variables
16   USE trd_oce        ! trends: ocean variables
17!!gm   USE dynhpg          ! hydrostatic pressure gradient   
18   USE zdfbfr         ! bottom friction
19   USE ldftra_oce     ! ocean active tracers lateral physics
20   USE sbc_oce        ! surface boundary condition: ocean
21   USE phycst         ! physical constants
22   USE trdvor         ! ocean vorticity trends
23   USE trdglo         ! trends:global domain averaged
24   USE trdmxl         ! ocean active mixed layer tracers trends
25   USE in_out_manager ! I/O manager
26   USE iom            ! I/O manager library
27   USE lib_mpp        ! MPP library
28   USE wrk_nemo       ! Memory allocation
29   USE ldfslp         ! Isopycnal slopes
30
31   USE yomhook, ONLY: lhook, dr_hook
32   USE parkind1, ONLY: jprb, jpim
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   trd_ken       ! called by trddyn module
38   PUBLIC   trd_ken_init  ! called by trdini module
39
40   INTEGER ::   nkstp       ! current time step
41
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   bu, bv   ! volume of u- and v-boxes
43   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   r1_bt    ! inverse of t-box volume
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48#  include "ldfeiv_substitute.h90"
49
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   INTEGER FUNCTION trd_ken_alloc()
58   INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
59   INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
60   REAL(KIND=jprb)               :: zhook_handle
61
62   CHARACTER(LEN=*), PARAMETER :: RoutineName='TRD_KEN_ALLOC'
63
64   IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
65
66      !!---------------------------------------------------------------------
67      !!                  ***  FUNCTION trd_ken_alloc  ***
68      !!---------------------------------------------------------------------
69      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
70      !
71      IF( lk_mpp             )   CALL mpp_sum ( trd_ken_alloc )
72      IF( trd_ken_alloc /= 0 )   CALL ctl_warn('trd_ken_alloc: failed to allocate arrays')
73   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
74   END FUNCTION trd_ken_alloc
75
76
77   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt )
78      !!---------------------------------------------------------------------
79      !!                  ***  ROUTINE trd_ken  ***
80      !!
81      !! ** Purpose :   output 3D Kinetic Energy trends using IOM
82      !!
83      !! ** Method  : - apply lbc to the input masked velocity trends
84      !!              - compute the associated KE trend:
85      !!          zke = 0.5 * (  mi-1[ un * putrd * bu ] + mj-1[ vn * pvtrd * bv]  ) / bt
86      !!      where bu, bv, bt are the volume of u-, v- and t-boxes.
87      !!              - vertical diffusion case (jpdyn_zdf):
88      !!          diagnose separately the KE trend associated with wind stress
89      !!              - bottom friction case (jpdyn_bfr):
90      !!          explicit case (ln_bfrimp=F): bottom trend put in the 1st level
91      !!                                       of putrd, pvtrd
92      !
93      !
94      !!----------------------------------------------------------------------
95      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
96      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
97      INTEGER                   , INTENT(in   ) ::   kt             ! time step
98      !
99      INTEGER ::   ji, jj, jk       ! dummy loop indices
100      INTEGER ::   ikbu  , ikbv     ! local integers
101      INTEGER ::   ikbum1, ikbvm1   !   -       -
102      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, zke2d   ! 2D workspace
103      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace
104      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
105      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
106      REAL(KIND=jprb)               :: zhook_handle
107
108      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRD_KEN'
109
110      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
111
112      !!----------------------------------------------------------------------
113      !
114      CALL wrk_alloc( jpi, jpj, jpk, zke )
115      !
116      CALL lbc_lnk( putrd, 'U', -1. )   ;   CALL lbc_lnk( pvtrd, 'V', -1. )      ! lateral boundary conditions
117      !
118      IF ( lk_vvl .AND. kt /= nkstp ) THEN   ! Variable volume: set box volume at the 1st call of kt time step
119         nkstp = kt
120         DO jk = 1, jpkm1
121            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk)
122            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk)
123            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) ) * tmask(:,:,jk)
124         END DO
125      ENDIF
126      !
127      zke(:,:,jpk) = 0._wp
128      zke(1,:, : ) = 0._wp
129      zke(:,1, : ) = 0._wp
130      DO jk = 1, jpkm1
131         DO jj = 2, jpj
132            DO ji = 2, jpi
133               zke(ji,jj,jk) = 0.5_wp * rau0 *( un(ji  ,jj,jk) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
134                  &                           + un(ji-1,jj,jk) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
135                  &                           + vn(ji,jj  ,jk) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
136                  &                           + vn(ji,jj-1,jk) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
137            END DO
138         END DO
139      END DO
140      !
141      SELECT CASE( ktrd )
142        CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg", zke )    ! hydrostatic pressure gradient
143        CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg", zke )    ! surface pressure gradient
144        CASE( jpdyn_spgexp );   CALL iom_put( "ketrd_spgexp", zke ) ! surface pressure gradient (explicit)
145        CASE( jpdyn_spgflt );   CALL iom_put( "ketrd_spgflt", zke ) ! surface pressure gradient (filter)
146        CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo", zke )    ! planetary vorticity
147        CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo", zke )    ! relative  vorticity     (or metric term)
148        CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg", zke )    ! Kinetic Energy gradient (or had)
149        CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad", zke )    ! vertical   advection
150        CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf", zke )    ! lateral diffusion
151        CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf", zke )    ! vertical diffusion
152                                 !                                   ! wind stress trends
153                                CALL wrk_alloc( jpi, jpj, z2dx, z2dy, zke2d )
154                     z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1u(:,:) * e2u(:,:) * umask(:,:,1)
155                     z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1v(:,:) * e2v(:,:) * vmask(:,:,1)
156                     zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
157                     DO jj = 2, jpj
158                         DO ji = 2, jpi
159                           zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
160                            &                         + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
161                         END DO
162                     END DO
163                                CALL iom_put( "ketrd_tau", zke2d )
164                                CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, zke2d )
165        CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr", zke )    ! bottom friction (explicit case)
166!!gm TO BE DONE properly
167!!gm only valid if ln_bfrimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
168!         IF(.NOT. ln_bfrimp) THEN
169!            DO jj = 1, jpj    !   
170!               DO ji = 1, jpi
171!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
172!                  ikbv = mbkv(ji,jj)   
173!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu)
174!                  z2dy(ji,jj) = vn(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv)
175!               END DO
176!            END DO
177!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
178!            DO jj = 2, jpj
179!               DO ji = 2, jpi
180!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
181!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
182!               END DO
183!            END DO
184!                              CALL iom_put( "ketrd_bfr", zke2d )    ! bottom friction (explicit case)
185!         ENDIF
186!!gm end
187        CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf", zke )    ! asselin filter trends
188!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdyn_bfr avant le swap !!!
189!! reflechir a une possible sauvegarde du "vrai" un,vn pour le calcul de atf....
190!
191!         IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case)
192!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
193!               DO ji = 1, jpi
194!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
195!                  ikbv = mbkv(ji,jj)
196!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) / fse3u(ji,jj,ikbu)
197!                  z2dy(ji,jj) = un(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) / fse3v(ji,jj,ikbv)
198!               END DO
199!            END DO
200!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
201!            DO jj = 2, jpj
202!               DO ji = 2, jpi
203!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
204!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   )
205!               END DO
206!            END DO
207!                              CALL iom_put( "ketrd_bfri", zke2d )
208!         ENDIF
209        CASE( jpdyn_ken )   ;   ! kinetic energy
210                    ! called in dynnxt.F90 before asselin time filter
211                    ! with putrd=ua and pvtrd=va
212                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
213                    CALL iom_put( "KE", zke )
214                    !
215                    CALL ken_p2k( kt , zke )
216                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
217# if defined key_ldfslp || key_esopa
218        CASE( jpdyn_eivke )
219            ! CMIP6 diagnostic tknebto = tendency of KE from
220            ! parameterized mesoscale eddy advection
221            ! = vertical_integral( k (N S)^2 ) rho dz
222            ! rho = reference density
223            ! S = isoneutral slope.
224            ! Most terms are on W grid so work on this grid
225            CALL wrk_alloc( jpi, jpj, zke2d )
226            zke2d(:,:) = 0._wp
227            DO jk = 1,jpk
228               DO ji = 1,jpi
229                  DO jj = 1,jpj
230                     zke2d(ji,jj) = zke2d(ji,jj) +  rau0 * fsaeiw(ji, jj, jk)               &
231                          &                      * ( wslpi(ji, jj, jk) * wslpi(ji,jj,jk)    &
232                          &                      +   wslpj(ji, jj, jk) * wslpj(ji,jj,jk) )  &
233                          &                      *   rn2(ji,jj,jk) * fse3w(ji, jj, jk)
234                  ENDDO
235               ENDDO
236            ENDDO
237            CALL iom_put("ketrd_eiv", zke2d)
238            CALL wrk_dealloc( jpi, jpj, zke2d )
239#endif
240         !
241      END SELECT
242      !
243      CALL wrk_dealloc( jpi, jpj, jpk, zke )
244      !
245      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
246   END SUBROUTINE trd_ken
247
248
249   SUBROUTINE ken_p2k( kt , pconv )
250      !!---------------------------------------------------------------------
251      !!                 ***  ROUTINE ken_p2k  ***
252      !!                   
253      !! ** Purpose :   compute rate of conversion from potential to kinetic energy
254      !!
255      !! ** Method  : - compute conv defined as -rau*g*w on T-grid points
256      !!
257      !! ** Work only for full steps and partial steps (ln_hpg_zco or ln_hpg_zps)
258      !!----------------------------------------------------------------------
259      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
260      !!
261      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   pconv
262      !
263      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
264      INTEGER  ::   iku, ikv                         ! temporary integers
265      REAL(wp) ::   zcoef                            ! temporary scalars
266      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zconv  ! temporary conv on W-grid
267      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
268      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
269      REAL(KIND=jprb)               :: zhook_handle
270
271      CHARACTER(LEN=*), PARAMETER :: RoutineName='KEN_P2K'
272
273      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
274
275      !!----------------------------------------------------------------------
276      !
277      CALL wrk_alloc( jpi,jpj,jpk, zconv )
278      !
279      ! Local constant initialization
280      zcoef = - rau0 * grav * 0.5_wp     
281     
282      !  Surface value (also valid in partial step case)
283      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * wn(:,:,1) * fse3w(:,:,1)
284
285      ! interior value (2=<jk=<jpkm1)
286      DO jk = 2, jpk
287         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * wn(:,:,jk) * fse3w(:,:,jk)
288      END DO
289
290      ! conv value on T-point
291      DO jk = 1, jpkm1
292         DO jj = 1, jpj
293            DO ji = 1, jpi
294               zcoef = 0.5_wp / fse3t(ji,jj,jk)
295               pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
296            END DO
297         END DO
298      END DO
299      !
300      CALL wrk_dealloc( jpi,jpj,jpk, zconv )     
301      !
302      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
303   END SUBROUTINE ken_p2k
304
305
306   SUBROUTINE trd_ken_init
307      !!---------------------------------------------------------------------
308      !!                  ***  ROUTINE trd_ken_init  ***
309      !!
310      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
311      !!----------------------------------------------------------------------
312      INTEGER  ::   ji, jj, jk   ! dummy loop indices
313      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
314      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
315      REAL(KIND=jprb)               :: zhook_handle
316
317      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRD_KEN_INIT'
318
319      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
320
321      !!----------------------------------------------------------------------
322      !
323      IF(lwp) THEN
324         WRITE(numout,*)
325         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends'
326         WRITE(numout,*) '~~~~~~~~~~~~~'
327      ENDIF
328      !                           ! allocate box volume arrays
329      IF ( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
330      !
331!!gm      IF( .NOT. (ln_hpg_zco.OR.ln_hpg_zps) )   &
332!!gm         &   CALL ctl_stop('trd_ken_init : only full and partial cells are coded for conversion rate')
333      !
334      IF ( .NOT.lk_vvl ) THEN     ! constant volume: bu, bv, 1/bt computed one for all
335         DO jk = 1, jpkm1
336            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk)
337            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk)
338            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) )
339         END DO
340      ENDIF
341      !
342      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
343   END SUBROUTINE trd_ken_init
344
345   !!======================================================================
346END MODULE trdken
Note: See TracBrowser for help on using the repository browser.