source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRD/trdglo.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp → rn_atfp (use namelist parameter everywhere)
rdtbt → rDt_e
nn_baro → nn_e
rn_scal_load → rn_load
rau0 → rho0

  • Property svn:keywords set to Id
File size: 28.1 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   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt, Kmm )
62      !!---------------------------------------------------------------------
63      !!                  ***  ROUTINE trd_glo  ***
64      !!
65      !! ** Purpose :   compute and print global domain averaged trends for
66      !!              T, T^2, momentum, KE, and KE<->PE
67      !!
68      !!----------------------------------------------------------------------
69      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
70      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
71      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
72      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
73      INTEGER                   , INTENT(in   ) ::   kt      ! time step
74      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
75      !!
76      INTEGER ::   ji, jj, jk      ! dummy loop indices
77      INTEGER ::   ikbu, ikbv      ! local integers
78      REAL(wp)::   zvm, zvt, zvs, z1_2rho0   ! local scalars
79      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
80      !!----------------------------------------------------------------------
81      !
82      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
83         !
84         SELECT CASE( ctype )
85         !
86         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
87            DO_3D_11_11( 1, jpkm1 )
88               zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
89               zvt = ptrdx(ji,jj,jk) * zvm
90               zvs = ptrdy(ji,jj,jk) * zvm
91               tmo(ktrd) = tmo(ktrd) + zvt   
92               smo(ktrd) = smo(ktrd) + zvs
93               t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
94               s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
95            END_3D
96            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
97            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
98               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
99               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
100               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
101               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
102            ENDIF
103            !
104            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
105               !
106               CALL glo_tra_wri( kt )             ! print the results in ocean.output
107               !               
108               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
109               smo(:) = 0._wp
110               t2 (:) = 0._wp
111               s2 (:) = 0._wp
112               !
113            ENDIF
114            !
115         CASE( 'DYN' )          !==  Momentum and KE  ==!       
116            DO_3D_10_10( 1, jpkm1 )
117               zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
118                  &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
119               zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
120                  &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
121               umo(ktrd) = umo(ktrd) + zvt
122               vmo(ktrd) = vmo(ktrd) + zvs
123               hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
124            END_3D
125            !                 
126            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
127               z1_2rho0 = 0.5_wp / rho0
128               DO_2D_10_10
129                  zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
130                     &                                                     * z1_2rho0       * e1e2u(ji,jj)
131                  zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
132                     &                                                     * z1_2rho0       * e1e2v(ji,jj)
133                  umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
134                  vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
135                  hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
136               END_2D
137            ENDIF
138            !                       
139!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
140!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
141!               !
142!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
143!                  z1_2rho0 = 0.5_wp / rho0
144!                  DO jj = 1, jpjm1
145!                     DO ji = 1, jpim1
146!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
147!                        ikbv = mbkv(ji,jj)
148!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
149!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
150!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
151!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
152!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
153!                     END DO
154!                  END DO
155!               ENDIF
156!
157!!gm top drag case is missing
158!
159!               !
160!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
161!               !               
162!               umo(:) = 0._wp                         ! reset for the next time step
163!               vmo(:) = 0._wp
164!               hke(:) = 0._wp
165!               !
166!            ENDIF
167!!gm end
168            !
169         END SELECT
170         !
171      ENDIF
172      !
173   END SUBROUTINE trd_glo
174
175
176   SUBROUTINE glo_dyn_wri( kt, Kmm )
177      !!---------------------------------------------------------------------
178      !!                  ***  ROUTINE glo_dyn_wri  ***
179      !!
180      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
181      !!----------------------------------------------------------------------
182      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
183      INTEGER, INTENT(in) ::   Kmm  ! time level index
184      !
185      INTEGER  ::   ji, jj, jk   ! dummy loop indices
186      REAL(wp) ::   zcof         ! local scalar
187      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
188      !!----------------------------------------------------------------------
189
190      ! I. Momentum trends
191      ! -------------------
192
193      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
194
195         ! I.1 Conversion potential energy - kinetic energy
196         ! --------------------------------------------------
197         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
198         zkx  (:,:,:) = 0._wp
199         zky  (:,:,:) = 0._wp
200         zkz  (:,:,:) = 0._wp
201         zkepe(:,:,:) = 0._wp
202   
203         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
204
205         zcof = 0.5_wp / rho0             ! Density flux at w-point
206         zkz(:,:,1) = 0._wp
207         DO jk = 2, jpk
208            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
209         END DO
210         
211         zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
212         DO_3D_10_10( 1, jpkm1 )
213            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
214            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
215         END_3D
216         
217         DO_3D_00_00( 1, jpkm1 )
218            zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
219               &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
220               &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
221               &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
222         END_3D
223
224         ! I.2 Basin averaged kinetic energy trend
225         ! ----------------------------------------
226         peke = 0._wp
227         DO jk = 1, jpkm1
228            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) * e3t(:,:,jk,Kmm) )
229         END DO
230         peke = grav * peke
231
232         ! I.3 Sums over the global domain
233         ! ---------------------------------
234         IF( lk_mpp ) THEN
235            CALL mpp_sum( 'trdglo', peke )
236            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
237            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
238            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
239         ENDIF
240
241         ! I.2 Print dynamic trends in the ocean.output file
242         ! --------------------------------------------------
243
244         IF(lwp) THEN
245            WRITE (numout,*)
246            WRITE (numout,*)
247            WRITE (numout,9500) kt
248            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
249            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
250            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
251            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
252            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
253            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
254            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
255            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
256            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
257            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
258            WRITE (numout,9511)
259            WRITE (numout,9512)                                                 &
260            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
261            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
262            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
263            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
264            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
265            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
266            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
267!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
268         ENDIF
269
270 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
271 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
272 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
273 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
274 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
275 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
276 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
277 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
278 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
279 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
280 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
281 9511    FORMAT(' -----------------------------------------------------------------------------')
282 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
283 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
284 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
285
286         IF(lwp) THEN
287            WRITE (numout,*)
288            WRITE (numout,*)
289            WRITE (numout,9520) kt
290            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
291            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
292            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
293            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
294            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
295            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
296            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
297            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
298            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
299            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
300            WRITE (numout,9531)
301            WRITE (numout,9532)   &
302            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
303            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
304            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
305            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
306!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
307         ENDIF
308
309 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
310 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
311 9522    FORMAT(' ke gradient               u2= ', e20.13)
312 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
313 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
314 9525    FORMAT(' vertical advection        u2= ', e20.13)
315 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
316 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
317 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
318 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
319 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
320 9531    FORMAT(' --------------------------------------------------')
321 9532    FORMAT(' total trend               u2= ', e20.13)
322 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
323 9534    FORMAT('       bottom stress       u2= ', e20.13)
324
325         IF(lwp) THEN
326            WRITE (numout,*)
327            WRITE (numout,*)
328            WRITE (numout,9540) kt
329            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
330            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
331            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
332            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
333            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
334            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
335            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
336            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
337            WRITE (numout,*)
338            WRITE (numout,*)
339         ENDIF
340
341 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
342 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
343 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
344 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
345 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
346 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
347 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
348 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
349 9548    FORMAT(' pressure gradient u2 = - 1/rho0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
350         !
351         ! Save potential to kinetic energy conversion for next time step
352         rpktrd = peke
353         !
354      ENDIF
355      !
356   END SUBROUTINE glo_dyn_wri
357
358
359   SUBROUTINE glo_tra_wri( kt )
360      !!---------------------------------------------------------------------
361      !!                  ***  ROUTINE glo_tra_wri  ***
362      !!
363      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
364      !!----------------------------------------------------------------------
365      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
366      !
367      INTEGER  ::   jk   ! loop indices
368      !!----------------------------------------------------------------------
369
370      ! I. Tracers trends
371      ! -----------------
372
373      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
374
375         ! I.1 Sums over the global domain
376         ! -------------------------------
377         IF( lk_mpp ) THEN
378            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
379            CALL mpp_sum( 'trdglo', smo, jptot_tra )
380            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
381            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
382         ENDIF
383
384         ! I.2 Print tracers trends in the ocean.output file
385         ! --------------------------------------------------
386         
387         IF(lwp) THEN
388            WRITE (numout,*)
389            WRITE (numout,*)
390            WRITE (numout,9400) kt
391            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
392            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
393            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
394            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
395            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
396            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
397            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
398            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
399            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
400            WRITE (numout,9409) 
401            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
402            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
403            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
404            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
405         ENDIF
406
4079400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
408              '              salinity',/' ============================')
4099401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4109411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4119402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4129403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4139404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4149405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4159406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4169407     FORMAT(' penetrative qsr             ',e20.13)
4179408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4189409     FORMAT(' -------------------------------------------------------------------------')
4199410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
420
421
422         IF(lwp) THEN
423            WRITE (numout,*)
424            WRITE (numout,*)
425            WRITE (numout,9420) kt
426            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
427            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
428            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
429            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
430            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
431            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
432            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
433            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
434            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
435            WRITE (numout,9429)
436            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
437            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
438            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
439            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
440         ENDIF
441
4429420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
443            '               salinity', /, ' ===============================')
4449421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4459431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4469422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4479423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4489424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4499425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4509426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4519427     FORMAT(' penetrative qsr           * t   ', e20.13)
4529428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4539429     FORMAT(' -----------------------------------------------------------------------------')
4549430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
455
456
457         IF(lwp) THEN
458            WRITE (numout,*)
459            WRITE (numout,*)
460            WRITE (numout,9440) kt
461            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
462            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
463            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
464            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
465            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
466            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
467            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
468            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
469            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
470            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
471            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
472         ENDIF
473
4749440     FORMAT(' tracer consistency at it= ',i6,   &
475            ' :         temperature','                salinity',/,   &
476            ' ==================================')
4779441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4789442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4799443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4809444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4819445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4829446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4839447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
4849448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
4859449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
486         !
487      ENDIF
488      !
489   END SUBROUTINE glo_tra_wri
490
491
492   SUBROUTINE trd_glo_init( Kmm )
493      !!---------------------------------------------------------------------
494      !!                  ***  ROUTINE trd_glo_init  ***
495      !!
496      !! ** Purpose :   Read the namtrd namelist
497      !!----------------------------------------------------------------------
498      INTEGER, INTENT(in) ::   Kmm   ! time level index
499      INTEGER  ::   ji, jj, jk   ! dummy loop indices
500      !!----------------------------------------------------------------------
501
502      IF(lwp) THEN
503         WRITE(numout,*)
504         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
505         WRITE(numout,*) '~~~~~~~~~~~~~'
506      ENDIF
507
508      ! Total volume at t-points:
509      tvolt = 0._wp
510      DO jk = 1, jpkm1
511         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
512      END DO
513      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
514
515      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
516
517      ! Initialization of potential to kinetic energy conversion
518      rpktrd = 0._wp
519
520      ! Total volume at u-, v- points:
521!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
522      tvolu = 0._wp
523      tvolv = 0._wp
524
525      DO_3D_00_00( 1, jpk )
526         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
527         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
528      END_3D
529      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
530      CALL mpp_sum( 'trdglo', tvolv )
531
532      IF(lwp) THEN
533         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
534         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
535      ENDIF
536      !
537   END SUBROUTINE trd_glo_init
538
539   !!======================================================================
540END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.