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.
trdicp.F90 in trunk/NEMO/OPA_SRC/TRD – NEMO

source: trunk/NEMO/OPA_SRC/TRD/trdicp.F90 @ 1708

Last change on this file since 1708 was 1708, checked in by rblod, 14 years ago

Stability for explicit bottom friction, see ticket #588
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.6 KB
Line 
1MODULE trdicp
2   !!======================================================================
3   !!                       ***  MODULE  trdicp  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic trends
5   !!=====================================================================
6   !! History :       !  91-12 (G. Madec)
7   !!                 !  92-06 (M. Imbard) add time step frequency
8   !!                 !  96-01 (G. Madec)  terrain following coordinates
9   !!            8.5  !  02-06 (G. Madec)  F90: Free form and module
10   !!            9.0  !  04-08 (C. Talandier) New trends organization
11   !!----------------------------------------------------------------------
12#if  defined key_trdtra   ||   defined key_trddyn   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_trdtra'  or                  active tracers trends diagnostics
15   !!   'key_trddyn'                            momentum trends diagnostics
16   !!----------------------------------------------------------------------
17   !!----------------------------------------------------------------------
18   !!   trd_icp          : compute the basin averaged properties for tra/dyn
19   !!   trd_dwr          : print dynmaic trends in ocean.output file
20   !!   trd_twr          : print tracers trends in ocean.output file
21   !!   trd_icp_init     : initialization step
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain variables
25   USE trdmod_oce      ! ocean variables trends
26   USE ldftra_oce      ! ocean active tracers: lateral physics
27   USE ldfdyn_oce      ! ocean dynamics: lateral physics
28   USE zdf_oce         ! ocean vertical physics
29   USE in_out_manager  ! I/O manager
30   USE lib_mpp         ! distibuted memory computing library
31   USE eosbn2          ! equation of state
32   USE phycst          ! physical constants
33
34   IMPLICIT NONE
35   PRIVATE
36
37   INTERFACE trd_icp
38      MODULE PROCEDURE trd_2d, trd_3d
39   END INTERFACE
40
41   PUBLIC   trd_icp       ! called by trdmod.F90
42   PUBLIC   trd_dwr       ! called by step.F90
43   PUBLIC   trd_twr       ! called by step.F90
44   PUBLIC   trd_icp_init  ! called by opa.F90
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !!   OPA 9.0 , LOCEAN-IPSL (2005)
51   !! $Id$
52   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )
58      !!---------------------------------------------------------------------
59      !!                  ***  ROUTINE trd_2d  ***
60      !!
61      !! ** Purpose : verify the basin averaged properties of tracers and/or
62      !!              momentum equations at every time step frequency nn_trd.
63      !!----------------------------------------------------------------------
64      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dx             ! Temperature or U trend
65      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptrd2dy             ! Salinity    or V trend
66      INTEGER                     , INTENT(in   ) ::   ktrd                ! tracer trend index
67      CHARACTER(len=3)            , INTENT(in   ) ::   ctype               ! momentum ('DYN') or tracers ('TRA') trends
68      CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   clpas     ! number of passage
69      !!
70      INTEGER  ::   ji, jj                                                 ! loop indices
71      CHARACTER(len=3) ::   cpas                                           ! number of passage
72      REAL(wp) ::   zmsku, zbtu, zbt                                       ! temporary scalars
73      REAL(wp) ::   zmskv, zbtv                                            !    "         "
74      !!----------------------------------------------------------------------
75
76      ! Control of optional arguments
77      cpas = 'fst'
78      IF( PRESENT(clpas) )  cpas = clpas
79
80      ! 1. Mask trends
81      ! --------------
82
83      SELECT CASE( ctype )
84      !
85      CASE( 'DYN' )              ! Momentum
86         DO jj = 1, jpjm1
87            DO ji = 1, jpim1
88               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,1)
89               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,1)
90               ptrd2dx(ji,jj) = ptrd2dx(ji,jj) * zmsku
91               ptrd2dy(ji,jj) = ptrd2dy(ji,jj) * zmskv
92            END DO
93         END DO
94         ptrd2dx(jpi, : ) = 0.e0      ;      ptrd2dy(jpi, : ) = 0.e0
95         ptrd2dx( : ,jpj) = 0.e0      ;      ptrd2dy( : ,jpj) = 0.e0
96         !
97      CASE( 'TRA' )              ! Tracers
98         ptrd2dx(:,:) = ptrd2dx(:,:) * tmask_i(:,:)
99         ptrd2dy(:,:) = ptrd2dy(:,:) * tmask_i(:,:)
100         !
101      END SELECT
102     
103      ! 2. Basin averaged tracer/momentum trends
104      ! ----------------------------------------
105
106      SELECT CASE( ctype )
107      !
108      CASE( 'DYN' )              ! Momentum
109         umo(ktrd) = 0.e0
110         vmo(ktrd) = 0.e0
111         !
112         SELECT CASE( ktrd )
113         !
114         CASE( jpdyn_trd_swf )         ! surface forcing
115            DO jj = 1, jpj
116               DO ji = 1, jpi
117                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
118                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
119               END DO
120            END DO
121            !
122         END SELECT
123         !
124      CASE( 'TRA' )              ! Tracers
125         IF( cpas == 'fst' )   THEN
126            tmo(ktrd) = 0.e0
127            smo(ktrd) = 0.e0
128         ENDIF
129         DO jj = 1, jpj
130            DO ji = 1, jpi
131               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
132               tmo(ktrd) =  tmo(ktrd) + ptrd2dx(ji,jj) * zbt
133               smo(ktrd) =  smo(ktrd) + ptrd2dy(ji,jj) * zbt
134            END DO
135         END DO
136         !
137      END SELECT
138     
139      ! 3. Basin averaged tracer/momentum square trends
140      ! ----------------------------------------------
141      ! c a u t i o n: field now
142     
143      SELECT CASE( ctype )
144      !
145      CASE( 'DYN' )              ! Momentum
146         hke(ktrd) = 0.e0
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
150               zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
151               hke(ktrd) = hke(ktrd)   &
152               &   + un(ji,jj,1) * ptrd2dx(ji,jj) * zbtu &
153               &   + vn(ji,jj,1) * ptrd2dy(ji,jj) * zbtv
154            END DO
155         END DO
156         !
157      CASE( 'TRA' )              ! Tracers
158         IF( cpas == 'fst' )   THEN
159            t2(ktrd) = 0.e0
160            s2(ktrd) = 0.e0
161         ENDIF
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1)
165               t2(ktrd) = t2(ktrd) + ptrd2dx(ji,jj) * zbt * tn(ji,jj,1)
166               s2(ktrd) = s2(ktrd) + ptrd2dy(ji,jj) * zbt * sn(ji,jj,1)
167            END DO
168         END DO
169         !     
170      END SELECT
171      !
172   END SUBROUTINE trd_2d
173
174
175   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd, ctype, clpas )
176      !!---------------------------------------------------------------------
177      !!                  ***  ROUTINE trd_3d  ***
178      !!
179      !! ** Purpose : verify the basin averaged properties of tracers and/or
180      !!              momentum equations at every time step frequency nn_trd.
181      !!----------------------------------------------------------------------
182      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dx            ! Temperature or U trend
183      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptrd3dy            ! Salinity    or V trend
184      INTEGER,                          INTENT(in   ) ::   ktrd               ! momentum or tracer trend index
185      CHARACTER(len=3),                 INTENT(in   ) ::   ctype              ! momentum ('DYN') or tracers ('TRA') trends
186      CHARACTER(len=3),                 INTENT(in   ), OPTIONAL ::   clpas    ! number of passage
187      !!
188      INTEGER ::   ji, jj, jk
189      CHARACTER(len=3) ::   cpas                                              ! number of passage
190      REAL(wp) ::   zbt, zbtu, zbtv, zmsku, zmskv                             ! temporary scalars
191      !!----------------------------------------------------------------------
192
193      ! Control of optional arguments
194      cpas = 'fst'
195      IF( PRESENT(clpas) )  cpas = clpas
196
197      ! 1. Mask the trends
198      ! ------------------
199
200      SELECT CASE( ctype )
201      !
202      CASE( 'DYN' )              ! Momentum       
203         DO jk = 1, jpk
204            DO jj = 1, jpjm1
205               DO ji = 1, jpim1
206                  zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
207                  zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
208                  ptrd3dx(ji,jj,jk) = ptrd3dx(ji,jj,jk) * zmsku
209                  ptrd3dy(ji,jj,jk) = ptrd3dy(ji,jj,jk) * zmskv
210               END DO
211            END DO
212         END DO
213         ptrd3dx(jpi, : ,:) = 0.e0      ;      ptrd3dy(jpi, : ,:) = 0.e0
214         ptrd3dx( : ,jpj,:) = 0.e0      ;      ptrd3dy( : ,jpj,:) = 0.e0
215         !
216      CASE( 'TRA' )              ! Tracers
217         DO jk = 1, jpk
218            ptrd3dx(:,:,jk) = ptrd3dx(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
219            ptrd3dy(:,:,jk) = ptrd3dy(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
220         END DO
221         !
222      END SELECT   
223
224      ! 2. Basin averaged tracer/momentum trends
225      ! ----------------------------------------
226     
227      SELECT CASE( ctype )
228      !
229      CASE( 'DYN' )              ! Momentum
230         umo(ktrd) = 0.e0
231         vmo(ktrd) = 0.e0
232         DO jk = 1, jpk
233            DO jj = 1, jpj
234               DO ji = 1, jpi
235                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
236                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
237                  umo(ktrd) = umo(ktrd) + ptrd3dx(ji,jj,jk) * zbtu
238                  vmo(ktrd) = vmo(ktrd) + ptrd3dy(ji,jj,jk) * zbtv
239               END DO
240            END DO
241         END DO
242         !
243      CASE( 'TRA' )              ! Tracers
244         IF( cpas == 'fst' )   THEN
245            tmo(ktrd) = 0.e0
246            smo(ktrd) = 0.e0
247         ENDIF
248         DO jk = 1, jpkm1
249            DO jj = 1, jpj
250               DO ji = 1, jpi
251                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
252                  tmo(ktrd) =  tmo(ktrd) + ptrd3dx(ji,jj,jk) * zbt
253                  smo(ktrd) =  smo(ktrd) + ptrd3dy(ji,jj,jk) * zbt
254               END DO
255            END DO
256         END DO
257         !
258      END SELECT
259
260      ! 3. Basin averaged tracer/momentum square trends
261      ! -----------------------------------------------
262      ! c a u t i o n: field now
263     
264      SELECT CASE( ctype )
265      !
266      CASE( 'DYN' )              ! Momentum
267         hke(ktrd) = 0.e0
268         DO jk = 1, jpk
269            DO jj = 1, jpj
270               DO ji = 1, jpi
271                  zbtu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
272                  zbtv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
273                  hke(ktrd) = hke(ktrd)   &
274                  &   + un(ji,jj,jk) * ptrd3dx(ji,jj,jk) * zbtu &
275                  &   + vn(ji,jj,jk) * ptrd3dy(ji,jj,jk) * zbtv
276               END DO
277            END DO
278         END DO
279         !
280      CASE( 'TRA' )              ! Tracers
281         IF( cpas == 'fst' )   THEN
282            t2(ktrd) = 0.e0
283            s2(ktrd) = 0.e0
284         ENDIF
285         DO jk = 1, jpk
286            DO jj = 1, jpj
287               DO ji = 1, jpi
288                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
289                  t2(ktrd) = t2(ktrd) + ptrd3dx(ji,jj,jk) * zbt * tn(ji,jj,jk)
290                  s2(ktrd) = s2(ktrd) + ptrd3dy(ji,jj,jk) * zbt * sn(ji,jj,jk)
291               END DO
292            END DO
293         END DO
294         !
295      END SELECT
296      !
297   END SUBROUTINE trd_3d
298
299
300
301   SUBROUTINE trd_icp_init
302      !!---------------------------------------------------------------------
303      !!                  ***  ROUTINE trd_icp_init  ***
304      !!
305      !! ** Purpose :   Read the namtrd namelist
306      !!----------------------------------------------------------------------
307      INTEGER  ::   ji, jj, jk
308      REAL(wp) ::   zmskt
309#if  defined key_trddyn
310      REAL(wp) ::   zmsku, zmskv
311#endif
312      !!----------------------------------------------------------------------
313
314      IF(lwp) THEN
315         WRITE(numout,*)
316         WRITE(numout,*) 'trd_icp_init : integral constraints properties trends'
317         WRITE(numout,*) '~~~~~~~~~~~~~'
318      ENDIF
319
320      ! Total volume at t-points:
321      tvolt = 0.e0
322      DO jk = 1, jpkm1
323         DO jj = 2, jpjm1
324            DO ji = fs_2, fs_jpim1   ! vector opt.
325               zmskt = tmask(ji,jj,jk) * tmask_i(ji,jj)
326               tvolt = tvolt + zmskt * e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk)
327            END DO
328         END DO
329      END DO
330      IF( lk_mpp )   CALL mpp_sum( tvolt )   ! sum over the global domain
331
332      IF(lwp) WRITE(numout,*) '                total ocean volume at T-point   tvolt = ',tvolt
333
334#if  defined key_trddyn
335      ! Initialization of potential to kinetic energy conversion
336      rpktrd = 0.e0
337
338      ! Total volume at u-, v- points:
339      tvolu = 0.e0
340      tvolv = 0.e0
341
342      DO jk = 1, jpk
343         DO jj = 2, jpjm1
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
346               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
347               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
348               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
349            END DO
350         END DO
351      END DO
352      IF( lk_mpp )   CALL mpp_sum( tvolu )   ! sums over the global domain
353      IF( lk_mpp )   CALL mpp_sum( tvolv )
354
355      IF(lwp) THEN
356         WRITE(numout,*) '                total ocean volume at U-point   tvolu = ',tvolu
357         WRITE(numout,*) '                total ocean volume at V-point   tvolv = ',tvolv
358      ENDIF
359#endif
360      !
361   END SUBROUTINE trd_icp_init
362
363
364   SUBROUTINE trd_dwr( kt )
365      !!---------------------------------------------------------------------
366      !!                  ***  ROUTINE trd_dwr  ***
367      !!
368      !! ** Purpose :  write dynamic trends in ocean.output
369      !!----------------------------------------------------------------------
370      INTEGER, INTENT(in) ::   kt                                  ! ocean time-step index
371      !!
372      INTEGER  ::   ji, jj, jk
373      REAL(wp) ::   ze1e2w, zcof, zbe1ru, zbe2rv, zbtr, ztz, zth   !    "      scalars
374      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkepe, zkx, zky, zkz   ! temporary arrays
375      !!----------------------------------------------------------------------
376
377      ! I. Momentum trends
378      ! -------------------
379
380      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
381
382         ! I.1 Conversion potential energy - kinetic energy
383         ! --------------------------------------------------
384         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
385
386         zkx(:,:,:)   = 0.e0
387         zky(:,:,:)   = 0.e0
388         zkz(:,:,:)   = 0.e0
389         zkepe(:,:,:) = 0.e0
390   
391         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
392
393         ! Density flux at w-point
394         DO jk = 2, jpk
395            DO jj = 1, jpj
396               DO ji = 1, jpi
397                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
398                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
399               END DO
400            END DO
401         END DO
402         zkz(:,:,1) = 0.e0
403         
404         ! Density flux at u and v-points
405         DO jk = 1, jpk
406            DO jj = 1, jpjm1
407               DO ji = 1, jpim1
408                  zcof   = 0.5 / rau0
409                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
410                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
411                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
412                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
413               END DO
414            END DO
415         END DO
416         
417         ! Density flux divergence at t-point
418         DO jk = 1, jpkm1
419            DO jj = 2, jpjm1
420               DO ji = 2, jpim1
421                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
422                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
423                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
424                    &             + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
425                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
426               END DO
427            END DO
428         END DO
429         zkepe( : , : ,jpk) = 0.e0
430         zkepe( : ,jpj, : ) = 0.e0
431         zkepe(jpi, : , : ) = 0.e0
432
433         ! I.2 Basin averaged kinetic energy trend
434         ! ----------------------------------------
435         peke = 0.e0
436         DO jk = 1,jpk
437            DO jj = 1, jpj
438               DO ji = 1, jpi
439                  peke = peke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk)   &
440                     &                     * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
441               END DO
442            END DO
443         END DO
444
445         ! I.3 Sums over the global domain
446         ! ---------------------------------
447         IF( lk_mpp ) THEN
448            CALL mpp_sum( peke )
449            CALL mpp_sum( umo , jptot_dyn )
450            CALL mpp_sum( vmo , jptot_dyn )
451            CALL mpp_sum( hke , jptot_dyn )
452         ENDIF
453
454         ! I.2 Print dynamic trends in the ocean.output file
455         ! --------------------------------------------------
456
457         IF(lwp) THEN
458            WRITE (numout,*)
459            WRITE (numout,*)
460            WRITE (numout,9500) kt
461            WRITE (numout,9501) umo(jpicpd_hpg) / tvolu, vmo(jpicpd_hpg) / tvolv
462            WRITE (numout,9502) umo(jpicpd_keg) / tvolu, vmo(jpicpd_keg) / tvolv
463            WRITE (numout,9503) umo(jpicpd_rvo) / tvolu, vmo(jpicpd_rvo) / tvolv
464            WRITE (numout,9504) umo(jpicpd_pvo) / tvolu, vmo(jpicpd_pvo) / tvolv
465            WRITE (numout,9505) umo(jpicpd_ldf) / tvolu, vmo(jpicpd_ldf) / tvolv
466            WRITE (numout,9506) umo(jpicpd_had) / tvolu, vmo(jpicpd_had) / tvolv
467            WRITE (numout,9507) umo(jpicpd_zad) / tvolu, vmo(jpicpd_zad) / tvolv
468            WRITE (numout,9508) umo(jpicpd_zdf) / tvolu, vmo(jpicpd_zdf) / tvolv
469            WRITE (numout,9509) umo(jpicpd_spg) / tvolu, vmo(jpicpd_spg) / tvolv
470            WRITE (numout,9510) umo(jpicpd_swf) / tvolu, vmo(jpicpd_swf) / tvolv
471            WRITE (numout,9511) umo(jpicpd_dat) / tvolu, vmo(jpicpd_dat) / tvolv
472            WRITE (numout,9512) umo(jpicpd_bfr) / tvolu, vmo(jpicpd_bfr) / tvolv
473            WRITE (numout,9513)
474            WRITE (numout,9514)                                                 &
475            &     (  umo(jpicpd_hpg) + umo(jpicpd_keg) + umo(jpicpd_rvo) + umo(jpicpd_pvo) + umo(jpicpd_ldf)   &
476            &      + umo(jpicpd_had) + umo(jpicpd_zad) + umo(jpicpd_zdf) + umo(jpicpd_spg) + umo(jpicpd_dat)   &
477            &      + umo(jpicpd_swf) + umo(jpicpd_bfr) ) / tvolu,   &
478            &     (  vmo(jpicpd_hpg) + vmo(jpicpd_keg) + vmo(jpicpd_rvo) + vmo(jpicpd_pvo) + vmo(jpicpd_ldf)   &
479            &      + vmo(jpicpd_had) + vmo(jpicpd_zad) + vmo(jpicpd_zdf) + vmo(jpicpd_spg) + vmo(jpicpd_dat)   &
480            &      + vmo(jpicpd_swf) + vmo(jpicpd_bfr) ) / tvolv
481         ENDIF
482
483 9500    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
484 9501    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
485 9502    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
486 9503    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
487 9504    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
488 9505    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
489 9506    FORMAT(' horizontal advection       u= ', e20.13, '    v= ', e20.13)
490 9507    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
491 9508    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
492 9509    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
493 9510    FORMAT(' surface wind forcing       u= ', e20.13, '    v= ', e20.13)
494 9511    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
495 9512    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
496 9513    FORMAT(' -----------------------------------------------------------------------------')
497 9514    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
498
499         IF(lwp) THEN
500            WRITE (numout,*)
501            WRITE (numout,*)
502            WRITE (numout,9520) kt
503            WRITE (numout,9521) hke(jpicpd_hpg) / tvolt
504            WRITE (numout,9522) hke(jpicpd_keg) / tvolt
505            WRITE (numout,9523) hke(jpicpd_rvo) / tvolt
506            WRITE (numout,9524) hke(jpicpd_pvo) / tvolt
507            WRITE (numout,9525) hke(jpicpd_ldf) / tvolt
508            WRITE (numout,9526) hke(jpicpd_had) / tvolt
509            WRITE (numout,9527) hke(jpicpd_zad) / tvolt
510            WRITE (numout,9528) hke(jpicpd_zdf) / tvolt
511            WRITE (numout,9529) hke(jpicpd_spg) / tvolt
512            WRITE (numout,9530) hke(jpicpd_swf) / tvolt
513            WRITE (numout,9531) hke(jpicpd_dat) / tvolt
514            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt
515            WRITE (numout,9533)
516            WRITE (numout,9534)   &
517            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   &
518            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   &
519            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt
520         ENDIF
521
522 9520    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
523 9521    FORMAT(' pressure gradient         u2= ', e20.13)
524 9522    FORMAT(' ke gradient               u2= ', e20.13)
525 9523    FORMAT(' relative vorticity term   u2= ', e20.13)
526 9524    FORMAT(' coriolis term             u2= ', e20.13)
527 9525    FORMAT(' horizontal diffusion      u2= ', e20.13)
528 9526    FORMAT(' horizontal advection      u2= ', e20.13)
529 9527    FORMAT(' vertical advection        u2= ', e20.13)
530 9528    FORMAT(' vertical diffusion        u2= ', e20.13)
531 9529    FORMAT(' surface pressure gradient u2= ', e20.13)
532 9530    FORMAT(' surface wind forcing      u2= ', e20.13)
533 9531    FORMAT(' dampimg term              u2= ', e20.13)
534 9532    FORMAT(' bottom flux               u2= ', e20.13)
535 9533    FORMAT(' --------------------------------------------------')
536 9534    FORMAT(' total trend               u2= ', e20.13)
537
538         IF(lwp) THEN
539            WRITE (numout,*)
540            WRITE (numout,*)
541            WRITE (numout,9540) kt
542            WRITE (numout,9541) ( hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
543            WRITE (numout,9542) ( hke(jpicpd_keg) + hke(jpicpd_had) + hke(jpicpd_zad) ) / tvolt
544            WRITE (numout,9543) ( hke(jpicpd_pvo) ) / tvolt
545            WRITE (numout,9544) ( hke(jpicpd_rvo) ) / tvolt
546            WRITE (numout,9545) ( hke(jpicpd_spg) ) / tvolt
547            WRITE (numout,9546) ( hke(jpicpd_ldf) ) / tvolt
548            WRITE (numout,9547) ( hke(jpicpd_zdf) ) / tvolt
549            WRITE (numout,9548) ( hke(jpicpd_hpg) ) / tvolt, rpktrd / tvolt
550            WRITE (numout,*)
551            WRITE (numout,*)
552         ENDIF
553
554 9540    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
555 9541    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
556 9542    FORMAT(' 0 = ke gradient + horizontal + vertical advection         : ', e20.13)
557 9543    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
558 9544    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)        : ', e20.13)
559 9545    FORMAT(' 0 = surface pressure gradient                             : ', e20.13)
560 9546    FORMAT(' 0 > horizontal diffusion                                  : ', e20.13)
561 9547    FORMAT(' 0 > vertical diffusion                                    : ', e20.13)
562 9548    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)                : ', e20.13, '  u.dz(rhop) =', e20.13)
563         !
564         ! Save potential to kinetic energy conversion for next time step
565         rpktrd = peke
566         !
567      ENDIF
568      !
569   END SUBROUTINE trd_dwr
570
571
572   SUBROUTINE trd_twr( kt )
573      !!---------------------------------------------------------------------
574      !!                  ***  ROUTINE trd_twr  ***
575      !!
576      !! ** Purpose :  write active tracers trends in ocean.output
577      !!----------------------------------------------------------------------
578      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
579      !!----------------------------------------------------------------------
580
581      ! I. Tracers trends
582      ! -----------------
583
584      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
585
586         ! I.1 Sums over the global domain
587         ! -------------------------------
588         IF( lk_mpp ) THEN
589            CALL mpp_sum( tmo, jptot_tra )   
590            CALL mpp_sum( smo, jptot_tra )
591            CALL mpp_sum( t2 , jptot_tra )
592            CALL mpp_sum( s2 , jptot_tra )
593         ENDIF
594
595         ! I.2 Print tracers trends in the ocean.output file
596         ! --------------------------------------------------
597         
598         IF(lwp) THEN
599            WRITE (numout,*)
600            WRITE (numout,*)
601            WRITE (numout,9400) kt
602            WRITE (numout,9401) (tmo(jpicpt_xad)+tmo(jpicpt_yad))/ tvolt, (smo(jpicpt_xad)+smo(jpicpt_yad))/ tvolt
603            WRITE (numout,9402)  tmo(jpicpt_zad) / tvolt, smo(jpicpt_zad) / tvolt
604            WRITE (numout,9403)  tmo(jpicpt_ldf) / tvolt, smo(jpicpt_ldf) / tvolt
605            WRITE (numout,9404)  tmo(jpicpt_zdf) / tvolt, smo(jpicpt_zdf) / tvolt
606            WRITE (numout,9405)  tmo(jpicpt_npc) / tvolt, smo(jpicpt_npc) / tvolt
607            WRITE (numout,9406)  tmo(jpicpt_dmp) / tvolt, smo(jpicpt_dmp) / tvolt
608            WRITE (numout,9407)  tmo(jpicpt_qsr) / tvolt
609            WRITE (numout,9408)  tmo(jpicpt_nsr) / tvolt, smo(jpicpt_nsr) / tvolt
610            WRITE (numout,9409) 
611            WRITE (numout,9410) (  tmo(jpicpt_xad) + tmo(jpicpt_yad) + tmo(jpicpt_zad) + tmo(jpicpt_ldf) + tmo(jpicpt_zdf)   &
612            &                    + tmo(jpicpt_npc) + tmo(jpicpt_dmp) + tmo(jpicpt_qsr) + tmo(jpicpt_nsr) ) / tvolt,   &
613            &                   (  smo(jpicpt_xad) + smo(jpicpt_yad) + smo(jpicpt_zad) + smo(jpicpt_ldf) + smo(jpicpt_zdf)   &
614            &                    + smo(jpicpt_npc) + smo(jpicpt_dmp)                   + smo(jpicpt_nsr) ) / tvolt
615         ENDIF
616
6179400     FORMAT(' tracer trend at it= ',i6,' :     temperature',   &
618              '              salinity',/' ============================')
6199401     FORMAT(' horizontal advection        ',e20.13,'     ',e20.13)
6209402     FORMAT(' vertical advection          ',e20.13,'     ',e20.13)
6219403     FORMAT(' horizontal diffusion        ',e20.13,'     ',e20.13)
6229404     FORMAT(' vertical diffusion          ',e20.13,'     ',e20.13)
6239405     FORMAT(' static instability mixing   ',e20.13,'     ',e20.13)
6249406     FORMAT(' damping term                ',e20.13,'     ',e20.13)
6259407     FORMAT(' penetrative qsr             ',e20.13)
6269408     FORMAT(' non solar radiation         ',e20.13,'     ',e20.13)
6279409     FORMAT(' -------------------------------------------------------------------------')
6289410     FORMAT(' total trend                 ',e20.13,'     ',e20.13)
629
630
631         IF(lwp) THEN
632            WRITE (numout,*)
633            WRITE (numout,*)
634            WRITE (numout,9420) kt
635            WRITE (numout,9421) ( t2(jpicpt_xad)+t2(jpicpt_yad) )/ tvolt, ( s2(jpicpt_xad)+s2(jpicpt_yad) )/ tvolt
636            WRITE (numout,9422)   t2(jpicpt_zad) / tvolt, s2(jpicpt_zad) / tvolt
637            WRITE (numout,9423)   t2(jpicpt_ldf) / tvolt, s2(jpicpt_ldf) / tvolt
638            WRITE (numout,9424)   t2(jpicpt_zdf) / tvolt, s2(jpicpt_zdf) / tvolt
639            WRITE (numout,9425)   t2(jpicpt_npc) / tvolt, s2(jpicpt_npc) / tvolt
640            WRITE (numout,9426)   t2(jpicpt_dmp) / tvolt, s2(jpicpt_dmp) / tvolt
641            WRITE (numout,9427)   t2(jpicpt_qsr) / tvolt
642            WRITE (numout,9428)   t2(jpicpt_nsr) / tvolt, s2(jpicpt_nsr) / tvolt
643            WRITE (numout,9429)
644            WRITE (numout,9430) (  t2(jpicpt_xad) + t2(jpicpt_yad) + t2(jpicpt_zad) + t2(jpicpt_ldf) + t2(jpicpt_zdf)   &
645            &                    + t2(jpicpt_npc) + t2(jpicpt_dmp) + t2(jpicpt_qsr) + t2(jpicpt_nsr) ) / tvolt,   &
646            &                   (  s2(jpicpt_xad) + s2(jpicpt_yad) + s2(jpicpt_zad) + s2(jpicpt_ldf) + s2(jpicpt_zdf)   &
647            &                    + s2(jpicpt_npc) + s2(jpicpt_dmp)                  + s2(jpicpt_nsr) ) / tvolt
648         ENDIF
649
6509420     FORMAT(' tracer**2 trend at it= ', i6, ' :      temperature',   &
651            '               salinity', /, ' ===============================')
6529421     FORMAT(' horizontal advection      * t   ', e20.13, '     ', e20.13)
6539422     FORMAT(' vertical advection        * t   ', e20.13, '     ', e20.13)
6549423     FORMAT(' horizontal diffusion      * t   ', e20.13, '     ', e20.13)
6559424     FORMAT(' vertical diffusion        * t   ', e20.13, '     ', e20.13)
6569425     FORMAT(' static instability mixing * t   ', e20.13, '     ', e20.13)
6579426     FORMAT(' damping term              * t   ', e20.13, '     ', e20.13)
6589427     FORMAT(' penetrative qsr           * t   ', e20.13)
6599428     FORMAT(' non solar radiation       * t   ', e20.13, '     ', e20.13)
6609429     FORMAT(' -----------------------------------------------------------------------------')
6619430     FORMAT(' total trend                *t = ', e20.13, '  *s = ', e20.13)
662
663
664         IF(lwp) THEN
665            WRITE (numout,*)
666            WRITE (numout,*)
667            WRITE (numout,9440) kt
668            WRITE (numout,9441) ( tmo(jpicpt_xad)+tmo(jpicpt_yad)+tmo(jpicpt_zad) )/tvolt,   &
669            &                   ( smo(jpicpt_xad)+smo(jpicpt_yad)+smo(jpicpt_zad) )/tvolt
670            WRITE (numout,9442)   tmo(jpicpt_zl1)/tvolt,  smo(jpicpt_zl1)/tvolt
671            WRITE (numout,9443)   tmo(jpicpt_ldf)/tvolt,  smo(jpicpt_ldf)/tvolt
672            WRITE (numout,9444)   tmo(jpicpt_zdf)/tvolt,  smo(jpicpt_zdf)/tvolt
673            WRITE (numout,9445)   tmo(jpicpt_npc)/tvolt,  smo(jpicpt_npc)/tvolt
674            WRITE (numout,9446) ( t2(jpicpt_xad)+t2(jpicpt_yad)+t2(jpicpt_zad) )/tvolt,    &
675            &                   ( s2(jpicpt_xad)+s2(jpicpt_yad)+s2(jpicpt_zad) )/tvolt
676            WRITE (numout,9447)   t2(jpicpt_ldf)/tvolt,   s2(jpicpt_ldf)/tvolt
677            WRITE (numout,9448)   t2(jpicpt_zdf)/tvolt,   s2(jpicpt_zdf)/tvolt
678            WRITE (numout,9449)   t2(jpicpt_npc)/tvolt,   s2(jpicpt_npc)/tvolt
679         ENDIF
680
6819440     FORMAT(' tracer consistency at it= ',i6,   &
682            ' :         temperature','                salinity',/,   &
683            ' ==================================')
6849441     FORMAT(' 0 = horizontal+vertical advection +    ',e20.13,'       ',e20.13)
6859442     FORMAT('     1st lev vertical advection         ',e20.13,'       ',e20.13)
6869443     FORMAT(' 0 = horizontal diffusion               ',e20.13,'       ',e20.13)
6879444     FORMAT(' 0 = vertical diffusion                 ',e20.13,'       ',e20.13)
6889445     FORMAT(' 0 = static instability mixing          ',e20.13,'       ',e20.13)
6899446     FORMAT(' 0 = horizontal+vertical advection * t  ',e20.13,'       ',e20.13)
6909447     FORMAT(' 0 > horizontal diffusion          * t  ',e20.13,'       ',e20.13)
6919448     FORMAT(' 0 > vertical diffusion            * t  ',e20.13,'       ',e20.13)
6929449     FORMAT(' 0 > static instability mixing     * t  ',e20.13,'       ',e20.13)
693         !
694      ENDIF
695      !
696   END SUBROUTINE trd_twr
697
698#   else
699   !!----------------------------------------------------------------------
700   !!   Default case :                                         Empty module
701   !!----------------------------------------------------------------------
702   INTERFACE trd_icp
703      MODULE PROCEDURE trd_2d, trd_3d
704   END INTERFACE
705
706CONTAINS
707   SUBROUTINE trd_2d( ptrd2dx, ptrd2dy, ktrd , ctype, clpas )       ! Empty routine
708      REAL, DIMENSION(:,:) ::   ptrd2dx, ptrd2dy
709      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
710      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
711      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
712      WRITE(*,*) 'trd_2d: You should not have seen this print! error ?', &
713          &       ptrd2dx(1,1), ptrd2dy(1,1), ktrd, ctype, clpas
714   END SUBROUTINE trd_2d
715   SUBROUTINE trd_3d( ptrd3dx, ptrd3dy, ktrd , ctype, clpas )       ! Empty routine
716      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy
717      INTEGER                     , INTENT(in   ) ::   ktrd         ! tracer trend index
718      CHARACTER(len=3)            , INTENT(in   ) ::   ctype        ! momentum ('DYN') or tracers ('TRA') trends
719      CHARACTER(len=3), INTENT(in), OPTIONAL ::   clpas             ! number of passage
720      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', &
721          &       ptrd3dx(1,1,1), ptrd3dy(1,1,1), ktrd, ctype, clpas
722   END SUBROUTINE trd_3d
723   SUBROUTINE trd_icp_init               ! Empty routine
724   END SUBROUTINE trd_icp_init
725   SUBROUTINE trd_dwr( kt )          ! Empty routine
726      WRITE(*,*) 'trd_dwr: You should not have seen this print! error ?', kt
727   END SUBROUTINE trd_dwr
728   SUBROUTINE trd_twr( kt )          ! Empty routine
729      WRITE(*,*) 'trd_twr: You should not have seen this print! error ?', kt
730   END SUBROUTINE trd_twr
731#endif
732
733   !!======================================================================
734END MODULE trdicp
Note: See TracBrowser for help on using the repository browser.