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.
trdglo.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRD/trdglo.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 28.2 KB
Line 
1MODULE trdglo
2   !!======================================================================
3   !!                       ***  MODULE  trdglo  ***
4   !! Ocean diagnostics:  global domain averaged tracer and momentum trends
5   !!=====================================================================
6   !! History :  1.0  !  2004-08  (C. Talandier) New trends organization
7   !!            3.5  !  2012-02  (G. Madec)  add 3D tracer zdf trend output using iom
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   trd_glo       : domain averaged budget of trends (including kinetic energy and T^2 trends)
12   !!   glo_dyn_wri   : print dynamic trends in ocean.output file
13   !!   glo_tra_wri   : print global T & T^2 trends in ocean.output file
14   !!   trd_glo_init  : initialization step
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean space and time domain variables
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE trd_oce        ! trends: ocean variables
20   USE phycst         ! physical constants
21   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
22   USE ldfdyn         ! ocean dynamics: lateral physics
23   USE zdf_oce        ! ocean vertical physics
24!!gm   USE zdfdrg         ! ocean vertical physics: bottom friction
25   USE zdfddm         ! ocean vertical physics: double diffusion
26   USE eosbn2         ! equation of state
27   USE phycst         ! physical constants
28   !
29   USE lib_mpp        ! distibuted memory computing library
30   USE in_out_manager ! I/O manager
31   USE iom            ! I/O manager library
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   trd_glo       ! called by trdtra and trddyn modules
37   PUBLIC   trd_glo_init  ! called by trdini module
38
39   !                     !!! Variables used for diagnostics
40   REAL(wp) ::   tvolt    ! volume of the whole ocean computed at t-points
41   REAL(wp) ::   tvolu    ! volume of the whole ocean computed at u-points
42   REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
43   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
44   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
45
46   !                     !!! domain averaged trends
47   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
48   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
49   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
50   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
51
52   !! * Substitutions
53#  include "vectopt_loop_substitute.h90"
54#  include "do_loop_substitute.h90"
55   !!----------------------------------------------------------------------
56   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
57   !! $Id$
58   !! Software governed by the CeCILL license (see ./LICENSE)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt, Kmm )
63      !!---------------------------------------------------------------------
64      !!                  ***  ROUTINE trd_glo  ***
65      !!
66      !! ** Purpose :   compute and print global domain averaged trends for
67      !!              T, T^2, momentum, KE, and KE<->PE
68      !!
69      !!----------------------------------------------------------------------
70      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
72      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
73      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
74      INTEGER                   , INTENT(in   ) ::   kt      ! time step
75      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
76      !!
77      INTEGER ::   ji, jj, jk      ! dummy loop indices
78      INTEGER ::   ikbu, ikbv      ! local integers
79      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
80      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
81      !!----------------------------------------------------------------------
82      !
83      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
84         !
85         SELECT CASE( ctype )
86         !
87         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
88            DO_3D_11_11( 1, jpkm1 )
89               zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
90               zvt = ptrdx(ji,jj,jk) * zvm
91               zvs = ptrdy(ji,jj,jk) * zvm
92               tmo(ktrd) = tmo(ktrd) + zvt   
93               smo(ktrd) = smo(ktrd) + zvs
94               t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
95               s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
96            END_3D
97            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
98            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
99               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
100               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
101               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
102               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
103            ENDIF
104            !
105            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
106               !
107               CALL glo_tra_wri( kt )             ! print the results in ocean.output
108               !               
109               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
110               smo(:) = 0._wp
111               t2 (:) = 0._wp
112               s2 (:) = 0._wp
113               !
114            ENDIF
115            !
116         CASE( 'DYN' )          !==  Momentum and KE  ==!       
117            DO_3D_10_10( 1, jpkm1 )
118               zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
119                  &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
120               zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
121                  &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
122               umo(ktrd) = umo(ktrd) + zvt
123               vmo(ktrd) = vmo(ktrd) + zvs
124               hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
125            END_3D
126            !                 
127            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
128               z1_2rau0 = 0.5_wp / rau0
129               DO_2D_10_10
130                  zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
131                     &                                                     * z1_2rau0       * e1e2u(ji,jj)
132                  zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
133                     &                                                     * z1_2rau0       * e1e2v(ji,jj)
134                  umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
135                  vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
136                  hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
137               END_2D
138            ENDIF
139            !                       
140!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
141!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
142!               !
143!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
144!                  z1_2rau0 = 0.5_wp / rau0
145!                  DO jj = 1, jpjm1
146!                     DO ji = 1, jpim1
147!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
148!                        ikbv = mbkv(ji,jj)
149!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
150!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
151!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
152!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
153!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
154!                     END DO
155!                  END DO
156!               ENDIF
157!
158!!gm top drag case is missing
159!
160!               !
161!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
162!               !               
163!               umo(:) = 0._wp                         ! reset for the next time step
164!               vmo(:) = 0._wp
165!               hke(:) = 0._wp
166!               !
167!            ENDIF
168!!gm end
169            !
170         END SELECT
171         !
172      ENDIF
173      !
174   END SUBROUTINE trd_glo
175
176
177   SUBROUTINE glo_dyn_wri( kt, Kmm )
178      !!---------------------------------------------------------------------
179      !!                  ***  ROUTINE glo_dyn_wri  ***
180      !!
181      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
182      !!----------------------------------------------------------------------
183      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
184      INTEGER, INTENT(in) ::   Kmm  ! time level index
185      !
186      INTEGER  ::   ji, jj, jk   ! dummy loop indices
187      REAL(wp) ::   zcof         ! local scalar
188      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
189      !!----------------------------------------------------------------------
190
191      ! I. Momentum trends
192      ! -------------------
193
194      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
195
196         ! I.1 Conversion potential energy - kinetic energy
197         ! --------------------------------------------------
198         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
199         zkx  (:,:,:) = 0._wp
200         zky  (:,:,:) = 0._wp
201         zkz  (:,:,:) = 0._wp
202         zkepe(:,:,:) = 0._wp
203   
204         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
205
206         zcof = 0.5_wp / rau0             ! Density flux at w-point
207         zkz(:,:,1) = 0._wp
208         DO jk = 2, jpk
209            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
210         END DO
211         
212         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
213         DO_3D_10_10( 1, jpkm1 )
214            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
215            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
216         END_3D
217         
218         DO_3D_00_00( 1, jpkm1 )
219            zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
220               &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
221               &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
222               &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
223         END_3D
224
225         ! I.2 Basin averaged kinetic energy trend
226         ! ----------------------------------------
227         peke = 0._wp
228         DO jk = 1, jpkm1
229            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) * e3t(:,:,jk,Kmm) )
230         END DO
231         peke = grav * peke
232
233         ! I.3 Sums over the global domain
234         ! ---------------------------------
235         IF( lk_mpp ) THEN
236            CALL mpp_sum( 'trdglo', peke )
237            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
238            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
239            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
240         ENDIF
241
242         ! I.2 Print dynamic trends in the ocean.output file
243         ! --------------------------------------------------
244
245         IF(lwp) THEN
246            WRITE (numout,*)
247            WRITE (numout,*)
248            WRITE (numout,9500) kt
249            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
250            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
251            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
252            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
253            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
254            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
255            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
256            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
257            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
258            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
259            WRITE (numout,9511)
260            WRITE (numout,9512)                                                 &
261            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
262            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
263            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
264            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
265            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
266            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
267            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
268!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
269         ENDIF
270
271 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
272 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
273 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
274 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
275 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
276 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
277 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
278 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
279 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
280 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
281 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
282 9511    FORMAT(' -----------------------------------------------------------------------------')
283 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
284 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
285 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
286
287         IF(lwp) THEN
288            WRITE (numout,*)
289            WRITE (numout,*)
290            WRITE (numout,9520) kt
291            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
292            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
293            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
294            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
295            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
296            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
297            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
298            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
299            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
300            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
301            WRITE (numout,9531)
302            WRITE (numout,9532)   &
303            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
304            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
305            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
306            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
307!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
308         ENDIF
309
310 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
311 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
312 9522    FORMAT(' ke gradient               u2= ', e20.13)
313 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
314 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
315 9525    FORMAT(' vertical advection        u2= ', e20.13)
316 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
317 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
318 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
319 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
320 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
321 9531    FORMAT(' --------------------------------------------------')
322 9532    FORMAT(' total trend               u2= ', e20.13)
323 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
324 9534    FORMAT('       bottom stress       u2= ', e20.13)
325
326         IF(lwp) THEN
327            WRITE (numout,*)
328            WRITE (numout,*)
329            WRITE (numout,9540) kt
330            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
331            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
332            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
333            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
334            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
335            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
336            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
337            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
338            WRITE (numout,*)
339            WRITE (numout,*)
340         ENDIF
341
342 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
343 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
344 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
345 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
346 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
347 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
348 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
349 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
350 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
351         !
352         ! Save potential to kinetic energy conversion for next time step
353         rpktrd = peke
354         !
355      ENDIF
356      !
357   END SUBROUTINE glo_dyn_wri
358
359
360   SUBROUTINE glo_tra_wri( kt )
361      !!---------------------------------------------------------------------
362      !!                  ***  ROUTINE glo_tra_wri  ***
363      !!
364      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
365      !!----------------------------------------------------------------------
366      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
367      !
368      INTEGER  ::   jk   ! loop indices
369      !!----------------------------------------------------------------------
370
371      ! I. Tracers trends
372      ! -----------------
373
374      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
375
376         ! I.1 Sums over the global domain
377         ! -------------------------------
378         IF( lk_mpp ) THEN
379            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
380            CALL mpp_sum( 'trdglo', smo, jptot_tra )
381            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
382            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
383         ENDIF
384
385         ! I.2 Print tracers trends in the ocean.output file
386         ! --------------------------------------------------
387         
388         IF(lwp) THEN
389            WRITE (numout,*)
390            WRITE (numout,*)
391            WRITE (numout,9400) kt
392            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
393            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
394            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
395            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
396            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
397            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
398            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
399            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
400            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
401            WRITE (numout,9409) 
402            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
403            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
404            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
405            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
406         ENDIF
407
4089400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
409              '              salinity',/' ============================')
4109401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4119411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4129402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4139403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4149404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4159405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4169406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4179407     FORMAT(' penetrative qsr             ',e20.13)
4189408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4199409     FORMAT(' -------------------------------------------------------------------------')
4209410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
421
422
423         IF(lwp) THEN
424            WRITE (numout,*)
425            WRITE (numout,*)
426            WRITE (numout,9420) kt
427            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
428            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
429            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
430            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
431            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
432            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
433            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
434            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
435            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
436            WRITE (numout,9429)
437            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
438            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
439            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
440            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
441         ENDIF
442
4439420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
444            '               salinity', /, ' ===============================')
4459421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4469431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4479422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4489423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4499424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4509425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4519426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4529427     FORMAT(' penetrative qsr           * t   ', e20.13)
4539428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4549429     FORMAT(' -----------------------------------------------------------------------------')
4559430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
456
457
458         IF(lwp) THEN
459            WRITE (numout,*)
460            WRITE (numout,*)
461            WRITE (numout,9440) kt
462            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
463            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
464            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
465            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
466            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
467            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
468            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
469            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
470            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
471            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
472            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
473         ENDIF
474
4759440     FORMAT(' tracer consistency at it= ',i6,   &
476            ' :         temperature','                salinity',/,   &
477            ' ==================================')
4789441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4799442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4809443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4819444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4829445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4839446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4849447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
4859448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
4869449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
487         !
488      ENDIF
489      !
490   END SUBROUTINE glo_tra_wri
491
492
493   SUBROUTINE trd_glo_init( Kmm )
494      !!---------------------------------------------------------------------
495      !!                  ***  ROUTINE trd_glo_init  ***
496      !!
497      !! ** Purpose :   Read the namtrd namelist
498      !!----------------------------------------------------------------------
499      INTEGER, INTENT(in) ::   Kmm   ! time level index
500      INTEGER  ::   ji, jj, jk   ! dummy loop indices
501      !!----------------------------------------------------------------------
502
503      IF(lwp) THEN
504         WRITE(numout,*)
505         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
506         WRITE(numout,*) '~~~~~~~~~~~~~'
507      ENDIF
508
509      ! Total volume at t-points:
510      tvolt = 0._wp
511      DO jk = 1, jpkm1
512         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
513      END DO
514      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
515
516      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
517
518      ! Initialization of potential to kinetic energy conversion
519      rpktrd = 0._wp
520
521      ! Total volume at u-, v- points:
522!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
523      tvolu = 0._wp
524      tvolv = 0._wp
525
526      DO_3D_00_00( 1, jpk )
527         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
528         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
529      END_3D
530      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
531      CALL mpp_sum( 'trdglo', tvolv )
532
533      IF(lwp) THEN
534         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
535         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
536      ENDIF
537      !
538   END SUBROUTINE trd_glo_init
539
540   !!======================================================================
541END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.