New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdicp.F90 in branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 4606

Last change on this file since 4606 was 4606, checked in by gm, 10 years ago

#1238 - v3.4_STABLE: correct the volume calculation for tendency diagnostics

  • Property svn:keywords set to Id
File size: 29.9 KB
RevLine 
[215]1MODULE trdicp
2   !!======================================================================
3   !!                       ***  MODULE  trdicp  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic trends
5   !!=====================================================================
[2528]6   !! History :  1.0  !  2004-08 (C. Talandier) New trends organization
[503]7   !!----------------------------------------------------------------------
[215]8#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
9   !!----------------------------------------------------------------------
10   !!   'key_trdtra'  or                  active tracers trends diagnostics
11   !!   'key_trddyn'                            momentum trends diagnostics
12   !!----------------------------------------------------------------------
[503]13   !!   trd_icp          : compute the basin averaged properties for tra/dyn
[215]14   !!   trd_dwr          : print dynmaic trends in ocean.output file
15   !!   trd_twr          : print tracers trends in ocean.output file
16   !!   trd_icp_init     : initialization step
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE trdmod_oce      ! ocean variables trends
21   USE ldftra_oce      ! ocean active tracers: lateral physics
22   USE ldfdyn_oce      ! ocean dynamics: lateral physics
23   USE zdf_oce         ! ocean vertical physics
24   USE in_out_manager  ! I/O manager
25   USE lib_mpp         ! distibuted memory computing library
26   USE eosbn2          ! equation of state
27   USE phycst          ! physical constants
[3294]28   USE wrk_nemo        ! Memory allocation
[215]29
[3294]30
[215]31   IMPLICIT NONE
32   PRIVATE
33
[503]34   INTERFACE trd_icp
[215]35      MODULE PROCEDURE trd_2d, trd_3d
36   END INTERFACE
37
[503]38   PUBLIC   trd_icp       ! called by trdmod.F90
39   PUBLIC   trd_dwr       ! called by step.F90
40   PUBLIC   trd_twr       ! called by step.F90
41   PUBLIC   trd_icp_init  ! called by opa.F90
[215]42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
[2528]47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
[2715]49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[215]50   !!----------------------------------------------------------------------
51CONTAINS
52
[2528]53   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )
[215]54      !!---------------------------------------------------------------------
55      !!                  ***  ROUTINE trd_2d  ***
56      !!
57      !! ** Purpose : verify the basin averaged properties of tracers and/or
[1601]58      !!              momentum equations at every time step frequency nn_trd.
[503]59      !!----------------------------------------------------------------------
[2715]60      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx   ! Temperature or U trend
61      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy   ! Salinity    or V trend
62      INTEGER                     , INTENT(in   ) ::   ktrd      ! tracer trend index
63      CHARACTER(len=3)            , INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
[215]64      !!
[2715]65      INTEGER ::   ji, jj   ! loop indices
[215]66      !!----------------------------------------------------------------------
67
[2715]68      SELECT CASE( ctype )    !==  Mask trends  ==!
[503]69      !
[2715]70      CASE( 'DYN' )                    ! Momentum
[215]71         DO jj = 1, jpjm1
72            DO ji = 1, jpim1
[2715]73               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
74               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
[215]75            END DO
76         END DO
[2715]77         ptrd2dx(jpi, : ) = 0._wp      ;      ptrd2dy(jpi, : ) = 0._wp
78         ptrd2dx( : ,jpj) = 0._wp      ;      ptrd2dy( : ,jpj) = 0._wp
[503]79         !
[2715]80      CASE( 'TRA' )                    ! Tracers
[215]81         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
82         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
[503]83         !
[215]84      END SELECT
85     
[2715]86      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
[503]87      !
[2715]88      CASE( 'DYN' )                    ! Momentum
89         umo(ktrd) = 0._wp
90         vmo(ktrd) = 0._wp
[503]91         !
92         SELECT CASE( ktrd )
93         CASE( jpdyn_trd_swf )         ! surface forcing
[2715]94            umo(ktrd) = SUM( ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1) )
95            vmo(ktrd) = SUM( ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1) )
[215]96         END SELECT
[503]97         !
98      CASE( 'TRA' )              ! Tracers
[2715]99         tmo(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
100         smo(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) )
[215]101      END SELECT
102     
[2715]103      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
[503]104      !
105      CASE( 'DYN' )              ! Momentum
[2715]106         hke(ktrd) = SUM(   un(:,:,1) * ptrd2dx(:,:) * e1u(:,:) * e2u(:,:) * fse3u(:,:,1)   &
107            &             + vn(:,:,1) * ptrd2dy(:,:) * e1v(:,:) * e2v(:,:) * fse3v(:,:,1)   )
[503]108         !
109      CASE( 'TRA' )              ! Tracers
[3294]110         t2(ktrd) = SUM( ptrd2dx(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_tem) )
111         s2(ktrd) = SUM( ptrd2dy(:,:) * e1e2t(:,:) * fse3t(:,:,1) * tsn(:,:,1,jp_sal) )
[503]112         !     
[215]113      END SELECT
[503]114      !
[215]115   END SUBROUTINE trd_2d
116
117
[2528]118   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype )
[215]119      !!---------------------------------------------------------------------
120      !!                  ***  ROUTINE trd_3d  ***
121      !!
122      !! ** Purpose : verify the basin averaged properties of tracers and/or
[1601]123      !!              momentum equations at every time step frequency nn_trd.
[503]124      !!----------------------------------------------------------------------
[2715]125      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx   ! Temperature or U trend
126      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy   ! Salinity    or V trend
127      INTEGER,                          INTENT(in   ) ::   ktrd      ! momentum or tracer trend index
128      CHARACTER(len=3),                 INTENT(in   ) ::   ctype     ! momentum ('DYN') or tracers ('TRA') trends
[215]129      !!
[2715]130      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[215]131      !!----------------------------------------------------------------------
132
[2715]133      SELECT CASE( ctype )    !==  Mask the trends  ==!
[503]134      !
135      CASE( 'DYN' )              ! Momentum       
[2715]136         DO jk = 1, jpkm1
[215]137            DO jj = 1, jpjm1
138               DO ji = 1, jpim1
[2715]139                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
140                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
[503]141               END DO
142            END DO
143         END DO
[2715]144         ptrd3dx(jpi, : ,:) = 0._wp      ;      ptrd3dy(jpi, : ,:) = 0._wp
145         ptrd3dx( : ,jpj,:) = 0._wp      ;      ptrd3dy( : ,jpj,:) = 0._wp
[503]146         !
147      CASE( 'TRA' )              ! Tracers
[2715]148         DO jk = 1, jpkm1
[215]149            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
150            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
[503]151         END DO
152         !
[215]153      END SELECT   
154
[2715]155      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum trends  ==!
[503]156      !
157      CASE( 'DYN' )              ! Momentum
[2715]158         umo(ktrd) = 0._wp
159         vmo(ktrd) = 0._wp
160         DO jk = 1, jpkm1
161            umo(ktrd) = umo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk) )
162            vmo(ktrd) = vmo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk) )
[215]163         END DO
[503]164         !
165      CASE( 'TRA' )              ! Tracers
[2715]166         tmo(ktrd) = 0._wp
167         smo(ktrd) = 0._wp
[215]168         DO jk = 1, jpkm1
[2715]169            tmo(ktrd) = tmo(ktrd) + SUM( ptrd3dx(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
170            smo(ktrd) = smo(ktrd) + SUM( ptrd3dy(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]171         END DO
[503]172         !
[215]173      END SELECT
174
[2715]175      SELECT CASE( ctype )    !==  Basin averaged tracer/momentum square trends  ==!   (now field)
[503]176      !
177      CASE( 'DYN' )              ! Momentum
[2715]178         hke(ktrd) = 0._wp
179         DO jk = 1, jpkm1
180            hke(ktrd) = hke(ktrd) + SUM(   un(:,:,jk) * ptrd3dx(:,:,jk) * e1u(:,:) * e2u(:,:) * fse3u(:,:,jk)   &
181               &                         + vn(:,:,jk) * ptrd3dy(:,:,jk) * e1v(:,:) * e2v(:,:) * fse3v(:,:,jk)   )
[215]182         END DO
[503]183         !
184      CASE( 'TRA' )              ! Tracers
[2715]185         t2(ktrd) = 0._wp
186         s2(ktrd) = 0._wp
187         DO jk = 1, jpkm1
[3294]188            t2(ktrd) = t2(ktrd) + SUM( ptrd3dx(:,:,jk) * tsn(:,:,jk,jp_tem) * e1e2t(:,:) * fse3t(:,:,jk) )
189            s2(ktrd) = s2(ktrd) + SUM( ptrd3dy(:,:,jk) * tsn(:,:,jk,jp_sal) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]190         END DO
[503]191         !
[215]192      END SELECT
[503]193      !
[215]194   END SUBROUTINE trd_3d
195
196
197   SUBROUTINE trd_icp_init
198      !!---------------------------------------------------------------------
199      !!                  ***  ROUTINE trd_icp_init  ***
200      !!
[503]201      !! ** Purpose :   Read the namtrd namelist
[215]202      !!----------------------------------------------------------------------
[2715]203      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[215]204      !!----------------------------------------------------------------------
205
206      IF(lwp) THEN
207         WRITE(numout,*)
208         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
209         WRITE(numout,*) '~~~~~~~~~~~~~'
210      ENDIF
211
212      ! Total volume at t-points:
[2715]213      tvolt = 0._wp
[215]214      DO jk = 1, jpkm1
[4606]215         tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
[215]216      END DO
217      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
218
[503]219      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
[215]220
221#if  defined key_trddyn
222      ! Initialization of potential to kinetic energy conversion
[2715]223      rpktrd = 0._wp
[215]224
225      ! Total volume at u-, v- points:
[2715]226      tvolu = 0._wp
227      tvolv = 0._wp
[215]228      DO jk = 1, jpk
229         DO jj = 2, jpjm1
230            DO ji = fs_2, fs_jpim1   ! vector opt.
[2715]231               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
232               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
[215]233            END DO
234         END DO
235      END DO
236      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
237      IF( lk_mpp )   CALL mpp_sum( tvolv )
238
239      IF(lwp) THEN
[503]240         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
241         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
[215]242      ENDIF
243#endif
[503]244      !
[215]245   END SUBROUTINE trd_icp_init
246
247
248   SUBROUTINE trd_dwr( kt )
249      !!---------------------------------------------------------------------
250      !!                  ***  ROUTINE trd_dwr  ***
251      !!
252      !! ** Purpose :  write dynamic trends in ocean.output
[503]253      !!----------------------------------------------------------------------
[2715]254      !
255      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
256      !
257      INTEGER  ::   ji, jj, jk   ! dummy loop indices
258      REAL(wp) ::   zcof         ! local scalar
[3294]259      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
[215]260      !!----------------------------------------------------------------------
261
[3294]262      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
[2715]263
[215]264      ! I. Momentum trends
265      ! -------------------
266
[1601]267      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
[215]268
269         ! I.1 Conversion potential energy - kinetic energy
270         ! --------------------------------------------------
271         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
[2715]272         zkx  (:,:,:) = 0._wp
273         zky  (:,:,:) = 0._wp
274         zkz  (:,:,:) = 0._wp
275         zkepe(:,:,:) = 0._wp
[215]276   
[2528]277         CALL eos( tsn, rhd, rhop )       ! now potential and in situ densities
[215]278
[2715]279         zcof = 0.5_wp / rau0             ! Density flux at w-point
280         zkz(:,:,1) = 0._wp
[215]281         DO jk = 2, jpk
[2715]282            zkz(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
[215]283         END DO
284         
[2715]285         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
286         DO jk = 1, jpkm1
[215]287            DO jj = 1, jpjm1
288               DO ji = 1, jpim1
[2715]289                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
290                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
[215]291               END DO
292            END DO
293         END DO
294         
[2715]295         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
[215]296            DO jj = 2, jpjm1
297               DO ji = 2, jpim1
[2715]298                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
299                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
300                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
301                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
[215]302               END DO
303            END DO
304         END DO
305
306         ! I.2 Basin averaged kinetic energy trend
307         ! ----------------------------------------
[2715]308         peke = 0._wp
309         DO jk = 1, jpkm1
310            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
[215]311         END DO
[2715]312         peke = grav * peke
[215]313
314         ! I.3 Sums over the global domain
315         ! ---------------------------------
316         IF( lk_mpp ) THEN
[503]317            CALL mpp_sum( peke )
318            CALL mpp_sum( umo , jptot_dyn )
319            CALL mpp_sum( vmo , jptot_dyn )
320            CALL mpp_sum( hke , jptot_dyn )
321         ENDIF
[215]322
323         ! I.2 Print dynamic trends in the ocean.output file
324         ! --------------------------------------------------
325
326         IF(lwp) THEN
327            WRITE (numout,*)
328            WRITE (numout,*)
329            WRITE (numout,9500) kt
[503]330            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
331            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
332            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
333            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
334            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
[1129]335            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
336            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
337            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
338            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
339            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
340            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
341            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
342            WRITE (numout,9513)
343            WRITE (numout,9514)                                                 &
[503]344            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
[1129]345            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
346            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
[503]347            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
[1129]348            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
349            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
[215]350         ENDIF
351
352 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
353 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
354 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
355 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
356 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
357 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
[1129]358 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
359 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
360 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
361 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
362 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
363 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
364 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
365 9513    FORMAT(' -----------------------------------------------------------------------------')
366 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
[215]367
368         IF(lwp) THEN
369            WRITE (numout,*)
370            WRITE (numout,*)
371            WRITE (numout,9520) kt
[503]372            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
373            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
374            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
375            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
376            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
[1129]377            WRITE (numout,9526) hke(jpicpd_had) / tvolt
378            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
379            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
380            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
381            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
382            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
[1708]383            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
384            WRITE (numout,9533)
385            WRITE (numout,9534)   &
[503]386            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
[1129]387            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
[1708]388            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
[215]389         ENDIF
390
391 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
392 9521    FORMAT(' pressure gradient         u2= ', e20.13)
393 9522    FORMAT(' ke gradient               u2= ', e20.13)
394 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
395 9524    FORMAT(' coriolis term             u2= ', e20.13)
396 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
[1129]397 9526    FORMAT(' horizontal advection      u2= ', e20.13)
398 9527    FORMAT(' vertical advection        u2= ', e20.13)
399 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
400 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
401 9530    FORMAT(' surface wind forcing      u2= ', e20.13)
402 9531    FORMAT(' dampimg term              u2= ', e20.13)
[1708]403 9532    FORMAT(' bottom flux               u2= ', e20.13)
404 9533    FORMAT(' --------------------------------------------------')
405 9534    FORMAT(' total trend               u2= ', e20.13)
[215]406
407         IF(lwp) THEN
408            WRITE (numout,*)
409            WRITE (numout,*)
410            WRITE (numout,9540) kt
[1129]411            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
412            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
[503]413            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
414            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
415            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
416            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
417            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
418            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
419            WRITE (numout,*)
420            WRITE (numout,*)
[215]421         ENDIF
422
423 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
424 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
[1129]425 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
[215]426 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
[503]427 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
428 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
429 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
430 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
431 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
432         !
[215]433         ! Save potential to kinetic energy conversion for next time step
434         rpktrd = peke
[503]435         !
[215]436      ENDIF
[503]437      !
[3294]438      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
[2715]439      !
[215]440   END SUBROUTINE trd_dwr
441
442
443   SUBROUTINE trd_twr( kt )
444      !!---------------------------------------------------------------------
445      !!                  ***  ROUTINE trd_twr  ***
446      !!
447      !! ** Purpose :  write active tracers trends in ocean.output
448      !!----------------------------------------------------------------------
[503]449      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[215]450      !!----------------------------------------------------------------------
451
452      ! I. Tracers trends
453      ! -----------------
454
[1601]455      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
[215]456
457         ! I.1 Sums over the global domain
458         ! -------------------------------
459         IF( lk_mpp ) THEN
[503]460            CALL mpp_sum( tmo, jptot_tra )   
461            CALL mpp_sum( smo, jptot_tra )
462            CALL mpp_sum( t2 , jptot_tra )
463            CALL mpp_sum( s2 , jptot_tra )
[215]464         ENDIF
465
466         ! I.2 Print tracers trends in the ocean.output file
467         ! --------------------------------------------------
468         
469         IF(lwp) THEN
470            WRITE (numout,*)
471            WRITE (numout,*)
472            WRITE (numout,9400) kt
[2528]473            WRITE (numout,9401)  tmo(jpicpt_xad) / tvolt, smo(jpicpt_xad) / tvolt
474            WRITE (numout,9411)  tmo(jpicpt_yad) / tvolt, smo(jpicpt_yad) / tvolt
[503]475            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
476            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
477            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
478            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
479            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
480            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
481            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
482            WRITE (numout,9409) 
483            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
484            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
485            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
486            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
[215]487         ENDIF
488
4899400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
490              '              salinity',/' ============================')
[2528]4919401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4929411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
[215]4939402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4949403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4959404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
[503]4969405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
[215]4979406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
[503]4989407     FORMAT(' penetrative qsr             ',e20.13)
4999408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
[215]5009409     FORMAT(' -------------------------------------------------------------------------')
5019410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
502
503
504         IF(lwp) THEN
505            WRITE (numout,*)
506            WRITE (numout,*)
507            WRITE (numout,9420) kt
[2528]508            WRITE (numout,9421)   t2(jpicpt_xad) / tvolt, s2(jpicpt_xad) / tvolt
509            WRITE (numout,9431)   t2(jpicpt_yad) / tvolt, s2(jpicpt_yad) / tvolt
[503]510            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
511            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
512            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
513            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
514            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
515            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
516            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
[215]517            WRITE (numout,9429)
[503]518            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
519            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
520            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
521            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
[215]522         ENDIF
523
5249420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
525            '               salinity', /, ' ===============================')
[2528]5269421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
5279431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
[215]5289422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
5299423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
5309424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
[503]5319425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
[215]5329426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
[503]5339427     FORMAT(' penetrative qsr           * t   ', e20.13)
5349428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
[215]5359429     FORMAT(' -----------------------------------------------------------------------------')
5369430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
537
538
539         IF(lwp) THEN
540            WRITE (numout,*)
541            WRITE (numout,*)
542            WRITE (numout,9440) kt
[503]543            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
544            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
545            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
546            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
547            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
548            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
549            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
550            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
551            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
552            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
553            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
[215]554         ENDIF
555
5569440     FORMAT(' tracer consistency at it= ',i6,   &
557            ' :         temperature','                salinity',/,   &
558            ' ==================================')
[503]5599441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5609442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5619443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5629444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5639445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5649446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5659447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5669448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5679449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
568         !
[215]569      ENDIF
[503]570      !
[215]571   END SUBROUTINE trd_twr
572
573#   else
574   !!----------------------------------------------------------------------
575   !!   Default case :                                         Empty module
576   !!----------------------------------------------------------------------
[601]577   INTERFACE trd_icp
578      MODULE PROCEDURE trd_2d, trd_3d
579   END INTERFACE
580
[215]581CONTAINS
[2528]582   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype )       ! Empty routine
[503]583      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
[601]584      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
585      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
[521]586      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
[2528]587          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype
[215]588   END SUBROUTINE trd_2d
[2528]589   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype )       ! Empty routine
[503]590      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
[601]591      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
592      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
[521]593      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
[2528]594          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype
[215]595   END SUBROUTINE trd_3d
596   SUBROUTINE trd_icp_init               ! Empty routine
597   END SUBROUTINE trd_icp_init
598   SUBROUTINE trd_dwr( kt )          ! Empty routine
599      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
600   END SUBROUTINE trd_dwr
601   SUBROUTINE trd_twr( kt )          ! Empty routine
602      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
603   END SUBROUTINE trd_twr
604#endif
605
606   !!======================================================================
607END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.