! VIPER_IGE - modeled after IMPATH() !23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 ! ! ROUTINE TO COMPUTE THE TRAJECTORIES OF A COLLECTION OF POINT VORTICES ! ! modified to reorder YZ data; first dimension set to two. C SUBROUTINE VIPER_IGE(T2,DTI,BZ,VZIN,ACSPEED,i_hwa,NTOUT, & TB,YPB,ZPB,GPB,YSB,ZSB,GSB) PARAMETER (MAXSAV=20000,NPTMAX=160,NYZMAX=2*NPTMAX,NXIMAX=721) PARAMETER (NSTPMX=5001,NVARBL=10,NZMAX=401,NNUT=8) C REAL TB(*), YPB(*), ZPB(*), GPB(*), YSB(*), ZSB(*), GSB(*) C c Tom - The air properties gamma_gas and R_gas have been added to common block c cvalb. COMMON /CVALB/GRAV,TZ,ADLAPS,C1,C2,C3,ZZ,ZIM,ZGE,YOVER,ZDOWN, & gamma_gas, R_gas, dis COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) X ,QZ(NZMAX),Q(NZMAX) X ,UCRZ(NZMAX),UCR(NZMAX) X ,UHWZ(NZMAX),UHW(NZMAX) & ,PRESS(NZMAX),THETA(NZMAX),BVFSQ(NZMAX),RHO(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /ATMB2/DUCR(NZMAX),DDUCR(NZMAX) COMMON /PATH1/DTSAV,TP(MAXSAV),KMAX,KOUNT COMMON /PATH2/YZP(2,NPTMAX,MAXSAV) COMMON /NGAM/GAM(NPTMAX),CGAM(NPTMAX),NPTS COMMON /GAMPB/GAMP(NPTMAX,MAXSAV),KG COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC COMMON /IZTOPB/ZLAST,QB,QB1,QB4 & ,IZTOPT,IZTOPQ,IZTOPU,IZTOPHW,NUTTOP COMMON /IZTPRB/IZTOPR COMMON /QB58B/QB5,QB8,BLAST,IZ2PDU COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3,i_hw common /ige/ DTFAC,EPS,H1,HMIN,CEJECTA,CKE,R0p,R0s, & TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF,ZGF, & KMAXA REAL YZ(2,NPTMAX) REAL R0p,R0s,CEJECTA,CKE,TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF C REAL PPY(MAXSAV),PPZ(MAXSAV),PPG(MAXSAV) REAL PPYI(NXIMAX),PPZI(NXIMAX),PPGI(NXIMAX) INTEGER NOK, NBAD, i_hw, i_hwa COMMON /GTOGB/YZPI(2,NPTMAX,NXIMAX),GAMI(NPTMAX,NXIMAX) & ,XI(NXIMAX),NXI C ! logical debug1, debug2, debug4, debug6, debug_ige common /debug/ debug1, debug2, debug4, debug6, debug_ige NOK = 0 NBAD = 0 tan3 = tan(3.0*pi/180.0) if( acspeed .eq. 0.0 ) then acspeedi = 0.0 else acspeedi = 1.0/acspeed end if i_hw = i_hwa c h1 = dti C FILL /THB/ PI = 4.*ATAN(1.) TWOPI = 2.*PI R2PI = 1./TWOPI PIS4 = .25*PI PIS3 = PI/3. PIS2 = .5*PI PI3S2 = 1.5*PI THBOT = PIS4 THTOP = TWOPI + PIS4 THMIN1 = PI + PIS4 THMIN2 = THMIN1 + PIS2 THMAX1 = PIS4 THMAX2 = PIS4 DELTH1 = THMIN1 - THMAX2 DELTH2 = THTOP - THMIN2 FACMIN = 0. FACMAX = 1. DELFAC = FACMAX - FACMIN ! IZTOPT = NZTDEG IZTOPQ = NZQ IZTOPU = NZUCR IZ2PDU = NZUCR IZTOPHW = NZUHW NUTTOP = NNUT ZLAST = ZZ BLAST = BZ QB1 = 0. QB4 = 0. QB5 = 0. QB8 = 0. ! IF(NZUCR .GT. 2) THEN NZUCRM = NZUCR - 1 DO IZ=2,NZUCRM DUCR(IZ) = (UCR(IZ+1)-UCR(IZ-1)) X /(UCRZ(IZ+1)-UCRZ(IZ-1)) DDUCR(IZ) = 4.*(UCR(IZ+1)+UCR(IZ-1) X -2.*UCR(IZ))/ X (UCRZ(IZ+1)-UCRZ(IZ-1))**2 ENDDO DUCR(1) = (UCR(2)-UCR(1))/(UCRZ(2)-UCRZ(1)) DUCR(NZUCR) = (UCR(NZUCR)-UCR(NZUCRM)) X /(UCRZ(NZUCR)-UCRZ(NZUCRM)) DDUCR(1) = DDUCR(2) DDUCR(NZUCR) = DDUCR(NZUCRM) ELSE DUCR(1) = (UCR(2)-UCR(1))/(UCRZ(2)-UCRZ(1)) DUCR(2) = DUCR(1) DDUCR(1) = 0. DDUCR(2) = 0. ENDIF KFLAG = 0 KUCR = 0 ! T1 = 0.0 KMAX = KMAXA NPTS = 2 GAM( 1 ) = GPB(NTOUT) GAMP(1,1) = GPB(NTOUT) GAM( 2 ) = GSB(NTOUT) GAMP(2,1) = GSB(NTOUT) YZ( 1,1) = YPB(NTOUT) YZ( 1,2) = YSB(NTOUT) YZ( 2,1) = ZPB(NTOUT) YZ( 2,2) = ZSB(NTOUT) BZ2 = sqrt( (YSB(NTOUT)-YPB(NTOUT))**2 + & (ZSB(NTOUT)-ZPB(NTOUT))**2 ) ZGE = ZGF*BZ2 ! UNLIKE VIPER OGE, VIPER IGE EVOLVES THE SIGNED CIRCULATIONS. C ! if ZIM and ZGE are the same altitude, go right to ZGE ! (and get out of this routine) IF(ABS(ZIM-ZGE) .LT. .001*ZIM) THEN KOUNT = 1 YZP(1,1,1) = YZ(1,1) YZP(2,1,1) = YZ(2,1) YZP(1,2,1) = YZ(1,2) YZP(2,2,1) = YZ(2,2) TLAST = 0. RETURN ENDIF C DTSAV=DTFAC*(T2-T1)/FLOAT(MAXSAV) CC PRINT 9999,(YZ(1,IY),IY=1,8),(YZ(2,IZ),IZ=1,8),(GAM(IG),IG=1,8) CC 9999 FORMAT('Y = ',8E13.4/'Z = ',8E13.4/'GAM = ',8E13.4) KG = 1 ! ! advance vortex field due to mutual advection ! use NPTSA to pass the number of vortices initially, because ! NPTS appears in a common block in IMIN, and that could lead ! to some confusion. ! c *** Not sure why, but imint does not reliably solve out to time t2. c *** Thus we add five additional dt. The solution is interpolated c *** onto a different time grid below, so doing an extra few steps c *** is of no lasting consequence. NPTSA = NPTS H1 = 0.5*dti CALL IMINT(YZ,NPTSA,T1,T2+5*dti,EPS,H1,HMIN,NOK,NBAD & ,CEJECTA,CKE,R0p,R0s,BZ,VZIN & ,TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF) c write(6,*)"YZ(1,1)=",YZ(1,1) c write(6,*)"YZ(2,1)=",YZ(2,1) c write(6,*)"YZ(1,2)=",YZ(1,2) c write(6,*)"YZ(2,2)=",YZ(2,2) c write(6,*)"NPTSA=",NPTSA c write(6,*)"T1=",T1 c write(6,*)"T2=",T2 c write(6,*)"EPS=",EPS c write(6,*)"H1=",H1 c write(6,*)"HMIN=",HMIN c write(6,*)"NOK=",NOK c write(6,*)"NBAD=",NBAD c write(6,*)"KG=",KG c write(6,*)"KOUNT=",KOUNT KMINI = MIN(KG,KOUNT) KG=KMINI-1 KOUNT=KMINI-1 IF(KG.NE.KOUNT) THEN WRITE(2,999) KG,KOUNT 999 FORMAT(/1X,'KG.NE.KOUNT'/1X,'KG,KOUNT = ',2I10) CLOSE(2) CCC CLOSE(3) STOP ENDIF C c print *, 'kount, npts, dt = ', kount, npts, tp(kount)-tp(kount-1) if( debug_ige ) then write(6,*)"KOUNT=",KOUNT end if c *** The procedure for sampling the solution is rediculous. I needed to c *** add one step to nxi to make things work. c print *, 'kount, tlast, t2 = ', kount, TP(KOUNT), t2 TLAST = min( TP(KOUNT), t2 ) c nxi = int((tlast-t1)/dti+1.0e-4)+1 nxi = int((tlast-t1)/dti+1.0e-4) + 1 dxi = dti DO 100 I=1,NXI 100 XI(I) = DXI*FLOAT(I-1) C CCC DO 103 J=1,KOUNT CCC WRITE(3,131) TP(J),(GAMP(I,J),YZP(1,I,J),YZP(2,I,J),I=1,2) CCC 103 CONTINUE C !JW DO 120 IPT=1,NPTSA DO 120 IPT=1,NPTS C DO 105 J=1,KOUNT PPG(J) = GAMP(IPT,J) PPY(J) = YZP(1,IPT,J) PPZ(J) = YZP(2,IPT,J) 105 CONTINUE C CALL GTOGL(KOUNT,TP,PPG,NXI,XI,PPGI) CALL GTOGL(KOUNT,TP,PPY,NXI,XI,PPYI) CALL GTOGL(KOUNT,TP,PPZ,NXI,XI,PPZI) C DO 110 J=1,NXI GAMI( IPT,J) = PPGI(J) YZPI(1,IPT,J) = PPYI(J) YZPI(2,IPT,J) = PPZI(J) 110 CONTINUE 120 CONTINUE C CC PRINT 121 CC 121 FORMAT(/////) DO 125 J=1,NXI 125 XI(J) = XI(J) + TB(NTOUT) ! write(6,*)"NTOUT,NXI=",NTOUT,NXI DO 130 J=1,NXI JJ = NTOUT + J CCC WRITE(3,131) XI(J),(GAMI(I,J),YZPI(1,I,J),YZPI(2,I,J),I=1,2) TB( JJ) = XI(J) YPB(JJ) = YZPI(1,1,J) ZPB(JJ) = YZPI(2,1,J) GPB(JJ) = GAMI( 1,J) YSB(JJ) = YZPI(1,2,J) ZSB(JJ) = YZPI(2,2,J) GSB(JJ) = GAMI( 2,J) 130 CONTINUE NTOUT = NTOUT + NXI if( debug2 ) then WRITE(2,*) 'NTOUT = ',NTOUT end if c write(40,40) tb(ntout), npts c40 format(f10.4,i5) CCC 131 FORMAT(F8.3,6F10.3) C CC PRINT 81, NOK,NBAD,KOUNT CC 81 FORMAT(' NOK,NBAD,KOUNT = ',3I10) CC DO 95 J=1,KOUNT CC PRINT 91, TP(J),(YZP(I,J),YZP(NPTS+I,J),I=1,NPTS) CC 91 FORMAT(F8.3,8F9.3) CC 95 CONTINUE C IF(TLAST.GE.T2) TLAST = -TLAST RETURN END C ! IMINT - it is not clear that i needed to modify this routine ! to identify 2 as the first dimension for all the Y variables. ! SUBROUTINE IMINT(YSTART,NPTSA,X1,X2,EPS,H1,HMIN,NOK,NBAD & ,CEJECTA,CKE,R0p,R0s,BZ,VZIN & ,TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF) PARAMETER (MAXSTP=20000,TWO=2.0,ZERO=0.0,TINY=1.E-20) PARAMETER (MAXSAV=20000,NPTMAX=160,NYZMAX=2*NPTMAX,NZMAX=401) PARAMETER (NMAX=NPTMAX) C COMMON /KFLAGB/KUCR,KFLAG C COMMON /PATH1/ DXSAV,XP(MAXSAV),KMAX,KOUNT COMMON /PATH2/ YP(2,NPTMAX,MAXSAV) c Tom - The air properties gamma_gas and R_gas have been added to common block c cvalb. COMMON /CVALB/GRAV,TZ,ADLAPS,C1,C2,C3,ZZ,ZIM,ZGE,YOVER,ZDOWN, & gamma_gas, R_gas, dis COMMON /NGAM/GAM(NPTMAX),CGAM(NPTMAX),NPTS COMMON /DGAMB/DGAM(NPTMAX),DELGAM COMMON /GAMPB/GAMP(NPTMAX,MAXSAV),KG COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) X ,QZ(NZMAX),Q(NZMAX) X ,UCRZ(NZMAX),UCR(NZMAX) X ,UHWZ(NZMAX),UHW(NZMAX) & ,PRESS(NZMAX),THETA(NZMAX),BVFSQ(NZMAX),RHO(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /ATMB2/DUCR(NZMAX),DDUCR(NZMAX) COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3,i_hw DIMENSION YSTART(2,NPTSA),YSCAL(2,NMAX),Y(2,NMAX),DYDX(2,NMAX), X THFAC(NPTMAX),DGAMPRIMEINIT(2),DGAMF(NPTMAX) & ,IGEPHASE(NPTMAX) DIMENSION DGAMC(2,NPTMAX),TIMEC(2,2,NPTMAX),ZMIN(NPTMAX) & ,GAM0(NPTMAX) INTEGER :: P2ND, S2ND, YCFLAG, IGEPHASE, PNEAR(NPTMAX) REAL :: GAMSMALL, DXEJECT1, DXEJECT2, DXEJECTP, DXEJECTS REAL :: RADII(NPTMAX), AREA(NPTMAX), DDUK(NPTMAX) REAL :: YC, R0p, R0s, CEJECTA, CKE REAL :: TEJECT1, TEJECT2, ZEJECT, SECDK, SECDIFF REAL :: DGAMPRIMEINIT,DGAMF,DGAMC,TIMEC REAL :: ZMIN,frac,GAM0 logical debug1, debug2, debug4, debug6, debug_ige common /debug/ debug1, debug2, debug4, debug6, debug_ige ! input variables: ! YSTART - (y,z) for all primary point vortices (no images) ! NPTSA - number of point vortices initially ! X1,X2 - initial and final times for this routine (internal clock) ! EPS - error tolerance for controlling timestep size passed to RKQC() ! H1 - time step (seconds) ! HMIN - minimum time step ! KOUNT - number of time steps completed in this routine ! R0p - radius of port vortex ! R0s - radius of starboard vortex ! DXEJECTP - time between port b.l. eruptions ! DXEJECTS - time between starbd b.l. eruptions ! for 1st eruption DXEJECTP is 1.5*[2.pi.b0.b0/GAMP(1,1)] ! subsequent eruptions it is 0.1*[8.pi*ZP*ZP/GAMP(1,t)] ! ! internal variables: ! Y - primary vortices + boundary-layer-ejection vortices ! MAXSTP - upper limit on time steps to take ! NMAX - number of vortices allowed in internal variables ! MAXSAV - maximum number of times output data is saved to disk ! NPTMAX - maximum number of vortices allowed ! NYZMAX - ????? ! NVAR - 2*NPTS, number of variables advanced, (2 for Y & Z) ! XSAV - latest time ! frac - fraction of secondary vortex energy resulting from primary ! DGAMC - circulation decay resulting from new-vortex creation ! DGAMC(1,i) - port vortex decay rate resulting from ith vortex ! DGAMC(2,i) - starbd vortex decay rate resulting from ith vortex ! TIMEC(1,1,i) - start time for applying DGAMC(1,i) to port vortex ! TIMEC(2,1,i) - end time for applying DGAMC(1,i) to port vortex ! TIMEC(1,2,i) - start time for applying DGAMC(1,i) to starbd vortex ! TIMEC(2,2,i) - end time for applying DGAMC(1,i) to starbd vortex ! IGEPHASE(i) - ige phase (0=oge, 1=nge, 2=ige1, 3=ige2) ! ! output variables: ! NOK - good data points ????? ! NBAD - bad data points ????? ! NPTS=NPTSA NVAR=NPTS*2 !! NPTS=NVAR/2 CC NPAIRS=NPTS/2 X=X1 H=SIGN(H1,X2-X1) NOK=0 NBAD=0 KOUNT=0 if( debug_ige ) then write(6,*)"* * * * * * * * * * * * * * * * * * * *" write(6,*)" UCR(1) = ",UCR(1) write(6,*)"* * * * * * * * * * * * * * * * * * * *" end if ! variables to hold the first port secondary vortex and the ! first starboard secondary vortex. P2ND = 0 S2ND = 0 ! initialize creation-drag circulation decay do ii=1,nptmax PNEAR(ii) = 0 ZMIN(ii) = 0.0 IGEPHASE(ii) = 1 do i=1,2 DGAMC( i,ii) = 0.0 TIMEC(1,i,ii) = 0.0 TIMEC(2,i,ii) = 0.0 enddo enddo ! initialize GAMSMALL GAMSMALL = 3.0/(45.5*2.)*(BZ*VZIN) * 1.0 ! initial and subsequent b.l. eruption build-up time ! DXEJECT1 = 1.3*0 + 0.5 ! DXEJECT2 = 0.02 DXEJECT1 = TEJECT1 DXEJECT2 = TEJECT2 ! initialize port/starboard eruption-time variables XEJECTP = 9999.99 XEJECTS = 9999.99 RADII(1) = R0p RADII(2) = R0s AREA(1) = PI * R0p*R0p AREA(2) = PI * R0s*R0s frac = 0.25 ! this is needed to interpolate DDUK IZ2PDU = NZUCR IZTOPU = NZUCR iztopq = nzq ! store the initial circulation decay DGAMPRIMEINIT(1) = DGAM(1) DGAMPRIMEINIT(2) = DGAM(2) ! ! DIS - kinematic viscosity, normally called nu [m^2/s] ! ! CEJECTA - effective area over which to integrate ejected boundary- ! layer vorticity to determine the initial circulation of an ejected ! vortex. ! ! CKE - geometrical factor for kinetic energy decay due to boundary ! scouring by primary vortex on close approach ! ! XEJECTP - time at which last port secondary vortex was ! ejected from the boundary layer. ! ! XEJECTS - time at which last starboard secondary vortex was ! ejected from the boundary layer. ! ! Tom - CEJECTA and CKE are now set in viper.F and passed into viper_ige. ! ! THFAC(I) is related to the angle made by the secondary vortices ! advected around the primaries. it is initiated to 1.0 for the ! first 4 points, the two primary vortices and their two images. ! however, we don't need it for the two images, since we are now ! treating primary vortices and images together. i've reduced the ! do-loop limit from 4 to 2 to address this. ! ! initialize THFAC for the primary vortices DO 5 I=1,2 5 THFAC(I) = 1. ! initiate the primary vortex positions from YSTART. DO 11 I=1,NPTS Y(1,I)=YSTART(1,I) Y(2,I)=YSTART(2,I) 11 CONTINUE ! XSAV = last time saved (for later output) ! DXSAV = a very small time step below which is too small ! YSCAL = solution scale factor for error evaluation ! NVORT = total number of vortices. initially it is just 2. ! NVORT = NPTS ! NSTP is the time-step number XSAV=X-DXSAV*TWO DO 16 NSTP=1,MAXSTP ! CGAM is used in the IGDRVS() routine. DO 140 I=1,NPTS 140 CGAM(I) = GAM(I) ! correct Z too low do i=1,npts if (y(2,i).lt.0.) then y(2,i) = amax1(y(2,i),zmin(i)) if( debug_ige ) then write(6,*)"i,z(i)=",i,y(2,i) end if endif enddo ! compute time derivatives resulting from multi-vortex advection CALL IGDRVS(Y,DYDX) ! IF((X+H-X2)*(X+H-X1).GT.ZERO) H=X2-X ! ! if new time will exceed X2, adjust time step to finish on X2 IF(X+H.GT.X2) H=X2-X ! YSCAL is used in RKQC() DO 12 I=1,NPTS YSCAL(1,I)=ABS(Y(1,I))+ABS(H*DYDX(1,I))+TINY YSCAL(2,I)=ABS(Y(2,I))+ABS(H*DYDX(2,I))+TINY 12 CONTINUE IF(KMAX.GT.0)THEN ! if X-XSAV > small value IF(ABS(X-XSAV).GE.ABS(DXSAV)) THEN ! and if our counter is less than the last step ... IF(KOUNT.LT.KMAX)THEN ! store output data KOUNT=KOUNT+1 XP(KOUNT)=X DO 13 I=1,NPTS YP(1,I,KOUNT)=Y(1,I) YP(2,I,KOUNT)=Y(2,I) 13 CONTINUE XSAV=X ! endif:KOUNTdxsav ENDIF ! endif:kmax>0 ENDIF ! update solution vectors CALL RKQC(Y,DYDX,NVAR,X,H,EPS,YSCAL,HDID,HNEXT & ,GAM,GAMSMALL,ZMIN) ! if you time advance more than DXSAV past XSAV, explore ! adding a new vortex. ! this is the same cretierion as KOUNT. why? IF(ABS(X-XSAV).GE.ABS(DXSAV)) THEN KG = KG + 1 !JW write(6,*)"NSTP,time,GAM(1),GAM(2)=",NSTP,X,GAM(1),GAM(2) ! you need to check if criteria is met to eject vortices. ! try: if surface-flow-convergence location is persistent. ! how persistent? ! within a fixed lateral distance range of the primary vortex ! for the time required to grow the b.l. by a factor of 5,6,7(???) ! note: as the primary vortices drop in altitude, the eruption ! point migrates towards the stage/player. so another height ! could work as a surrogate criterion on this one. ! ! ! update GAM(I) and store in GAMP(I,KG) ! we should update continuously in RK routine DO 150 I=1,NPTS ! lin-model OGE decay rate: -(0.1020 + 0.687 EDR*)*GAM0*V0/B0 ! lin-model NGE decay rate: -(0.0715 + 0.108 EDR*)*GAM0*V0/B0 ! use this for 2ndary vortices until you have better theory QK = FK( Y(2,I) , QZ, Q, NZQ, iztopq ) ! QSTAR = 0.0 QSTAR = (QK*BZ)**0.333333/VZIN ! get crosswind 2nd derivative at vortex height DDUK(I) = FK( Y(2,I),UCRZ,DDUCR,NZUCR,IZ2PDU ) if (abs(GAM(I)).lt.GAMSMALL) then GAM(I) = 0.0 DGAM(I) = 0.0 else if (I.le.2) then ! linear-model nge decay ! linear-model ige2 decay ! DGAM(I) = - sign( 0.0400 + 0.020*QSTAR , GAM(I) ) ! DGAM(I) = - sign( 0.0472 + 0.400*QSTAR , GAM(I) ) !!! if (IGEPHASE(I) .le. 1) then !!! DGAM(I) = - sign( 0.0715 + 0.108*QSTAR , GAM(I) ) !!! & * VZIN*VZIN*TWOPI !!! else !!! DGAM(I) = - sign( 0.0472 + 0.400*QSTAR , GAM(I) ) !!! & * VZIN*VZIN*TWOPI !!! endif !!! if (Y(2,I) .le. 0.70*BZ) then !!! IGEPHASE(I) = 2 !!! DGAM(I) = - sign( 0.0591 + 0.915*QSTAR , GAM(I) ) !!! & * VZIN*VZIN*TWOPI !!! endif ! exponential decay of primary vortices ! DGAM(I) = - ( 0.1500 + 0.000*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.1200 + 0.000*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.1000 + 0.000*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.2000 + 0.000*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.1500 + 0.000*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.1250 + 0.000*QSTAR )*GAM(I)*VZIN/BZ DGAM(I) = - ( 0.02117+0.6964*QSTAR )*GAM(I)*VZIN/BZ ! DGAM(I) = - ( 0.071 + 0.421*QSTAR )*GAM(I)*VZIN/BZ ! surface-friction drag contribution to circulation decay call gamdotfriction(DGAMF(I),GAM(I),Y(2,I),DYDX(2,I) & ,DIS,CKE,RADII(I)) ! circulation decay for primary vortices ! DGAM(I) = DGAMF(I) + DGAMPRIMEINIT(I) DGAM(I) = DGAMF(I)*1. + DGAM(I) & + 1.0*AREA(I)*DYDX(2,I)*DDUK(I) ! write(6,*)" DGAMF(",I,")=",DGAMF(I) ! ! compute the diffusion interaction between vortex II ! and a primary vortex I. use DGAMC array for this. ! do II=3,NPTS SGII = sign( 1.0 , GAM(II)*GAM(I) ) RGII = (Y(1,II)-Y(1,I))*(Y(1,II)-Y(1,I)) & + (Y(2,II)-Y(2,I))*(Y(2,II)-Y(2,I)) ! ! proximity leads to decay for the secondary vortex ! regardless its sign. the circulation lost diffuses ! with the primary vortex, either increasing it or ! decreasing it. ! ! add this to the secondary vortex and subtract it ! from the primary. ! ! oops! this is actually the location of minimum ! velocity. ! DGAMC(I,II) = - 0.5 * sign(1.0,GAM(II))*( ! & abs(abs(GAM(I))**.333-SGII*abs(GAM(II))**.333)**3. ! & / (TWOPI*RGII) ) ! ! change to diffusion proportional to straining field ! of primary vortex at the site of the secondary vortex ! DGAMC(I,II) = - 0.20 * GAM(II) ! DGAMC(I,II) = - 0.05 * GAM(II) ! DGAMC(I,II) = - 0.00 * GAM(II) DGAMC(I,II) = - SECDIFF * GAM(II) & * abs(GAM(I)) / (TWOPI*RGII) !testing testing ! DGAMC(I,II) = 0.0 DGAM(I) = DGAM(I) - DGAMC(I,II) enddo !!! do II=3,NPTS !! if ( X.ge.TIMEC(1,I,II) .and. !! & X.le.TIMEC(2,I,II) ) then !!! if (PNEAR(II).eq.I) then !!! RR = (Y(1,I)-Y(1,II))**2 + (Y(2,I)-Y(2,II))**2 !!! if (RR.lt.BZ*BZ*1.0*1.0) then !!! DGAM(I) = DGAM(I) + GAM0(II)*GAM(II)*0.1 !! DGAM(I) = DGAM(I) + DGAMC(I,II)*0.0 ! write(6,*)" DGAMC(",I,II,")=",DGAMC(I,II) !!! endif !!! endif !!! enddo else ! circulation decay for ejected secondary vortices ! linear-model nge deay ! DGAM(I) = - sign( 0.0715 + 0.108*QSTAR , GAM(I) ) ! linear-model ige2 decay ! DGAM(I) = - sign( 0.0472 + 0.400*QSTAR , GAM(I) ) ! DGAM(I) = - sign( 0.0400 + 0.020*QSTAR , GAM(I) ) ! run_viper_7 DGAM(I) = - sign( 0.1000 + 0.000*QSTAR , GAM(I) ) ! run_viper_11 DGAM(I) = - sign( 0.1000 + 0.000*QSTAR , GAM(I) ) ! run_viper_8 DGAM(I) = - sign( 0.0100 + 0.000*QSTAR , GAM(I) ) ! run_viper_14* DGAM(I) = - sign( 0.0100 + 0.000*QSTAR , GAM(I) ) ! & * VZIN*VZIN*TWOPI ! DGAM(I) = - ( 0.0300 + 0.100*QSTAR )*GAM(I) ! & * abs(GAM(I))/abs(GAM0(I)) ! testing testing DGAM(I) = -GAM0(I)*GAM(I) * 1.0 DGAM(I) = DGAM(I) + DGAMC(1,I) + DGAMC(2,I) ! surface-friction drag contribution to circulation decay ! call gamdotfriction(DGAMF(I),GAM(I),Y(2,I),DYDX(2,I) ! & ,DIS,CKE,RADII(I)) ! call gamdotfriction(DGAMF(I),GAM(I),Y(2,I),DYDX(2,I) ! & ,DIS,CKE,0.10 ) ! ! DGAM(I) = DGAMF(I)*0. + DGAM(I) ! & + 2.*AREA(I)*DYDX(2,I)*DDUK(I) ! DGAM(I) = + DGAM(I) & + 1.*AREA(I)*DYDX(2,I)*DDUK(I) endif if( i_hw .eq. 1 ) then UHWK = FK( Y(2,I),UHWZ,UHW,NZUHW,IZTOPU) DGAM(I) = DGAM(I) * ( 1 - UHWK*ACSPEEDI ) end if GAMOLD = GAM(I) GAM(I) = GAMOLD + DGAM(I)*(x-xsav) ! GAM(I) = GAMOLD + DGAM(I)*(h) ! if GAM(I) changes sign, set it to zero. if ( GAM(I)*GAMOLD .lt. 0.0 ) then GAM(I) = 0.0 endif endif GAMP(I,KG) = GAM(I) 150 CONTINUE ! if either primary circulation hits zero magnitude ... ! GAM(1) is negative; it is port ! GAM(2) is positive; it is starboard ! IF(GAM(1) .GE. 0. .OR. GAM(2) .LE. 0.) THEN !JWJW IF(GAM(1).GE.0. .OR. GAM(2).LE.0. ) THEN ! ... set all circulation magnitude's to zero and set flag ! ! do we really want to bail on everything if only one of ! the primary vortices peters out? we could still allow ! the other to continue. test this, but then try without. ! !JWJW DO IG=1,NPTS !JWJW GAM(IG) = 0. !JWJW DGAM(IG) = 0. !JWJW GAMP(IG,KG) = 0. !JWJW KFLAG = 1 !JWJW ENDDO ! endif GAM1 or GAM2 decay to zero !JWJW ENDIF ! endif for X-XSAV > DXSAV ENDIF !JOE - check this out IF(HDID.EQ.H)THEN NOK=NOK+1 ELSE NBAD=NBAD+1 ENDIF ZTEST = AMIN1(Y(2,1),Y(2,2)) ! port vortex drops into IGE if ( Y(2,1) .le. ZGE ) then ! begin buildup of first ejected port vortex if ( XEJECTP .eq. 9999.99 ) then XEJECTP = X DXEJECTP = DXEJECT1*2.5*BZ & /ABS(GAMP(1,1)*1./(PI*BZ)+1.*UCR(1)) !!! DXEJECTP = DXEJECT1*4.*BZ !!! & /ABS(GAMP(1,1)*2./(PI*BZ)+0.*UCR(1)) if( debug_ige ) then write(6,*)"time,DXEJECTP = ",X,DXEJECTP end if endif ! finish buildup of new port-side b.l. vortex: at least DXEJECTP ! since last port vortex was ejected; not reached maximum ! vortices; neither primary vortex has decayed away. if ( (X-XEJECTP.GE.ABS(DXEJECTP)) & .AND. NPTS.LT.NPTMAX & .AND. KFLAG.ne.1 ! & .AND. 0.eq.1 & ) then ! locate boundary-layer convergence point YC near port vortex call blconvergencepoint(Y,GAM,RADII,2,1,UCR,NZUCR,YC,YCFLAG) ! if convergence point was found, eject boundary layer vortex if (YCFLAG.eq.0) then if( debug_ige ) then write(6,*)"NPTS, time, port YC=",NPTS+1,X,YC end if ! reset DXEJECTP DXEJECTP = DXEJECT2 * 8.*Y(2,1)/ & ABS(GAM(1)/(PI*Y(2,1))+1.*UCR(1)) if( debug_ige ) then write(6,*)"DXEJECTP = ",DXEJECTP end if ! add a new boundary-layer vortex nearby call addvortex(YC,YNEW,ZNEW,GAMNEW,DGAMNEW,GAMNEW2 & ,X,Y(1,1),Y(2,1),GAM(1),NPTS,NVAR,P2ND & , 1.,DIS,DXEJECTP,CEJECTA,BZ,VZIN,XEJECTP & ,ZEJECT,SECDK,RADII(NPTS+1),AREA(NPTS+1)) Y(1, NPTS) = YNEW Y(2, NPTS) = ZNEW ZMIN( NPTS) = ZNEW GAM( NPTS) = GAMNEW GAM0( NPTS) = DGAMNEW DGAM( NPTS) =-DGAMNEW*GAM(NPTS) PNEAR(NPTS) = 1 if( debug_ige ) then write(6,*)"NPTS,time,GAM,DGAM,ZNEW=",NPTS,X,GAMNEW & ,DGAMNEW,ZNEW end if ! set YSCAL for RKQC() and CGAM for IGDRVS() YSCAL(1,NPTS) = ABS(YNEW)+TINY YSCAL(2,NPTS) = ABS(ZNEW)+TINY CGAM(NPTS) = GAMNEW ! store GAMNEW in output array GAMP(NPTS,KG) = GAMNEW ! compute contribution to primary vortex circulation decay call creationdrain(dgamcnew,Y(1,1),Y(2,1),GAM(1),R0p & ,YNEW,ZNEW,GAMNEW ,R0p,DXEJECTP,frac) DGAMC( 1,NPTS) = dgamcnew TIMEC(1,1,NPTS) = X TIMEC(2,1,NPTS) = X+DXEJECTP*1. ! write(6,*)" DGAMC(1,",NPTS,")=",dgamcnew ! end found boundary-layer convergence endif ! end port vortex boundary-layer ejection buildup endif ! end port vortex in ground effect endif ! starboard vortex drops into IGE if ( Y(2,2) .le. ZGE ) then ! begin buildup of first ejected starboard vortex if ( XEJECTS .eq. 9999.99 ) then XEJECTS = X DXEJECTS = DXEJECT1*2.5*BZ & /ABS(GAMP(2,1)*1./(PI*BZ)+1.*UCR(1)) !!! DXEJECTS = DXEJECT1*4.*BZ !!! & /ABS(GAMP(2,1)*2./(PI*BZ)+0.*UCR(1)) if( debug_ige ) then write(6,*)"time,DXEJECTS = ",X,DXEJECTS end if endif ! finish buildup of new port-side b.l. vortex: at least DXEJECTS ! since last port vortex was ejected; not reached maximum ! vortices; neither primary vortex has decayed away. if ( (X-XEJECTS.GE.ABS(DXEJECTS)) & .AND. NPTS.LT.NPTMAX & .AND. KFLAG.ne.1 ! & .AND. 0.eq.1 & ) then ! locate boundary-layer convergence point YC near starbrd vortex call blconvergencepoint(Y,GAM,RADII,2,2,UCR,NZUCR,YC,YCFLAG) ! if convergence point was found, eject boundary layer vortex if (YCFLAG.eq.0) then if( debug_ige ) then write(6,*)"NPTS, time, starboard YC=",NPTS+1,X,YC end if ! reset DXEJECTS DXEJECTS = DXEJECT2 * 8.*Y(2,2)/ & ABS(GAM(2)/(PI*Y(2,2))+1.*UCR(1)) if( debug_ige ) then write(6,*)"DXEJECTS = ",DXEJECTS end if ! add a new boundary-layer vortex nearby call addvortex(YC,YNEW,ZNEW,GAMNEW,DGAMNEW,GAMNEW2 & ,X,Y(1,2),Y(2,2),GAM(2),NPTS,NVAR,S2ND & ,-1.,DIS,DXEJECTS,CEJECTA,BZ,VZIN,XEJECTS & ,ZEJECT,SECDK,RADII(NPTS+1),AREA(NPTS+1)) Y(1, NPTS) = YNEW Y(2, NPTS) = ZNEW ZMIN( NPTS) = ZNEW GAM( NPTS) = GAMNEW GAM0( NPTS) = DGAMNEW DGAM( NPTS) =-DGAMNEW*GAM(NPTS) PNEAR(NPTS) = 2 if( debug_ige ) then write(6,*)"NPTS,time,GAM,DGAM,ZNEW=",NPTS,X,GAMNEW & ,DGAMNEW,ZNEW end if ! set YSCAL for RKQC() and CGAM for IGDRVS() YSCAL(1,NPTS) = ABS(YNEW)+TINY YSCAL(2,NPTS) = ABS(ZNEW)+TINY CGAM(NPTS) = GAMNEW ! store GAMNEW in output array GAMP(NPTS,KG) = GAMNEW ! compute contribution to primary vortex circulation decay call creationdrain(dgamcnew,Y(1,2),Y(2,2),GAM(2),R0s & ,YNEW,ZNEW,GAMNEW ,R0s,DXEJECTS,frac) DGAMC( 2,NPTS) = dgamcnew TIMEC(1,2,NPTS) = X TIMEC(2,2,NPTS) = X+DXEJECTS*1. ! write(6,*)" DGAMC(2,",NPTS,")=",dgamcnew ! end found boundary-layer convergence endif ! end starboard vortex boundary-layer ejection buildup endif ! end starboard vortex in ground effect endif ! needed just before existing loop IF ( KMAX.NE.0 & .AND. (KFLAG.eq.1 .OR. & NSTP.eq.MAXSTP .OR. & H.eq.0) ) THEN KOUNT=KOUNT+1 XP(KOUNT)=X DO 15 I=1,NPTS YP(1,I,KOUNT)=Y(1,I) YP(2,I,KOUNT)=Y(2,I) 15 CONTINUE ENDIF IF( (X-X2)*(X2-X1).GE.ZERO .OR. & ZTEST.LE.ZGE*0.15 .OR. & H.le.0. ) THEN if( debug_ige ) then if ( (X-X2)*(X2-X1).GE.ZERO ) then write(6,*)" time > T2, X2=",X2 endif if (ZTEST.le.ZGE*0.15) then write(6,*)"ZTEST < 0.15 ZGE, ZTEST=",ZTEST endif if (H.le.0) then write(6,*)"H < or = 0, H=",H endif end if RETURN ENDIF IF ( Y(2,1).gt.4.0*YSTART(2,1) .OR. & Y(2,2).gt.4.0*YSTART(2,2) ) THEN if( debug_ige ) then write(6,*)" higher than we started!" end if RETURN ENDIF IF(GAM(1) .GE. 0. .AND. GAM(2) .LE. 0.) THEN DO IG=1,NPTS GAM(IG) = 0. DGAM(IG) = 0. GAMP(IG,KG) = 0. KFLAG = 1 ENDDO ENDIF C hnext = sign(max(abs(hnext),hmin),hnext) IF(ABS(HNEXT).LT.HMIN) THEN WRITE(2,991) 991 FORMAT(/'Stepsize smaller than minimum.') CLOSE(2) CCC CLOSE(3) STOP ENDIF H=HNEXT ! time stepping loop !JW write(6,*)"NSTP = ",NSTP 16 CONTINUE if( debug_ige ) then write(6,*)"hi joe, done looping" end if WRITE(2,993) 993 FORMAT(/'Too many steps.') CLOSE(2) CCC CLOSE(3) STOP END !======================================================================= subroutine blconvergencepoint(YZ,G,R,N,J,U,NU,YC,YCFLAG) !----------------------------------------------------------------------- ! ! NAME ! blconvergencepoint: locate surface-velocity convergent point on ! boundary layer for flow field induced by N ! vortices in the frame of primary vortex J ! ! USAGE ! call blconvergencepoint(YZ,G,R,N,J,U,NU,YC,YCFLAG) ! ! INPUT ! YZ(1,N) = y coordinates for N vortices ! YZ(2,N) = z coordinates for N vortices ! G( N) = circulations for N vortices ! R( N) = core radii for N vortices ! J = index of primary vortex associated with convergence pt ! U( NU) = cross wind ! ! OUTPUT ! YC = y coordinate for convergence point ! YCFLAG = 0 (convergence) or 1 (no convergence) ! ! EXAMPLE ! ! NOTES ! vortices are assumed to have Burnham-Hallock velocity profiles: ! ! r^2 Gam0 r ! Gam(r) = Gam0 --------- or V_theta(r) = ---- --------- ! r^2 + R^2 2.pi r^2 + R^2 ! ! this routine returns the location of flow convergence along the ! bottom boundary of the induced flow field for N Burnham-Hallock ! vorticies in the frame of vortex J. ! ! in co-moving frame, for convergence location, U(z) is ignored ! ! AUTHOR ! Joe Werne 2012 December ! !----------------------------------------------------------------------- implicit none integer :: N,NU integer :: J,I,II,ITER,YCFLAG integer, parameter :: MAXITER=1000 real :: YZ(2,N),G(N),R(N),U(NU),YC real :: RHS,RHS1,PII,B,D,Y,YOLD,SURF,V_J & ,COEFFM,COEFFP,ERR ! iteration tolerance ERR = 1e-8 ! 1/pi PII = 1./(4.0*atan(1.0)) ! separation of first two primary vorticies B = sqrt((YZ(1,1)-YZ(1,2))**2+(YZ(2,1)-YZ(2,2))**2) ! initialize horizontal velocity of J vortex V_J = 0.0 ! initialize YCFLAG=1 (no convergence) YCFLAG = 1 ! step through all vortices do I=1,N ! these coefficients are needed to compute the contributions to ! the velocity of the jth vortex by the ith vortex and its image COEFFM = 0.5*PII*G(I) & / (R(I)**2+(YZ(1,j)-YZ(1,i))**2+(YZ(2,j)-YZ(2,i))**2) COEFFP = 0.5*PII*G(I) & / (R(I)**2+(YZ(1,j)-YZ(1,i))**2+(YZ(2,j)+YZ(2,i))**2) ! lateral avection of j-vortex (neglecting cross wind) V_J = V_J + COEFFM*( YZ(2,I) - YZ(2,J) ) & + COEFFP*( YZ(2,I) + YZ(2,J) ) enddo ! guess location Y of the convergence point along the BL D = sqrt(3.)*(abs(YZ(2,J))*2./B)**1.57 if ( G(J) .gt. 0. ) then Y = YZ(1,J) + D else Y = YZ(1,J) - D endif ! iterate to find convergence point do ITER=1,MAXITER ! store Y in YOLD YOLD = Y c write(6,*)"BLCONPT - YOLD,Y = ",YOLD,Y ! initialize surface velocity in j-vortex's co-moving frame SURF = -V_J ! I cycles through all indices except J, starting with I=J+1 do II=1,N-1 I = mod(II-1+J,N) + 1 ! surface flow in j-vortex's co-moving frame SURF = SURF & + G(I)*PII*YZ(2,I)/(R(I)**2+(Y-YZ(1,I))**2+(YZ(2,I)**2)) enddo ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! note: if we added the contribution from the jth vortex, we'd have the ! following: ! ! SURF = SURF + G(J)/PI*Z(J) / ( R*R + (Y-Y(J))**2 + Z(J)**2 ) ! ! but we want to find the point Y where SURF=0; therefore, ! ! 0 = SURF + G(J)/PI*Z(J) / ( R*R + (Y-Y(J))**2 + Z(J)**2 ) ! G(J)/PI*Z(J) / ( R*R + (Y-Y(J))**2 + Z(J)**2 ) = -SURF ! R*R + (Y-Y(J))**2 + Z(J)**2 = -G(J)/PI*Z(J)/SURF ! ! note that the LHS is positive definite, so the RHS must be positive. ! if it's not, you have a problem. ! ! (Y-Y(J))**2 = - G(J)/PI*Z(J)/SURF - R*R - Z(J)*Z(J) ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! RHS1 = - G(J)*PII*YZ(2,J)/SURF - R(I)*R(I) - YZ(2,J)*YZ(2,J) if (RHS1 .gt. 0.) then D = sqrt(RHS1) if (G(J) .gt. 0.) then Y = YZ(1,J) + D else Y = YZ(1,J) - D endif else c write(6,*)"RHS1 = ",RHS1 c write(6,*)"RHS1 was negative - no convergence point." Y=-999999.99 YCFLAG=1 goto 86 endif ! if ( abs(Y-YOLD) .lt. abs(YOLD)*ERR ) then ! iteration converged YCFLAG=0 goto 86 endif ! ! end ITER loop enddo ! ! break out of ITER loop 86 continue ! ! save convergent point YC = Y ! return end !======================================================================= subroutine gamdotfriction(DGAMF,GAM,Z,W,DIS,CKE,R0) !----------------------------------------------------------------------- real DGAMF,GAM,Z,W,DIS,CKE,R0,USURF COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC ! estimate surface velocity beneath primary vortices USURF = GAM/(Z*PI) ! estimate drag coefficient for primary vortices. ! also check skin friction for roughness length. ! specify viscosity elsewhere and pass in. DIS=nu ! define variables. ! ! check appropriate value for vortex core radius R ! DIS = kinematic viscosity, normally called nu [m^2/s] ! RE = max streamwise surface reynolds # under vortex ! CD = skin-friction Drag coefficient associated with RE ! R0 = vortex core radius [m] ! RE = ABS(GAM)/(PI*DIS) CD = 0.032 * RE**(-1./7.) ! fully rough ... CD = 0.032 / 0.027 * 0.017 ! add additional friction drag to primary vortices ! ! dGam = - Cke.CF.Gam.|Gam|/(Z.PI)^2 - Gam.(Z_dot/Z) ! ---- ------------------------------------------- ! dt 1 + 2.ln(2.Z/R0) ! ! we will initially neglect the second term. ! note: CF should blow up when ejections occur. ! we can handle this by tripling CKE during ejection ! events. note, CKE is set to 1.0 at the top of this ! routine. we can improve on this handling of ejection- ! enhanced circulation decay by including pressure-drag ! physics, but we'll save that for a later version. ! DGAMF = - GAM * ( CKE*CD*ABS(GAM)/(Z*PI)**2 + max(0.,W/Z) ) & /(1.+2.*alog(2.*Z/R0)) ! DGAMF = - GAM * ( CKE*CD*ABS(GAM)/(Z*PI)**2 + W/Z ) ! & /(1.+2.*alog(2.0*2.*Z/R0)) ! ! DGAMF = - ( CKE*CD*GAM*ABS(GAM/(Z*PI)**2) + GAM*W/Z ) ! & /(1.+2.*alog(2.*Z/R0)) ! RETURN END !======================================================================= subroutine addvortex(YC,YNEW,ZNEW,GAMNEW,DGAMNEW,GAMNEW2 & ,time,Y,Z,GAM,NPTS,NVAR,PS2ND & ,S,DIS,DXEJECT,CEJECTA,B0,V0,XEJECT,ZEJECT & ,SECDK,RADIUS,AREA) !----------------------------------------------------------------------- PARAMETER (NZMAX=401) COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) & ,QZ(NZMAX),Q(NZMAX) & ,UCRZ(NZMAX),UCR(NZMAX) & ,UHWZ(NZMAX),UHW(NZMAX) & ,PRESS(NZMAX),THETA(NZMAX),BVFSQ(NZMAX),RHO(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /ATMB2/DUCR(NZMAX),DDUCR(NZMAX) COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, & THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, & DELTH1,DELTH2,FACMIN,FACMAX,DELFAC INTEGER NPTS,NVAR,PS2ND REAL YC,YNEW,ZNEW,GAMNEW,DGAMNEW,GAMNEW2 & ,time,Y,Z,GAM,S & ,DIS,DXEJECT,CEJECTA,B0,V0,XEJECT,ZEJECT & ,SECDK & ,QK & ,RADIUS,AREA ! add a new boundary-layer vortex nearby NPTS = NPTS + 1 NVAR = NPTS*2 ! store the 1st port or 1st starboard 2ndary vortex number in P2ND IF (PS2ND.EQ.0) THEN PS2ND=NPTS ENDIF ! vortex location rule: left of left vortex, midway between ! the left vortex and the convergence point for flow along the ! boundary, evaluated in the laterally moving frame of the ! left vortex and 5% up from the ground. 5% is arbitrary. base ! on the height of the boundary layer in a later version. ! ! ejection occurs midway between convergence point and primary YNEW = 0.5*( Y + YC ) ! if (YC.lt.Y) then ! YNEW = Y-0.5*sqrt(3.)*Z ! else ! YNEW = Y+0.5*sqrt(3.)*Z ! endif ! ZNEW = .05* Z ! ZNEW = .10* Z ! ZNEW = Z * ZEJECT !JOE - LOOK HERE ! ZNEW = 0.2* Z ! ZNEW = 0.5* Z ! rule for ejected vortex circulation: ! ! GAMMA_dot = -0.5(Usurf)|Usurf| 0.166 Re^(-1/7) Ceject ! ! b.l. vorticity is Usurf/delta ! ! b.l. height is delta = z 0.166 Re^(-1/7) ! ! rate of circulation production deduced from Doligalski ! and Walker (84). ! ! old thoughts below ! ! rule for ejected vortex circulation: the surface vorticity ! is dU/dz = -u*|u*|/nu, where u* is the friction velocity, ! u*=Usurf.sqrt(Cf/2)=>dU/dz=-Usurf |Usurf| Cf/(2.nu). note ! that Usurf is the surface velocity in the lab frame, and ! the sign of the ejected vortex is opposite the surface-flow ! direciton in the lab frame. ! ! the circulation is the area integral of the vorticity, so ! we multiply the surface-vorticity dU/dz by the area defined ! by: the boundary-layer depth, times the length the ejected ! vortex travels in the time between vortex ejections, times ! a constant geometrical factor CEJECTA. ! ! estimate surface velocity beneath primary port vortex in ! order to estimate circulation strength of ejected vortex. ! cross wind enhances USURF1. ! ! note: at Re=10^7, circulation production (Gamma_dot) is ! Gamma_dot = -Usurf|Usurf| 0.166 Re^(-1/7) Ceject/2 ! = -Usurf|Usurf| Ceject/120 ! ! and we set Ceject=50. ! -Usurf ! another estimate for Gamma_dot = ----- |Usurf| d Ceject ! d ! Gamma_dot = -Usurf |Usurf| Ceject ! ! which is going to result in Ceject ~ O(1) ! ! Usurf is maximum surface velocity beneath primary vortex USURF = GAM*Z/((Z*Z+0.*(YNEW-Y)*0.*(YNEW-Y))*PI) + UCR(1) RE = ABS( GAM/(PI*DIS) + UCR(1)*Z/DIS ) BLDELZ = Z*0.166*RE**(-1./7.) ! fully rough ... BLDELZ = Z*0.1 AREA = BLDELZ*ABS(USURF)*(time-XEJECT) * CEJECTA RADIUS = SQRT(AREA/PI) ZNEW = 3.*BLDELZ ZNEW = 0.05 * Z ZNEW = 0.30 * Z ZNEW = ABS(Y-YNEW) * 0.5 ZNEW = 3.0 * BLDELZ ZNEW = 0.10 * Z ZNEW = 0.20 * Z GAMNEW = -USURF*ABS(USURF)*(time-XEJECT)*CEJECTA ! & * ( abs(UCR(1)/1.43) )**0.0 ! ! GAMNEW = -USURF*ABS(USURF)*0.5 * 0.166*RE**(-1./7.) *CEJECTA ! & *(time-XEJECT) ! ! & *DXEJECT ! ! CF = 0.0277 * RE**(-1./7.) ! DUDZ = -USURF*ABS(USURF)*CF/(2.*DIS) ! AREA = Z*0.1643*RE**(-1./7.)*ABS(USURF*DXEJECT) ! GAMNEW = S*max(AREA*DUDZ*CEJECTA*S,-GAM*0.0*S) ! if (ABS(GAMNEW).lt.GAMSMALL) then ! GAMNEW = GAMSMALL*SIGN(1.,GAMNEW) ! endif USURF2 = GAM*Z/((Z*Z+0.*(YNEW-Y)*0.*(YNEW-Y))*PI) RE = ABS( GAM/(PI*DIS) ) GAMNEW2= -USURF2*ABS(USURF2)*0.5*0.166*RE**(-1./7.)*CEJECTA & *(time-XEJECT) ! & *DXEJECT ! ! CF = 0.0277 * RE**(-1./7.) ! DUDZ = -USURF*ABS(USURF)*CF/(2.*DIS) ! AREA = Z*0.1643*RE**(-1./7.)*ABS(USURF*DXEJECT) ! GAMNEW2= S*max(AREA*DUDZ*CEJECTA*S,-GAM*0.0*S) ! lin-model OGE decay rate: -(0.1020 + 0.687 EDR*)*GAM0*V0/B0 ! lin-model NGE decay rate: -(0.0715 + 0.108 EDR*)*GAM0*V0/B0 ! use this until you have better theory iztopq = nzq QK = FK( ZNEW , QZ, Q, NZQ, iztopq ) QSTAR = (QK*B0)**0.333333/V0 ! DGAMNEW = - sign( 0.0715 + 0.108*QSTAR , GAMNEW )*V0*V0*TWOPI ! DGAMNEW = - sign( 0.0472 + 0.400*QSTAR , GAMNEW )*V0*V0*TWOPI ! DGAMNEW = - sign( 0.0400 + 0.020*QSTAR , GAMNEW )*V0*V0*TWOPI * 1. ! DGAMNEW = - sign( 0.1000 + 0.000*QSTAR , GAMNEW )*V0*V0*TWOPI * 1. ! DGAMNEW = - sign( 0.0100 + 0.000*QSTAR , GAMNEW )*V0*V0*TWOPI * 1. DGAMNEW = ( 0.0045 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.0090 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.0008 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.0135 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.0070 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.0035 + 0.000*QSTAR ) * ABS(USURF)/BLDELZ * 1. DGAMNEW = ( 0.02117+ 0.6964*QSTAR) * ABS(USURF)/BLDELZ*SECDK & / ( 0.15 + 0.6*abs(UCR(1)) ) ! & / ( abs(UCR(1)/1.43) )**0.7 if (USURF*GAMNEW.gt.0.) then DGAMNEW = DGAMNEW * 1.04 else DGAMNEW = DGAMNEW * 0.96 endif ! store the time of the last secondary vortex ejected XEJECT = time return end !======================================================================= subroutine creationdrain(dgam,y,z,gam,r0,y2,z2,gam2,r02,delt,frac) !----------------------------------------------------------------------- ! ! NAME ! creationdrain: compute the circulation decay of a primary vortex ! resulting from the creation of a boundary-layer ! vortex. ! ! USAGE ! call creationdrain(dgam,y,z,gam,r0,y2,z2,gam2,r02,delt,frac) ! ! INPUT ! y = primary vortex lateral position ! z = primary vortex vertical position ! gam = primary vortex circulation ! r0 = primary vortex core radius ! y2 = secondary vortex lateral position ! z2 = secondary vortex vertical position ! gam2 = secondary vortex circulation ! r02 = secondary vortex core radius ! delt = duration of secondary vortex creation ! frac = fraction of creation energy offered by primary vortex ! ! OUTPUT ! dgam - primary vortex circulation decay required to conserve KE ! ! EXAMPLE ! ! NOTES ! vortices are assumed to have Burnham-Hallock velocity profiles: ! ! r^2 Gam0 r ! Gam(r) = Gam0 --------- or V_theta(r) = ---- --------- ! r^2 + R^2 2.pi r^2 + R^2 ! ! KE of a single primary vortex and it's image ... ! ! KE1 = pi.k1.k1 + 2.pi.k1.k1.ln(2.z1/R1) = pi.k1.k1(1+2.ln(2.z1/R1)) ! ! KE2 = pi.k2.k2 + 2.pi.k2.k2.ln(2.z2/R2) = pi.k2.k2(1+2.ln(2.z2/R2)) ! ! KE1' = KE1 - f.KE2 ! ! pi.k1'.k1'(1+2.ln(2.z1/R1) = pi.k1.k1(1+2.ln(2.z1/R1)) ! - f.pi.k2.k2(1+2.ln(2.z2/R2)) ! ! k1'.k1' = k1.k1 - f.k2.k2(1+2.ln(2.z2/R2))/(1.+2.ln(2.z1/R1)) ! ! AUTHOR ! Joe Werne 2013 August ! !----------------------------------------------------------------------- implicit none real :: dgam,y,z,gam,r0,y2,z2,gam2,r02,delt,frac real :: gamnew,fraction,coeffb,coeffc,a,d,ke1,ke2,ke3 !! a = 0.5*sqrt( (z-z2)*(z-z2) + (y-y2)*(y-y2) ) !! d = 0.5*sqrt( (z+z2)*(z+z2) + (y-y2)*(y-y2) ) !! coeffb = 4.*gam2*alog(a/d)/(1.+2.*alog(2.*z/r0)) !! coeffc = gam2*gam2*(1.+2.*alog(2.*z2/r02))/(1.+2.*alog(2.*z/r0)) !! & - gam*gam !!! ke1 = gam *gam *(1.+(r0 *r0 +2.*z *z )/(z *sqrt(r0 *r0 +z *z )) !!! & *alog(z /r0 +sqrt(1.+z *z /(r0 *r0 )))) !!! !!! ke2 = gam2*gam2*(1.+(r02*r02+2.*z2*z2)/(z2*sqrt(r02*r02+z2*z2)) !!! & *alog(z2/r02+sqrt(1.+z2*z2/(r02*r02)))) !!! !!! ke3 = ke1 - frac*ke2 !!! write(6,*) "ke1,2,3=",ke1,ke2,ke3 !!! gamnew = sqrt( ke3/(1.+(r0 *r0 +2.*z *z )/(z *sqrt(r0 *r0 +z *z )) !!! & *alog(z /r0 +sqrt(1.+z *z /(r0 *r0 )))) ) !!! ke1 = gam*gam ke2 = gam2*gam2 ke3 = ke1 - frac*ke2 gamnew = sign( sqrt(amax1(ke3,0.0)) , gam ) ! write(6,*) "ke1,2,3,gamnew=",ke1,ke2,ke3,gamnew ! fraction of secondary-vortex creation taken from primary vortex ! (it could come from the mean wind too.) dgam = ( gamnew - gam ) / delt return end