[4] | 1 | !> \file slope_surf.f90 |
---|
| 2 | !! File to compute the surface slope by finite difference |
---|
| 3 | !< |
---|
| 4 | |
---|
| 5 | !> SUBROUTINE: slope_surf |
---|
| 6 | !! Compute the surface slope by finite difference |
---|
| 7 | !! \author Catritz |
---|
| 8 | !! \date ... |
---|
| 9 | !! @note Used module |
---|
| 10 | !! @note - use module3D_phy |
---|
| 11 | !< |
---|
| 12 | |
---|
| 13 | subroutine slope_surf |
---|
| 14 | |
---|
| 15 | use module3D_phy |
---|
| 16 | implicit none |
---|
| 17 | |
---|
| 18 | real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx) |
---|
| 19 | inv_4dx = 1./(4.*dx) |
---|
| 20 | |
---|
| 21 | |
---|
| 22 | |
---|
| 23 | do j=1,ny ! slope along x on node > |
---|
| 24 | do i=2,nx |
---|
| 25 | sdx(i,j)=(s(i,j)-s(i-1,j))/dx |
---|
| 26 | end do |
---|
| 27 | end do |
---|
| 28 | |
---|
| 29 | do j=2,ny ! slope along y on node ^ |
---|
| 30 | do i=1,nx |
---|
| 31 | sdy(i,j)=(s(i,j)-s(i,j-1))/dy |
---|
| 32 | end do |
---|
| 33 | end do |
---|
| 34 | |
---|
| 35 | do j=2,ny-1 ! slope amplitude on node o |
---|
| 36 | do i=2,nx-1 |
---|
| 37 | slope(i,j)=(s(i+1,j)-s(i-1,j))**2 + (s(i,j+1)-s(i,j-1))**2 |
---|
| 38 | slope(i,j)=slope(i,j)**0.5 / dx /2. |
---|
| 39 | end do |
---|
| 40 | end do |
---|
| 41 | |
---|
| 42 | |
---|
| 43 | do j=2,ny ! slope along x on node ^ (average) |
---|
| 44 | do i=2,nx-1 |
---|
| 45 | sdxmy(i,j)= ((s(i+1,j)-s(i-1,j))+s(i+1,j-1)-s(i-1,j-1))*inv_4dx |
---|
| 46 | end do |
---|
| 47 | end do |
---|
| 48 | |
---|
| 49 | do j=2,ny-1 ! slope along y on node > (average) |
---|
| 50 | do i=2,nx |
---|
| 51 | sdymx(i,j)= ((s(i,j+1)-s(i,j-1))+s(i-1,j+1)-s(i-1,j-1))*inv_4dx |
---|
| 52 | enddo |
---|
| 53 | enddo |
---|
| 54 | |
---|
| 55 | slope2mx(:,:)=sdx(:,:)**2+sdymx(:,:)**2 ! slope amplitude on node > |
---|
| 56 | slope2my(:,:)=sdy(:,:)**2+sdxmy(:,:)**2 ! slope amplitude on node ^ |
---|
| 57 | |
---|
| 58 | ! conditions on the grid boundaries |
---|
| 59 | |
---|
| 60 | slope(1,:) = 0. |
---|
| 61 | slope(:,1) = 0. |
---|
| 62 | slope(nx,:) = 0. |
---|
| 63 | slope(:,ny) = 0. |
---|
| 64 | |
---|
| 65 | slope2mx(1,:) = 0. |
---|
| 66 | slope2mx(:,1) = 0. |
---|
| 67 | slope2mx(:,ny) = 0. |
---|
| 68 | |
---|
| 69 | slope2my(1,:) = 0. |
---|
| 70 | slope2my(:,1) = 0. |
---|
| 71 | slope2my(nx,:) = 0. |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | ! limitation to 2e-2 |
---|
| 75 | sdx(:,:)=max(sdx(:,:),-2.e-2) |
---|
| 76 | sdx(:,:)=min(sdx(:,:),2.e-2) |
---|
| 77 | sdxmy(:,:)=max(sdxmy(:,:),-2.e-2) |
---|
| 78 | sdxmy(:,:)=min(sdxmy(:,:),2.e-2) |
---|
| 79 | slope2mx(:,:)=min(slope2mx(:,:),4.e-4) |
---|
| 80 | |
---|
| 81 | sdy(:,:)=max(sdy(:,:),-2.e-2) |
---|
| 82 | sdy(:,:)=min(sdy(:,:),2.e-2) |
---|
| 83 | sdymx(:,:)=max(sdymx(:,:),-2.e-2) |
---|
| 84 | sdymx(:,:)=min(sdymx(:,:),2.e-2) |
---|
| 85 | slope2my(:,:)=min(slope2my(:,:),4.e-4) |
---|
| 86 | |
---|
| 87 | slope(:,:)=min(slope(:,:),2.e-2) |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | debug_3D(:,:,21)=sdx(:,:) |
---|
| 91 | debug_3D(:,:,22)=sdy(:,:) |
---|
| 92 | end subroutine slope_surf |
---|