|
XBeach
|
00001 module spaceparams 00002 use spaceparamsdef 00003 use xmpi_module 00004 use mnemmodule 00005 implicit none 00006 save 00007 ! 00008 00009 #ifdef USEMPI 00010 00011 ! This interface an cause complains by memcache because the it is conditional on uninitialized variables. 00012 ! I don't think this is an error 00013 interface space_distribute 00014 module procedure space_distribute_matrix_real8 00015 module procedure space_distribute_matrix_integer 00016 module procedure space_distribute_block_real8 00017 module procedure space_distribute_block_integer 00018 module procedure space_distribute_block4_real8 00019 module procedure space_distribute_vector 00020 module procedure space_distribute_block_vector 00021 end interface space_distribute 00022 00023 interface space_shift_borders 00024 module procedure space_shift_borders_matrix_real8 00025 module procedure space_shift_borders_block_real8 00026 end interface space_shift_borders 00027 00028 interface space_collect 00029 module procedure space_collect_block_real8 00030 module procedure space_collect_block_integer 00031 module procedure space_collect_block4_real8 00032 module procedure space_collect_block4_integer 00033 module procedure space_collect_matrix_real8 00034 module procedure space_collect_matrix_integer 00035 end interface space_collect 00036 00037 #endif 00038 00039 contains 00040 00041 subroutine indextos(s,index,t) 00042 ! 00043 ! given s and index, this subroutine returns in t 00044 ! a pointer to the requested array 00045 ! 00046 use mnemmodule 00047 use logging_module 00048 use spaceparamsdef 00049 implicit none 00050 type (spacepars), intent(in),target :: s 00051 integer, intent(in) :: index 00052 type(arraytype), intent(out) :: t 00053 00054 if (index .lt. 1 .or. index .gt. numvars) then 00055 call writelog('els','(a,i3,a)','invalid index ',index,' in indextos. Program will stop') 00056 call halt_program 00057 endif 00058 00059 select case(index) 00060 include 'indextos.inc' 00061 end select 00062 00063 end subroutine indextos 00064 00065 subroutine index_allocate(s,par,index,choice) 00066 ! allocates, deallocates reallocates in type s, based on index 00067 ! choice: 00068 ! 'a': allocate 00069 ! 'd': deallocate 00070 ! 'r': reallocate 00071 use params 00072 use logging_module 00073 use spaceparamsdef 00074 implicit none 00075 type (spacepars), intent(inout) :: s 00076 type(parameters), intent(in) :: par 00077 integer, intent(in) :: index 00078 character(*) :: choice 00079 00080 ! the .inc files contain code for allocatable entities 00081 ! scalars will be skipped silently 00082 select case(choice(1:1)) 00083 case('a','A') 00084 select case(index) 00085 case default 00086 continue 00087 include 'index_allocate.inc' 00088 end select 00089 case ('d','D') 00090 select case(index) 00091 case default 00092 continue 00093 include 'index_deallocate.inc' 00094 end select 00095 case ('r','R') 00096 select case(index) 00097 case default 00098 continue 00099 include 'index_reallocate.inc' 00100 end select 00101 end select 00102 end subroutine index_allocate 00103 00104 logical function index_allocated(s,index) 00105 use logging_module 00106 use spaceparamsdef 00107 implicit none 00108 type (spacepars), intent(in) :: s 00109 integer, intent(in) :: index 00110 00111 logical :: r 00112 00113 r = .true. ! for scalars: function will return always .true. 00114 select case(index) 00115 case default 00116 continue 00117 include 'index_allocated.inc' 00118 end select 00119 00120 index_allocated = r 00121 00122 end function index_allocated 00123 00124 ! Generated subroutine to allocate the scalars in s 00125 subroutine space_alloc_scalars(s) 00126 use mnemmodule 00127 implicit none 00128 type(spacepars),intent(inout) :: s 00129 00130 include 'space_alloc_scalars.inc' 00131 00132 end subroutine space_alloc_scalars 00133 00134 ! Generated subroutine to allocate all arrays in s 00135 subroutine space_alloc_arrays(s,par) 00136 use mnemmodule 00137 use params 00138 implicit none 00139 type(spacepars),intent(inout) :: s 00140 type(parameters),intent(in) :: par 00141 00142 include 'space_alloc_arrays.inc' 00143 00144 end subroutine space_alloc_arrays 00145 00146 ! Generated subroutine to allocate dummy address for all arrays in s 00147 subroutine space_alloc_arrays_dummies(s,par) 00148 use mnemmodule 00149 use params 00150 implicit none 00151 type(spacepars),intent(inout) :: s 00152 type(parameters),intent(in) :: par 00153 00154 include 'space_alloc_arrays_dummies.inc' 00155 00156 end subroutine space_alloc_arrays_dummies 00157 00158 #ifdef USEMPI 00159 00160 ! copies scalars from sg to sl on xmaster, and distributes 00161 ! them 00162 subroutine space_copy_scalars(sg,sl) 00163 use mnemmodule 00164 implicit none 00165 type(spacepars),intent(inout) :: sg,sl 00166 00167 type(arraytype) :: tg,tl 00168 integer :: j 00169 00170 do j = 1,numvars 00171 call indextos(sg,j,tg) 00172 if (tg%rank .eq. 0) then 00173 call indextos(sl,j,tl) 00174 select case (tg%type) 00175 case('i') 00176 tl%i0 = tg%i0 00177 case('r') 00178 tl%r0 = tg%r0 00179 end select 00180 endif 00181 enddo 00182 end subroutine space_copy_scalars 00183 00184 ! The scalars are needed for allocating the arrays, 00185 ! we distribute them here: 00186 subroutine space_distribute_scalars(sl) 00187 use mnemmodule 00188 use xmpi_module 00189 type (spacepars) :: sl 00190 00191 type (arraytype) :: tl 00192 integer :: i 00193 logical, parameter :: toall = .true. 00194 00195 do i=1,numvars 00196 call indextos(sl,i,tl) 00197 if (tl%rank .eq. 0) then 00198 if (tl%type .eq. 'i') then 00199 call xmpi_bcast(tl%i0,toall) 00200 else 00201 call xmpi_bcast(tl%r0,toall) 00202 endif 00203 endif 00204 enddo 00205 00206 end subroutine space_distribute_scalars 00207 !________________________________________________________________________ 00208 subroutine space_distribute_matrix_real8(sl,a,b,toall) 00209 use xmpi_module 00210 use general_mpi_module 00211 implicit none 00212 type (spacepars), intent(inout) :: sl 00213 real*8, dimension(:,:), intent(in) :: a 00214 real*8, dimension(:,:), intent(out) :: b 00215 include 'space_distribute.inc' 00216 call matrix_distr(a,b,sl%is,sl%lm,sl%js,sl%ln,src,comm) 00217 end subroutine space_distribute_matrix_real8 00218 00219 !________________________________________________________________________ 00220 subroutine space_distribute_matrix_integer(sl,a,b,toall) 00221 use xmpi_module 00222 use general_mpi_module 00223 implicit none 00224 type (spacepars), intent(inout) :: sl 00225 integer, dimension(:,:), intent(in) :: a 00226 integer, dimension(:,:), intent(out) :: b 00227 include 'space_distribute.inc' 00228 00229 call matrix_distr(a,b,sl%is,sl%lm,sl%js,sl%ln,src,comm) 00230 00231 end subroutine space_distribute_matrix_integer 00232 00233 !________________________________________________________________________ 00234 subroutine space_distribute_block_real8(sl,a,b,toall) 00235 use xmpi_module 00236 use general_mpi_module 00237 implicit none 00238 type (spacepars), intent(inout) :: sl 00239 real*8, dimension(:,:,:), intent(in) :: a 00240 real*8, dimension(:,:,:), intent(out) :: b 00241 00242 real*8, dimension(1,1) :: dum 00243 integer :: i 00244 include 'space_distribute.inc' 00245 00246 do i=1,size(b,3) ! assuming that b is allocated on all processes 00247 if (iamsrc) then 00248 call matrix_distr(a(:,:,i),b(:,:,i),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00249 else 00250 call matrix_distr(dum, b(:,:,i),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00251 endif 00252 enddo 00253 00254 end subroutine space_distribute_block_real8 00255 !________________________________________________________________________ 00256 00257 subroutine space_distribute_block_integer(sl,a,b,toall) 00258 use xmpi_module 00259 use general_mpi_module 00260 implicit none 00261 type (spacepars), intent(inout) :: sl 00262 integer, dimension(:,:,:), intent(in) :: a 00263 integer, dimension(:,:,:), intent(out) :: b 00264 00265 integer, dimension(1,1) :: dum 00266 integer :: i 00267 include 'space_distribute.inc' 00268 00269 do i=1,size(b,3) ! assuming that b is allocated on all processes 00270 if(iamsrc) then 00271 call matrix_distr(a(:,:,i),b(:,:,i),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00272 else 00273 call matrix_distr(dum, b(:,:,i),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00274 endif 00275 00276 enddo 00277 00278 end subroutine space_distribute_block_integer 00279 !________________________________________________________________________ 00280 00281 subroutine space_distribute_block4_real8(sl,a,b,toall) 00282 use xmpi_module 00283 use general_mpi_module 00284 implicit none 00285 type (spacepars), intent(inout) :: sl 00286 real*8, dimension(:,:,:,:), intent(in) :: a 00287 real*8, dimension(:,:,:,:), intent(out) :: b 00288 00289 integer :: i,j 00290 real*8, dimension(1,1) :: dum 00291 include 'space_distribute.inc' 00292 00293 do i=1,size(b,3) 00294 do j=1,size(b,4) 00295 if(iamsrc) then 00296 call matrix_distr(a(:,:,i,j),b(:,:,i,j),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00297 else 00298 call matrix_distr(dum, b(:,:,i,j),sl%is,sl%lm,sl%js,sl%ln,src,comm) 00299 endif 00300 enddo 00301 enddo 00302 00303 end subroutine space_distribute_block4_real8 00304 !________________________________________________________________________ 00305 00306 subroutine space_distribute_vector(xy,sl,a,b,toall) 00307 use xmpi_module 00308 use general_mpi_module 00309 implicit none 00310 ! Not sure what this xy is doing here. 00311 ! Ok, the master process holds a vector of length ny+1 or nx+1 00312 ! This vector will be distributed among all processes. 00313 ! if xy .eq. 'x', the distribution is done according 00314 ! to the values in is and lm, i.e. the global vector 00315 ! has a size equal to the size of the first dimension of the 00316 ! global area: nx+1. 00317 ! In the other case, the distribution is done according to 00318 ! the values in js and ln, i.e. the global vector has a size 00319 ! equal to the second dimension of the global area: ny+1 00320 ! 00321 character, intent(in) :: xy 00322 type (spacepars), intent(inout), target :: sl 00323 real*8, dimension(:), intent(in) :: a 00324 real*8, dimension(:), intent(out) :: b 00325 00326 integer, dimension(:), pointer :: ijs,lmn 00327 include 'space_distribute.inc' 00328 00329 if(xy .eq.'x') then 00330 ijs => sl%is 00331 lmn => sl%lm 00332 else 00333 ijs => sl%js 00334 lmn => sl%ln 00335 endif 00336 call vector_distr(a,b,ijs,lmn,src,comm) 00337 00338 end subroutine space_distribute_vector 00339 !________________________________________________________________________ 00340 00341 subroutine space_distribute_block_vector(xy,sl,a,b,toall) 00342 use xmpi_module 00343 use general_mpi_module 00344 implicit none 00345 character, intent(in) :: xy 00346 type (spacepars), intent(in) :: sl 00347 real*8, dimension(:,:), intent(in) :: a 00348 real*8, dimension(:,:), intent(out) :: b 00349 00350 ! integer :: i 00351 include 'space_distribute.inc' 00352 00353 !DF do i=1,sl%ntheta 00354 ! do i=1,size(b,2) 00355 ! call space_distribute(xy,sl,a(:,i),b(:,i)) 00356 ! enddo 00357 00358 select case(xy) 00359 case('y') 00360 call vector_distr(a,b,sl%js,sl%ln,src,comm) 00361 case default 00362 print *,'error in space_distribute_block_vector, other than "y" not tested' 00363 call xmpi_abort() 00364 end select 00365 00366 end subroutine space_distribute_block_vector 00367 !________________________________________________________________________ 00368 00369 00370 subroutine space_distribute_space(sg,sl,par) 00371 use xmpi_module 00372 use logging_module 00373 use general_mpi_module 00374 use params 00375 use mnemmodule 00376 00377 implicit none 00378 type(spacepars), intent(inout) :: sg 00379 type(spacepars), intent(inout) :: sl 00380 type(parameters) :: par 00381 00382 integer :: i,j,lid,eid, wid 00383 real*8, pointer, dimension(:) :: vectorg, vectorl 00384 type (arraytype) :: tg, tl 00385 00386 logical, parameter :: toall = .true. 00387 integer, parameter :: nbord = 2 00388 character(1000) :: txt 00389 00390 ! 00391 ! This subroutine takes care that all contents of the global 00392 ! space is distributed to the local space. 00393 ! 00394 00395 ! allocate scalars 00396 00397 call space_alloc_scalars(sl) 00398 00399 ! copy scalars to sl, only on master 00400 ! distributing will take place later 00401 00402 00403 if(xmaster) then 00404 call space_copy_scalars(sg,sl) 00405 endif 00406 00407 ! copy scalars to all processes, nx and ny will be adapted later 00408 call space_distribute_scalars(sl) 00409 call space_distribute_scalars(sg) 00410 00411 ! Also, the isleft, isright, istop and isbot logicals from 00412 ! the xmpi module are put in sg and sl. 00413 ! 00414 00415 ! 00416 ! Distribute is,js,ln,lm,isleft,isright,istop,isbot 00417 ! 00418 00419 if(xmaster) then 00420 allocate(sg%is(xmpi_size)) 00421 allocate(sg%js(xmpi_size)) 00422 allocate(sg%lm(xmpi_size)) 00423 allocate(sg%ln(xmpi_size)) 00424 allocate(sg%isleft(xmpi_size)) 00425 allocate(sg%isright(xmpi_size)) 00426 allocate(sg%istop(xmpi_size)) 00427 allocate(sg%isbot(xmpi_size)) 00428 endif 00429 00430 00431 allocate(sl%is(xmpi_size)) 00432 allocate(sl%js(xmpi_size)) 00433 allocate(sl%lm(xmpi_size)) 00434 allocate(sl%ln(xmpi_size)) 00435 allocate(sl%isleft(xmpi_size)) 00436 allocate(sl%isright(xmpi_size)) 00437 allocate(sl%istop(xmpi_size)) 00438 allocate(sl%isbot(xmpi_size)) 00439 00440 if(xmaster) then 00441 call det_submatrices(sg%nx+1, sg%ny+1, xmpi_m, xmpi_n, & 00442 sg%is, sg%lm, sg%js, sg%ln, & 00443 sg%isleft, sg%isright, sg%istop, sg%isbot) 00444 call writelog('l','','--------------------------------') 00445 call writelog('l','','MPI implementation: ') 00446 call writelog('sl','','Distribution of matrix on processors') 00447 call writelog('sl','',' proc is lm js ln') 00448 do i=1,xmpi_size 00449 call writelog('sl','(i5,i5,i5,i5,i5)',i-1,sg%is(i),sg%lm(i),sg%js(i),sg%ln(i)) 00450 enddo 00451 call writelog('ls','',' proc left right top bot') 00452 do i=1,xmpi_size 00453 call writelog('ls','',i-1,sg%isleft(i),sg%isright(i),sg%istop(i),sg%isbot(i)) 00454 enddo 00455 call writelog('l','','--------------------------------') 00456 endif 00457 00458 ! determine the computational regions, i.e. the parts of the matrices 00459 ! that are computed in the processes. 00460 ! In general, the 2 borders left, right, top and bottom are not 00461 ! computed but received from the neighbours. 00462 ! Exceptions for the rows and colomns of left, right, top or 00463 ! bottom processes. 00464 ! 00465 ! icgs, icge: start and end indices in global coordinates 1st dimension 00466 ! jcgs: jcge: start and end indices in global coordinates 2nd dimension 00467 ! icls, icle: start and end indices in local coordinates 1st dimension 00468 ! jcls, jcle: start and end indices in local coordinates 2nd dimension 00469 00470 if (xmaster) then 00471 allocate (sg%icgs(xmpi_size)) 00472 allocate (sg%icge(xmpi_size)) 00473 allocate (sg%jcgs(xmpi_size)) 00474 allocate (sg%jcge(xmpi_size)) 00475 allocate (sg%icls(xmpi_size)) 00476 allocate (sg%icle(xmpi_size)) 00477 allocate (sg%jcls(xmpi_size)) 00478 allocate (sg%jcle(xmpi_size)) 00479 endif 00480 00481 allocate (sl%icgs(xmpi_size)) 00482 allocate (sl%icge(xmpi_size)) 00483 allocate (sl%jcgs(xmpi_size)) 00484 allocate (sl%jcge(xmpi_size)) 00485 allocate (sl%icls(xmpi_size)) 00486 allocate (sl%icle(xmpi_size)) 00487 allocate (sl%jcls(xmpi_size)) 00488 allocate (sl%jcle(xmpi_size)) 00489 00490 if (xmaster) then 00491 do i = 1,xmpi_size 00492 sg%icgs(i) = sg%is(i) + nbord 00493 sg%icge(i) = sg%is(i) + sg%lm(i) - 1 - nbord 00494 sg%jcgs(i) = sg%js(i) + nbord 00495 sg%jcge(i) = sg%js(i) + sg%ln(i) - 1 - nbord 00496 sg%icls(i) = nbord + 1 00497 sg%icle(i) = sg%lm(i) - nbord 00498 sg%jcls(i) = nbord + 1 00499 sg%jcle(i) = sg%ln(i) - nbord 00500 00501 if (sg%istop(i)) then 00502 sg%icgs(i) = sg%icgs(i) - nbord 00503 sg%icls(i) = sg%icls(i) - nbord 00504 endif 00505 00506 if (sg%isbot(i)) then 00507 sg%icge(i) = sg%icge(i) + nbord 00508 sg%icle(i) = sg%icle(i) + nbord 00509 endif 00510 00511 if (sg%isleft(i)) then 00512 sg%jcgs(i) = sg%jcgs(i) - nbord 00513 sg%jcls(i) = sg%jcls(i) - nbord 00514 endif 00515 00516 if (sg%isright(i)) then 00517 sg%jcge(i) = sg%jcge(i) + nbord 00518 sg%jcle(i) = sg%jcle(i) + nbord 00519 endif 00520 00521 enddo 00522 endif 00523 00524 if (xmaster) then 00525 sl%is = sg%is 00526 sl%js = sg%js 00527 sl%lm = sg%lm 00528 sl%ln = sg%ln 00529 sl%isleft = sg%isleft 00530 sl%isright = sg%isright 00531 sl%istop = sg%istop 00532 sl%isbot = sg%isbot 00533 00534 sl%icgs = sg%icgs 00535 sl%icge = sg%icge 00536 sl%jcgs = sg%jcgs 00537 sl%jcge = sg%jcge 00538 sl%icls = sg%icls 00539 sl%icle = sg%icle 00540 sl%jcls = sg%jcls 00541 sl%jcle = sg%jcle 00542 endif 00543 00544 call xmpi_bcast(sl%is,toall) 00545 call xmpi_bcast(sl%js,toall) 00546 call xmpi_bcast(sl%lm,toall) 00547 call xmpi_bcast(sl%ln,toall) 00548 call xmpi_bcast(sl%isleft,toall) 00549 call xmpi_bcast(sl%isright,toall) 00550 call xmpi_bcast(sl%istop,toall) 00551 call xmpi_bcast(sl%isbot,toall) 00552 00553 call xmpi_bcast(sl%icgs,toall) 00554 call xmpi_bcast(sl%icge,toall) 00555 call xmpi_bcast(sl%jcgs,toall) 00556 call xmpi_bcast(sl%jcge,toall) 00557 call xmpi_bcast(sl%icls,toall) 00558 call xmpi_bcast(sl%icle,toall) 00559 call xmpi_bcast(sl%jcls,toall) 00560 call xmpi_bcast(sl%jcle,toall) 00561 00562 if (xmaster) then 00563 call writelog('l','','--------------------------------') 00564 call writelog('sl','','computational domains on processors') 00565 call writelog('sl','',' proc icgs icge jcgs jcge icls icle jcls jcle') 00566 do i=1,xmpi_size 00567 write(txt,'(9i7)')i-1,sg%icgs(i),sg%icge(i),sg%jcgs(i),sg%jcge(i), & 00568 sg%icls(i),sg%icle(i),sg%jcls(i),sg%jcle(i) 00569 call writelog('sl','',txt) 00570 enddo 00571 call writelog('l','','--------------------------------') 00572 endif 00573 00574 ! 00575 ! compute the values for local nx and ny 00576 ! 00577 00578 if(xcompute) then 00579 sl%nx = sl%lm(xmpi_rank+1) - 1 00580 sl%ny = sl%ln(xmpi_rank+1) - 1 00581 endif 00582 00583 00584 ! 00585 ! allocate all arrays in sl 00586 ! 00587 00588 call space_alloc_arrays(sl,par) 00589 00590 ! for each variable in sg, find out how to distribute it 00591 ! and distribute 00592 ! 00593 00594 do i = 1,numvars 00595 00596 call indextos(sg,i,tg) 00597 call indextos(sl,i,tl) 00598 select case (tl%btype) 00599 case('b') ! have to broadcast this 00600 select case (tl%type) 00601 case('i') 00602 select case(tl%rank) 00603 case(0) ! do nothing, scalars are already in place 00604 ! Have to be prudent here. all scalars are broadcasted, exept nx 00605 ! and ny 00606 !if (tl%name .ne. mnem_nx .and. tl%name .ne. mnem_ny) then 00607 ! if(xmaster) then 00608 ! tl%i0 = tg%i0 00609 ! endif 00610 ! call xmpi_bcast(tl%i0) 00611 !endif 00612 case(1) 00613 if(xmaster) then 00614 tl%i1 = tg%i1 00615 endif 00616 call xmpi_bcast(tl%i1,toall) 00617 case(2) ! wwvv : try to fix error with integer pdish 00618 if(xmaster) then 00619 tl%i2 = tg%i2 00620 endif 00621 call xmpi_bcast(tl%i2,toall) 00622 case default 00623 goto 100 00624 end select ! rank 00625 case('r') 00626 select case(tl%rank) 00627 case(0) ! do nothing here, scalars are already in place 00628 ! Have to be prudent here. all scalars are broadcasted, exept nx 00629 ! and ny 00630 ! here only real*8 is broadcasted, so no check necessary here 00631 !if(xmaster) then 00632 ! tl%r0 = tg%r0 00633 !endif 00634 !call xmpi_bcast(tl%r0) 00635 case(1) 00636 if(xmaster) then 00637 tl%r1 = tg%r1 00638 endif 00639 call xmpi_bcast(tl%r1,toall) 00640 case(2) 00641 if(xmaster) then 00642 tl%r2 = tg%r2 00643 endif 00644 call xmpi_bcast(tl%r2,toall) 00645 case default 00646 goto 100 00647 end select ! rank 00648 case default 00649 goto 100 00650 end select ! type 00651 case('d') ! have to distribute this 00652 select case(tl%type) 00653 case ('i') 00654 select case(tl%rank) 00655 case(1) ! wwvv 2013: superfluous, there is a default case 00656 goto 100 00657 case(2) 00658 call space_distribute(sl,tg%i2,tl%i2) 00659 case(3) 00660 call space_distribute(sl,tg%i3,tl%i3) 00661 case default 00662 goto 100 00663 end select ! rank 00664 case ('r') 00665 select case(tl%rank) 00666 case(1) 00667 ! Robert: these don't exist anymore? 00668 ! case(mnem_xz,mnem_xu) 00669 ! call space_distribute_vector('x',sl,tg%r1,tl%r1) 00670 ! case(mnem_yz, mnem_yv, mnem_bi) 00671 select case(tl%name) 00672 case(mnem_bi) 00673 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00674 case(mnem_runup) 00675 ! Not sure why vector expects a name.... 00676 ! This name is x or y, it relates to the length of the vector. 00677 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00678 case(mnem_hrunup) 00679 ! Not sure why vector expects a name.... 00680 ! This name is x or y, it relates to the length of the vector. 00681 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00682 case(mnem_xhrunup) 00683 ! Not sure why vector expects a name.... 00684 ! This name is x or y, it relates to the length of the vector. 00685 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00686 case(mnem_strucslope) 00687 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00688 case(mnem_istruct) 00689 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00690 case(mnem_iwl) 00691 call space_distribute_vector('y',sl,tg%r1,tl%r1) 00692 case default 00693 goto 100 00694 end select 00695 case(2) 00696 call space_distribute(sl,tg%r2,tl%r2) 00697 case(3) 00698 call space_distribute(sl,tg%r3,tl%r3) 00699 case(4) 00700 call space_distribute(sl,tg%r4,tl%r4) 00701 case default 00702 continue 00703 end select ! rank 00704 case default 00705 goto 100 00706 end select ! type 00707 case('2') 00708 ! the umean case 00709 ! dimension 2,s%ny+1 00710 if(xmaster) then 00711 allocate(vectorg(sg%ny+1)) 00712 else 00713 allocate(vectorg(1)) 00714 vectorg = 0 00715 endif 00716 allocate(vectorl(sl%ny+1)) 00717 ! Let's keep this 2 case explicit. (because tg%r2 is only known on master) 00718 ! This one is distributed as 2 vectors. 00719 do j=1,2 00720 if (xmaster) then 00721 vectorg = tg%r2(j,:) 00722 endif 00723 call space_distribute("y", sl,vectorg,vectorl) 00724 tl%r2(j,:) = vectorl 00725 enddo 00726 if(xmaster) then 00727 deallocate(vectorg) 00728 endif 00729 deallocate(vectorl) 00730 case default 00731 goto 100 00732 end select ! btype 00733 enddo ! numvars 00734 return 00735 00736 100 continue 00737 call writelog('sel','',xmpi_rank,': Error in space_distribute_space, trying to distribute:') 00738 call get_logfileid(lid,eid, wid) 00739 call printvar(tl,lid,eid, wid) 00740 call halt_program 00741 return 00742 00743 end subroutine space_distribute_space 00744 00745 subroutine space_who_has(sl,i,j,p) 00746 ! determine which process contains element(i,j) 00747 ! sl: local spacepars 00748 ! i,j: global indices 00749 ! p: process in xmpi_ocomm which contains the element 00750 00751 ! icgs, icge: start and end indices in global coordinates 1st dimension 00752 ! jcgs: jcge: start and end indices in global coordinates 2nd dimension 00753 implicit none 00754 type(spacepars), intent(in) :: sl 00755 integer, intent(in) :: i,j 00756 integer, intent(out) :: p 00757 00758 integer k 00759 00760 p = -123 00761 00762 do k = 1,xmpi_size 00763 if (i .ge. sl%icgs(k) .and. i .le. sl%icge(k) .and. & 00764 & j .ge. sl%jcgs(k) .and. j .le. sl%jcge(k)) then 00765 p = xmpi_rank_to_orank(k-1) 00766 exit 00767 endif 00768 enddo 00769 00770 end subroutine space_who_has 00771 00772 subroutine space_global_to_local(sl,ig,jg,il,jl) 00773 ! given global coordinates (ig,jg), compute local coordinates (il,jl) 00774 implicit none 00775 type(spacepars), intent(in) :: sl 00776 integer, intent(in) :: ig,jg 00777 integer, intent(out) :: il,jl 00778 00779 il = ig - sl%is(xmpi_rank+1) + 1 00780 jl = jg - sl%js(xmpi_rank+1) + 1 00781 00782 end subroutine space_global_to_local 00783 00784 subroutine space_local_to_global(sl,il,jl,ig,jg) 00785 ! given local coordinates (il,jl), compute global coordinates (ig,jg) 00786 implicit none 00787 type(spacepars), intent(in) :: sl 00788 integer, intent(in) :: il,jl 00789 integer, intent(out) :: ig,jg 00790 00791 ig = il + sl%is(xmpi_rank+1) - 1 00792 jg = jl + sl%js(xmpi_rank+1) - 1 00793 00794 end subroutine space_local_to_global 00795 00796 subroutine space_shift_borders_matrix_real8(a) 00797 use general_mpi_module 00798 use xmpi_module 00799 implicit none 00800 real*8, intent(inout),dimension(:,:) :: a 00801 call shift_borders_matrix_real8(a,xmpi_left,xmpi_right, & 00802 xmpi_top,xmpi_bot,xmpi_comm) 00803 end subroutine space_shift_borders_matrix_real8 00804 00805 subroutine space_shift_borders_block_real8(a) 00806 use general_mpi_module 00807 use xmpi_module 00808 implicit none 00809 real*8, intent(inout),dimension(:,:,:) :: a 00810 00811 integer i 00812 00813 do i=1,size(a,3) 00814 call shift_borders_matrix_real8(a(:,:,i),xmpi_left,xmpi_right,& 00815 xmpi_top,xmpi_bot,xmpi_comm) 00816 enddo 00817 end subroutine space_shift_borders_block_real8 00818 00819 ! 00820 ! wwvv a subtle point with the collect subroutines: the second 00821 ! argument: the matrix wherein the submatrices are to be collected, 00822 ! does not have to be available on the non-master processes, so 00823 ! the dimensions are not defined. 00824 ! The third argument is always defined, on master and non-master 00825 ! processes, so its dimensions (notably the 3rd in the block subroutines 00826 ! are available 00827 ! 00828 ! parameters of the space_collect routines: 00829 ! s: spacepars: LOCAL s 00830 ! a: output: in this matrix the submatrices are collected 00831 ! b: input: the local submatrix 00832 00833 subroutine space_collect_block_real8(s,a,b) 00834 use general_mpi_module 00835 use xmpi_module 00836 implicit none 00837 type(spacepars), intent(in) :: s 00838 real*8, dimension(:,:,:), intent(out) :: a 00839 real*8, dimension(:,:,:), intent(in) :: b 00840 00841 integer i 00842 real*8, dimension(1,1) :: dum 00843 00844 do i = 1,size(b,3) 00845 if (xomaster) then 00846 call matrix_coll(a(:,:,i),b(:,:,i),s%is,s%lm,s%js,s%ln, & 00847 s%isleft,s%isright,s%istop,s%isbot, & 00848 xmpi_omaster,xmpi_ocomm) 00849 else 00850 call matrix_coll(dum, b(:,:,i),s%is,s%lm,s%js,s%ln, & 00851 s%isleft,s%isright,s%istop,s%isbot, & 00852 xmpi_omaster,xmpi_ocomm) 00853 endif 00854 enddo 00855 00856 end subroutine space_collect_block_real8 00857 00858 subroutine space_collect_block_integer(s,a,b) 00859 use general_mpi_module 00860 use xmpi_module 00861 implicit none 00862 type(spacepars), intent(in) :: s 00863 integer, dimension(:,:,:), intent(out) :: a 00864 integer, dimension(:,:,:), intent(in) :: b 00865 00866 integer i 00867 00868 !real*8, dimension(:,:,:), allocatable :: ra,rb 00869 ! wwvv changed this into 00870 ! wwvv huh? TODO 00871 integer, dimension(:,:,:), allocatable :: ra,rb 00872 integer :: m,n,o 00873 integer, dimension(1,1) :: dum 00874 00875 m = size(b,1) 00876 n = size(b,2) 00877 o = size(b,3) 00878 00879 allocate(rb(m,n,o)) 00880 00881 if (xomaster) then 00882 m = size(a,1) 00883 n = size(a,2) 00884 o = size(a,3) 00885 allocate(ra(m,n,o)) 00886 else 00887 allocate(ra(1,1,1)) 00888 endif 00889 00890 rb = b 00891 00892 do i = 1,o 00893 if(xomaster) then 00894 call matrix_coll(ra(:,:,i),rb(:,:,i),s%is,s%lm,s%js,s%ln, & 00895 s%isleft,s%isright,s%istop,s%isbot, & 00896 xmpi_omaster,xmpi_ocomm) 00897 else 00898 call matrix_coll(dum, rb(:,:,i),s%is,s%lm,s%js,s%ln, & 00899 s%isleft,s%isright,s%istop,s%isbot, & 00900 xmpi_omaster,xmpi_ocomm) 00901 endif 00902 enddo 00903 00904 if (xomaster) then 00905 a = ra 00906 endif 00907 00908 deallocate(ra,rb) 00909 00910 end subroutine space_collect_block_integer 00911 00912 subroutine space_collect_block4_real8(s,a,b) 00913 use general_mpi_module 00914 use xmpi_module 00915 implicit none 00916 type(spacepars), intent(in) :: s 00917 real*8, dimension(:,:,:,:), intent(out) :: a 00918 real*8, dimension(:,:,:,:), intent(in) :: b 00919 00920 integer i,j 00921 real*8, dimension(1,1) :: dum 00922 00923 do j = 1,size(b,4) 00924 do i = 1,size(b,3) 00925 if (xomaster) then 00926 call matrix_coll(a(:,:,i,j),b(:,:,i,j),s%is,s%lm,s%js,s%ln, & 00927 s%isleft,s%isright,s%istop,s%isbot, & 00928 xmpi_omaster,xmpi_ocomm) 00929 else 00930 call matrix_coll(dum, b(:,:,i,j),s%is,s%lm,s%js,s%ln, & 00931 s%isleft,s%isright,s%istop,s%isbot, & 00932 xmpi_omaster,xmpi_ocomm) 00933 endif 00934 enddo 00935 enddo 00936 00937 end subroutine space_collect_block4_real8 00938 00939 subroutine space_collect_block4_integer(s,a,b) 00940 use general_mpi_module 00941 use xmpi_module 00942 implicit none 00943 type(spacepars), intent(in) :: s 00944 integer, dimension(:,:,:,:), intent(out) :: a 00945 integer, dimension(:,:,:,:), intent(in) :: b 00946 00947 integer i,j 00948 integer :: dum(1,1) 00949 00950 do j = 1,size(b,4) 00951 do i = 1,size(b,3) 00952 if(xomaster) then 00953 call matrix_coll(a(:,:,i,j),b(:,:,i,j),s%is,s%lm,s%js,s%ln, & 00954 s%isleft,s%isright,s%istop,s%isbot, & 00955 xmpi_omaster,xmpi_ocomm) 00956 else 00957 call matrix_coll(dum, b(:,:,i,j),s%is,s%lm,s%js,s%ln, & 00958 s%isleft,s%isright,s%istop,s%isbot, & 00959 xmpi_omaster,xmpi_ocomm) 00960 endif 00961 enddo 00962 enddo 00963 00964 end subroutine space_collect_block4_integer 00965 00966 subroutine space_collect_matrix_real8(s,a,b) 00967 use general_mpi_module 00968 use xmpi_module 00969 implicit none 00970 type(spacepars), intent(in) :: s 00971 real*8, dimension(:,:), intent(out) :: a 00972 real*8, dimension(:,:), intent(in) :: b 00973 00974 call matrix_coll(a,b,s%is,s%lm,s%js,s%ln, & 00975 s%isleft,s%isright,s%istop,s%isbot, & 00976 xmpi_omaster,xmpi_ocomm) 00977 00978 end subroutine space_collect_matrix_real8 00979 00980 subroutine space_collect_vector_real8(loc,s,a,b) 00981 use general_mpi_module 00982 use xmpi_module 00983 implicit none 00984 type(spacepars), intent(in) :: s 00985 real*8, dimension(:,:), intent(out) :: a 00986 real*8, dimension(:,:), intent(in) :: b 00987 character*1, intent(in) :: loc 00988 select case (loc) 00989 case ('t') 00990 call matrix_coll(a,b,s%is,s%lm,s%js,s%ln, & 00991 s%isleft,s%isright,s%istop,s%isbot, & 00992 xmpi_omaster,xmpi_ocomm) 00993 case ('b') 00994 call matrix_coll(a,b,s%is,s%lm,s%js,s%ln, & 00995 s%isleft,s%isright,s%istop,s%isbot, & 00996 xmpi_omaster,xmpi_ocomm) 00997 case ('l') 00998 call matrix_coll(a,b,s%is,s%lm,s%js,s%ln, & 00999 s%isleft,s%isright,s%istop,s%isbot, & 01000 xmpi_omaster,xmpi_ocomm) 01001 case ('r') 01002 call matrix_coll(a,b,s%is,s%lm,s%js,s%ln, & 01003 s%isleft,s%isright,s%istop,s%isbot, & 01004 xmpi_omaster,xmpi_ocomm) 01005 end select 01006 01007 end subroutine space_collect_vector_real8 01008 01009 subroutine space_collect_matrix_integer(s,a,b) 01010 use general_mpi_module 01011 use xmpi_module 01012 implicit none 01013 type(spacepars), intent(in) :: s 01014 integer, dimension(:,:), intent(out) :: a 01015 integer, dimension(:,:), intent(in) :: b 01016 01017 ! not used often, so we convert the integers to real*8, 01018 ! collect and convert back. 01019 ! if this routine becomes heavily used, than a special 01020 ! matrix_coll has to be made. (Now we understand why 01021 ! C++ has templates) 01022 01023 !real*8, dimension(:,:), allocatable :: ra,rb 01024 ! wwvv changed this into 01025 integer, dimension(:,:), allocatable :: ra,rb 01026 01027 integer :: m,n 01028 01029 m = size(b,1) 01030 n = size(b,2) 01031 01032 allocate(rb(m,n)) 01033 01034 if (xomaster) then 01035 m = size(a,1) 01036 n = size(a,2) 01037 allocate(ra(m,n)) 01038 else 01039 allocate(ra(1,1)) 01040 endif 01041 01042 rb = b 01043 01044 call matrix_coll(ra,rb,s%is,s%lm,s%js,s%ln, & 01045 s%isleft,s%isright,s%istop,s%isbot, & 01046 xmpi_omaster,xmpi_ocomm) 01047 01048 if (xomaster) then 01049 a = ra 01050 endif 01051 01052 deallocate(ra,rb) 01053 01054 end subroutine space_collect_matrix_integer 01055 01056 #endif 01057 01058 01059 #ifdef USEMPI 01060 ! 01061 ! collects data from processes in master 01062 ! using the index number of the variable to be 01063 ! collected 01064 ! 01065 subroutine space_collect_index(sg,sl,par,index) 01066 use xmpi_module 01067 use mnemmodule 01068 use logging_module 01069 use params 01070 type(spacepars) :: sg 01071 type(spacepars), intent(in) :: sl 01072 type(parameters) :: par 01073 integer, intent(in) :: index 01074 01075 integer :: lid,eid, wid,ierr 01076 type(arraytype) :: tg,tl 01077 logical :: iscollected 01078 01079 #ifdef USEMPE 01080 call MPE_Log_event(event_coll_start,0,'cstart') 01081 #endif 01082 01083 if (xomaster) then 01084 iscollected = sg%collected(index) 01085 endif 01086 01087 call MPI_Bcast(iscollected,1,MPI_LOGICAL,xmpi_omaster,xmpi_ocomm,ierr) 01088 01089 ! return if this alreay has been collected 01090 if(iscollected) then 01091 return 01092 endif 01093 01094 #ifdef USEMPI 01095 if (xomaster) then 01096 call index_allocate(sg,par,index,'r') 01097 endif 01098 #endif 01099 01100 call indextos(sl,index,tl) 01101 call indextos(sg,index,tg) 01102 01103 select case(tl%type) 01104 case('i') 01105 select case(tl%rank) 01106 case(0) ! nothing to do 01107 case(2) 01108 call space_collect(sl, tg%i2, tl%i2) 01109 case(3) 01110 call space_collect(sl, tg%i3, tl%i3) 01111 case default ! case 1 and 4 are not handled 01112 goto 100 01113 end select ! rank 01114 case('r') 01115 select case(tl%rank) 01116 case(0) ! nothing to do 01117 case(2) 01118 !if (tl%name .eq. mnem_umean) then 01119 ! goto 100 01120 !endif 01121 call space_collect(sl,tg%r2,tl%r2) 01122 case(3) 01123 call space_collect(sl,tg%r3,tl%r3) 01124 case(4) 01125 call space_collect(sl,tg%r4,tl%r4) 01126 case default 01127 end select ! rank 01128 case default 01129 end select ! type 01130 01131 #ifdef USEMPE 01132 call MPE_Log_event(event_coll_start,0,'cend') 01133 #endif 01134 01135 sg%collected(index) = .true. 01136 return 01137 01138 100 continue 01139 call writelog('lse','','Problem in space_collect_index with variable ',trim(tg%name ) ) 01140 call writelog('lse','','Don''t know how to collect that on the masternode') 01141 call get_logfileid(lid,eid, wid) 01142 call printvar(tl,lid,eid, wid) 01143 call halt_program 01144 01145 end subroutine space_collect_index 01146 01147 subroutine space_collect_mnem(sg,sl,par,mnem) 01148 use params 01149 type(spacepars) :: sg 01150 type(spacepars), intent(in) :: sl 01151 type(parameters) :: par 01152 character(*) :: mnem 01153 01154 integer index 01155 index = chartoindex(mnem) 01156 call space_collect_index(sg,sl,par,index) 01157 01158 end subroutine space_collect_mnem 01159 #endif 01160 01161 subroutine gridprops (s) 01162 01163 IMPLICIT NONE 01164 01165 type(spacepars),target :: s 01166 ! Temporary local arrays and variables 01167 integer :: i, j 01168 real*8,dimension(:,:),allocatable :: xc ! x-coordinate c-points 01169 real*8,dimension(:,:),allocatable :: yc ! y-coordinate c-points 01170 real*8 :: dsdnu ! surface of cell centered around u-point 01171 real*8 :: dsdnv ! surface of cell centered around v-point 01172 real*8 :: dsdnz ! surface of cell centered around z-point 01173 real*8 :: x1,y1,x2,y2,x3,y3,x4,y4 01174 01175 allocate (xc(s%nx+1,s%ny+1)) 01176 allocate (yc(s%nx+1,s%ny+1)) 01177 01178 ! x and y were read in grid_bathy.f90 (initialize.f90) and can be either (cartesian) world coordinates or 01179 ! XBeach coordinates (xori, yori and alfa are nonzero). Here x and y are transformed to cartesian world coordinates. 01180 ! XBeach performs all computations (after implementation of curvi-lineair option) world coordinate grid. 01181 01182 ! world coordinates of z-points 01183 s%xz=s%xori+s%x*cos(s%alfa)-s%y*sin(s%alfa) 01184 s%yz=s%yori+s%x*sin(s%alfa)+s%y*cos(s%alfa) 01185 01186 ! world coordinates of u-points 01187 do j=1,s%ny+1 01188 do i=1,s%nx 01189 s%xu(i,j)=.5d0*(s%xz(i,j)+s%xz(i+1,j)) 01190 s%yu(i,j)=.5d0*(s%yz(i,j)+s%yz(i+1,j)) 01191 enddo 01192 s%xu(s%nx+1,j)=1.5d0*s%xz(s%nx+1,j)-0.5d0*s%xz(s%nx,j) 01193 s%yu(s%nx+1,j)=1.5d0*s%yz(s%nx+1,j)-0.5d0*s%yz(s%nx,j) 01194 enddo 01195 01196 ! world coordinates of v-points 01197 if (s%ny>0) then 01198 do i=1,s%nx+1 01199 do j=1,s%ny 01200 s%xv(i,j)=.5d0*(s%xz(i,j)+s%xz(i,j+1)) 01201 s%yv(i,j)=.5d0*(s%yz(i,j)+s%yz(i,j+1)) 01202 enddo 01203 s%xv(i,s%ny+1)=1.5d0*s%xz(i,s%ny+1)-0.5d0*s%xz(i,s%ny) 01204 s%yv(i,s%ny+1)=1.5d0*s%yz(i,s%ny+1)-0.5d0*s%yz(i,s%ny) 01205 enddo 01206 else 01207 s%xv=s%xz 01208 s%yv=s%yz 01209 endif 01210 01211 ! world coordinates of corner points 01212 if (s%ny>0) then 01213 do j=1,s%ny 01214 do i=1,s%nx 01215 xc(i,j)=.25d0*(s%xz(i,j)+s%xz(i+1,j)+s%xz(i,j+1)+s%xz(i+1,j+1)) 01216 yc(i,j)=.25d0*(s%yz(i,j)+s%yz(i+1,j)+s%yz(i,j+1)+s%yz(i+1,j+1)) 01217 enddo 01218 xc(s%nx+1,j)=0.5d0*(s%xu(s%nx+1,j)+s%xu(s%nx+1,j+1)) 01219 yc(s%nx+1,j)=0.5d0*(s%yu(s%nx+1,j)+s%yu(s%nx+1,j+1)) 01220 enddo 01221 do i=1,s%nx 01222 xc(i,s%ny+1)=0.5d0*(s%xv(i,s%ny+1)+s%xv(i+1,s%ny+1)) 01223 yc(i,s%ny+1)=0.5d0*(s%yv(i,s%ny+1)+s%yv(i+1,s%ny+1)) 01224 enddo 01225 xc(s%nx+1,s%ny+1)=1.5d0*s%xu(s%nx+1,s%ny+1)-0.5*s%xu(s%nx,s%ny+1) 01226 yc(s%nx+1,s%ny+1)=1.5d0*s%yu(s%nx+1,s%ny+1)-0.5*s%yu(s%nx,s%ny+1) 01227 else 01228 xc=s%xu 01229 yc=s%yu 01230 endif 01231 01232 ! s%dsu 01233 01234 do j=1,s%ny+1 01235 do i=1,s%nx 01236 s%dsu(i,j)=((s%xz(i+1,j)-s%xz(i,j))**2+(s%yz(i+1,j)-s%yz(i,j))**2)**(0.5d0) 01237 enddo 01238 s%dsu(s%nx+1,j)=s%dsu(s%nx,j) 01239 enddo 01240 01241 ! s%dsz 01242 01243 do j=1,s%ny+1 01244 do i=2,s%nx+1 01245 s%dsz(i,j)=((s%xu(i,j)-s%xu(i-1,j))**2+(s%yu(i,j)-s%yu(i-1,j))**2)**(0.5d0) 01246 enddo 01247 s%dsz(1,j)=s%dsz(2,j) ! Robert -> was i=s%nx+1 (calculated in loop), so i=1 more likely 01248 enddo 01249 01250 ! s%dsv 01251 01252 if (s%ny>0) then 01253 do j=1,s%ny+1 01254 do i=2,s%nx+1 01255 s%dsv(i,j)=((xc(i,j)-xc(i-1,j))**2+(yc(i,j)-yc(i-1,j))**2)**(0.5d0) 01256 enddo 01257 s%dsv(1,j)=s%dsv(2,j) ! Robert, need to have this for wave advec 01258 s%dsv(s%nx+1,j)=s%dsv(s%nx,j) 01259 enddo 01260 else 01261 s%dsv=s%dsz 01262 endif 01263 01264 ! s%dsc 01265 01266 if (s%ny>0) then 01267 do j=1,s%ny+1 01268 do i=1,s%nx 01269 s%dsc(i,j)=((s%xv(i+1,j)-s%xv(i,j))**2+(s%yv(i+1,j)-s%yv(i,j))**2)**(0.5d0) 01270 enddo 01271 s%dsc(s%nx+1,j)=s%dsc(s%nx,j) 01272 enddo 01273 else 01274 s%dsc=s%dsu 01275 endif 01276 01277 ! s%dnu 01278 01279 if (s%ny>0) then 01280 do j=2,s%ny+1 01281 do i=1,s%nx+1 01282 s%dnu(i,j)=((xc(i,j)-xc(i,j-1))**2+(yc(i,j)-yc(i,j-1))**2)**(0.5d0) 01283 enddo 01284 enddo 01285 s%dnu(:,1)=s%dnu(:,2) 01286 else 01287 s%dnu=100.d0 01288 endif 01289 01290 ! s%dnz 01291 01292 if (s%ny>0) then 01293 do j=2,s%ny+1 01294 do i=1,s%nx+1 01295 s%dnz(i,j)=((s%xv(i,j)-s%xv(i,j-1))**2+(s%yv(i,j)-s%yv(i,j-1))**2)**(0.5d0) 01296 enddo 01297 enddo 01298 s%dnz(:,1)=s%dnz(:,2) 01299 else 01300 s%dnz=100.d0 01301 endif 01302 01303 ! s%dnv 01304 01305 if (s%ny>0) then 01306 do j=1,s%ny 01307 do i=1,s%nx+1 01308 s%dnv(i,j)=((s%xz(i,j+1)-s%xz(i,j))**2+(s%yz(i,j+1)-s%yz(i,j))**2)**(0.5d0) 01309 enddo 01310 enddo 01311 s%dnv(:,s%ny+1)=s%dnv(:,s%ny) 01312 else 01313 s%dnv=100.d0 01314 endif 01315 01316 ! s%dnc 01317 01318 if (s%ny>0) then 01319 do j=1,s%ny 01320 do i=1,s%nx+1 01321 s%dnc(i,j)=((s%xu(i,j+1)-s%xu(i,j))**2+(s%yu(i,j+1)-s%yu(i,j))**2)**(0.5d0) 01322 enddo 01323 enddo 01324 s%dnc(:,s%ny+1)=s%dnc(:,s%ny) 01325 else 01326 s%dnc=100.d0 01327 endif 01328 01329 01330 if (s%ny>0) then 01331 01332 ! dsdnu 01333 01334 do j=2,s%ny+1 01335 do i=1,s%nx 01336 x1=s%xv(i ,j ) - s%xv(i ,j-1) 01337 x3=s%xv(i+1,j-1) - s%xv(i ,j-1) 01338 x2=s%xv(i+1,j ) - s%xv(i+1,j-1) 01339 x4=s%xv(i+1,j ) - s%xv(i ,j ) 01340 y1=s%yv(i ,j ) - s%yv(i ,j-1) 01341 y3=s%yv(i+1,j-1) - s%yv(i ,j-1) 01342 y2=s%yv(i+1,j ) - s%yv(i+1,j-1) 01343 y4=s%yv(i+1,j ) - s%yv(i ,j ) 01344 dsdnu=0.5d0*(abs(x1*y3-x3*y1)+abs(x2*y4-x4*y2)) 01345 s%dsdnui(i,j)=1.d0/dsdnu 01346 enddo 01347 enddo 01348 s%dsdnui(:,1)=s%dsdnui(:,2) 01349 s%dsdnui(s%nx+1,:)=s%dsdnui(s%nx,:) 01350 01351 ! dsdnv 01352 01353 do j=1,s%ny 01354 do i=2,s%nx+1 01355 x1=s%xu(i-1,j+1) - s%xu(i-1,j ) 01356 x3=s%xu(i ,j ) - s%xu(i-1,j ) 01357 x2=s%xu(i ,j+1) - s%xu(i ,j ) 01358 x4=s%xu(i ,j+1) - s%xu(i-1,j+1) 01359 y1=s%yu(i-1,j+1) - s%yu(i-1,j ) 01360 y3=s%yu(i ,j ) - s%yu(i-1,j ) 01361 y2=s%yu(i ,j+1) - s%yu(i ,j ) 01362 y4=s%yu(i ,j+1) - s%yu(i-1,j+1) 01363 dsdnv=0.5d0*(abs(x1*y3-x3*y1)+abs(x2*y4-x4*y2)) 01364 s%dsdnvi(i,j)=1.d0/dsdnv 01365 enddo 01366 enddo 01367 s%dsdnvi(:,s%ny+1)=s%dsdnvi(:,s%ny) 01368 s%dsdnvi(1,:)=s%dsdnvi(2,:) 01369 01370 ! dsdnz 01371 01372 do j=2,s%ny+1 01373 do i=2,s%nx+1 01374 x1=xc(i-1,j ) - xc(i-1,j-1) 01375 x3=xc(i ,j-1) - xc(i-1,j-1) 01376 x2=xc(i ,j ) - xc(i ,j-1) 01377 x4=xc(i ,j ) - xc(i-1,j ) 01378 y1=yc(i-1,j ) - yc(i-1,j-1) 01379 y3=yc(i ,j-1) - yc(i-1,j-1) 01380 y2=yc(i ,j ) - yc(i ,j-1) 01381 y4=yc(i ,j ) - yc(i-1,j ) 01382 dsdnz=0.5d0*(abs(x1*y3-x3*y1)+abs(x2*y4-x4*y2)) 01383 s%dsdnzi(i,j)=1.d0/dsdnz 01384 enddo 01385 enddo 01386 s%dsdnzi(:,1)=s%dsdnzi(:,2) 01387 s%dsdnzi(1,:)=s%dsdnzi(2,:) 01388 01389 else 01390 01391 s%dsdnui=1.d0/(s%dsu*s%dnu) 01392 s%dsdnvi=1.d0/(s%dsv*s%dnv) 01393 s%dsdnzi=1.d0/(s%dsz*s%dnz) 01394 01395 endif 01396 01397 ! s%alfaz, grid orientation in z-points 01398 01399 do j=1,s%ny+1 01400 do i=2,s%nx 01401 s%alfaz(i,j)=atan2(s%yz(i+1,j)-s%yz(i-1,j),s%xz(i+1,j)-s%xz(i-1,j)) 01402 enddo 01403 s%alfaz(1,j)=s%alfaz(2,j) 01404 s%alfaz(s%nx+1,j)=s%alfaz(s%nx,j) 01405 enddo 01406 01407 ! s%alfau, grid orientation in u-points 01408 01409 do j=1,s%ny+1 01410 do i=1,s%nx 01411 s%alfau(i,j)=atan2(s%yz(i+1,j)-s%yz(i,j),s%xz(i+1,j)-s%xz(i,j)) 01412 enddo 01413 s%alfau(s%nx+1,j)=s%alfau(s%nx,j) 01414 enddo 01415 01416 ! s%alfav, grid orientation in v-points 01417 01418 if (s%ny>0) then 01419 do i=1,s%nx+1 01420 do j=1,s%ny 01421 s%alfav(i,j)=atan2(s%yz(i,j+1)-s%yz(i,j),s%xz(i,j+1)-s%xz(i,j)) 01422 enddo 01423 s%alfav(i,s%ny+1)=s%alfav(i,s%ny) 01424 enddo 01425 else 01426 s%alfav=s%alfaz 01427 endif 01428 01429 do j=1,s%ny+1 01430 s%sdist(1,j)=0 01431 do i=2,s%nx+1 01432 s%sdist(i,j)=s%sdist(i-1,j)+s%dsu(i-1,j) 01433 enddo 01434 enddo 01435 01436 do i=1,s%nx+1 01437 s%ndist(i,1)=0 01438 do j=2,s%ny+1 01439 s%ndist(i,j)=s%ndist(i,j-1)+s%dnv(i,j-1) 01440 enddo 01441 enddo 01442 01443 deallocate (xc) 01444 deallocate (yc) 01445 01446 end subroutine gridprops 01447 01448 subroutine ranges_init(s) 01449 use xmpi_module 01450 implicit none 01451 type(spacepars) :: s 01452 01453 imin_ee = 2 01454 imax_ee = s%nx 01455 if (s%ny>0) then 01456 jmin_ee = 2 01457 jmax_ee = s%ny 01458 else 01459 jmin_ee = 1 01460 jmax_ee = 1 01461 endif 01462 01463 imin_uu = 2 01464 imax_uu = s%nx-1 01465 if (s%ny>0) then 01466 jmin_uu = 2 01467 jmax_uu = s%ny 01468 else 01469 jmin_uu = 1 01470 jmax_uu = 1 01471 endif 01472 01473 imin_vv = 2 01474 imax_vv = s%nx 01475 if (s%ny>0) then 01476 jmin_vv = 2 01477 jmax_vv = s%ny-1 01478 else 01479 jmin_vv = 1 01480 jmax_vv = 1 01481 endif 01482 01483 imin_zs = 2 01484 imax_zs = s%nx 01485 if (s%ny>0) then 01486 jmin_zs = 2 01487 jmax_zs = s%ny 01488 else 01489 jmin_zs = 1 01490 jmax_zs = 1 01491 endif 01492 #ifdef USEMPI 01493 01494 if (.not. xmpi_istop) then 01495 imin_ee = 3 01496 !imin_uu = 2 ! in combination with shift_ee 01497 imin_uu = 3 01498 imin_vv = 3 01499 imin_zs = 3 01500 endif 01501 01502 if (.not. xmpi_isbot) then 01503 imax_ee = s%nx-1 01504 !imax_uu = s%nx-2 01505 imax_uu = s%nx-1 ! in combination with shift_ee 01506 imax_vv = s%nx-1 01507 imax_zs = s%nx-1 01508 endif 01509 01510 if (.not. xmpi_isleft) then 01511 jmin_ee = 3 01512 jmin_uu = 3 01513 ! jmin_vv = 2 01514 jmin_vv = 3 ! in combination with shift_ee 01515 jmin_zs = 3 01516 endif 01517 01518 if (.not. xmpi_isright) then 01519 jmax_ee = s%ny-1 01520 jmax_uu = s%ny-1 01521 ! jmax_vv = s%ny-2 01522 jmax_vv = s%ny-1 ! in combination with shift_ee 01523 jmax_zs = s%ny-1 01524 endif 01525 01526 #endif 01527 end subroutine ranges_init 01528 01529 end module spaceparams
1.7.6.1
