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/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90 @ 10047

Last change on this file since 10047 was 10047, checked in by jpalmier, 6 years ago

merge with GO6_package_branch 9385-10020 ; plus debug OMIP_DIC

File size: 14.7 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      ENDIF
290      !                           ! allocate box volume arrays
291      IF ( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
292      !
293!!gm      IF( .NOT. (ln_hpg_zco.OR.ln_hpg_zps) )   &
294!!gm         &   CALL ctl_stop('trd_ken_init : only full and partial cells are coded for conversion rate')
295      !
296      IF ( .NOT.lk_vvl ) THEN     ! constant volume: bu, bv, 1/bt computed one for all
297         DO jk = 1, jpkm1
298            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk)
299            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk)
300            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) )
301         END DO
302      ENDIF
303      !
304   END SUBROUTINE trd_ken_init
305
306   !!======================================================================
307END MODULE trdken
Note: See TracBrowser for help on using the repository browser.