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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • 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          ! lateral diffusion: eddy diffusivity & EIV coeff.
22   USE ldfdyn          ! 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   !
29   USE lib_mpp         ! distibuted memory computing library
30   USE in_out_manager  ! I/O manager
31   USE iom             ! I/O manager library
32   USE wrk_nemo        ! Memory allocation
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   trd_glo       ! called by trdtra and trddyn modules
38   PUBLIC   trd_glo_init  ! called by trdini module
39
40   !                     !!! Variables used for diagnostics
41   REAL(wp) ::   tvolt    ! volume of the whole ocean computed at t-points
42   REAL(wp) ::   tvolu    ! volume of the whole ocean computed at u-points
43   REAL(wp) ::   tvolv    ! volume of the whole ocean computed at v-points
44   REAL(wp) ::   rpktrd   ! potential to kinetic energy conversion
45   REAL(wp) ::   peke     ! conversion potential energy - kinetic energy trend
46
47   !                     !!! domain averaged trends
48   REAL(wp), DIMENSION(jptot_tra) ::   tmo, smo   ! temperature and salinity trends
49   REAL(wp), DIMENSION(jptot_tra) ::   t2 , s2    ! T^2 and S^2 trends
50   REAL(wp), DIMENSION(jptot_dyn) ::   umo, vmo   ! momentum trends
51   REAL(wp), DIMENSION(jptot_dyn) ::   hke        ! kinetic energy trends (u^2+v^2)
52
53   !! * Substitutions
54#  include "domzgr_substitute.h90"
55#  include "vectopt_loop_substitute.h90"
56#  include "zdfddm_substitute.h90"
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE trd_glo( ptrdx, ptrdy, ktrd, ctype, kt )
65      !!---------------------------------------------------------------------
66      !!                  ***  ROUTINE trd_glo  ***
67      !!
68      !! ** Purpose :   compute and print global domain averaged trends for
69      !!              T, T^2, momentum, KE, and KE<->PE
70      !!
71      !!----------------------------------------------------------------------
72      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
73      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
74      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
75      CHARACTER(len=3)          , INTENT(in   ) ::   ctype   ! momentum or tracers trends type (='DYN'/'TRA')
76      INTEGER                   , INTENT(in   ) ::   kt      ! time step
77      !!
78      INTEGER ::   ji, jj, jk      ! dummy loop indices
79      INTEGER ::   ikbu, ikbv      ! local integers
80      REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
81      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace
82      !!----------------------------------------------------------------------
83
84      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
85
86      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
87         !
88         SELECT CASE( ctype )
89         !
90         CASE( 'TRA' )          !==  Tracers (T & S)  ==!
91            DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask)
92               DO jj = 1, jpj
93                  DO ji = 1, jpi       
94                     zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
95                     zvt = ptrdx(ji,jj,jk) * zvm
96                     zvs = ptrdy(ji,jj,jk) * zvm
97                     tmo(ktrd) = tmo(ktrd) + zvt   
98                     smo(ktrd) = smo(ktrd) + zvs
99                     t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem)
100                     s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal)
101                  END DO
102               END DO
103            END DO
104            !                       ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
105            IF( .NOT.lk_vvl .AND. ktrd == jptra_zad ) THEN 
106               tmo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * e1e2t(:,:) )
107               smo(jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
108               t2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) * e1e2t(:,:)  )
109               s2 (jptra_sad) = SUM( wn(:,:,1) * tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) * e1e2t(:,:)  )
110            ENDIF
111            !
112            IF( ktrd == jptra_atf ) THEN     ! last trend (asselin time filter)
113               !
114               CALL glo_tra_wri( kt )             ! print the results in ocean.output
115               !               
116               tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero)
117               smo(:) = 0._wp
118               t2 (:) = 0._wp
119               s2 (:) = 0._wp
120               !
121            ENDIF
122            !
123         CASE( 'DYN' )          !==  Momentum and KE  ==!       
124            DO jk = 1, jpkm1
125               DO jj = 1, jpjm1
126                  DO ji = 1, jpim1
127                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
128                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * fse3u(ji,jj,jk)
129                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
130                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
131                     umo(ktrd) = umo(ktrd) + zvt
132                     vmo(ktrd) = vmo(ktrd) + zvs
133                     hke(ktrd) = hke(ktrd) + un(ji,jj,jk) * zvt + vn(ji,jj,jk) * zvs
134                  END DO
135               END DO
136            END DO
137            !                 
138            IF( ktrd == jpdyn_zdf ) THEN      ! zdf trend: compute separately the surface forcing trend
139               z1_2rau0 = 0.5_wp / rau0
140               DO jj = 1, jpjm1
141                  DO ji = 1, jpim1
142                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
143                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
144                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
145                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * fse3u(ji,jj,jk)
146                     umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
147                     vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
148                     hke(jpdyn_tau) = hke(jpdyn_tau) + un(ji,jj,1) * zvt + vn(ji,jj,1) * zvs
149                  END DO
150               END DO
151            ENDIF
152            !                       
153            IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
154               !
155               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction
156                  z1_2rau0 = 0.5_wp / rau0
157                  DO jj = 1, jpjm1
158                     DO ji = 1, jpim1
159                        ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
160                        ikbv = mbkv(ji,jj)
161                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
162                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
163                        umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
164                        vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
165                        hke(jpdyn_bfri) = hke(jpdyn_bfri) + un(ji,jj,ikbu) * zvt + vn(ji,jj,ikbv) * zvs
166                     END DO
167                  END DO
168               ENDIF
169               !
170               CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
171               !               
172               umo(:) = 0._wp                         ! reset for the next time step
173               vmo(:) = 0._wp
174               hke(:) = 0._wp
175               !
176            ENDIF
177            !
178         END SELECT
179         !
180      ENDIF
181      !
182      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
183      !
184   END SUBROUTINE trd_glo
185
186
187   SUBROUTINE glo_dyn_wri( kt )
188      !!---------------------------------------------------------------------
189      !!                  ***  ROUTINE glo_dyn_wri  ***
190      !!
191      !! ** Purpose :  write global averaged U, KE, PE<->KE trends in ocean.output
192      !!----------------------------------------------------------------------
193      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
194      !
195      INTEGER  ::   ji, jj, jk   ! dummy loop indices
196      REAL(wp) ::   zcof         ! local scalar
197      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe 
198      !!----------------------------------------------------------------------
199
200      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
201
202      ! I. Momentum trends
203      ! -------------------
204
205      IF( MOD( kt, nn_trd ) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
206
207         ! I.1 Conversion potential energy - kinetic energy
208         ! --------------------------------------------------
209         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
210         zkx  (:,:,:) = 0._wp
211         zky  (:,:,:) = 0._wp
212         zkz  (:,:,:) = 0._wp
213         zkepe(:,:,:) = 0._wp
214   
215         CALL eos( tsn, rhd, rhop )       ! now potential density
216
217         zcof = 0.5_wp / rau0             ! Density flux at w-point
218         zkz(:,:,1) = 0._wp
219         DO jk = 2, jpk
220            zkz(:,:,jk) = zcof * e1e2t(:,:) * wn(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
221         END DO
222         
223         zcof   = 0.5_wp / rau0           ! Density flux at u and v-points
224         DO jk = 1, jpkm1
225            DO jj = 1, jpjm1
226               DO ji = 1, jpim1
227                  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) )
228                  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) )
229               END DO
230            END DO
231         END DO
232         
233         DO jk = 1, jpkm1                 ! Density flux divergence at t-point
234            DO jj = 2, jpjm1
235               DO ji = 2, jpim1
236                  zkepe(ji,jj,jk) = - (  zkz(ji,jj,jk) - zkz(ji  ,jj  ,jk+1)               &
237                     &                 + zkx(ji,jj,jk) - zkx(ji-1,jj  ,jk  )               &
238                     &                 + zky(ji,jj,jk) - zky(ji  ,jj-1,jk  )   )           &
239                     &              / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) * tmask_i(ji,jj)
240               END DO
241            END DO
242         END DO
243
244         ! I.2 Basin averaged kinetic energy trend
245         ! ----------------------------------------
246         peke = 0._wp
247         DO jk = 1, jpkm1
248            peke = peke + SUM( zkepe(:,:,jk) * fsdept(:,:,jk) * e1e2t(:,:) * fse3t(:,:,jk) )
249         END DO
250         peke = grav * peke
251
252         ! I.3 Sums over the global domain
253         ! ---------------------------------
254         IF( lk_mpp ) THEN
255            CALL mpp_sum( peke )
256            CALL mpp_sum( umo , jptot_dyn )
257            CALL mpp_sum( vmo , jptot_dyn )
258            CALL mpp_sum( hke , jptot_dyn )
259         ENDIF
260
261         ! I.2 Print dynamic trends in the ocean.output file
262         ! --------------------------------------------------
263
264         IF(lwp) THEN
265            WRITE (numout,*)
266            WRITE (numout,*)
267            WRITE (numout,9500) kt
268            WRITE (numout,9501) umo(jpdyn_hpg) / tvolu, vmo(jpdyn_hpg) / tvolv
269            WRITE (numout,9502) umo(jpdyn_keg) / tvolu, vmo(jpdyn_keg) / tvolv
270            WRITE (numout,9503) umo(jpdyn_rvo) / tvolu, vmo(jpdyn_rvo) / tvolv
271            WRITE (numout,9504) umo(jpdyn_pvo) / tvolu, vmo(jpdyn_pvo) / tvolv
272            WRITE (numout,9505) umo(jpdyn_zad) / tvolu, vmo(jpdyn_zad) / tvolv
273            WRITE (numout,9506) umo(jpdyn_ldf) / tvolu, vmo(jpdyn_ldf) / tvolv
274            WRITE (numout,9507) umo(jpdyn_zdf) / tvolu, vmo(jpdyn_zdf) / tvolv
275            WRITE (numout,9508) umo(jpdyn_spg) / tvolu, vmo(jpdyn_spg) / tvolv
276            WRITE (numout,9509) umo(jpdyn_bfr) / tvolu, vmo(jpdyn_bfr) / tvolv
277            WRITE (numout,9510) umo(jpdyn_atf) / tvolu, vmo(jpdyn_atf) / tvolv
278            WRITE (numout,9511)
279            WRITE (numout,9512)                                                 &
280            &     (  umo(jpdyn_hpg) + umo(jpdyn_keg) + umo(jpdyn_rvo) + umo(jpdyn_pvo)   &
281            &      + umo(jpdyn_zad) + umo(jpdyn_ldf) + umo(jpdyn_zdf) + umo(jpdyn_spg)   &
282            &      + umo(jpdyn_bfr) + umo(jpdyn_atf) ) / tvolu,   &
283            &     (  vmo(jpdyn_hpg) + vmo(jpdyn_keg) + vmo(jpdyn_rvo) + vmo(jpdyn_pvo)   &
284            &      + vmo(jpdyn_zad) + vmo(jpdyn_ldf) + vmo(jpdyn_zdf) + vmo(jpdyn_spg)   &
285            &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
286            WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
287            IF( ln_bfrimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
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         ENDIF
328
329 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
330 9521    FORMAT(' hydro pressure gradient   u2= ', e20.13)
331 9522    FORMAT(' ke gradient               u2= ', e20.13)
332 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
333 9524    FORMAT(' planetary vorticity term  u2= ', e20.13)
334 9525    FORMAT(' vertical advection        u2= ', e20.13)
335 9526    FORMAT(' horizontal diffusion      u2= ', e20.13)
336 9527    FORMAT(' vertical diffusion        u2= ', e20.13)
337 9528    FORMAT(' surface pressure gradient u2= ', e20.13)
338 9529    FORMAT(' explicit bottom friction  u2= ', e20.13)
339 9530    FORMAT(' Asselin time filter       u2= ', e20.13)
340 9531    FORMAT(' --------------------------------------------------')
341 9532    FORMAT(' total trend               u2= ', e20.13)
342 9533    FORMAT(' incl. surface wind stress u2= ', e20.13)
343 9534    FORMAT('       bottom stress       u2= ', e20.13)
344
345         IF(lwp) THEN
346            WRITE (numout,*)
347            WRITE (numout,*)
348            WRITE (numout,9540) kt
349            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt
350            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt
351            WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt
352            WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt
353            WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt
354            WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt
355            WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt
356            WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt
357            WRITE (numout,*)
358            WRITE (numout,*)
359         ENDIF
360
361 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
362 9541    FORMAT(' 0 = non linear term (true if KE conserved)                : ', e20.13)
363 9542    FORMAT(' 0 = ke gradient + vertical advection                      : ', e20.13)
364 9543    FORMAT(' 0 = coriolis term  (true if KE conserving scheme)         : ', e20.13)
365 9544    FORMAT(' 0 = vorticity term (true if KE conserving scheme)         : ', e20.13)
366 9545    FORMAT(' 0 = surface pressure gradient  ???                        : ', e20.13)
367 9546    FORMAT(' 0 < horizontal diffusion                                  : ', e20.13)
368 9547    FORMAT(' 0 < vertical diffusion                                    : ', e20.13)
369 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
370         !
371         ! Save potential to kinetic energy conversion for next time step
372         rpktrd = peke
373         !
374      ENDIF
375      !
376      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
377      !
378   END SUBROUTINE glo_dyn_wri
379
380
381   SUBROUTINE glo_tra_wri( kt )
382      !!---------------------------------------------------------------------
383      !!                  ***  ROUTINE glo_tra_wri  ***
384      !!
385      !! ** Purpose :  write global domain averaged of T and T^2 trends in ocean.output
386      !!----------------------------------------------------------------------
387      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
388      !
389      INTEGER  ::   jk   ! loop indices
390      !!----------------------------------------------------------------------
391
392      ! I. Tracers trends
393      ! -----------------
394
395      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
396
397         ! I.1 Sums over the global domain
398         ! -------------------------------
399         IF( lk_mpp ) THEN
400            CALL mpp_sum( tmo, jptot_tra )   
401            CALL mpp_sum( smo, jptot_tra )
402            CALL mpp_sum( t2 , jptot_tra )
403            CALL mpp_sum( s2 , jptot_tra )
404         ENDIF
405
406         ! I.2 Print tracers trends in the ocean.output file
407         ! --------------------------------------------------
408         
409         IF(lwp) THEN
410            WRITE (numout,*)
411            WRITE (numout,*)
412            WRITE (numout,9400) kt
413            WRITE (numout,9401)  tmo(jptra_xad) / tvolt, smo(jptra_xad) / tvolt
414            WRITE (numout,9411)  tmo(jptra_yad) / tvolt, smo(jptra_yad) / tvolt
415            WRITE (numout,9402)  tmo(jptra_zad) / tvolt, smo(jptra_zad) / tvolt
416            WRITE (numout,9403)  tmo(jptra_ldf) / tvolt, smo(jptra_ldf) / tvolt
417            WRITE (numout,9404)  tmo(jptra_zdf) / tvolt, smo(jptra_zdf) / tvolt
418            WRITE (numout,9405)  tmo(jptra_npc) / tvolt, smo(jptra_npc) / tvolt
419            WRITE (numout,9406)  tmo(jptra_dmp) / tvolt, smo(jptra_dmp) / tvolt
420            WRITE (numout,9407)  tmo(jptra_qsr) / tvolt
421            WRITE (numout,9408)  tmo(jptra_nsr) / tvolt, smo(jptra_nsr) / tvolt
422            WRITE (numout,9409) 
423            WRITE (numout,9410) (  tmo(jptra_xad) + tmo(jptra_yad) + tmo(jptra_zad) + tmo(jptra_ldf) + tmo(jptra_zdf)   &
424            &                    + tmo(jptra_npc) + tmo(jptra_dmp) + tmo(jptra_qsr) + tmo(jptra_nsr) ) / tvolt,   &
425            &                   (  smo(jptra_xad) + smo(jptra_yad) + smo(jptra_zad) + smo(jptra_ldf) + smo(jptra_zdf)   &
426            &                    + smo(jptra_npc) + smo(jptra_dmp)                   + smo(jptra_nsr) ) / tvolt
427         ENDIF
428
4299400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
430              '              salinity',/' ============================')
4319401     FORMAT(' zonal      advection        ',e20.13,'     ',e20.13)
4329411     FORMAT(' meridional advection        ',e20.13,'     ',e20.13)
4339402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
4349403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
4359404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
4369405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
4379406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
4389407     FORMAT(' penetrative qsr             ',e20.13)
4399408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
4409409     FORMAT(' -------------------------------------------------------------------------')
4419410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
442
443
444         IF(lwp) THEN
445            WRITE (numout,*)
446            WRITE (numout,*)
447            WRITE (numout,9420) kt
448            WRITE (numout,9421)   t2(jptra_xad) / tvolt, s2(jptra_xad) / tvolt
449            WRITE (numout,9431)   t2(jptra_yad) / tvolt, s2(jptra_yad) / tvolt
450            WRITE (numout,9422)   t2(jptra_zad) / tvolt, s2(jptra_zad) / tvolt
451            WRITE (numout,9423)   t2(jptra_ldf) / tvolt, s2(jptra_ldf) / tvolt
452            WRITE (numout,9424)   t2(jptra_zdf) / tvolt, s2(jptra_zdf) / tvolt
453            WRITE (numout,9425)   t2(jptra_npc) / tvolt, s2(jptra_npc) / tvolt
454            WRITE (numout,9426)   t2(jptra_dmp) / tvolt, s2(jptra_dmp) / tvolt
455            WRITE (numout,9427)   t2(jptra_qsr) / tvolt
456            WRITE (numout,9428)   t2(jptra_nsr) / tvolt, s2(jptra_nsr) / tvolt
457            WRITE (numout,9429)
458            WRITE (numout,9430) (  t2(jptra_xad) + t2(jptra_yad) + t2(jptra_zad) + t2(jptra_ldf) + t2(jptra_zdf)   &
459            &                    + t2(jptra_npc) + t2(jptra_dmp) + t2(jptra_qsr) + t2(jptra_nsr) ) / tvolt,   &
460            &                   (  s2(jptra_xad) + s2(jptra_yad) + s2(jptra_zad) + s2(jptra_ldf) + s2(jptra_zdf)   &
461            &                    + s2(jptra_npc) + s2(jptra_dmp)                  + s2(jptra_nsr) ) / tvolt
462         ENDIF
463
4649420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
465            '               salinity', /, ' ===============================')
4669421     FORMAT(' zonal      advection      * t   ', e20.13, '     ', e20.13)
4679431     FORMAT(' meridional advection      * t   ', e20.13, '     ', e20.13)
4689422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
4699423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
4709424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
4719425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
4729426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
4739427     FORMAT(' penetrative qsr           * t   ', e20.13)
4749428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
4759429     FORMAT(' -----------------------------------------------------------------------------')
4769430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
477
478
479         IF(lwp) THEN
480            WRITE (numout,*)
481            WRITE (numout,*)
482            WRITE (numout,9440) kt
483            WRITE (numout,9441) ( tmo(jptra_xad)+tmo(jptra_yad)+tmo(jptra_zad) )/tvolt,   &
484            &                   ( smo(jptra_xad)+smo(jptra_yad)+smo(jptra_zad) )/tvolt
485            WRITE (numout,9442)   tmo(jptra_sad)/tvolt,  smo(jptra_sad)/tvolt
486            WRITE (numout,9443)   tmo(jptra_ldf)/tvolt,  smo(jptra_ldf)/tvolt
487            WRITE (numout,9444)   tmo(jptra_zdf)/tvolt,  smo(jptra_zdf)/tvolt
488            WRITE (numout,9445)   tmo(jptra_npc)/tvolt,  smo(jptra_npc)/tvolt
489            WRITE (numout,9446) ( t2(jptra_xad)+t2(jptra_yad)+t2(jptra_zad) )/tvolt,    &
490            &                   ( s2(jptra_xad)+s2(jptra_yad)+s2(jptra_zad) )/tvolt
491            WRITE (numout,9447)   t2(jptra_ldf)/tvolt,   s2(jptra_ldf)/tvolt
492            WRITE (numout,9448)   t2(jptra_zdf)/tvolt,   s2(jptra_zdf)/tvolt
493            WRITE (numout,9449)   t2(jptra_npc)/tvolt,   s2(jptra_npc)/tvolt
494         ENDIF
495
4969440     FORMAT(' tracer consistency at it= ',i6,   &
497            ' :         temperature','                salinity',/,   &
498            ' ==================================')
4999441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
5009442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
5019443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
5029444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
5039445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
5049446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
5059447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
5069448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
5079449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
508         !
509      ENDIF
510      !
511   END SUBROUTINE glo_tra_wri
512
513
514   SUBROUTINE trd_glo_init
515      !!---------------------------------------------------------------------
516      !!                  ***  ROUTINE trd_glo_init  ***
517      !!
518      !! ** Purpose :   Read the namtrd namelist
519      !!----------------------------------------------------------------------
520      INTEGER  ::   ji, jj, jk   ! dummy loop indices
521      !!----------------------------------------------------------------------
522
523      IF(lwp) THEN
524         WRITE(numout,*)
525         WRITE(numout,*) 'trd_glo_init : integral constraints properties trends'
526         WRITE(numout,*) '~~~~~~~~~~~~~'
527      ENDIF
528
529      ! Total volume at t-points:
530      tvolt = 0._wp
531      DO jk = 1, jpkm1
532         tvolt = tvolt + SUM( e1e2t(:,:) * fse3t(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) )
533      END DO
534      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
535
536      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
537
538      ! Initialization of potential to kinetic energy conversion
539      rpktrd = 0._wp
540
541      ! Total volume at u-, v- points:
542!!gm :  bug?  je suis quasi sur que le produit des tmask_i ne correspond pas exactement au umask_i et vmask_i !
543      tvolu = 0._wp
544      tvolv = 0._wp
545
546      DO jk = 1, jpk
547         DO jj = 2, jpjm1
548            DO ji = fs_2, fs_jpim1   ! vector opt.
549               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)
550               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)
551            END DO
552         END DO
553      END DO
554      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
555      IF( lk_mpp )   CALL mpp_sum( tvolv )
556
557      IF(lwp) THEN
558         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
559         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
560      ENDIF
561      !
562   END SUBROUTINE trd_glo_init
563
564   !!======================================================================
565END MODULE trdglo
Note: See TracBrowser for help on using the repository browser.