1 subroutine genini (iBC, BC, y, ac, iper, ilwork, 2 & ifath, velbar, nsons,x, 3 & shp, shgl, shpb, shglb ) 4c 5c---------------------------------------------------------------------- 6c This routine reads the initial values in primitive form (density, 7c velocity and temperature), satisfies the boundary conditions and 8c converts them to Y-variables. 9c 10c input: 11c iBC (nshg) : boundary condition code 12c BC (nshg,ndofBC) : boundary condition constrain data 13c x (numnp,nsd) : locations of nodes, numnp-> # of node 14c nsd-> space dimension, 1=x, 2=y, 3=z 15C 16c 17c output: 18c y (nshg,ndof) : initial values of Y variables 19c 20c 21c Farzin Shakib, Winter 1986. 22c Zdenek Johan, Winter 1991. (Fortran 90) 23c---------------------------------------------------------------------- 24c 25 use specialBC ! gets itvn from here 26 use convolImpFlow ! brings in ntimeptpT and other variables 27 use convolRCRFlow ! brings RCR variables 28 include "common.h" 29 include "mpif.h" 30 include "auxmpi.h" 31c 32 dimension iBC(nshg), iper(nshg), 33 & BC(nshg,ndofBC), y(nshg,ndof), 34 & ac(nshg,ndof), x(numnp,nsd), 35 & shp(MAXTOP,maxsh,MAXQPT), 36 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 37 & shpb(MAXTOP,maxsh,MAXQPT), 38 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 39 40c 41 dimension ilwork(nlwork) 42c 43c stuff for dynamic model s.w.avg and wall model 44c 45 dimension ifath(numnp), velbar(nfath,nflow), 46 & nsons(nfath) 47 48 character*20 fname1 49 character*10 cname2 50 character*5 cname 51c 52c.... --------------------------> Restart <--------------------------- 53c 54c.... read q from [RESTAR.INP], reset LSTEP 55c 56 call restar ('in ', y, ac) 57c 58 if((itwmod.gt.0) 59 & .or. (nsonmax.eq.1 .and. iLES.gt.0) ) then 60 call rwvelb('in ',velbar,ifail) 61c 62c if the read failed calculate velbar 63c 64 if(ifail.eq.1) then 65 call getvel (y, ilwork, iBC, 66 & nsons, ifath, velbar) 67 endif 68 69 endif ! for the itwmod or irscale 70c 71c.... time varying boundary conditions as set from file bct.dat and impt.dat 72c (see function for format in file bctint.f) 73c 74 if (itvn .gt. 0 ) then !for inlet velocities 75 call initBCt( x, iBC, BC) 76 call BCint(lstep*Delt(1),shp,shgl,shpb,shglb,x, BC, iBC) 77 endif 78 if (impfile .gt. 0 ) then !for impedance BC 79 if(myrank.eq.master) then 80 write(*,*) 'reading Qhistor.dat' 81 endif 82 open(unit=816, file='Qhistor.dat',status='old') 83 read (816,*) ntimeptpT 84 allocate (QHistImp(ntimeptpT+1,numImpSrfs)) 85 allocate (QHistTry(ntimeptpT,numImpSrfs)) 86 allocate (QHistTryF(ntimeptpT,numImpSrfs)) 87 do j=1,ntimeptpT+1 88 read(816,*) (QHistImp(j,n),n=1,numImpSrfs) !read flow history 89 enddo 90 close(816) 91 call initImpt() !read impedance data and initialize begin/end values 92 do i=2,ntimeptpT 93 call Impint((ntimeptpT-i+1)*Delt(1),i) !return Imp values in reverse order ZN->Z0 94 enddo 95 96 allocate (poldImp(0:MAXSURF)) !for pressure part that depends on the history only 97 endif 98 if (ircrfile .gt. 0 ) then !for RCR BC 99 call initRCRt() !read RCR data 100 dtRCR(:) = Delt(1)/(ValueListRCR(2,:)*ValueListRCR(3,:)) 101 allocate (RCRConvCoef(nstep(1)+lstep+2,numRCRSrfs)) !for convolution coeff 102 !last one needed only to have array of same size as imp BC 103 allocate (QHistRCR(nstep(1)+lstep+1,numRCRSrfs)) !for flow history 104 QHistRCR = zero 105 if(lstep.ne.0) then 106 fname1 = 'Qhistor' 107 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 108 fname1 = trim(fname1) 109 if(myrank.eq.master) then 110 write(*,*) 'reading ', fname1 111 endif 112 open(unit=816, file=fname1, status='old') 113 read(816,*) k 114 do j=1,lstep+1 115 read(816,*) (QHistRCR(j,n),n=1,numRCRSrfs) !read flow history 116 enddo 117 close(816) 118 endif 119 allocate (poldRCR(0:MAXSURF)) !for pressure part that depends on the history only 120 allocate (HopRCR(0:MAXSURF)) !for H operator contribution 121 endif 122 123 124ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 125c this subroutine initBCprofileScale is called to generate the correct 126c BC array, including the siuation of Take BC from IC for inlet. 127c so this subroutine must be called before BCs are applied to ICs 128c as those following lines do 129c call INIBCprofile(BC,y,x) 130c call MPI_BARRIER(MPI_COMM_WORLD,ierr) 131ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 132c 133c.... satisfy the boundary conditions 134c 135 call itrBC (y, ac, iBC, BC, iper, ilwork) 136 137 itempr=mod(impl(1),2) ! tempr solve if impl odd 138 if(itempr.eq.1) then 139 isclr=0 140 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 141 endif 142 do isclr=1,nsclr 143 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 144 enddo 145 146 if((irscale.ge.0) .and. (myrank.eq.master)) then 147 call genscale(y, x, iBC) 148 endif 149c 150c.... ---------------------------> Echo <---------------------------- 151c 152c.... echo the initial data 153c 154 if ((necho .lt. 0).and.(myrank.eq.master)) then 155 do n = 1, nshg 156 if (mod(n,50) .eq. 1) write(iecho,1000) ititle,(i,i=1,ndof) 157 write (iecho,1100) n, (y(n,i),i=1,ndof) 158 enddo 159 endif 160c 161c.... return 162c 163 return 164c 1651000 format(a80,//, 166 & ' I n i t i a l V a l u e s ',//, 167 & ' Node ',/, 168 & ' Number ',6x,6('dof',i1,:,10x)) 1691100 format(1p,2x,i5,5x,5(e12.5,2x)) 170c 171 end 172