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

source: trunk/NEMO/OPA_SRC/TRD/trddyn.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.5 KB
Line 
1MODULE trddyn
2   !!======================================================================
3   !!                       ***  MODULE  trddyn  ***
4   !! Ocean trends :  ocean momentum trends
5   !!=====================================================================
6#if   defined key_trddyn   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_trddyn'                                   momentum trend diag.
9   !!----------------------------------------------------------------------
10   !!   trd_dyn      : verify the basin averaged properties for momentum
11   !!   trd_dyn_init :
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE phycst          ! physical constants
17   USE trdtra_oce      ! ocean active tracer trend
18   USE trddyn_oce      ! ocean dynamics trend
19   USE zdf_oce         ! ocean vertical physics
20   USE ldftra_oce      ! ocean active tracers: lateral physics
21   USE ldfdyn_oce      ! ocean dynamics: lateral physics
22   USE sol_oce         ! solver variables
23   USE obc_oce         ! ocean lateral open boundary condition
24   USE eosbn2          ! equation of state (eos routine)
25   USE in_out_manager  ! I/O manager
26   USE lib_mpp         ! distributed memory computing library
27
28   IMPLICIT NONE
29   PRIVATE
30
31   !! * Routine accessibility
32   PUBLIC trd_dyn      ! called by step.F90
33   PUBLIC trd_dyn_init ! called by step.F90
34
35   !! * Shared module vaiables
36   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .TRUE.   ! momentum trend flag
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !!   OPA 9.0 , LODYC-IPSL  (2003)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46
47   SUBROUTINE trd_dyn( kt )
48      !!---------------------------------------------------------------------
49      !!                  ***  ROUTINE trd_dyn  ***
50      !!
51      !! ** Purpose :   verify the basin averaged properties of the momentum
52      !!      equation at every time step frequency ntrd.
53      !!      momentum equation:
54      !!       * non linear momentum trend (relative vorticity trend + hori-
55      !!      zontal kinetic energy gradient trend + vertical advection
56      !!      trend) :
57      !!      1. conserve the momentum
58      !!      2. conserve the potential relative vorticity (rotn/e3f)
59      !!         (except for the vertical advection)
60      !!      3. conserve the potential enstrophy (except for the
61      !!         vertical advection trend) IF no cpp key activated
62      !!        or #defined key_vorcombined
63      !!      4. conserve the horizontal kinetic energy if #defined
64      !!        key_vorenergy
65      !!
66      !! ** Method :
67      !!      damping trend
68      !!      horizontal pressure gradient trend
69      !!      horizontal kinetic energy gradient trend
70      !!      relative vorticity term trend
71      !!      coriolis trend
72      !!
73      !! History :
74      !!        !  91-12  (G. Madec)  Original code
75      !!        !  92-06  (M. Imbard)  add time step frequency
76      !!        !  96-01  (G. Madec)  Statement function for e3
77      !!                              suppression of common work arrays
78      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
79      !!---------------------------------------------------------------------
80      !! * Arguments
81      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
82
83      !! * Local declarations
84      INTEGER ::   ji, jj, jk, jn        ! dummy loop indices
85
86      REAL(wp) ::   ze1e2w,zbe1ru,zbe2rv,zbtr,zth,ztz,zpeke,zcof
87      REAL(wp) ::   zmsku, zmskv
88      REAL(wp) ::   zumo(11),zvmo(11),zhke(10)
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
90         zkepe, zkx, zky, zkz            ! workspace
91     
92      NAMELIST/namtrd/ ntrd, nctls
93      !!----------------------------------------------------------------------
94
95
96      ! 1. forcing trend and mask trends
97      ! --------------------------------
98     
99      IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
100
101         ! Mask surface forcing and bottom friction fluxes
102         DO jj = 1, jpjm1
103            DO ji = 1, jpim1
104               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
105               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
106               tautrd(ji,jj,1) = tautrd(ji,jj,1) * zmsku
107               tautrd(ji,jj,2) = tautrd(ji,jj,2) * zmskv
108               tautrd(ji,jj,3) = tautrd(ji,jj,3) * zmsku
109               tautrd(ji,jj,4) = tautrd(ji,jj,4) * zmskv
110            END DO
111         END DO
112         tautrd( : ,jpj,1:4) = 0.e0
113         tautrd(jpi, : ,1:4) = 0.e0
114         
115         ! Mask the trends
116         DO jn = 1,9
117            DO jk = 1,jpk
118               DO jj = 1, jpjm1
119                  DO ji = 1, jpim1
120                     zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
121                     zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
122                     utrd(ji,jj,jk,jn) = utrd(ji,jj,jk,jn) * zmsku
123                     vtrd(ji,jj,jk,jn) = vtrd(ji,jj,jk,jn) * zmskv
124                  END DO
125               END DO
126            END DO
127         END DO
128         utrd( : ,jpj,:,1:9) = 0.e0      ;      vtrd( : ,jpj,:,1:9) = 0.e0
129         utrd(jpi, : ,:,1:9) = 0.e0      ;      vtrd(jpi, : ,:,1:9) = 0.e0
130         
131         ! Conversion potential energy - kinetic energy
132         ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
133
134         zkx(:,:,:) = 0.e0
135         zky(:,:,:) = 0.e0
136         zkz(:,:,:) = 0.e0
137     
138         CALL eos( tn, sn, rhd, rhop )       ! now potential and in situ densities
139
140         ! Density flux at w-point
141         DO jk = 2, jpk
142            DO jj = 1, jpj
143               DO ji = 1, jpi
144                  ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj)
145                  zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) )
146               END DO
147            END DO
148         END DO
149         zkz  (:,:, 1 ) = 0.e0
150         
151         ! Density flux at u and v-points
152         DO jk = 1, jpk
153            DO jj = 1, jpjm1
154               DO ji = 1, jpim1
155                  zcof   = 0.5 / rau0
156                  zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
157                  zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
158                  zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
159                  zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
160               END DO
161            END DO
162         END DO
163         
164         ! Density flux divergence at t-point
165         DO jk = 1, jpkm1
166            DO jj = 1, jpjm1
167               DO ji = 1, jpim1
168                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
169                  ztz = - zbtr * (    zkz(ji,jj,jk) - zkz(ji,jj,jk+1) )
170                  zth = - zbtr * (  ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) )   &
171                                  + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) )  )
172                  zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj)
173               END DO
174            END DO
175         END DO
176         zkepe( : , : ,jpk) = 0.e0
177         zkepe( : ,jpj, : ) = 0.e0
178         zkepe(jpi, : , : ) = 0.e0
179
180
181         ! 2. Basin averaged momentum trend
182         ! --------------------------------
183         
184         DO jn = 1,9
185            zumo(jn) = 0.
186            zvmo(jn) = 0.
187            DO jk = 1, jpk
188               DO jj = 1, jpj
189                  DO ji = 1, jpi
190                     zumo(jn) = zumo(jn) + utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
191                     zvmo(jn) = zvmo(jn) + vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
192                  END DO
193               END DO
194            END DO
195         END DO
196
197         zumo(10) = 0.
198         zvmo(10) = 0.
199         zumo(11) = 0.
200         zvmo(11) = 0.
201         DO jj = 1, jpj
202            DO ji = 1, jpi
203               zumo(10) = zumo(10) + tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)
204               zvmo(10) = zvmo(10) + tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
205               zumo(11) = zumo(11) + tautrd(ji,jj,3)
206               zvmo(11) = zvmo(11) + tautrd(ji,jj,4)
207            END DO
208         END DO
209         
210         
211         ! 3. Basin averaged kinetic energy trend
212         ! --------------------------------------
213         
214         ! Field before, because after the array swap
215         
216         DO jn = 1,9
217            zhke(jn) = 0.e0
218            DO jk = 1, jpk
219               DO jj = 1, jpj
220                  DO ji = 1, jpi
221                     zhke(jn) = zhke(jn)   &
222                     &   + ub(ji,jj,jk) * utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)   &
223                     &   + vb(ji,jj,jk) * vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
224                  END DO
225               END DO
226            END DO
227         END DO
228         
229         zhke(10) = 0.e0
230         DO jj = 1, jpj
231            DO ji = 1, jpi
232               zhke(10) = zhke(10)   &
233               &   + ub(ji,jj,1) * tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1)   &
234               &   + vb(ji,jj,1) * tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1)
235            END DO
236         END DO
237         
238         zpeke = 0.e0
239         DO jk = 1,jpk
240            DO jj = 1, jpj
241               DO ji = 1, jpi
242                  zpeke    = zpeke + zkepe(ji,jj,jk) * g * fsdept(ji,jj,jk)   &
243                  &                        * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
244               END DO
245            END DO
246         END DO
247         
248# if defined key_mpp
249         CALL mpp_sum( zpeke )
250         CALL mpp_sum( zumo , 11 )
251         CALL mpp_sum( zvmo , 11 )
252         CALL mpp_sum( zhke , 10 )
253# endif
254
255
256         ! 4. Print
257         ! --------
258
259         IF(lwp) THEN
260            WRITE (numout,*)
261            WRITE (numout,*)
262            WRITE (numout,9400) kt
263            WRITE (numout,9401) zumo( 1) / tvolu, zvmo( 1) / tvolv
264            WRITE (numout,9402) zumo( 2) / tvolu, zvmo( 2) / tvolv
265            WRITE (numout,9403) zumo( 3) / tvolu, zvmo( 3) / tvolv
266            WRITE (numout,9404) zumo( 4) / tvolu, zvmo( 4) / tvolv
267            WRITE (numout,9405) zumo( 5) / tvolu, zvmo( 5) / tvolv
268            WRITE (numout,9406) zumo( 6) / tvolu, zvmo( 6) / tvolv
269            WRITE (numout,9407) zumo( 7) / tvolu, zvmo( 7) / tvolv
270            WRITE (numout,9408) zumo( 8) / tvolu, zvmo( 8) / tvolv
271            WRITE (numout,9409) zumo(10) / tvolu, zvmo(10) / tvolv
272            WRITE (numout,9410) zumo( 9) / tvolu, zvmo( 9) / tvolv
273            WRITE (numout,9411) zumo(11) / tvolu, zvmo(11) / tvolv
274            WRITE (numout,9412)
275            WRITE (numout,9413)                                                     &
276            &     (  zumo(1) + zumo(2) + zumo(3) + zumo(4) + zumo(5) + zumo(6)      &
277            &      + zumo(7) + zumo(8) + zumo(9) + zumo(10) + zumo(11) ) / tvolu,   &
278            &     (  zvmo(1) + zvmo(2) + zvmo(3) + zvmo(4) + zvmo(5) + zvmo(6)      &
279            &      + zvmo(7) + zvmo(8) + zvmo(9) + zvmo(10) + zvmo(11) ) / tvolv
280         ENDIF
281
282 9400    FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================')
283 9401    FORMAT(' pressure gradient          u= ', e20.13, '    v= ', e20.13)
284 9402    FORMAT(' ke gradient                u= ', e20.13, '    v= ', e20.13)
285 9403    FORMAT(' relative vorticity term    u= ', e20.13, '    v= ', e20.13)
286 9404    FORMAT(' coriolis term              u= ', e20.13, '    v= ', e20.13)
287 9405    FORMAT(' horizontal diffusion       u= ', e20.13, '    v= ', e20.13)
288 9406    FORMAT(' vertical advection         u= ', e20.13, '    v= ', e20.13)
289 9407    FORMAT(' vertical diffusion         u= ', e20.13, '    v= ', e20.13)
290 9408    FORMAT(' surface pressure gradient  u= ', e20.13, '    v= ', e20.13)
291 9409    FORMAT(' forcing term               u= ', e20.13, '    v= ', e20.13)
292 9410    FORMAT(' dampimg term               u= ', e20.13, '    v= ', e20.13)
293 9411    FORMAT(' bottom flux                u= ', e20.13, '    v= ', e20.13)
294 9412    FORMAT(' -----------------------------------------------------------------------------')
295 9413    FORMAT(' total trend                u= ', e20.13, '    v= ', e20.13)
296
297         IF(lwp) THEN
298            WRITE (numout,*)
299            WRITE (numout,*)
300            WRITE (numout,9420) kt
301            WRITE (numout,9421) zhke( 1) / tvols
302            WRITE (numout,9422) zhke( 2) / tvols
303            WRITE (numout,9423) zhke( 3) / tvols
304            WRITE (numout,9424) zhke( 4) / tvols
305            WRITE (numout,9425) zhke( 5) / tvols
306            WRITE (numout,9426) zhke( 6) / tvols
307            WRITE (numout,9427) zhke( 7) / tvols
308            WRITE (numout,9428) zhke( 8) / tvols
309            WRITE (numout,9429) zhke(10) / tvols
310            WRITE (numout,9430) zhke( 9) / tvols
311            WRITE (numout,9431)
312            WRITE (numout,9432)   &
313            &     (  zhke(1) + zhke(2) + zhke(3) + zhke(4) + zhke(5) + zhke(6)   &
314            &      + zhke(7) + zhke(8) + zhke(9) + zhke(10) ) / tvols
315         ENDIF
316
317 9420    FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================')
318 9421    FORMAT(' pressure gradient         u2= ', e20.13)
319 9422    FORMAT(' ke gradient               u2= ', e20.13)
320 9423    FORMAT(' relative vorticity term   u2= ', e20.13)
321 9424    FORMAT(' coriolis term             u2= ', e20.13)
322 9425    FORMAT(' horizontal diffusion      u2= ', e20.13)
323 9426    FORMAT(' vertical advection        u2= ', e20.13)
324 9427    FORMAT(' vertical diffusion        u2= ', e20.13)
325 9428    FORMAT(' surface pressure gradient u2= ', e20.13)
326 9429    FORMAT(' forcing term              u2= ', e20.13)
327 9430    FORMAT(' dampimg term              u2= ', e20.13)
328 9431    FORMAT(' --------------------------------------------------')
329 9432    FORMAT(' total trend               u2= ', e20.13)
330
331         IF(lwp) THEN
332            WRITE (numout,*)
333            WRITE (numout,*)
334            WRITE (numout,9440) kt
335            WRITE (numout,9441) ( zhke(2) + zhke(3) + zhke(6) ) / tvols
336            WRITE (numout,9442) ( zhke(2) + zhke(6) ) / tvols
337            WRITE (numout,9443) ( zhke(4) ) / tvols
338            WRITE (numout,9444) ( zhke(3) ) / tvols
339            WRITE (numout,9445) ( zhke(8) ) / tvols
340            WRITE (numout,9446) ( zhke(5) ) / tvols
341            WRITE (numout,9447) ( zhke(7) ) / tvols
342            WRITE (numout,9448) ( zhke(1) ) / tvols, rpktrd / tvols
343         ENDIF
344
345 9440    FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================')
346 9441    FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13)
347 9442    FORMAT(' 0 = ke gradient + vertical advection              : ', e20.13)
348 9443    FORMAT(' 0 = coriolis term  (true if key_vorenergy or key_combined): ', e20.13)
349 9444    FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.)    : ', e20.13)
350 9445    FORMAT(' 0 = surface pressure gradient                     : ', e20.13)
351 9446    FORMAT(' 0 > horizontal diffusion                          : ', e20.13)
352 9447    FORMAT(' 0 > vertical diffusion                            : ', e20.13)
353 9448    FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop)        : ', e20.13, '  u.dz(rhop) =', e20.13)
354
355         ! Save potential to kinetic energy conversion for next time step
356         
357         rpktrd = zpeke
358         
359      ENDIF
360
361   END SUBROUTINE trd_dyn
362
363
364   SUBROUTINE trd_dyn_init
365      !!---------------------------------------------------------------------
366      !!                  ***  ROUTINE trd_dyn_init  ***
367      !!
368      !! ** Purpose :   
369      !!
370      !! ** Method :
371      !!
372      !! History :
373      !!   9.0  !  03-09  (G. Madec)  Original code
374      !!---------------------------------------------------------------------
375      !! * Local declarations
376      INTEGER ::   ji, jj, jk        ! dummy loop indices
377
378      REAL(wp) ::   zmskt,zmsku,zmskv
379     
380      NAMELIST/namtrd/ ntrd, nctls
381      !!----------------------------------------------------------------------
382
383      ! 0. Initialization
384      ! -----------------
385
386      ! set the trends to zero
387      utrd  (:,:,:,:) = 0.e0
388      vtrd  (:,:,:,:) = 0.e0
389      tautrd(:,:,  :) = 0.e0
390
391      ! namelist namtrd : ocean momentum trend
392      REWIND( numnam )
393      READ  ( numnam, namtrd )
394
395      IF(lwp) THEN
396         WRITE(numout,*) 'trd_dyn : read namelist namtrd'
397         WRITE(numout,*) '~~~~~~~'
398         WRITE(numout,*) ' time step frequency trend       ntrd  = ', ntrd
399         WRITE(numout,*) ' '
400      ENDIF
401
402#if defined key_s_coord || defined key_partial_steps
403      ! dynamics trend diagnostics not yet implemented
404      IF(lwp) WRITE(numout,*) 'trd_dyn : dynamics trend diagnostics not yet implemented'
405      nstop = nstop + 1
406#endif
407
408      ! total volume at u-, v- and t-points:
409      tvols = 0.
410      tvolu = 0.
411      tvolv = 0.
412
413      DO jk = 1, jpk
414         DO jj = 2, jpjm1
415            DO ji = fs_2, fs_jpim1   ! vector opt.
416               zmskt =                      tmask_i(ji,jj) * tmask(ji,jj,jk)
417               zmsku = tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)
418               zmskv = tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)
419               tvols = tvols + zmskt * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
420               tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
421               tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
422            END DO
423         END DO
424      END DO
425# if defined key_mpp
426      CALL mpp_sum( tvols )
427      CALL mpp_sum( tvolu )
428      CALL mpp_sum( tvolv )
429# endif
430
431      IF(lwp) THEN
432         WRITE(numout,*)
433         WRITE(numout,*) 'trd_dyn : frequency ntrd = ', ntrd
434         WRITE(numout,*) '~~~~~~~   tvols = ', tvols
435         WRITE(numout,*) '          tvolu = ', tvolu,' tvolv = ', tvolv
436      ENDIF
437
438      ! 0.2 Initialization of potential to kinetic energy conversion
439     
440      rpktrd = 0.e0
441
442   END SUBROUTINE trd_dyn_init
443
444#   else 
445   !!----------------------------------------------------------------------
446   !!   Default option :                      NO mementum trend diagnostics
447   !!----------------------------------------------------------------------
448   LOGICAL, PUBLIC, PARAMETER ::   lk_trddyn = .FALSE.   ! momentum trend flag
449CONTAINS
450   SUBROUTINE trd_dyn( kt )        ! Empty routine
451      WRITE(*,*) kt
452   END SUBROUTINE trd_dyn
453   SUBROUTINE trd_dyn_init         ! Empty routine
454   END SUBROUTINE trd_dyn_init
455#endif
456
457   !!======================================================================
458END MODULE trddyn
Note: See TracBrowser for help on using the repository browser.