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/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRD – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRD/trdglo.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

Mixed precision version, tested up to 30 years on ORCA2.
  • Property svn:keywords set to Id
File size: 28.6 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(dp) ::   tvolt    ! volume of the whole ocean computed at t-points
41   REAL(dp) ::   tvolu    ! volume of the whole ocean computed at u-points
42   REAL(dp) ::   tvolv    ! volume of the whole ocean computed at v-points
43   REAL(dp) ::   rpktrd   ! potential to kinetic energy conversion
44   REAL(dp) ::   peke     ! conversion potential energy - kinetic energy trend
45
46   !                     !!! domain averaged trends
47   REAL(dp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
48   REAL(dp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
49   REAL(dp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
50   REAL(dp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
51
52   !! * Substitutions
53#  include "do_loop_substitute.h90"
54#  include "domzgr_substitute.h90"
55#  include "single_precision_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt, Kmm )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE trd_glo  ***
66      !!
67      !! ** Purpose :   compute and print global domain averaged trends for
68      !!              T, T^2, momentum, KE, and KE<->PE
69      !!
70      !!----------------------------------------------------------------------
71      REAL(dp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
72      REAL(dp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
73      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
74      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
75      INTEGER                   , INTENT(in   ) ::   kt      ! time step
76      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
77      !!
78      INTEGER ::   ji, jj, jk      ! dummy loop indices
79      INTEGER ::   ikbu, ikbv      ! local integers
80      REAL(dp)::   zvm, zvt, zvs, z1_2rho0   ! local scalars
81      REAL(dp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
82      !!----------------------------------------------------------------------
83      !
84      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
85         !
86         SELECT CASE( ctype )
87         !
88         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
89            DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! global sum of mask volume trend and trend*T (including interior mask)
90               zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
91               zvt = ptrdx(ji,jj,jk) * zvm
92               zvs = ptrdy(ji,jj,jk) * zvm
93               tmo(ktrd) = tmo(ktrd) + zvt   
94               smo(ktrd) = smo(ktrd) + zvs
95               t2 (ktrd) = t2(ktrd)  + zvt * ts(ji,jj,jk,jp_tem,Kmm)
96               s2 (ktrd) = s2(ktrd)  + zvs * ts(ji,jj,jk,jp_sal,Kmm)
97            END_3D
98            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
99            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
100               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
101               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
102               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
103               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
104            ENDIF
105            !
106            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
107               !
108               CALL glo_tra_wri( kt )             ! print the results in ocean.output
109               !               
110               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
111               smo(:) = 0._wp
112               t2 (:) = 0._wp
113               s2 (:) = 0._wp
114               !
115            ENDIF
116            !
117         CASE( 'DYN' )          !==  Momentum and KE  ==!       
118            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
119               zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
120                  &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
121               zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
122                  &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
123               umo(ktrd) = umo(ktrd) + zvt
124               vmo(ktrd) = vmo(ktrd) + zvs
125               hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
126            END_3D
127            !                 
128            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
129               z1_2rho0 = 0.5_wp / rho0
130               DO_2D( 1, 0, 1, 0 )
131                  zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
132                     &                                                     * z1_2rho0       * e1e2u(ji,jj)
133                  zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
134                     &                                                     * z1_2rho0       * e1e2v(ji,jj)
135                  umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
136                  vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
137                  hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
138               END_2D
139            ENDIF
140            !                       
141!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
142!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
143!               !
144!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
145!                  z1_2rho0 = 0.5_wp / rho0
146!                  DO jj = 1, jpjm1
147!                     DO ji = 1, jpim1
148!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
149!                        ikbv = mbkv(ji,jj)
150!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
151!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
152!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
153!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
154!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
155!                     END DO
156!                  END DO
157!               ENDIF
158!
159!!gm top drag case is missing
160!
161!               !
162!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
163!               !               
164!               umo(:) = 0._wp                         ! reset for the next time step
165!               vmo(:) = 0._wp
166!               hke(:) = 0._wp
167!               !
168!            ENDIF
169!!gm end
170            !
171         END SELECT
172         !
173      ENDIF
174      !
175   END SUBROUTINE trd_glo
176
177
178   SUBROUTINE glo_dyn_wri( kt, Kmm )
179      !!---------------------------------------------------------------------
180      !!                  ***  ROUTINE glo_dyn_wri  ***
181      !!
182      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
183      !!----------------------------------------------------------------------
184      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
185      INTEGER, INTENT(in) ::   Kmm  ! time level index
186      !
187      INTEGER  ::   ji, jj, jk   ! dummy loop indices
188      REAL(dp) ::   zcof         ! local scalar
189      REAL(dp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
190      !!----------------------------------------------------------------------
191
192      ! I. Momentum trends
193      ! -------------------
194
195      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
196
197         ! I.1 Conversion potential energy - kinetic energy
198         ! --------------------------------------------------
199         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
200         zkx  (:,:,:) = 0._wp
201         zky  (:,:,:) = 0._wp
202         zkz  (:,:,:) = 0._wp
203         zkepe(:,:,:) = 0._wp
204   
205         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
206
207         zcof = 0.5_wp / rho0             ! Density flux at w-point
208         zkz(:,:,1) = 0._wp
209         DO jk = 2, jpk
210            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
211         END DO
212         
213         zcof   = 0.5_wp / rho0           ! Density flux at u and v-points
214         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
215            zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
216               &                              *  uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
217            zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
218               &                              *  vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
219         END_3D
220         
221         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Density flux divergence at t-point
222            zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
223               &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
224               &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
225               &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
226         END_3D
227
228         ! I.2 Basin averaged kinetic energy trend
229         ! ----------------------------------------
230         peke = 0._wp
231         DO jk = 1, jpkm1
232            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:)   &
233               &                               * e3t(:,:,jk,Kmm) )
234         END DO
235         peke = grav * peke
236
237         ! I.3 Sums over the global domain
238         ! ---------------------------------
239         IF( lk_mpp ) THEN
240            CALL mpp_sum( 'trdglo', peke )
241            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
242            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
243            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
244         ENDIF
245
246         ! I.2 Print dynamic trends in the ocean.output file
247         ! --------------------------------------------------
248
249         IF(lwp) THEN
250            WRITE (numout,*)
251            WRITE (numout,*)
252            WRITE (numout,9500) kt
253            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
254            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
255            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
256            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
257            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
258            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
259            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
260            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
261            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
262            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
263            WRITE (numout,9511)
264            WRITE (numout,9512)                                                 &
265            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
266            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
267            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
268            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
269            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
270            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
271            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
272!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
273         ENDIF
274
275 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
276 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
277 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
278 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
279 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
280 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
281 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
282 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
283 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
284 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
285 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
286 9511    FORMAT(' -----------------------------------------------------------------------------')
287 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
288 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
289 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
290
291         IF(lwp) THEN
292            WRITE (numout,*)
293            WRITE (numout,*)
294            WRITE (numout,9520) kt
295            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
296            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
297            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
298            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
299            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
300            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
301            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
302            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
303            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
304            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
305            WRITE (numout,9531)
306            WRITE (numout,9532)   &
307            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
308            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
309            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
310            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
311!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
312         ENDIF
313
314 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
315 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
316 9522    FORMAT(' ke gradient               u2= ', e20.13)
317 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
318 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
319 9525    FORMAT(' vertical advection        u2= ', e20.13)
320 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
321 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
322 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
323 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
324 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
325 9531    FORMAT(' --------------------------------------------------')
326 9532    FORMAT(' total trend               u2= ', e20.13)
327 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
328 9534    FORMAT('       bottom stress       u2= ', e20.13)
329
330         IF(lwp) THEN
331            WRITE (numout,*)
332            WRITE (numout,*)
333            WRITE (numout,9540) kt
334            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
335            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
336            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
337            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
338            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
339            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
340            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
341            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
342            WRITE (numout,*)
343            WRITE (numout,*)
344         ENDIF
345
346 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
347 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
348 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
349 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
350 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
351 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
352 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
353 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
354 9548    FORMAT(' pressure gradient u2 = - 1/rho0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
355         !
356         ! Save potential to kinetic energy conversion for next time step
357         rpktrd = peke
358         !
359      ENDIF
360      !
361   END SUBROUTINE glo_dyn_wri
362
363
364   SUBROUTINE glo_tra_wri( kt )
365      !!---------------------------------------------------------------------
366      !!                  ***  ROUTINE glo_tra_wri  ***
367      !!
368      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
369      !!----------------------------------------------------------------------
370      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
371      !
372      INTEGER  ::   jk   ! loop indices
373      !!----------------------------------------------------------------------
374
375      ! I. Tracers trends
376      ! -----------------
377
378      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
379
380         ! I.1 Sums over the global domain
381         ! -------------------------------
382         IF( lk_mpp ) THEN
383            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
384            CALL mpp_sum( 'trdglo', smo, jptot_tra )
385            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
386            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
387         ENDIF
388
389         ! I.2 Print tracers trends in the ocean.output file
390         ! --------------------------------------------------
391         
392         IF(lwp) THEN
393            WRITE (numout,*)
394            WRITE (numout,*)
395            WRITE (numout,9400) kt
396            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
397            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
398            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
399            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
400            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
401            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
402            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
403            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
404            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
405            WRITE (numout,9409) 
406            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
407            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
408            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
409            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
410         ENDIF
411
4129400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
413              '              salinity',/' ============================')
4149401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4159411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4169402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4179403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4189404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4199405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4209406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4219407     FORMAT(' penetrative qsr             ',e20.13)
4229408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4239409     FORMAT(' -------------------------------------------------------------------------')
4249410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
425
426
427         IF(lwp) THEN
428            WRITE (numout,*)
429            WRITE (numout,*)
430            WRITE (numout,9420) kt
431            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
432            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
433            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
434            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
435            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
436            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
437            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
438            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
439            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
440            WRITE (numout,9429)
441            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
442            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
443            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
444            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
445         ENDIF
446
4479420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
448            '               salinity', /, ' ===============================')
4499421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4509431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4519422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4529423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4539424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4549425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4559426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4569427     FORMAT(' penetrative qsr           * t   ', e20.13)
4579428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4589429     FORMAT(' -----------------------------------------------------------------------------')
4599430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
460
461
462         IF(lwp) THEN
463            WRITE (numout,*)
464            WRITE (numout,*)
465            WRITE (numout,9440) kt
466            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
467            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
468            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
469            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
470            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
471            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
472            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
473            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
474            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
475            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
476            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
477         ENDIF
478
4799440     FORMAT(' tracer consistency at it= ',i6,   &
480            ' :         temperature','                salinity',/,   &
481            ' ==================================')
4829441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4839442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4849443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4859444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4869445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4879446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4889447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
4899448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
4909449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
491         !
492      ENDIF
493      !
494   END SUBROUTINE glo_tra_wri
495
496
497   SUBROUTINE trd_glo_init( Kmm )
498      !!---------------------------------------------------------------------
499      !!                  ***  ROUTINE trd_glo_init  ***
500      !!
501      !! ** Purpose :   Read the namtrd namelist
502      !!----------------------------------------------------------------------
503      INTEGER, INTENT(in) ::   Kmm   ! time level index
504      INTEGER  ::   ji, jj, jk   ! dummy loop indices
505      !!----------------------------------------------------------------------
506
507      IF(lwp) THEN
508         WRITE(numout,*)
509         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
510         WRITE(numout,*) '~~~~~~~~~~~~~'
511      ENDIF
512
513      ! Total volume at t-points:
514      tvolt = 0._wp
515      DO jk = 1, jpkm1
516         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
517      END DO
518      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
519
520      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
521
522      ! Initialization of potential to kinetic energy conversion
523      rpktrd = 0._wp
524
525      ! Total volume at u-, v- points:
526!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
527      tvolu = 0._wp
528      tvolv = 0._wp
529
530      DO_3D( 0, 0, 0, 0, 1, jpk )
531         tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)   &
532            &                                       * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
533         tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   &
534            &                                       * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
535      END_3D
536      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
537      CALL mpp_sum( 'trdglo', tvolv )
538
539      IF(lwp) THEN
540         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
541         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
542      ENDIF
543      !
544   END SUBROUTINE trd_glo_init
545
546   !!======================================================================
547END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.