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.
trdicp.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/trdicp.F90 @ 3316

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

Ediag branche: #927 add 3D output for dyn & tracer trends
  • Property svn:keywords set to Id
File size: 31.8 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#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   'key_trdtra'  or                  active tracers trends diagnostics
12   !!   'key_trddyn'                            momentum trends diagnostics
13   !!----------------------------------------------------------------------
14   !!   trd_icp          : compute the basin averaged properties for tra/dyn
15   !!   trd_dwr          : print dynmaic trends in ocean.output file
16   !!   trd_twr          : print tracers trends in ocean.output file
17   !!   trd_icp_init     : initialization step
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
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 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   INTERFACE trd_icp
37      MODULE PROCEDURE trd_2d, trd_3d
38   END INTERFACE
39
40   PUBLIC   trd_icp       ! called by trdmod.F90
41   PUBLIC   trd_dwr       ! called by step.F90
42   PUBLIC   trd_twr       ! called by step.F90
43   PUBLIC   trd_icp_init  ! called by opa.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48#  include "zdfddm_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )
57      !!---------------------------------------------------------------------
58      !!                  ***  ROUTINE trd_2d  ***
59      !!
60      !! ** Purpose :   compute and print the domain averaged properties of tracers
61      !!              and/or momentum equations at each nn_trd time step.
62      !!----------------------------------------------------------------------
63      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ptrd2dx   ! Temperature or U trend
64      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ptrd2dy   ! Salinity    or V trend
65      INTEGER                 , INTENT(in   ) ::   ktrd      ! tracer trend index
66      CHARACTER(len=3)        , INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
67      !!
68      INTEGER ::   ji, jj   ! loop indices
69      !!----------------------------------------------------------------------
70
71      SELECT CASE( ctype )    !==  Mask trends  ==!
72      !
73      CASE( 'DYN' )                    ! Momentum
74         DO jj = 1, jpjm1
75            DO ji = 1, jpim1
76               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
77               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
78            END DO
79         END DO
80         ptrd2dx(jpi, : ) = 0._wp      ;      ptrd2dy(jpi, : ) = 0._wp
81         ptrd2dx( : ,jpj) = 0._wp      ;      ptrd2dy( : ,jpj) = 0._wp
82         !
83      CASE( 'TRA' )                    ! Tracers
84         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
85         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
86         !
87      END SELECT
88     
89      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
90      !
91      CASE( 'DYN' )                    ! Momentum
92         umo(ktrd) = 0._wp
93         vmo(ktrd) = 0._wp
94         !
95         SELECT CASE( ktrd )
96         CASE( jpdyn_trd_swf )         ! surface forcing
97            umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
98            vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
99         END SELECT
100         !
101      CASE( 'TRA' )              ! Tracers
102         tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
103         smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
104      END SELECT
105     
106      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
107      !
108      CASE( 'DYN' )              ! Momentum
109         hke(ktrd) = SUM(   un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1)   &
110            &             + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1)   )
111         !
112      CASE( 'TRA' )              ! Tracers
113         t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_tem) )
114         s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_sal) )
115         !     
116      END SELECT
117      !
118   END SUBROUTINE trd_2d
119
120
121   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype )
122      !!---------------------------------------------------------------------
123      !!                  ***  ROUTINE trd_3d  ***
124      !!
125      !! ** Purpose : verify the basin averaged properties of tracers and/or
126      !!              momentum equations at every time step frequency nn_trd.
127      !!----------------------------------------------------------------------
128      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrd3dx   ! Temperature or U trend
129      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrd3dy   ! Salinity    or V trend
130      INTEGER,                    INTENT(in   ) ::   ktrd      ! momentum or tracer trend index
131      CHARACTER(len=3),           INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
132      !!
133      INTEGER  ::   ji, jj, jk   ! dummy loop indices
134      !!----------------------------------------------------------------------
135
136      SELECT CASE( ctype )    !==  Mask the trends  ==!
137      !
138      CASE( 'DYN' )              ! Momentum       
139         DO jk = 1, jpkm1
140            DO jj = 1, jpjm1
141               DO ji = 1, jpim1
142                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
143                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
144               END DO
145            END DO
146         END DO
147         ptrd3dx(jpi, : ,:) = 0._wp      ;      ptrd3dy(jpi, : ,:) = 0._wp
148         ptrd3dx( : ,jpj,:) = 0._wp      ;      ptrd3dy( : ,jpj,:) = 0._wp
149         !
150      CASE( 'TRA' )              ! Tracers
151         DO jk = 1, jpkm1
152            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
153            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
154         END DO
155         !
156      END SELECT   
157
158      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
159      !
160      CASE( 'DYN' )              ! Momentum
161         umo(ktrd) = 0._wp
162         vmo(ktrd) = 0._wp
163         DO jk = 1, jpkm1
164            umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
165            vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
166         END DO
167         !
168      CASE( 'TRA' )              ! Tracers
169         tmo(ktrd) = 0._wp
170         smo(ktrd) = 0._wp
171         DO jk = 1, jpkm1
172            tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
173            smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
174         END DO
175         !
176      END SELECT
177
178      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
179      !
180      CASE( 'DYN' )              ! Momentum
181         hke(ktrd) = 0._wp
182         DO jk = 1, jpkm1
183            hke(ktrd) = hke(ktrd) + SUM(   un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk)   &
184               &                         + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk)   )
185         END DO
186         !
187      CASE( 'TRA' )              ! Tracers
188         t2(ktrd) = 0._wp
189         s2(ktrd) = 0._wp
190         DO jk = 1, jpkm1
191            t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tsn(:,:,jk,jp_tem) * e1e2t(:,:) * fse3t(:,:,jk) )
192            s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * tsn(:,:,jk,jp_sal) * e1e2t(:,:) * fse3t(:,:,jk) )
193         END DO
194         !
195      END SELECT
196      !
197   END SUBROUTINE trd_3d
198
199
200   SUBROUTINE trd_icp_init
201      !!---------------------------------------------------------------------
202      !!                  ***  ROUTINE trd_icp_init  ***
203      !!
204      !! ** Purpose :   Read the namtrd namelist
205      !!----------------------------------------------------------------------
206      INTEGER  ::   ji, jj, jk   ! dummy loop indices
207      !!----------------------------------------------------------------------
208
209      IF(lwp) THEN
210         WRITE(numout,*)
211         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
212         WRITE(numout,*) '~~~~~~~~~~~~~'
213      ENDIF
214
215      ! Total volume at t-points:
216      tvolt = 0._wp
217      DO jk = 1, jpkm1
218         tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
219      END DO
220      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
221
222      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
223
224#if  defined key_trddyn
225      ! Initialization of potential to kinetic energy conversion
226      rpktrd = 0._wp
227
228      ! Total volume at u-, v- points:
229      tvolu = 0._wp
230      tvolv = 0._wp
231
232      DO jk = 1, jpk
233         DO jj = 2, jpjm1
234            DO ji = fs_2, fs_jpim1   ! vector opt.
235               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)
236               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)
237            END DO
238         END DO
239      END DO
240      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
241      IF( lk_mpp )   CALL mpp_sum( tvolv )
242
243      IF(lwp) THEN
244         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
245         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
246      ENDIF
247#endif
248      !
249   END SUBROUTINE trd_icp_init
250
251
252   SUBROUTINE trd_dwr( kt )
253      !!---------------------------------------------------------------------
254      !!                  ***  ROUTINE trd_dwr  ***
255      !!
256      !! ** Purpose :  write dynamic trends in ocean.output
257      !!----------------------------------------------------------------------
258      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
259      !
260      INTEGER  ::   ji, jj, jk   ! dummy loop indices
261      REAL(wp) ::   zcof         ! local scalar
262      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
263      !!----------------------------------------------------------------------
264
265      IF( ln_3D_trd_t .AND. ln_3D_trd_d )   RETURN            ! do nothing if 3D output with IOM
266
267      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
268
269      ! I. Momentum trends
270      ! -------------------
271
272      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
273
274         ! I.1 Conversion potential energy - kinetic energy
275         ! --------------------------------------------------
276         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
277         zkx  (:,:,:) = 0._wp
278         zky  (:,:,:) = 0._wp
279         zkz  (:,:,:) = 0._wp
280         zkepe(:,:,:) = 0._wp
281   
282         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
283
284         zcof = 0.5_wp / rau0             ! Density flux at w-point
285         zkz(:,:,1) = 0._wp
286         DO jk = 2, jpk
287            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
288         END DO
289         
290         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
291         DO jk = 1, jpkm1
292            DO jj = 1, jpjm1
293               DO ji = 1, jpim1
294                  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) )
295                  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) )
296               END DO
297            END DO
298         END DO
299         
300         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
301            DO jj = 2, jpjm1
302               DO ji = 2, jpim1
303                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
304                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
305                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
306                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
307               END DO
308            END DO
309         END DO
310
311         ! I.2 Basin averaged kinetic energy trend
312         ! ----------------------------------------
313         peke = 0._wp
314         DO jk = 1, jpkm1
315            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
316         END DO
317         peke = grav * peke
318
319         ! I.3 Sums over the global domain
320         ! ---------------------------------
321         IF( lk_mpp ) THEN
322            CALL mpp_sum( peke )
323            CALL mpp_sum( umo , jptot_dyn )
324            CALL mpp_sum( vmo , jptot_dyn )
325            CALL mpp_sum( hke , jptot_dyn )
326         ENDIF
327
328         ! I.2 Print dynamic trends in the ocean.output file
329         ! --------------------------------------------------
330
331         IF(lwp) THEN
332            WRITE (numout,*)
333            WRITE (numout,*)
334            WRITE (numout,9500) kt
335            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
336            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
337            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
338            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
339            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
340            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
341            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
342            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
343            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
344            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
345            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
346            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
347            WRITE (numout,9513)
348            WRITE (numout,9514)                                                 &
349            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
350            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
351            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
352            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
353            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
354            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
355         ENDIF
356
357 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
358 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
359 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
360 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
361 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
362 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
363 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
364 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
365 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
366 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
367 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
368 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
369 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
370 9513    FORMAT(' -----------------------------------------------------------------------------')
371 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
372
373         IF(lwp) THEN
374            WRITE (numout,*)
375            WRITE (numout,*)
376            WRITE (numout,9520) kt
377            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
378            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
379            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
380            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
381            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
382            WRITE (numout,9526) hke(jpicpd_had) / tvolt
383            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
384            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
385            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
386            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
387            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
388            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
389            WRITE (numout,9533)
390            WRITE (numout,9534)   &
391            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
392            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
393            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
394         ENDIF
395
396 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
397 9521    FORMAT(' pressure gradient         u2= ', e20.13)
398 9522    FORMAT(' ke gradient               u2= ', e20.13)
399 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
400 9524    FORMAT(' coriolis term             u2= ', e20.13)
401 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
402 9526    FORMAT(' horizontal advection      u2= ', e20.13)
403 9527    FORMAT(' vertical advection        u2= ', e20.13)
404 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
405 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
406 9530    FORMAT(' surface wind forcing      u2= ', e20.13)
407 9531    FORMAT(' dampimg term              u2= ', e20.13)
408 9532    FORMAT(' bottom flux               u2= ', e20.13)
409 9533    FORMAT(' --------------------------------------------------')
410 9534    FORMAT(' total trend               u2= ', e20.13)
411
412         IF(lwp) THEN
413            WRITE (numout,*)
414            WRITE (numout,*)
415            WRITE (numout,9540) kt
416            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
417            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
418            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
419            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
420            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
421            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
422            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
423            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
424            WRITE (numout,*)
425            WRITE (numout,*)
426         ENDIF
427
428 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
429 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
430 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
431 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
432 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
433 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
434 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
435 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
436 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
437         !
438         ! Save potential to kinetic energy conversion for next time step
439         rpktrd = peke
440         !
441      ENDIF
442      !
443      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
444      !
445   END SUBROUTINE trd_dwr
446
447
448   SUBROUTINE trd_twr( kt )
449      !!---------------------------------------------------------------------
450      !!                  ***  ROUTINE trd_twr  ***
451      !!
452      !! ** Purpose :  write active tracers trends in ocean.output
453      !!----------------------------------------------------------------------
454      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
455      !
456      INTEGER  ::   jk   ! loop indices
457      REAL(wp), POINTER, DIMENSION(:,:,:)  ::   zwt, zws, ztrdt, ztrds   ! 3D workspace
458      !!----------------------------------------------------------------------
459
460
461      IF( ln_3D_trd_t ) THEN      ! 3D output: treat the vertical diffusion trends (if iso)
462         !
463         CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt, ztrds )
464         !
465         IF( ln_traldf_iso ) THEN      ! iso-neutral diffusion : re-compute the PURE vertical diffusive trend
466            !                                 !  zdf trends using now field (called after the swap)
467            zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
468            zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
469            DO jk = 2, jpk
470               zwt(:,:,jk) =   avt(:,:,jk) * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
471               zws(:,:,jk) = fsavs(:,:,jk) * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk)
472            END DO
473            !
474            ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
475            DO jk = 1, jpkm1
476               ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
477               ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk)
478            END DO
479            CALL iom_put( "ttrd_zdfp", ztrdt )        ! PURE vertical diffusion (no isoneutral contribution)
480            CALL iom_put( "strd_zdfp", ztrds )
481         ENDIF
482         !
483         CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt, ztrds )
484         !
485         RETURN                     ! do nothing else if 3D output with IOM
486         !
487      ENDIF
488
489
490      ! I. Tracers trends
491      ! -----------------
492
493      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
494
495         ! I.1 Sums over the global domain
496         ! -------------------------------
497         IF( lk_mpp ) THEN
498            CALL mpp_sum( tmo, jptot_tra )   
499            CALL mpp_sum( smo, jptot_tra )
500            CALL mpp_sum( t2 , jptot_tra )
501            CALL mpp_sum( s2 , jptot_tra )
502         ENDIF
503
504         ! I.2 Print tracers trends in the ocean.output file
505         ! --------------------------------------------------
506         
507         IF(lwp) THEN
508            WRITE (numout,*)
509            WRITE (numout,*)
510            WRITE (numout,9400) kt
511            WRITE (numout,9401)  tmo(jpicpt_xad) / tvolt, smo(jpicpt_xad) / tvolt
512            WRITE (numout,9411)  tmo(jpicpt_yad) / tvolt, smo(jpicpt_yad) / tvolt
513            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
514            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
515            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
516            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
517            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
518            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
519            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
520            WRITE (numout,9409) 
521            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
522            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
523            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
524            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
525         ENDIF
526
5279400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
528              '              salinity',/' ============================')
5299401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
5309411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
5319402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
5329403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
5339404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
5349405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
5359406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
5369407     FORMAT(' penetrative qsr             ',e20.13)
5379408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
5389409     FORMAT(' -------------------------------------------------------------------------')
5399410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
540
541
542         IF(lwp) THEN
543            WRITE (numout,*)
544            WRITE (numout,*)
545            WRITE (numout,9420) kt
546            WRITE (numout,9421)   t2(jpicpt_xad) / tvolt, s2(jpicpt_xad) / tvolt
547            WRITE (numout,9431)   t2(jpicpt_yad) / tvolt, s2(jpicpt_yad) / tvolt
548            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
549            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
550            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
551            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
552            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
553            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
554            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
555            WRITE (numout,9429)
556            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
557            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
558            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
559            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
560         ENDIF
561
5629420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
563            '               salinity', /, ' ===============================')
5649421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
5659431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
5669422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
5679423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
5689424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
5699425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
5709426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
5719427     FORMAT(' penetrative qsr           * t   ', e20.13)
5729428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
5739429     FORMAT(' -----------------------------------------------------------------------------')
5749430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
575
576
577         IF(lwp) THEN
578            WRITE (numout,*)
579            WRITE (numout,*)
580            WRITE (numout,9440) kt
581            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
582            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
583            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
584            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
585            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
586            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
587            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
588            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
589            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
590            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
591            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
592         ENDIF
593
5949440     FORMAT(' tracer consistency at it= ',i6,   &
595            ' :         temperature','                salinity',/,   &
596            ' ==================================')
5979441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5989442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5999443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
6009444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
6019445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
6029446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
6039447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
6049448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
6059449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
606         !
607      ENDIF
608      !
609   END SUBROUTINE trd_twr
610
611#   else
612   !!----------------------------------------------------------------------
613   !!   Default case :                                         Empty module
614   !!----------------------------------------------------------------------
615   INTERFACE trd_icp
616      MODULE PROCEDURE trd_2d, trd_3d
617   END INTERFACE
618
619CONTAINS
620   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )       ! Empty routine
621      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
622      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
623      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
624      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
625          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
626   END SUBROUTINE trd_2d
627   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype )       ! Empty routine
628      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
629      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
630      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
631      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
632          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
633   END SUBROUTINE trd_3d
634   SUBROUTINE trd_icp_init               ! Empty routine
635   END SUBROUTINE trd_icp_init
636   SUBROUTINE trd_dwr( kt )          ! Empty routine
637      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
638   END SUBROUTINE trd_dwr
639   SUBROUTINE trd_twr( kt )          ! Empty routine
640      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
641   END SUBROUTINE trd_twr
642#endif
643
644   !!======================================================================
645END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.