1 subroutine sonfath (ncorp, vnumnpl, ifath, nsons, tfath, 2 & numpe, numnp, maxnode) 3 4c This program writes to geom.dat.<1..numpe> the number sons and 5c fathers required for averaging in the z-direction (Turbulent 6c Boundary Layer), (x,y)-space (Turbulent Channel Flow), 7c or (x,y,z)-space (Isotropic Turbulence) 8 9 integer inform,n1m,n2m,n3m,numpe,numpe1,numnp,maxnode, tfath, 10 & counter,nc8,nc4,node,ifil,itmp 11 integer gfath(0:numnp),ifath(numpe,maxnode),ncorp(numpe,maxnode), 12 & nsons(tfath),dummy(numnp),iper(numnp), vnumnpl(numpe) 13 integer, allocatable :: imap(:),invmap(:) 14 allocate (imap(numnp)) 15 allocate (invmap(numnp)) 16 17 write(*,*)'Enter number of nodes in each direction' 18 read(*,*)n1m,n2m,n3m 19 20 write(*,*)'For Z averaging enter 1 For X-Y averaging enter 2' 21 write(*,*)'For X-Y-Z averaging enter 3 W=0 SPEBC enter 4' 22 write(*,*)'fully local model 5' 23 read(*,*)inform 24 25c 26c.... read and generate map 27c 28 kmap = 15 29 30 open(kmap,file='geom.kmap',status='old',err=183) 31 32 do i=1,numnp 33 read(kmap,*,err=183) iorig 34 imap(i)=iorig+1 ! imap takes arguement of original and returns new 35 invmap(iorig+1)=i ! invmap takes arguemnt of new and returns original 36 enddo 37 goto 184 38 183 write(*,*) 'geom.kmap not read properly...assuming inform=5' 39 inform=5 40 184 continue 41 42 nsoncor=1 43 if(inform.eq.4) then 44 inform=1 45 nsoncor=0 46 endif 47 48 gfath(:) = 0 49 50 nc4=n3m*(n2m-1)+1 51 nc8=n3m*n2m*n1m+1-n3m 52 53 if (numpe.eq.1)then 54 numnp = n1m*n2m*n3m 55 maxnode = numnp 56 endif 57 58 if (inform.eq.1) then ! For averaging over Z - lines 59 60 do i = 1, tfath 61 nsons(i) = n3m-nsoncor 62 enddo 63 64 counter = 1 65 do i = 1, nc8, n3m 66 do j = i, i+n3m-1 67 gfath(j) = counter 68 enddo 69 counter = counter + 1 70 enddo 71 72 if(numpe.gt.1)then 73 do j = 1, numpe 74 do i = 1, maxnode 75 ifath(j,i) = gfath(imap(ncorp(j,i))) 76 enddo 77 enddo 78 else 79 do i = 1, maxnode 80 ifath(1,i) = gfath(imap(i)) 81 enddo 82 endif 83 84 endif 85 86 if (inform.eq.2) then ! For averaging over X-Z planes 87 88 do i = 1, tfath 89 nsons(i) = n1m*n3m - (n1m+n3m-1) 90 enddo 91 92 counter = 1 93 do j = 1, nc4, n3m 94 do i = 1, n1m 95 do k = j, j+n3m-1 96 node = (i-1)*n2m*n3m + k 97 gfath(node) = counter 98 enddo 99 enddo 100 counter = counter + 1 101 enddo 102 103 if(numpe.gt.1)then 104 do j = 1, numpe 105 do i = 1, maxnode 106 ifath(j,i) = gfath(imap(ncorp(j,i))) 107 enddo 108 enddo 109 else 110 do i = 1, maxnode 111 ifath(1,i) = gfath(imap(i)) 112 enddo 113 endif 114 115 endif 116 117 if( inform.eq.3 )then !For averaging over X-Y-Z 118 119 open (unit=96,FILE='geom.iper',form='unformatted', 120 & status='old') 121 read(96) (iper(l),l=1,numnp) 122 close(96) 123 124 dummy(:) = 1 125 126 do i = 1, numnp 127 if(iper(i).eq.0)then 128 iper(i)=i 129 endif 130 enddo 131 do j = 1, numnp 132 i = iper(j) 133 if (i .ne. j) then 134 dummy(j) = 0 135 endif 136 enddo 137 itmp=sum(dummy) 138 nsons(tfath) = itmp 139 write(*,*)'nsons(1)=',nsons(1) 140 ifath=1 141 endif 142 143 if(inform.eq.5) then ! local clipping 144 145 do i = 1, numpe 146 ifil = 99 147 148 nsons=1 149 do j=1,vnumnpl(i) 150 ifath(i,j)=j 151 enddo 152 153 enddo 154 endif 155 156 return 157 end 158