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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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