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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch
  • Property svn:keywords set to Id
File size: 28.5 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   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 "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
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      !!
75      INTEGER ::   ji, jj, jk      ! dummy loop indices
76      INTEGER ::   ikbu, ikbv      ! local integers
77      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
78      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
79      !!----------------------------------------------------------------------
80      !
81      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
82         !
83         SELECT CASE( ctype )
84         !
85         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
86            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
87               DO jj = 1, jpj
88                  DO ji = 1, jpi       
89                     zvm = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * 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 * tsn(ji,jj,jk,jp_tem)
95                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
96                  END DO
97               END DO
98            END DO
99            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
100            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
101               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:) )
102               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
103               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:)  )
104               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
105            ENDIF
106            !
107            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
108               !
109               CALL glo_tra_wri( kt )             ! print the results in ocean.output
110               !               
111               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
112               smo(:) = 0._wp
113               t2 (:) = 0._wp
114               s2 (:) = 0._wp
115               !
116            ENDIF
117            !
118         CASE( 'DYN' )          !==  Momentum and KE  ==!       
119            DO jk = 1, jpkm1
120               DO jj = 1, jpjm1
121                  DO ji = 1, jpim1
122                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
123                        &                                     * e1e2u  (ji,jj) * e3u_n(ji,jj,jk)
124                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
125                        &                                     * e1e2v  (ji,jj) * e3u_n(ji,jj,jk)
126                     umo(ktrd) = umo(ktrd) + zvt
127                     vmo(ktrd) = vmo(ktrd) + zvs
128                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
129                  END DO
130               END DO
131            END DO
132            !                 
133            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
134               z1_2rau0 = 0.5_wp / rau0
135               DO jj = 1, jpjm1
136                  DO ji = 1, jpim1
137                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
138                        &                                                     * z1_2rau0       * e1e2u(ji,jj)
139                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
140                        &                                                     * z1_2rau0       * e1e2v(ji,jj)
141                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
142                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
143                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
144                  END DO
145               END DO
146            ENDIF
147            !                       
148            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
149               !
150               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
151                  z1_2rau0 = 0.5_wp / rau0
152                  DO jj = 1, jpjm1
153                     DO ji = 1, jpim1
154                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
155                        ikbv = mbkv(ji,jj)
156                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * un(ji,jj,ikbu) * e1e2u(ji,jj)
157                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vn(ji,jj,ikbv) * e1e2v(ji,jj)
158                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
159                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
160                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
161                     END DO
162                  END DO
163               ENDIF
164!!gm top drag case is missing
165               !
166               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
167               !               
168               umo(:) = 0._wp                         ! reset for the next time step
169               vmo(:) = 0._wp
170               hke(:) = 0._wp
171               !
172            ENDIF
173            !
174         END SELECT
175         !
176      ENDIF
177      !
178   END SUBROUTINE trd_glo
179
180
181   SUBROUTINE glo_dyn_wri( kt )
182      !!---------------------------------------------------------------------
183      !!                  ***  ROUTINE glo_dyn_wri  ***
184      !!
185      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
186      !!----------------------------------------------------------------------
187      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
188      !
189      INTEGER  ::   ji, jj, jk   ! dummy loop indices
190      REAL(wp) ::   zcof         ! local scalar
191      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
192      !!----------------------------------------------------------------------
193
194      ! I. Momentum trends
195      ! -------------------
196
197      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
198
199         ! I.1 Conversion potential energy - kinetic energy
200         ! --------------------------------------------------
201         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
202         zkx  (:,:,:) = 0._wp
203         zky  (:,:,:) = 0._wp
204         zkz  (:,:,:) = 0._wp
205         zkepe(:,:,:) = 0._wp
206   
207         CALL eos( tsn, rhd, rhop )       ! now potential density
208
209         zcof = 0.5_wp / rau0             ! Density flux at w-point
210         zkz(:,:,1) = 0._wp
211         DO jk = 2, jpk
212            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
213         END DO
214         
215         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
216         DO jk = 1, jpkm1
217            DO jj = 1, jpjm1
218               DO ji = 1, jpim1
219                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
220                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
221               END DO
222            END DO
223         END DO
224         
225         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
226            DO jj = 2, jpjm1
227               DO ji = 2, jpim1
228                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
229                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
230                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
231                     &              / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
232               END DO
233            END DO
234         END DO
235
236         ! I.2 Basin averaged kinetic energy trend
237         ! ----------------------------------------
238         peke = 0._wp
239         DO jk = 1, jpkm1
240            peke = peke + SUM( zkepe(:,:,jk) * gdept_n(:,:,jk) * e1e2t(:,:) * e3t_n(:,:,jk) )
241         END DO
242         peke = grav * peke
243
244         ! I.3 Sums over the global domain
245         ! ---------------------------------
246         IF( lk_mpp ) THEN
247            CALL mpp_sum( peke )
248            CALL mpp_sum( umo , jptot_dyn )
249            CALL mpp_sum( vmo , jptot_dyn )
250            CALL mpp_sum( hke , jptot_dyn )
251         ENDIF
252
253         ! I.2 Print dynamic trends in the ocean.output file
254         ! --------------------------------------------------
255
256         IF(lwp) THEN
257            WRITE (numout,*)
258            WRITE (numout,*)
259            WRITE (numout,9500) kt
260            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
261            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
262            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
263            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
264            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
265            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
266            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
267            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
268            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
269            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
270            WRITE (numout,9511)
271            WRITE (numout,9512)                                                 &
272            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
273            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
274            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
275            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
276            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
277            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
278            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
279            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
280         ENDIF
281
282 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
283 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
284 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
285 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
286 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
287 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
288 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
289 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
290 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
291 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
292 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
293 9511    FORMAT(' -----------------------------------------------------------------------------')
294 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
295 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
296 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
297
298         IF(lwp) THEN
299            WRITE (numout,*)
300            WRITE (numout,*)
301            WRITE (numout,9520) kt
302            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
303            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
304            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
305            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
306            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
307            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
308            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
309            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
310            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
311            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
312            WRITE (numout,9531)
313            WRITE (numout,9532)   &
314            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
315            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
316            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
317            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
318            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
319         ENDIF
320
321 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
322 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
323 9522    FORMAT(' ke gradient               u2= ', e20.13)
324 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
325 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
326 9525    FORMAT(' vertical advection        u2= ', e20.13)
327 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
328 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
329 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
330 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
331 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
332 9531    FORMAT(' --------------------------------------------------')
333 9532    FORMAT(' total trend               u2= ', e20.13)
334 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
335 9534    FORMAT('       bottom stress       u2= ', e20.13)
336
337         IF(lwp) THEN
338            WRITE (numout,*)
339            WRITE (numout,*)
340            WRITE (numout,9540) kt
341            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
342            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
343            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
344            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
345            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
346            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
347            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
348            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
349            WRITE (numout,*)
350            WRITE (numout,*)
351         ENDIF
352
353 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
354 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
355 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
356 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
357 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
358 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
359 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
360 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
361 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
362         !
363         ! Save potential to kinetic energy conversion for next time step
364         rpktrd = peke
365         !
366      ENDIF
367      !
368   END SUBROUTINE glo_dyn_wri
369
370
371   SUBROUTINE glo_tra_wri( kt )
372      !!---------------------------------------------------------------------
373      !!                  ***  ROUTINE glo_tra_wri  ***
374      !!
375      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
376      !!----------------------------------------------------------------------
377      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
378      !
379      INTEGER  ::   jk   ! loop indices
380      !!----------------------------------------------------------------------
381
382      ! I. Tracers trends
383      ! -----------------
384
385      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
386
387         ! I.1 Sums over the global domain
388         ! -------------------------------
389         IF( lk_mpp ) THEN
390            CALL mpp_sum( tmo, jptot_tra )   
391            CALL mpp_sum( smo, jptot_tra )
392            CALL mpp_sum( t2 , jptot_tra )
393            CALL mpp_sum( s2 , jptot_tra )
394         ENDIF
395
396         ! I.2 Print tracers trends in the ocean.output file
397         ! --------------------------------------------------
398         
399         IF(lwp) THEN
400            WRITE (numout,*)
401            WRITE (numout,*)
402            WRITE (numout,9400) kt
403            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
404            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
405            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
406            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
407            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
408            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
409            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
410            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
411            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
412            WRITE (numout,9409) 
413            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
414            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
415            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
416            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
417         ENDIF
418
4199400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
420              '              salinity',/' ============================')
4219401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4229411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4239402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4249403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4259404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4269405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4279406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4289407     FORMAT(' penetrative qsr             ',e20.13)
4299408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4309409     FORMAT(' -------------------------------------------------------------------------')
4319410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
432
433
434         IF(lwp) THEN
435            WRITE (numout,*)
436            WRITE (numout,*)
437            WRITE (numout,9420) kt
438            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
439            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
440            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
441            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
442            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
443            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
444            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
445            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
446            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
447            WRITE (numout,9429)
448            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
449            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
450            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
451            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
452         ENDIF
453
4549420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
455            '               salinity', /, ' ===============================')
4569421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4579431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4589422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4599423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4609424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4619425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4629426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4639427     FORMAT(' penetrative qsr           * t   ', e20.13)
4649428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4659429     FORMAT(' -----------------------------------------------------------------------------')
4669430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
467
468
469         IF(lwp) THEN
470            WRITE (numout,*)
471            WRITE (numout,*)
472            WRITE (numout,9440) kt
473            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
474            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
475            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
476            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
477            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
478            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
479            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
480            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
481            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
482            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
483            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
484         ENDIF
485
4869440     FORMAT(' tracer consistency at it= ',i6,   &
487            ' :         temperature','                salinity',/,   &
488            ' ==================================')
4899441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4909442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4919443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4929444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4939445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4949446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4959447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
4969448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
4979449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
498         !
499      ENDIF
500      !
501   END SUBROUTINE glo_tra_wri
502
503
504   SUBROUTINE trd_glo_init
505      !!---------------------------------------------------------------------
506      !!                  ***  ROUTINE trd_glo_init  ***
507      !!
508      !! ** Purpose :   Read the namtrd namelist
509      !!----------------------------------------------------------------------
510      INTEGER  ::   ji, jj, jk   ! dummy loop indices
511      !!----------------------------------------------------------------------
512
513      IF(lwp) THEN
514         WRITE(numout,*)
515         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
516         WRITE(numout,*) '~~~~~~~~~~~~~'
517      ENDIF
518
519      ! Total volume at t-points:
520      tvolt = 0._wp
521      DO jk = 1, jpkm1
522         tvolt = tvolt + SUM( e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
523      END DO
524      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
525
526      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
527
528      ! Initialization of potential to kinetic energy conversion
529      rpktrd = 0._wp
530
531      ! Total volume at u-, v- points:
532!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
533      tvolu = 0._wp
534      tvolv = 0._wp
535
536      DO jk = 1, jpk
537         DO jj = 2, jpjm1
538            DO ji = fs_2, fs_jpim1   ! vector opt.
539               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u_n(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
540               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v_n(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
541            END DO
542         END DO
543      END DO
544      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
545      IF( lk_mpp )   CALL mpp_sum( tvolv )
546
547      IF(lwp) THEN
548         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
549         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
550      ENDIF
551      !
552   END SUBROUTINE trd_glo_init
553
554   !!======================================================================
555END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.