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.
trdglo.F90 in NEMO/trunk/src/OCE/TRD – NEMO

source: NEMO/trunk/src/OCE/TRD/trdglo.F90 @ 13286

Last change on this file since 13286 was 13237, checked in by smasson, 4 years ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3
  • Property svn:keywords set to Id
File size: 28.4 KB
Line 
1MODULE trdglo
2   !!======================================================================
3   !!                       ***  MODULE  trdglo  ***
4   !! Ocean diagnostics:  global domain averaged tracer and momentum trends
5   !!=====================================================================
6   !! History :  1.0  !  2004-08  (C. Talandier) New trends organization
7   !!            3.5  !  2012-02  (G. Madec)  add 3D tracer zdf trend output using iom
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   trd_glo       : domain averaged budget of trends (including kinetic energy and T^2 trends)
12   !!   glo_dyn_wri   : print dynamic trends in ocean.output file
13   !!   glo_tra_wri   : print global T & T^2 trends in ocean.output file
14   !!   trd_glo_init  : initialization step
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE trd_oce        ! trends: ocean variables
20   USE phycst         ! physical constants
21   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
22   USE ldfdyn         ! ocean dynamics: lateral physics
23   USE zdf_oce        ! ocean vertical physics
24!!gm   USE zdfdrg         ! ocean vertical physics: bottom friction
25   USE zdfddm         ! ocean vertical physics: double diffusion
26   USE eosbn2         ! equation of state
27   USE phycst         ! physical constants
28   !
29   USE lib_mpp        ! distibuted memory computing library
30   USE in_out_manager ! I/O manager
31   USE iom            ! I/O manager library
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   trd_glo       ! called by trdtra and trddyn modules
37   PUBLIC   trd_glo_init  ! called by trdini module
38
39   !                     !!! Variables used for diagnostics
40   REAL(wp) ::   tvolt    ! volume of the whole ocean computed at t-points
41   REAL(wp) ::   tvolu    ! volume of the whole ocean computed at u-points
42   REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
43   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
44   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
45
46   !                     !!! domain averaged trends
47   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
48   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
49   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
50   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
51
52   !! * Substitutions
53#  include "do_loop_substitute.h90"
54#  include "domzgr_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
57   !! $Id$
58   !! Software governed by the CeCILL license (see ./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt, Kmm )
63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE trd_glo  ***
65      !!
66      !! ** Purpose :   compute and print global domain averaged trends for
67      !!              T, T^2, momentum, KE, and KE<->PE
68      !!
69      !!----------------------------------------------------------------------
70      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
72      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
73      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
74      INTEGER                   , INTENT(in   ) ::   kt      ! time step
75      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
76      !!
77      INTEGER ::   ji, jj, jk      ! dummy loop indices
78      INTEGER ::   ikbu, ikbv      ! local integers
79      REAL(wp)::   zvm, zvt, zvs, z1_2rho0   ! local scalars
80      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
81      !!----------------------------------------------------------------------
82      !
83      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
84         !
85         SELECT CASE( ctype )
86         !
87         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
88            DO_3D_11_11( 1, jpkm1 )
89               zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
90               zvt = ptrdx(ji,jj,jk) * zvm
91               zvs = ptrdy(ji,jj,jk) * zvm
92               tmo(ktrd) = tmo(ktrd) + zvt   
93               smo(ktrd) = smo(ktrd) + zvs
94               t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
95               s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
96            END_3D
97            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
98            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
99               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
100               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
101               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
102               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
103            ENDIF
104            !
105            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
106               !
107               CALL glo_tra_wri( kt )             ! print the results in ocean.output
108               !               
109               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
110               smo(:) = 0._wp
111               t2 (:) = 0._wp
112               s2 (:) = 0._wp
113               !
114            ENDIF
115            !
116         CASE( 'DYN' )          !==  Momentum and KE  ==!       
117            DO_3D_10_10( 1, jpkm1 )
118               zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
119                  &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
120               zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
121                  &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
122               umo(ktrd) = umo(ktrd) + zvt
123               vmo(ktrd) = vmo(ktrd) + zvs
124               hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
125            END_3D
126            !                 
127            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
128               z1_2rho0 = 0.5_wp / rho0
129               DO_2D_10_10
130                  zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
131                     &                                                     * z1_2rho0       * e1e2u(ji,jj)
132                  zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
133                     &                                                     * z1_2rho0       * e1e2v(ji,jj)
134                  umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
135                  vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
136                  hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
137               END_2D
138            ENDIF
139            !                       
140!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
141!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
142!               !
143!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
144!                  z1_2rho0 = 0.5_wp / rho0
145!                  DO jj = 1, jpjm1
146!                     DO ji = 1, jpim1
147!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
148!                        ikbv = mbkv(ji,jj)
149!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
150!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
151!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
152!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
153!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
154!                     END DO
155!                  END DO
156!               ENDIF
157!
158!!gm top drag case is missing
159!
160!               !
161!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
162!               !               
163!               umo(:) = 0._wp                         ! reset for the next time step
164!               vmo(:) = 0._wp
165!               hke(:) = 0._wp
166!               !
167!            ENDIF
168!!gm end
169            !
170         END SELECT
171         !
172      ENDIF
173      !
174   END SUBROUTINE trd_glo
175
176
177   SUBROUTINE glo_dyn_wri( kt, Kmm )
178      !!---------------------------------------------------------------------
179      !!                  ***  ROUTINE glo_dyn_wri  ***
180      !!
181      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
182      !!----------------------------------------------------------------------
183      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
184      INTEGER, INTENT(in) ::   Kmm  ! time level index
185      !
186      INTEGER  ::   ji, jj, jk   ! dummy loop indices
187      REAL(wp) ::   zcof         ! local scalar
188      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
189      !!----------------------------------------------------------------------
190
191      ! I. Momentum trends
192      ! -------------------
193
194      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
195
196         ! I.1 Conversion potential energy - kinetic energy
197         ! --------------------------------------------------
198         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
199         zkx  (:,:,:) = 0._wp
200         zky  (:,:,:) = 0._wp
201         zkz  (:,:,:) = 0._wp
202         zkepe(:,:,:) = 0._wp
203   
204         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
205
206         zcof = 0.5_wp / rho0             ! Density flux at w-point
207         zkz(:,:,1) = 0._wp
208         DO jk = 2, jpk
209            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
210         END DO
211         
212         zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
213         DO_3D_10_10( 1, jpkm1 )
214            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
215               &                              *  uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
216            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
217               &                              *  vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
218         END_3D
219         
220         DO_3D_00_00( 1, jpkm1 )
221            zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
222               &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
223               &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
224               &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
225         END_3D
226
227         ! I.2 Basin averaged kinetic energy trend
228         ! ----------------------------------------
229         peke = 0._wp
230         DO jk = 1, jpkm1
231            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:)   &
232               &                               * e3t(:,:,jk,Kmm) )
233         END DO
234         peke = grav * peke
235
236         ! I.3 Sums over the global domain
237         ! ---------------------------------
238         IF( lk_mpp ) THEN
239            CALL mpp_sum( 'trdglo', peke )
240            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
241            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
242            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
243         ENDIF
244
245         ! I.2 Print dynamic trends in the ocean.output file
246         ! --------------------------------------------------
247
248         IF(lwp) THEN
249            WRITE (numout,*)
250            WRITE (numout,*)
251            WRITE (numout,9500) kt
252            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
253            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
254            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
255            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
256            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
257            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
258            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
259            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
260            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
261            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
262            WRITE (numout,9511)
263            WRITE (numout,9512)                                                 &
264            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
265            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
266            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
267            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
268            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
269            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
270            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
271!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
272         ENDIF
273
274 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
275 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
276 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
277 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
278 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
279 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
280 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
281 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
282 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
283 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
284 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
285 9511    FORMAT(' -----------------------------------------------------------------------------')
286 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
287 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
288 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
289
290         IF(lwp) THEN
291            WRITE (numout,*)
292            WRITE (numout,*)
293            WRITE (numout,9520) kt
294            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
295            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
296            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
297            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
298            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
299            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
300            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
301            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
302            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
303            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
304            WRITE (numout,9531)
305            WRITE (numout,9532)   &
306            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
307            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
308            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
309            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
310!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
311         ENDIF
312
313 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
314 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
315 9522    FORMAT(' ke gradient               u2= ', e20.13)
316 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
317 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
318 9525    FORMAT(' vertical advection        u2= ', e20.13)
319 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
320 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
321 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
322 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
323 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
324 9531    FORMAT(' --------------------------------------------------')
325 9532    FORMAT(' total trend               u2= ', e20.13)
326 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
327 9534    FORMAT('       bottom stress       u2= ', e20.13)
328
329         IF(lwp) THEN
330            WRITE (numout,*)
331            WRITE (numout,*)
332            WRITE (numout,9540) kt
333            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
334            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
335            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
336            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
337            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
338            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
339            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
340            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
341            WRITE (numout,*)
342            WRITE (numout,*)
343         ENDIF
344
345 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
346 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
347 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
348 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
349 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
350 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
351 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
352 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
353 9548    FORMAT(' pressure gradient u2 = - 1/rho0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
354         !
355         ! Save potential to kinetic energy conversion for next time step
356         rpktrd = peke
357         !
358      ENDIF
359      !
360   END SUBROUTINE glo_dyn_wri
361
362
363   SUBROUTINE glo_tra_wri( kt )
364      !!---------------------------------------------------------------------
365      !!                  ***  ROUTINE glo_tra_wri  ***
366      !!
367      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
368      !!----------------------------------------------------------------------
369      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
370      !
371      INTEGER  ::   jk   ! loop indices
372      !!----------------------------------------------------------------------
373
374      ! I. Tracers trends
375      ! -----------------
376
377      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
378
379         ! I.1 Sums over the global domain
380         ! -------------------------------
381         IF( lk_mpp ) THEN
382            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
383            CALL mpp_sum( 'trdglo', smo, jptot_tra )
384            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
385            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
386         ENDIF
387
388         ! I.2 Print tracers trends in the ocean.output file
389         ! --------------------------------------------------
390         
391         IF(lwp) THEN
392            WRITE (numout,*)
393            WRITE (numout,*)
394            WRITE (numout,9400) kt
395            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
396            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
397            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
398            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
399            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
400            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
401            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
402            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
403            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
404            WRITE (numout,9409) 
405            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
406            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
407            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
408            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
409         ENDIF
410
4119400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
412              '              salinity',/' ============================')
4139401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4149411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4159402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4169403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4179404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4189405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4199406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4209407     FORMAT(' penetrative qsr             ',e20.13)
4219408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4229409     FORMAT(' -------------------------------------------------------------------------')
4239410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
424
425
426         IF(lwp) THEN
427            WRITE (numout,*)
428            WRITE (numout,*)
429            WRITE (numout,9420) kt
430            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
431            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
432            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
433            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
434            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
435            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
436            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
437            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
438            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
439            WRITE (numout,9429)
440            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
441            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
442            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
443            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
444         ENDIF
445
4469420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
447            '               salinity', /, ' ===============================')
4489421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4499431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4509422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4519423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4529424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4539425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4549426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4559427     FORMAT(' penetrative qsr           * t   ', e20.13)
4569428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4579429     FORMAT(' -----------------------------------------------------------------------------')
4589430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
459
460
461         IF(lwp) THEN
462            WRITE (numout,*)
463            WRITE (numout,*)
464            WRITE (numout,9440) kt
465            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
466            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
467            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
468            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
469            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
470            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
471            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
472            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
473            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
474            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
475            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
476         ENDIF
477
4789440     FORMAT(' tracer consistency at it= ',i6,   &
479            ' :         temperature','                salinity',/,   &
480            ' ==================================')
4819441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4829442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4839443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4849444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4859445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4869446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4879447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
4889448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
4899449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
490         !
491      ENDIF
492      !
493   END SUBROUTINE glo_tra_wri
494
495
496   SUBROUTINE trd_glo_init( Kmm )
497      !!---------------------------------------------------------------------
498      !!                  ***  ROUTINE trd_glo_init  ***
499      !!
500      !! ** Purpose :   Read the namtrd namelist
501      !!----------------------------------------------------------------------
502      INTEGER, INTENT(in) ::   Kmm   ! time level index
503      INTEGER  ::   ji, jj, jk   ! dummy loop indices
504      !!----------------------------------------------------------------------
505
506      IF(lwp) THEN
507         WRITE(numout,*)
508         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
509         WRITE(numout,*) '~~~~~~~~~~~~~'
510      ENDIF
511
512      ! Total volume at t-points:
513      tvolt = 0._wp
514      DO jk = 1, jpkm1
515         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
516      END DO
517      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
518
519      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
520
521      ! Initialization of potential to kinetic energy conversion
522      rpktrd = 0._wp
523
524      ! Total volume at u-, v- points:
525!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
526      tvolu = 0._wp
527      tvolv = 0._wp
528
529      DO_3D_00_00( 1, jpk )
530         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
531            &                                       * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
532         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
533            &                                       * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
534      END_3D
535      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
536      CALL mpp_sum( 'trdglo', tvolv )
537
538      IF(lwp) THEN
539         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
540         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
541      ENDIF
542      !
543   END SUBROUTINE trd_glo_init
544
545   !!======================================================================
546END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.