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 branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 ​svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package
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_oce      ! ocean active tracers: lateral physics
22   USE ldfdyn_oce      ! ocean dynamics: lateral physics
23   USE zdf_oce         ! ocean vertical physics
24   USE zdfbfr          ! bottom friction
25   USE zdfddm          ! ocean vertical physics: double diffusion
26   USE eosbn2          ! equation of state
27   USE phycst          ! physical constants
28   USE lib_mpp         ! distibuted memory computing library
29   USE in_out_manager  ! I/O manager
30   USE iom             ! I/O manager library
31   USE wrk_nemo        ! Memory allocation
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 "domzgr_substitute.h90"
54#  include "vectopt_loop_substitute.h90"
55#  include "zdfddm_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE trd_glo  ***
66      !!
67      !! ** Purpose :   compute and print global domain averaged trends for
68      !!              T, T^2, momentum, KE, and KE<->PE
69      !!
70      !!----------------------------------------------------------------------
71      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
73      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
74      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
75      INTEGER                   , INTENT(in   ) ::   kt      ! time step
76      !!
77      INTEGER ::   ji, jj, jk      ! dummy loop indices
78      INTEGER ::   ikbu, ikbv      ! local integers
79      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
80      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
81      !!----------------------------------------------------------------------
82
83      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
84
85      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
86         !
87         SELECT CASE( ctype )
88         !
89         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
90            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
91               DO jj = 1, jpj
92                  DO ji = 1, jpi       
93                     zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
94                     zvt = ptrdx(ji,jj,jk) * zvm
95                     zvs = ptrdy(ji,jj,jk) * zvm
96                     tmo(ktrd) = tmo(ktrd) + zvt   
97                     smo(ktrd) = smo(ktrd) + zvs
98                     t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
99                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
100                  END DO
101               END DO
102            END DO
103            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
104            IF( .NOT.lk_vvl .AND. ktrd == jptra_zad ) THEN 
105               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
106               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
107               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
108               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
109            ENDIF
110            !
111            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
112               !
113               CALL glo_tra_wri( kt )             ! print the results in ocean.output
114               !               
115               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
116               smo(:) = 0._wp
117               t2 (:) = 0._wp
118               s2 (:) = 0._wp
119               !
120            ENDIF
121            !
122         CASE( 'DYN' )          !==  Momentum and KE  ==!       
123            DO jk = 1, jpkm1
124               DO jj = 1, jpjm1
125                  DO ji = 1, jpim1
126                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
127                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
128                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
129                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
130                     umo(ktrd) = umo(ktrd) + zvt
131                     vmo(ktrd) = vmo(ktrd) + zvs
132                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
133                  END DO
134               END DO
135            END DO
136            !                 
137            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
138               z1_2rau0 = 0.5_wp / rau0
139               DO jj = 1, jpjm1
140                  DO ji = 1, jpim1
141                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
142                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
143                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
144                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
145                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
146                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
147                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
148                  END DO
149               END DO
150            ENDIF
151            !                       
152            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
153               !
154               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction
155                  z1_2rau0 = 0.5_wp / rau0
156                  DO jj = 1, jpjm1
157                     DO ji = 1, jpim1
158                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
159                        ikbv = mbkv(ji,jj)
160                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
161                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
162                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
163                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
164                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
165                     END DO
166                  END DO
167               ENDIF
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            !
177         END SELECT
178         !
179      ENDIF
180      !
181      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
182      !
183   END SUBROUTINE trd_glo
184
185
186   SUBROUTINE glo_dyn_wri( kt )
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      !
194      INTEGER  ::   ji, jj, jk   ! dummy loop indices
195      REAL(wp) ::   zcof         ! local scalar
196      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
197      !!----------------------------------------------------------------------
198
199      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
200
201      ! I. Momentum trends
202      ! -------------------
203
204      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
205
206         ! I.1 Conversion potential energy - kinetic energy
207         ! --------------------------------------------------
208         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
209         zkx  (:,:,:) = 0._wp
210         zky  (:,:,:) = 0._wp
211         zkz  (:,:,:) = 0._wp
212         zkepe(:,:,:) = 0._wp
213   
214         CALL eos( tsn, rhd, rhop )       ! now potential density
215
216         zcof = 0.5_wp / rau0             ! Density flux at w-point
217         zkz(:,:,1) = 0._wp
218         DO jk = 2, jpk
219            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
220         END DO
221         
222         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
223         DO jk = 1, jpkm1
224            DO jj = 1, jpjm1
225               DO ji = 1, jpim1
226                  zkx(ji,jj,jk) = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
227                  zky(ji,jj,jk) = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
228               END DO
229            END DO
230         END DO
231         
232         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
233            DO jj = 2, jpjm1
234               DO ji = 2, jpim1
235                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
236                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
237                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
238                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
239               END DO
240            END DO
241         END DO
242
243         ! I.2 Basin averaged kinetic energy trend
244         ! ----------------------------------------
245         peke = 0._wp
246         DO jk = 1, jpkm1
247            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
248         END DO
249         peke = grav * peke
250
251         ! I.3 Sums over the global domain
252         ! ---------------------------------
253         IF( lk_mpp ) THEN
254            CALL mpp_sum( peke )
255            CALL mpp_sum( umo , jptot_dyn )
256            CALL mpp_sum( vmo , jptot_dyn )
257            CALL mpp_sum( hke , jptot_dyn )
258         ENDIF
259
260         ! I.2 Print dynamic trends in the ocean.output file
261         ! --------------------------------------------------
262
263         IF(lwp) THEN
264            WRITE (numout,*)
265            WRITE (numout,*)
266            WRITE (numout,9500) kt
267            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
268            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
269            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
270            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
271            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
272            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
273            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
274            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
275            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
276            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
277            WRITE (numout,9511)
278            WRITE (numout,9512)                                                 &
279            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
280            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
281            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
282            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
283            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
284            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
285            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
286            IF( ln_bfrimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
287            IF(lflush) CALL flush(numout)
288         ENDIF
289
290 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
291 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
292 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
293 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
294 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
295 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
296 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
297 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
298 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
299 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
300 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
301 9511    FORMAT(' -----------------------------------------------------------------------------')
302 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
303 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
304 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
305
306         IF(lwp) THEN
307            WRITE (numout,*)
308            WRITE (numout,*)
309            WRITE (numout,9520) kt
310            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
311            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
312            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
313            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
314            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
315            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
316            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
317            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
318            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
319            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
320            WRITE (numout,9531)
321            WRITE (numout,9532)   &
322            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
323            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
324            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
325            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
326            IF( ln_bfrimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
327            IF(lflush) CALL flush(numout)
328         ENDIF
329
330 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
331 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
332 9522    FORMAT(' ke gradient               u2= ', e20.13)
333 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
334 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
335 9525    FORMAT(' vertical advection        u2= ', e20.13)
336 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
337 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
338 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
339 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
340 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
341 9531    FORMAT(' --------------------------------------------------')
342 9532    FORMAT(' total trend               u2= ', e20.13)
343 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
344 9534    FORMAT('       bottom stress       u2= ', e20.13)
345
346         IF(lwp) THEN
347            WRITE (numout,*)
348            WRITE (numout,*)
349            WRITE (numout,9540) kt
350            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
351            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
352            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
353            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
354            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
355            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
356            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
357            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
358            WRITE (numout,*)
359            WRITE (numout,*)
360            IF(lflush) CALL flush(numout)
361         ENDIF
362
363 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
364 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
365 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
366 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
367 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
368 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
369 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
370 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
371 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
372         !
373         ! Save potential to kinetic energy conversion for next time step
374         rpktrd = peke
375         !
376      ENDIF
377      !
378      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
379      !
380   END SUBROUTINE glo_dyn_wri
381
382
383   SUBROUTINE glo_tra_wri( kt )
384      !!---------------------------------------------------------------------
385      !!                  ***  ROUTINE glo_tra_wri  ***
386      !!
387      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
388      !!----------------------------------------------------------------------
389      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
390      !
391      INTEGER  ::   jk   ! loop indices
392      !!----------------------------------------------------------------------
393
394      ! I. Tracers trends
395      ! -----------------
396
397      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
398
399         ! I.1 Sums over the global domain
400         ! -------------------------------
401         IF( lk_mpp ) THEN
402            CALL mpp_sum( tmo, jptot_tra )   
403            CALL mpp_sum( smo, jptot_tra )
404            CALL mpp_sum( t2 , jptot_tra )
405            CALL mpp_sum( s2 , jptot_tra )
406         ENDIF
407
408         ! I.2 Print tracers trends in the ocean.output file
409         ! --------------------------------------------------
410         
411         IF(lwp) THEN
412            WRITE (numout,*)
413            WRITE (numout,*)
414            WRITE (numout,9400) kt
415            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
416            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
417            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
418            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
419            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
420            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
421            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
422            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
423            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
424            WRITE (numout,9409) 
425            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
426            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
427            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
428            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
429            IF(lflush) CALL flush(numout)
430         ENDIF
431
4329400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
433              '              salinity',/' ============================')
4349401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4359411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4369402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4379403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4389404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4399405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4409406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4419407     FORMAT(' penetrative qsr             ',e20.13)
4429408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4439409     FORMAT(' -------------------------------------------------------------------------')
4449410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
445
446
447         IF(lwp) THEN
448            WRITE (numout,*)
449            WRITE (numout,*)
450            WRITE (numout,9420) kt
451            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
452            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
453            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
454            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
455            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
456            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
457            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
458            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
459            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
460            WRITE (numout,9429)
461            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
462            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
463            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
464            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
465            IF(lflush) CALL flush(numout)
466         ENDIF
467
4689420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
469            '               salinity', /, ' ===============================')
4709421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4719431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4729422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4739423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4749424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4759425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4769426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4779427     FORMAT(' penetrative qsr           * t   ', e20.13)
4789428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4799429     FORMAT(' -----------------------------------------------------------------------------')
4809430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
481
482
483         IF(lwp) THEN
484            WRITE (numout,*)
485            WRITE (numout,*)
486            WRITE (numout,9440) kt
487            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
488            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
489            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
490            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
491            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
492            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
493            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
494            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
495            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
496            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
497            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
498            IF(lflush) CALL flush(numout)
499         ENDIF
500
5019440     FORMAT(' tracer consistency at it= ',i6,   &
502            ' :         temperature','                salinity',/,   &
503            ' ==================================')
5049441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5059442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5069443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5079444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5089445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5099446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5109447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5119448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5129449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
513         !
514      ENDIF
515      !
516   END SUBROUTINE glo_tra_wri
517
518
519   SUBROUTINE trd_glo_init
520      !!---------------------------------------------------------------------
521      !!                  ***  ROUTINE trd_glo_init  ***
522      !!
523      !! ** Purpose :   Read the namtrd namelist
524      !!----------------------------------------------------------------------
525      INTEGER  ::   ji, jj, jk   ! dummy loop indices
526      !!----------------------------------------------------------------------
527
528      IF(lwp) THEN
529         WRITE(numout,*)
530         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
531         WRITE(numout,*) '~~~~~~~~~~~~~'
532         IF(lflush) CALL flush(numout)
533      ENDIF
534
535      ! Total volume at t-points:
536      tvolt = 0._wp
537      DO jk = 1, jpkm1
538         tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
539      END DO
540      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
541
542      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
543      IF(lwp .AND. lflush) CALL flush(numout)
544
545      ! Initialization of potential to kinetic energy conversion
546      rpktrd = 0._wp
547
548      ! Total volume at u-, v- points:
549!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
550      tvolu = 0._wp
551      tvolv = 0._wp
552
553      DO jk = 1, jpk
554         DO jj = 2, jpjm1
555            DO ji = fs_2, fs_jpim1   ! vector opt.
556               tvolu = tvolu + e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
557               tvolv = tvolv + e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
558            END DO
559         END DO
560      END DO
561      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
562      IF( lk_mpp )   CALL mpp_sum( tvolv )
563
564      IF(lwp) THEN
565         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
566         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
567         IF(lflush) CALL flush(numout)
568      ENDIF
569      !
570   END SUBROUTINE trd_glo_init
571
572   !!======================================================================
573END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.