source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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