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

source: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 790

Last change on this file since 790 was 790, checked in by gm, 17 years ago

dev_001_GM - complete theprevious comit with omitted routines

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.9 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, trd_u2d, trd_u3d
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   !! Variables used for diagnostics
47   REAL(wp) ::   tvolt        !: volume of the whole ocean computed at t-points
48   REAL(wp) ::   tvolu        !: volume of the whole ocean computed at u-points
49   REAL(wp) ::   tvolv        !: volume of the whole ocean computed at v-points
50
51   !! Active Tracer trend diagnostics variables
52   REAL(wp), DIMENSION(jpt_trd,2) ::   tsmo         !: tracers trends average
53   REAL(wp), DIMENSION(jpt_trd,2) ::   ts2          !: tracers square trends average
54
55   !! Momentum trends diagnostics variables
56   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo         !: momentum trends average
57   REAL(wp), DIMENSION(jptot_dyn) ::   hke              !: momentum square trends average
58   REAL(wp) ::   rpktrd   !: potential to kinetic energy conversion
59   REAL(wp) ::   peke     !: conversion potential energy - kinetic energy trend
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90"
63#  include "vectopt_loop_substitute.h90"
64   !!----------------------------------------------------------------------
65   !!   OPA 9.0 , LOCEAN-IPSL (2005)
66   !! $Header$
67   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
68   !!----------------------------------------------------------------------
69
70CONTAINS
71
72   SUBROUTINE trd_2d( ptrd, ktra, ktrd , ctype, clpas )
73      !!---------------------------------------------------------------------
74      !!                  ***  ROUTINE trd_2d  ***
75      !!
76      !! ** Purpose : verify the basin averaged properties of tracers and/or
77      !!              momentum equations at every time step frequency ntrd.
78      !!----------------------------------------------------------------------
79      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptrd      ! tracer or u trend
80      INTEGER                     , INTENT(in   )           ::   ktra      ! tracer index
81      INTEGER                     , INTENT(in   )           ::   ktrd      ! tracer trend index
82      CHARACTER(len=3)            , INTENT(in   )           ::   ctype     ! tracer type (='TRA' or 'TRC')
83      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
84      !!
85      CHARACTER(len=3) ::   cpas   ! number of passage
86      REAL(wp)         ::   zsum   ! temporary scalars
87      !!----------------------------------------------------------------------
88
89      cpas = 'fst'                             ! Control of optional arguments
90      IF( PRESENT(clpas) )  cpas = clpas
91
92      ! 1. Mask volumic trends
93      ptrd(:,:) = e1t(:,:) * e2t(:,:) * fse3t(:,:,1) * ptrd(:,:) * tmask_i(:,:)
94     
95      ! 2. Basin averaged tracer trends
96      SELECT CASE( ctype )
97      CASE( 'TRA' )              ! Tracers
98         zsum = SUM( ptrd(:,:) )
99         IF( cpas == 'fst' ) THEN   ;   tsmo(ktrd,ktra) =                   zsum
100         ELSE                       ;   tsmo(ktrd,ktra) = tsmo(ktrd,ktra) + zsum
101         ENDIF
102      CASE( 'TRC' )              ! Passive tracers
103         !      ....  to be done
104      END SELECT
105     
106      ! 3. Basin averaged tracer square trends (i.e. trd * now tracer field)
107      SELECT CASE( ctype )
108      CASE( 'TRA' )              ! Tracers
109         IF( ktra == jp_tem )   zsum = SUM(  ptrd(:,:) * tn(:,:,1) )
110         IF( ktra == jp_sal )   zsum = SUM(  ptrd(:,:) * sn(:,:,1) )
111         IF( cpas == 'fst' ) THEN   ;   ts2(ktrd,ktra) =                  zsum
112         ELSE                       ;   ts2(ktrd,ktra) = ts2(ktrd,ktra) + zsum
113         ENDIF
114      CASE( 'TRC' )              ! Passive tracers
115         !      ....  to be done
116      END SELECT
117      !
118   END SUBROUTINE trd_2d
119
120
121   SUBROUTINE trd_3d( ptrd, ktra, ktrd , ctype, clpas )
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 ntrd.
127      !!----------------------------------------------------------------------
128      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   ptrd      ! tracer or u trend
129      INTEGER                         , INTENT(in   )           ::   ktra      ! tracer index
130      INTEGER                         , INTENT(in   )           ::   ktrd      ! tracer trend index
131      CHARACTER(len=3)                , INTENT(in   )           ::   ctype     ! tracer type (='TRA' or 'TRC')
132      CHARACTER(len=3)                , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
133      !!
134      INTEGER          ::   jk     ! dummy loop indices
135      CHARACTER(len=3) ::   cpas   ! number of passage
136      REAL(wp)         ::   zsum   ! temporary scalars
137      REAL(wp), DIMENSION(jpi,jpj) ::   zsurf          ! 2D workspace array
138      !!----------------------------------------------------------------------
139
140      cpas = 'fst'                             ! Control of optional arguments
141      IF( PRESENT(clpas) )  cpas = clpas
142
143      ! 1. Mask volumic trends
144      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
145      DO jk = 1, jpk
146         ptrd(:,:,jk) = zsurf(:,:) * fse3t(:,:,jk) * ptrd(:,:,jk) * tmask(:,:,jk)
147      END DO
148
149      ! 2. Basin averaged tracer trends
150      SELECT CASE( ctype )
151      CASE( 'TRA' )              ! Tracers
152         zsum = SUM( ptrd(:,:,:) )
153         IF( cpas == 'fst' ) THEN   ;   tsmo(ktrd,ktra) =                   zsum
154         ELSE                       ;   tsmo(ktrd,ktra) = tsmo(ktrd,ktra) + zsum
155         ENDIF
156      CASE( 'TRC' )              ! Passive tracers
157         !      ....  to be done
158      END SELECT
159
160      ! 3. Basin averaged tracer square trends (i.e. trd * now tracer field)
161      SELECT CASE( ctype )
162      CASE( 'TRA' )              ! Tracers
163         IF( ktra == jp_tem )   zsum = SUM(  ptrd(:,:,:) * tn(:,:,:) )
164         IF( ktra == jp_sal )   zsum = SUM(  ptrd(:,:,:) * sn(:,:,:) )
165         IF( cpas == 'fst' ) THEN   ;   ts2(ktrd,ktra) =                  zsum
166         ELSE                       ;   ts2(ktrd,ktra) = ts2(ktrd,ktra) + zsum
167         ENDIF
168      CASE( 'TRC' )              ! Passive tracers
169         !      ....  to be done
170      END SELECT
171      !
172   END SUBROUTINE trd_3d
173
174
175   SUBROUTINE trd_u2d( ptrdu, ptrdv, ktrd , ctype, clpas )
176      !!---------------------------------------------------------------------
177      !!                  ***  ROUTINE trd_2d  ***
178      !!
179      !! ** Purpose : verify the basin averaged properties of tracers and/or
180      !!              momentum equations at every time step frequency ntrd.
181      !!----------------------------------------------------------------------
182      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)           ::   ptrdu, ptrdv   ! U and V momentum trends
183      INTEGER                     , INTENT(in   )           ::   ktrd           ! tracer trend index
184      CHARACTER(len=3)            , INTENT(in   )           ::   ctype          ! momentum ('DYN')
185      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas          ! number of passage
186      !!
187      INTEGER          ::   ji, jj        ! loop indices
188      CHARACTER(len=3) ::   cpas          ! number of passage
189      REAL(wp)         ::   zmsku, zbtu   ! temporary scalars
190      REAL(wp)         ::   zmskv, zbtv   !    "         "
191      !!----------------------------------------------------------------------
192
193      cpas = 'fst'                           ! Control of optional arguments
194      IF( PRESENT(clpas) )  cpas = clpas
195
196      SELECT CASE( ctype )
197      !
198      CASE( 'DYN' )              ! Momentum
199
200      ! 1. Mask trends
201      ! --------------
202      DO jj = 1, jpjm1
203         DO ji = 1, jpim1
204            zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
205            zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
206            ptrdu(ji,jj) = ptrdu(ji,jj) * zmsku
207            ptrdv(ji,jj) = ptrdv(ji,jj) * zmskv
208         END DO
209      END DO
210      ptrdu(jpi, : ) = 0.e0      ;      ptrdv(jpi, : ) = 0.e0
211      ptrdu( : ,jpj) = 0.e0      ;      ptrdv( : ,jpj) = 0.e0
212     
213      ! 2. Basin averaged tracer/momentum trends
214      ! ----------------------------------------
215      umo(ktrd) = 0.e0
216      vmo(ktrd) = 0.e0
217      !
218      SELECT CASE( ktrd )
219      !
220      CASE( jpdyn_trd_swf )         ! surface forcing
221         DO jj = 1, jpj
222            DO ji = 1, jpi
223               umo(ktrd) = umo(ktrd) + ptrdu(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
224               vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
225            END DO
226         END DO
227            !
228         CASE( jpdyn_trd_bfr )         ! bottom friction fluxes
229            DO jj = 1, jpj
230               DO ji = 1, jpi
231                  umo(ktrd) = umo(ktrd) + ptrdu(ji,jj)
232                  vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj)
233               END DO
234            END DO
235            !
236      END SELECT
237      !   
238      ! 3. Basin averaged tracer/momentum square trends
239      ! ----------------------------------------------
240      ! c a u t i o n: field now
241     
242         hke(ktrd) = 0.e0
243         DO jj = 1, jpj
244            DO ji = 1, jpi
245               zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
246               zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
247               hke(ktrd) = hke(ktrd) + un(ji,jj,1) * ptrdu(ji,jj) * zbtu   &
248               &                     + vn(ji,jj,1) * ptrdv(ji,jj) * zbtv
249            END DO
250         END DO
251         !
252      END SELECT
253      !
254   END SUBROUTINE trd_u2d
255
256
257   SUBROUTINE trd_u3d( ptrdu, ptrdv, ktrd, ctype, clpas )
258      !!---------------------------------------------------------------------
259      !!                  ***  ROUTINE trd_3d  ***
260      !!
261      !! ** Purpose : verify the basin averaged properties of tracers and/or
262      !!              momentum equations at every time step frequency ntrd.
263      !!----------------------------------------------------------------------
264      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::   ptrdu, ptrdv   ! U and V momentum trends
265      INTEGER,                          INTENT(in   )           ::   ktrd           ! momentum trend index
266      CHARACTER(len=3),                 INTENT(in   )           ::   ctype          ! momentum (='DYN')
267      CHARACTER(len=3),                 INTENT(in   ), OPTIONAL ::   clpas          ! number of passage
268      !!
269      INTEGER ::   ji, jj, jk
270      CHARACTER(len=3) ::   cpas               ! number of passage
271      REAL(wp) ::   zbtu, zbtv, zmsku, zmskv   ! temporary scalars
272      !!----------------------------------------------------------------------
273
274      cpas = 'fst'                            ! Control of optional arguments
275      IF( PRESENT(clpas) )  cpas = clpas
276
277      SELECT CASE( ctype )
278      !
279      CASE( 'DYN' )              ! Momentum       
280
281         ! 1. Mask the trends
282         ! ------------------
283         DO jk = 1, jpk
284            DO jj = 1, jpjm1
285               DO ji = 1, jpim1
286                  zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
287                  zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
288                  ptrdu(ji,jj,jk) = ptrdu(ji,jj,jk) * zmsku
289                  ptrdv(ji,jj,jk) = ptrdv(ji,jj,jk) * zmskv
290               END DO
291            END DO
292         END DO
293         ptrdu(jpi, : ,:) = 0.e0      ;      ptrdv(jpi, : ,:) = 0.e0
294         ptrdu( : ,jpj,:) = 0.e0      ;      ptrdv( : ,jpj,:) = 0.e0
295         !
296
297      ! 2. Basin averaged tracer/momentum trends
298      ! ----------------------------------------
299     
300         umo(ktrd) = 0.e0
301         vmo(ktrd) = 0.e0
302         DO jk = 1, jpk
303            DO jj = 1, jpj
304               DO ji = 1, jpi
305                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
306                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
307                  umo(ktrd) = umo(ktrd) + ptrdu(ji,jj,jk) * zbtu
308                  vmo(ktrd) = vmo(ktrd) + ptrdv(ji,jj,jk) * zbtv
309               END DO
310            END DO
311         END DO
312         !
313
314      ! 3. Basin averaged tracer/momentum square trends
315      ! -----------------------------------------------
316      ! c a u t i o n: field now
317     
318         hke(ktrd) = 0.e0
319         DO jk = 1, jpk
320            DO jj = 1, jpj
321               DO ji = 1, jpi
322                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
323                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
324                  hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * ptrdu(ji,jj,jk) * zbtu   &
325                  &                     + vn(ji,jj,jk) * ptrdv(ji,jj,jk) * zbtv
326               END DO
327            END DO
328         END DO
329         !
330      END SELECT
331      !
332   END SUBROUTINE trd_u3d
333
334
335   SUBROUTINE trd_icp_init
336      !!---------------------------------------------------------------------
337      !!                  ***  ROUTINE trd_icp_init  ***
338      !!
339      !! ** Purpose :   Read the namtrd namelist
340      !!----------------------------------------------------------------------
341      INTEGER  ::   jk
342#if  defined key_trddyn
343      INTEGER  ::   ji, jj
344      REAL(wp) ::   zmsku, zmskv
345#endif
346      REAL(wp), DIMENSION(jpi,jpj) ::   zsurf          ! 2D workspace array
347      !!----------------------------------------------------------------------
348
349      IF(lwp) THEN
350         WRITE(numout,*)
351         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
352         WRITE(numout,*) '~~~~~~~~~~~~~'
353      ENDIF
354
355      ! Total volume at t-points:
356      tvolt = 0.e0
357      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)
358      DO jk = 1, jpkm1
359            tvolt = tvolt + SUM( tmask(:,:,jk) * zsurf(:,:) * fse3t(:,:,jk) )
360      END DO
361      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
362
363      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
364
365#if  defined key_trddyn
366      ! Initialization of potential to kinetic energy conversion
367      rpktrd = 0.e0
368
369      ! Total volume at u-, v- points:
370      tvolu = 0.e0
371      tvolv = 0.e0
372
373      DO jk = 1, jpk
374         DO jj = 2, jpjm1
375            DO ji = fs_2, fs_jpim1   ! vector opt.
376               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
377               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
378               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
379               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
380            END DO
381         END DO
382      END DO
383      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
384      IF( lk_mpp )   CALL mpp_sum( tvolv )
385
386      IF(lwp) THEN
387         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
388         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
389      ENDIF
390#endif
391      !
392   END SUBROUTINE trd_icp_init
393
394
395   SUBROUTINE trd_dwr( kt )
396      !!---------------------------------------------------------------------
397      !!                  ***  ROUTINE trd_dwr  ***
398      !!
399      !! ** Purpose :  write dynamic trends in ocean.output
400      !!----------------------------------------------------------------------
401      INTEGER, INTENT(in) ::   kt                                  ! ocean time-step index
402      !!
403      INTEGER  ::   ji, jj, jk
404      REAL(wp) ::   zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
405      REAL(wp), DIMENSION(jpi,jpj)     ::   zsurf          ! 2D workspace array
406      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkepe, zkx, zky, zkz   ! temporary arrays
407      !!----------------------------------------------------------------------
408
409      ! I. Momentum trends
410      ! -------------------
411
412      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
413
414         ! I.1 Conversion potential energy - kinetic energy
415         ! --------------------------------------------------
416         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
417
418         zkx(:,:,:)   = 0.e0
419         zky(:,:,:)   = 0.e0
420         zkz(:,:,:)   = 0.e0
421         zkepe(:,:,:) = 0.e0
422   
423         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
424
425         ! Density flux at w-point
426         zsurf(:,:) = 0.5 * e1t(:,:) * e2t(:,:) * tmask_i(:,:) /rau0
427!!gm better use bn2....
428         zkz(:,:,1) = 0.e0
429         DO jk = 2, jpk
430            zkz(:,:,jk) = zsurf(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) )
431         END DO
432         
433         ! Density flux at u and v-points
434         DO jk = 1, jpk
435            DO jj = 1, jpjm1
436               DO ji = 1, jpim1
437                  zcof   = 0.5 / rau0
438                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
439                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
440                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
441                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
442               END DO
443            END DO
444         END DO
445         
446         ! Density flux divergence at t-point
447         DO jk = 1, jpkm1
448            DO jj = 2, jpjm1
449               DO ji = 2, jpim1
450                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
451                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
452                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
453                    &             + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
454                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
455               END DO
456            END DO
457         END DO
458         zkepe( : , : ,jpk) = 0.e0
459         zkepe( : ,jpj, : ) = 0.e0
460         zkepe(jpi, : , : ) = 0.e0
461
462         ! I.2 Basin averaged kinetic energy trend
463         ! ----------------------------------------
464         peke = 0.e0
465         DO jk = 1,jpk
466            DO jj = 1, jpj
467               DO ji = 1, jpi
468                  peke = peke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
469                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
470               END DO
471            END DO
472         END DO
473
474         ! I.3 Sums over the global domain
475         ! ---------------------------------
476         IF( lk_mpp ) THEN
477            CALL mpp_sum( peke )
478            CALL mpp_sum( umo , jptot_dyn )
479            CALL mpp_sum( vmo , jptot_dyn )
480            CALL mpp_sum( hke , jptot_dyn )
481         ENDIF
482
483         ! I.2 Print dynamic trends in the ocean.output file
484         ! --------------------------------------------------
485
486         IF(lwp) THEN
487            WRITE (numout,*)
488            WRITE (numout,*)
489            WRITE (numout,9500) kt
490            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
491            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
492            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
493            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
494            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
495            WRITE (numout,9506) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
496            WRITE (numout,9507) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
497            WRITE (numout,9508) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
498            WRITE (numout,9509) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
499            WRITE (numout,9510) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
500            WRITE (numout,9511) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
501            WRITE (numout,9512)
502            WRITE (numout,9513)                                                 &
503            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
504            &      + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat) + umo(jpicpd_swf)   &
505            &      + umo(jpicpd_bfr) ) / tvolu,   &
506            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
507            &      + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat) + vmo(jpicpd_swf)   &
508            &      + vmo(jpicpd_bfr) ) / tvolv
509         ENDIF
510
511 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
512 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
513 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
514 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
515 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
516 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
517 9506    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
518 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
519 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
520 9509    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
521 9510    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
522 9511    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
523 9512    FORMAT(' -----------------------------------------------------------------------------')
524 9513    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
525
526         IF(lwp) THEN
527            WRITE (numout,*)
528            WRITE (numout,*)
529            WRITE (numout,9520) kt
530            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
531            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
532            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
533            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
534            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
535            WRITE (numout,9526) hke(jpicpd_zad) / tvolt
536            WRITE (numout,9527) hke(jpicpd_zdf) / tvolt
537            WRITE (numout,9528) hke(jpicpd_spg) / tvolt
538            WRITE (numout,9529) hke(jpicpd_swf) / tvolt
539            WRITE (numout,9530) hke(jpicpd_dat) / tvolt
540            WRITE (numout,9531)
541            WRITE (numout,9532)   &
542            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
543            &      + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat) + hke(jpicpd_swf) ) / tvolt
544         ENDIF
545
546 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
547 9521    FORMAT(' pressure gradient         u2= ', e20.13)
548 9522    FORMAT(' ke gradient               u2= ', e20.13)
549 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
550 9524    FORMAT(' coriolis term             u2= ', e20.13)
551 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
552 9526    FORMAT(' vertical advection        u2= ', e20.13)
553 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
554 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
555 9529    FORMAT(' surface wind forcing      u2= ', e20.13)
556 9530    FORMAT(' dampimg term              u2= ', e20.13)
557 9531    FORMAT(' --------------------------------------------------')
558 9532    FORMAT(' total trend               u2= ', e20.13)
559
560         IF(lwp) THEN
561            WRITE (numout,*)
562            WRITE (numout,*)
563            WRITE (numout,9540) kt
564            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_zad) ) / tvolt
565            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_zad) ) / tvolt
566            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
567            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
568            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
569            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
570            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
571            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
572            WRITE (numout,*)
573            WRITE (numout,*)
574         ENDIF
575
576 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
577 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
578 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
579 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
580 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
581 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
582 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
583 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
584 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
585         !
586         ! Save potential to kinetic energy conversion for next time step
587         rpktrd = peke
588         !
589      ENDIF
590      !
591   END SUBROUTINE trd_dwr
592
593
594   SUBROUTINE trd_twr( kt )
595      !!---------------------------------------------------------------------
596      !!                  ***  ROUTINE trd_twr  ***
597      !!
598      !! ** Purpose :  write active tracers trends in ocean.output
599      !!----------------------------------------------------------------------
600      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
601      !!----------------------------------------------------------------------
602
603      ! I. Tracers trends
604      ! -----------------
605
606      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
607
608         ! I.1 Sums over the global domain
609         ! -------------------------------
610         IF( lk_mpp ) THEN
611            CALL mpp_sum( tsmo, jpt_trd )   
612            CALL mpp_sum( ts2 , jpt_trd )
613         ENDIF
614
615         ! I.2 Print tracers trends in the ocean.output file
616         ! --------------------------------------------------
617         
618         IF(lwp) THEN
619            WRITE (numout,*)
620            WRITE (numout,*)
621            WRITE (numout,9400) kt
622            WRITE (numout,9401) ( tsmo(jpt_trd_xad,jp_tem) + tsmo(jpt_trd_yad,jp_tem)  ) / tvolt,   &
623               &                ( tsmo(jpt_trd_xad,jp_sal) + tsmo(jpt_trd_yad,jp_sal)  ) / tvolt
624            WRITE (numout,9402)   tsmo(jpt_trd_zad,jp_tem) / tvolt,   tsmo(jpt_trd_zad,jp_sal) / tvolt
625            WRITE (numout,9403)   tsmo(jpt_trd_ldf,jp_tem) / tvolt,   tsmo(jpt_trd_ldf,jp_sal) / tvolt
626            WRITE (numout,9404)   tsmo(jpt_trd_zdf,jp_tem) / tvolt,   tsmo(jpt_trd_zdf,jp_sal) / tvolt
627            WRITE (numout,9405)   tsmo(jpt_trd_npc,jp_tem) / tvolt,   tsmo(jpt_trd_npc,jp_sal) / tvolt
628            WRITE (numout,9406)   tsmo(jpt_trd_dmp,jp_tem) / tvolt,   tsmo(jpt_trd_dmp,jp_sal) / tvolt
629            WRITE (numout,9407)   tsmo(jpt_trd_qsr,jp_tem) / tvolt
630            WRITE (numout,9408)   tsmo(jpt_trd_qns,jp_tem) / tvolt,   tsmo(jpt_trd_qns,jp_sal) / tvolt
631            WRITE (numout,9409) 
632            WRITE (numout,9410) (  tsmo(jpt_trd_xad,jp_tem) + tsmo(jpt_trd_yad,jp_tem) + tsmo(jpt_trd_zad,jp_tem)     &
633               &                 + tsmo(jpt_trd_ldf,jp_tem) + tsmo(jpt_trd_zdf,jp_tem) + tsmo(jpt_trd_npc,jp_tem)     &
634               &                 + tsmo(jpt_trd_dmp,jp_tem) + tsmo(jpt_trd_qsr,jp_tem) + tsmo(jpt_trd_qns,jp_tem) )   &
635               &                / tvolt,   &
636               &                (  tsmo(jpt_trd_xad,jp_sal) + tsmo(jpt_trd_yad,jp_sal) + tsmo(jpt_trd_zad,jp_sal)     &
637               &                 + tsmo(jpt_trd_ldf,jp_sal) + tsmo(jpt_trd_zdf,jp_sal) + tsmo(jpt_trd_npc,jp_sal)     &
638               &                 + tsmo(jpt_trd_dmp,jp_sal)                           + tsmo(jpt_trd_qns,jp_sal) )   &
639               &                / tvolt
640
641
6429400     FORMAT(' tracer trend at it= ',i6,' :     temperature              salinity', /,   &
643                ' ============================')
6449401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
6459402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
6469403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
6479404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
6489405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
6499406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
6509407     FORMAT(' penetrative qsr             ',e20.13)
6519408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
6529409     FORMAT(' -------------------------------------------------------------------------')
6539410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
654
655
656            WRITE (numout,*)
657            WRITE (numout,*)
658            WRITE (numout,9420) kt
659            WRITE (numout,9421) ( ts2(jpt_trd_xad,jp_tem) + ts2(jpt_trd_yad,jp_tem) ) / tvolt,   &
660               &                ( ts2(jpt_trd_xad,jp_sal) + ts2(jpt_trd_yad,jp_sal) ) / tvolt
661            WRITE (numout,9422)   ts2(jpt_trd_zad,jp_tem) / tvolt, ts2(jpt_trd_zad,jp_sal) / tvolt
662            WRITE (numout,9423)   ts2(jpt_trd_ldf,jp_tem) / tvolt, ts2(jpt_trd_ldf,jp_sal) / tvolt
663            WRITE (numout,9424)   ts2(jpt_trd_zdf,jp_tem) / tvolt, ts2(jpt_trd_zdf,jp_sal) / tvolt
664            WRITE (numout,9425)   ts2(jpt_trd_npc,jp_tem) / tvolt, ts2(jpt_trd_npc,jp_sal) / tvolt
665            WRITE (numout,9426)   ts2(jpt_trd_dmp,jp_tem) / tvolt, ts2(jpt_trd_dmp,jp_sal) / tvolt
666            WRITE (numout,9427)   ts2(jpt_trd_qsr,jp_tem) / tvolt
667            WRITE (numout,9428)   ts2(jpt_trd_qns,jp_tem) / tvolt, ts2(jpt_trd_qns,jp_sal) / tvolt
668            WRITE (numout,9429)
669            WRITE (numout,9430) (  ts2(jpt_trd_xad,jp_tem) + ts2(jpt_trd_yad,jp_tem) + ts2(jpt_trd_zad,jp_tem)     &
670               &                 + ts2(jpt_trd_ldf,jp_tem) + ts2(jpt_trd_zdf,jp_tem) + ts2(jpt_trd_npc,jp_tem)     &
671               &                 + ts2(jpt_trd_dmp,jp_tem) + ts2(jpt_trd_qsr,jp_tem) + ts2(jpt_trd_qns,jp_tem) )   &
672               &                / tvolt,   &
673               &                (  ts2(jpt_trd_xad,jp_sal) + ts2(jpt_trd_yad,jp_sal) + ts2(jpt_trd_zad,jp_sal)     &
674               &                 + ts2(jpt_trd_ldf,jp_sal) + ts2(jpt_trd_zdf,jp_sal) + ts2(jpt_trd_npc,jp_sal)     &
675               &                 + ts2(jpt_trd_dmp,jp_sal)                          + ts2(jpt_trd_qns,jp_sal) )   &
676               &                / tvolt
677
6789420     FORMAT(' tracer**2 trend at it= ',i6,' :     temperature              salinity', /,   &
679                ' ============================')
6809421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
6819422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
6829423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
6839424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
6849425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
6859426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
6869427     FORMAT(' penetrative qsr           * t   ', e20.13)
6879428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
6889429     FORMAT(' -----------------------------------------------------------------------------')
6899430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
690
691
692            WRITE (numout,*)
693            WRITE (numout,*)
694            WRITE (numout,9440) kt
695            WRITE (numout,9441) ( tsmo(jpt_trd_xad,jp_tem)+tsmo(jpt_trd_yad,jp_tem)+tsmo(jpt_trd_zad,jp_tem) )/tvolt,  &
696            &                   ( tsmo(jpt_trd_xad,jp_sal)+tsmo(jpt_trd_yad,jp_sal)+tsmo(jpt_trd_zad,jp_sal) )/tvolt
697            WRITE (numout,9442)   tsmo(jpt_trd_zl1,jp_tem)/tvolt,  tsmo(jpt_trd_zl1,jp_sal)/tvolt
698            WRITE (numout,9443)   tsmo(jpt_trd_ldf,jp_tem)/tvolt,  tsmo(jpt_trd_ldf,jp_sal)/tvolt
699            WRITE (numout,9444)   tsmo(jpt_trd_zdf,jp_tem)/tvolt,  tsmo(jpt_trd_zdf,jp_sal)/tvolt
700            WRITE (numout,9445)   tsmo(jpt_trd_npc,jp_tem)/tvolt,  tsmo(jpt_trd_npc,jp_sal)/tvolt
701            WRITE (numout,9446) ( ts2(jpt_trd_xad,jp_tem)+ts2(jpt_trd_yad,jp_tem)+ts2(jpt_trd_zad,jp_tem) )/tvolt,    &
702            &                   ( ts2(jpt_trd_xad,jp_sal)+ts2(jpt_trd_yad,jp_sal)+ts2(jpt_trd_zad,jp_sal) )/tvolt
703            WRITE (numout,9447)   ts2(jpt_trd_ldf,jp_tem)/tvolt,   ts2(jpt_trd_ldf,jp_sal)/tvolt
704            WRITE (numout,9448)   ts2(jpt_trd_zdf,jp_tem)/tvolt,   ts2(jpt_trd_zdf,jp_sal)/tvolt
705            WRITE (numout,9449)   ts2(jpt_trd_npc,jp_tem)/tvolt,   ts2(jpt_trd_npc,jp_sal)/tvolt
706         ENDIF
707
7089440     FORMAT(' tracer consistency at it= ',i6, ' :         temperature','                salinity',   &
709                /, ' ==================================')
7109441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
7119442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
7129443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
7139444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
7149445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
7159446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
7169447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
7179448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
7189449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
719         !
720      ENDIF
721      !
722   END SUBROUTINE trd_twr
723
724#   else
725   !!----------------------------------------------------------------------
726   !!   Default case :                                         Empty module
727   !!----------------------------------------------------------------------
728   INTERFACE trd_icp
729      MODULE PROCEDURE trd_2d, trd_3d
730   END INTERFACE
731
732CONTAINS
733   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )       ! Empty routine
734      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
735      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
736      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
737      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
738      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
739          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype, clpas
740   END SUBROUTINE trd_2d
741   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype, clpas )       ! Empty routine
742      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
743      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
744      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
745      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
746      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
747          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype, clpas
748   END SUBROUTINE trd_3d
749   SUBROUTINE trd_icp_init               ! Empty routine
750   END SUBROUTINE trd_icp_init
751   SUBROUTINE trd_dwr( kt )          ! Empty routine
752      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
753   END SUBROUTINE trd_dwr
754   SUBROUTINE trd_twr( kt )          ! Empty routine
755      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
756   END SUBROUTINE trd_twr
757#endif
758
759   !!======================================================================
760END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.