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 @ 1388

Last change on this file since 1388 was 1152, checked in by rblod, 16 years ago

Convert cvs header to svn Id, step II

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.8 KB
RevLine 
[215]1MODULE trdicp
2   !!======================================================================
3   !!                       ***  MODULE  trdicp  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic trends
5   !!=====================================================================
[503]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   !!----------------------------------------------------------------------
[215]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   !!----------------------------------------------------------------------
[503]18   !!   trd_icp          : compute the basin averaged properties for tra/dyn
[215]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
[503]37   INTERFACE trd_icp
[215]38      MODULE PROCEDURE trd_2d, trd_3d
39   END INTERFACE
40
[503]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
[215]45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
[247]50   !!   OPA 9.0 , LOCEAN-IPSL (2005)
[1152]51   !! $Id$
[503]52   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[215]53   !!----------------------------------------------------------------------
54
55CONTAINS
56
[503]57   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )
[215]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.
[503]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
[215]69      !!
[503]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                                            !    "         "
[215]74      !!----------------------------------------------------------------------
75
[503]76      ! Control of optional arguments
77      cpas = 'fst'
78      IF( PRESENT(clpas) )  cpas = clpas
[215]79
[503]80      ! 1. Mask trends
81      ! --------------
[215]82
[503]83      SELECT CASE( ctype )
84      !
85      CASE( 'DYN' )              ! Momentum
[215]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
[503]96         !
97      CASE( 'TRA' )              ! Tracers
[215]98         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
99         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
[503]100         !
[215]101      END SELECT
102     
[503]103      ! 2. Basin averaged tracer/momentum trends
104      ! ----------------------------------------
[215]105
[503]106      SELECT CASE( ctype )
107      !
108      CASE( 'DYN' )              ! Momentum
[215]109         umo(ktrd) = 0.e0
110         vmo(ktrd) = 0.e0
[503]111         !
112         SELECT CASE( ktrd )
113         !
114         CASE( jpdyn_trd_swf )         ! surface forcing
[215]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
[503]121            !
122         CASE( jpdyn_trd_bfr )         ! bottom friction fluxes
[215]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
[503]129            !
[215]130         END SELECT
[503]131         !
132      CASE( 'TRA' )              ! Tracers
133         IF( cpas == 'fst' )   THEN
134            tmo(ktrd) = 0.e0
135            smo(ktrd) = 0.e0
136         ENDIF
[215]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
[503]144         !
[215]145      END SELECT
146     
[503]147      ! 3. Basin averaged tracer/momentum square trends
148      ! ----------------------------------------------
[215]149      ! c a u t i o n: field now
150     
[503]151      SELECT CASE( ctype )
152      !
153      CASE( 'DYN' )              ! Momentum
[215]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
[503]164         !
165      CASE( 'TRA' )              ! Tracers
166         IF( cpas == 'fst' )   THEN
167            t2(ktrd) = 0.e0
168            s2(ktrd) = 0.e0
169         ENDIF
[215]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
[503]177         !     
[215]178      END SELECT
[503]179      !
[215]180   END SUBROUTINE trd_2d
181
182
[503]183   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype, clpas )
[215]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.
[503]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
[215]195      !!
196      INTEGER ::   ji, jj, jk
[503]197      CHARACTER(len=3) ::   cpas                                              ! number of passage
198      REAL(wp) ::   zbt, zbtu, zbtv, zmsku, zmskv                             ! temporary scalars
[215]199      !!----------------------------------------------------------------------
200
[503]201      ! Control of optional arguments
202      cpas = 'fst'
203      IF( PRESENT(clpas) )  cpas = clpas
[215]204
[503]205      ! 1. Mask the trends
206      ! ------------------
[215]207
[503]208      SELECT CASE( ctype )
209      !
210      CASE( 'DYN' )              ! Momentum       
[215]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
[503]218               END DO
219            END DO
220         END DO
[215]221         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
222         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
[503]223         !
224      CASE( 'TRA' )              ! Tracers
[215]225         DO jk = 1, jpk
226            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
227            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
[503]228         END DO
229         !
[215]230      END SELECT   
231
[503]232      ! 2. Basin averaged tracer/momentum trends
233      ! ----------------------------------------
[215]234     
[503]235      SELECT CASE( ctype )
236      !
237      CASE( 'DYN' )              ! Momentum
[215]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
[503]250         !
251      CASE( 'TRA' )              ! Tracers
252         IF( cpas == 'fst' )   THEN
253            tmo(ktrd) = 0.e0
254            smo(ktrd) = 0.e0
255         ENDIF
[215]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
[503]265         !
[215]266      END SELECT
267
[503]268      ! 3. Basin averaged tracer/momentum square trends
269      ! -----------------------------------------------
[215]270      ! c a u t i o n: field now
271     
[503]272      SELECT CASE( ctype )
273      !
274      CASE( 'DYN' )              ! Momentum
[215]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
[503]287         !
288      CASE( 'TRA' )              ! Tracers
289         IF( cpas == 'fst' )   THEN
290            t2(ktrd) = 0.e0
291            s2(ktrd) = 0.e0
292         ENDIF
[215]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
[503]302         !
[215]303      END SELECT
[503]304      !
[215]305   END SUBROUTINE trd_3d
306
307
308
309   SUBROUTINE trd_icp_init
310      !!---------------------------------------------------------------------
311      !!                  ***  ROUTINE trd_icp_init  ***
312      !!
[503]313      !! ** Purpose :   Read the namtrd namelist
[215]314      !!----------------------------------------------------------------------
[503]315      INTEGER  ::   ji, jj, jk
[215]316      REAL(wp) ::   zmskt
317#if  defined key_trddyn
[503]318      REAL(wp) ::   zmsku, zmskv
[215]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
[503]340      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
[215]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
[503]364         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
365         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
[215]366      ENDIF
367#endif
[503]368      !
[215]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
[503]377      !!----------------------------------------------------------------------
378      INTEGER, INTENT(in) ::   kt                                  ! ocean time-step index
[215]379      !!
[503]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
[215]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
[503]394         zkx(:,:,:)   = 0.e0
395         zky(:,:,:)   = 0.e0
396         zkz(:,:,:)   = 0.e0
[215]397         zkepe(:,:,:) = 0.e0
398   
399         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
400
[503]401         ! Density flux at w-point
[215]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
[503]410         zkz(:,:,1) = 0.e0
[215]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
[503]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
[215]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
[503]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
[1129]474            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
475            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
476            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
477            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
478            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
479            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
480            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
481            WRITE (numout,9513)
482            WRITE (numout,9514)                                                 &
[503]483            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
[1129]484            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
485            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
[503]486            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
[1129]487            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
488            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
[215]489         ENDIF
490
491 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
492 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
493 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
494 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
495 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
496 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
[1129]497 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
498 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
499 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
500 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
501 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
502 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
503 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
504 9513    FORMAT(' -----------------------------------------------------------------------------')
505 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
[215]506
507         IF(lwp) THEN
508            WRITE (numout,*)
509            WRITE (numout,*)
510            WRITE (numout,9520) kt
[503]511            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
512            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
513            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
514            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
515            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
[1129]516            WRITE (numout,9526) hke(jpicpd_had) / tvolt
517            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
518            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
519            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
520            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
521            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
522            WRITE (numout,9532)
523            WRITE (numout,9533)   &
[503]524            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
[1129]525            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
526            &      + hke(jpicpd_swf) ) / tvolt
[215]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)
[1129]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(' --------------------------------------------------')
542 9533    FORMAT(' total trend               u2= ', e20.13)
[215]543
544         IF(lwp) THEN
545            WRITE (numout,*)
546            WRITE (numout,*)
547            WRITE (numout,9540) kt
[1129]548            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
549            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
[503]550            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
551            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
552            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
553            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
554            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
555            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
556            WRITE (numout,*)
557            WRITE (numout,*)
[215]558         ENDIF
559
560 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
561 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
[1129]562 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
[215]563 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
[503]564 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
565 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
566 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
567 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
568 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
569         !
[215]570         ! Save potential to kinetic energy conversion for next time step
571         rpktrd = peke
[503]572         !
[215]573      ENDIF
[503]574      !
[215]575   END SUBROUTINE trd_dwr
576
577
578   SUBROUTINE trd_twr( kt )
579      !!---------------------------------------------------------------------
580      !!                  ***  ROUTINE trd_twr  ***
581      !!
582      !! ** Purpose :  write active tracers trends in ocean.output
583      !!----------------------------------------------------------------------
[503]584      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[215]585      !!----------------------------------------------------------------------
586
587      ! I. Tracers trends
588      ! -----------------
589
590      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
591
592         ! I.1 Sums over the global domain
593         ! -------------------------------
594         IF( lk_mpp ) THEN
[503]595            CALL mpp_sum( tmo, jptot_tra )   
596            CALL mpp_sum( smo, jptot_tra )
597            CALL mpp_sum( t2 , jptot_tra )
598            CALL mpp_sum( s2 , jptot_tra )
[215]599         ENDIF
600
601         ! I.2 Print tracers trends in the ocean.output file
602         ! --------------------------------------------------
603         
604         IF(lwp) THEN
605            WRITE (numout,*)
606            WRITE (numout,*)
607            WRITE (numout,9400) kt
[503]608            WRITE (numout,9401) (tmo(jpicpt_xad)+tmo(jpicpt_yad))/ tvolt, (smo(jpicpt_xad)+smo(jpicpt_yad))/ tvolt
609            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
610            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
611            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
612            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
613            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
614            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
615            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
616            WRITE (numout,9409) 
617            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
618            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
619            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
620            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
[215]621         ENDIF
622
6239400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
624              '              salinity',/' ============================')
6259401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
6269402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
6279403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
6289404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
[503]6299405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
[215]6309406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
[503]6319407     FORMAT(' penetrative qsr             ',e20.13)
6329408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
[215]6339409     FORMAT(' -------------------------------------------------------------------------')
6349410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
635
636
637         IF(lwp) THEN
638            WRITE (numout,*)
639            WRITE (numout,*)
640            WRITE (numout,9420) kt
[503]641            WRITE (numout,9421) ( t2(jpicpt_xad)+t2(jpicpt_yad) )/ tvolt, ( s2(jpicpt_xad)+s2(jpicpt_yad) )/ tvolt
642            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
643            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
644            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
645            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
646            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
647            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
648            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
[215]649            WRITE (numout,9429)
[503]650            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
651            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
652            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
653            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
[215]654         ENDIF
655
6569420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
657            '               salinity', /, ' ===============================')
6589421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
6599422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
6609423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
6619424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
[503]6629425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
[215]6639426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
[503]6649427     FORMAT(' penetrative qsr           * t   ', e20.13)
6659428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
[215]6669429     FORMAT(' -----------------------------------------------------------------------------')
6679430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
668
669
670         IF(lwp) THEN
671            WRITE (numout,*)
672            WRITE (numout,*)
673            WRITE (numout,9440) kt
[503]674            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
675            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
676            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
677            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
678            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
679            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
680            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
681            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
682            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
683            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
684            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
[215]685         ENDIF
686
6879440     FORMAT(' tracer consistency at it= ',i6,   &
688            ' :         temperature','                salinity',/,   &
689            ' ==================================')
[503]6909441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
6919442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
6929443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
6939444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
6949445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
6959446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
6969447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
6979448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
6989449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
699         !
[215]700      ENDIF
[503]701      !
[215]702   END SUBROUTINE trd_twr
703
704#   else
705   !!----------------------------------------------------------------------
706   !!   Default case :                                         Empty module
707   !!----------------------------------------------------------------------
[601]708   INTERFACE trd_icp
709      MODULE PROCEDURE trd_2d, trd_3d
710   END INTERFACE
711
[215]712CONTAINS
[521]713   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )       ! Empty routine
[503]714      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
[601]715      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
716      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
717      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
[521]718      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
719          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype, clpas
[215]720   END SUBROUTINE trd_2d
[521]721   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype, clpas )       ! Empty routine
[503]722      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
[601]723      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
724      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
725      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
[521]726      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
727          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype, clpas
[215]728   END SUBROUTINE trd_3d
729   SUBROUTINE trd_icp_init               ! Empty routine
730   END SUBROUTINE trd_icp_init
731   SUBROUTINE trd_dwr( kt )          ! Empty routine
732      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
733   END SUBROUTINE trd_dwr
734   SUBROUTINE trd_twr( kt )          ! Empty routine
735      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
736   END SUBROUTINE trd_twr
737#endif
738
739   !!======================================================================
740END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.