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) 57 58ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 59 if(matflg(1,1).eq.0)then ! compressible code 60 call INIprofile(BC,y,x) 61 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 62 endif 63cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 64c 65 if((itwmod.gt.0) 66 & .or. (nsonmax.eq.1 .and. iLES.gt.0) ) then 67 call rwvelb('in ',velbar,ifail) 68c 69c if the read failed calculate velbar 70c 71 if(ifail.eq.1) then 72 call getvel (y, ilwork, iBC, 73 & nsons, ifath, velbar) 74 endif 75 76 endif ! for the itwmod or irscale 77c 78c.... time varying boundary conditions as set from file bct.dat and impt.dat 79c (see function for format in file bctint.f) 80c 81 if (itvn .gt. 0 ) then !for inlet velocities 82 call initBCt( x, iBC, BC) 83 call BCint(lstep*Delt(1),shp,shgl,shpb,shglb,x, BC, iBC) 84 endif 85 if (impfile .gt. 0 ) then !for impedance BC 86 if(myrank.eq.master) then 87 write(*,*) 'reading Qhistor.dat' 88 endif 89 open(unit=816, file='Qhistor.dat',status='old') 90 read (816,*) ntimeptpT 91 allocate (QHistImp(ntimeptpT+1,numImpSrfs)) 92 allocate (QHistTry(ntimeptpT,numImpSrfs)) 93 allocate (QHistTryF(ntimeptpT,numImpSrfs)) 94 do j=1,ntimeptpT+1 95 read(816,*) (QHistImp(j,n),n=1,numImpSrfs) !read flow history 96 enddo 97 close(816) 98 call initImpt() !read impedance data and initialize begin/end values 99 do i=2,ntimeptpT 100 call Impint((ntimeptpT-i+1)*Delt(1),i) !return Imp values in reverse order ZN->Z0 101 enddo 102 103 allocate (poldImp(0:MAXSURF)) !for pressure part that depends on the history only 104 endif 105 if (ircrfile .gt. 0 ) then !for RCR BC 106 call initRCRt() !read RCR data 107 dtRCR(:) = Delt(1)/(ValueListRCR(2,:)*ValueListRCR(3,:)) 108 allocate (RCRConvCoef(nstep(1)+lstep+2,numRCRSrfs)) !for convolution coeff 109 !last one needed only to have array of same size as imp BC 110 allocate (QHistRCR(nstep(1)+lstep+1,numRCRSrfs)) !for flow history 111 QHistRCR = zero 112 if(lstep.ne.0) then 113 fname1 = 'Qhistor' 114 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 115 fname1 = trim(fname1) 116 if(myrank.eq.master) then 117 write(*,*) 'reading ', fname1 118 endif 119 open(unit=816, file=fname1, status='old') 120 read(816,*) k 121 do j=1,lstep+1 122 read(816,*) (QHistRCR(j,n),n=1,numRCRSrfs) !read flow history 123 enddo 124 close(816) 125 endif 126 allocate (poldRCR(0:MAXSURF)) !for pressure part that depends on the history only 127 allocate (HopRCR(0:MAXSURF)) !for H operator contribution 128 endif 129 130 131ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 132c this subroutine initBCprofileScale is called to generate the correct 133c BC array, including the siuation of Take BC from IC for inlet. 134c so this subroutine must be called before BCs are applied to ICs 135c as those following lines do 136c call INIBCprofile(BC,y,x) 137c call MPI_BARRIER(MPI_COMM_WORLD,ierr) 138ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 139c 140c.... satisfy the boundary conditions 141c 142 call itrBC (y, ac, iBC, BC, iper, ilwork) 143 144 itempr=mod(impl(1),2) ! tempr solve if impl odd 145 if(itempr.eq.1) then 146 isclr=0 147 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 148 endif 149 do isclr=1,nsclr 150 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 151 enddo 152 153 if((irscale.ge.0) .and. (myrank.eq.master)) then 154 call genscale(y, x, iBC) 155 endif 156c 157c.... ---------------------------> Echo <---------------------------- 158c 159c.... echo the initial data 160c 161 if ((necho .lt. 0).and.(myrank.eq.master)) then 162 do n = 1, nshg 163 if (mod(n,50) .eq. 1) write(iecho,1000) ititle,(i,i=1,ndof) 164 write (iecho,1100) n, (y(n,i),i=1,ndof) 165 enddo 166 endif 167c 168c.... return 169c 170 return 171c 1721000 format(a80,//, 173 & ' I n i t i a l V a l u e s ',//, 174 & ' Node ',/, 175 & ' Number ',6x,6('dof',i1,:,10x)) 1761100 format(1p,2x,i5,5x,5(e12.5,2x)) 177c 178 end 179