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

source: trunk/NEMO/OPA_SRC/TRD/trdicp.F90 @ 703

Last change on this file since 703 was 699, checked in by smasson, 17 years ago

insert revision Id
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.4 KB
Line 
1MODULE trdicp
2   !!======================================================================
3   !!                       ***  MODULE  trdicp  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic trends
5   !!=====================================================================
6   !! History :       !  91-12 (G. Madec)
7   !!                 !  92-06 (M. Imbard) add time step frequency
8   !!                 !  96-01 (G. Madec)  terrain following coordinates
9   !!            8.5  !  02-06 (G. Madec)  F90: Free form and module
10   !!            9.0  !  04-08 (C. Talandier) New trends organization
11   !!----------------------------------------------------------------------
12#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_trdtra'  or                  active tracers trends diagnostics
15   !!   'key_trddyn'                            momentum trends diagnostics
16   !!----------------------------------------------------------------------
17   !!----------------------------------------------------------------------
18   !!   trd_icp          : compute the basin averaged properties for tra/dyn
19   !!   trd_dwr          : print dynmaic trends in ocean.output file
20   !!   trd_twr          : print tracers trends in ocean.output file
21   !!   trd_icp_init     : initialization step
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain variables
25   USE trdmod_oce      ! ocean variables trends
26   USE ldftra_oce      ! ocean active tracers: lateral physics
27   USE ldfdyn_oce      ! ocean dynamics: lateral physics
28   USE zdf_oce         ! ocean vertical physics
29   USE in_out_manager  ! I/O manager
30   USE lib_mpp         ! distibuted memory computing library
31   USE eosbn2          ! equation of state
32   USE phycst          ! physical constants
33
34   IMPLICIT NONE
35   PRIVATE
36
37   INTERFACE trd_icp
38      MODULE PROCEDURE trd_2d, trd_3d
39   END INTERFACE
40
41   PUBLIC   trd_icp       ! called by trdmod.F90
42   PUBLIC   trd_dwr       ! called by step.F90
43   PUBLIC   trd_twr       ! called by step.F90
44   PUBLIC   trd_icp_init  ! called by opa.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !!   OPA 9.0 , LOCEAN-IPSL (2005)
51   !! $Id$
52   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )
58      !!---------------------------------------------------------------------
59      !!                  ***  ROUTINE trd_2d  ***
60      !!
61      !! ** Purpose : verify the basin averaged properties of tracers and/or
62      !!              momentum equations at every time step frequency ntrd.
63      !!----------------------------------------------------------------------
64      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx             ! Temperature or U trend
65      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy             ! Salinity    or V trend
66      INTEGER                     , INTENT(in   ) ::   ktrd                ! tracer trend index
67      CHARACTER(len=3)            , INTENT(in   ) ::   ctype               ! momentum ('DYN') or tracers ('TRA') trends
68      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
69      !!
70      INTEGER  ::   ji, jj                                                 ! loop indices
71      CHARACTER(len=3) ::   cpas                                           ! number of passage
72      REAL(wp) ::   zmsku, zbtu, zbt                                       ! temporary scalars
73      REAL(wp) ::   zmskv, zbtv                                            !    "         "
74      !!----------------------------------------------------------------------
75
76      ! Control of optional arguments
77      cpas = 'fst'
78      IF( PRESENT(clpas) )  cpas = clpas
79
80      ! 1. Mask trends
81      ! --------------
82
83      SELECT CASE( ctype )
84      !
85      CASE( 'DYN' )              ! Momentum
86         DO jj = 1, jpjm1
87            DO ji = 1, jpim1
88               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
89               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
90               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * zmsku
91               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * zmskv
92            END DO
93         END DO
94         ptrd2dx(jpi, : ) = 0.e0      ;      ptrd2dy(jpi, : ) = 0.e0
95         ptrd2dx( : ,jpj) = 0.e0      ;      ptrd2dy( : ,jpj) = 0.e0
96         !
97      CASE( 'TRA' )              ! Tracers
98         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
99         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
100         !
101      END SELECT
102     
103      ! 2. Basin averaged tracer/momentum trends
104      ! ----------------------------------------
105
106      SELECT CASE( ctype )
107      !
108      CASE( 'DYN' )              ! Momentum
109         umo(ktrd) = 0.e0
110         vmo(ktrd) = 0.e0
111         !
112         SELECT CASE( ktrd )
113         !
114         CASE( jpdyn_trd_swf )         ! surface forcing
115            DO jj = 1, jpj
116               DO ji = 1, jpi
117                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
118                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
119               END DO
120            END DO
121            !
122         CASE( jpdyn_trd_bfr )         ! bottom friction fluxes
123            DO jj = 1, jpj
124               DO ji = 1, jpi
125                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj)
126                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj)
127               END DO
128            END DO
129            !
130         END SELECT
131         !
132      CASE( 'TRA' )              ! Tracers
133         IF( cpas == 'fst' )   THEN
134            tmo(ktrd) = 0.e0
135            smo(ktrd) = 0.e0
136         ENDIF
137         DO jj = 1, jpj
138            DO ji = 1, jpi
139               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
140               tmo(ktrd) =  tmo(ktrd) + ptrd2dx(ji,jj) * zbt
141               smo(ktrd) =  smo(ktrd) + ptrd2dy(ji,jj) * zbt
142            END DO
143         END DO
144         !
145      END SELECT
146     
147      ! 3. Basin averaged tracer/momentum square trends
148      ! ----------------------------------------------
149      ! c a u t i o n: field now
150     
151      SELECT CASE( ctype )
152      !
153      CASE( 'DYN' )              ! Momentum
154         hke(ktrd) = 0.e0
155         DO jj = 1, jpj
156            DO ji = 1, jpi
157               zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
158               zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
159               hke(ktrd) = hke(ktrd)   &
160               &   + un(ji,jj,1) * ptrd2dx(ji,jj) * zbtu &
161               &   + vn(ji,jj,1) * ptrd2dy(ji,jj) * zbtv
162            END DO
163         END DO
164         !
165      CASE( 'TRA' )              ! Tracers
166         IF( cpas == 'fst' )   THEN
167            t2(ktrd) = 0.e0
168            s2(ktrd) = 0.e0
169         ENDIF
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
173               t2(ktrd) = t2(ktrd) + ptrd2dx(ji,jj) * zbt * tn(ji,jj,1)
174               s2(ktrd) = s2(ktrd) + ptrd2dy(ji,jj) * zbt * sn(ji,jj,1)
175            END DO
176         END DO
177         !     
178      END SELECT
179      !
180   END SUBROUTINE trd_2d
181
182
183   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype, clpas )
184      !!---------------------------------------------------------------------
185      !!                  ***  ROUTINE trd_3d  ***
186      !!
187      !! ** Purpose : verify the basin averaged properties of tracers and/or
188      !!              momentum equations at every time step frequency ntrd.
189      !!----------------------------------------------------------------------
190      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx            ! Temperature or U trend
191      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy            ! Salinity    or V trend
192      INTEGER,                          INTENT(in   ) ::   ktrd               ! momentum or tracer trend index
193      CHARACTER(len=3),                 INTENT(in   ) ::   ctype              ! momentum ('DYN') or tracers ('TRA') trends
194      CHARACTER(len=3),                 INTENT(in   ), OPTIONAL ::   clpas    ! number of passage
195      !!
196      INTEGER ::   ji, jj, jk
197      CHARACTER(len=3) ::   cpas                                              ! number of passage
198      REAL(wp) ::   zbt, zbtu, zbtv, zmsku, zmskv                             ! temporary scalars
199      !!----------------------------------------------------------------------
200
201      ! Control of optional arguments
202      cpas = 'fst'
203      IF( PRESENT(clpas) )  cpas = clpas
204
205      ! 1. Mask the trends
206      ! ------------------
207
208      SELECT CASE( ctype )
209      !
210      CASE( 'DYN' )              ! Momentum       
211         DO jk = 1, jpk
212            DO jj = 1, jpjm1
213               DO ji = 1, jpim1
214                  zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
215                  zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
216                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * zmsku
217                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * zmskv
218               END DO
219            END DO
220         END DO
221         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
222         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
223         !
224      CASE( 'TRA' )              ! Tracers
225         DO jk = 1, jpk
226            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
227            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
228         END DO
229         !
230      END SELECT   
231
232      ! 2. Basin averaged tracer/momentum trends
233      ! ----------------------------------------
234     
235      SELECT CASE( ctype )
236      !
237      CASE( 'DYN' )              ! Momentum
238         umo(ktrd) = 0.e0
239         vmo(ktrd) = 0.e0
240         DO jk = 1, jpk
241            DO jj = 1, jpj
242               DO ji = 1, jpi
243                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
244                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
245                  umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * zbtu
246                  vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * zbtv
247               END DO
248            END DO
249         END DO
250         !
251      CASE( 'TRA' )              ! Tracers
252         IF( cpas == 'fst' )   THEN
253            tmo(ktrd) = 0.e0
254            smo(ktrd) = 0.e0
255         ENDIF
256         DO jk = 1, jpkm1
257            DO jj = 1, jpj
258               DO ji = 1, jpi
259                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
260                  tmo(ktrd) =  tmo(ktrd) + ptrd3dx(ji,jj,jk) * zbt
261                  smo(ktrd) =  smo(ktrd) + ptrd3dy(ji,jj,jk) * zbt
262               END DO
263            END DO
264         END DO
265         !
266      END SELECT
267
268      ! 3. Basin averaged tracer/momentum square trends
269      ! -----------------------------------------------
270      ! c a u t i o n: field now
271     
272      SELECT CASE( ctype )
273      !
274      CASE( 'DYN' )              ! Momentum
275         hke(ktrd) = 0.e0
276         DO jk = 1, jpk
277            DO jj = 1, jpj
278               DO ji = 1, jpi
279                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
280                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
281                  hke(ktrd) = hke(ktrd)   &
282                  &   + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * zbtu &
283                  &   + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * zbtv
284               END DO
285            END DO
286         END DO
287         !
288      CASE( 'TRA' )              ! Tracers
289         IF( cpas == 'fst' )   THEN
290            t2(ktrd) = 0.e0
291            s2(ktrd) = 0.e0
292         ENDIF
293         DO jk = 1, jpk
294            DO jj = 1, jpj
295               DO ji = 1, jpi
296                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
297                  t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * zbt * tn(ji,jj,jk)
298                  s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * zbt * sn(ji,jj,jk)
299               END DO
300            END DO
301         END DO
302         !
303      END SELECT
304      !
305   END SUBROUTINE trd_3d
306
307
308
309   SUBROUTINE trd_icp_init
310      !!---------------------------------------------------------------------
311      !!                  ***  ROUTINE trd_icp_init  ***
312      !!
313      !! ** Purpose :   Read the namtrd namelist
314      !!----------------------------------------------------------------------
315      INTEGER  ::   ji, jj, jk
316      REAL(wp) ::   zmskt
317#if  defined key_trddyn
318      REAL(wp) ::   zmsku, zmskv
319#endif
320      !!----------------------------------------------------------------------
321
322      IF(lwp) THEN
323         WRITE(numout,*)
324         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
325         WRITE(numout,*) '~~~~~~~~~~~~~'
326      ENDIF
327
328      ! Total volume at t-points:
329      tvolt = 0.e0
330      DO jk = 1, jpkm1
331         DO jj = 2, jpjm1
332            DO ji = fs_2, fs_jpim1   ! vector opt.
333               zmskt = tmask(ji,jj,jk) * tmask_i(ji,jj)
334               tvolt = tvolt + zmskt * e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)
335            END DO
336         END DO
337      END DO
338      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
339
340      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
341
342#if  defined key_trddyn
343      ! Initialization of potential to kinetic energy conversion
344      rpktrd = 0.e0
345
346      ! Total volume at u-, v- points:
347      tvolu = 0.e0
348      tvolv = 0.e0
349
350      DO jk = 1, jpk
351         DO jj = 2, jpjm1
352            DO ji = fs_2, fs_jpim1   ! vector opt.
353               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
354               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
355               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
356               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
357            END DO
358         END DO
359      END DO
360      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
361      IF( lk_mpp )   CALL mpp_sum( tvolv )
362
363      IF(lwp) THEN
364         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
365         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
366      ENDIF
367#endif
368      !
369   END SUBROUTINE trd_icp_init
370
371
372   SUBROUTINE trd_dwr( kt )
373      !!---------------------------------------------------------------------
374      !!                  ***  ROUTINE trd_dwr  ***
375      !!
376      !! ** Purpose :  write dynamic trends in ocean.output
377      !!----------------------------------------------------------------------
378      INTEGER, INTENT(in) ::   kt                                  ! ocean time-step index
379      !!
380      INTEGER  ::   ji, jj, jk
381      REAL(wp) ::   ze1e2w, zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
382      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkepe, zkx, zky, zkz   ! temporary arrays
383      !!----------------------------------------------------------------------
384
385      ! I. Momentum trends
386      ! -------------------
387
388      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
389
390         ! I.1 Conversion potential energy - kinetic energy
391         ! --------------------------------------------------
392         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
393
394         zkx(:,:,:)   = 0.e0
395         zky(:,:,:)   = 0.e0
396         zkz(:,:,:)   = 0.e0
397         zkepe(:,:,:) = 0.e0
398   
399         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
400
401         ! Density flux at w-point
402         DO jk = 2, jpk
403            DO jj = 1, jpj
404               DO ji = 1, jpi
405                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
406                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
407               END DO
408            END DO
409         END DO
410         zkz(:,:,1) = 0.e0
411         
412         ! Density flux at u and v-points
413         DO jk = 1, jpk
414            DO jj = 1, jpjm1
415               DO ji = 1, jpim1
416                  zcof   = 0.5 / rau0
417                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
418                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
419                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
420                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
421               END DO
422            END DO
423         END DO
424         
425         ! Density flux divergence at t-point
426         DO jk = 1, jpkm1
427            DO jj = 2, jpjm1
428               DO ji = 2, jpim1
429                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
430                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
431                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
432                    &             + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
433                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
434               END DO
435            END DO
436         END DO
437         zkepe( : , : ,jpk) = 0.e0
438         zkepe( : ,jpj, : ) = 0.e0
439         zkepe(jpi, : , : ) = 0.e0
440
441         ! I.2 Basin averaged kinetic energy trend
442         ! ----------------------------------------
443         peke = 0.e0
444         DO jk = 1,jpk
445            DO jj = 1, jpj
446               DO ji = 1, jpi
447                  peke = peke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
448                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
449               END DO
450            END DO
451         END DO
452
453         ! I.3 Sums over the global domain
454         ! ---------------------------------
455         IF( lk_mpp ) THEN
456            CALL mpp_sum( peke )
457            CALL mpp_sum( umo , jptot_dyn )
458            CALL mpp_sum( vmo , jptot_dyn )
459            CALL mpp_sum( hke , jptot_dyn )
460         ENDIF
461
462         ! I.2 Print dynamic trends in the ocean.output file
463         ! --------------------------------------------------
464
465         IF(lwp) THEN
466            WRITE (numout,*)
467            WRITE (numout,*)
468            WRITE (numout,9500) kt
469            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
470            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
471            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
472            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
473            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
474            WRITE (numout,9506) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
475            WRITE (numout,9507) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
476            WRITE (numout,9508) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
477            WRITE (numout,9509) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
478            WRITE (numout,9510) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
479            WRITE (numout,9511) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
480            WRITE (numout,9512)
481            WRITE (numout,9513)                                                 &
482            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
483            &      + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat) + umo(jpicpd_swf)   &
484            &      + umo(jpicpd_bfr) ) / tvolu,   &
485            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
486            &      + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat) + vmo(jpicpd_swf)   &
487            &      + 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(' vertical advection         u= ', e20.13, '    v= ', e20.13)
497 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
498 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
499 9509    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
500 9510    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
501 9511    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
502 9512    FORMAT(' -----------------------------------------------------------------------------')
503 9513    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
504
505         IF(lwp) THEN
506            WRITE (numout,*)
507            WRITE (numout,*)
508            WRITE (numout,9520) kt
509            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
510            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
511            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
512            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
513            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
514            WRITE (numout,9526) hke(jpicpd_zad) / tvolt
515            WRITE (numout,9527) hke(jpicpd_zdf) / tvolt
516            WRITE (numout,9528) hke(jpicpd_spg) / tvolt
517            WRITE (numout,9529) hke(jpicpd_swf) / tvolt
518            WRITE (numout,9530) hke(jpicpd_dat) / tvolt
519            WRITE (numout,9531)
520            WRITE (numout,9532)   &
521            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
522            &      + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat) + hke(jpicpd_swf) ) / tvolt
523         ENDIF
524
525 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
526 9521    FORMAT(' pressure gradient         u2= ', e20.13)
527 9522    FORMAT(' ke gradient               u2= ', e20.13)
528 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
529 9524    FORMAT(' coriolis term             u2= ', e20.13)
530 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
531 9526    FORMAT(' vertical advection        u2= ', e20.13)
532 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
533 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
534 9529    FORMAT(' surface wind forcing      u2= ', e20.13)
535 9530    FORMAT(' dampimg term              u2= ', e20.13)
536 9531    FORMAT(' --------------------------------------------------')
537 9532    FORMAT(' total trend               u2= ', e20.13)
538
539         IF(lwp) THEN
540            WRITE (numout,*)
541            WRITE (numout,*)
542            WRITE (numout,9540) kt
543            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_zad) ) / tvolt
544            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_zad) ) / tvolt
545            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
546            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
547            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
548            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
549            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
550            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
551            WRITE (numout,*)
552            WRITE (numout,*)
553         ENDIF
554
555 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
556 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
557 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
558 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
559 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
560 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
561 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
562 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
563 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
564         !
565         ! Save potential to kinetic energy conversion for next time step
566         rpktrd = peke
567         !
568      ENDIF
569      !
570   END SUBROUTINE trd_dwr
571
572
573   SUBROUTINE trd_twr( kt )
574      !!---------------------------------------------------------------------
575      !!                  ***  ROUTINE trd_twr  ***
576      !!
577      !! ** Purpose :  write active tracers trends in ocean.output
578      !!----------------------------------------------------------------------
579      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
580      !!----------------------------------------------------------------------
581
582      ! I. Tracers trends
583      ! -----------------
584
585      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
586
587         ! I.1 Sums over the global domain
588         ! -------------------------------
589         IF( lk_mpp ) THEN
590            CALL mpp_sum( tmo, jptot_tra )   
591            CALL mpp_sum( smo, jptot_tra )
592            CALL mpp_sum( t2 , jptot_tra )
593            CALL mpp_sum( s2 , jptot_tra )
594         ENDIF
595
596         ! I.2 Print tracers trends in the ocean.output file
597         ! --------------------------------------------------
598         
599         IF(lwp) THEN
600            WRITE (numout,*)
601            WRITE (numout,*)
602            WRITE (numout,9400) kt
603            WRITE (numout,9401) (tmo(jpicpt_xad)+tmo(jpicpt_yad))/ tvolt, (smo(jpicpt_xad)+smo(jpicpt_yad))/ tvolt
604            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
605            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
606            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
607            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
608            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
609            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
610            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
611            WRITE (numout,9409) 
612            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
613            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
614            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
615            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
616         ENDIF
617
6189400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
619              '              salinity',/' ============================')
6209401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
6219402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
6229403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
6239404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
6249405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
6259406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
6269407     FORMAT(' penetrative qsr             ',e20.13)
6279408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
6289409     FORMAT(' -------------------------------------------------------------------------')
6299410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
630
631
632         IF(lwp) THEN
633            WRITE (numout,*)
634            WRITE (numout,*)
635            WRITE (numout,9420) kt
636            WRITE (numout,9421) ( t2(jpicpt_xad)+t2(jpicpt_yad) )/ tvolt, ( s2(jpicpt_xad)+s2(jpicpt_yad) )/ tvolt
637            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
638            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
639            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
640            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
641            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
642            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
643            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
644            WRITE (numout,9429)
645            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
646            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
647            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
648            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
649         ENDIF
650
6519420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
652            '               salinity', /, ' ===============================')
6539421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
6549422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
6559423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
6569424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
6579425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
6589426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
6599427     FORMAT(' penetrative qsr           * t   ', e20.13)
6609428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
6619429     FORMAT(' -----------------------------------------------------------------------------')
6629430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
663
664
665         IF(lwp) THEN
666            WRITE (numout,*)
667            WRITE (numout,*)
668            WRITE (numout,9440) kt
669            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
670            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
671            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
672            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
673            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
674            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
675            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
676            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
677            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
678            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
679            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
680         ENDIF
681
6829440     FORMAT(' tracer consistency at it= ',i6,   &
683            ' :         temperature','                salinity',/,   &
684            ' ==================================')
6859441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
6869442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
6879443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
6889444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
6899445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
6909446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
6919447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
6929448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
6939449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
694         !
695      ENDIF
696      !
697   END SUBROUTINE trd_twr
698
699#   else
700   !!----------------------------------------------------------------------
701   !!   Default case :                                         Empty module
702   !!----------------------------------------------------------------------
703   INTERFACE trd_icp
704      MODULE PROCEDURE trd_2d, trd_3d
705   END INTERFACE
706
707CONTAINS
708   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )       ! Empty routine
709      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
710      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
711      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
712      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
713      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
714          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype, clpas
715   END SUBROUTINE trd_2d
716   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype, clpas )       ! Empty routine
717      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
718      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
719      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
720      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
721      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
722          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype, clpas
723   END SUBROUTINE trd_3d
724   SUBROUTINE trd_icp_init               ! Empty routine
725   END SUBROUTINE trd_icp_init
726   SUBROUTINE trd_dwr( kt )          ! Empty routine
727      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
728   END SUBROUTINE trd_dwr
729   SUBROUTINE trd_twr( kt )          ! Empty routine
730      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
731   END SUBROUTINE trd_twr
732#endif
733
734   !!======================================================================
735END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.