modupdatebasic.F90 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434
  1. !
  2. ! $Id: modupdatebasic.F90 4779 2014-09-19 14:21:37Z rblod $
  3. !
  4. ! AGRIF (Adaptive Grid Refinement In Fortran)
  5. !
  6. ! Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
  7. ! Christophe Vouland (Christophe.Vouland@imag.fr)
  8. !
  9. ! This program is free software; you can redistribute it and/or modify
  10. ! it under the terms of the GNU General Public License as published by
  11. ! the Free Software Foundation; either version 2 of the License, or
  12. ! (at your option) any later version.
  13. !
  14. ! This program is distributed in the hope that it will be useful,
  15. ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  16. ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  17. ! GNU General Public License for more details.
  18. !
  19. ! You should have received a copy of the GNU General Public License
  20. ! along with this program; if not, write to the Free Software
  21. ! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  22. !
  23. !
  24. !
  25. !> Module containing different procedures of update (copy, average, full_weighting)
  26. !! used in the #Agrif_Update module.
  27. !===================================================================================================
  28. !
  29. module Agrif_UpdateBasic
  30. !
  31. use Agrif_Types
  32. implicit none
  33. integer, dimension(:,:), allocatable :: indchildcopy
  34. integer, dimension(:,:), allocatable :: indchildaverage
  35. !
  36. contains
  37. !
  38. !===================================================================================================
  39. ! subroutine Agrif_basicupdate_copy1d
  40. !
  41. !> Carries out a copy on a parent grid (vector x) from its child grid (vector y).
  42. !---------------------------------------------------------------------------------------------------
  43. subroutine Agrif_basicupdate_copy1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
  44. !---------------------------------------------------------------------------------------------------
  45. real, dimension(np), intent(out) :: x !< Coarse output data to parent
  46. real, dimension(nc), intent(in) :: y !< Fine input data from child
  47. integer, intent(in) :: np !< Length of parent array
  48. integer, intent(in) :: nc !< Length of child array
  49. real, intent(in) :: s_parent !< Parent grid position (s_root = 0)
  50. real, intent(in) :: s_child !< Child grid position (s_root = 0)
  51. real, intent(in) :: ds_parent !< Parent grid dx (ds_root = 1)
  52. real, intent(in) :: ds_child !< Child grid dx (ds_root = 1)
  53. !---------------------------------------------------------------------------------------------------
  54. integer :: i, locind_child_left, coeffraf
  55. !
  56. coeffraf = nint(ds_parent/ds_child)
  57. locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
  58. !
  59. if ( coeffraf == 1 ) then
  60. !CDIR ALTCODE
  61. x(1:np) = y(locind_child_left:locind_child_left+np-1)
  62. return
  63. endif
  64. !
  65. !CDIR ALTCODE
  66. do i = 1,np
  67. x(i) = y(locind_child_left)
  68. locind_child_left = locind_child_left + coeffraf
  69. enddo
  70. !---------------------------------------------------------------------------------------------------
  71. end subroutine Agrif_basicupdate_copy1d
  72. !===================================================================================================
  73. !
  74. !===================================================================================================
  75. ! subroutine Agrif_basicupdate_copy1d_before
  76. !
  77. !> Precomputes index for a copy on a parent grid (vector x) from its child grid (vector y).
  78. !---------------------------------------------------------------------------------------------------
  79. subroutine Agrif_basicupdate_copy1d_before ( nc2, np, nc, s_parent, s_child, ds_parent, ds_child, dir )
  80. !---------------------------------------------------------------------------------------------------
  81. integer, intent(in) :: nc2 !< Length of parent array
  82. integer, intent(in) :: np !< Length of parent array
  83. integer, intent(in) :: nc !< Length of child array
  84. real, intent(in) :: s_parent !< Parent grid position (s_root = 0)
  85. real, intent(in) :: s_child !< Child grid position (s_root = 0)
  86. real, intent(in) :: ds_parent !< Parent grid dx (ds_root = 1)
  87. real, intent(in) :: ds_child !< Child grid dx (ds_root = 1)
  88. integer, intent(in) :: dir !< Direction
  89. !---------------------------------------------------------------------------------------------------
  90. integer, dimension(:,:), allocatable :: indchildcopy_tmp
  91. integer :: i, n_old, locind_child_left, coeffraf
  92. !
  93. coeffraf = nint(ds_parent/ds_child)
  94. !
  95. locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
  96. if ( .not.allocated(indchildcopy) ) then
  97. allocate(indchildcopy(np*nc2, 3))
  98. else
  99. n_old = size(indchildcopy,1)
  100. if ( n_old < np*nc2 ) then
  101. allocate( indchildcopy_tmp(1:n_old, 3))
  102. indchildcopy_tmp = indchildcopy
  103. deallocate(indchildcopy)
  104. allocate(indchildcopy(np*nc2, 3))
  105. indchildcopy(1:n_old,:) = indchildcopy_tmp
  106. deallocate(indchildcopy_tmp)
  107. endif
  108. endif
  109. !
  110. do i = 1,np
  111. indchildcopy(i,dir) = locind_child_left
  112. locind_child_left = locind_child_left + coeffraf
  113. enddo
  114. !
  115. do i = 2,nc2
  116. indchildcopy(1+(i-1)*np:i*np,dir) = indchildcopy(1+(i-2)*np:(i-1)*np,dir) + nc
  117. enddo
  118. !---------------------------------------------------------------------------------------------------
  119. end subroutine Agrif_basicupdate_copy1d_before
  120. !===================================================================================================
  121. !
  122. !===================================================================================================
  123. ! subroutine Agrif_basicupdate_copy1d_after
  124. !
  125. !> Carries out a copy on a parent grid (vector x) from its child grid (vector y)
  126. !! using precomputed index.
  127. !---------------------------------------------------------------------------------------------------
  128. subroutine Agrif_basicupdate_copy1d_after ( x, y, np, nc, dir )
  129. !---------------------------------------------------------------------------------------------------
  130. real, dimension(np), intent(out) :: x !< Coarse output data to parent
  131. real, dimension(nc), intent(in) :: y !< Fine input data from child
  132. integer, intent(in) :: np !< Length of parent array
  133. integer, intent(in) :: nc !< Length of child array
  134. integer, intent(in) :: dir !< Direction
  135. !---------------------------------------------------------------------------------------------------
  136. integer :: i
  137. !
  138. !CDIR ALTCODE
  139. do i = 1,np
  140. x(i) = y(indchildcopy(i,dir))
  141. enddo
  142. !---------------------------------------------------------------------------------------------------
  143. end subroutine Agrif_basicupdate_copy1d_after
  144. !===================================================================================================
  145. !
  146. !===================================================================================================
  147. ! subroutine Agrif_basicupdate_average1d
  148. !
  149. !> Carries out an update by average on a parent grid (vector x)from its child grid (vector y).
  150. !---------------------------------------------------------------------------------------------------
  151. subroutine Agrif_basicupdate_average1d ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
  152. !---------------------------------------------------------------------------------------------------
  153. REAL, DIMENSION(np), intent(out) :: x
  154. REAL, DIMENSION(nc), intent(in) :: y
  155. INTEGER, intent(in) :: np,nc
  156. REAL, intent(in) :: s_parent, s_child
  157. REAL, intent(in) :: ds_parent, ds_child
  158. !
  159. INTEGER :: i, ii, locind_child_left, coeffraf
  160. REAL :: xpos, invcoeffraf
  161. INTEGER :: nbnonnuls
  162. INTEGER :: diffmod
  163. !
  164. coeffraf = nint(ds_parent/ds_child)
  165. invcoeffraf = 1./coeffraf
  166. !
  167. if (coeffraf == 1) then
  168. locind_child_left = 1 + nint((s_parent - s_child)/ds_child)
  169. x(1:np) = y(locind_child_left:locind_child_left+np-1)
  170. return
  171. endif
  172. !
  173. xpos = s_parent
  174. x = 0.
  175. !
  176. diffmod = 0
  177. !
  178. IF ( mod(coeffraf,2) == 0 ) diffmod = 1
  179. !
  180. locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child)
  181. !
  182. IF (Agrif_UseSpecialValueInUpdate) THEN
  183. do i = 1,np
  184. nbnonnuls = 0
  185. !CDIR NOVECTOR
  186. do ii = -coeffraf/2+locind_child_left+diffmod, &
  187. coeffraf/2+locind_child_left
  188. IF (y(ii) /= Agrif_SpecialValueFineGrid) THEN
  189. ! nbnonnuls = nbnonnuls + 1
  190. x(i) = x(i) + y(ii)
  191. ENDIF
  192. enddo
  193. ! IF (nbnonnuls /= 0) THEN
  194. ! x(i) = x(i)/nbnonnuls
  195. ! ELSE
  196. ! x(i) = Agrif_SpecialValueFineGrid
  197. ! ENDIF
  198. locind_child_left = locind_child_left + coeffraf
  199. enddo
  200. ELSE
  201. !
  202. !CDIR ALTCODE
  203. do i = 1,np
  204. !CDIR NOVECTOR
  205. do ii = -coeffraf/2+locind_child_left+diffmod, &
  206. coeffraf/2+locind_child_left
  207. x(i) = x(i) + y(ii)
  208. enddo
  209. ! x(i) = x(i)*invcoeffraf
  210. locind_child_left = locind_child_left + coeffraf
  211. enddo
  212. IF (.not.Agrif_Update_Weights) THEN
  213. x = x * invcoeffraf
  214. ENDIF
  215. ENDIF
  216. !---------------------------------------------------------------------------------------------------
  217. end subroutine Agrif_basicupdate_average1d
  218. !===================================================================================================
  219. !
  220. !===================================================================================================
  221. ! subroutine Average1dPrecompute
  222. !
  223. !> Carries out an update by average on a parent grid (vector x)from its child grid (vector y).
  224. !---------------------------------------------------------------------------------------------------
  225. subroutine Average1dPrecompute ( nc2, np, nc, s_parent, s_child, ds_parent, ds_child, dir )
  226. !---------------------------------------------------------------------------------------------------
  227. INTEGER, intent(in) :: nc2, np, nc
  228. REAL, intent(in) :: s_parent, s_child
  229. REAL, intent(in) :: ds_parent, ds_child
  230. INTEGER, intent(in) :: dir
  231. !
  232. INTEGER, DIMENSION(:,:), ALLOCATABLE :: indchildaverage_tmp
  233. INTEGER :: i, locind_child_left, coeffraf
  234. REAL :: xpos
  235. INTEGER :: diffmod
  236. !
  237. coeffraf = nint(ds_parent/ds_child)
  238. xpos = s_parent
  239. diffmod = 0
  240. !
  241. IF ( mod(coeffraf,2) == 0 ) diffmod = 1
  242. !
  243. locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child)
  244. !
  245. if (.not.allocated(indchildaverage)) then
  246. allocate(indchildaverage(np*nc2,3))
  247. else
  248. if (size(indchildaverage,1)<np*nc2) then
  249. allocate( indchildaverage_tmp(size(indchildaverage,1),size(indchildaverage,2)))
  250. indchildaverage_tmp = indchildaverage
  251. deallocate(indchildaverage)
  252. allocate(indchildaverage(np*nc2,3))
  253. indchildaverage(1:size(indchildaverage_tmp,1),1:size(indchildaverage_tmp,2)) = indchildaverage_tmp
  254. deallocate(indchildaverage_tmp)
  255. endif
  256. endif
  257. !
  258. do i = 1,np
  259. indchildaverage(i,dir)= -coeffraf/2+locind_child_left+diffmod
  260. locind_child_left = locind_child_left + coeffraf
  261. enddo
  262. !
  263. do i = 2,nc2
  264. indchildaverage(1+(i-1)*np:i*np,dir) = indchildaverage(1+(i-2)*np:(i-1)*np,dir) + nc
  265. enddo
  266. !---------------------------------------------------------------------------------------------------
  267. end subroutine Average1dPrecompute
  268. !===================================================================================================
  269. !
  270. !===================================================================================================
  271. ! subroutine Average1dAfterCompute
  272. !
  273. !> Carries out an update by average on a parent grid (vector x) from its child grid (vector y).
  274. !---------------------------------------------------------------------------------------------------
  275. subroutine Average1dAfterCompute ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child, dir )
  276. !---------------------------------------------------------------------------------------------------
  277. REAL, DIMENSION(np), intent(inout) :: x
  278. REAL, DIMENSION(nc), intent(in) :: y
  279. INTEGER, intent(in) :: np, nc
  280. REAL, intent(in) :: s_parent, s_child
  281. REAL, intent(in) :: ds_parent, ds_child
  282. INTEGER, intent(in) :: dir
  283. !
  284. REAL :: invcoeffraf
  285. INTEGER :: i, j, coeffraf
  286. INTEGER, DIMENSION(np) :: nbnonnuls
  287. REAL, DIMENSION(0:7), parameter :: invcoeff = (/1.,1.,0.5,1./3.,0.25,0.2,1./6.,1./7./)
  288. !
  289. coeffraf = nint(ds_parent/ds_child)
  290. invcoeffraf = 1./coeffraf
  291. !
  292. IF (Agrif_UseSpecialValueInUpdate) THEN
  293. !
  294. ! nbnonnuls = 0
  295. do j = 1,coeffraf
  296. do i = 1,np
  297. IF (y(indchildaverage(i,dir) + j -1) /= Agrif_SpecialValueFineGrid) THEN
  298. ! nbnonnuls(i) = nbnonnuls(i) + 1
  299. x(i) = x(i) + y(indchildaverage(i,dir) + j-1 )
  300. ENDIF
  301. enddo
  302. enddo
  303. do i=1,np
  304. ! x(i) = x(i)*invcoeff(nbnonnuls(i))
  305. ! if (nbnonnuls(i) == 0) x(i) = Agrif_SpecialValueFineGrid
  306. enddo
  307. !
  308. ELSE
  309. !
  310. !CDIR NOLOOPCHG
  311. do j = 1,coeffraf
  312. !CDIR VECTOR
  313. do i= 1,np
  314. x(i) = x(i) + y(indchildaverage(i,dir) + j-1 )
  315. enddo
  316. enddo
  317. IF (.not.Agrif_Update_Weights) THEN
  318. x = x * invcoeffraf
  319. ENDIF
  320. !
  321. ENDIF
  322. !---------------------------------------------------------------------------------------------------
  323. end subroutine Average1dAfterCompute
  324. !===================================================================================================
  325. !
  326. !===================================================================================================
  327. ! subroutine Agrif_basicupdate_full_weighting1D
  328. !
  329. !> Carries out an update by full_weighting on a parent grid (vector x) from its child grid (vector y).
  330. !---------------------------------------------------------------------------------------------------
  331. subroutine Agrif_basicupdate_full_weighting1D ( x, y, np, nc, s_parent, s_child, ds_parent, ds_child )
  332. !---------------------------------------------------------------------------------------------------
  333. real, dimension(np), intent(out) :: x
  334. real, dimension(nc), intent(in) :: y
  335. integer, intent(in) :: np, nc
  336. real, intent(in) :: s_parent, s_child
  337. real, intent(in) :: ds_parent, ds_child
  338. !---------------------------------------------------------------------------------------------------
  339. REAL :: xpos, xposfin
  340. INTEGER :: i, ii, diffmod
  341. INTEGER :: it1, it2
  342. INTEGER :: i1, i2
  343. INTEGER :: coeffraf
  344. INTEGER :: locind_child_left
  345. REAL :: sumweight, invsumweight
  346. REAL :: weights(-Agrif_MaxRaff:Agrif_MaxRaff)
  347. coeffraf = nint(ds_parent/ds_child)
  348. locind_child_left = 1 + agrif_int((s_parent-s_child)/ds_child)
  349. !
  350. if (coeffraf == 1) then
  351. x(1:np) = y(locind_child_left:locind_child_left+np-1)
  352. return
  353. endif
  354. !
  355. xpos = s_parent
  356. x = 0.
  357. !
  358. xposfin = s_child + ds_child * (locind_child_left - 1)
  359. IF (abs(xposfin - xpos) < 0.001) THEN
  360. diffmod = 0
  361. ELSE
  362. diffmod = 1
  363. ENDIF
  364. !
  365. if (diffmod == 1) THEN
  366. invsumweight=1./(2.*coeffraf**2)
  367. do i = -coeffraf,-1
  368. weights(i) = invsumweight*(2*(coeffraf+i)+1)
  369. enddo
  370. do i = 0,coeffraf-1
  371. weights(i) = weights(-(i+1))
  372. enddo
  373. it1 = -coeffraf
  374. i1 = -(coeffraf-1)+locind_child_left
  375. i2 = 2*coeffraf - 1
  376. else
  377. invsumweight=1./coeffraf**2
  378. do i = -(coeffraf-1),0
  379. weights(i) = invsumweight*(coeffraf + i)
  380. enddo
  381. do i=1,coeffraf-1
  382. weights(i) = invsumweight*(coeffraf - i)
  383. enddo
  384. it1 = -(coeffraf-1)
  385. i1 = -(coeffraf-1)+locind_child_left
  386. i2 = 2*coeffraf - 2
  387. endif
  388. !
  389. sumweight = 0.
  390. MYLOOP : do i = 1,np
  391. !
  392. it2 = it1
  393. ! sumweight = 0.
  394. do ii = i1,i1+i2
  395. !
  396. IF (Agrif_UseSpecialValueInUpdate) THEN
  397. IF (y(ii) /= Agrif_SpecialValueFineGrid) THEN
  398. x(i) = x(i) + weights(it2)*y(ii)
  399. ! sumweight = sumweight+weights(it2)
  400. ELSE
  401. x(i) = Agrif_SpecialValueFineGrid
  402. i1=i1+coeffraf
  403. CYCLE MYLOOP
  404. ENDIF
  405. ELSE
  406. x(i) = x(i) + weights(it2)*y(ii)
  407. ENDIF
  408. it2 = it2+1
  409. !
  410. enddo
  411. !
  412. i1 = i1 + coeffraf
  413. !
  414. enddo MYLOOP
  415. IF (Agrif_UseSpecialValueInUpdate) THEN
  416. x = x * coeffraf ! x will be divided by coeffraf later in modupdate.F90
  417. ENDIF
  418. !---------------------------------------------------------------------------------------------------
  419. end subroutine Agrif_basicupdate_full_weighting1D
  420. !===================================================================================================
  421. !
  422. end module Agrif_UpdateBasic