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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90 @ 3326

Last change on this file since 3326 was 3326, checked in by gm, 12 years ago

Ediag branche: #927 add Potential Energy trend diagnostics (trdpen.F90)

  • Property svn:keywords set to Id
File size: 10.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   USE zdfbfr         ! bottom friction
18   USE ldftra_oce     ! ocean active tracers lateral physics
19   USE sbc_oce        ! surface boundary condition: ocean
20   USE phycst         ! physical constants
21   USE trdvor         ! ocean vorticity trends
22   USE trdglo         ! trends:global domain averaged
23   USE trdmld         ! ocean active mixed layer tracers trends
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O manager library
26   USE lib_mpp        ! MPP library
27   USE wrk_nemo       ! Memory allocation
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   trd_ken       ! called by trddyn module
33   PUBLIC   trd_ken_init  ! called by trdini module
34
35   INTEGER ::   nkstp       ! current time step
36   REAL(wp)::   r1_2_rau0   ! = 0.5 / rau0
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 "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   INTEGER FUNCTION trd_ken_alloc()
52      !!---------------------------------------------------------------------
53      !!                  ***  FUNCTION trd_ken_alloc  ***
54      !!---------------------------------------------------------------------
55      ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
56      !
57      IF( lk_mpp             )   CALL mpp_sum ( trd_ken_alloc )
58      IF( trd_ken_alloc /= 0 )   CALL ctl_warn('trd_ken_alloc: failed to allocate arrays')
59   END FUNCTION trd_ken_alloc
60
61
62   SUBROUTINE trd_ken( putrd, pvtrd, ktrd, kt )
63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE trd_ken  ***
65      !!
66      !! ** Purpose :   output 3D Kinetic Energy trends using IOM
67      !!
68      !! ** Method  : - apply lbc to the input masked velocity trends
69      !!              - compute the associated KE trend:
70      !!          zke = 0.5 * (  mi-1[ un * putrd * bu ] + mj-1[ vn * pvtrd * bv]  ) / bt
71      !!      where bu, bv, bt are the volume of u-, v- and t-boxes.
72      !!              - vertical diffusion case (jpdyn_zdf):
73      !!          diagnose separately the KE trend associated with wind stress
74      !!              - bottom friction case (jpdyn_bfr):
75      !!          explicit case (ln_bfrimp=F): bottom trend put in the 1st level
76      !!                                       of putrd, pvtrd
77      !
78      !
79      !!----------------------------------------------------------------------
80      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V masked trends
81      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
82      INTEGER                   , INTENT(in   ) ::   kt             ! time step
83      !
84      INTEGER ::   ji, jj, jk       ! dummy loop indices
85      INTEGER ::   ikbu  , ikbv     ! local integers
86      INTEGER ::   ikbum1, ikbvm1   !   -       -
87      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, zke2d   ! 2D workspace
88      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace
89      !!----------------------------------------------------------------------
90      !
91      CALL wrk_alloc( jpi, jpj     , z2dx, z2dy, zke2d )
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      IF ( lk_vvl .AND. kt /= nkstp ) THEN   ! Variable volume: set box volume at the 1st call of kt time step
97         nkstp = kt
98         DO jk = 1, jpkm1
99            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk)
100            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk)
101            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) ) * tmask(:,:,jk)
102         END DO
103      ENDIF
104      !
105      zke(:,:,jpk) = 0._wp
106      zke(1,:, : ) = 0._wp
107      zke(:,1, : ) = 0._wp
108      DO jk = 1, jpkm1
109         DO jj = 2, jpj
110            DO ji = 2, jpj
111               zke(ji,jj,jk) = 0.5_wp * (   un(ji  ,jj,jk) * putrd(ji  ,jj,jk) * bu(ji  ,jj,jk)  &
112                  &                       + un(ji-1,jj,jk) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk)  &
113                  &                       + vn(ji,jj  ,jk) * pvtrd(ji,jj  ,jk) * bv(ji,jj  ,jk)  &
114                  &                       + vn(ji,jj-1,jk) * pvtrd(ji,jj-1,jk) * bv(ji,jj-1,jk)  ) * r1_bt(ji,jj,jk)
115            END DO
116         END DO
117      END DO
118      !
119      SELECT CASE( ktrd )
120      CASE( jpdyn_hpg )   ;   CALL iom_put( "ketrd_hpg", zke )    ! hydrostatic pressure gradient
121      CASE( jpdyn_spg )   ;   CALL iom_put( "ketrd_spg", zke )    ! surface pressure gradient
122      CASE( jpdyn_pvo )   ;   CALL iom_put( "ketrd_pvo", zke )    ! planetary vorticity
123      CASE( jpdyn_rvo )   ;   CALL iom_put( "ketrd_rvo", zke )    ! relative  vorticity     (or metric term)
124      CASE( jpdyn_keg )   ;   CALL iom_put( "ketrd_keg", zke )    ! Kinetic Energy gradient (or had)
125      CASE( jpdyn_zad )   ;   CALL iom_put( "ketrd_zad", zke )    ! vertical   advection
126      CASE( jpdyn_ldf )   ;   CALL iom_put( "ketrd_ldf", zke )    ! lateral diffusion
127      CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf", zke )    ! vertical diffusion
128                              !                                   ! wind stress trends
129         z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1u(:,:) * e2u(:,:) * umask(:,:,1)
130         z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1v(:,:) * e2v(:,:) * vmask(:,:,1)
131         zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
132         DO jj = 2, jpj
133            DO ji = 2, jpj
134               zke2d(ji,jj) = r1_2_rau0 * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
135                  &                         + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,1)
136            END DO            !
137         END DO               !
138                              CALL iom_put( "ketrd_tau", zke2d )
139      CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr", zke )    ! bottom friction (explicit)
140!!gm TO BE DONE properly
141!         IF( .NOT.ln_bfrimp ) THEN
142!            DO jj = 1, jpj    !   
143!               DO ji = 1, jpi
144!                  ikbu = mbku(ji,jj)         ! deepest ocean u- & v-levels
145!                  ikbv = mbkv(ji,jj)   
146!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu)
147!                  z2dy(ji,jj) = vn(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv)
148!               END DO
149!            END DO
150!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
151!            DO jj = 2, jpj
152!               DO ji = 2, jpj
153!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
154!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   ) * r1_bt(ji,jj,  BEURK!!!
155!               END DO
156!            END DO
157!                              CALL iom_put( "ketrd_bfr", zke2d )    ! bottom friction (explicit case)
158!!gm only valid if ln_bfrimp=T otherwise the bottom stress as to be recomputed at the end of the compuation....
159!         ENDIF
160!!gm end
161      CASE( jpdyn_atf )   ;   CALL iom_put( "ketrd_atf", zke )    ! asselin filter trends
162
163
164!! a faire !!!!  idee changer dynnxt pour avoir un appel a jpdynbfr avant le swap !!!
165!! reflechir à une possible sauvegarde du "vrai" un,vn pour le calcul de atf....
166!
167!         IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case)
168!            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
169!               DO ji = 1, jpi
170!                  ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
171!                  ikbv = mbkv(ji,jj)
172!                  z2dx(ji,jj) = un(ji,jj,ikbu) * bfrua(ji,jj) * un(ji,jj,ikbu) / fse3u(ji,jj,ikbu)
173!                  z2dy(ji,jj) = un(ji,jj,ikbu) * bfrva(ji,jj) * vn(ji,jj,ikbv) / fse3v(ji,jj,ikbv)
174!               END DO
175!            END DO
176!            zke2d(1,:) = 0._wp   ;   zke2d(:,1) = 0._wp
177!            DO jj = 2, jpj
178!               DO ji = 2, jpj
179!                  zke2d(ji,jj) = 0.5_wp * (   z2dx(ji,jj) + z2dx(ji-1,jj)   &
180!                     &                      + z2dy(ji,jj) + z2dy(ji,jj-1)   )
181!               END DO
182!            END DO
183!                              CALL iom_put( "ketrd_bfr", zke2d )
184!         ENDIF
185         !
186      END SELECT
187      !
188      CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, zke2d )
189      CALL wrk_dealloc( jpi, jpj, jpk, zke )
190      !
191   END SUBROUTINE trd_ken
192
193
194   SUBROUTINE trd_ken_init
195      !!---------------------------------------------------------------------
196      !!                  ***  ROUTINE trd_ken_init  ***
197      !!
198      !! ** Purpose :   initialisation of 3D Kinetic Energy trend diagnostic
199      !!----------------------------------------------------------------------
200      INTEGER  ::   ji, jj, jk   ! dummy loop indices
201      !!----------------------------------------------------------------------
202      !
203      IF(lwp) THEN
204         WRITE(numout,*)
205         WRITE(numout,*) 'trd_ken_init : 3D Kinetic Energy trends'
206         WRITE(numout,*) '~~~~~~~~~~~~~'
207      ENDIF
208      !                           ! allocate box volume arrays
209      IF ( trd_ken_alloc() /= 0 )   CALL ctl_stop('trd_ken_alloc: failed to allocate arrays')
210      !
211      IF ( .NOT.lk_vvl ) THEN     ! constant volume: bu, bv, 1/bt computed one for all
212         DO jk = 1, jpkm1
213            bu   (:,:,jk) =  e1u(:,:) * e2u(:,:) * fse3u_n(:,:,jk)
214            bv   (:,:,jk) =  e1v(:,:) * e2v(:,:) * fse3v_n(:,:,jk)
215            r1_bt(:,:,jk) = 1._wp / ( e1e2t(:,:) * fse3t_n(:,:,jk) )
216         END DO
217      ENDIF
218      !
219      nkstp     = 0
220      r1_2_rau0 = 0.5_wp / rau0
221      !
222   END SUBROUTINE trd_ken_init
223
224   !!======================================================================
225END MODULE trdken
Note: See TracBrowser for help on using the repository browser.