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/trunk/src/OCE/TRD – NEMO

source: NEMO/trunk/src/OCE/TRD/trdken.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

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