shapiro.F90 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113
  1. MODULE shapiro
  2. !!==============================================================================
  3. !! *** MODULE shapiro ***
  4. !! spatial filtering of input field
  5. !!==============================================================================
  6. !! History : ! 09-08 (S. Cailleau ) from N. Ferry
  7. ! 09-09 (C. Regnier ) corrections
  8. ! 04-10 (J.M. Molines) module and nemo standard
  9. !!----------------------------------------------------------------------
  10. !! * Modules used
  11. USE in_out_manager
  12. USE dom_oce ! ocean space and time domain
  13. USE timing ! Timing
  14. USE lbclnk
  15. IMPLICIT NONE
  16. PRIVATE
  17. PUBLIC Shapiro_1D ! use by sbcblk_core and sbcssr
  18. CONTAINS
  19. SUBROUTINE Shapiro_1D(ptabin, kiter, cd_overlap, ptabout) !GIG
  20. !!----------------------------------------------------------------------
  21. !! *** routine Shapiro_1D ***
  22. !!
  23. !! ** Purpose : Multiple application (kiter) of a shapiro filter
  24. !! on ptabin to produce ptabout.
  25. !!
  26. !! ** Method :
  27. !!
  28. !! ** Action : ptabout filtered output from ptabin
  29. !!
  30. !!----------------------------------------------------------------------
  31. INTEGER, INTENT(IN) :: kiter ! number of iterations to perform
  32. REAL(wp), DIMENSION(:,:), INTENT(IN) :: ptabin ! input array
  33. CHARACTER(len=*), INTENT(IN) :: cd_overlap ! = one of MERCA_GLOB, REGULAR_GLOB, ORCA_GLOB (??)
  34. REAL(wp), DIMENSION(:,:), INTENT(OUT) :: ptabout ! output array
  35. ! * Local variable
  36. INTEGER :: ji, jj, jn ! dummy loop index
  37. INTEGER :: jpim1, jpjm1 ! for compatibility ...
  38. REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zvarout ! working array
  39. REAL(wp), PARAMETER :: rp_aniso_diff_XY=2.25 ! anisotrope case (???)
  40. ! Empirical value for 140 iterations
  41. ! for an anisotropic ratio of 1.5.
  42. ! (re ???)
  43. REAL(wp) :: zalphax ! weight coeficient (x direction)
  44. REAL(wp) :: zalphay ! weight coeficient (y direction)
  45. REAL(wp) :: znum ! numerator
  46. REAL(wp) :: zden ! denominator
  47. !
  48. !------------------------------------------------------------------------------
  49. !
  50. IF( ln_timing ) CALL timing_start('Shapiro')
  51. !
  52. ALLOCATE(zvarout(jpi,jpj) )
  53. ! Global ocean case
  54. IF (( cd_overlap == 'MERCA_GLOB' ) .OR. &
  55. ( cd_overlap == 'REGULAR_GLOB' ) .OR. &
  56. ( cd_overlap == 'ORCA_GLOB' )) THEN
  57. ptabout(:,:) = ptabin(:,:)
  58. zvarout(:,:) = ptabout(:,:) ! ptabout intent out ???
  59. zalphax=1./2.
  60. zalphay=1./2.
  61. ! Dx/Dy=rp_aniso_diff_XY , D_ = vitesse de diffusion
  62. ! 140 passes du fitre, Lx/Ly=1.5, le rp_aniso_diff_XY correspondant est:
  63. IF ( rp_aniso_diff_XY >= 1. ) zalphay=zalphay/rp_aniso_diff_XY
  64. IF ( rp_aniso_diff_XY < 1. ) zalphax=zalphax*rp_aniso_diff_XY
  65. jpim1=jpi - 1
  66. jpjm1=jpj - 1
  67. DO jn = 1,kiter
  68. DO jj = 2,jpjm1
  69. DO ji = 2,jpim1
  70. ! We crop on the coast
  71. znum = zvarout(ji,jj) &
  72. + 0.25*zalphax*(zvarout(ji-1,jj )-zvarout(ji,jj))*tmask(ji-1,jj ,1) &
  73. + 0.25*zalphax*(zvarout(ji+1,jj )-zvarout(ji,jj))*tmask(ji+1,jj ,1) &
  74. + 0.25*zalphay*(zvarout(ji ,jj-1)-zvarout(ji,jj))*tmask(ji ,jj-1,1) &
  75. + 0.25*zalphay*(zvarout(ji ,jj+1)-zvarout(ji,jj))*tmask(ji ,jj+1,1)
  76. ptabout(ji,jj)=znum*tmask(ji,jj,1)+ptabin(ji,jj)*(1.-tmask(ji,jj,1))
  77. ENDDO ! end loop jj
  78. ENDDO ! end loop ji
  79. !
  80. !
  81. ! Periodical condition in case of cd_overlap (global ocean)
  82. ! - on a mercator projection grid we consider that singular point at poles
  83. ! are a mean of the values at points of the previous latitude
  84. ! - on ORCA and regular grid we copy the values at points of the previous latitude
  85. IF ( cd_overlap == 'MERCAT_GLOB' ) THEN
  86. !GIG case unchecked ! JMM for sure not valid in MPP (BUG)
  87. ptabout(1,1) = SUM(ptabout(:,2)) / jpi
  88. ptabout(jpi,jpj) = SUM(ptabout(:,jpj-1)) / jpi
  89. ELSE
  90. CALL lbc_lnk('shapiro',ptabout, 'T', 1.) ! Boundary condition
  91. ENDIF
  92. zvarout(:,:) = ptabout(:,:)
  93. ENDDO ! end loop jn
  94. ENDIF
  95. DEALLOCATE(zvarout)
  96. IF( ln_timing ) CALL timing_stop('Shapiro')
  97. !
  98. END SUBROUTINE Shapiro_1D
  99. END MODULE shapiro