source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 18 months ago

The Dr Hook changes from my perl code.

File size: 30.2 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   USE yomhook, ONLY: lhook, dr_hook
34   USE parkind1, ONLY: jprb, jpim
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   trd_glo       ! called by trdtra and trddyn modules
40   PUBLIC   trd_glo_init  ! called by trdini module
41
42   !                     !!! Variables used for diagnostics
43   REAL(wp) ::   tvolt    ! volume of the whole ocean computed at t-points
44   REAL(wp) ::   tvolu    ! volume of the whole ocean computed at u-points
45   REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
46   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
47   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
48
49   !                     !!! domain averaged trends
50   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
51   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
52   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
53   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "vectopt_loop_substitute.h90"
58#  include "zdfddm_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
67      !!---------------------------------------------------------------------
68      !!                  ***  ROUTINE trd_glo  ***
69      !!
70      !! ** Purpose :   compute and print global domain averaged trends for
71      !!              T, T^2, momentum, KE, and KE<->PE
72      !!
73      !!----------------------------------------------------------------------
74      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
75      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
76      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
77      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
78      INTEGER                   , INTENT(in   ) ::   kt      ! time step
79      !!
80      INTEGER ::   ji, jj, jk      ! dummy loop indices
81      INTEGER ::   ikbu, ikbv      ! local integers
82      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
83      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
84      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
85      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
86      REAL(KIND=jprb)               :: zhook_handle
87
88      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRD_GLO'
89
90      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
91
92      !!----------------------------------------------------------------------
93
94      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
95
96      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
97         !
98         SELECT CASE( ctype )
99         !
100         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
101            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
102               DO jj = 1, jpj
103                  DO ji = 1, jpi       
104                     zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
105                     zvt = ptrdx(ji,jj,jk) * zvm
106                     zvs = ptrdy(ji,jj,jk) * zvm
107                     tmo(ktrd) = tmo(ktrd) + zvt   
108                     smo(ktrd) = smo(ktrd) + zvs
109                     t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
110                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
111                  END DO
112               END DO
113            END DO
114            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
115            IF( .NOT.lk_vvl .AND. ktrd == jptra_zad ) THEN 
116               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
117               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
118               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
119               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
120            ENDIF
121            !
122            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
123               !
124               CALL glo_tra_wri( kt )             ! print the results in ocean.output
125               !               
126               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
127               smo(:) = 0._wp
128               t2 (:) = 0._wp
129               s2 (:) = 0._wp
130               !
131            ENDIF
132            !
133         CASE( 'DYN' )          !==  Momentum and KE  ==!       
134            DO jk = 1, jpkm1
135               DO jj = 1, jpjm1
136                  DO ji = 1, jpim1
137                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
138                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
139                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
140                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
141                     umo(ktrd) = umo(ktrd) + zvt
142                     vmo(ktrd) = vmo(ktrd) + zvs
143                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
144                  END DO
145               END DO
146            END DO
147            !                 
148            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
149               z1_2rau0 = 0.5_wp / rau0
150               DO jj = 1, jpjm1
151                  DO ji = 1, jpim1
152                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
153                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
154                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
155                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
156                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
157                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
158                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
159                  END DO
160               END DO
161            ENDIF
162            !                       
163            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
164               !
165               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction
166                  z1_2rau0 = 0.5_wp / rau0
167                  DO jj = 1, jpjm1
168                     DO ji = 1, jpim1
169                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
170                        ikbv = mbkv(ji,jj)
171                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
172                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
173                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
174                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
175                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
176                     END DO
177                  END DO
178               ENDIF
179               !
180               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
181               !               
182               umo(:) = 0._wp                         ! reset for the next time step
183               vmo(:) = 0._wp
184               hke(:) = 0._wp
185               !
186            ENDIF
187            !
188         END SELECT
189         !
190      ENDIF
191      !
192      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
193      !
194      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
195   END SUBROUTINE trd_glo
196
197
198   SUBROUTINE glo_dyn_wri( kt )
199      !!---------------------------------------------------------------------
200      !!                  ***  ROUTINE glo_dyn_wri  ***
201      !!
202      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
203      !!----------------------------------------------------------------------
204      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
205      !
206      INTEGER  ::   ji, jj, jk   ! dummy loop indices
207      REAL(wp) ::   zcof         ! local scalar
208      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
209      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
210      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
211      REAL(KIND=jprb)               :: zhook_handle
212
213      CHARACTER(LEN=*), PARAMETER :: RoutineName='GLO_DYN_WRI'
214
215      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
216
217      !!----------------------------------------------------------------------
218
219      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
220
221      ! I. Momentum trends
222      ! -------------------
223
224      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
225
226         ! I.1 Conversion potential energy - kinetic energy
227         ! --------------------------------------------------
228         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
229         zkx  (:,:,:) = 0._wp
230         zky  (:,:,:) = 0._wp
231         zkz  (:,:,:) = 0._wp
232         zkepe(:,:,:) = 0._wp
233   
234         CALL eos( tsn, rhd, rhop )       ! now potential density
235
236         zcof = 0.5_wp / rau0             ! Density flux at w-point
237         zkz(:,:,1) = 0._wp
238         DO jk = 2, jpk
239            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
240         END DO
241         
242         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
243         DO jk = 1, jpkm1
244            DO jj = 1, jpjm1
245               DO ji = 1, jpim1
246                  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) )
247                  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) )
248               END DO
249            END DO
250         END DO
251         
252         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
253            DO jj = 2, jpjm1
254               DO ji = 2, jpim1
255                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
256                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
257                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
258                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
259               END DO
260            END DO
261         END DO
262
263         ! I.2 Basin averaged kinetic energy trend
264         ! ----------------------------------------
265         peke = 0._wp
266         DO jk = 1, jpkm1
267            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
268         END DO
269         peke = grav * peke
270
271         ! I.3 Sums over the global domain
272         ! ---------------------------------
273         IF( lk_mpp ) THEN
274            CALL mpp_sum( peke )
275            CALL mpp_sum( umo , jptot_dyn )
276            CALL mpp_sum( vmo , jptot_dyn )
277            CALL mpp_sum( hke , jptot_dyn )
278         ENDIF
279
280         ! I.2 Print dynamic trends in the ocean.output file
281         ! --------------------------------------------------
282
283         IF(lwp) THEN
284            WRITE (numout,*)
285            WRITE (numout,*)
286            WRITE (numout,9500) kt
287            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
288            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
289            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
290            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
291            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
292            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
293            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
294            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
295            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
296            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
297            WRITE (numout,9511)
298            WRITE (numout,9512)                                                 &
299            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
300            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
301            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
302            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
303            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
304            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
305            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
306            IF( ln_bfrimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
307         ENDIF
308
309 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
310 9501    FORMAT(' hydro pressure gradient    u= ', e20.13, '    v= ', e20.13)
311 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
312 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
313 9504    FORMAT(' planetary vorticity term   u= ', e20.13, '    v= ', e20.13)
314 9505    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
315 9506    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
316 9507    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
317 9508    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
318 9509    FORMAT(' explicit bottom friction   u= ', e20.13, '    v= ', e20.13)
319 9510    FORMAT(' Asselin time filter        u= ', e20.13, '    v= ', e20.13)
320 9511    FORMAT(' -----------------------------------------------------------------------------')
321 9512    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
322 9513    FORMAT(' incl. surface wind stress  u= ', e20.13, '    v= ', e20.13)
323 9514    FORMAT('       bottom stress        u= ', e20.13, '    v= ', e20.13)
324
325         IF(lwp) THEN
326            WRITE (numout,*)
327            WRITE (numout,*)
328            WRITE (numout,9520) kt
329            WRITE (numout,9521) hke(jpdyn_hpg) / tvolt
330            WRITE (numout,9522) hke(jpdyn_keg) / tvolt
331            WRITE (numout,9523) hke(jpdyn_rvo) / tvolt
332            WRITE (numout,9524) hke(jpdyn_pvo) / tvolt
333            WRITE (numout,9525) hke(jpdyn_zad) / tvolt
334            WRITE (numout,9526) hke(jpdyn_ldf) / tvolt
335            WRITE (numout,9527) hke(jpdyn_zdf) / tvolt
336            WRITE (numout,9528) hke(jpdyn_spg) / tvolt
337            WRITE (numout,9529) hke(jpdyn_bfr) / tvolt
338            WRITE (numout,9530) hke(jpdyn_atf) / tvolt
339            WRITE (numout,9531)
340            WRITE (numout,9532)   &
341            &     (  hke(jpdyn_hpg) + hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_pvo)   &
342            &      + hke(jpdyn_zad) + hke(jpdyn_ldf) + hke(jpdyn_zdf) + hke(jpdyn_spg)   &
343            &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
344            WRITE (numout,9533) hke(jpdyn_tau) / tvolt
345            IF( ln_bfrimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
346         ENDIF
347
348 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
349 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
350 9522    FORMAT(' ke gradient               u2= ', e20.13)
351 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
352 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
353 9525    FORMAT(' vertical advection        u2= ', e20.13)
354 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
355 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
356 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
357 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
358 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
359 9531    FORMAT(' --------------------------------------------------')
360 9532    FORMAT(' total trend               u2= ', e20.13)
361 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
362 9534    FORMAT('       bottom stress       u2= ', e20.13)
363
364         IF(lwp) THEN
365            WRITE (numout,*)
366            WRITE (numout,*)
367            WRITE (numout,9540) kt
368            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
369            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
370            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
371            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
372            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
373            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
374            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
375            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
376            WRITE (numout,*)
377            WRITE (numout,*)
378         ENDIF
379
380 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
381 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
382 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
383 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
384 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
385 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
386 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
387 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
388 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
389         !
390         ! Save potential to kinetic energy conversion for next time step
391         rpktrd = peke
392         !
393      ENDIF
394      !
395      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
396      !
397      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
398   END SUBROUTINE glo_dyn_wri
399
400
401   SUBROUTINE glo_tra_wri( kt )
402      !!---------------------------------------------------------------------
403      !!                  ***  ROUTINE glo_tra_wri  ***
404      !!
405      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
406      !!----------------------------------------------------------------------
407      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
408      !
409      INTEGER  ::   jk   ! loop indices
410      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
411      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
412      REAL(KIND=jprb)               :: zhook_handle
413
414      CHARACTER(LEN=*), PARAMETER :: RoutineName='GLO_TRA_WRI'
415
416      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
417
418      !!----------------------------------------------------------------------
419
420      ! I. Tracers trends
421      ! -----------------
422
423      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
424
425         ! I.1 Sums over the global domain
426         ! -------------------------------
427         IF( lk_mpp ) THEN
428            CALL mpp_sum( tmo, jptot_tra )   
429            CALL mpp_sum( smo, jptot_tra )
430            CALL mpp_sum( t2 , jptot_tra )
431            CALL mpp_sum( s2 , jptot_tra )
432         ENDIF
433
434         ! I.2 Print tracers trends in the ocean.output file
435         ! --------------------------------------------------
436         
437         IF(lwp) THEN
438            WRITE (numout,*)
439            WRITE (numout,*)
440            WRITE (numout,9400) kt
441            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
442            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
443            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
444            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
445            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
446            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
447            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
448            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
449            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
450            WRITE (numout,9409) 
451            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
452            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
453            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
454            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
455         ENDIF
456
4579400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
458              '              salinity',/' ============================')
4599401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4609411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4619402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4629403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4639404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4649405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4659406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4669407     FORMAT(' penetrative qsr             ',e20.13)
4679408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4689409     FORMAT(' -------------------------------------------------------------------------')
4699410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
470
471
472         IF(lwp) THEN
473            WRITE (numout,*)
474            WRITE (numout,*)
475            WRITE (numout,9420) kt
476            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
477            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
478            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
479            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
480            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
481            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
482            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
483            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
484            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
485            WRITE (numout,9429)
486            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
487            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
488            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
489            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
490         ENDIF
491
4929420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
493            '               salinity', /, ' ===============================')
4949421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4959431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4969422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4979423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4989424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4999425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
5009426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
5019427     FORMAT(' penetrative qsr           * t   ', e20.13)
5029428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
5039429     FORMAT(' -----------------------------------------------------------------------------')
5049430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
505
506
507         IF(lwp) THEN
508            WRITE (numout,*)
509            WRITE (numout,*)
510            WRITE (numout,9440) kt
511            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
512            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
513            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
514            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
515            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
516            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
517            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
518            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
519            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
520            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
521            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
522         ENDIF
523
5249440     FORMAT(' tracer consistency at it= ',i6,   &
525            ' :         temperature','                salinity',/,   &
526            ' ==================================')
5279441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5289442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5299443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5309444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5319445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5329446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5339447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5349448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5359449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
536         !
537      ENDIF
538      !
539      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
540   END SUBROUTINE glo_tra_wri
541
542
543   SUBROUTINE trd_glo_init
544      !!---------------------------------------------------------------------
545      !!                  ***  ROUTINE trd_glo_init  ***
546      !!
547      !! ** Purpose :   Read the namtrd namelist
548      !!----------------------------------------------------------------------
549      INTEGER  ::   ji, jj, jk   ! dummy loop indices
550      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
551      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
552      REAL(KIND=jprb)               :: zhook_handle
553
554      CHARACTER(LEN=*), PARAMETER :: RoutineName='TRD_GLO_INIT'
555
556      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
557
558      !!----------------------------------------------------------------------
559
560      IF(lwp) THEN
561         WRITE(numout,*)
562         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
563         WRITE(numout,*) '~~~~~~~~~~~~~'
564      ENDIF
565
566      ! Total volume at t-points:
567      tvolt = 0._wp
568      DO jk = 1, jpkm1
569         tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
570      END DO
571      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
572
573      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
574
575      ! Initialization of potential to kinetic energy conversion
576      rpktrd = 0._wp
577
578      ! Total volume at u-, v- points:
579!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
580      tvolu = 0._wp
581      tvolv = 0._wp
582
583      DO jk = 1, jpk
584         DO jj = 2, jpjm1
585            DO ji = fs_2, fs_jpim1   ! vector opt.
586               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)
587               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)
588            END DO
589         END DO
590      END DO
591      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
592      IF( lk_mpp )   CALL mpp_sum( tvolv )
593
594      IF(lwp) THEN
595         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
596         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
597      ENDIF
598      !
599      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
600   END SUBROUTINE trd_glo_init
601
602   !!======================================================================
603END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.