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 branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

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