source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

  • Property svn:keywords set to Id
File size: 34.3 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   !!----------------------------------------------------------------------
8#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
9   !!----------------------------------------------------------------------
10   !!   'key_trdtra'  or                  active tracers trends diagnostics
11   !!   'key_trddyn'                            momentum trends diagnostics
12   !!----------------------------------------------------------------------
13   !!   trd_icp          : compute the basin averaged properties for tra/dyn
14   !!   trd_dwr          : print dynmaic trends in ocean.output file
15   !!   trd_twr          : print tracers trends in ocean.output file
16   !!   trd_icp_init     : initialization step
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE trdmod_oce      ! ocean variables trends
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 in_out_manager  ! I/O manager
25   USE lib_mpp         ! distibuted memory computing library
26   USE eosbn2          ! equation of state
27   USE phycst          ! physical constants
28
29   IMPLICIT NONE
30   PRIVATE
31
32   INTERFACE trd_icp
33      MODULE PROCEDURE trd_2d, trd_3d
34   END INTERFACE
35
36   PUBLIC   trd_icp       ! called by trdmod.F90
37   PUBLIC   trd_dwr       ! called by step.F90
38   PUBLIC   trd_twr       ! called by step.F90
39   PUBLIC   trd_icp_init  ! called by opa.F90
40
41   !! * Control permutation of array indices
42#  include "oce_ftrans.h90"
43#  include "dom_oce_ftrans.h90"
44#  include "trdmld_oce_ftrans.h90"
45#  include "ldftra_oce_ftrans.h90"
46#  include "ldfdyn_oce_ftrans.h90"
47#  include "zdf_oce_ftrans.h90"
48
49   !! * Substitutions
50#  include "domzgr_substitute.h90"
51#  include "vectopt_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )
60      !!---------------------------------------------------------------------
61      !!                  ***  ROUTINE trd_2d  ***
62      !!
63      !! ** Purpose : verify the basin averaged properties of tracers and/or
64      !!              momentum equations at every time step frequency nn_trd.
65      !!----------------------------------------------------------------------
66      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx   ! Temperature or U trend
67      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy   ! Salinity    or V trend
68      INTEGER                     , INTENT(in   ) ::   ktrd      ! tracer trend index
69      CHARACTER(len=3)            , INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
70      !!
71      INTEGER ::   ji, jj   ! loop indices
72      !!----------------------------------------------------------------------
73
74      SELECT CASE( ctype )    !==  Mask trends  ==!
75      !
76      CASE( 'DYN' )                    ! Momentum
77         DO jj = 1, jpjm1
78            DO ji = 1, jpim1
79               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
80               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
81            END DO
82         END DO
83         ptrd2dx(jpi, : ) = 0._wp      ;      ptrd2dy(jpi, : ) = 0._wp
84         ptrd2dx( : ,jpj) = 0._wp      ;      ptrd2dy( : ,jpj) = 0._wp
85         !
86      CASE( 'TRA' )                    ! Tracers
87         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
88         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
89         !
90      END SELECT
91     
92      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
93      !
94      CASE( 'DYN' )                    ! Momentum
95         umo(ktrd) = 0._wp
96         vmo(ktrd) = 0._wp
97         !
98         SELECT CASE( ktrd )
99         CASE( jpdyn_trd_swf )         ! surface forcing
100            umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
101            vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
102         END SELECT
103         !
104      CASE( 'TRA' )              ! Tracers
105         tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
106         smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
107      END SELECT
108     
109      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
110      !
111      CASE( 'DYN' )              ! Momentum
112         hke(ktrd) = SUM(   un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1)   &
113            &             + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1)   )
114         !
115      CASE( 'TRA' )              ! Tracers
116         t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tn(:,:,1) )
117         s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * sn(:,:,1) )
118         !     
119      END SELECT
120      !
121   END SUBROUTINE trd_2d
122
123
124   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype )
125      !!---------------------------------------------------------------------
126      !!                  ***  ROUTINE trd_3d  ***
127      !!
128      !! ** Purpose : verify the basin averaged properties of tracers and/or
129      !!              momentum equations at every time step frequency nn_trd.
130      !!----------------------------------------------------------------------
131
132      !! DCSE_NEMO: This style defeats ftrans
133!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx   ! Temperature or U trend
134!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy   ! Salinity    or V trend
135
136!FTRANS ptrd3dx ptrd3dy :I :I :z
137      REAL(wp), INTENT(inout) ::   ptrd3dx(jpi,jpj,jpkorig)   ! Temperature or U trend
138      REAL(wp), INTENT(inout) ::   ptrd3dy(jpi,jpj,jpkorig)   ! Salinity    or V trend
139
140      INTEGER,                          INTENT(in   ) ::   ktrd      ! momentum or tracer trend index
141      CHARACTER(len=3),                 INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
142      !!
143      INTEGER  ::   ji, jj, jk   ! dummy loop indices
144      !!----------------------------------------------------------------------
145
146      SELECT CASE( ctype )    !==  Mask the trends  ==!
147      !
148      CASE( 'DYN' )              ! Momentum       
149#if defined key_z_first
150         DO jj = 1, jpjm1
151            DO ji = 1, jpim1
152               DO jk = 1, jpkm1
153#else
154         DO jk = 1, jpkm1
155            DO jj = 1, jpjm1
156               DO ji = 1, jpim1
157#endif
158                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
159                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
160               END DO
161            END DO
162         END DO
163         ptrd3dx(jpi, : ,:) = 0._wp      ;      ptrd3dy(jpi, : ,:) = 0._wp
164         ptrd3dx( : ,jpj,:) = 0._wp      ;      ptrd3dy( : ,jpj,:) = 0._wp
165         !
166      CASE( 'TRA' )              ! Tracers
167#if defined key_z_first
168         DO jj = 1, jpj
169            DO ji = 1, jpi
170               DO jk = 1, jpkm1
171                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
172                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
173               END DO
174            END DO
175         END DO
176#else
177         DO jk = 1, jpkm1
178            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
179            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
180         END DO
181#endif
182         !
183      END SELECT   
184
185      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
186      !
187      CASE( 'DYN' )              ! Momentum
188         umo(ktrd) = 0._wp
189         vmo(ktrd) = 0._wp
190#if defined key_z_first
191         !! DCSE_NEMO: this changes the order of summation
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194               DO jk = 1, jpkm1
195                  umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
196                  vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
197               END DO
198            END DO
199         END DO
200#else
201         DO jk = 1, jpkm1
202            umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
203            vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
204         END DO
205#endif
206         !
207      CASE( 'TRA' )              ! Tracers
208         tmo(ktrd) = 0._wp
209         smo(ktrd) = 0._wp
210#if defined key_z_first
211         !! DCSE_NEMO: this changes the order of summation
212         DO jj = 1, jpj
213            DO ji = 1, jpi
214               DO jk = 1, jpkm1
215                  tmo(ktrd) = tmo(ktrd) + ptrd3dx(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk) )
216                  smo(ktrd) = smo(ktrd) + ptrd3dy(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk) )
217               END DO
218            END DO
219         END DO
220#else
221         DO jk = 1, jpkm1
222            tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
223            smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
224         END DO
225#endif
226         !
227      END SELECT
228
229      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
230      !
231      CASE( 'DYN' )              ! Momentum
232         hke(ktrd) = 0._wp
233#if defined key_z_first
234         !! DCSE_NEMO: this changes the order of summation
235         DO jj = 1, jpj
236            DO ji = 1, jpi
237               DO jk = 1, jpkm1
238                  hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)   &
239                     &                  + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
240               END DO
241            END DO
242         END DO
243#else
244         DO jk = 1, jpkm1
245            hke(ktrd) = hke(ktrd) + SUM(   un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk)   &
246               &                         + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk)   )
247         END DO
248#endif
249         !
250      CASE( 'TRA' )              ! Tracers
251         t2(ktrd) = 0._wp
252         s2(ktrd) = 0._wp
253#if defined key_z_first
254         !! DCSE_NEMO: this changes the order of summation
255         DO jj = 1, jpj
256            DO ji = 1, jpi
257               DO jk = 1, jpkm1
258                  t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * tn(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
259                  s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * sn(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
260               END DO
261            END DO
262         END DO
263#else
264         DO jk = 1, jpkm1
265            t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
266            s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * sn(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
267         END DO
268#endif
269         !
270      END SELECT
271      !
272   END SUBROUTINE trd_3d
273
274
275   SUBROUTINE trd_icp_init
276      !!---------------------------------------------------------------------
277      !!                  ***  ROUTINE trd_icp_init  ***
278      !!
279      !! ** Purpose :   Read the namtrd namelist
280      !!----------------------------------------------------------------------
281      INTEGER  ::   ji, jj, jk   ! dummy loop indices
282      !!----------------------------------------------------------------------
283
284      IF(lwp) THEN
285         WRITE(numout,*)
286         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
287         WRITE(numout,*) '~~~~~~~~~~~~~'
288      ENDIF
289
290      ! Total volume at t-points:
291      tvolt = 0._wp
292#if defined key_z_first
293      DO jj = 1, jpj
294         DO ji = 1, jpi
295            DO jk = 1, jpkm1
296               tvolt = tvolt + e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
297            END DO
298         END DO
299      END DO
300#else
301      DO jk = 1, jpkm1
302      !! DCSE_NEMO: This looks plain wrong
303!        tvolt = SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
304         tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
305      END DO
306#endif
307      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
308
309      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
310
311#if  defined key_trddyn
312      ! Initialization of potential to kinetic energy conversion
313      rpktrd = 0._wp
314
315      ! Total volume at u-, v- points:
316      tvolu = 0._wp
317      tvolv = 0._wp
318
319#if defined key_z_first
320      DO jj = 2, jpjm1
321         DO ji = 2, jpim1
322            DO jk = 1, jpk
323#else
324      DO jk = 1, jpk
325         DO jj = 2, jpjm1
326            DO ji = fs_2, fs_jpim1   ! vector opt.
327#endif
328               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)
329               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)
330            END DO
331         END DO
332      END DO
333      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
334      IF( lk_mpp )   CALL mpp_sum( tvolv )
335
336      IF(lwp) THEN
337         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
338         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
339      ENDIF
340#endif
341      !
342   END SUBROUTINE trd_icp_init
343
344
345   SUBROUTINE trd_dwr( kt )
346      !!---------------------------------------------------------------------
347      !!                  ***  ROUTINE trd_dwr  ***
348      !!
349      !! ** Purpose :  write dynamic trends in ocean.output
350      !!----------------------------------------------------------------------
351      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
352      USE wrk_nemo, ONLY:   zkepe => wrk_3d_1 , zkx => wrk_3d_2   ! 3D workspace
353      USE wrk_nemo, ONLY:   zky   => wrk_3d_3 , zkz => wrk_3d_4   !  -      -
354
355      !! DCSE_NEMO: need additional directives for renamed module variables
356!FTRANS zkepe zkx zky zkz :I :I :z
357
358      !
359      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
360      !
361      INTEGER  ::   ji, jj, jk   ! dummy loop indices
362      REAL(wp) ::   zcof         ! local scalar
363      !!----------------------------------------------------------------------
364
365      IF( wrk_in_use(3, 1,2,3,4) ) THEN
366         CALL ctl_stop('trd_dwr: requested workspace arrays unavailable')   ;   RETURN
367      ENDIF
368
369      ! I. Momentum trends
370      ! -------------------
371
372      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
373
374         ! I.1 Conversion potential energy - kinetic energy
375         ! --------------------------------------------------
376         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
377         zkx  (:,:,:) = 0._wp
378         zky  (:,:,:) = 0._wp
379         zkz  (:,:,:) = 0._wp
380         zkepe(:,:,:) = 0._wp
381   
382         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
383
384         zcof = 0.5_wp / rau0             ! Density flux at w-point
385#if defined key_z_first
386         DO jj = 1, jpj
387            DO ji = 1, jpi
388               zkz(ji,jj,1) = 0._wp
389               DO jk = 2, jpk
390                  zkz(ji,jj,jk) = e1e2t(ji,jj) * wn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) * tmask_i(ji,jj)
391               END DO
392            END DO
393         END DO
394#else
395         zkz(:,:,1) = 0._wp
396         DO jk = 2, jpk
397            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
398         END DO
399#endif
400         
401         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
402#if defined key_z_first
403         DO jj = 1, jpjm1
404            DO ji = 1, jpim1
405               DO jk = 1, jpkm1
406#else
407         DO jk = 1, jpkm1
408            DO jj = 1, jpjm1
409               DO ji = 1, jpim1
410#endif
411                  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) )
412                  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) )
413               END DO
414            END DO
415         END DO
416         
417#if defined key_z_first
418         DO jj = 2, jpjm1                 ! Density flux divergence at t-point
419            DO ji = 2, jpim1
420               DO jk = 1, jpkm1
421#else
422         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
423            DO jj = 2, jpjm1
424               DO ji = 2, jpim1
425#endif
426                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
427                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
428                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
429                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
430               END DO
431            END DO
432         END DO
433
434         ! I.2 Basin averaged kinetic energy trend
435         ! ----------------------------------------
436         peke = 0._wp
437#if defined key_z_first
438         DO jj = 1, jpj
439            DO ji = 1, jpi
440               DO jk = 1, jpkm1
441                  peke = peke + zkepe(ji,jj,jk) * fsdept(ji,jj,jk) * e1e2t(ji,jj) * fse3t(ji,jj,jk)
442               END DO
443            END DO
444         END DO
445#else
446         DO jk = 1, jpkm1
447            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
448         END DO
449#endif
450         peke = grav * peke
451
452         ! I.3 Sums over the global domain
453         ! ---------------------------------
454         IF( lk_mpp ) THEN
455            CALL mpp_sum( peke )
456            CALL mpp_sum( umo , jptot_dyn )
457            CALL mpp_sum( vmo , jptot_dyn )
458            CALL mpp_sum( hke , jptot_dyn )
459         ENDIF
460
461         ! I.2 Print dynamic trends in the ocean.output file
462         ! --------------------------------------------------
463
464         IF(lwp) THEN
465            WRITE (numout,*)
466            WRITE (numout,*)
467            WRITE (numout,9500) kt
468            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
469            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
470            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
471            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
472            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
473            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
474            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
475            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
476            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
477            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
478            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
479            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
480            WRITE (numout,9513)
481            WRITE (numout,9514)                                                 &
482            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
483            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
484            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
485            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
486            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
487            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
488         ENDIF
489
490 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
491 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
492 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
493 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
494 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
495 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
496 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
497 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
498 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
499 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
500 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
501 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
502 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
503 9513    FORMAT(' -----------------------------------------------------------------------------')
504 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
505
506         IF(lwp) THEN
507            WRITE (numout,*)
508            WRITE (numout,*)
509            WRITE (numout,9520) kt
510            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
511            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
512            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
513            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
514            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
515            WRITE (numout,9526) hke(jpicpd_had) / tvolt
516            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
517            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
518            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
519            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
520            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
521            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
522            WRITE (numout,9533)
523            WRITE (numout,9534)   &
524            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
525            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
526            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
527         ENDIF
528
529 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
530 9521    FORMAT(' pressure gradient         u2= ', e20.13)
531 9522    FORMAT(' ke gradient               u2= ', e20.13)
532 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
533 9524    FORMAT(' coriolis term             u2= ', e20.13)
534 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
535 9526    FORMAT(' horizontal advection      u2= ', e20.13)
536 9527    FORMAT(' vertical advection        u2= ', e20.13)
537 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
538 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
539 9530    FORMAT(' surface wind forcing      u2= ', e20.13)
540 9531    FORMAT(' dampimg term              u2= ', e20.13)
541 9532    FORMAT(' bottom flux               u2= ', e20.13)
542 9533    FORMAT(' --------------------------------------------------')
543 9534    FORMAT(' total trend               u2= ', e20.13)
544
545         IF(lwp) THEN
546            WRITE (numout,*)
547            WRITE (numout,*)
548            WRITE (numout,9540) kt
549            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
550            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
551            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
552            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
553            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
554            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
555            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
556            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
557            WRITE (numout,*)
558            WRITE (numout,*)
559         ENDIF
560
561 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
562 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
563 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
564 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
565 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
566 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
567 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
568 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
569 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
570         !
571         ! Save potential to kinetic energy conversion for next time step
572         rpktrd = peke
573         !
574      ENDIF
575      !
576      IF( wrk_not_released(3, 1,2,3,4) )   CALL ctl_stop('trd_dwr: failed to release workspace arrays')
577      !
578   END SUBROUTINE trd_dwr
579
580
581   SUBROUTINE trd_twr( kt )
582      !!---------------------------------------------------------------------
583      !!                  ***  ROUTINE trd_twr  ***
584      !!
585      !! ** Purpose :  write active tracers trends in ocean.output
586      !!----------------------------------------------------------------------
587      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
588      !!----------------------------------------------------------------------
589
590      ! I. Tracers trends
591      ! -----------------
592
593      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
594
595         ! I.1 Sums over the global domain
596         ! -------------------------------
597         IF( lk_mpp ) THEN
598            CALL mpp_sum( tmo, jptot_tra )   
599            CALL mpp_sum( smo, jptot_tra )
600            CALL mpp_sum( t2 , jptot_tra )
601            CALL mpp_sum( s2 , jptot_tra )
602         ENDIF
603
604         ! I.2 Print tracers trends in the ocean.output file
605         ! --------------------------------------------------
606         
607         IF(lwp) THEN
608            WRITE (numout,*)
609            WRITE (numout,*)
610            WRITE (numout,9400) kt
611            WRITE (numout,9401)  tmo(jpicpt_xad) / tvolt, smo(jpicpt_xad) / tvolt
612            WRITE (numout,9411)  tmo(jpicpt_yad) / tvolt, smo(jpicpt_yad) / tvolt
613            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
614            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
615            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
616            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
617            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
618            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
619            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
620            WRITE (numout,9409) 
621            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
622            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
623            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
624            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
625         ENDIF
626
6279400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
628              '              salinity',/' ============================')
6299401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
6309411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
6319402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
6329403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
6339404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
6349405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
6359406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
6369407     FORMAT(' penetrative qsr             ',e20.13)
6379408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
6389409     FORMAT(' -------------------------------------------------------------------------')
6399410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
640
641
642         IF(lwp) THEN
643            WRITE (numout,*)
644            WRITE (numout,*)
645            WRITE (numout,9420) kt
646            WRITE (numout,9421)   t2(jpicpt_xad) / tvolt, s2(jpicpt_xad) / tvolt
647            WRITE (numout,9431)   t2(jpicpt_yad) / tvolt, s2(jpicpt_yad) / tvolt
648            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
649            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
650            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
651            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
652            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
653            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
654            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
655            WRITE (numout,9429)
656            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
657            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
658            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
659            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
660         ENDIF
661
6629420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
663            '               salinity', /, ' ===============================')
6649421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
6659431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
6669422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
6679423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
6689424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
6699425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
6709426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
6719427     FORMAT(' penetrative qsr           * t   ', e20.13)
6729428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
6739429     FORMAT(' -----------------------------------------------------------------------------')
6749430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
675
676
677         IF(lwp) THEN
678            WRITE (numout,*)
679            WRITE (numout,*)
680            WRITE (numout,9440) kt
681            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
682            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
683            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
684            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
685            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
686            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
687            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
688            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
689            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
690            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
691            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
692         ENDIF
693
6949440     FORMAT(' tracer consistency at it= ',i6,   &
695            ' :         temperature','                salinity',/,   &
696            ' ==================================')
6979441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
6989442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
6999443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
7009444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
7019445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
7029446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
7039447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
7049448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
7059449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
706         !
707      ENDIF
708      !
709   END SUBROUTINE trd_twr
710
711#   else
712   !!----------------------------------------------------------------------
713   !!   Default case :                                         Empty module
714   !!----------------------------------------------------------------------
715   INTERFACE trd_icp
716      MODULE PROCEDURE trd_2d, trd_3d
717   END INTERFACE
718
719CONTAINS
720   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )       ! Empty routine
721      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
722      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
723      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
724      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
725          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
726   END SUBROUTINE trd_2d
727   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype )       ! Empty routine
728      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
729      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
730      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
731      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
732          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
733   END SUBROUTINE trd_3d
734   SUBROUTINE trd_icp_init               ! Empty routine
735   END SUBROUTINE trd_icp_init
736   SUBROUTINE trd_dwr( kt )          ! Empty routine
737      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
738   END SUBROUTINE trd_dwr
739   SUBROUTINE trd_twr( kt )          ! Empty routine
740      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
741   END SUBROUTINE trd_twr
742#endif
743
744   !!======================================================================
745END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.