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/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/TRD – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/TRD/trdglo.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush
File size: 29.1 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   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
62      !!---------------------------------------------------------------------
63      !!                  ***  ROUTINE trd_glo  ***
64      !!
65      !! ** Purpose :   compute and print global domain averaged trends for
66      !!              T, T^2, momentum, KE, and KE<->PE
67      !!
68      !!----------------------------------------------------------------------
69      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
70      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
71      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
72      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
73      INTEGER                   , INTENT(in   ) ::   kt      ! time step
74      !!
75      INTEGER ::   ji, jj, jk      ! dummy loop indices
76      INTEGER ::   ikbu, ikbv      ! local integers
77      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
78      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
79      !!----------------------------------------------------------------------
80      !
81      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
82         !
83         SELECT CASE( ctype )
84         !
85         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
86            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
87               DO jj = 1, jpj
88                  DO ji = 1, jpi       
89                     zvm = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * 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 * tsn(ji,jj,jk,jp_tem)
95                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
96                  END DO
97               END DO
98            END DO
99            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
100            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
101               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:) )
102               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
103               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:) * tmask_i(:,:)  )
104               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:) * tmask_i(:,:)  )
105            ENDIF
106            !
107            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
108               !
109               CALL glo_tra_wri( kt )             ! print the results in ocean.output
110               !               
111               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
112               smo(:) = 0._wp
113               t2 (:) = 0._wp
114               s2 (:) = 0._wp
115               !
116            ENDIF
117            !
118         CASE( 'DYN' )          !==  Momentum and KE  ==!       
119            DO jk = 1, jpkm1
120               DO jj = 1, jpjm1
121                  DO ji = 1, jpim1
122                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
123                        &                                     * e1e2u  (ji,jj) * e3u_n(ji,jj,jk)
124                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
125                        &                                     * e1e2v  (ji,jj) * e3u_n(ji,jj,jk)
126                     umo(ktrd) = umo(ktrd) + zvt
127                     vmo(ktrd) = vmo(ktrd) + zvs
128                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
129                  END DO
130               END DO
131            END DO
132            !                 
133            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
134               z1_2rau0 = 0.5_wp / rau0
135               DO jj = 1, jpjm1
136                  DO ji = 1, jpim1
137                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
138                        &                                                     * z1_2rau0       * e1e2u(ji,jj)
139                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
140                        &                                                     * z1_2rau0       * e1e2v(ji,jj)
141                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
142                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
143                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
144                  END DO
145               END DO
146            ENDIF
147            !                       
148!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
149!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
150!               !
151!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
152!                  z1_2rau0 = 0.5_wp / rau0
153!                  DO jj = 1, jpjm1
154!                     DO ji = 1, jpim1
155!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
156!                        ikbv = mbkv(ji,jj)
157!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * un(ji,jj,ikbu) * e1e2u(ji,jj)
158!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vn(ji,jj,ikbv) * e1e2v(ji,jj)
159!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
160!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
161!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
162!                     END DO
163!                  END DO
164!               ENDIF
165!
166!!gm top drag case is missing
167!
168!               !
169!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
170!               !               
171!               umo(:) = 0._wp                         ! reset for the next time step
172!               vmo(:) = 0._wp
173!               hke(:) = 0._wp
174!               !
175!            ENDIF
176!!gm end
177            !
178         END SELECT
179         !
180      ENDIF
181      !
182   END SUBROUTINE trd_glo
183
184
185   SUBROUTINE glo_dyn_wri( kt )
186      !!---------------------------------------------------------------------
187      !!                  ***  ROUTINE glo_dyn_wri  ***
188      !!
189      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
190      !!----------------------------------------------------------------------
191      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
192      !
193      INTEGER  ::   ji, jj, jk   ! dummy loop indices
194      REAL(wp) ::   zcof         ! local scalar
195      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
196      !!----------------------------------------------------------------------
197
198      ! I. Momentum trends
199      ! -------------------
200
201      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
202
203         ! I.1 Conversion potential energy - kinetic energy
204         ! --------------------------------------------------
205         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
206         zkx  (:,:,:) = 0._wp
207         zky  (:,:,:) = 0._wp
208         zkz  (:,:,:) = 0._wp
209         zkepe(:,:,:) = 0._wp
210   
211         CALL eos( tsn, rhd, rhop )       ! now potential density
212
213         zcof = 0.5_wp / rau0             ! Density flux at w-point
214         zkz(:,:,1) = 0._wp
215         DO jk = 2, jpk
216            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
217         END DO
218         
219         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
220         DO jk = 1, jpkm1
221            DO jj = 1, jpjm1
222               DO ji = 1, jpim1
223                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
224                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
225               END DO
226            END DO
227         END DO
228         
229         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
230            DO jj = 2, jpjm1
231               DO ji = 2, jpim1
232                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
233                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
234                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
235                     &              / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
236               END DO
237            END DO
238         END DO
239
240         ! I.2 Basin averaged kinetic energy trend
241         ! ----------------------------------------
242         peke = 0._wp
243         DO jk = 1, jpkm1
244            peke = peke + SUM( zkepe(:,:,jk) * gdept_n(:,:,jk) * e1e2t(:,:) * e3t_n(:,:,jk) )
245         END DO
246         peke = grav * peke
247
248         ! I.3 Sums over the global domain
249         ! ---------------------------------
250         IF( lk_mpp ) THEN
251            CALL mpp_sum( 'trdglo', peke )
252            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
253            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
254            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
255         ENDIF
256
257         ! I.2 Print dynamic trends in the ocean.output file
258         ! --------------------------------------------------
259
260         IF(lwp) THEN
261            WRITE (numout,*)
262            WRITE (numout,*)
263            WRITE (numout,9500) kt
264            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
265            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
266            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
267            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
268            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
269            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
270            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
271            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
272            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
273            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
274            WRITE (numout,9511)
275            WRITE (numout,9512)                                                 &
276            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
277            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
278            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
279            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
280            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
281            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
282            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
283!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
284            IF(lflush) CALL FLUSH(numout)
285         ENDIF
286
287 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
288 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
289 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
290 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
291 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
292 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
293 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
294 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
295 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
296 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
297 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
298 9511    FORMAT(' -----------------------------------------------------------------------------')
299 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
300 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
301 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
302
303         IF(lwp) THEN
304            WRITE (numout,*)
305            WRITE (numout,*)
306            WRITE (numout,9520) kt
307            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
308            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
309            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
310            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
311            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
312            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
313            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
314            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
315            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
316            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
317            WRITE (numout,9531)
318            WRITE (numout,9532)   &
319            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
320            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
321            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
322            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
323!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
324            IF(lflush) CALL FLUSH(numout)
325         ENDIF
326
327 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
328 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
329 9522    FORMAT(' ke gradient               u2= ', e20.13)
330 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
331 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
332 9525    FORMAT(' vertical advection        u2= ', e20.13)
333 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
334 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
335 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
336 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
337 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
338 9531    FORMAT(' --------------------------------------------------')
339 9532    FORMAT(' total trend               u2= ', e20.13)
340 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
341 9534    FORMAT('       bottom stress       u2= ', e20.13)
342
343         IF(lwp) THEN
344            WRITE (numout,*)
345            WRITE (numout,*)
346            WRITE (numout,9540) kt
347            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
348            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
349            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
350            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
351            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
352            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
353            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
354            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
355            WRITE (numout,*)
356            WRITE (numout,*)
357            IF(lflush) CALL FLUSH(numout)
358         ENDIF
359
360 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
361 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
362 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
363 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
364 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
365 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
366 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
367 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
368 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
369         !
370         ! Save potential to kinetic energy conversion for next time step
371         rpktrd = peke
372         !
373      ENDIF
374      !
375   END SUBROUTINE glo_dyn_wri
376
377
378   SUBROUTINE glo_tra_wri( kt )
379      !!---------------------------------------------------------------------
380      !!                  ***  ROUTINE glo_tra_wri  ***
381      !!
382      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
383      !!----------------------------------------------------------------------
384      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
385      !
386      INTEGER  ::   jk   ! loop indices
387      !!----------------------------------------------------------------------
388
389      ! I. Tracers trends
390      ! -----------------
391
392      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
393
394         ! I.1 Sums over the global domain
395         ! -------------------------------
396         IF( lk_mpp ) THEN
397            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
398            CALL mpp_sum( 'trdglo', smo, jptot_tra )
399            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
400            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
401         ENDIF
402
403         ! I.2 Print tracers trends in the ocean.output file
404         ! --------------------------------------------------
405         
406         IF(lwp) THEN
407            WRITE (numout,*)
408            WRITE (numout,*)
409            WRITE (numout,9400) kt
410            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
411            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
412            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
413            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
414            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
415            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
416            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
417            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
418            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
419            WRITE (numout,9409) 
420            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
421            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
422            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
423            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
424            IF(lflush) CALL FLUSH(numout)
425         ENDIF
426
4279400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
428              '              salinity',/' ============================')
4299401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4309411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4319402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4329403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4339404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4349405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4359406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4369407     FORMAT(' penetrative qsr             ',e20.13)
4379408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4389409     FORMAT(' -------------------------------------------------------------------------')
4399410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
440
441
442         IF(lwp) THEN
443            WRITE (numout,*)
444            WRITE (numout,*)
445            WRITE (numout,9420) kt
446            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
447            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
448            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
449            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
450            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
451            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
452            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
453            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
454            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
455            WRITE (numout,9429)
456            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
457            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
458            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
459            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
460            IF(lflush) CALL FLUSH(numout)
461         ENDIF
462
4639420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
464            '               salinity', /, ' ===============================')
4659421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4669431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4679422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4689423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4699424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4709425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4719426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4729427     FORMAT(' penetrative qsr           * t   ', e20.13)
4739428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4749429     FORMAT(' -----------------------------------------------------------------------------')
4759430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
476
477
478         IF(lwp) THEN
479            WRITE (numout,*)
480            WRITE (numout,*)
481            WRITE (numout,9440) kt
482            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
483            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
484            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
485            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
486            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
487            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
488            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
489            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
490            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
491            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
492            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
493            IF(lflush) CALL FLUSH(numout)
494         ENDIF
495
4969440     FORMAT(' tracer consistency at it= ',i6,   &
497            ' :         temperature','                salinity',/,   &
498            ' ==================================')
4999441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5009442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5019443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5029444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5039445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5049446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5059447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5069448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5079449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
508         !
509      ENDIF
510      !
511   END SUBROUTINE glo_tra_wri
512
513
514   SUBROUTINE trd_glo_init
515      !!---------------------------------------------------------------------
516      !!                  ***  ROUTINE trd_glo_init  ***
517      !!
518      !! ** Purpose :   Read the namtrd namelist
519      !!----------------------------------------------------------------------
520      INTEGER  ::   ji, jj, jk   ! dummy loop indices
521      !!----------------------------------------------------------------------
522
523      IF(lwp) THEN
524         WRITE(numout,*)
525         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
526         WRITE(numout,*) '~~~~~~~~~~~~~'
527         IF(lflush) CALL FLUSH(numout)
528      ENDIF
529
530      ! Total volume at t-points:
531      tvolt = 0._wp
532      DO jk = 1, jpkm1
533         tvolt = tvolt + SUM( e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
534      END DO
535      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
536
537      IF(lwp) THEN
538         WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
539         IF(lflush) CALL FLUSH(numout)
540      ENDIF
541
542      ! Initialization of potential to kinetic energy conversion
543      rpktrd = 0._wp
544
545      ! Total volume at u-, v- points:
546!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
547      tvolu = 0._wp
548      tvolv = 0._wp
549
550      DO jk = 1, jpk
551         DO jj = 2, jpjm1
552            DO ji = fs_2, fs_jpim1   ! vector opt.
553               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u_n(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
554               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v_n(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
555            END DO
556         END DO
557      END DO
558      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
559      CALL mpp_sum( 'trdglo', tvolv )
560
561      IF(lwp) THEN
562         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
563         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
564         IF(lflush) CALL FLUSH(numout)
565      ENDIF
566      !
567   END SUBROUTINE trd_glo_init
568
569   !!======================================================================
570END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.