Changeset 2839 for LMDZ5/branches/testing/libf/misc
- Timestamp:
- Mar 30, 2017, 4:16:38 PM (8 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 3 deleted
- 2 edited
- 2 copied
-
. (modified) (1 prop)
-
libf/misc/regr1_conserv_m.F90 (deleted)
-
libf/misc/regr1_lint_m.F90 (deleted)
-
libf/misc/regr3_lint_m.F90 (deleted)
-
libf/misc/regr_conserv_m.F90 (copied) (copied from LMDZ5/trunk/libf/misc/regr_conserv_m.F90)
-
libf/misc/regr_lint_m.F90 (copied) (copied from LMDZ5/trunk/libf/misc/regr_lint_m.F90)
-
libf/misc/slopes_m.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 2786,2788-2790,2792-2814,2816-2838
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/misc/slopes_m.F90
r2471 r2839 1 moduleslopes_m1 MODULE slopes_m 2 2 3 3 ! Author: Lionel GUEZ 4 5 implicit none 6 7 interface slopes 8 ! This generic function computes second order slopes with Van 9 ! Leer slope-limiting, given cell averages. Reference: Dukowicz, 10 ! 1987, SIAM Journal on Scientific and Statistical Computing, 8, 11 ! 305. 12 13 ! The only difference between the specific functions is the rank 14 ! of the first argument and the equal rank of the result. 15 16 ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1 17 ! real, intent(in):: x(:) ! (n + 1) cell edges 18 ! real slopes, same shape as f ! (n, ...) 19 20 module procedure slopes1, slopes2, slopes3, slopes4 21 end interface 22 23 private 24 public slopes 25 26 contains 27 28 pure function slopes1(f, x) 29 30 real, intent(in):: f(:) 31 real, intent(in):: x(:) 32 real slopes1(size(f)) 33 34 ! Local: 35 integer n, i 36 real xc(size(f)) ! (n) cell centers 37 real h 38 39 !------------------------------------------------------ 40 41 n = size(f) 42 forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. 43 slopes1(1) = 0. 44 slopes1(n) = 0. 45 46 do i = 2, n - 1 47 if (f(i) >= max(f(i - 1), f(i + 1)) .or. f(i) & 48 <= min(f(i - 1), f(i + 1))) then 49 ! Local extremum 50 slopes1(i) = 0. 51 else 52 ! (f(i - 1), f(i), f(i + 1)) strictly monotonous 53 54 ! Second order slope: 55 slopes1(i) = (f(i + 1) - f(i - 1)) / (xc(i + 1) - xc(i - 1)) 56 57 ! Slope limitation: 58 h = abs(x(i + 1) - xc(i)) 59 slopes1(i) = sign(min(abs(slopes1(i)), abs(f(i + 1) - f(i)) / h, & 60 abs(f(i) - f(i - 1)) / h), slopes1(i)) 61 end if 62 end do 63 64 end function slopes1 65 66 !************************************************************* 67 68 pure function slopes2(f, x) 69 70 real, intent(in):: f(:, :) 71 real, intent(in):: x(:) 72 real slopes2(size(f, 1), size(f, 2)) 73 74 ! Local: 75 integer n, i, j 76 real xc(size(f, 1)) ! (n) cell centers 77 real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) 78 79 !------------------------------------------------------ 80 81 n = size(f, 1) 82 forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. 83 84 forall (i = 2:n - 1) 85 h(i) = abs(x(i + 1) - xc(i)) 86 delta_xc(i) = xc(i + 1) - xc(i - 1) 87 end forall 88 89 do j = 1, size(f, 2) 90 slopes2(1, j) = 0. 91 92 do i = 2, n - 1 93 if (f(i, j) >= max(f(i - 1, j), f(i + 1, j)) .or. & 94 f(i, j) <= min(f(i - 1, j), f(i + 1, j))) then 95 ! Local extremum 96 slopes2(i, j) = 0. 97 else 98 ! (f(i - 1, j), f(i, j), f(i + 1, j)) 99 ! strictly monotonous 100 101 ! Second order slope: 102 slopes2(i, j) = (f(i + 1, j) - f(i - 1, j)) / delta_xc(i) 103 104 ! Slope limitation: 105 slopes2(i, j) = sign(min(abs(slopes2(i, j)), & 106 abs(f(i + 1, j) - f(i, j)) / h(i), & 107 abs(f(i, j) - f(i - 1, j)) / h(i)), slopes2(i, j)) 108 end if 109 end do 110 111 slopes2(n, j) = 0. 112 end do 113 114 end function slopes2 115 116 !************************************************************* 117 118 pure function slopes3(f, x) 119 120 real, intent(in):: f(:, :, :) 121 real, intent(in):: x(:) 122 real slopes3(size(f, 1), size(f, 2), size(f, 3)) 123 124 ! Local: 125 integer n, i, j, k 126 real xc(size(f, 1)) ! (n) cell centers 127 real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) 128 129 !------------------------------------------------------ 130 131 n = size(f, 1) 132 forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. 133 134 forall (i = 2:n - 1) 135 h(i) = abs(x(i + 1) - xc(i)) 136 delta_xc(i) = xc(i + 1) - xc(i - 1) 137 end forall 138 139 do k = 1, size(f, 3) 140 do j = 1, size(f, 2) 141 slopes3(1, j, k) = 0. 142 143 do i = 2, n - 1 144 if (f(i, j, k) >= max(f(i - 1, j, k), f(i + 1, j, k)) .or. & 145 f(i, j, k) <= min(f(i - 1, j, k), f(i + 1, j, k))) then 146 ! Local extremum 147 slopes3(i, j, k) = 0. 148 else 149 ! (f(i - 1, j, k), f(i, j, k), f(i + 1, j, k)) 150 ! strictly monotonous 151 152 ! Second order slope: 153 slopes3(i, j, k) = (f(i + 1, j, k) - f(i - 1, j, k)) & 154 / delta_xc(i) 155 156 ! Slope limitation: 157 slopes3(i, j, k) = sign(min(abs(slopes3(i, j, k)), & 158 abs(f(i + 1, j, k) - f(i, j, k)) / h(i), & 159 abs(f(i, j, k) - f(i - 1, j, k)) / h(i)), slopes3(i, j, k)) 160 end if 161 end do 162 163 slopes3(n, j, k) = 0. 164 end do 165 end do 166 167 end function slopes3 168 169 !************************************************************* 170 171 pure function slopes4(f, x) 172 173 real, intent(in):: f(:, :, :, :) 174 real, intent(in):: x(:) 175 real slopes4(size(f, 1), size(f, 2), size(f, 3), size(f, 4)) 176 177 ! Local: 178 integer n, i, j, k, l 179 real xc(size(f, 1)) ! (n) cell centers 180 real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) 181 182 !------------------------------------------------------ 183 184 n = size(f, 1) 185 forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. 186 187 forall (i = 2:n - 1) 188 h(i) = abs(x(i + 1) - xc(i)) 189 delta_xc(i) = xc(i + 1) - xc(i - 1) 190 end forall 191 192 do l = 1, size(f, 4) 193 do k = 1, size(f, 3) 194 do j = 1, size(f, 2) 195 slopes4(1, j, k, l) = 0. 196 197 do i = 2, n - 1 198 if (f(i, j, k, l) >= max(f(i - 1, j, k, l), f(i + 1, j, k, l)) & 199 .or. f(i, j, k, l) & 200 <= min(f(i - 1, j, k, l), f(i + 1, j, k, l))) then 201 ! Local extremum 202 slopes4(i, j, k, l) = 0. 203 else 204 ! (f(i - 1, j, k, l), f(i, j, k, l), f(i + 1, j, k, l)) 205 ! strictly monotonous 206 207 ! Second order slope: 208 slopes4(i, j, k, l) = (f(i + 1, j, k, l) & 209 - f(i - 1, j, k, l)) / delta_xc(i) 210 211 ! Slope limitation: 212 slopes4(i, j, k, l) = sign(min(abs(slopes4(i, j, k, l)), & 213 abs(f(i + 1, j, k, l) - f(i, j, k, l)) / h(i), & 214 abs(f(i, j, k, l) - f(i - 1, j, k, l)) / h(i)), & 215 slopes4(i, j, k, l)) 216 end if 217 end do 218 219 slopes4(n, j, k, l) = 0. 220 end do 221 end do 222 end do 223 224 end function slopes4 225 226 end module slopes_m 4 ! Extension / factorisation: David CUGNET 5 6 IMPLICIT NONE 7 8 ! Those generic function computes second order slopes with Van 9 ! Leer slope-limiting, given cell averages. Reference: Dukowicz, 10 ! 1987, SIAM Journal on Scientific and Statistical Computing, 8, 11 ! 305. 12 13 ! The only difference between the specific functions is the rank 14 ! of the first argument and the equal rank of the result. 15 16 ! slope(ix,...) acts on ix th dimension. 17 18 ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1 19 ! real, intent(in):: x(:) ! (n + 1) cell edges 20 ! real slopes, same shape as f ! (n, ...) 21 INTERFACE slopes 22 MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5 23 END INTERFACE 24 25 PRIVATE 26 PUBLIC :: slopes 27 28 CONTAINS 29 30 !------------------------------------------------------------------------------- 31 ! 32 PURE FUNCTION slopes1(ix, f, x) 33 ! 34 !------------------------------------------------------------------------------- 35 ! Arguments: 36 INTEGER, INTENT(IN) :: ix 37 REAL, INTENT(IN) :: f(:) 38 REAL, INTENT(IN) :: x(:) 39 REAL :: slopes1(SIZE(f,1)) 40 !------------------------------------------------------------------------------- 41 ! Local: 42 INTEGER :: n, i, j, sta(2), sto(2) 43 REAL :: xc(SIZE(f,1)) ! (n) cell centers 44 REAL :: h(2:SIZE(f,1)-1), delta_xc(2:SIZE(f,1)-1) ! (2:n-1) 45 REAL :: fm, ff, fp, dx 46 !------------------------------------------------------------------------------- 47 n=SIZE(f,ix) 48 FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. 49 FORALL(i=2:n-1) 50 h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) 51 END FORALL 52 slopes1(:)=0. 53 DO i=2,n-1 54 ff=f(i); fm=f(i-1); fp=f(i+1) 55 IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN 56 slopes1(i)=0.; CYCLE !--- Local extremum 57 !--- 2nd order slope ; (fm, ff, fp) strictly monotonous 58 slopes1(i)=(fp-fm)/delta_xc(i) 59 !--- Slope limitation 60 slopes1(i) = SIGN(MIN(ABS(slopes1(i)), & 61 ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes1(i) ) 62 END IF 63 END DO 64 65 END FUNCTION slopes1 66 ! 67 !------------------------------------------------------------------------------- 68 69 70 !------------------------------------------------------------------------------- 71 ! 72 PURE FUNCTION slopes2(ix, f, x) 73 ! 74 !------------------------------------------------------------------------------- 75 ! Arguments: 76 INTEGER, INTENT(IN) :: ix 77 REAL, INTENT(IN) :: f(:, :) 78 REAL, INTENT(IN) :: x(:) 79 REAL :: slopes2(SIZE(f,1),SIZE(f,2)) 80 !------------------------------------------------------------------------------- 81 ! Local: 82 INTEGER :: n, i, j, sta(2), sto(2) 83 REAL, ALLOCATABLE :: xc(:) ! (n) cell centers 84 REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) 85 REAL :: fm, ff, fp, dx 86 !------------------------------------------------------------------------------- 87 n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) 88 FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. 89 FORALL(i=2:n-1) 90 h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) 91 END FORALL 92 slopes2(:,:)=0. 93 sta=[1,1]; sta(ix)=2 94 sto=SHAPE(f); sto(ix)=n-1 95 DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) 96 DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) 97 ff=f(i,j) 98 SELECT CASE(ix) 99 CASE(1); fm=f(i-1,j); fp=f(i+1,j) 100 CASE(2); fm=f(i,j-1); fp=f(i,j+1) 101 END SELECT 102 IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN 103 slopes2(i,j)=0.; CYCLE !--- Local extremum 104 !--- 2nd order slope ; (fm, ff, fp) strictly monotonous 105 slopes2(i,j)=(fp-fm)/dx 106 !--- Slope limitation 107 slopes2(i,j) = SIGN(MIN(ABS(slopes2(i,j)), & 108 ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes2(i,j) ) 109 END IF 110 END DO 111 END DO 112 DEALLOCATE(xc,h,delta_xc) 113 114 END FUNCTION slopes2 115 ! 116 !------------------------------------------------------------------------------- 117 118 119 !------------------------------------------------------------------------------- 120 ! 121 PURE FUNCTION slopes3(ix, f, x) 122 ! 123 !------------------------------------------------------------------------------- 124 ! Arguments: 125 INTEGER, INTENT(IN) :: ix 126 REAL, INTENT(IN) :: f(:, :, :) 127 REAL, INTENT(IN) :: x(:) 128 REAL :: slopes3(SIZE(f,1),SIZE(f,2),SIZE(f,3)) 129 !------------------------------------------------------------------------------- 130 ! Local: 131 INTEGER :: n, i, j, k, sta(3), sto(3) 132 REAL, ALLOCATABLE :: xc(:) ! (n) cell centers 133 REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) 134 REAL :: fm, ff, fp, dx 135 !------------------------------------------------------------------------------- 136 n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) 137 FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. 138 FORALL(i=2:n-1) 139 h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) 140 END FORALL 141 slopes3(:,:,:)=0. 142 sta=[1,1,1]; sta(ix)=2 143 sto=SHAPE(f); sto(ix)=n-1 144 DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) 145 DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) 146 DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) 147 ff=f(i,j,k) 148 SELECT CASE(ix) 149 CASE(1); fm=f(i-1,j,k); fp=f(i+1,j,k) 150 CASE(2); fm=f(i,j-1,k); fp=f(i,j+1,k) 151 CASE(3); fm=f(i,j,k-1); fp=f(i,j,k+1) 152 END SELECT 153 IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN 154 slopes3(i,j,k)=0.; CYCLE !--- Local extremum 155 !--- 2nd order slope ; (fm, ff, fp) strictly monotonous 156 slopes3(i,j,k)=(fp-fm)/dx 157 !--- Slope limitation 158 slopes3(i,j,k) = SIGN(MIN(ABS(slopes3(i,j,k)), & 159 ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes3(i,j,k) ) 160 END IF 161 END DO 162 END DO 163 END DO 164 DEALLOCATE(xc,h,delta_xc) 165 166 END FUNCTION slopes3 167 ! 168 !------------------------------------------------------------------------------- 169 170 171 !------------------------------------------------------------------------------- 172 ! 173 PURE FUNCTION slopes4(ix, f, x) 174 ! 175 !------------------------------------------------------------------------------- 176 ! Arguments: 177 INTEGER, INTENT(IN) :: ix 178 REAL, INTENT(IN) :: f(:, :, :, :) 179 REAL, INTENT(IN) :: x(:) 180 REAL :: slopes4(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4)) 181 !------------------------------------------------------------------------------- 182 ! Local: 183 INTEGER :: n, i, j, k, l, m, sta(4), sto(4) 184 REAL, ALLOCATABLE :: xc(:) ! (n) cell centers 185 REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) 186 REAL :: fm, ff, fp, dx 187 !------------------------------------------------------------------------------- 188 n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) 189 FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. 190 FORALL(i=2:n-1) 191 h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) 192 END FORALL 193 slopes4(:,:,:,:)=0. 194 sta=[1,1,1,1]; sta(ix)=2 195 sto=SHAPE(f); sto(ix)=n-1 196 DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) 197 DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) 198 DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) 199 DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) 200 ff=f(i,j,k,l) 201 SELECT CASE(ix) 202 CASE(1); fm=f(i-1,j,k,l); fp=f(i+1,j,k,l) 203 CASE(2); fm=f(i,j-1,k,l); fp=f(i,j+1,k,l) 204 CASE(3); fm=f(i,j,k-1,l); fp=f(i,j,k+1,l) 205 CASE(4); fm=f(i,j,k,l-1); fp=f(i,j,k,l+1) 206 END SELECT 207 IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN 208 slopes4(i,j,k,l)=0.; CYCLE !--- Local extremum 209 !--- 2nd order slope ; (fm, ff, fp) strictly monotonous 210 slopes4(i,j,k,l)=(fp-fm)/dx 211 !--- Slope limitation 212 slopes4(i,j,k,l) = SIGN(MIN(ABS(slopes4(i,j,k,l)), & 213 ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes4(i,j,k,l) ) 214 END IF 215 END DO 216 END DO 217 END DO 218 END DO 219 DEALLOCATE(xc,h,delta_xc) 220 221 END FUNCTION slopes4 222 ! 223 !------------------------------------------------------------------------------- 224 225 226 !------------------------------------------------------------------------------- 227 ! 228 PURE FUNCTION slopes5(ix, f, x) 229 ! 230 !------------------------------------------------------------------------------- 231 ! Arguments: 232 INTEGER, INTENT(IN) :: ix 233 REAL, INTENT(IN) :: f(:, :, :, :, :) 234 REAL, INTENT(IN) :: x(:) 235 REAL :: slopes5(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4),SIZE(f,5)) 236 !------------------------------------------------------------------------------- 237 ! Local: 238 INTEGER :: n, i, j, k, l, m, sta(5), sto(5) 239 REAL, ALLOCATABLE :: xc(:) ! (n) cell centers 240 REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) 241 REAL :: fm, ff, fp, dx 242 !------------------------------------------------------------------------------- 243 n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1)) 244 FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2. 245 FORALL(i=2:n-1) 246 h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1) 247 END FORALL 248 slopes5(:,:,:,:,:)=0. 249 sta=[1,1,1,1,1]; sta(ix)=2 250 sto=SHAPE(f); sto(ix)=n-1 251 DO m=sta(5),sto(5); IF(ix==5) dx=delta_xc(m) 252 DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) 253 DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) 254 DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) 255 DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) 256 ff=f(i,j,k,l,m) 257 SELECT CASE(ix) 258 CASE(1); fm=f(i-1,j,k,l,m); fp=f(i+1,j,k,l,m) 259 CASE(2); fm=f(i,j-1,k,l,m); fp=f(i,j+1,k,l,m) 260 CASE(3); fm=f(i,j,k-1,l,m); fp=f(i,j,k+1,l,m) 261 CASE(4); fm=f(i,j,k,l-1,m); fp=f(i,j,k,l+1,m) 262 CASE(5); fm=f(i,j,k,l,m-1); fp=f(i,j,k,l,m+1) 263 END SELECT 264 IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN 265 slopes5(i,j,k,l,m)=0.; CYCLE !--- Local extremum 266 !--- 2nd order slope ; (fm, ff, fp) strictly monotonous 267 slopes5(i,j,k,l,m)=(fp-fm)/dx 268 !--- Slope limitation 269 slopes5(i,j,k,l,m) = SIGN(MIN(ABS(slopes5(i,j,k,l,m)), & 270 ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes5(i,j,k,l,m) ) 271 END IF 272 END DO 273 END DO 274 END DO 275 END DO 276 END DO 277 DEALLOCATE(xc,h,delta_xc) 278 279 END FUNCTION slopes5 280 ! 281 !------------------------------------------------------------------------------- 282 283 END MODULE slopes_m
Note: See TracChangeset
for help on using the changeset viewer.
