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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdicp.F90 @ 2733

Last change on this file since 2733 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

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