[4] | 1 | subroutine compact_ddx(data,deriv,flag,end_values,m1,m2,m3) |
---|
| 2 | ! Routine to compute the x derivative (1st dimension) |
---|
| 3 | ! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming |
---|
| 4 | ! a uniformly spaced grid mesh with dx=h. |
---|
| 5 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 6 | ! inputs |
---|
| 7 | ! data m1,m2,m3 array of input values |
---|
| 8 | ! m1,m2,m3 dimensions of field to be differentiated |
---|
| 9 | ! h x grid spacing |
---|
| 10 | ! flag flag which identifies periodicity |
---|
| 11 | ! end_values the end values for the derivatives |
---|
| 12 | ! (required but now obsolete) |
---|
| 13 | ! |
---|
| 14 | ! output |
---|
| 15 | ! deriv m1,m2,m3 array containing the optimized |
---|
| 16 | ! 6th order accurate d/dx(data) |
---|
| 17 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 18 | |
---|
| 19 | use grid_info |
---|
| 20 | use boundary_information |
---|
| 21 | use user_parameters, only: LAPACK_FLAG |
---|
| 22 | implicit none |
---|
| 23 | |
---|
| 24 | integer, intent(in) :: flag,m1,m2,m3 |
---|
| 25 | real, dimension(m1,m2,m3) :: data |
---|
| 26 | real, dimension(m1,m2,m3) :: deriv |
---|
| 27 | real, dimension(2) :: end_values |
---|
| 28 | real, dimension(:,:), allocatable :: work |
---|
| 29 | integer :: i,j,k,ierr,level,maxlevels |
---|
| 30 | real :: h |
---|
| 31 | parameter(maxlevels=12) |
---|
| 32 | |
---|
| 33 | real, allocatable :: kx(:),in(:),out(:) |
---|
| 34 | real :: pi |
---|
| 35 | integer*8 :: plan_f, plan_i !size must match pointer size, 8 for DEC alphas |
---|
| 36 | integer :: N |
---|
| 37 | include 'fftw_kw.inc' ! predefined FFTW constants |
---|
| 38 | |
---|
| 39 | if( boundary_type2(1,1)=='periodic' ) then ! spectral differentiation in x, use FFTW |
---|
| 40 | |
---|
| 41 | ! preliminaries, setup fftw "plans" |
---|
| 42 | pi=4.*atan(1.) |
---|
| 43 | N=m1-1 ! m1 includes both boundary points, N only 1 |
---|
| 44 | allocate( kx(N),in(N),out(N) ) |
---|
| 45 | |
---|
| 46 | call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE) |
---|
| 47 | call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE) |
---|
| 48 | do i=1,N/2+1 |
---|
| 49 | kx(i)=2.*pi*(i-1.) ! FFTW real to half-complex storage |
---|
| 50 | enddo |
---|
| 51 | do i=2,N/2 |
---|
| 52 | kx(N-i+2)=kx(i) ! FFTW real to half-complex storage |
---|
| 53 | enddo |
---|
| 54 | |
---|
| 55 | do j=1,m2 |
---|
| 56 | do k=1,m3 |
---|
| 57 | |
---|
| 58 | in(1:N)=data(1:N,j,k) |
---|
| 59 | call rfftw_f77_one(plan_f,in,out) ! FORWARD TRANSFORM |
---|
| 60 | |
---|
| 61 | ! multiply by wavenumber array & normalize by N |
---|
| 62 | out(1:N) = kx(1:N)*out(1:N)/float(N) |
---|
| 63 | ! multiply by i (fold and negate) |
---|
| 64 | do i=2,N/2 |
---|
| 65 | in(i)=-out(N-i+2) ! real <-- (-imaginary) |
---|
| 66 | in(N-i+2)=out(i) ! imag <-- ( real ) |
---|
| 67 | enddo |
---|
| 68 | in(1)=0. ! dc |
---|
| 69 | in(N/2+1)=0. !nyquist |
---|
| 70 | |
---|
| 71 | call rfftw_f77_one(plan_i,in,out) ! INVERSE TRANSFORM |
---|
| 72 | deriv(1:N,j,k)=out(1:N) ! store results in deriv |
---|
| 73 | deriv(m1,j,k)=out(1) ! carry 2nd periodic end point along |
---|
| 74 | |
---|
| 75 | enddo |
---|
| 76 | enddo |
---|
| 77 | |
---|
| 78 | call rfftw_f77_destroy_plan(plan_f) |
---|
| 79 | call rfftw_f77_destroy_plan(plan_i) |
---|
| 80 | deallocate( kx,in,out ) |
---|
| 81 | return |
---|
| 82 | |
---|
| 83 | endif |
---|
| 84 | |
---|
| 85 | !************* nonperiodic, use compact differentiation ******************** |
---|
| 86 | allocate( work(m1,m2) ) |
---|
| 87 | |
---|
| 88 | ! determine an appropriate grid level from m1 |
---|
| 89 | do level=0,maxlevels |
---|
| 90 | if( m1 == Grid(level)%n1 ) goto 999 |
---|
| 91 | enddo |
---|
| 92 | 999 continue ! can use ddx and piv arrays from grid level |
---|
| 93 | h = Grid(level)%dxi |
---|
| 94 | |
---|
| 95 | !loop through each xy horizontal plane |
---|
| 96 | |
---|
| 97 | do k=1,m3 |
---|
| 98 | |
---|
| 99 | !Create rhs vectors, store in work(:,:) |
---|
| 100 | work(1,:)=( -(274.)*data(1,:,k)+(600.)*data(2,:,k) & |
---|
| 101 | -(600.)*data(3,:,k)+(400.)*data(4,:,k) & |
---|
| 102 | -(150.)*data(5,:,k)+ (24.)*data(6,:,k) )/(120.*h) |
---|
| 103 | ! work(2,:)=( -(5./9.)*data(1,:,k)-(1./2.)*data(2,:,k) & |
---|
| 104 | ! +(1.)*data(3,:,k)+(1./18.)*data(4,:,k) )/h |
---|
| 105 | work(2,:)=( -(24.)*data(1,:,k)-(130.)*data(2,:,k) & |
---|
| 106 | +(240.)*data(3,:,k)-(120.)*data(4,:,k) & |
---|
| 107 | +(40.)*data(5,:,k)- (6.)*data(6,:,k) )/(120.*h) |
---|
| 108 | work(3,:)= ( (1./36.)*(data(5,:,k)-data(1,:,k)) & |
---|
| 109 | + (7./9.)*(data(4,:,k)-data(2,:,k)) )/h |
---|
| 110 | |
---|
| 111 | do i=4,m1-3 |
---|
| 112 | work(i,:)= ( (c1/6.)*(data(i+3,:,k)-data(i-3,:,k)) & |
---|
| 113 | + (b1/4.)*(data(i+2,:,k)-data(i-2,:,k)) & |
---|
| 114 | + (a1/2.)*(data(i+1,:,k)-data(i-1,:,k)) )/h |
---|
| 115 | enddo |
---|
| 116 | |
---|
| 117 | work(m1-2,:)= (-(1./36.)*(data(m1-4,:,k)-data(m1,:,k)) & |
---|
| 118 | - (7./9.)*(data(m1-3,:,k)-data(m1-1,:,k)) )/h |
---|
| 119 | ! work(m1-1,:)=( (5./9.)*data(m1,:,k)+(1./2.)*data(m1-1,:,k) & |
---|
| 120 | ! -(1)*data(m1-2,:,k)-(1./18.)*data(m1-3,:,k) )/h |
---|
| 121 | work(m1-1,:)=( (24.)*data(m1,:,k)+(130.)*data(m1-1,:,k) & |
---|
| 122 | -(240.)*data(m1-2,:,k)+(120.)*data(m1-3,:,k) & |
---|
| 123 | -(40.)*data(m1-4,:,k)+ (6.)*data(m1-5,:,k) )/(120.*h) |
---|
| 124 | work(m1,:)=( (274.)*data(m1,:,k)-(600.)*data(m1-1,:,k) & |
---|
| 125 | +(600.)*data(m1-2,:,k)-(400.)*data(m1-3,:,k) & |
---|
| 126 | +(150.)*data(m1-4,:,k)- (24.)*data(m1-5,:,k) )/(120.*h) |
---|
| 127 | |
---|
| 128 | !Compute the x derivatives for plane k, overwriting the array work |
---|
| 129 | |
---|
| 130 | if( LAPACK_FLAG == 'double' ) then |
---|
| 131 | !!! call DGBTRS('N',m1,2,2,m2,Grid(level)%ddx,7,Grid(level)%piv_dx,work(1,1),m1,ierr) |
---|
| 132 | write(0,*) 'LAPACK_FLAG = double instead of single' |
---|
| 133 | stop |
---|
| 134 | elseif( LAPACK_FLAG == 'single' ) then |
---|
| 135 | call SGBTRS('N',m1,2,2,m2,Grid(level)%ddx,7,Grid(level)%piv_dx,work(1,1),m1,ierr) |
---|
| 136 | endif |
---|
| 137 | |
---|
| 138 | !If deriv. end conditions are 'periodic' (flag = 5) then |
---|
| 139 | !overwrite the endpts with nonlocal estimate |
---|
| 140 | if( flag .eq. 5 ) then |
---|
| 141 | |
---|
| 142 | work(1,:) = ( 2.*data(m1-2,:,k) -16.*data(m1-1,:,k) & |
---|
| 143 | +16.*data(2,:,k) -2.*data(3,:,k) )/(24.*h) |
---|
| 144 | work(m1,:)=work(1,:) |
---|
| 145 | |
---|
| 146 | work(2,:) = ( 2.*data(m1-1,:,k) -16.*data(m1,:,k) & |
---|
| 147 | +16.*data(3,:,k) -2.*data(4,:,k) )/(24.*h) |
---|
| 148 | work(m1-1,:) = ( 2.*data(m1-3,:,k) -16.*data(m1-2,:,k) & |
---|
| 149 | +16.*data(1,:,k) -2.*data(2,:,k) )/(24.*h) |
---|
| 150 | |
---|
| 151 | work(3,:) = ( 2.*data(m1,:,k) -16.*data(2,:,k) & |
---|
| 152 | +16.*data(4,:,k) -2.*data(5,:,k) )/(24.*h) |
---|
| 153 | work(m1-2,:) = ( 2.*data(m1-4,:,k) -16.*data(m1-3,:,k) & |
---|
| 154 | +16.*data(m1-1,:,k) -2.*data(m1,:,k) )/(24.*h) |
---|
| 155 | |
---|
| 156 | endif |
---|
| 157 | |
---|
| 158 | !store the results in deriv(:,:,k) |
---|
| 159 | deriv(:,:,k)=work(:,:) |
---|
| 160 | |
---|
| 161 | !Go on to the next k plane |
---|
| 162 | enddo |
---|
| 163 | |
---|
| 164 | deallocate( work ) |
---|
| 165 | end subroutine compact_ddx |
---|
| 166 | |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | |
---|
| 171 | subroutine compact_ddy(data,deriv,flag,end_values,m1,m2,m3) |
---|
| 172 | ! Routine to compute the y derivative (2nd dimension) |
---|
| 173 | ! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming |
---|
| 174 | ! a uniformly spaced grid mesh with dy=h. |
---|
| 175 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 176 | ! inputs |
---|
| 177 | ! data m1,m2,m3 array of input values |
---|
| 178 | ! m1,m2,m3 dimensions of field to be differentiated |
---|
| 179 | ! h y grid spacing |
---|
| 180 | ! flag flag which identifies periodicity |
---|
| 181 | ! end_values the end values for the derivatives |
---|
| 182 | ! (required but obsolete) |
---|
| 183 | ! |
---|
| 184 | ! output |
---|
| 185 | ! deriv m1,m2,m3 array containing the optimized |
---|
| 186 | ! 6th order accurate d/dy(data) |
---|
| 187 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 188 | |
---|
| 189 | use grid_info |
---|
| 190 | use boundary_information |
---|
| 191 | use user_parameters, only: LAPACK_FLAG |
---|
| 192 | implicit none |
---|
| 193 | |
---|
| 194 | integer, intent(in) :: flag,m1,m2,m3 |
---|
| 195 | real, dimension(m1,m2,m3), intent(in) :: data |
---|
| 196 | real, dimension(m1,m2,m3), intent(out) :: deriv |
---|
| 197 | real, dimension(2), intent(in) :: end_values |
---|
| 198 | !!! double precision, dimension(:,:), allocatable :: work |
---|
| 199 | real, dimension(:,:), allocatable :: work |
---|
| 200 | integer :: i,j,k,ierr,level,maxlevels |
---|
| 201 | real :: h |
---|
| 202 | parameter(maxlevels=12) |
---|
| 203 | |
---|
| 204 | real, allocatable :: ky(:),in(:),out(:) |
---|
| 205 | real :: pi |
---|
| 206 | integer*8 :: plan_f, plan_i !size must match pointer size, 8 for DEC alphas |
---|
| 207 | integer :: N |
---|
| 208 | include 'fftw_kw.inc' ! predefined FFTW constants |
---|
| 209 | |
---|
| 210 | if( boundary_type4(1,1)=='periodic' ) then ! spectral differentiation in y, use FFTW |
---|
| 211 | |
---|
| 212 | ! preliminaries, setup fftw "plans" |
---|
| 213 | pi=4.*atan(1.) |
---|
| 214 | N=m2-1 ! m2 includes both boundary points, N only 1 |
---|
| 215 | allocate( ky(N),in(N),out(N) ) |
---|
| 216 | |
---|
| 217 | call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE) |
---|
| 218 | call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE) |
---|
| 219 | do j=1,N/2+1 |
---|
| 220 | ky(j)=2.*pi*(j-1.) ! FFTW real to half-complex storage |
---|
| 221 | enddo |
---|
| 222 | do j=2,N/2 |
---|
| 223 | ky(N-j+2)=ky(j) ! FFTW real to half-complex storage |
---|
| 224 | enddo |
---|
| 225 | |
---|
| 226 | do i=1,m1 |
---|
| 227 | do k=1,m3 |
---|
| 228 | |
---|
| 229 | in(1:N)=data(i,1:N,k) |
---|
| 230 | call rfftw_f77_one(plan_f,in,out) ! FORWARD TRANSFORM |
---|
| 231 | |
---|
| 232 | ! multiply by wavenumber array & normalize by N |
---|
| 233 | out(1:N) = ky(1:N)*out(1:N)/float(N) |
---|
| 234 | ! multiply by i (fold and negate) |
---|
| 235 | do j=2,N/2 |
---|
| 236 | in(j)=-out(N-j+2) ! real <-- (-imaginary) |
---|
| 237 | in(N-j+2)=out(j) ! imag <-- ( real ) |
---|
| 238 | enddo |
---|
| 239 | in(1)=0. ! dc |
---|
| 240 | in(N/2+1)=0. !nyquist |
---|
| 241 | |
---|
| 242 | call rfftw_f77_one(plan_i,in,out) ! INVERSE TRANSFORM |
---|
| 243 | deriv(i,1:N,k)=out(1:N) ! store results in deriv |
---|
| 244 | deriv(i,m2,k)=out(1) ! carry 2nd periodic end point along |
---|
| 245 | |
---|
| 246 | enddo |
---|
| 247 | enddo |
---|
| 248 | |
---|
| 249 | call rfftw_f77_destroy_plan(plan_f) |
---|
| 250 | call rfftw_f77_destroy_plan(plan_i) |
---|
| 251 | deallocate( ky,in,out ) |
---|
| 252 | return |
---|
| 253 | |
---|
| 254 | endif |
---|
| 255 | |
---|
| 256 | !************* nonperiodic, use compact differentiation ******************** |
---|
| 257 | allocate( work(m2,m1) ) |
---|
| 258 | |
---|
| 259 | ! determine an appropriate grid level from m2 |
---|
| 260 | do level=0,maxlevels-1 |
---|
| 261 | if( m2 == Grid(level)%n2 ) goto 999 |
---|
| 262 | enddo |
---|
| 263 | 999 continue ! can use ddx and piv arrays from grid level |
---|
| 264 | h = Grid(level)%deta |
---|
| 265 | |
---|
| 266 | !loop through each xy horizontal plane |
---|
| 267 | |
---|
| 268 | do k=1,m3 |
---|
| 269 | |
---|
| 270 | ! create rhs vectors, store in work(:,:,1) |
---|
| 271 | work(1,:)=( -(274.)*data(:,1,k)+(600.)*data(:,2,k) & |
---|
| 272 | -(600.)*data(:,3,k)+(400.)*data(:,4,k) & |
---|
| 273 | -(150.)*data(:,5,k)+ (24.)*data(:,6,k) )/(120.*h) |
---|
| 274 | ! work(2,:)=( -(5./9.)*data(:,1,k)-(1./2.)*data(:,2,k) & |
---|
| 275 | ! +(1.)*data(:,3,k)+(1./18.)*data(:,4,k) )/h |
---|
| 276 | work(2,:)=( -(24.)*data(:,1,k)-(130.)*data(:,2,k) & |
---|
| 277 | +(240.)*data(:,3,k)-(120.)*data(:,4,k) & |
---|
| 278 | +(40.)*data(:,5,k)- (6.)*data(:,6,k) )/(120.*h) |
---|
| 279 | work(3,:) = ( (1./36.)*( data(:,5,k) - data(:,1,k) ) & |
---|
| 280 | + (7./9.)*( data(:,4,k) - data(:,2,k) ) )/h |
---|
| 281 | do j=4,m2-3 |
---|
| 282 | work(j,:) = ( (c1/6.)*( data(:,j+3,k) - data(:,j-3,k) ) & |
---|
| 283 | + (b1/4.)*( data(:,j+2,k) - data(:,j-2,k) ) & |
---|
| 284 | + (a1/2.)*( data(:,j+1,k) - data(:,j-1,k) ) )/h |
---|
| 285 | enddo |
---|
| 286 | |
---|
| 287 | work(m2-2,:) = (-(1./36.)*( data(:,m2-4,k) - data(:,m2,k) ) & |
---|
| 288 | - (7./9.)*( data(:,m2-3,k) - data(:,m2-1,k) ) )/h |
---|
| 289 | ! work(m2-1,:)= ( (5./9.)*data(:,m2,k)+(1./2.)*data(:,m2-1,k) & |
---|
| 290 | ! -(1.)*data(:,m2-2,k)-(1./18.)*data(:,m2-3,k) )/h |
---|
| 291 | work(m2-1,:)=( (24.)*data(:,m2,k)+(130.)*data(:,m2-1,k) & |
---|
| 292 | -(240.)*data(:,m2-2,k)+(120.)*data(:,m2-3,k) & |
---|
| 293 | -(40.)*data(:,m2-4,k)+ (6.)*data(:,m2-5,k) )/(120.*h) |
---|
| 294 | work(m2,:)=( (274.)*data(:,m2,k)-(600.)*data(:,m2-1,k) & |
---|
| 295 | +(600.)*data(:,m2-2,k)-(400.)*data(:,m2-3,k) & |
---|
| 296 | +(150.)*data(:,m2-4,k)- (24.)*data(:,m2-5,k) )/(120.*h) |
---|
| 297 | |
---|
| 298 | !Compute the y derivatives for plane k, overwriting the array work |
---|
| 299 | if( LAPACK_FLAG == 'double' ) then |
---|
| 300 | !!! call DGBTRS('N',m2,2,2,m1,Grid(level)%ddy,7,Grid(level)%piv_dy,work(1,1),m2,ierr) |
---|
| 301 | write(0,*) 'LAPACK_FLAG = double instead of single' |
---|
| 302 | stop |
---|
| 303 | elseif( LAPACK_FLAG == 'single' ) then |
---|
| 304 | call SGBTRS('N',m2,2,2,m1,Grid(level)%ddy,7,Grid(level)%piv_dy,work(1,1),m2,ierr) |
---|
| 305 | endif |
---|
| 306 | |
---|
| 307 | !If deriv. end conditions are 'periodic' (flag = 5) then |
---|
| 308 | !overwrite the endpts with nonlocal estimate |
---|
| 309 | if( flag .eq. 5 ) then |
---|
| 310 | work(1,:) = ( 2.*data(:,m2-2,k) -16.*data(:,m2-1,k) & |
---|
| 311 | +16.*data(:,2,k) -2.*data(:,3,k) )/(24.*h) |
---|
| 312 | work(m2,:)=work(1,:) |
---|
| 313 | |
---|
| 314 | work(2,:) = ( 2.*data(:,m2-1,k) -16.*data(:,m2,k) & |
---|
| 315 | +16.*data(:,3,k) -2.*data(:,4,k) )/(24.*h) |
---|
| 316 | work(m2-1,:) = ( 2.*data(:,m2-3,k) -16.*data(:,m2-2,k) & |
---|
| 317 | +16.*data(:,1,k) -2.*data(:,2,k) )/(24.*h) |
---|
| 318 | |
---|
| 319 | work(3,:) = ( 2.*data(:,m2,k) -16.*data(:,2,k) & |
---|
| 320 | +16.*data(:,4,k) -2.*data(:,5,k) )/(24.*h) |
---|
| 321 | work(m2-2,:) = ( 2.*data(:,m2-4,k) -16.*data(:,m2-3,k) & |
---|
| 322 | +16.*data(:,m2-1,k) -2.*data(:,m2,k) )/(24.*h) |
---|
| 323 | |
---|
| 324 | endif |
---|
| 325 | |
---|
| 326 | !Store the results in deriv(:,:,k) |
---|
| 327 | deriv(1:m1,1:m2,k)=transpose( work(1:m2,1:m1) ) |
---|
| 328 | |
---|
| 329 | !Go on to the next k plane |
---|
| 330 | enddo |
---|
| 331 | |
---|
| 332 | deallocate( work ) |
---|
| 333 | end subroutine compact_ddy |
---|
| 334 | |
---|
| 335 | |
---|
| 336 | |
---|
| 337 | |
---|
| 338 | subroutine compact_ddz(data,deriv,flag,end_values,m1,m2,m3) |
---|
| 339 | ! Routine to compute the x derivative (3rd dimension) |
---|
| 340 | ! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming |
---|
| 341 | ! a uniformly spaced grid mesh with dz=h. |
---|
| 342 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 343 | ! inputs |
---|
| 344 | ! data m1,m2,m3 array of input values |
---|
| 345 | ! m1,m2,m3 dimensions of field to be differentiated |
---|
| 346 | ! h z grid spacing |
---|
| 347 | ! flag flag which identifies end point scheme |
---|
| 348 | ! end_values the end values for the derivatives |
---|
| 349 | ! (required but obsolete) |
---|
| 350 | ! |
---|
| 351 | ! output |
---|
| 352 | ! deriv m1,m2,m3 array containing the optimized |
---|
| 353 | ! 6th order accurate d/dz(data) |
---|
| 354 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 355 | |
---|
| 356 | use grid_info |
---|
| 357 | use user_parameters, only: LAPACK_FLAG |
---|
| 358 | implicit none |
---|
| 359 | |
---|
| 360 | integer, intent(in) :: flag,m1,m2,m3 |
---|
| 361 | real, dimension(m1,m2,m3), intent(in) :: data |
---|
| 362 | real, dimension(m1,m2,m3), intent(out) :: deriv |
---|
| 363 | real, dimension(2), intent(in) :: end_values |
---|
| 364 | !!! double precision, dimension(:,:), allocatable :: work |
---|
| 365 | real, dimension(:,:), allocatable :: work |
---|
| 366 | integer :: i,j,k,ierr,level,maxlevels |
---|
| 367 | real :: h |
---|
| 368 | parameter(maxlevels=12) |
---|
| 369 | |
---|
| 370 | allocate( work(m3,m1) ) |
---|
| 371 | |
---|
| 372 | ! determine an appropriate grid level from m3 |
---|
| 373 | do level=0,maxlevels-1 |
---|
| 374 | if( m3 == Grid(level)%n3 ) goto 999 |
---|
| 375 | enddo |
---|
| 376 | 999 continue ! can use ddx and piv arrays from grid level |
---|
| 377 | h = Grid(level)%dzeta |
---|
| 378 | |
---|
| 379 | !loop through each xz plane |
---|
| 380 | |
---|
| 381 | do j=1,m2 |
---|
| 382 | |
---|
| 383 | !create rhs vectors, store in work(:,:,1) |
---|
| 384 | work(1,:)=( -(274.)*data(:,j,1)+(600.)*data(:,j,2) & |
---|
| 385 | -(600.)*data(:,j,3)+(400.)*data(:,j,4) & |
---|
| 386 | -(150.)*data(:,j,5)+(24.)*data(:,j,6) )/(120.*h) |
---|
| 387 | ! work(2,:)=(-(5./9.)*data(:,j,1)-(1./2.)*data(:,j,2) & |
---|
| 388 | ! +(1.)*data(:,j,3)+(1./18.)*data(:,j,4) )/h |
---|
| 389 | work(2,:)=( -(24.)*data(:,j,1)-(130.)*data(:,j,2) & |
---|
| 390 | +(240.)*data(:,j,3)-(120.)*data(:,j,4) & |
---|
| 391 | +(40.)*data(:,j,5)-(6.)*data(:,j,6) )/(120.*h) |
---|
| 392 | work(3,:)=( (1./36.)* ( data(:,j,5) - data(:,j,1) ) & |
---|
| 393 | + (7./9.)* ( data(:,j,4) - data(:,j,2) ) )/h |
---|
| 394 | |
---|
| 395 | do k=4,m3-3 |
---|
| 396 | work(k,:)=( (c1/6.)* ( data(:,j,k+3) - data(:,j,k-3) ) & |
---|
| 397 | + (b1/4.)* ( data(:,j,k+2) - data(:,j,k-2) ) & |
---|
| 398 | + (a1/2.)* ( data(:,j,k+1) - data(:,j,k-1) ) )/h |
---|
| 399 | enddo |
---|
| 400 | |
---|
| 401 | work(m3-2,:)=(-(1./36.)* ( data(:,j,m3-4) - data(:,j,m3) ) & |
---|
| 402 | - (7./9.)* ( data(:,j,m3-3) - data(:,j,m3-1) ) )/h |
---|
| 403 | ! work(m3-1,:)= ( (5./9.)*data(:,j,m3)+(1./2.)*data(:,j,m3-1) & |
---|
| 404 | ! -(1.)*data(:,j,m3-2)-(1./18.)*data(:,j,m3-3) )/h |
---|
| 405 | work(m3-1,:)=( (24.)*data(:,j,m3)+(130.)*data(:,j,m3-1) & |
---|
| 406 | -(240.)*data(:,j,m3-2)+(120.)*data(:,j,m3-3) & |
---|
| 407 | -(40.)*data(:,j,m3-4)+(6.)*data(:,j,m3-5) )/(120.*h) |
---|
| 408 | work(m3,:)=( (274.)*data(:,j,m3)-(600.)*data(:,j,m3-1) & |
---|
| 409 | +(600.)*data(:,j,m3-2)-(400.)*data(:,j,m3-3) & |
---|
| 410 | +(150.)*data(:,j,m3-4)-(24.)*data(:,j,m3-5) )/(120.*h) |
---|
| 411 | |
---|
| 412 | !Compute the z derivatives for plane j, overwriting the array work |
---|
| 413 | if( LAPACK_FLAG == 'double' ) then |
---|
| 414 | !!! call DGBTRS('N',m3,2,2,m1,Grid(level)%ddz,7,Grid(level)%piv_dz,work(1,1),m3,ierr) |
---|
| 415 | write(0,*) 'LAPACK_FLAG = double instead of single' |
---|
| 416 | stop |
---|
| 417 | elseif( LAPACK_FLAG == 'single' ) then |
---|
| 418 | !!! write(0,*) '********** ddz = ', Grid(level)%ddz |
---|
| 419 | !!! write(0,*) '********** piv_dz = ', Grid(level)%piv_dz |
---|
| 420 | call SGBTRS('N',m3,2,2,m1,Grid(level)%ddz,7,Grid(level)%piv_dz,work(1,1),m3,ierr) |
---|
| 421 | endif |
---|
| 422 | |
---|
| 423 | !!! do i=1,m1 |
---|
| 424 | !!! do k=1,m3 |
---|
| 425 | !!! if(abs(work(k,i)) .lt. 1.e-20) work(k,i) = 0.0 |
---|
| 426 | !!! enddo |
---|
| 427 | !!! enddo |
---|
| 428 | |
---|
| 429 | !Store the results in deriv(:,:,k) |
---|
| 430 | do i=1,m1 |
---|
| 431 | deriv(i,j,:)=work(:,i) |
---|
| 432 | enddo |
---|
| 433 | |
---|
| 434 | !If deriv. end conditions are 'periodic' (flag = 5) then |
---|
| 435 | !overwrite the endpts w/ the avg of the endpts |
---|
| 436 | if( flag .eq. 5 ) then |
---|
| 437 | work(1,:)=(deriv(:,j,1)+deriv(:,j,m3))/2. |
---|
| 438 | deriv(:,j,1)=work(1,:) |
---|
| 439 | deriv(:,j,m3)=work(1,:) |
---|
| 440 | endif |
---|
| 441 | |
---|
| 442 | !Go on to the next j plane |
---|
| 443 | enddo |
---|
| 444 | |
---|
| 445 | deallocate ( work ) |
---|
| 446 | end subroutine compact_ddz |
---|
| 447 | |
---|
| 448 | |
---|
| 449 | |
---|
| 450 | |
---|
| 451 | |
---|
| 452 | subroutine compact_ddz_mpi(data,deriv,flag,end_values,m1,m2,m3) |
---|
| 453 | ! This version computes the z derivative of a block of 3d data. |
---|
| 454 | ! There are three ways in which this routine can be used, depending |
---|
| 455 | ! on the values of nz, m3 and nzderiv. The cases are denoted CASE1 |
---|
| 456 | ! CASE2 and CASE3 as described below. |
---|
| 457 | ! |
---|
| 458 | ! CASE1: nprocs=1 and nzderiv=m3=nz |
---|
| 459 | ! Serial execution. No communication calls necessary. |
---|
| 460 | ! Derivatives implemented with single call to |
---|
| 461 | ! compact_ddz with a vector length of nzderiv. |
---|
| 462 | |
---|
| 463 | ! CASE2: nprocs>1 and nzderiv=m3 < nz |
---|
| 464 | ! Parallel execution but no data exchanges between processors |
---|
| 465 | ! required. Each processor computes z derivatives over its own |
---|
| 466 | ! data space, i.e. over vector lengths of m3. |
---|
| 467 | ! Derivatives implemented by having each processor execute |
---|
| 468 | ! compact_ddz with a vector length of nzderiv. |
---|
| 469 | ! |
---|
| 470 | ! CASE3: nprocs>1 and nzderiv=nz > m3. |
---|
| 471 | ! No processor possesses the required data locally and so |
---|
| 472 | ! data exchanges between processors are required both before |
---|
| 473 | ! and after differentiation. The algorithm is described below. |
---|
| 474 | ! On entry, each process possesses a "slice" of the total data space, |
---|
| 475 | ! i.e. a block with dimensions (m1,m2,m3). Since the compact schemes |
---|
| 476 | ! are nonlocal, i.e. all nz elements are required to compute a z derivative. |
---|
| 477 | ! In this implementation, each process will compute z derivatives within |
---|
| 478 | ! a "block of columns", i.e. a block of data of size (m1/nprocs,m2,nz) |
---|
| 479 | ! or ( (m1-1)/nprocs,m2,n2 ) if ffts are being used in x&y |
---|
| 480 | ! Since the basic data distribution consists of "blocks of rows", exchange |
---|
| 481 | ! of data between processes is required. |
---|
| 482 | ! |
---|
| 483 | ! Each process lacks (nprocs-1) blocks of size (locnx,m2,m3) |
---|
| 484 | ! needed to make up its "block of columns". On entry, each process |
---|
| 485 | ! has (nprocs-1) data "blocks" of size (locnx,m2,m3) that are |
---|
| 486 | ! not needed locally but required by other processes to complete their |
---|
| 487 | ! "block of columns". |
---|
| 488 | ! |
---|
| 489 | ! This algorithm thus proceeds as follows: |
---|
| 490 | ! (1) ship data not needed locally and receive data required to |
---|
| 491 | ! complete a "block of columns". |
---|
| 492 | ! (2) compute z derivatives over the processes "block of columns". |
---|
| 493 | ! (3) send results to other processes, receive their results |
---|
| 494 | ! (4) store the results in appropriate block locations of the original |
---|
| 495 | ! "slice" of the total data space. |
---|
| 496 | !******************************************************* |
---|
| 497 | use mpi_parameters |
---|
| 498 | use boundary_information |
---|
| 499 | use columnslab |
---|
| 500 | |
---|
| 501 | implicit none |
---|
| 502 | |
---|
| 503 | integer, intent(in) :: flag,m1,m2,m3 |
---|
| 504 | real, dimension(m1,m2,m3), intent(in) :: data |
---|
| 505 | real, dimension(m1,m2,m3), intent(out) :: deriv |
---|
| 506 | real, dimension(2), intent(in) :: end_values |
---|
| 507 | integer :: i,j,k,nzderiv,locnx |
---|
| 508 | |
---|
| 509 | include 'mpif.h' |
---|
| 510 | integer :: tag,ierr,status_array(MPI_STATUS_SIZE,2*nprocs),reqs(2*nprocs) |
---|
| 511 | integer :: dest,source,istart,iend,kstart,kend,nz |
---|
| 512 | integer :: sizeofreal |
---|
| 513 | integer :: iproc,numwords |
---|
| 514 | real, allocatable :: block_of_cols(:,:,:) |
---|
| 515 | |
---|
| 516 | if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then |
---|
| 517 | locnx = (m1-1)/nprocs |
---|
| 518 | else |
---|
| 519 | locnx = m1/nprocs |
---|
| 520 | endif |
---|
| 521 | |
---|
| 522 | |
---|
| 523 | nz=m3*nprocs ! global # of points |
---|
| 524 | nzderiv=nz ! always do global differentiation |
---|
| 525 | |
---|
| 526 | ! Both CASE1 and CASE2 require simple calls to compact_ddz |
---|
| 527 | if( nprocs .eq. 1 .and. m3 .eq. nz .and. nzderiv .eq. nz ) then |
---|
| 528 | goto 999 ! CASE1 |
---|
| 529 | elseif(nprocs .gt. 1 .and. nzderiv .eq. m3 .and. m3 .lt. nz) then |
---|
| 530 | goto 999 ! CASE2 |
---|
| 531 | elseif(nprocs .gt. 1 .and. nzderiv .eq. nz .and. nz .gt. m3) then |
---|
| 532 | ! CASE3 logic |
---|
| 533 | |
---|
| 534 | allocate ( block_of_cols(locnx,m2,nprocs*m3) ) |
---|
| 535 | block_of_cols=0. |
---|
| 536 | |
---|
| 537 | ! Create a derived data type for ease of swaps |
---|
| 538 | call MPI_TYPE_EXTENT(MPI_REAL,sizeofreal,ierr) |
---|
| 539 | ! data type for 1d section of basic data slice: m1/nprocs |
---|
| 540 | call MPI_TYPE_VECTOR(locnx,1,1,MPI_REAL,oneslice,ierr) |
---|
| 541 | ! data type for 2d section of basic data slice: m2=ny 1d sections |
---|
| 542 | call MPI_TYPE_HVECTOR(m2,1,m1*sizeofreal,oneslice,twoslice,ierr) |
---|
| 543 | ! data type for 3d section of basic data slice:locnz=nz/numprocs 2d sections |
---|
| 544 | call MPI_TYPE_HVECTOR(m3,1,m1*m2*sizeofreal,twoslice,subslice,ierr) |
---|
| 545 | ! declare data type |
---|
| 546 | call MPI_TYPE_COMMIT(subslice,ierr) |
---|
| 547 | |
---|
| 548 | !!******************************************************************************* |
---|
| 549 | !!************* DO THE DATA EXCHANGES **************************************** |
---|
| 550 | !!******************************************************************************* |
---|
| 551 | !! (1) Ship portions of my slice needed by other processes and |
---|
| 552 | !! receive data needed to fill in "block_of_cols" |
---|
| 553 | |
---|
| 554 | call slabs_to_columns(data,block_of_cols) |
---|
| 555 | |
---|
| 556 | !! (2) Compute z derivatives over my "block of columns". |
---|
| 557 | !! Store results in the array "col_of_derivs". |
---|
| 558 | |
---|
| 559 | call compact_ddz(block_of_cols,block_of_cols,flag,end_values,locnx,m2,nz) |
---|
| 560 | |
---|
| 561 | !! (3) Send derivative results to other processes |
---|
| 562 | !! and receive their results |
---|
| 563 | |
---|
| 564 | call columns_to_slabs(block_of_cols,deriv) |
---|
| 565 | |
---|
| 566 | !!******************************************************************************* |
---|
| 567 | !!************* END DATA EXCHANGES **************************************** |
---|
| 568 | !!******************************************************************************* |
---|
| 569 | |
---|
| 570 | if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then |
---|
| 571 | deriv(m1,:,:)=deriv(1,:,:) ! we didn't compute the extra bdry point |
---|
| 572 | endif |
---|
| 573 | |
---|
| 574 | |
---|
| 575 | call MPI_BARRIER(comm,ierr) |
---|
| 576 | |
---|
| 577 | !!! free data type |
---|
| 578 | call MPI_TYPE_FREE(oneslice,ierr) |
---|
| 579 | call MPI_TYPE_FREE(twoslice,ierr) |
---|
| 580 | call MPI_TYPE_FREE(subslice,ierr) |
---|
| 581 | |
---|
| 582 | deallocate(block_of_cols) |
---|
| 583 | return ! end CASE3 |
---|
| 584 | |
---|
| 585 | else |
---|
| 586 | write(0,*) ' size logic error using compact_ddz_mpi ' |
---|
| 587 | write(0,*) ' nz : ',nz |
---|
| 588 | write(0,*) ' locnz : ',m3 |
---|
| 589 | write(0,*) ' nzderiv : ',nzderiv |
---|
| 590 | write(0,*) ' nprocs : ',nprocs |
---|
| 591 | endif |
---|
| 592 | |
---|
| 593 | 999 continue ! CASE1 & CASE2 |
---|
| 594 | call compact_ddz(data,deriv,flag,end_values,m1,m2,m3) |
---|
| 595 | |
---|
| 596 | end subroutine compact_ddz_mpi |
---|
| 597 | |
---|