New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdglo.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdglo.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 29.0 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, Kmm )
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      INTEGER                   , INTENT(in   ) ::   Kmm     ! time level index
75      !!
76      INTEGER ::   ji, jj, jk      ! dummy loop indices
77      INTEGER ::   ikbu, ikbv      ! local integers
78      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
79      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
80      !!----------------------------------------------------------------------
81      !
82      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
83         !
84         SELECT CASE( ctype )
85         !
86         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
87            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
88               DO jj = 1, jpj
89                  DO ji = 1, jpi       
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 DO
98               END DO
99            END DO
100            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
101            IF( ln_linssh .AND. ktrd == jptra_zad ) THEN 
102               tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
103               smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
104               t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
105               s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:)  )
106            ENDIF
107            !
108            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
109               !
110               CALL glo_tra_wri( kt )             ! print the results in ocean.output
111               !               
112               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
113               smo(:) = 0._wp
114               t2 (:) = 0._wp
115               s2 (:) = 0._wp
116               !
117            ENDIF
118            !
119         CASE( 'DYN' )          !==  Momentum and KE  ==!       
120            DO jk = 1, jpkm1
121               DO jj = 1, jpjm1
122                  DO ji = 1, jpim1
123                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
124                        &                                     * e1e2u  (ji,jj) * e3u(ji,jj,jk,Kmm)
125                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
126                        &                                     * e1e2v  (ji,jj) * e3u(ji,jj,jk,Kmm)
127                     umo(ktrd) = umo(ktrd) + zvt
128                     vmo(ktrd) = vmo(ktrd) + zvs
129                     hke(ktrd) = hke(ktrd) + uu(ji,jj,jk,Kmm) * zvt + vv(ji,jj,jk,Kmm) * zvs
130                  END DO
131               END DO
132            END DO
133            !                 
134            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
135               z1_2rau0 = 0.5_wp / rau0
136               DO jj = 1, jpjm1
137                  DO ji = 1, jpim1
138                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
139                        &                                                     * z1_2rau0       * e1e2u(ji,jj)
140                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
141                        &                                                     * z1_2rau0       * e1e2v(ji,jj)
142                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
143                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
144                     hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
145                  END DO
146               END DO
147            ENDIF
148            !                       
149!!gm  miss placed calculation   ===>>>> to be done in dynzdf.F90
150!            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
151!               !
152!               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction
153!                  z1_2rau0 = 0.5_wp / rau0
154!                  DO jj = 1, jpjm1
155!                     DO ji = 1, jpim1
156!                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
157!                        ikbv = mbkv(ji,jj)
158!                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu(ji,jj,ikbu,Kmm) * e1e2u(ji,jj)
159!                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv(ji,jj,ikbv,Kmm) * e1e2v(ji,jj)
160!                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
161!                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
162!                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + uu(ji,jj,ikbu,Kmm) * zvt + vv(ji,jj,ikbv,Kmm) * zvs
163!                     END DO
164!                  END DO
165!               ENDIF
166!
167!!gm top drag case is missing
168!
169!               !
170!               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
171!               !               
172!               umo(:) = 0._wp                         ! reset for the next time step
173!               vmo(:) = 0._wp
174!               hke(:) = 0._wp
175!               !
176!            ENDIF
177!!gm end
178            !
179         END SELECT
180         !
181      ENDIF
182      !
183   END SUBROUTINE trd_glo
184
185
186   SUBROUTINE glo_dyn_wri( kt, Kmm )
187      !!---------------------------------------------------------------------
188      !!                  ***  ROUTINE glo_dyn_wri  ***
189      !!
190      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
191      !!----------------------------------------------------------------------
192      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
193      INTEGER, INTENT(in) ::   Kmm  ! time level index
194      !
195      INTEGER  ::   ji, jj, jk   ! dummy loop indices
196      REAL(wp) ::   zcof         ! local scalar
197      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe 
198      !!----------------------------------------------------------------------
199
200      ! I. Momentum trends
201      ! -------------------
202
203      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
204
205         ! I.1 Conversion potential energy - kinetic energy
206         ! --------------------------------------------------
207         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
208         zkx  (:,:,:) = 0._wp
209         zky  (:,:,:) = 0._wp
210         zkz  (:,:,:) = 0._wp
211         zkepe(:,:,:) = 0._wp
212   
213         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
214
215         zcof = 0.5_wp / rau0             ! Density flux at w-point
216         zkz(:,:,1) = 0._wp
217         DO jk = 2, jpk
218            zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
219         END DO
220         
221         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
222         DO jk = 1, jpkm1
223            DO jj = 1, jpjm1
224               DO ji = 1, jpim1
225                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
226                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
227               END DO
228            END DO
229         END DO
230         
231         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
232            DO jj = 2, jpjm1
233               DO ji = 2, jpim1
234                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
235                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
236                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
237                     &              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
238               END DO
239            END DO
240         END DO
241
242         ! I.2 Basin averaged kinetic energy trend
243         ! ----------------------------------------
244         peke = 0._wp
245         DO jk = 1, jpkm1
246            peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) * e3t(:,:,jk,Kmm) )
247         END DO
248         peke = grav * peke
249
250         ! I.3 Sums over the global domain
251         ! ---------------------------------
252         IF( lk_mpp ) THEN
253            CALL mpp_sum( 'trdglo', peke )
254            CALL mpp_sum( 'trdglo', umo , jptot_dyn )
255            CALL mpp_sum( 'trdglo', vmo , jptot_dyn )
256            CALL mpp_sum( 'trdglo', hke , jptot_dyn )
257         ENDIF
258
259         ! I.2 Print dynamic trends in the ocean.output file
260         ! --------------------------------------------------
261
262         IF(lwp) THEN
263            WRITE (numout,*)
264            WRITE (numout,*)
265            WRITE (numout,9500) kt
266            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
267            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
268            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
269            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
270            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
271            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
272            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
273            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
274            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
275            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
276            WRITE (numout,9511)
277            WRITE (numout,9512)                                                 &
278            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
279            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
280            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
281            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
282            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
283            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
284            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
285!!gm            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
286         ENDIF
287
288 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
289 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
290 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
291 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
292 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
293 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
294 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
295 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
296 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
297 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
298 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
299 9511    FORMAT(' -----------------------------------------------------------------------------')
300 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
301 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
302 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
303
304         IF(lwp) THEN
305            WRITE (numout,*)
306            WRITE (numout,*)
307            WRITE (numout,9520) kt
308            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
309            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
310            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
311            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
312            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
313            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
314            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
315            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
316            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
317            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
318            WRITE (numout,9531)
319            WRITE (numout,9532)   &
320            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
321            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
322            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
323            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
324!!gm            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
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         ENDIF
358
359 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
360 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
361 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
362 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
363 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
364 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
365 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
366 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
367 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
368         !
369         ! Save potential to kinetic energy conversion for next time step
370         rpktrd = peke
371         !
372      ENDIF
373      !
374   END SUBROUTINE glo_dyn_wri
375
376
377   SUBROUTINE glo_tra_wri( kt )
378      !!---------------------------------------------------------------------
379      !!                  ***  ROUTINE glo_tra_wri  ***
380      !!
381      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
382      !!----------------------------------------------------------------------
383      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
384      !
385      INTEGER  ::   jk   ! loop indices
386      !!----------------------------------------------------------------------
387
388      ! I. Tracers trends
389      ! -----------------
390
391      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
392
393         ! I.1 Sums over the global domain
394         ! -------------------------------
395         IF( lk_mpp ) THEN
396            CALL mpp_sum( 'trdglo', tmo, jptot_tra )   
397            CALL mpp_sum( 'trdglo', smo, jptot_tra )
398            CALL mpp_sum( 'trdglo', t2 , jptot_tra )
399            CALL mpp_sum( 'trdglo', s2 , jptot_tra )
400         ENDIF
401
402         ! I.2 Print tracers trends in the ocean.output file
403         ! --------------------------------------------------
404         
405         IF(lwp) THEN
406            WRITE (numout,*)
407            WRITE (numout,*)
408            WRITE (numout,9400) kt
409            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
410            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
411            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
412            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
413            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
414            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
415            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
416            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
417            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
418            WRITE (numout,9409) 
419            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
420            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
421            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
422            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
423         ENDIF
424
4259400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
426              '              salinity',/' ============================')
4279401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4289411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4299402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4309403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4319404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4329405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4339406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4349407     FORMAT(' penetrative qsr             ',e20.13)
4359408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4369409     FORMAT(' -------------------------------------------------------------------------')
4379410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
438
439
440         IF(lwp) THEN
441            WRITE (numout,*)
442            WRITE (numout,*)
443            WRITE (numout,9420) kt
444            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
445            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
446            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
447            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
448            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
449            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
450            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
451            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
452            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
453            WRITE (numout,9429)
454            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
455            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
456            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
457            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
458         ENDIF
459
4609420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
461            '               salinity', /, ' ===============================')
4629421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4639431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4649422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4659423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4669424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4679425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4689426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4699427     FORMAT(' penetrative qsr           * t   ', e20.13)
4709428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4719429     FORMAT(' -----------------------------------------------------------------------------')
4729430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
473
474
475         IF(lwp) THEN
476            WRITE (numout,*)
477            WRITE (numout,*)
478            WRITE (numout,9440) kt
479            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
480            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
481            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
482            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
483            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
484            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
485            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
486            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
487            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
488            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
489            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
490         ENDIF
491
4929440     FORMAT(' tracer consistency at it= ',i6,   &
493            ' :         temperature','                salinity',/,   &
494            ' ==================================')
4959441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
4969442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
4979443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
4989444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
4999445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5009446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5019447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5029448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5039449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
504         !
505      ENDIF
506      !
507   END SUBROUTINE glo_tra_wri
508
509
510   SUBROUTINE trd_glo_init( Kmm )
511      !!---------------------------------------------------------------------
512      !!                  ***  ROUTINE trd_glo_init  ***
513      !!
514      !! ** Purpose :   Read the namtrd namelist
515      !!----------------------------------------------------------------------
516      INTEGER, INTENT(in) ::   Kmm   ! time level index
517      INTEGER  ::   ji, jj, jk   ! dummy loop indices
518      !!----------------------------------------------------------------------
519
520      IF(lwp) THEN
521         WRITE(numout,*)
522         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
523         WRITE(numout,*) '~~~~~~~~~~~~~'
524      ENDIF
525
526      ! Total volume at t-points:
527      tvolt = 0._wp
528      DO jk = 1, jpkm1
529         tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
530      END DO
531      CALL mpp_sum( 'trdglo', tvolt )   ! sum over the global domain
532
533      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
534
535      ! Initialization of potential to kinetic energy conversion
536      rpktrd = 0._wp
537
538      ! Total volume at u-, v- points:
539!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
540      tvolu = 0._wp
541      tvolv = 0._wp
542
543      DO jk = 1, jpk
544         DO jj = 2, jpjm1
545            DO ji = fs_2, fs_jpim1   ! vector opt.
546               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
547               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
548            END DO
549         END DO
550      END DO
551      CALL mpp_sum( 'trdglo', tvolu )   ! sums over the global domain
552      CALL mpp_sum( 'trdglo', tvolv )
553
554      IF(lwp) THEN
555         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
556         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
557      ENDIF
558      !
559   END SUBROUTINE trd_glo_init
560
561   !!======================================================================
562END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.