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/dev_r9922_GC3_cpl_pkg/src/OCE/TRD – NEMO

source: NEMO/branches/UKMO/dev_r9922_GC3_cpl_pkg/src/OCE/TRD/trdglo.F90 @ 9947

Last change on this file since 9947 was 9947, checked in by timgraham, 6 years ago

Clear svn keywords

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(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 licence     (./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( peke )
252            CALL mpp_sum( umo , jptot_dyn )
253            CALL mpp_sum( vmo , jptot_dyn )
254            CALL mpp_sum( 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         ENDIF
285
286 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
287 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
288 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
289 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
290 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
291 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
292 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
293 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
294 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
295 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
296 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
297 9511    FORMAT(' -----------------------------------------------------------------------------')
298 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
299 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
300 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
301
302         IF(lwp) THEN
303            WRITE (numout,*)
304            WRITE (numout,*)
305            WRITE (numout,9520) kt
306            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
307            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
308            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
309            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
310            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
311            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
312            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
313            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
314            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
315            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
316            WRITE (numout,9531)
317            WRITE (numout,9532)   &
318            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
319            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
320            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
321            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
322!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
323         ENDIF
324
325 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
326 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
327 9522    FORMAT(' ke gradient               u2= ', e20.13)
328 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
329 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
330 9525    FORMAT(' vertical advection        u2= ', e20.13)
331 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
332 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
333 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
334 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
335 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
336 9531    FORMAT(' --------------------------------------------------')
337 9532    FORMAT(' total trend               u2= ', e20.13)
338 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
339 9534    FORMAT('       bottom stress       u2= ', e20.13)
340
341         IF(lwp) THEN
342            WRITE (numout,*)
343            WRITE (numout,*)
344            WRITE (numout,9540) kt
345            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
346            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
347            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
348            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
349            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
350            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
351            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
352            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
353            WRITE (numout,*)
354            WRITE (numout,*)
355         ENDIF
356
357 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
358 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
359 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
360 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
361 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
362 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
363 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
364 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
365 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
366         !
367         ! Save potential to kinetic energy conversion for next time step
368         rpktrd = peke
369         !
370      ENDIF
371      !
372   END SUBROUTINE glo_dyn_wri
373
374
375   SUBROUTINE glo_tra_wri( kt )
376      !!---------------------------------------------------------------------
377      !!                  ***  ROUTINE glo_tra_wri  ***
378      !!
379      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
380      !!----------------------------------------------------------------------
381      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
382      !
383      INTEGER  ::   jk   ! loop indices
384      !!----------------------------------------------------------------------
385
386      ! I. Tracers trends
387      ! -----------------
388
389      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
390
391         ! I.1 Sums over the global domain
392         ! -------------------------------
393         IF( lk_mpp ) THEN
394            CALL mpp_sum( tmo, jptot_tra )   
395            CALL mpp_sum( smo, jptot_tra )
396            CALL mpp_sum( t2 , jptot_tra )
397            CALL mpp_sum( s2 , jptot_tra )
398         ENDIF
399
400         ! I.2 Print tracers trends in the ocean.output file
401         ! --------------------------------------------------
402         
403         IF(lwp) THEN
404            WRITE (numout,*)
405            WRITE (numout,*)
406            WRITE (numout,9400) kt
407            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
408            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
409            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
410            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
411            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
412            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
413            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
414            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
415            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
416            WRITE (numout,9409) 
417            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
418            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
419            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
420            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
421         ENDIF
422
4239400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
424              '              salinity',/' ============================')
4259401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4269411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4279402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4289403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4299404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4309405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4319406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4329407     FORMAT(' penetrative qsr             ',e20.13)
4339408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4349409     FORMAT(' -------------------------------------------------------------------------')
4359410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
436
437
438         IF(lwp) THEN
439            WRITE (numout,*)
440            WRITE (numout,*)
441            WRITE (numout,9420) kt
442            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
443            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
444            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
445            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
446            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
447            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
448            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
449            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
450            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
451            WRITE (numout,9429)
452            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
453            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
454            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
455            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
456         ENDIF
457
4589420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
459            '               salinity', /, ' ===============================')
4609421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4619431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4629422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4639423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4649424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4659425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4669426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4679427     FORMAT(' penetrative qsr           * t   ', e20.13)
4689428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4699429     FORMAT(' -----------------------------------------------------------------------------')
4709430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
471
472
473         IF(lwp) THEN
474            WRITE (numout,*)
475            WRITE (numout,*)
476            WRITE (numout,9440) kt
477            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
478            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
479            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
480            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
481            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
482            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
483            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
484            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
485            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
486            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
487            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
488         ENDIF
489
4909440     FORMAT(' tracer consistency at it= ',i6,   &
491            ' :         temperature','                salinity',/,   &
492            ' ==================================')
4939441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4949442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4959443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4969444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4979445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
4989446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
4999447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5009448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5019449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
502         !
503      ENDIF
504      !
505   END SUBROUTINE glo_tra_wri
506
507
508   SUBROUTINE trd_glo_init
509      !!---------------------------------------------------------------------
510      !!                  ***  ROUTINE trd_glo_init  ***
511      !!
512      !! ** Purpose :   Read the namtrd namelist
513      !!----------------------------------------------------------------------
514      INTEGER  ::   ji, jj, jk   ! dummy loop indices
515      !!----------------------------------------------------------------------
516
517      IF(lwp) THEN
518         WRITE(numout,*)
519         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
520         WRITE(numout,*) '~~~~~~~~~~~~~'
521      ENDIF
522
523      ! Total volume at t-points:
524      tvolt = 0._wp
525      DO jk = 1, jpkm1
526         tvolt = tvolt + SUM( e1e2t(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
527      END DO
528      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
529
530      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
531
532      ! Initialization of potential to kinetic energy conversion
533      rpktrd = 0._wp
534
535      ! Total volume at u-, v- points:
536!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
537      tvolu = 0._wp
538      tvolv = 0._wp
539
540      DO jk = 1, jpk
541         DO jj = 2, jpjm1
542            DO ji = fs_2, fs_jpim1   ! vector opt.
543               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)
544               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)
545            END DO
546         END DO
547      END DO
548      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
549      IF( lk_mpp )   CALL mpp_sum( tvolv )
550
551      IF(lwp) THEN
552         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
553         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
554      ENDIF
555      !
556   END SUBROUTINE trd_glo_init
557
558   !!======================================================================
559END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.