source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 3317

Last change on this file since 3317 was 3317, checked in by gm, 9 years ago

Ediag branche: #927 restructuration of the trdicp computation - part I

  • Property svn:keywords set to Id
File size: 29.2 KB
Line 
1MODULE trdicp
2   !!======================================================================
3   !!                       ***  MODULE  trdicp  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic 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_budget     : domain averaged budget of trends (including kinetic energy and T^2 trends)
12   !!   trd_icp        : compute the basin averaged properties for tra/dyn
13   !!   trd_dwr        : print dynmaic trends in ocean.output file
14   !!   trd_twr        : print tracers trends in ocean.output file
15   !!   trd_icp_init   : initialization step
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain variables
19   USE sbc_oce         ! surface boundary condition: ocean
20   USE phycst          ! physical constants
21   USE trdmod_oce      ! ocean variables trends
22   USE ldftra_oce      ! ocean active tracers: lateral physics
23   USE ldfdyn_oce      ! ocean dynamics: lateral physics
24   USE zdf_oce         ! ocean vertical physics
25   USE zdfbfr          ! bottom friction
26   USE zdfddm          ! ocean vertical physics: double diffusion
27   USE eosbn2          ! equation of state
28   USE phycst          ! physical constants
29   USE lib_mpp         ! distibuted memory computing library
30   USE in_out_manager  ! I/O manager
31   USE iom             ! I/O manager library
32   USE wrk_nemo        ! Memory allocation
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   trd_budget    ! called by trdmod.F90
38   PUBLIC   trd_dwr       ! called by step.F90
39   PUBLIC   trd_twr       ! called by step.F90
40   PUBLIC   trd_icp_init  ! called by opa.F90
41
42   !                     !!! Variables used for diagnostics
43   REAL(wp) ::   tvolt    ! volume of the whole ocean computed at t-points
44   REAL(wp) ::   tvolu    ! volume of the whole ocean computed at u-points
45   REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
46   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
47   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
48
49   !                     !!! domain averaged trends
50   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
51   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
52   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
53   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58#  include "zdfddm_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE trd_budget( ptrdx, ptrdy, ktrd, ctype, kt )
67      !!---------------------------------------------------------------------
68      !!                  ***  ROUTINE trd_budget  ***
69      !!
70      !! ** Purpose : integral constraint diagnostics for momentum and/or tracer trends
71      !!             
72      !!----------------------------------------------------------------------
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
75      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
76      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type 'DYN'/'TRA'
77      INTEGER                   , INTENT(in   ) ::   kt      ! time step
78      !!
79      INTEGER ::   ji, jj, jk      ! dummy loop indices
80      INTEGER ::   ikbu, ikbv      ! local integers
81      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
82      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
83      !!----------------------------------------------------------------------
84
85      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
86
87      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
88         !
89         SELECT CASE( ctype )
90         !
91         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
92            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
93               DO jj = 1, jpj
94                  DO ji = 1, jpi       
95                     zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
96                     zvt = ptrdx(ji,jj,jk) * zvm
97                     zvs = ptrdy(ji,jj,jk) * zvm
98                     tmo(ktrd) = tmo(ktrd) + zvt   
99                     smo(ktrd) = smo(ktrd) + zvs
100                     t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
101                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
102                  END DO
103               END DO
104            END DO
105            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
106            IF( .NOT.lk_vvl .AND. ktrd == jptra_trd_zad ) THEN 
107               tmo(jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
108               smo(jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
109               t2 (jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
110               s2 (jptra_trd_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
111            ENDIF
112            !
113            IF( ktrd == jptra_trd_atf ) THEN     ! last trend (asselin time filter)
114               !
115               CALL trd_twr( kt )   ! print the results in ocean.output
116               !               
117               tmo(:) = 0._wp       ! prepare the next time step (domain averaged array reset to zero)
118               smo(:) = 0._wp
119               t2 (:) = 0._wp
120               s2 (:) = 0._wp
121               !
122            ENDIF
123            !
124         CASE( 'DYN' )          !==  Momentum and KE  ==!       
125            DO jk = 1, jpkm1
126               DO jj = 1, jpjm1
127                  DO ji = 1, jpim1
128                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
129                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
130                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
131                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
132                     umo(ktrd) = umo(ktrd) + zvt
133                     vmo(ktrd) = vmo(ktrd) + zvs
134                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
135                  END DO
136               END DO
137            END DO
138            !                 
139            IF( ktrd == jpdyn_trd_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
140               z1_2rau0 = 0.5_wp / rau0
141               DO jj = 1, jpjm1
142                  DO ji = 1, jpim1
143                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
144                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
145                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
146                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
147                     umo(jpdyn_trd_tau) = umo(jpdyn_trd_tau) + zvt
148                     vmo(jpdyn_trd_tau) = vmo(jpdyn_trd_tau) + zvs
149                     hke(jpdyn_trd_tau) = hke(jpdyn_trd_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
150                  END DO
151               END DO
152            ENDIF
153            !                       
154            IF( ktrd == jpdyn_trd_atf ) THEN     ! last trend (asselin time filter)
155               !
156               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction
157                  z1_2rau0 = 0.5_wp / rau0
158                  DO jj = 1, jpjm1
159                     DO ji = 1, jpim1
160                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
161                        ikbv = mbkv(ji,jj)
162                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
163                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
164                        umo(jpdyn_trd_bfr) = umo(jpdyn_trd_bfr) + zvt
165                        vmo(jpdyn_trd_bfr) = vmo(jpdyn_trd_bfr) + zvs
166                        hke(jpdyn_trd_bfr) = hke(jpdyn_trd_bfr) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
167                     END DO
168                  END DO
169               ENDIF
170               !
171               CALL trd_dwr( kt )                     ! print the results in ocean.output
172               !               
173               umo(:) = 0._wp                         ! reset for the next time step
174               vmo(:) = 0._wp
175               hke(:) = 0._wp
176               !
177            ENDIF
178            !
179         END SELECT
180         !
181      ENDIF
182      !
183      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
184      !
185   END SUBROUTINE trd_budget
186
187
188   SUBROUTINE trd_icp_init
189      !!---------------------------------------------------------------------
190      !!                  ***  ROUTINE trd_icp_init  ***
191      !!
192      !! ** Purpose :   Read the namtrd namelist
193      !!----------------------------------------------------------------------
194      INTEGER  ::   ji, jj, jk   ! dummy loop indices
195      !!----------------------------------------------------------------------
196
197      IF(lwp) THEN
198         WRITE(numout,*)
199         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
200         WRITE(numout,*) '~~~~~~~~~~~~~'
201      ENDIF
202
203      ! Total volume at t-points:
204      tvolt = 0._wp
205      DO jk = 1, jpkm1
206         tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
207      END DO
208      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
209
210      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
211
212#if  defined key_trddyn
213      ! Initialization of potential to kinetic energy conversion
214      rpktrd = 0._wp
215
216      ! Total volume at u-, v- points:
217      tvolu = 0._wp
218      tvolv = 0._wp
219
220      DO jk = 1, jpk
221         DO jj = 2, jpjm1
222            DO ji = fs_2, fs_jpim1   ! vector opt.
223               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)
224               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)
225            END DO
226         END DO
227      END DO
228      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
229      IF( lk_mpp )   CALL mpp_sum( tvolv )
230
231      IF(lwp) THEN
232         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
233         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
234      ENDIF
235#endif
236      !
237   END SUBROUTINE trd_icp_init
238
239
240   SUBROUTINE trd_dwr( kt )
241      !!---------------------------------------------------------------------
242      !!                  ***  ROUTINE trd_dwr  ***
243      !!
244      !! ** Purpose :  write dynamic trends in ocean.output
245      !!----------------------------------------------------------------------
246      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
247      !
248      INTEGER  ::   ji, jj, jk   ! dummy loop indices
249      REAL(wp) ::   zcof         ! local scalar
250      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
251      !!----------------------------------------------------------------------
252
253      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
254
255      ! I. Momentum trends
256      ! -------------------
257
258      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
259
260         ! I.1 Conversion potential energy - kinetic energy
261         ! --------------------------------------------------
262         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
263         zkx  (:,:,:) = 0._wp
264         zky  (:,:,:) = 0._wp
265         zkz  (:,:,:) = 0._wp
266         zkepe(:,:,:) = 0._wp
267   
268         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
269
270         zcof = 0.5_wp / rau0             ! Density flux at w-point
271         zkz(:,:,1) = 0._wp
272         DO jk = 2, jpk
273            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
274         END DO
275         
276         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
277         DO jk = 1, jpkm1
278            DO jj = 1, jpjm1
279               DO ji = 1, jpim1
280                  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) )
281                  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) )
282               END DO
283            END DO
284         END DO
285         
286         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
287            DO jj = 2, jpjm1
288               DO ji = 2, jpim1
289                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
290                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
291                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
292                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
293               END DO
294            END DO
295         END DO
296
297         ! I.2 Basin averaged kinetic energy trend
298         ! ----------------------------------------
299         peke = 0._wp
300         DO jk = 1, jpkm1
301            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
302         END DO
303         peke = grav * peke
304
305         ! I.3 Sums over the global domain
306         ! ---------------------------------
307         IF( lk_mpp ) THEN
308            CALL mpp_sum( peke )
309            CALL mpp_sum( umo , jptot_dyn )
310            CALL mpp_sum( vmo , jptot_dyn )
311            CALL mpp_sum( hke , jptot_dyn )
312         ENDIF
313
314         ! I.2 Print dynamic trends in the ocean.output file
315         ! --------------------------------------------------
316
317         IF(lwp) THEN
318            WRITE (numout,*)
319            WRITE (numout,*)
320            WRITE (numout,9500) kt
321            WRITE (numout,9501) umo(jpdyn_trd_hpg) / tvolu, vmo(jpdyn_trd_hpg) / tvolv
322            WRITE (numout,9509) umo(jpdyn_trd_spg) / tvolu, vmo(jpdyn_trd_spg) / tvolv
323            WRITE (numout,9502) umo(jpdyn_trd_keg) / tvolu, vmo(jpdyn_trd_keg) / tvolv
324            WRITE (numout,9503) umo(jpdyn_trd_rvo) / tvolu, vmo(jpdyn_trd_rvo) / tvolv
325            WRITE (numout,9504) umo(jpdyn_trd_pvo) / tvolu, vmo(jpdyn_trd_pvo) / tvolv
326            WRITE (numout,9507) umo(jpdyn_trd_zad) / tvolu, vmo(jpdyn_trd_zad) / tvolv
327            WRITE (numout,9505) umo(jpdyn_trd_ldf) / tvolu, vmo(jpdyn_trd_ldf) / tvolv
328            WRITE (numout,9508) umo(jpdyn_trd_zdf) / tvolu, vmo(jpdyn_trd_zdf) / tvolv
329            WRITE (numout,9510) umo(jpdyn_trd_tau) / tvolu, vmo(jpdyn_trd_tau) / tvolv
330            WRITE (numout,9511) umo(jpdyn_trd_bfr) / tvolu, vmo(jpdyn_trd_bfr) / tvolv
331            WRITE (numout,9512) umo(jpdyn_trd_atf) / tvolu, vmo(jpdyn_trd_atf) / tvolv
332            WRITE (numout,9513)
333            WRITE (numout,9514)                                                 &
334            &     (  umo(jpdyn_trd_hpg) + umo(jpdyn_trd_spg) + umo(jpdyn_trd_keg) + umo(jpdyn_trd_rvo)   &
335            &      + umo(jpdyn_trd_pvo) + umo(jpdyn_trd_zad) + umo(jpdyn_trd_ldf) + umo(jpdyn_trd_zdf)   &
336            &      + umo(jpdyn_trd_tau) + umo(jpdyn_trd_bfr) + umo(jpdyn_trd_atf) ) / tvolu,   &
337            &     (  vmo(jpdyn_trd_hpg) + vmo(jpdyn_trd_spg) + vmo(jpdyn_trd_keg) + vmo(jpdyn_trd_rvo)   &
338            &      + vmo(jpdyn_trd_pvo) + vmo(jpdyn_trd_zad) + vmo(jpdyn_trd_ldf) + vmo(jpdyn_trd_zdf)   &
339            &      + vmo(jpdyn_trd_tau) + vmo(jpdyn_trd_bfr) + vmo(jpdyn_trd_atf) ) / tvolv
340         ENDIF
341
342 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
343 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
344 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
345 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
346 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
347 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
348 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
349 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
350 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
351 9510    FORMAT(' surface wind stress        u= ', e20.13, '    v= ', e20.13)
352 9511    FORMAT(' bottom friction            u= ', e20.13, '    v= ', e20.13)
353 9512    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
354 9513    FORMAT(' -----------------------------------------------------------------------------')
355 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
356
357         IF(lwp) THEN
358            WRITE (numout,*)
359            WRITE (numout,*)
360            WRITE (numout,9520) kt
361            WRITE (numout,9521) hke(jpdyn_trd_hpg) / tvolt
362            WRITE (numout,9529) hke(jpdyn_trd_spg) / tvolt
363            WRITE (numout,9522) hke(jpdyn_trd_keg) / tvolt
364            WRITE (numout,9523) hke(jpdyn_trd_rvo) / tvolt
365            WRITE (numout,9524) hke(jpdyn_trd_pvo) / tvolt
366            WRITE (numout,9527) hke(jpdyn_trd_zad) / tvolt
367            WRITE (numout,9525) hke(jpdyn_trd_ldf) / tvolt
368            WRITE (numout,9528) hke(jpdyn_trd_zdf) / tvolt
369            WRITE (numout,9530) hke(jpdyn_trd_tau) / tvolt
370            WRITE (numout,9531) hke(jpdyn_trd_bfr) / tvolt
371            WRITE (numout,9532) hke(jpdyn_trd_atf) / tvolt
372            WRITE (numout,9533)
373            WRITE (numout,9534)   &
374            &     (  hke(jpdyn_trd_hpg) + hke(jpdyn_trd_spg) + hke(jpdyn_trd_keg) + hke(jpdyn_trd_rvo)   &
375            &      + hke(jpdyn_trd_pvo) + hke(jpdyn_trd_zad) + hke(jpdyn_trd_ldf) + hke(jpdyn_trd_zdf)   &
376            &      + hke(jpdyn_trd_tau) + hke(jpdyn_trd_bfr) + hke(jpdyn_trd_atf)  ) / tvolt
377         ENDIF
378
379 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
380 9521    FORMAT(' pressure gradient         u2= ', e20.13)
381 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
382 9522    FORMAT(' ke gradient               u2= ', e20.13)
383 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
384 9524    FORMAT(' coriolis term             u2= ', e20.13)
385 9527    FORMAT(' vertical advection        u2= ', e20.13)
386 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
387 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
388 9530    FORMAT(' surface wind stress       u2= ', e20.13)
389 9531    FORMAT(' bottom friction           u2= ', e20.13)
390 9532    FORMAT(' Asselin time filter       u2= ', e20.13)
391 9533    FORMAT(' --------------------------------------------------')
392 9534    FORMAT(' total trend               u2= ', e20.13)
393
394         IF(lwp) THEN
395            WRITE (numout,*)
396            WRITE (numout,*)
397            WRITE (numout,9540) kt
398            WRITE (numout,9541) ( hke(jpdyn_trd_keg) + hke(jpdyn_trd_rvo) + hke(jpdyn_trd_zad) ) / tvolt
399            WRITE (numout,9542) ( hke(jpdyn_trd_keg) + hke(jpdyn_trd_zad) ) / tvolt
400            WRITE (numout,9543) ( hke(jpdyn_trd_pvo) ) / tvolt
401            WRITE (numout,9544) ( hke(jpdyn_trd_rvo) ) / tvolt
402            WRITE (numout,9545) ( hke(jpdyn_trd_spg) ) / tvolt
403            WRITE (numout,9546) ( hke(jpdyn_trd_ldf) ) / tvolt
404            WRITE (numout,9547) ( hke(jpdyn_trd_zdf) ) / tvolt
405            WRITE (numout,9548) ( hke(jpdyn_trd_hpg) ) / tvolt, rpktrd / tvolt
406            WRITE (numout,*)
407            WRITE (numout,*)
408         ENDIF
409
410 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
411 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
412 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
413 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
414 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
415 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
416 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
417 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
418 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
419         !
420         ! Save potential to kinetic energy conversion for next time step
421         rpktrd = peke
422         !
423      ENDIF
424      !
425      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
426      !
427   END SUBROUTINE trd_dwr
428
429
430   SUBROUTINE trd_twr( kt )
431      !!---------------------------------------------------------------------
432      !!                  ***  ROUTINE trd_twr  ***
433      !!
434      !! ** Purpose :  write active tracers trends in ocean.output
435      !!----------------------------------------------------------------------
436      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
437      !
438      INTEGER  ::   jk   ! loop indices
439      !!----------------------------------------------------------------------
440
441      ! I. Tracers trends
442      ! -----------------
443
444      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
445
446         ! I.1 Sums over the global domain
447         ! -------------------------------
448         IF( lk_mpp ) THEN
449            CALL mpp_sum( tmo, jptot_tra )   
450            CALL mpp_sum( smo, jptot_tra )
451            CALL mpp_sum( t2 , jptot_tra )
452            CALL mpp_sum( s2 , jptot_tra )
453         ENDIF
454
455         ! I.2 Print tracers trends in the ocean.output file
456         ! --------------------------------------------------
457         
458         IF(lwp) THEN
459            WRITE (numout,*)
460            WRITE (numout,*)
461            WRITE (numout,9400) kt
462            WRITE (numout,9401)  tmo(jptra_trd_xad) / tvolt, smo(jptra_trd_xad) / tvolt
463            WRITE (numout,9411)  tmo(jptra_trd_yad) / tvolt, smo(jptra_trd_yad) / tvolt
464            WRITE (numout,9402)  tmo(jptra_trd_zad) / tvolt, smo(jptra_trd_zad) / tvolt
465            WRITE (numout,9403)  tmo(jptra_trd_ldf) / tvolt, smo(jptra_trd_ldf) / tvolt
466            WRITE (numout,9404)  tmo(jptra_trd_zdf) / tvolt, smo(jptra_trd_zdf) / tvolt
467            WRITE (numout,9405)  tmo(jptra_trd_npc) / tvolt, smo(jptra_trd_npc) / tvolt
468            WRITE (numout,9406)  tmo(jptra_trd_dmp) / tvolt, smo(jptra_trd_dmp) / tvolt
469            WRITE (numout,9407)  tmo(jptra_trd_qsr) / tvolt
470            WRITE (numout,9408)  tmo(jptra_trd_nsr) / tvolt, smo(jptra_trd_nsr) / tvolt
471            WRITE (numout,9409) 
472            WRITE (numout,9410) (  tmo(jptra_trd_xad) + tmo(jptra_trd_yad) + tmo(jptra_trd_zad) + tmo(jptra_trd_ldf) + tmo(jptra_trd_zdf)   &
473            &                    + tmo(jptra_trd_npc) + tmo(jptra_trd_dmp) + tmo(jptra_trd_qsr) + tmo(jptra_trd_nsr) ) / tvolt,   &
474            &                   (  smo(jptra_trd_xad) + smo(jptra_trd_yad) + smo(jptra_trd_zad) + smo(jptra_trd_ldf) + smo(jptra_trd_zdf)   &
475            &                    + smo(jptra_trd_npc) + smo(jptra_trd_dmp)                   + smo(jptra_trd_nsr) ) / tvolt
476         ENDIF
477
4789400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
479              '              salinity',/' ============================')
4809401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4819411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4829402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4839403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4849404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4859405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4869406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4879407     FORMAT(' penetrative qsr             ',e20.13)
4889408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4899409     FORMAT(' -------------------------------------------------------------------------')
4909410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
491
492
493         IF(lwp) THEN
494            WRITE (numout,*)
495            WRITE (numout,*)
496            WRITE (numout,9420) kt
497            WRITE (numout,9421)   t2(jptra_trd_xad) / tvolt, s2(jptra_trd_xad) / tvolt
498            WRITE (numout,9431)   t2(jptra_trd_yad) / tvolt, s2(jptra_trd_yad) / tvolt
499            WRITE (numout,9422)   t2(jptra_trd_zad) / tvolt, s2(jptra_trd_zad) / tvolt
500            WRITE (numout,9423)   t2(jptra_trd_ldf) / tvolt, s2(jptra_trd_ldf) / tvolt
501            WRITE (numout,9424)   t2(jptra_trd_zdf) / tvolt, s2(jptra_trd_zdf) / tvolt
502            WRITE (numout,9425)   t2(jptra_trd_npc) / tvolt, s2(jptra_trd_npc) / tvolt
503            WRITE (numout,9426)   t2(jptra_trd_dmp) / tvolt, s2(jptra_trd_dmp) / tvolt
504            WRITE (numout,9427)   t2(jptra_trd_qsr) / tvolt
505            WRITE (numout,9428)   t2(jptra_trd_nsr) / tvolt, s2(jptra_trd_nsr) / tvolt
506            WRITE (numout,9429)
507            WRITE (numout,9430) (  t2(jptra_trd_xad) + t2(jptra_trd_yad) + t2(jptra_trd_zad) + t2(jptra_trd_ldf) + t2(jptra_trd_zdf)   &
508            &                    + t2(jptra_trd_npc) + t2(jptra_trd_dmp) + t2(jptra_trd_qsr) + t2(jptra_trd_nsr) ) / tvolt,   &
509            &                   (  s2(jptra_trd_xad) + s2(jptra_trd_yad) + s2(jptra_trd_zad) + s2(jptra_trd_ldf) + s2(jptra_trd_zdf)   &
510            &                    + s2(jptra_trd_npc) + s2(jptra_trd_dmp)                  + s2(jptra_trd_nsr) ) / tvolt
511         ENDIF
512
5139420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
514            '               salinity', /, ' ===============================')
5159421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
5169431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
5179422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
5189423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
5199424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
5209425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
5219426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
5229427     FORMAT(' penetrative qsr           * t   ', e20.13)
5239428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
5249429     FORMAT(' -----------------------------------------------------------------------------')
5259430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
526
527
528         IF(lwp) THEN
529            WRITE (numout,*)
530            WRITE (numout,*)
531            WRITE (numout,9440) kt
532            WRITE (numout,9441) ( tmo(jptra_trd_xad)+tmo(jptra_trd_yad)+tmo(jptra_trd_zad) )/tvolt,   &
533            &                   ( smo(jptra_trd_xad)+smo(jptra_trd_yad)+smo(jptra_trd_zad) )/tvolt
534            WRITE (numout,9442)   tmo(jptra_trd_sad)/tvolt,  smo(jptra_trd_sad)/tvolt
535            WRITE (numout,9443)   tmo(jptra_trd_ldf)/tvolt,  smo(jptra_trd_ldf)/tvolt
536            WRITE (numout,9444)   tmo(jptra_trd_zdf)/tvolt,  smo(jptra_trd_zdf)/tvolt
537            WRITE (numout,9445)   tmo(jptra_trd_npc)/tvolt,  smo(jptra_trd_npc)/tvolt
538            WRITE (numout,9446) ( t2(jptra_trd_xad)+t2(jptra_trd_yad)+t2(jptra_trd_zad) )/tvolt,    &
539            &                   ( s2(jptra_trd_xad)+s2(jptra_trd_yad)+s2(jptra_trd_zad) )/tvolt
540            WRITE (numout,9447)   t2(jptra_trd_ldf)/tvolt,   s2(jptra_trd_ldf)/tvolt
541            WRITE (numout,9448)   t2(jptra_trd_zdf)/tvolt,   s2(jptra_trd_zdf)/tvolt
542            WRITE (numout,9449)   t2(jptra_trd_npc)/tvolt,   s2(jptra_trd_npc)/tvolt
543         ENDIF
544
5459440     FORMAT(' tracer consistency at it= ',i6,   &
546            ' :         temperature','                salinity',/,   &
547            ' ==================================')
5489441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5499442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5509443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5519444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5529445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5539446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5549447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5559448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5569449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
557         !
558      ENDIF
559      !
560   END SUBROUTINE trd_twr
561
562   !!======================================================================
563END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.