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

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90 @ 8733

Last change on this file since 8733 was 8733, checked in by dancopsey, 6 years ago

Remove svn keywords.

File size: 12.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 sbc_oce        ! surface boundary condition: ocean
16   USE zdf_oce        ! ocean vertical physics variables
17   USE trd_oce        ! trends: ocean variables
18!!gm   USE dynhpg          ! hydrostatic pressure gradient   
19   USE zdfbfr         ! bottom friction
20   USE ldftra         ! ocean active tracers lateral physics
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   !
26   USE in_out_manager ! I/O manager
27   USE iom            ! I/O manager library
28   USE lib_mpp        ! MPP library
29   USE wrk_nemo       ! Memory allocation
30   USE ldfslp         ! Isopycnal slopes
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   trd_ken       ! called by trddyn module
36   PUBLIC   trd_ken_init  ! called by trdini module
37
38   INTEGER ::   nkstp       ! current time step
39
40   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   bu, bv   ! volume of u- and v-boxes
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   r1_bt    ! inverse of t-box volume
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION trd_ken_alloc()
53      !!---------------------------------------------------------------------
54      !!                  ***  FUNCTION trd_ken_alloc  ***
55      !!---------------------------------------------------------------------
56      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
57      !
58      IF( lk_mpp             )   CALL mpp_sum ( trd_ken_alloc )
59      IF( trd_ken_alloc /= 0 )   CALL ctl_warn('trd_ken_alloc: failed to allocate arrays')
60   END FUNCTION trd_ken_alloc
61
62
63   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE trd_ken  ***
66      !!
67      !! ** Purpose :   output 3D Kinetic Energy trends using IOM
68      !!
69      !! ** Method  : - apply lbc to the input masked velocity trends
70      !!              - compute the associated KE trend:
71      !!          zke = 0.5 * (  mi-1[ un * putrd * bu ] + mj-1[ vn * pvtrd * bv]  ) / bt
72      !!      where bu, bv, bt are the volume of u-, v- and t-boxes.
73      !!              - vertical diffusion case (jpdyn_zdf):
74      !!          diagnose separately the KE trend associated with wind stress
75      !!              - bottom friction case (jpdyn_bfr):
76      !!          explicit case (ln_bfrimp=F): bottom trend put in the 1st level
77      !!                                       of putrd, pvtrd
78      !
79      !
80      !!----------------------------------------------------------------------
81      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
82      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
83      INTEGER                   , INTENT(in   ) ::   kt             ! time step
84      !
85      INTEGER ::   ji, jj, jk       ! dummy loop indices
86      INTEGER ::   ikbu  , ikbv     ! local integers
87      INTEGER ::   ikbum1, ikbvm1   !   -       -
88      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, zke2d   ! 2D workspace
89      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace
90      !!----------------------------------------------------------------------
91      !
92      CALL wrk_alloc( jpi, jpj, jpk, zke )
93      !
94      CALL lbc_lnk( putrd, 'U', -1. )   ;   CALL lbc_lnk( pvtrd, 'V', -1. )      ! lateral boundary conditions
95      !
96      nkstp = kt
97      DO jk = 1, jpkm1
98         bu   (:,:,jk) =    e1e2u(:,:) * e3u_n(:,:,jk)
99         bv   (:,:,jk) =    e1e2v(:,:) * e3v_n(:,:,jk)
100         r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t_n(:,:,jk) * tmask(:,:,jk)
101      END DO
102      !
103      zke(:,:,jpk) = 0._wp
104      zke(1,:, : ) = 0._wp
105      zke(:,1, : ) = 0._wp
106      DO jk = 1, jpkm1
107         DO jj = 2, jpj
108            DO ji = 2, jpi
109               zke(ji,jj,jk) = 0.5_wp * rau0 *( un(ji  ,jj,jk) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
110                  &                           + un(ji-1,jj,jk) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
111                  &                           + vn(ji,jj  ,jk) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
112                  &                           + vn(ji,jj-1,jk) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
113            END DO
114         END DO
115      END DO
116      !
117      SELECT CASE( ktrd )
118         CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg"   , zke )    ! hydrostatic pressure gradient
119         CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg"   , zke )    ! surface pressure gradient
120         CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo"   , zke )    ! planetary vorticity
121         CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo"   , zke )    ! relative  vorticity     (or metric term)
122         CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg"   , zke )    ! Kinetic Energy gradient (or had)
123         CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad"   , zke )    ! vertical   advection
124         CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf"   , zke )    ! lateral diffusion
125         CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion
126         !                   !                                          ! wind stress trends
127                                 CALL wrk_alloc( jpi, jpj, z2dx, z2dy, zke2d )
128                           z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
129                           z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
130                           zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
131                           DO jj = 2, jpj
132                              DO ji = 2, jpi
133                                 zke2d(ji,jj) = r1_rau0 * 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
134                                 &                                   + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
135                              END DO
136                           END DO
137                                 CALL iom_put( "ketrd_tau"   , zke2d )  !
138                                 CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, zke2d )
139         CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case)
140!!gm TO BE DONE properly
141!!gm only valid if ln_bfrimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
142!         IF(.NOT. ln_bfrimp) THEN
143!            DO jj = 1, jpj    !   
144!               DO ji = 1, jpi
145!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
146!                  ikbv = mbkv(ji,jj)   
147!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu)
148!                  z2dy(ji,jj) = vn(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv)
149!               END DO
150!            END DO
151!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
152!            DO jj = 2, jpj
153!               DO ji = 2, jpi
154!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
155!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
156!               END DO
157!            END DO
158!                                    CALL iom_put( "ketrd_bfr"  , zke2d )   ! bottom friction (explicit case)
159!         ENDIF
160!!gm end
161         CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf"   , zke )    ! asselin filter trends
162!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdyn_bfr avant le swap !!!
163!! reflechir a une possible sauvegarde du "vrai" un,vn pour le calcul de atf....
164!
165!         IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case)
166!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
167!               DO ji = 1, jpi
168!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
169!                  ikbv = mbkv(ji,jj)
170!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
171!                  z2dy(ji,jj) = un(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
172!               END DO
173!            END DO
174!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
175!            DO jj = 2, jpj
176!               DO ji = 2, jpi
177!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
178!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   )
179!               END DO
180!            END DO
181!                              CALL iom_put( "ketrd_bfri", zke2d )
182!         ENDIF
183        CASE( jpdyn_ken )   ;   ! kinetic energy
184                    ! called in dynnxt.F90 before asselin time filter
185                    ! with putrd=ua and pvtrd=va
186                    zke(:,:,:) = 0.5_wp * zke(:,:,:)
187                    CALL iom_put( "KE", zke )
188                    !
189                    CALL ken_p2k( kt , zke )
190                      CALL iom_put( "ketrd_convP2K", zke )     ! conversion -rau*g*w
191         !
192      END SELECT
193      !
194      CALL wrk_dealloc( jpi, jpj, jpk, zke )
195      !
196   END SUBROUTINE trd_ken
197
198
199   SUBROUTINE ken_p2k( kt , pconv )
200      !!---------------------------------------------------------------------
201      !!                 ***  ROUTINE ken_p2k  ***
202      !!                   
203      !! ** Purpose :   compute rate of conversion from potential to kinetic energy
204      !!
205      !! ** Method  : - compute conv defined as -rau*g*w on T-grid points
206      !!
207      !! ** Work only for full steps and partial steps (ln_hpg_zco or ln_hpg_zps)
208      !!----------------------------------------------------------------------
209      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
210      !!
211      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   pconv
212      !
213      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
214      INTEGER  ::   iku, ikv                         ! temporary integers
215      REAL(wp) ::   zcoef                            ! temporary scalars
216      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zconv  ! temporary conv on W-grid
217      !!----------------------------------------------------------------------
218      !
219      CALL wrk_alloc( jpi,jpj,jpk, zconv )
220      !
221      ! Local constant initialization
222      zcoef = - rau0 * grav * 0.5_wp     
223     
224      !  Surface value (also valid in partial step case)
225      zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * wn(:,:,1) * e3w_n(:,:,1)
226
227      ! interior value (2=<jk=<jpkm1)
228      DO jk = 2, jpk
229         zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * wn(:,:,jk) * e3w_n(:,:,jk)
230      END DO
231
232      ! conv value on T-point
233      DO jk = 1, jpkm1
234         DO jj = 1, jpj
235            DO ji = 1, jpi
236               zcoef = 0.5_wp / e3t_n(ji,jj,jk)
237               pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
238            END DO
239         END DO
240      END DO
241      !
242      CALL wrk_dealloc( jpi,jpj,jpk, zconv )     
243      !
244   END SUBROUTINE ken_p2k
245
246
247   SUBROUTINE trd_ken_init
248      !!---------------------------------------------------------------------
249      !!                  ***  ROUTINE trd_ken_init  ***
250      !!
251      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
252      !!----------------------------------------------------------------------
253      INTEGER  ::   ji, jj, jk   ! dummy loop indices
254      !!----------------------------------------------------------------------
255      !
256      IF(lwp) THEN
257         WRITE(numout,*)
258         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends'
259         WRITE(numout,*) '~~~~~~~~~~~~~'
260      ENDIF
261      !                           ! allocate box volume arrays
262      IF( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
263      !
264   END SUBROUTINE trd_ken_init
265
266   !!======================================================================
267END MODULE trdken
Note: See TracBrowser for help on using the repository browser.