MODULE trddyn !!====================================================================== !! *** MODULE trddyn *** !! Ocean trends : ocean momentum trends !!===================================================================== #if defined key_trddyn || defined key_esopa !!---------------------------------------------------------------------- !! 'key_trddyn' momentum trend diag. !!---------------------------------------------------------------------- !! trd_dyn : verify the basin averaged properties for momentum !! trd_dyn_init : !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE trdtra_oce ! ocean active tracer trend USE trddyn_oce ! ocean dynamics trend USE zdf_oce ! ocean vertical physics USE ldftra_oce ! ocean active tracers: lateral physics USE ldfdyn_oce ! ocean dynamics: lateral physics USE sol_oce ! solver variables USE obc_oce ! ocean lateral open boundary condition USE eosbn2 ! equation of state (eos routine) USE in_out_manager ! I/O manager USE lib_mpp ! distributed memory computing library IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC trd_dyn ! called by step.F90 PUBLIC trd_dyn_init ! called by step.F90 !! * Shared module vaiables LOGICAL, PUBLIC, PARAMETER :: lk_trddyn = .TRUE. !: momentum trend flag !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE trd_dyn( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE trd_dyn *** !! !! ** Purpose : verify the basin averaged properties of the momentum !! equation at every time step frequency ntrd. !! momentum equation: !! * non linear momentum trend (relative vorticity trend + hori- !! zontal kinetic energy gradient trend + vertical advection !! trend) : !! 1. conserve the momentum !! 2. conserve the potential relative vorticity (rotn/e3f) !! (except for the vertical advection) !! 3. conserve the potential enstrophy (except for the !! vertical advection trend) IF no cpp key activated !! or #defined key_vorcombined !! 4. conserve the horizontal kinetic energy if #defined !! key_vorenergy !! !! ** Method : !! damping trend !! horizontal pressure gradient trend !! horizontal kinetic energy gradient trend !! relative vorticity term trend !! coriolis trend !! !! History : !! ! 91-12 (G. Madec) Original code !! ! 92-06 (M. Imbard) add time step frequency !! ! 96-01 (G. Madec) Statement function for e3 !! suppression of common work arrays !! 8.5 ! 02-09 (G. Madec) F90: Free form and module !!--------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index !! * Local declarations INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: ze1e2w,zbe1ru,zbe2rv,zbtr,zth,ztz,zpeke,zcof REAL(wp) :: zmsku, zmskv REAL(wp) :: zumo(11),zvmo(11),zhke(10) REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zkepe, zkx, zky, zkz ! workspace NAMELIST/namtrd/ ntrd, nctls !!---------------------------------------------------------------------- ! 1. forcing trend and mask trends ! -------------------------------- IF( MOD(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN ! Mask surface forcing and bottom friction fluxes DO jj = 1, jpjm1 DO ji = 1, jpim1 zmsku = tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,jk) zmskv = tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) tautrd(ji,jj,1) = tautrd(ji,jj,1) * zmsku tautrd(ji,jj,2) = tautrd(ji,jj,2) * zmskv tautrd(ji,jj,3) = tautrd(ji,jj,3) * zmsku tautrd(ji,jj,4) = tautrd(ji,jj,4) * zmskv END DO END DO tautrd( : ,jpj,1:4) = 0.e0 tautrd(jpi, : ,1:4) = 0.e0 ! Mask the trends DO jn = 1,9 DO jk = 1,jpk DO jj = 1, jpjm1 DO ji = 1, jpim1 zmsku = tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,jk) zmskv = tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) utrd(ji,jj,jk,jn) = utrd(ji,jj,jk,jn) * zmsku vtrd(ji,jj,jk,jn) = vtrd(ji,jj,jk,jn) * zmskv END DO END DO END DO END DO utrd( : ,jpj,:,1:9) = 0.e0 ; vtrd( : ,jpj,:,1:9) = 0.e0 utrd(jpi, : ,:,1:9) = 0.e0 ; vtrd(jpi, : ,:,1:9) = 0.e0 ! Conversion potential energy - kinetic energy ! c a u t i o n here, trends are computed at kt+1 (now , but after the swap) zkx(:,:,:) = 0.e0 zky(:,:,:) = 0.e0 zkz(:,:,:) = 0.e0 CALL eos( tn, sn, rhd, rhop ) ! now potential and in situ densities ! Density flux at w-point DO jk = 2, jpk DO jj = 1, jpj DO ji = 1, jpi ze1e2w = 0.5 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) * tmask_i(ji,jj) zkz(ji,jj,jk) = ze1e2w / rau0 * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) END DO END DO END DO zkz (:,:, 1 ) = 0.e0 ! Density flux at u and v-points DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, jpim1 zcof = 0.5 / rau0 zbe1ru = zcof * e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) zbe2rv = zcof * e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) zkx(ji,jj,jk) = zbe1ru * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) ) zky(ji,jj,jk) = zbe2rv * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) ) END DO END DO END DO ! Density flux divergence at t-point DO jk = 1, jpkm1 DO jj = 1, jpjm1 DO ji = 1, jpim1 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) ztz = - zbtr * ( zkz(ji,jj,jk) - zkz(ji,jj,jk+1) ) zth = - zbtr * ( ( zkx(ji,jj,jk) - zkx(ji-1,jj,jk) ) & + ( zky(ji,jj,jk) - zky(ji,jj-1,jk) ) ) zkepe(ji,jj,jk) = (zth + ztz) * tmask(ji,jj,jk) * tmask_i(ji,jj) END DO END DO END DO zkepe( : , : ,jpk) = 0.e0 zkepe( : ,jpj, : ) = 0.e0 zkepe(jpi, : , : ) = 0.e0 ! 2. Basin averaged momentum trend ! -------------------------------- DO jn = 1,9 zumo(jn) = 0. zvmo(jn) = 0. DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zumo(jn) = zumo(jn) + utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) zvmo(jn) = zvmo(jn) + vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) END DO END DO END DO END DO zumo(10) = 0. zvmo(10) = 0. zumo(11) = 0. zvmo(11) = 0. DO jj = 1, jpj DO ji = 1, jpi zumo(10) = zumo(10) + tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) zvmo(10) = zvmo(10) + tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) zumo(11) = zumo(11) + tautrd(ji,jj,3) zvmo(11) = zvmo(11) + tautrd(ji,jj,4) END DO END DO ! 3. Basin averaged kinetic energy trend ! -------------------------------------- ! Field before, because after the array swap DO jn = 1,9 zhke(jn) = 0.e0 DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi zhke(jn) = zhke(jn) & & + ub(ji,jj,jk) * utrd(ji,jj,jk,jn) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & & + vb(ji,jj,jk) * vtrd(ji,jj,jk,jn) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) END DO END DO END DO END DO zhke(10) = 0.e0 DO jj = 1, jpj DO ji = 1, jpi zhke(10) = zhke(10) & & + ub(ji,jj,1) * tautrd(ji,jj,1) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) & & + vb(ji,jj,1) * tautrd(ji,jj,2) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) END DO END DO zpeke = 0.e0 DO jk = 1,jpk DO jj = 1, jpj DO ji = 1, jpi zpeke = zpeke + zkepe(ji,jj,jk) * grav * fsdept(ji,jj,jk) & & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) END DO END DO END DO IF( lk_mpp ) THEN CALL mpp_sum( zpeke ) CALL mpp_sum( zumo , 11 ) CALL mpp_sum( zvmo , 11 ) CALL mpp_sum( zhke , 10 ) ENDIF ! 4. Print ! -------- IF(lwp) THEN WRITE (numout,*) WRITE (numout,*) WRITE (numout,9400) kt WRITE (numout,9401) zumo( 1) / tvolu, zvmo( 1) / tvolv WRITE (numout,9402) zumo( 2) / tvolu, zvmo( 2) / tvolv WRITE (numout,9403) zumo( 3) / tvolu, zvmo( 3) / tvolv WRITE (numout,9404) zumo( 4) / tvolu, zvmo( 4) / tvolv WRITE (numout,9405) zumo( 5) / tvolu, zvmo( 5) / tvolv WRITE (numout,9406) zumo( 6) / tvolu, zvmo( 6) / tvolv WRITE (numout,9407) zumo( 7) / tvolu, zvmo( 7) / tvolv WRITE (numout,9408) zumo( 8) / tvolu, zvmo( 8) / tvolv WRITE (numout,9409) zumo(10) / tvolu, zvmo(10) / tvolv WRITE (numout,9410) zumo( 9) / tvolu, zvmo( 9) / tvolv WRITE (numout,9411) zumo(11) / tvolu, zvmo(11) / tvolv WRITE (numout,9412) WRITE (numout,9413) & & ( zumo(1) + zumo(2) + zumo(3) + zumo(4) + zumo(5) + zumo(6) & & + zumo(7) + zumo(8) + zumo(9) + zumo(10) + zumo(11) ) / tvolu, & & ( zvmo(1) + zvmo(2) + zvmo(3) + zvmo(4) + zvmo(5) + zvmo(6) & & + zvmo(7) + zvmo(8) + zvmo(9) + zvmo(10) + zvmo(11) ) / tvolv ENDIF 9400 FORMAT(' momentum trend at it= ', i6, ' :', /' ==============================') 9401 FORMAT(' pressure gradient u= ', e20.13, ' v= ', e20.13) 9402 FORMAT(' ke gradient u= ', e20.13, ' v= ', e20.13) 9403 FORMAT(' relative vorticity term u= ', e20.13, ' v= ', e20.13) 9404 FORMAT(' coriolis term u= ', e20.13, ' v= ', e20.13) 9405 FORMAT(' horizontal diffusion u= ', e20.13, ' v= ', e20.13) 9406 FORMAT(' vertical advection u= ', e20.13, ' v= ', e20.13) 9407 FORMAT(' vertical diffusion u= ', e20.13, ' v= ', e20.13) 9408 FORMAT(' surface pressure gradient u= ', e20.13, ' v= ', e20.13) 9409 FORMAT(' forcing term u= ', e20.13, ' v= ', e20.13) 9410 FORMAT(' dampimg term u= ', e20.13, ' v= ', e20.13) 9411 FORMAT(' bottom flux u= ', e20.13, ' v= ', e20.13) 9412 FORMAT(' -----------------------------------------------------------------------------') 9413 FORMAT(' total trend u= ', e20.13, ' v= ', e20.13) IF(lwp) THEN WRITE (numout,*) WRITE (numout,*) WRITE (numout,9420) kt WRITE (numout,9421) zhke( 1) / tvols WRITE (numout,9422) zhke( 2) / tvols WRITE (numout,9423) zhke( 3) / tvols WRITE (numout,9424) zhke( 4) / tvols WRITE (numout,9425) zhke( 5) / tvols WRITE (numout,9426) zhke( 6) / tvols WRITE (numout,9427) zhke( 7) / tvols WRITE (numout,9428) zhke( 8) / tvols WRITE (numout,9429) zhke(10) / tvols WRITE (numout,9430) zhke( 9) / tvols WRITE (numout,9431) WRITE (numout,9432) & & ( zhke(1) + zhke(2) + zhke(3) + zhke(4) + zhke(5) + zhke(6) & & + zhke(7) + zhke(8) + zhke(9) + zhke(10) ) / tvols ENDIF 9420 FORMAT(' kinetic energy trend at it= ', i6, ' :', /' ====================================') 9421 FORMAT(' pressure gradient u2= ', e20.13) 9422 FORMAT(' ke gradient u2= ', e20.13) 9423 FORMAT(' relative vorticity term u2= ', e20.13) 9424 FORMAT(' coriolis term u2= ', e20.13) 9425 FORMAT(' horizontal diffusion u2= ', e20.13) 9426 FORMAT(' vertical advection u2= ', e20.13) 9427 FORMAT(' vertical diffusion u2= ', e20.13) 9428 FORMAT(' surface pressure gradient u2= ', e20.13) 9429 FORMAT(' forcing term u2= ', e20.13) 9430 FORMAT(' dampimg term u2= ', e20.13) 9431 FORMAT(' --------------------------------------------------') 9432 FORMAT(' total trend u2= ', e20.13) IF(lwp) THEN WRITE (numout,*) WRITE (numout,*) WRITE (numout,9440) kt WRITE (numout,9441) ( zhke(2) + zhke(3) + zhke(6) ) / tvols WRITE (numout,9442) ( zhke(2) + zhke(6) ) / tvols WRITE (numout,9443) ( zhke(4) ) / tvols WRITE (numout,9444) ( zhke(3) ) / tvols WRITE (numout,9445) ( zhke(8) ) / tvols WRITE (numout,9446) ( zhke(5) ) / tvols WRITE (numout,9447) ( zhke(7) ) / tvols WRITE (numout,9448) ( zhke(1) ) / tvols, rpktrd / tvols ENDIF 9440 FORMAT(' energetic consistency at it= ', i6, ' :', /' =========================================') 9441 FORMAT(' 0 = non linear term(true if key_vorenergy or key_combined): ', e20.13) 9442 FORMAT(' 0 = ke gradient + vertical advection : ', e20.13) 9443 FORMAT(' 0 = coriolis term (true if key_vorenergy or key_combined): ', e20.13) 9444 FORMAT(' 0 = uh.( rot(u) x uh ) (true if enstrophy conser.) : ', e20.13) 9445 FORMAT(' 0 = surface pressure gradient : ', e20.13) 9446 FORMAT(' 0 > horizontal diffusion : ', e20.13) 9447 FORMAT(' 0 > vertical diffusion : ', e20.13) 9448 FORMAT(' pressure gradient u2 = - 1/rau0 u.dz(rhop) : ', e20.13, ' u.dz(rhop) =', e20.13) ! Save potential to kinetic energy conversion for next time step rpktrd = zpeke ENDIF END SUBROUTINE trd_dyn SUBROUTINE trd_dyn_init !!--------------------------------------------------------------------- !! *** ROUTINE trd_dyn_init *** !! !! ** Purpose : !! !! ** Method : !! !! History : !! 9.0 ! 03-09 (G. Madec) Original code !!--------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zmskt,zmsku,zmskv NAMELIST/namtrd/ ntrd, nctls !!---------------------------------------------------------------------- ! 0. Initialization ! ----------------- ! set the trends to zero utrd (:,:,:,:) = 0.e0 vtrd (:,:,:,:) = 0.e0 tautrd(:,:, :) = 0.e0 ! namelist namtrd : ocean momentum trend REWIND( numnam ) READ ( numnam, namtrd ) IF(lwp) THEN WRITE(numout,*) 'trd_dyn : read namelist namtrd' WRITE(numout,*) '~~~~~~~' WRITE(numout,*) ' time step frequency trend ntrd = ', ntrd WRITE(numout,*) ' ' ENDIF #if defined key_s_coord || defined key_partial_steps ! dynamics trend diagnostics not yet implemented IF(lwp) WRITE(numout,*) 'trd_dyn : dynamics trend diagnostics not yet implemented' nstop = nstop + 1 #endif ! total volume at u-, v- and t-points: tvols = 0. tvolu = 0. tvolv = 0. DO jk = 1, jpk DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zmskt = tmask_i(ji,jj) * tmask(ji,jj,jk) zmsku = tmask_i(ji+1,jj ) * tmask_i(ji,jj) * umask(ji,jj,jk) zmskv = tmask_i(ji ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) tvols = tvols + zmskt * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) tvolu = tvolu + zmsku * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) tvolv = tvolv + zmskv * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) END DO END DO END DO IF( lk_mpp ) CALL mpp_sum( tvols ) ! sums over the global domain IF( lk_mpp ) CALL mpp_sum( tvolu ) IF( lk_mpp ) CALL mpp_sum( tvolv ) IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'trd_dyn : frequency ntrd = ', ntrd WRITE(numout,*) '~~~~~~~ tvols = ', tvols WRITE(numout,*) ' tvolu = ', tvolu,' tvolv = ', tvolv ENDIF ! 0.2 Initialization of potential to kinetic energy conversion rpktrd = 0.e0 END SUBROUTINE trd_dyn_init # else !!---------------------------------------------------------------------- !! Default option : NO mementum trend diagnostics !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_trddyn = .FALSE. !: momentum trend flag CONTAINS SUBROUTINE trd_dyn( kt ) ! Empty routine WRITE(*,*) 'trd_dyn: You should not have seen this print! error?', kt END SUBROUTINE trd_dyn SUBROUTINE trd_dyn_init ! Empty routine END SUBROUTINE trd_dyn_init #endif !!====================================================================== END MODULE trddyn