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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 3325

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

Ediag branche: #927 add Kinetic Energy trend diagnostics (trdken.F90)
  • 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_oce      ! ocean active tracers: lateral physics
22   USE ldfdyn_oce      ! 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   USE lib_mpp         ! distibuted memory computing library
29   USE in_out_manager  ! I/O manager
30   USE iom             ! I/O manager library
31   USE wrk_nemo        ! Memory allocation
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 "domzgr_substitute.h90"
54#  include "vectopt_loop_substitute.h90"
55#  include "zdfddm_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE trd_glo  ***
66      !!
67      !! ** Purpose :   compute and print global domain averaged trends for
68      !!              T, T^2, momentum, KE, and KE<->PE
69      !!
70      !!----------------------------------------------------------------------
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
73      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
74      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
75      INTEGER                   , INTENT(in   ) ::   kt      ! time step
76      !!
77      INTEGER ::   ji, jj, jk      ! dummy loop indices
78      INTEGER ::   ikbu, ikbv      ! local integers
79      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
80      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
81      !!----------------------------------------------------------------------
82
83      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
84
85      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
86         !
87         SELECT CASE( ctype )
88         !
89         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
90            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
91               DO jj = 1, jpj
92                  DO ji = 1, jpi       
93                     zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
94                     zvt = ptrdx(ji,jj,jk) * zvm
95                     zvs = ptrdy(ji,jj,jk) * zvm
96                     tmo(ktrd) = tmo(ktrd) + zvt   
97                     smo(ktrd) = smo(ktrd) + zvs
98                     t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
99                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
100                  END DO
101               END DO
102            END DO
103            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
104            IF( .NOT.lk_vvl .AND. ktrd == jptra_zad ) THEN 
105               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
106               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
107               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
108               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
109            ENDIF
110            !
111            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
112               !
113               CALL glo_tra_wri( kt )             ! print the results in ocean.output
114               !               
115               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
116               smo(:) = 0._wp
117               t2 (:) = 0._wp
118               s2 (:) = 0._wp
119               !
120            ENDIF
121            !
122         CASE( 'DYN' )          !==  Momentum and KE  ==!       
123            DO jk = 1, jpkm1
124               DO jj = 1, jpjm1
125                  DO ji = 1, jpim1
126                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
127                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
128                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
129                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
130                     umo(ktrd) = umo(ktrd) + zvt
131                     vmo(ktrd) = vmo(ktrd) + zvs
132                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
133                  END DO
134               END DO
135            END DO
136            !                 
137            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
138               z1_2rau0 = 0.5_wp / rau0
139               DO jj = 1, jpjm1
140                  DO ji = 1, jpim1
141                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
142                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
143                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
144                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
145                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
146                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
147                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
148                  END DO
149               END DO
150            ENDIF
151            !                       
152            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
153               !
154               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction
155                  z1_2rau0 = 0.5_wp / rau0
156                  DO jj = 1, jpjm1
157                     DO ji = 1, jpim1
158                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
159                        ikbv = mbkv(ji,jj)
160                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
161                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
162                        umo(jpdyn_bfr) = umo(jpdyn_bfr) + zvt
163                        vmo(jpdyn_bfr) = vmo(jpdyn_bfr) + zvs
164                        hke(jpdyn_bfr) = hke(jpdyn_bfr) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
165                     END DO
166                  END DO
167               ENDIF
168               !
169               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
170               !               
171               umo(:) = 0._wp                         ! reset for the next time step
172               vmo(:) = 0._wp
173               hke(:) = 0._wp
174               !
175            ENDIF
176            !
177         END SELECT
178         !
179      ENDIF
180      !
181      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
182      !
183   END SUBROUTINE trd_glo
184
185
186   SUBROUTINE glo_dyn_wri( kt )
187      !!---------------------------------------------------------------------
188      !!                  ***  ROUTINE glo_dyn_wri  ***
189      !!
190      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
191      !!----------------------------------------------------------------------
192      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
193      !
194      INTEGER  ::   ji, jj, jk   ! dummy loop indices
195      REAL(wp) ::   zcof         ! local scalar
196      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
197      !!----------------------------------------------------------------------
198
199      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
200
201      ! I. Momentum trends
202      ! -------------------
203
204      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
205
206         ! I.1 Conversion potential energy - kinetic energy
207         ! --------------------------------------------------
208         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
209         zkx  (:,:,:) = 0._wp
210         zky  (:,:,:) = 0._wp
211         zkz  (:,:,:) = 0._wp
212         zkepe(:,:,:) = 0._wp
213   
214         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
215
216         zcof = 0.5_wp / rau0             ! Density flux at w-point
217         zkz(:,:,1) = 0._wp
218         DO jk = 2, jpk
219            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
220         END DO
221         
222         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
223         DO jk = 1, jpkm1
224            DO jj = 1, jpjm1
225               DO ji = 1, jpim1
226                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
227                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
228               END DO
229            END DO
230         END DO
231         
232         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
233            DO jj = 2, jpjm1
234               DO ji = 2, jpim1
235                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
236                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
237                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
238                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
239               END DO
240            END DO
241         END DO
242
243         ! I.2 Basin averaged kinetic energy trend
244         ! ----------------------------------------
245         peke = 0._wp
246         DO jk = 1, jpkm1
247            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
248         END DO
249         peke = grav * peke
250
251         ! I.3 Sums over the global domain
252         ! ---------------------------------
253         IF( lk_mpp ) THEN
254            CALL mpp_sum( peke )
255            CALL mpp_sum( umo , jptot_dyn )
256            CALL mpp_sum( vmo , jptot_dyn )
257            CALL mpp_sum( hke , jptot_dyn )
258         ENDIF
259
260         ! I.2 Print dynamic trends in the ocean.output file
261         ! --------------------------------------------------
262
263         IF(lwp) THEN
264            WRITE (numout,*)
265            WRITE (numout,*)
266            WRITE (numout,9500) kt
267            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
268            WRITE (numout,9509) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
269            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
270            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
271            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
272            WRITE (numout,9507) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
273            WRITE (numout,9505) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
274            WRITE (numout,9508) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
275            WRITE (numout,9510) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
276            WRITE (numout,9511) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
277            WRITE (numout,9512) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
278            WRITE (numout,9513)
279            WRITE (numout,9514)                                                 &
280            &     (  umo(jpdyn_hpg) + umo(jpdyn_spg) + umo(jpdyn_keg) + umo(jpdyn_rvo)   &
281            &      + umo(jpdyn_pvo) + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf)   &
282            &      + umo(jpdyn_tau) + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
283            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_spg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo)   &
284            &      + vmo(jpdyn_pvo) + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf)   &
285            &      + vmo(jpdyn_tau) + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
286         ENDIF
287
288 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
289 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
290 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
291 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
292 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
293 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
294 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
295 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
296 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
297 9510    FORMAT(' surface wind stress        u= ', e20.13, '    v= ', e20.13)
298 9511    FORMAT(' bottom friction            u= ', e20.13, '    v= ', e20.13)
299 9512    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
300 9513    FORMAT(' -----------------------------------------------------------------------------')
301 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
302
303         IF(lwp) THEN
304            WRITE (numout,*)
305            WRITE (numout,*)
306            WRITE (numout,9520) kt
307            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
308            WRITE (numout,9529) hke(jpdyn_spg) / tvolt
309            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
310            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
311            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
312            WRITE (numout,9527) hke(jpdyn_zad) / tvolt
313            WRITE (numout,9525) hke(jpdyn_ldf) / tvolt
314            WRITE (numout,9528) hke(jpdyn_zdf) / tvolt
315            WRITE (numout,9530) hke(jpdyn_tau) / tvolt
316            WRITE (numout,9531) hke(jpdyn_bfr) / tvolt
317            WRITE (numout,9532) hke(jpdyn_atf) / tvolt
318            WRITE (numout,9533)
319            WRITE (numout,9534)   &
320            &     (  hke(jpdyn_hpg) + hke(jpdyn_spg) + hke(jpdyn_keg) + hke(jpdyn_rvo)   &
321            &      + hke(jpdyn_pvo) + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf)   &
322            &      + hke(jpdyn_tau) + hke(jpdyn_bfr) + hke(jpdyn_atf)  ) / tvolt
323         ENDIF
324
325 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
326 9521    FORMAT(' pressure gradient         u2= ', e20.13)
327 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
328 9522    FORMAT(' ke gradient               u2= ', e20.13)
329 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
330 9524    FORMAT(' coriolis term             u2= ', e20.13)
331 9527    FORMAT(' vertical advection        u2= ', e20.13)
332 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
333 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
334 9530    FORMAT(' surface wind stress       u2= ', e20.13)
335 9531    FORMAT(' bottom friction           u2= ', e20.13)
336 9532    FORMAT(' Asselin time filter       u2= ', e20.13)
337 9533    FORMAT(' --------------------------------------------------')
338 9534    FORMAT(' total trend               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      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
372      !
373   END SUBROUTINE glo_dyn_wri
374
375
376   SUBROUTINE glo_tra_wri( kt )
377      !!---------------------------------------------------------------------
378      !!                  ***  ROUTINE glo_tra_wri  ***
379      !!
380      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
381      !!----------------------------------------------------------------------
382      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
383      !
384      INTEGER  ::   jk   ! loop indices
385      !!----------------------------------------------------------------------
386
387      ! I. Tracers trends
388      ! -----------------
389
390      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
391
392         ! I.1 Sums over the global domain
393         ! -------------------------------
394         IF( lk_mpp ) THEN
395            CALL mpp_sum( tmo, jptot_tra )   
396            CALL mpp_sum( smo, jptot_tra )
397            CALL mpp_sum( t2 , jptot_tra )
398            CALL mpp_sum( s2 , jptot_tra )
399         ENDIF
400
401         ! I.2 Print tracers trends in the ocean.output file
402         ! --------------------------------------------------
403         
404         IF(lwp) THEN
405            WRITE (numout,*)
406            WRITE (numout,*)
407            WRITE (numout,9400) kt
408            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
409            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
410            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
411            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
412            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
413            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
414            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
415            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
416            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
417            WRITE (numout,9409) 
418            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
419            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
420            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
421            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
422         ENDIF
423
4249400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
425              '              salinity',/' ============================')
4269401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4279411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4289402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4299403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4309404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4319405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4329406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4339407     FORMAT(' penetrative qsr             ',e20.13)
4349408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4359409     FORMAT(' -------------------------------------------------------------------------')
4369410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
437
438
439         IF(lwp) THEN
440            WRITE (numout,*)
441            WRITE (numout,*)
442            WRITE (numout,9420) kt
443            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
444            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
445            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
446            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
447            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
448            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
449            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
450            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
451            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
452            WRITE (numout,9429)
453            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
454            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
455            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
456            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
457         ENDIF
458
4599420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
460            '               salinity', /, ' ===============================')
4619421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4629431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4639422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4649423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4659424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4669425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4679426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4689427     FORMAT(' penetrative qsr           * t   ', e20.13)
4699428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4709429     FORMAT(' -----------------------------------------------------------------------------')
4719430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
472
473
474         IF(lwp) THEN
475            WRITE (numout,*)
476            WRITE (numout,*)
477            WRITE (numout,9440) kt
478            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
479            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
480            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
481            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
482            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
483            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
484            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
485            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
486            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
487            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
488            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
489         ENDIF
490
4919440     FORMAT(' tracer consistency at it= ',i6,   &
492            ' :         temperature','                salinity',/,   &
493            ' ==================================')
4949441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4959442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4969443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4979444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4989445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4999446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5009447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5019448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5029449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
503         !
504      ENDIF
505      !
506   END SUBROUTINE glo_tra_wri
507
508
509   SUBROUTINE trd_glo_init
510      !!---------------------------------------------------------------------
511      !!                  ***  ROUTINE trd_glo_init  ***
512      !!
513      !! ** Purpose :   Read the namtrd namelist
514      !!----------------------------------------------------------------------
515      INTEGER  ::   ji, jj, jk   ! dummy loop indices
516      !!----------------------------------------------------------------------
517
518      IF(lwp) THEN
519         WRITE(numout,*)
520         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
521         WRITE(numout,*) '~~~~~~~~~~~~~'
522      ENDIF
523
524      ! Total volume at t-points:
525      tvolt = 0._wp
526      DO jk = 1, jpkm1
527         tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
528      END DO
529      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
530
531      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
532
533      ! Initialization of potential to kinetic energy conversion
534      rpktrd = 0._wp
535
536      ! Total volume at u-, v- points:
537!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
538      tvolu = 0._wp
539      tvolv = 0._wp
540
541      DO jk = 1, jpk
542         DO jj = 2, jpjm1
543            DO ji = fs_2, fs_jpim1   ! vector opt.
544               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
545               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
546            END DO
547         END DO
548      END DO
549      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
550      IF( lk_mpp )   CALL mpp_sum( tvolv )
551
552      IF(lwp) THEN
553         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
554         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
555      ENDIF
556      !
557   END SUBROUTINE trd_glo_init
558
559   !!======================================================================
560END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.