dynkeg.F90 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162
  1. MODULE dynkeg
  2. !!======================================================================
  3. !! *** MODULE dynkeg ***
  4. !! Ocean dynamics: kinetic energy gradient trend
  5. !!======================================================================
  6. !! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
  7. !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
  8. !! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module
  9. !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
  10. !!----------------------------------------------------------------------
  11. !!----------------------------------------------------------------------
  12. !! dyn_keg : update the momentum trend with the horizontal tke
  13. !!----------------------------------------------------------------------
  14. USE oce ! ocean dynamics and tracers
  15. USE dom_oce ! ocean space and time domain
  16. USE trd_oce ! trends: ocean variables
  17. USE trddyn ! trend manager: dynamics
  18. !
  19. USE in_out_manager ! I/O manager
  20. USE lbclnk ! ocean lateral boundary conditions (or mpp link)
  21. USE lib_mpp ! MPP library
  22. USE prtctl ! Print control
  23. USE wrk_nemo ! Memory Allocation
  24. USE timing ! Timing
  25. IMPLICIT NONE
  26. PRIVATE
  27. PUBLIC dyn_keg ! routine called by step module
  28. INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme)
  29. INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
  30. !
  31. REAL(wp) :: r1_48 = 1._wp / 48._wp !: =1/(4*2*6)
  32. !! * Substitutions
  33. # include "vectopt_loop_substitute.h90"
  34. !!----------------------------------------------------------------------
  35. !! NEMO/OPA 3.6 , NEMO Consortium (2015)
  36. !! $Id: dynkeg.F90 4990 2014-12-15 16:42:49Z timgraham $
  37. !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
  38. !!----------------------------------------------------------------------
  39. CONTAINS
  40. SUBROUTINE dyn_keg( kt, kscheme )
  41. !!----------------------------------------------------------------------
  42. !! *** ROUTINE dyn_keg ***
  43. !!
  44. !! ** Purpose : Compute the now momentum trend due to the horizontal
  45. !! gradient of the horizontal kinetic energy and add it to the
  46. !! general momentum trend.
  47. !!
  48. !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that
  49. !! conserve kinetic energy. Compute the now horizontal kinetic energy
  50. !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
  51. !! * kscheme = nkeg_HW : Hollingsworth correction following
  52. !! Arakawa (2001). The now horizontal kinetic energy is given by:
  53. !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 )
  54. !! + mj-1( 2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2 ) ]
  55. !!
  56. !! Take its horizontal gradient and add it to the general momentum
  57. !! trend (ua,va).
  58. !! ua = ua - 1/e1u di[ zhke ]
  59. !! va = va - 1/e2v dj[ zhke ]
  60. !!
  61. !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
  62. !! - send this trends to trd_dyn (l_trddyn=T) for post-processing
  63. !!
  64. !! ** References : Arakawa, A., International Geophysics 2001.
  65. !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
  66. !!----------------------------------------------------------------------
  67. INTEGER, INTENT( in ) :: kt ! ocean time-step index
  68. INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme
  69. !
  70. INTEGER :: ji, jj, jk ! dummy loop indices
  71. REAL(wp) :: zu, zv ! temporary scalars
  72. REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
  73. REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
  74. !!----------------------------------------------------------------------
  75. !
  76. IF( nn_timing == 1 ) CALL timing_start('dyn_keg')
  77. !
  78. CALL wrk_alloc( jpi,jpj,jpk, zhke )
  79. !
  80. IF( kt == nit000 ) THEN
  81. IF(lwp) WRITE(numout,*)
  82. IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
  83. IF(lwp) WRITE(numout,*) '~~~~~~~'
  84. ENDIF
  85. IF( l_trddyn ) THEN ! Save ua and va trends
  86. CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
  87. ztrdu(:,:,:) = ua(:,:,:)
  88. ztrdv(:,:,:) = va(:,:,:)
  89. ENDIF
  90. zhke(:,:,jpk) = 0._wp
  91. SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==!
  92. !
  93. CASE ( nkeg_C2 ) !-- Standard scheme --!
  94. DO jk = 1, jpkm1
  95. DO jj = 2, jpj
  96. DO ji = fs_2, jpi ! vector opt.
  97. zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) &
  98. & + un(ji ,jj ,jk) * un(ji ,jj ,jk)
  99. zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &
  100. & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk)
  101. zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
  102. END DO
  103. END DO
  104. END DO
  105. !
  106. CASE ( nkeg_HW ) !-- Hollingsworth scheme --!
  107. DO jk = 1, jpkm1
  108. DO jj = 2, jpjm1
  109. DO ji = fs_2, jpim1 ! vector opt.
  110. zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) &
  111. & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) &
  112. & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) &
  113. & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) )
  114. !
  115. zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &
  116. & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) &
  117. & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) &
  118. & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) )
  119. zhke(ji,jj,jk) = r1_48 * ( zv + zu )
  120. END DO
  121. END DO
  122. END DO
  123. CALL lbc_lnk( zhke, 'T', 1. )
  124. !
  125. END SELECT
  126. !
  127. DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==!
  128. DO jj = 2, jpjm1
  129. DO ji = fs_2, fs_jpim1 ! vector opt.
  130. ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
  131. va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
  132. END DO
  133. END DO
  134. END DO
  135. !
  136. IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic
  137. ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
  138. ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
  139. CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
  140. CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
  141. ENDIF
  142. !
  143. IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' keg - Ua: ', mask1=umask, &
  144. & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
  145. !
  146. CALL wrk_dealloc( jpi,jpj,jpk, zhke )
  147. !
  148. IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')
  149. !
  150. END SUBROUTINE dyn_keg
  151. !!======================================================================
  152. END MODULE dynkeg