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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdken.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 12.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 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      CALL mpp_sum ( 'trdken', trd_ken_alloc )
57      IF( trd_ken_alloc /= 0 )   CALL ctl_stop( 'STOP', 'trd_ken_alloc: failed to allocate arrays' )
58   END FUNCTION trd_ken_alloc
59
60
61   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt, Kmm )
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[ uu(Kmm) * putrd * bu ] + mj-1[ vv(Kmm) * 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      INTEGER                   , INTENT(in   ) ::   Kmm            ! time level index
83      !
84      INTEGER ::   ji, jj, jk       ! dummy loop indices
85      INTEGER ::   ikbu  , ikbv     ! local integers
86      INTEGER ::   ikbum1, ikbvm1   !   -       -
87      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2dx, z2dy, zke2d   ! 2D workspace
88      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zke                 ! 3D workspace
89      !!----------------------------------------------------------------------
90      !
91      CALL lbc_lnk_multi( 'trdken', putrd, 'U', -1. , pvtrd, 'V', -1. )      ! lateral boundary conditions
92      !
93      nkstp = kt
94      DO jk = 1, jpkm1
95         bu   (:,:,jk) =    e1e2u(:,:) * e3u(:,:,jk,Kmm)
96         bv   (:,:,jk) =    e1e2v(:,:) * e3v(:,:,jk,Kmm)
97         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
98      END DO
99      !
100      zke(:,:,jpk) = 0._wp
101      zke(1,:, : ) = 0._wp
102      zke(:,1, : ) = 0._wp
103      DO jk = 1, jpkm1
104         DO jj = 2, jpj
105            DO ji = 2, jpi
106               zke(ji,jj,jk) = 0.5_wp * rau0 *( uu(ji  ,jj,jk,Kmm) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
107                  &                           + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
108                  &                           + vv(ji,jj  ,jk,Kmm) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
109                  &                           + vv(ji,jj-1,jk,Kmm) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
110            END DO
111         END DO
112      END DO
113      !
114      SELECT CASE( ktrd )
115         CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg"   , zke )    ! hydrostatic pressure gradient
116         CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg"   , zke )    ! surface pressure gradient
117         CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo"   , zke )    ! planetary vorticity
118         CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo"   , zke )    ! relative  vorticity     (or metric term)
119         CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg"   , zke )    ! Kinetic Energy gradient (or had)
120         CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad"   , zke )    ! vertical   advection
121         CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf"   , zke )    ! lateral diffusion
122         CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion
123         !                   !                                          ! wind stress trends
124                                 ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
125                           z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
126                           z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
127                           zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
128                           DO jj = 2, jpj
129                              DO ji = 2, jpi
130                                 zke2d(ji,jj) = r1_rau0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
131                                 &                                   + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
132                              END DO
133                           END DO
134                                 CALL iom_put( "ketrd_tau"   , zke2d )  !
135                                 DEALLOCATE( z2dx , z2dy , zke2d )
136         CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case)
137!!gm TO BE DONE properly
138!!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
139!         IF(.NOT. ln_drgimp) THEN
140!            DO jj = 1, jpj    !   
141!               DO ji = 1, jpi
142!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
143!                  ikbv = mbkv(ji,jj)   
144!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
145!                  z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
146!               END DO
147!            END DO
148!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
149!            DO jj = 2, jpj
150!               DO ji = 2, jpi
151!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
152!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
153!               END DO
154!            END DO
155!                                    CALL iom_put( "ketrd_bfr"  , zke2d )   ! bottom friction (explicit case)
156!         ENDIF
157!!gm end
158         CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf"   , zke )    ! asselin filter trends
159!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdyn_bfr avant le swap !!!
160!! reflechir a une possible sauvegarde du "vrai" uu(Kmm),vv(Kmm) pour le calcul de atf....
161!
162!         IF( ln_drgimp ) THEN                                          ! bottom friction (implicit case)
163!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
164!               DO ji = 1, jpi
165!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
166!                  ikbv = mbkv(ji,jj)
167!                  z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm) / e3u(ji,jj,ikbu,Kmm)
168!                  z2dy(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm) / e3v(ji,jj,ikbv,Kmm)
169!               END DO
170!            END DO
171!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
172!            DO jj = 2, jpj
173!               DO ji = 2, jpi
174!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
175!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   )
176!               END DO
177!            END DO
178!                              CALL iom_put( "ketrd_bfri", zke2d )
179!         ENDIF
180        CASE( jpdyn_ken )   ;   ! kinetic energy
181                    ! called in dynnxt.F90 before asselin time filter
182                    ! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
183                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
184                    CALL iom_put( "KE", zke )
185                    !
186                    CALL ken_p2k( kt , zke, Kmm )
187                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
188         !
189      END SELECT
190      !
191   END SUBROUTINE trd_ken
192
193
194   SUBROUTINE ken_p2k( kt , pconv, Kmm )
195      !!---------------------------------------------------------------------
196      !!                 ***  ROUTINE ken_p2k  ***
197      !!                   
198      !! ** Purpose :   compute rate of conversion from potential to kinetic energy
199      !!
200      !! ** Method  : - compute conv defined as -rau*g*w on T-grid points
201      !!
202      !! ** Work only for full steps and partial steps (ln_hpg_zco or ln_hpg_zps)
203      !!----------------------------------------------------------------------
204      INTEGER                   , INTENT(in   ) ::   kt      ! ocean time-step index
205      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
206      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pconv   !
207      !
208      INTEGER  ::   ji, jj, jk   ! dummy loop indices
209      INTEGER  ::   iku, ikv     ! local integers
210      REAL(wp) ::   zcoef        ! local scalars
211      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zconv  ! 3D workspace
212      !!----------------------------------------------------------------------
213      !
214      ! Local constant initialization
215      zcoef = - rau0 * grav * 0.5_wp     
216     
217      !  Surface value (also valid in partial step case)
218      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
219
220      ! interior value (2=<jk=<jpkm1)
221      DO jk = 2, jpk
222         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
223      END DO
224
225      ! conv value on T-point
226      DO jk = 1, jpkm1
227         DO jj = 1, jpj
228            DO ji = 1, jpi
229               zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
230               pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
231            END DO
232         END DO
233      END DO
234      !
235   END SUBROUTINE ken_p2k
236
237
238   SUBROUTINE trd_ken_init
239      !!---------------------------------------------------------------------
240      !!                  ***  ROUTINE trd_ken_init  ***
241      !!
242      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
243      !!----------------------------------------------------------------------
244      INTEGER  ::   ji, jj, jk   ! dummy loop indices
245      !!----------------------------------------------------------------------
246      !
247      IF(lwp) THEN
248         WRITE(numout,*)
249         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends'
250         WRITE(numout,*) '~~~~~~~~~~~~~'
251      ENDIF
252      !                           ! allocate box volume arrays
253      IF( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
254      !
255   END SUBROUTINE trd_ken_init
256
257   !!======================================================================
258END MODULE trdken
Note: See TracBrowser for help on using the repository browser.