! JAW_V01 ! C ALGORITHM FOR COMPUTING TRAJECTORIES OF TRAILING VORTICES C (BASED ON F. PROCTOR'S SHEAR MODEL) C ! JOE - new flag variables are being passed ! JOE - added headwind variables NHWA, HWZA, HWA ! JOE - hard code default values for the following control flags: ! FLAG_PARTICLE = 0.17 ! FLAG_GRAV = 1.0 ! FLAG_VCDOTSC = 1.0 ! FLAG_VBDOTSB = 1.0 ! FLAG_IGE_NGE = 1.0 * 0.0 ! FLAG_DYNBUOY = 0.0 ! FLAG_LINK = 0.0 ! ! JOE - if z < ZMFA * b0, go into viper_ige(), this is z-limit for NGE ! ! ! JOE - and reuse the control flags as follows: ! FLAG_ENTR_CIR_EDR = AFLAG_PARTICLE ! FLAG_ENTR_MOM_EDR = AFLAG_GRAV ! FLAG_ENTR_DEN_EDR = AFLAG_VCDOTSC ! FLAG_ENTR_CIR_N2 = AFLAG_VBDOTSB ! FLAG_ENTR_MOM_N2 = AFLAG_IGE_NGE ! FLAG_ENTR_DEN_N2 = AFLAG_DYNBUOY ! FLAG_CIR_BARO = AFLAG_LINK ! ! QZ,Q = (z,edr) ! SUBROUTINE VIPER( & YZA,ZZA,VZA,BZA,ZMFA,ZGFA,GMFA,GRFA,GNGA, X NTA,TZA,TA,rhoa,bvfsqa, & NUA,UZA,UA,NHWA,HWZA,HWA,NQA,QZA,QA, & ATLAG_FAC,ATIGE1_FAC, X AT1,AT2,NTOGA,AT1IM,AT2IM,NTIMGA, X AT1GE1,AT2GE1,NTGE1A,AT1GE2,AT2GE2,NTGE2A, X KMAXA,ADTFAC,AEPS,AH1,AHMIN, X AIRPORTALTA,RHOAIRA,ACSPEEDA, & R0p,R0s,CEJECTA,CKE, & TEJECT1,TEJECT2,ZEJECT, & SECDK, SECDIFF, & AFLAG_ENTR_CIR,AFLAG_ENTR_CIR_EDR,AFLAG_ENTR_CIR_N2, & AFLAG_ENTR_MOM,AFLAG_ENTR_MOM_EDR,AFLAG_ENTR_MOM_N2, & AFLAG_ENTR_DEN,AFLAG_ENTR_DEN_EDR,AFLAG_ENTR_DEN_N2, & AFLAG_GRAV, & AFLAG_VCDOTSC, & AFLAG_VBDOTSB, & AFLAG_CI, & AFLAG_AEFF_B, & AFLAG_AEFF_S, & AFLAG_LINK, & AFLAG_IGE_NGE, & AFLAG_DYNBUOY, & AFLAG_PARTICLE, & NB,TB, & YPB,ZPB,GPB, & YSB,ZSB,GSB,FLAG, & BASENAME, & ALINEAR,hw_avg, & ATLD,ATSSD) C c JOE - changed NVAR to 10 so i can include particle density variable PARAMETER (NSTPMX=5001,NVAR=10,NZMAX=401,NNUT=8, X MAXSAV=20000,NPTMAX=160,NYZMAX=2*NPTMAX & ,NXIMAX=721) C c JOE - add headwind variables, HWZA(*),HWA(*), c AIRPORTALT,RHOAIR,ACSPEED REAL YZA,ZZA,VZA,BZA,ZMFA,ZGFA,GMFA,GRFA,GNGA REAL NTA INTEGER NUA,NHWA,NQA,ALINEAR,LINEAR REAL TZA(*),TA(*), rhoa(*), bvfsqa(*), & UZA(*),UA(*),HWZA(*),HWA(*),QZA(*),QA(*) REAL TB(*),YPB(*),ZPB(*),GPB(*),YSB(*),ZSB(*),GSB(*) REAL ATLAG_FAC,ATIGE1_FAC, & AT1,AT2,NTOGA,AT1IM,AT2IM,NTIMGA, & AT1GE1,AT2GE1,NTGE1A,AT1GE2,AT2GE2,NTGE2A, & KMAXA,ADTFAC,AEPS,AH1,AHMIN, & AIRPORTALTA,RHOAIRA,ACSPEEDA, & R0p,R0s,CEJECTA,CKE, & TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF REAL DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 REAL TLREPORT, TSSREPORT, TL, TSS, TLD, TSSD C C JOE - a bunch of flag variables to control inclusion c of various terms in the evolution equations. REAL AFLAG_ENTR_CIR , AFLAG_ENTR_CIR_EDR , AFLAG_ENTR_CIR_N2 REAL AFLAG_ENTR_MOM , AFLAG_ENTR_MOM_EDR , AFLAG_ENTR_MOM_N2 REAL AFLAG_ENTR_DEN , AFLAG_ENTR_DEN_EDR , AFLAG_ENTR_DEN_N2 REAL AFLAG_GRAV REAL AFLAG_VCDOTSC REAL AFLAG_VBDOTSB REAL AFLAG_CI REAL AFLAG_AEFF_B REAL AFLAG_AEFF_S REAL AFLAG_LINK REAL AFLAG_IGE_NGE REAL AFLAG_DYNBUOY REAL AFLAG_PARTICLE REAL ATLD , ATSSD REAL FLAG INTEGER i_old CHARACTER(80) BASENAME C COMMON /IFLAGB/IFLAG COMMON /KFLAGB/KUCR,KFLAG C COMMON /OUTB/TOUT(5001), X YPOUT(5001),ZPOUT(5001),GPOUT(5001), X YSOUT(5001),ZSOUT(5001),GSOUT(5001),NTOUT C COMMON /VIPERB/GAMZ,BSV,FAC11,BETA1,ALPH1,FAC12,FAC13, X BETA2,BETA3,ALPH2,FAC21,FAC22,FAC23,THLF,FP3,CSH & , TLREPORT, TSSREPORT, TLD, TSSD, hw_fac, i_old COMMON /SARPB/C0,BZ,VZIN,R3,F0,VSB,VVSB,RVSB2,AFAC,DBFAC, X DIMT,BRATIO,GAMAVB,ARGTB,RHOINIT,TDEGINIT COMMON /SARPB2/BZ2,VVSB2 COMMON /TSTARB/TS,BS,DLOGBS,TCFAC,IFIRST C COMMON /PATH/TT(NSTPMX),VV(NVAR,NSTPMX),SEP(NSTPMX) COMMON /PATHPS/SEPPS(NSTPMX),SEPY(NSTPMX),SEPZ(NSTPMX) 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 COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) X ,QZ(NZMAX),Q(NZMAX) X ,UCRZ(NZMAX),UCR(NZMAX) X ,UHWZ(NZMAX),UHW(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /ATMB2/DUCR(NZMAX),DDUCR(NZMAX) COMMON /RHOB/BVFSQ(NZMAX),RHO(NZMAX) COMMON /STRAT/STRN COMMON /CDB/REFAC,CDFAC,REHI,RELO,CDHI,CDLO,XNUT(NNUT),XNU(NNUT) COMMON /IZTOPB/ZLAST,QB,QB1,QB4 & ,IZTOPT,IZTOPQ,IZTOPU,IZTOPHW,NUTTOP COMMON /IZTPRB/IZTOPR COMMON /QB58B/QB5,QB8,BLAST,IZ2PDU COMMON /GTOGB/YZPI(2,NPTMAX,NXIMAX),GAMI(NPTMAX,NXIMAX) & ,XI(NXIMAX),NXI C COMMON /PATH1/DTSAV,TP(MAXSAV),KMAX,KOUNT COMMON /PATH2/YZP(2,NPTMAX,MAXSAV) COMMON /NGAM/GAM(NPTMAX),CGAM(NPTMAX),NPTS COMMON /DGAMB/DGAM(NPTMAX),DELGAM COMMON /DBDTB/DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC C JOE - common block of flag variables to control inclusion c of various terms in the evolution equations. REAL FLAG_ENTR_CIR , FLAG_ENTR_CIR_EDR , FLAG_ENTR_CIR_N2 REAL FLAG_ENTR_MOM , FLAG_ENTR_MOM_EDR , FLAG_ENTR_MOM_N2 REAL FLAG_ENTR_DEN , FLAG_ENTR_DEN_EDR , FLAG_ENTR_DEN_N2 REAL FLAG_GRAV,FLAG_VCDOTSC,FLAG_VBDOTSB,FLAG_LINK REAL AEFF_S,FLAG_IGE_NGE,FLAG_DYNBUOY,FLAG_PARTICLE REAL AEFF_B,FLAG_CI COMMON /FLAGS/FLAG_ENTR_CIR,FLAG_ENTR_CIR_EDR,FLAG_ENTR_CIR_N2 & ,FLAG_ENTR_MOM,FLAG_ENTR_MOM_EDR,FLAG_ENTR_MOM_N2 & ,FLAG_ENTR_DEN,FLAG_ENTR_DEN_EDR,FLAG_ENTR_DEN_N2 X ,FLAG_GRAV,FLAG_VCDOTSC,FLAG_VBDOTSB,FLAG_LINK X ,AEFF_S,FLAG_IGE_NGE,FLAG_DYNBUOY,FLAG_PARTICLE X ,AEFF_B,FLAG_CI,LINEAR COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3 REAL VSTART(NVAR),XNUTVL(NNUT),XNUVAL(NNUT) REAL GAMHA(NSTPMX),GAMHB(NSTPMX), X GAMDA(NSTPMX),GAMDB(NSTPMX), X VVYP(NSTPMX),VVYS(NSTPMX),VVZP(NSTPMX),VVZS(NSTPMX) REAL GAOUT(NSTPMX),GBOUT(NSTPMX) REAL TDEGZ2(NZMAX),BVFSQ2(NZMAX),QNSQ2(NZMAX),QNSQ(NZMAX) c JOE - add variable to report Fr in DENSITY.DAT file REAL FROUT(NSTPMX),QSTAR(NSTPMX),COSANGALL(NSTPMX) C real mu, nu, mu_ref logical debug1, debug2, debug4, debug6, debug_ige common /debug/ debug1, debug2, debug4, debug6, debug_ige c *** The logical variables debug2 and debug4 control output sent c *** units 1, 2, 4, and 6. debug1 = .false. debug2 = .false. debug4 = .false. debug6 = .false. debug_ige = .false. ! if i_old=1, use code tom was using, else use updated code. i_old = 1 i_old = 0 ! initialize TLREPORT, TSSREPORT TLREPORT = 0.0 TSSREPORT = 0.0 if( i_old .eq. 1 ) then AFLAG_ENTR_CIR = 0.1285 AFLAG_ENTR_MOM = 0.0 AFLAG_ENTR_DEN = 0.2 AFLAG_ENTR_CIR_EDR = 0.5 AFLAG_ENTR_MOM_EDR = -0.1 AFLAG_ENTR_DEN_EDR = 0.0 AFLAG_ENTR_CIR_N2 = 0.8 AFLAG_ENTR_MOM_N2 = 1.0 AFLAG_ENTR_DEN_N2 = 2.0 AFLAG_AEFF_S = 0.1 endif DATA XNUTVL/-50.,0.,10.,15.,20.,30.,40.,60./, X XNUVAL/.092,.132,.141,.145,.15,.16,.169,.188/ C AIRPORTALT = AIRPORTALTA RHOAIR = RHOAIRA if (ACSPEEDA.ne.0.) then ACSPEEDI = 1./ACSPEEDA else ACSPEEDI = 0.0 stop endif !!!! ACSPEEDI = 0.0 IFLAG = 0 FLAG = 0 C YZ = YZA ZZ = ZZA VZ = VZA BZ = BZA ZIMFAC = ZMFA ZGEFAC = ZGFA GAMFAC = GMFA GERFAC = GRFA GEANG = GNGA NZTDEG = INT(NTA) NTAABS = IABS(NZTDEG) NZUCR = NUA NZUHW = NHWA NZQ = NQA WPOLD = -VZ WSOLD = -VZ WP = -VZ WS = -VZ tminp = 9999. tmins = 9999. LINEAR = ALINEAR tlag = ATLAG_FAC * BZ / VZ tige1 = ATIGE1_FAC * BZ / VZ C DO I=1,NTAABS TDEGZ(I) = TZA( I) TDEG( I) = TA( I) rho( i) = rhoa( i) bvfsq(i) = bvfsqa(i) ENDDO DO I=1,NUA UCRZ(I) = UZA(I) UCR(I) = UA(I) ENDDO ! joe - turning off headwind effect ! write(6,*)" warning: headwind effect turned off in viper.f" DO I=1,NHWA UHWZ(I) = HWZA(I) UHW(I) = HWA(I) ENDDO DO I=1,NQA QZ(I) = QZA(I) Q(I) = QA(I) ENDDO C NTOG = INT(NTOGA) NTIMG = INT(NTIMGA) NTGE1 = INT(NTGE1A) NTGE2 = INT(NTGE2A) KMAX = INT(KMAXA) C T1 = AT1 T2 = AT2 T1IM = AT1IM T2IM = AT2IM T1GE1 = AT1GE1 T2GE1 = AT2GE1 T1GE2 = AT1GE2 T2GE2 = AT2GE2 DTFAC = ADTFAC EPS = AEPS H1 = AH1 HMIN = AHMIN C JOE - set variables in FLAGS common block FLAG_ENTR_CIR = AFLAG_ENTR_CIR FLAG_ENTR_MOM = AFLAG_ENTR_MOM FLAG_ENTR_DEN = AFLAG_ENTR_DEN FLAG_ENTR_CIR_EDR = AFLAG_ENTR_CIR_EDR FLAG_ENTR_MOM_EDR = AFLAG_ENTR_MOM_EDR FLAG_ENTR_DEN_EDR = AFLAG_ENTR_DEN_EDR FLAG_ENTR_CIR_N2 = AFLAG_ENTR_CIR_N2 FLAG_ENTR_MOM_N2 = AFLAG_ENTR_MOM_N2 FLAG_ENTR_DEN_N2 = AFLAG_ENTR_DEN_N2 FLAG_GRAV = AFLAG_GRAV FLAG_VCDOTSC = AFLAG_VCDOTSC FLAG_VBDOTSB = AFLAG_VBDOTSB FLAG_CI = AFLAG_CI AEFF_B = AFLAG_AEFF_B AEFF_S = AFLAG_AEFF_S FLAG_LINK = AFLAG_LINK FLAG_IGE_NGE = AFLAG_IGE_NGE FLAG_DYNBUOY = AFLAG_DYNBUOY FLAG_PARTICLE = AFLAG_PARTICLE c JOE - and hard code default values for the variables you used FLAG_PARTICLE = 0.17 C C****************************** NEW ****************************** CCC READ(1,*) YZ,ZZ,VZ,BZ CCC YZ = ACYZ CCC ZZ = ACZZ CCC VZ = ACVZ CCC BZ = ACBZ C******************************************************************* CCC READ(1,*) ZIMFAC,ZGEFAC,GAMFAC,GERFAC,GEANG CCC READ(1,*) T1,T2,NT CCC READ(1,*) T1IM,T2IM,NTIM CCC READ(1,*) T1GE1,T2GE1,NTGE1 CCC READ(1,*) T1GE2,T2GE2,NTGE2 CCC READ(1,*) KMAX,DTFAC,EPS,H1,HMIN CCC CLOSE(1) C c JOE - shift all times by T0 so the vortices start where you c measure them to be c c override initial conditions in ACDATA c iover = 0 -- use ACDATA initial condition c iover = 1 -- override iniitial condition c iover is set below T0 = 0.0 C if( debug2 ) then OPEN(2,FILE='VIPER18.OUT',STATUS='UNKNOWN') WRITE(2,1) yz,zz,vz,bz, x zimfac,zgefac,gamfac,gerfac,geang, x t1,t2,ntog, x t1im,t2im,ntimg, x t1ge1,t2ge1,ntge1, x t1ge2,t2ge2,ntge2, x kmax,dtfac,eps,h1,hmin 1 FORMAT(/' yz,zz,vz,bz'/' zimfac,zgefac,gamfac,gerfac,geang'/ x ' t1,t2,ntog'/' t1im,t2im,ntimg'/ x ' t1ge1,t2ge1,ntge1'/' t1ge2,t2ge2,ntge2'/ x ' kmax,dtfac,eps,h1,hmin'/ x 1x,4(1pe15.5)/1x,5(1pe15.5)/ x 1x,2(1pe15.5),i15/1x,2(1pe15.5),i15/ x 1x,2(1pe15.5),i15/1x,2(1pe15.5),i15/ x 1x,i15,4(1pe15.5)) end if C NXIIM = IABS(NTIMG) NXIGE1 = NTGE1 NXIGE2 = NTGE2 C KFLAG = 0 KUCR = 0 IF(NTOG .LT. 0) THEN NTOG = -NTOG KUCR = 1 ENDIF C C CSH = 0.0 for Version 1.0 (v00 - use NTIMG > 0) CSH = 0.0 IF(NTIMG .LT. 0) THEN NTIMG = IABS(NTIMG) CSH = 0.25 ENDIF C BRATIO = 1. C VZIN = VZ R3 = 1./3. F0 = 0.7676**0.75 VSB = VZ/BZ BSV = BZ/VZ VVSB = VZ*VSB RVSB2 = .25/(VSB*VSB) C c JOE - what's up? ZIM = ZIMFAC*BZ IF(ZIM .GT. .99*ZZ) ZIM = 0.99*ZZ C if( debug2 ) then WRITE(2,3) VZ,ZIM 3 FORMAT(/1X,'VZ,ZIM = ',2(1PE15.5)) end if C C FILL /VIPERB/ PI = 4.*ATAN(1.) TWOPI = 2.*PI GAMLC = 3./8. BETA1 = 3./5. ALPH1 = 4./3. CC1 = 0.19 A1 = 0.42 BETA2 = 0.68 BETA3 = 1./4. ALPH2 = 15./8. CC2 = 0.22 A2 = 0.035 THLF = 0. FP3 = 1.27*ALOG(0.3) + 0.57 FAC11 = VVSB*GAMLC*BETA1 FAC12 = VSB*CC1 FAC13 = VSB*VSB*A1 FAC21 = PI*VZ*VZ FAC22 = VSB*CC2 FAC23 = TWOPI*VVSB*A2 CSH = 0.25 TAN3 = tan(3./180.*PI) hw_fac = 1.0/( 1.0 - hw_avg/ACSPEEDA ) C C FILL /THB/ 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 C if( debug2 ) then WRITE(2,1051) 1051 FORMAT(/1X,'GAMMA ROTATION DECAY (ANGLE,FACTOR)') DO 1055 I=1,25 ARG = TWOPI*FLOAT(I-1)/24. THFOUT = THFACF(ARG) ARG = ARG*180/PI WRITE(2,1053) ARG,THFOUT 1053 FORMAT(1X,2F20.5) 1055 CONTINUE end if C TZ = 273.15 ADLAPS = 0.00976 C C READ PROFILES CCC OPEN(1,FILE='TDATA') CCC READ(1,*) NZTDEG IF(NZTDEG.LT.0) THEN TZ = 0. ADLAPS = 0. NZTDEG = -NZTDEG ENDIF CCC DO 10 IZ=1,NZTDEG CCC 10 READ(1,*) TDEGZ(IZ),TDEG(IZ) CCC CLOSE(1) c Tom - Compute the kinematic viscosity at the surface (needed by c viper_ige) mu_ref = 1.7885e-5 T_ref = 288.15 T_s = 110.0 mu = mu_ref*(Tdeg(1)/T_ref)**1.5*(T_ref+T_s)/(Tdeg(1)+T_s) nu = mu/rho(1) DIS = nu c print *, 'T, nu = ', Tdeg(1), DIS C CCC OPEN(1,FILE='UDATA') CCC READ(1,*) NZUCR CCC DO 20 IZ=1,NZUCR CCC 20 READ(1,*) UCRZ(IZ),UCR(IZ) CCC CLOSE(1) C 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 C if( debug2 ) then WRITE(2,27) (iz,ucrz(iz),ucr(iz),ducr(iz),dducr(iz), X iz=1,nzucr) 27 FORMAT(/' udat,dudat,ddudat'/(1x,i10,4(1pe15.5))) end if C CCC OPEN(1,FILE='QDATA') CCC READ(1,*) NZQ CCC DO 40 IZ=1,NZQ CCC 40 READ(1,*) QZ(IZ),Q(IZ) CCC CLOSE(1) if( debug2 ) then WRITE(2,47) (iz,qz(iz),q(iz),iz=1,nzq) 47 FORMAT(/' qdata'/(1x,i10,2(1pe15.5))) WRITE(2,49) 49 FORMAT() end if C c JOE - add IZTOPHW IZTOPT = NZTDEG IZTOPQ = NZQ IZTOPU = NZUCR IZ2PDU = NZUCR IZTOPHW = NZUHW NUTTOP = NNUT ZLAST = ZZ BLAST = BZ QB1 = 0. QB4 = 0. QB5 = 0. QB8 = 0. C IZTOPR = NZTDEG C c JOE - override initial conditions in ACDATA c iover = 0 -- use ACDATA initial condition c iover = 1 -- override iniitial condition c VOS and VOP are speeds, hence positive (it appears) iover = 0 ioverx= 1 ioverx= 0 VOS = 1.7656 VOP = 2.38829 ZOS = 241.8 ZOP = 233.5 YOS = 86.7009 YOP = 39.1009 GAMZ = TWOPI * ( VZ*BZ*(1-iover) + 3.73607*48.3182*iover ) c JOE - case 1561 (~unstratified) VOS = 3.1 VOP = 3.1 ZOS =(168.3 + 5)*0.0 + 237.29*0.0 + 800. ZOP =(179.0 + 5)*0.0 + 234.27*0.0 + 800. YOS =( 49.6 )*0.0 + 28.175*0.0 + 22.75 YOP =( 25.4 )*0.0 + 6.165*0.0 - 22.75 GAMZ = TWOPI * ( VZ*BZ*(1-iover) + 4.55*26.46*iover ) GAMZ = TWOPI * ( VZ*BZ*(1-iover) + 3.73607*48.3182*iover ) GAMZ = TWOPI * ( VZ*BZ*(1-iover) + 4.55*26.46*iover ) c JOE - speed sign VSTART(1) = -VZ*(1-iover) - (VOP*0.+VZ*1.)*iover VSTART(1) = VZ*(1-iover) + (VOP*0.+VZ*1.)*iover VSTART(2) =(1000.0*0.+ZZ)*(1-ioverx) + ZOP*ioverx c JOE-THE PORT VORTEX IS ON THE LEFT AND THE STARBOARD VORTEX IS c ON THE RIGHT. VSTART(3) = (YZ - 0.5*BZ)*(1-ioverx) + YOP*ioverx c JOE-PORT AND STARBOARD CIRCULATION MAGNITUDES ARE EVOLVED BY THIS c CODE; I.E., THEY ARE BOTH POSITIVE. VSTART(4) = GAMZ*(1-iover) + (GAMZ*0.) * iover c JOE-VERTICAL VELOCITIES ARE NEGATIVE IF THE VORTICES ARE MOVING c DOWN. c JOE - speed sign VSTART(5) = -VZ*(1-iover) - (VOS*0.+VZ*1.)*iover VSTART(5) = VZ*(1-iover) + (VOS*0.+VZ*1.)*iover VSTART(6) =(1000.0*0.+ZZ)*(1-ioverx) + ZOS*ioverx VSTART(7) = (YZ + 0.5*BZ)*(1-ioverx) + YOS*ioverx VSTART(8) = GAMZ*(1-iover) + (GAMZ*0.) * iover ZAV0 = 0.5*(VSTART(2)+VSTART(6)) RHOINIT = FK(ZAV0,TDEGZ,RHO ,NZTDEG,IZTOPR) TDEGINIT = FK(ZAV0,TDEGZ,TDEG,NZTDEG,IZTOPT) VSTART(9) = RHOINIT c JOE particle density VSTART(10)= 1.0 C NT = NTOG NTIM = 0 if( debug_ige ) then write(6,*)" NT = ",NT end if ! ! ! T1,T2 are passed to RKDUMB2() and used to compute the ! timestep. T2 is not modified. the final time when ! exciting RKDUMB2() is identified by NTIM, which is the ! last timestep completed. ! !!!!! CALL RKDUMB2(VSTART,NVAR,T1,T2,NT,NTIM,TSS) CALL RKDUMB2(VSTART,NVAR,T1,T2,NT,NTIM,TL,TSS) IF(NTIM.GT.0) THEN IIM = 1 ELSE IIM = 0 NTIM = -NTIM ENDIF C DO 60 IT=1,NTIM SEPY(IT) = VV(7,IT) - VV(3,IT) SEPZ(IT) = VV(6,IT) - VV(2,IT) CCC SEPPS(IT) = VV(7,IT) - VV(3,IT) SEPPS(IT) = SQRT(SEPY(IT)**2 + SEPZ(IT)**2) HALFB = 0.5*SEPPS(IT) TWOPIB = TWOPI*SEPPS(IT) GAMHA(IT) = -VV(4,IT) GAMHB(IT) = VV(8,IT) GAMDA(IT) = TWOPIB*VV(5,IT) GAMDB(IT) = -TWOPIB*VV(1,IT) VVYP(IT) = VV(3,IT) VVYS(IT) = VV(7,IT) VVZP(IT) = VV(2,IT) VVZS(IT) = VV(6,IT) 60 CONTINUE C c JOE - hazard circ is written to traject.dat, so decent circ is c written here. if( debug2 ) then WRITE(2,607) (TT(IT),SEPY(IT),SEPZ(IT),SEPPS(IT), X GAMDA(IT),GAMDB(IT), IT=1,NTIM) 607 FORMAT(/1X,'SEPARATION HISTORY'//(6F12.3)) end if C DGMDTA = (GAMHA(NTIM)-GAMHA(NTIM-1))/(TT(NTIM)-TT(NTIM-1)) DGMDTB = (GAMHB(NTIM)-GAMHB(NTIM-1))/(TT(NTIM)-TT(NTIM-1)) if( debug2 ) then WRITE(2,611) TT(NTIM),DGMDTA,DGMDTB 611 FORMAT(/1X,'TIME,DGMDTA,DGMDTB AFTER PHASE I = ', X 3(1PE15.4)/) end if C NTMAX = NTOG IF(IIM.EQ.1) NTMAX = NTIM - 1 IF(IIM.EQ.0 .AND. NTIM.LE.NTMAX) NTMAX=NTIM C NTOUT = NTMAX c JOE write output here to DENSITY.DAT c what are we writing out? c 1. FR = 1/N* c 2. EDR* c 3. non-dim time c 4. non-dim ave height (ave of both vortices) c 5. non-dim density of vortex cell c 6. non-dim particle density of vortex cell c 7. non-dim circulation port vortex c 8. non-dim lateral position port vortex c 9. non-dim height port vortex c 10. non-dim velocity port vortex c 11. non-dim circulation starboard vortex c 12. non-dim lateral position starboard vortex c 13. non-dim height starboard vortex c 14. non-dim velocity starboard vortex do IT = 1 , NTMAX ZAV0 = 0.5*(VV(2,IT)+VV(6,IT)) BVFSQK = FK(ZAV0,TDEGZ,BVFSQ,NZTDEG,IZTOPR) QK = FK(ZAV0, QZ, Q, NZQ, NZQ) if (BVFSQK .lt. 0.) BVFSQK = 0.0 c Tom - We need to guard against N^2=0 here FROUT(IT) = VZ/BZ/sqrt(BVFSQK+1.0e-10) QSTAR(IT) = (QK*BZ)**0.333333/VZ TANANG = ABS(VV(6,IT)-VV(2,IT))/ABS(VV(7,IT)-VV(3,IT)) ANG = ATAN(TANANG) COSANGALL(IT) = COS(ANG) enddo !JW write normalized output to DENSITY.DAT if( debug4 ) then OPEN(4,FILE='DENSITY.DAT',STATUS='UNKNOWN') WRITE(4,62) (FROUT(IT) & ,QSTAR(IT) & ,TT(IT)*VZ/BZ !JW & ,TT(IT) & ,(VVZP(IT)+VVZS(IT))*0.5/BZ !JW & ,(VVZP(IT)+VVZS(IT))*0.5 & ,VV( 9,IT)/VV( 9,1) & ,VV(10,IT)/VV(10,1) !JW & ,VV( 9,IT) !JW & ,VV(10,IT) & ,VV( 4,IT)/GAMZ & ,VV( 3,IT)/BZ & ,VV( 2,IT)/BZ & ,VV( 1,IT)/VZ !JW & ,VV( 4,IT) !JW & ,VV( 3,IT) !JW & ,VV( 2,IT) !JW & ,VV( 1,IT) & ,VV( 8,IT)/GAMZ & ,VV( 7,IT)/BZ & ,VV( 6,IT)/BZ & ,VV( 5,IT)/VZ !JW & ,VV( 8,IT) !JW & ,VV( 7,IT) !JW & ,VV( 6,IT) !JW & ,VV( 5,IT) !JW & ,COSANGALL(IT) & ,IT=1,NTMAX) 62 FORMAT(14F10.3) CLOSE(4) !!-! & ,VV( 4,IT)/GAMZ !!-! & ,VV( 8,IT)/GAMZ !!-! & ,VV( 1,IT)/VZ !!-! & ,VV( 5,IT)/VZ endif ! set link time and short-wave instabilitiy time ! to pass to calling routine ATLD = TLD ATSSD = TSSD c do it=1, ntmax c write(87,87) tt(it), vv(1,it), vv(5,it) c87 format(1p,6e12.4) c end do c VP = VECY(1) speed, port vortex c ZP = VECY(2) height, port vortex c YP = VECY(3) horizontal position, port vortex c GAMP = VECY(4) hazard circulation, port vortex c VS = VECY(5) c ZS = VECY(6) c YS = VECY(7) c GAMS = VECY(8) c RHOC = VECY(9) c ptcls= VECY(10) C DO IT=1,NTMAX CCC GAOUT(IT) = GAMDA(IT) CCC GBOUT(IT) = GAMDB(IT) GAOUT(IT) = GAMHA(IT) GBOUT(IT) = GAMHB(IT) ENDDO C DO IT=1,NTMAX TOUT(IT) = TT(IT) YPOUT(IT) = VVYP(IT) ZPOUT(IT) = VVZP(IT) GPOUT(IT) = GAOUT(IT) YSOUT(IT) = VVYS(IT) ZSOUT(IT) = VVZS(IT) GSOUT(IT) = GBOUT(IT) ENDDO C CCC OPEN(3,FILE='TRAJEC.DAT',STATUS='UNKNOWN') CCC WRITE(3,61) (TT(IT),GAOUT(IT),VVYP(IT),VVZP(IT), CCC X GBOUT(IT),VVYS(IT),VVZS(IT),IT=1,NTMAX) CCC 61 FORMAT(F8.3,6F10.3) ! ! JOE - set IIM=0 if you want to skip IGE altogether ! IIM=0 IF(FLAG_IGE_NGE.eq.0) THEN IIM=0 if( debug6 ) then write(6,*)"skipping ige by setting FLAG_IGE_NGE=0 in viper.F" end if ENDIF ! this block has been modified to only write out the two primary ! vortices and not specify the images. images are taken into ! account in igsolv() ! IF(IIM.EQ.1) THEN if( debug6 ) then write(6,*)" descending into IGE ..." end if CCC GAM1 = GAMDA(NTIM) CCC GAM2 = GAMDB(NTIM) GAM1 = GAMHA(NTIM) GAM2 = GAMHB(NTIM) !! GAM3 = -GAM1 !! GAM4 = -GAM2 CCC DGFAC = DGMDT/GAM1 DGFACP = 0. DGFACS = 0. !JW these if-blocks are no longer ignored; see just below IF(GAM1 .NE. 0.) DGFACP = DGMDTA/GAM1 IF(GAM2 .NE. 0.) DGFACS = DGMDTB/GAM2 !JW DGFACP = DGMDTA/GAM1 !JW DGFACS = -DGMDTB/GAM2 ! why divide by GAM1,2 just to multiply by it again? ! DGAM(1,2) are time derivatives of GAM1,2 DGAM(1) = GAM1*DGFACP DGAM(2) = GAM2*DGFACS !! DGAM(3) = GAM3*DGFACP !! DGAM(4) = GAM4*DGFACS if( debug2 ) then WRITE(2,99) (DGAM(I),I=1,2) 99 FORMAT('DGAM = ',2(1PE15.4)/) !! WRITE(2,99) (DGAM(I),I=1,4) !! 99 FORMAT('DGAM = ',4(1PE15.4)/) end if Y13 = VVYP(NTIM) Y24 = VVYS(NTIM) CCC Z12 = VVZ(NTIM) CCC Z34 = -Z12 Z1 = VVZP(NTIM) Z2 = VVZS(NTIM) !! Z3 = -Z1 !! Z4 = -Z2 !! NPTS = 4 NPTS = 2 OPEN(1,FILE='CTRL.DAT',STATUS='UNKNOWN') WRITE(1,*) KMAX,DTFAC,T1IM,T2IM,EPS,H1,HMIN,NPTS,NXIIM WRITE(1,*) GAM1,Y13,Z1 WRITE(1,*) GAM2,Y24,Z2 !! WRITE(1,*) GAM3,Y13,Z3 !! WRITE(1,*) GAM4,Y24,Z4 CLOSE(1) C BZ2 = SEPPS(NTIM) VVAV = 0.5*(VV(1,NTIM) + VV(5,NTIM)) VVSB2 = VVAV**2/BZ2 C ZGE0 = ZGEFAC*BZ2 ZGE = ZGE0 IF(ZGE .GT. 0.99*ZIM) ZGE = 0.99*ZIM C GERAD0 = GERFAC*BZ2 IF(ZGE .LT. ZGE0) THEN GERAD = GERAD0*ZGE/ZGE0 RADFAC = GERAD/GERAD0 GAMFAC = GAMFAC*RADFAC**2 ELSE GERAD = GERAD0 ENDIF C GEARG = GEANG*PI/180. ZDOWN = GERAD*COS(GEARG) YOVER = GERAD*SIN(GEARG) C if( debug2 ) then WRITE(2,83) ZGE,GERAD,GAMFAC,ZDOWN,YOVER 83 FORMAT(/1X,'ZGE = ',1PE15.5/ X 1X,'GERAD,GAMFAC,ZDOWN,YOVER = ',4(1PE15.5)) end if C DGMDT = 0.5*(ABS(DGMDTA) + ABS(DGMDTB)) DELGAM = DGMDT if( debug2 ) then WRITE(2,91) 91 FORMAT(//1x,'CALL VIPER_IGE') end if if( debug6 ) then write(6,*) "calling viper_ige at t=",TT(NTIM) end if CALL VIPER_IGE(NTIM,TLAST,CEJECTA,CKE,R0p,R0s,DIS & ,TEJECT1,TEJECT2,ZEJECT,SECDK,SECDIFF) ! sketch code here: ! ! RKDUMB2 - when vortices drop low enough, mark as entering NGE, then exit ! ! VIPER_IGE - advect vortices, add when needed. C CC IF(TLAST .GE. 0. .AND. GAM(1) .LT. 0.) THEN ! ! this section adds a new vortex pair, then calls ge1pth(). ! originally the new pair's images were also added, but the ! core particle-advection algorithm has been modified to ! include image contributions with their associated primary ! vortex, so we only need to add the pair now. ! ! VIPER_IGE call ENDIF C FLAG = IFLAG NB = NTOUT ! 15m^2 * 0 + 20m^2 R15R15 = 15.*15.*0. + 20.*20.*1. + 25.*25.*0. ! if using non-dimensional units if (BZ.lt.2.) then R15R15 = 0.333333*0.333333 R15R15 = 0.444444*0.444444 endif DO I=1,NB TB(I) = TOUT(I) YPB(I) = YPOUT(I) ZPB(I) = ZPOUT(I) YSB(I) = YSOUT(I) ZSB(I) = ZSOUT(I) ! initialize port and starboard 5m-15m circulation values GPB(I) = GPOUT(I) GSB(I) = GSOUT(I) ! add to 5m-15m circulation if a 2ndary vortex is within 15m ! testing testing DO J=3,NPTS RPRP = (YPB(I)-YZPI(1,J,I))**2 + (ZPB(I)-YZPI(2,J,I))**2 RSRS = (YSB(I)-YZPI(1,J,I))**2 + (ZSB(I)-YZPI(2,J,I))**2 if (RPRP.lt.R15R15) then GPB(I) = GPB(I) + GAMI(J,I) endif if (RSRS.lt.R15R15) then GSB(I) = GSB(I) + GAMI(J,I) endif ENDDO ENDDO C NEWNB = NB !JW90 CONTINUE !JW IF(GPB(NEWNB) .GT. 0. .OR. GSB(NEWNB) .LT.0.) THEN !JW NEWNB = NEWNB - 1 !JW IF(NEWNB .GT. 1) GO TO 90 !JW ENDIF !JW IF(NEWNB .GT. 1) NB = NEWNB do I = 1 , NB ZAV0 = 0.5*(ZPB(I)+ZSB(I)) BVFSQK = FK(ZAV0,TDEGZ,BVFSQ,NZTDEG,IZTOPR) QK = FK(ZAV0, QZ, Q, NZQ, NZQ) if (BVFSQK .lt. 0.) BVFSQK = 0.0 c Tom - We need to guard against N^2=0 here FROUT(I) = VZ/BZ/sqrt(BVFSQK+1.0e-10) QSTAR(I) = (QK*BZ)**0.333333/VZ TANANG = ABS(ZPB(I)-ZSB(I))/ABS(YPB(I)-YSB(I)) ANG = ATAN(TANANG) COSANGALL(I) = COS(ANG) enddo ! write out the .DENSITY file for OGE+IGE if( debug4 ) then OPEN(4,FILE='Output/'//trim(basename)//'.DENSITY') WRITE(4,62) (FROUT(I) & ,QSTAR(I) & ,TB(I)*VZ/BZ & ,(ZPB(I)+ZSB(I))*0.5/BZ & ,VV( 9,I)/VV( 9,1) & ,VV(10,I)/VV(10,1) & ,GPB(I)/GAMZ & ,YPB(I)/BZ & ,ZPB(I)/BZ & ,VV( 1,I)/VZ & ,GSB(I)/GAMZ & ,YSB(I)/BZ & ,ZSB(I)/BZ & ,VV( 5,I)/VZ & ,I=1,NB) CLOSE(4) end if if( debug_ige ) then write(6,*)"NPTS,VZ,BZ=",NPTS,VZ,BZ end if ! write out the .DENSITY2 file for 2ndary IGE vortices if( debug4 ) then OPEN(4,FILE='Output/'//trim(basename)//'.DENSITY2') do J = 1 , NPTS do I = 1 , NXI ZAV0 = YZPI(2,J,I) BVFSQK = FK(ZAV0,TDEGZ,BVFSQ,NZTDEG,IZTOPR) QK = FK(ZAV0, QZ, Q, NZQ, NZQ) if (BVFSQK .lt. 0.) BVFSQK = 0.0 c Tom - We need to guard against N^2=0 here FROUT(I) = VZ/BZ/sqrt(BVFSQK+1.0e-10) QSTAR(I) = (QK*BZ)**0.333333/VZ enddo WRITE(4,63) (FROUT(I) & ,QSTAR(I) & ,XI( I)*VZ/BZ & ,YZPI(2,J,I)/BZ & ,0.0 & ,0.0 & ,GAMI( J,I)/GAMZ & ,YZPI(1,J,I)/BZ & ,YZPI(2,J,I)/BZ & ,0.0 & ,I=1,NXI) enddo CLOSE(4) 63 FORMAT(10F10.3) end if C if( debug2 ) then WRITE(2,*) 'NTOUT = ',NTOUT WRITE(2,999) 999 FORMAT(/' VIPER IS DONE.') CLOSE(2) end if CCC CLOSE(3) c JOE - write TSS to screen CCC write(6,*)"TSS = ",TSS CCC write(6,*)"T = ",TSS * BZ / VZ RETURN END C c JOE - added TSS to RKDUMB2() arg list !!!!! SUBROUTINE RKDUMB2(VSTART,NVAR,X1,X2,NSTEP,KIM,TSS) SUBROUTINE RKDUMB2(VSTART,NVAR,X1,X2,NSTEP,KIM,TL,TSS) c JOE - changed NVARMX to 10 to add particle density variable PARAMETER(NVARMX=10,NSTPMX=5001) COMMON /SARPB/C0,BZ,VZIN,R3,F0,VSB,VVSB,RVSB2,AFAC,DBFAC, X X,BRATIO,GAMAVB,ARGTB,RHOINIT,TDEGINIT COMMON /DBDTB/DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 REAL DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 COMMON /PATH/XX(NSTPMX),Y(NVARMX,NSTPMX),SEP(NSTPMX) COMMON /HB/H,HH,H6,XH 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 COMMON /STRAT/STRN DIMENSION VSTART(*),V(NVARMX),DV(NVARMX),VOUT(NVARMX) c JOE - define "time" REAL time, TSS, TL logical debug1, debug2, debug4, debug6, debug_ige common /debug/ debug1, debug2, debug4, debug6, debug_ige C DO 11 I=1,NVAR V(I) = VSTART(I) Y(I,1) = V(I) 11 CONTINUE XX(1) = X1 X = X1 H = (X2-X1)/FLOAT(NSTEP-1) IF(ZZ.LT.BZ) H = 0.1*H HH = 0.5*H H6=H/6. XH=X+HH C CCC KZERO = 0 KLIM = NSTEP - 1 KFIX48 = 0 CCC KFIX15 = 0 SEPTOT = BZ SEPOLD = BZ BDBDT = 0. DBDT1 = 0. BDBDT1 = 0. time = X1 !!!!! CALL DERIVS2(0.,V,DV,TSS,time) CALL DERIVS2(0.,V,DV,TL,TSS,time) if( debug2 ) then write(2,*) 'DIMT,SEPTOT,BRATIO = ',X,SEPTOT,BRATIO end if C*************************************** V2TSTB = ZZ C*************************************** DO 13 K=1,KLIM C WPOLD = WP WSOLD = WS time = XX(K) c if ( time*VSB .gt. 11.9 ) then if ( (time*VSB .gt. TSS*(6.9-2.6*STRN) .or. & (time*VSB .gt. TL+6)) .and. 1.eq.2 ) then vtfac = 1.0 !! write(6,*)" TSS = ",TSS !! write(6,*)" TL = ",TL !! write(6,*)" STRN= ",STRN KIM = -K VTFAC = 1. ! make sure K-1=0 is not specified for initial IGE runs IF (K.GT.1) THEN DO IV=1,NVAR Y(IV,K) = Y(IV,K-1) + VTFAC*(Y(IV,K)-Y(IV,K-1)) V(IV) = Y(IV,K) ENDDO XX(K) = XX(K-1) + VTFAC*(XX(K)-XX(K-1)) C DARG = VTFAC*H X = XX(K-1) CALL DERIVS2(DARG,V,DV,TL,TSS,time) SEP(K) = BZ/BRATIO ENDIF X = XX(K) RETURN endif IF (V(2).LE.ZIM .AND. V(6).LE.ZIM) THEN KIM = +K C**************************************** VTFAC = 1. V26MAX = AMAX1(V(2),V(6)) c IF(V26MAX .NE. V2TSTB) VTFAC = (ZIM-V2TSTB)/(V26MAX-V2TSTB) ! make sure K-1=0 is not specified for initial IGE runs IF (K.GT.1) THEN DO IV=1,NVAR Y(IV,K) = Y(IV,K-1) + VTFAC*(Y(IV,K)-Y(IV,K-1)) V(IV) = Y(IV,K) ENDDO XX(K) = XX(K-1) + VTFAC*(XX(K)-XX(K-1)) C DARG = VTFAC*H X = XX(K-1) !!!!! CALL DERIVS2(DARG,V,DV,TSS,time) CALL DERIVS2(DARG,V,DV,TL,TSS,time) SEP(K) = BZ/BRATIO ENDIF X = XX(K) C***************************************** RETURN ENDIF CCC IF(V(4) .LT. 0.01 .OR. V(8) .LT. 0.01) THEN CCC KIM = -(K-1) CCC XX(K) = 0. CCC RETURN CCC ENDIF CCC IF(DV(4) .GT. 0.01 .OR. DV(8) .GT. 0.01) THEN CCC KIM = -(K-1) CCC XX(K) = 0. CCC RETURN CCC ENDIF C C 16May09: Set cutoff when hazard circ <= 50 !!! IF((V(4) .LE. 25. .OR. V(8) .LE. 25.) .AND. KFIX48.EQ.0) THEN !!! KFIX48 = K !!! Y(4,K) = 25. !!! Y(8,K) = 25. !!! ENDIF !!! IF(KFIX48 .NE. 0) THEN !!! V(4) = 25. !!! V(8) = 25. !!! ENDIF IF((V(4) .LE. 0.) .AND. KFIX48.EQ.0) THEN KFIX48 = K Y(4,K) = 0. ENDIF IF((V(8) .LE. 0.) .AND. KFIX48.EQ.0) THEN KFIX48 = K Y(8,K) = 0. ENDIF c JOE-IF(KFIX48 .NE. 0) THEN c JOE- V(4) = 25. c JOE- V(8) = 25. c JOE-ENDIF C CCC IF((V(1) .GT. 0. .OR. V(5) .GT. 0.) .AND. KFIX15 .EQ. 0) THEN C FOLLOWING TWO IF TESTS ARE ALWAYS FALSE CCC IF((V(1) .GT. 0. .OR. V(5) .GT. 0.) .AND. KFIX15 .NE. 0) THEN CCC KFIX15 = K CCC Y(1,K) = 0. CCC Y(5,K) = 0. CCC ENDIF CCC IF(KFIX15 .NE. 0) THEN CCC V(1) = 0. CCC V(5) = 0. CCC ENDIF C c JOE - you passed TSS as a new last argument to DERIVS2(V,DV) time = XX(K) !!!!! CALL DERIVS2(0.,V,DV,TSS,time) CALL DERIVS2(0.,V,DV,TL,TSS,time) SEP(K) = BZ/BRATIO C C*************************************** V26MAX = AMAX1(V(2),V(6)) V2TSTB = V26MAX C*************************************** CALL RK42(V,DV,NVAR,VOUT,TL,TSS,time) c JOE - make sure V(10) - particle density - does not go negative if (V(10).lt.0.0) then V(10) = 0.0 endif C IF((X+H).EQ.X) THEN if( debug2 ) then WRITE(2,991) 991 FORMAT(/'X+H .EQ. X') CLOSE(2) end if CCC CLOSE(3) STOP ENDIF C ! ! if rebound has occured in NGE, reset tminp,tmins ! ! we are below the NGE altitude (ZIM=1.5*B0) if (V(2) .lt. ZIM) then ! the vertical velocity rebounded if ( (WPOLD.lt.0) .and. (WP.ge.0.) ) then ! tminp is still the default value of 9999. if (tminp.eq.9999.) then tminp = X - H*WPOLD/(WP-WPOLD) if( debug6 ) then write(6,*)"tminp = ",tminp end if endif endif endif ! ! we are below the NGE altitude (ZIM=1.5*B0) if (V(6) .lt. ZIM) then ! the vertical velocity rebounded if ( (WSOLD.lt.0) .and. (WS.ge.0.) ) then ! tmins is still the default value of 9999. if (tmins.eq.9999.) then tmins = X - H*WSOLD/(WS-WSOLD) if( debug6 ) then write(6,*)"tmins = ",tmins end if endif endif endif ! write(6,*)"ZP,ZS,ZIM = ",ZP,ZS,ZIM ! write(6,*)"tminp = ",tminp ! write(6,*)"tmins = ",tmins X = X + H XX(K+1) = X DO 12 I=1,NVAR Y(I,K+1) = VOUT(I) V(I) = VOUT(I) 12 CONTINUE ZP = V(2) YP = V(3) ZS = V(6) YS = V(7) ZPDOT = DV(2) YPDOT = DV(3) ZSDOT = DV(6) YSDOT = DV(7) SEPTOT = SQRT((YS-YP)**2 + (ZS-ZP)**2) CCC DBDT = (SEPTOT - SEPOLD) / H CCC BDBDT = (YP-YS)*(YPDOT-YSDOT) + (ZP-ZS)*(ZPDOT-ZSDOT) CCC SEPOLD = SEPTOT BRATIO = BZ / SEPTOT if( debug2 ) then write(2,*) 'DIMT,SEPTOT,BRATIO = ',X,SEPTOT,BRATIO end if C if( debug2 ) then write(2,881) X,(DV(I),I=1,4) 881 FORMAT('DIMT,DVP = ',5(1PE15.5)) write(2,882) X,(DV(I),I=5,8) 882 FORMAT('DIMT,DVS = ',5(1PE15.5)) end if C 13 CONTINUE C IF(V(2).LE.ZIM .AND. V(6).LE.ZIM) THEN KIM = NSTEP C**************************************** K = NSTEP VTFAC = 1. V26MAX = AMAX1(V(2),V(6)) c IF(V26MAX .NE. V2TSTB) VTFAC = (ZIM-V2TSTB)/(V26MAX-V2TSTB) DO IV=1,NVAR Y(IV,K) = Y(IV,K-1) + VTFAC*(Y(IV,K)-Y(IV,K-1)) V(IV) = Y(IV,K) ENDDO XX(K) = XX(K-1) + VTFAC*(XX(K)-XX(K-1)) DARG = VTFAC*H X = XX(K-1) time = XX(NSTEP) c JOE - you passed TSS to DERIVS2(V,DV) !!!!! CALL DERIVS2(DARG,V,DV,TSS,time) CALL DERIVS2(DARG,V,DV,TL,TSS,time) SEP(NSTEP) = BZ/BRATIO X = XX(K) C***************************************** ELSE KIM = -NSTEP ENDIF C RETURN END C C DERIVS C ROUTINE TO COMPUTE RHS OF MODEL DIFFEQ C c JOE - added TSS to DERIVS2() call !!!!! SUBROUTINE DERIVS2(DELT,VECY,DVECY,TSS,time) SUBROUTINE DERIVS2(DELT,VECY,DVECY,TL,TSS,time) REAL VECY(*),DVECY(*) C PARAMETER (NZMAX=401,NNUT=8) C CCC EXTERNAL FUNCD c JOE - what time is it? in "time" REAL time REAL N2STAR, N2STAR2 REAL ROLLUP, CI, CII REAL DVECY4, DVECY8 REAL TLREPORT, TSSREPORT, TL, TSS, TLD, TSSD, hw_fac INTEGER i_old C COMMON /KFLAGB/KUCR,KFLAG COMMON /VIPERB/GAMZ,BSV,FAC11,BETA1,ALPH1,FAC12,FAC13, X BETA2,BETA3,ALPH2,FAC21,FAC22,FAC23,THLF,FP3,CSH & , TLREPORT, TSSREPORT, TLD, TSSD, hw_fac, i_old COMMON /SARPB/C0,BZ,VZIN,R3,F0,VSB,VVSB,RVSB2,AFAC,DBFAC, X DIMT,BRATIO,GAMAVB,ARGTB,RHOINIT,TDEGINIT COMMON /DBDTB/DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 REAL DBDT,DBDT1,BDBDT,BDBDT1,SEPTOT,SEPOLD & ,WPOLD,WSOLD,WP,WS,tminp,tmins,tlag,tige1 COMMON /TSTARB/TS,BS,DLOGBS,TCFAC,IFIRST COMMON /THB/PI,TWOPI,R2PI,PIS4,PIS2,PI3S2, X THBOT,THTOP,THMIN1,THMIN2,THMAX1,THMAX2, X DELTH1,DELTH2,FACMIN,FACMAX,DELFAC 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 COMMON /ATMB/TDEGZ(NZMAX),TDEG(NZMAX) X ,QZ(NZMAX),Q(NZMAX) X ,UCRZ(NZMAX),UCR(NZMAX) X ,UHWZ(NZMAX),UHW(NZMAX) & ,NZTDEG,NZQ,NZUCR,NZUHW COMMON /ATMB2/DUCR(NZMAX),DDUCR(NZMAX) COMMON /RHOB/BVFSQ(NZMAX),RHO(NZMAX) COMMON /STRAT/STRN COMMON /CDB/REFAC,CDFAC,REHI,RELO,CDHI,CDLO,XNUT(NNUT),XNU(NNUT) COMMON /IZTOPB/ZLAST,QB,QB1,QB4 & ,IZTOPT,IZTOPQ,IZTOPU,IZTOPHW,NUTTOP COMMON /IZTPRB/IZTOPR COMMON /QB58B/QB5,QB8,BLAST,IZ2PDU COMMON /HB/H,HH,H6,XH COMMON /AIRPLANE/AIRPORTALT,RHOAIR,ACSPEEDI,TAN3 C JOE - common block of flag variables to control inclusion c of various terms in the evolution equations. C REAL FLAG_ENTR_CIR , FLAG_ENTR_CIR_EDR , FLAG_ENTR_CIR_N2 REAL FLAG_ENTR_MOM , FLAG_ENTR_MOM_EDR , FLAG_ENTR_MOM_N2 REAL FLAG_ENTR_DEN , FLAG_ENTR_DEN_EDR , FLAG_ENTR_DEN_N2 REAL FLAG_GRAV,FLAG_VCDOTSC,FLAG_VBDOTSB,FLAG_LINK REAL AEFF_S,FLAG_IGE_NGE,FLAG_DYNBUOY,FLAG_PARTICLE REAL AEFF_B,FLAG_CI INTEGER LINEAR COMMON /FLAGS/FLAG_ENTR_CIR,FLAG_ENTR_CIR_EDR,FLAG_ENTR_CIR_N2 & ,FLAG_ENTR_MOM,FLAG_ENTR_MOM_EDR,FLAG_ENTR_MOM_N2 & ,FLAG_ENTR_DEN,FLAG_ENTR_DEN_EDR,FLAG_ENTR_DEN_N2 X ,FLAG_GRAV,FLAG_VCDOTSC,FLAG_VBDOTSB,FLAG_LINK X ,AEFF_S,FLAG_IGE_NGE,FLAG_DYNBUOY,FLAG_PARTICLE X ,AEFF_B,FLAG_CI,LINEAR REAL RTE0,RTE1,RTE2,RTE3,FCIRBARO REAL RTE0BASE,RTE1BASE,RTE2BASE,RTE3BASE logical debug1, debug2, debug4, debug6, debug_ige common /debug/ debug1, debug2, debug4, debug6, debug_ige real w_avg C c TOM - Protected against overflow HSEC(X) = 2./(EXP(min(abs(X),35.0))+EXP(-min(abs(X),35.0))) C VP = VECY(1) ZP = VECY(2) YP = VECY(3) GAMP = VECY(4) VS = VECY(5) ZS = VECY(6) YS = VECY(7) GAMS = VECY(8) SEPSP = YS - YP ZAV = (ZS+ZP)*0.5 C RHOC = VECY(9) C UCRKP = FK(ZP,UCRZ,UCR,NZUCR,IZTOPU) UCRKS = FK(ZS,UCRZ,UCR,NZUCR,IZTOPU) c JOE - headwind UHWKP = FK(ZP ,UHWZ,UHW,NZUHW,IZTOPHW) UHWKS = FK(ZS ,UHWZ,UHW,NZUHW,IZTOPHW) UHWKA = FK(ZAV,UHWZ,UHW,NZUHW,IZTOPHW) QK = FK(ZAV,QZ,Q,NZQ,IZTOPQ) ARGT = DIMT + DELT C if( debug2 ) then write(2,991) argt,vp,vs,zp,zs,yp,ys,gamp,gams 991 format('argt,vp,vs,zp,zs,yp,ys,gamp,gams = ',9F12.3) end if C 16May09: Set cutoff when hazard circ <= 50 c JOE-IF(GAMP .LE. 25. .OR. GAMS .LE. 25. .OR. SEPSP .LE. 0.1*BZ) THEN !!! IF(GAMP .LE. 25. .OR. GAMS .LE. 25. ) THEN !!! DO I=1,8 !!! DVECY(I) = 0. !!! ENDDO !!! IF(KUCR.EQ.1) THEN !!! UCRKP = 0. !!! UCRKS = 0. !!! ENDIF !!! DVECY(3) = UCRKP !!! DVECY(7) = UCRKS !!! RETURN !!! ENDIF ! ! if primary circulation magnitude is zero (or less), ! vortices no longer exist - don't drift them with the ! mean cross wind. IF(GAMP .LE. 0. .AND. GAMS .LE. 0.) THEN DO I=1,8 DVECY(I) = 0. ENDDO IF(KUCR.EQ.1) THEN UCRKP = 0. UCRKS = 0. ENDIF DVECY(3) = UCRKP DVECY(7) = UCRKS RETURN ENDIF C SEPSQ = SEPSP*SEPSP CCC SEPTOT = SEPSP SEPTOTSQ = SEPSQ + (ZS-ZP)**2 SEPTOT = SQRT(SEPSQ + (ZS-ZP)**2) c JOE - this is the angle computed when delz and dely are both c positive. TANANG = ABS(ZS-ZP)/ABS(SEPSP) ANG = ATAN(TANANG) CCC DEGANG = ANG*180/3.14159 CCC print *, 'dimt,degang = ',dimt,degang SINANG = SIN(ANG) COSANG = COS(ANG) ABSVP = ABS(VP) ABSVS = ABS(VS) C SUMZ = ZP + ZS SUMZSQ = SUMZ*SUMZ c JOE - DENOM is for image vortices. that's why the sum c looks odd. don't worry about it for now. it's c for near-ground effects. DENOM = TWOPI*(SEPSQ + SUMZSQ) FOURPI = 2.*TWOPI CCC TWOPIB = TWOPI*SEPSP TWOPIB = TWOPI*SEPTOT c JOE - these gammas have signs and are estimated from V. c ZPCOR = -septot.Vp.(YS-YP)/((YS-YP)^2+(ZP+ZS)^2) c D's are for "descent" c JOE - speed sign c GAMDP = TWOPIB*VS c GAMDS = -TWOPIB*VP GAMDP = -TWOPIB*VS GAMDS = +TWOPIB*VP ZPCOR = GAMDS*SEPSP/DENOM ZSCOR = -GAMDP*SEPSP/DENOM YPCOR1 = GAMDS*SUMZ/DENOM YPCOR2 = GAMDP/(FOURPI*ZP) YSCOR1 = GAMDP*SUMZ/DENOM YSCOR2 = GAMDS/(FOURPI*ZS) CCC write(2,*) VP,ZP,YP,GAMP,VS,ZS,YS,GAMS CCC write(2,*) ZPCOR,YPCOR1,YPCOR2,ZSCOR,YSCOR1,YSCOR2 C ZAV = 0.5*(ZP + ZS) GAMAV = 0.5*(GAMP + GAMS) IZTOPT = NZTDEG TDEGK = FK(ZAV,TDEGZ,TDEG,NZTDEG,IZTOPT) c Tom - There used to be some statements here that attempted to c compute N^2. The equations did not appear correct to me and we c have already computed N^2 correctly in viper. We use the latter c values now. iztop = iztopt tdegk = fk(zav,Tdegz,tdeg ,nztdeg,iztop) bvfsqk = fk(zav,Tdegz,bvfsq,nztdeg,iztop) iztop = nztdeg bvfsqpsg = fk( zp,Tdegz,bvfsq,nztdeg,iztop)/grav iztop = nztdeg bvfsqssg = fk( zs,Tdegz,bvfsq,nztdeg,iztop)/grav ! N2STAR has sign N2STAR = BSV*BSV*BVFSQK !joe write(6,*)"N2STAR, BVFSQK = ",N2STAR,BVFSQK IF(BVFSQK .LT. 0.) BVFSQK = 0. !! try this: min( N2STAR ) = 0 because denver data looks strange !! still looks strange. perhaps when N2STAR is negative we have !! no real idea what the stratification is??? !! N2STAR = BSV*BSV*BVFSQK !!! IF(BVFSQK .LT. 0.) THEN !!! STRN = 0.0 !!! ELSE ! STRN has no sign; note: min(BVFSQK)=0 STRN = BSV*SQRT(ABS(BVFSQK)) !!! ENDIF STRN2 = STRN*STRN STRN2P = STRN**2.5 C QB1 = QB1 + (ZAV-ZLAST) * STRN2P QB4 = QB4 + (ZAV-ZLAST) * STRN2 QB5 = QB5 + (ZAV-ZLAST) * STRN2P QB8 = QB8 + (ZAV-ZLAST) * STRN2 ZLAST = ZAV C RHOK = FK(ZAV,TDEGZ,RHO ,NZTDEG,IZTOPR) C QK = FK(ZAV,QZ,Q,NZQ,IZTOPQ) EPSTAR = (BZ*QK)**R3/VZIN epstar = max( epstar , 1.0e-3 ) IF(EPSTAR .LT. 0.001) THEN TL = 9.0 ELSE IF(EPSTAR .LT. 0.0121 .AND. EPSTAR .GE. 0.001) THEN TL = 9.18 - 180.*EPSTAR ELSE IF(EPSTAR .LT. 0.2535 .AND. EPSTAR .GE. 0.0121) THEN C using Fred P's curve fit eqn. Oct 2008 C TL = RTSAFE(EPSTAR,FUNCD,1.,9.,1.E-5) TL = -1.5583*log(EPSTAR) + 0.1556 ELSE TL = F0*EPSTAR**(-0.75) ENDIF BIGT = ARGT*VSB C TSS IF(EPSTAR .LT. 0.25) THEN !! TSS = (0.57 - 1.27*ALOG(EPSTAR)) * EXP(-1.15*STRN) if (N2STAR .ge. 0.0) then TSS = (0.57 - 1.27*ALOG(EPSTAR))*(.16+.84*EXP(-3.*STRN)) else TSS = 1000. endif ELSE !! TSS = 0.825*EPSTAR**(-0.75) * EXP(-1.15*STRN) if (N2STAR .ge. 0.0) then TSS = 0.825*EPSTAR**(-0.75) *(.16+.84*EXP(-3.*STRN)) else TSS = 1000. endif ENDIF CCC IF(EPSTAR .LT. 0.3) THEN CCC TSS = -(1.27*ALOG(EPSTAR) + 0.57) * EXP(-1.15*STRN) CCC ELSE CCC TSS = -FP3 * EXP(-1.15*STRN) CCC ENDIF c *** Account for the time dialation due to headwind c print *, 'epstar, hw_fac, TL, TSS = ', epstar, hw_fac, TL, TSS TL = TL *hw_fac TSS = TSS*hw_fac if (TSSREPORT .eq. 0.0 .and. i_old.ne.1) then if (BIGT .gt. TSS) then c write(6,*)" t(shortwave) = ",TSS/VSB," sec, T* = ",TSS TSSD = TSS/VSB TSSREPORT = 1.0 endif endif if (TLREPORT .eq. 0.0 .and. i_old.ne.1) then if (BIGT .gt. TL + 0.1) then c write(6,*)" t(link) = ",TL/VSB," sec, T* = ",TL TLD = TL/VSB TLREPORT = 1.0 endif endif if (BIGT .gt. TL+0.1 .and. BIGT .lt. TL+3.1) then FFL = 1.0 + 0.35*TLREPORT*( 1. + 2.*EPSTAR ) & /(1.+0.0*(TL-TSS)**2.) else FFL = 1.0 endif TSSFAC = -(1.27*ALOG(EPSTAR) + 0.57) IF(EPSTAR .GT. 0.25) X TSSFAC = -(1.27*ALOG(EPSTAR) + 2.554*ALOG(1+EPSTAR)) c JOE - remove stratification effect on circulation ! JOE - use TSS above (don't overwrite TSS here) ! JOE TSS = TSSFAC*EXP(-1.15*STRN) !!!! TSS = TSSFAC*1.0 C C RAMPING FUNCTION, FT IF(THLF .EQ. 0. .AND. GAMAV .LE. 0.5*GAMZ) THEN RATIO = (ARGT-ARGTB)/(GAMAV-GAMAVB) THLF = ARGT + (0.5*GAMZ - GAMAV)*RATIO THLF = THLF*VSB ENDIF IF(THLF .EQ. 0.) THEN FT = 1. ELSE FT = 1. - 0.4*(BIGT-THLF) IF(FT .GT. 1.) FT = 1. IF(FT .LT. 0.) FT = 0. CCC FT = AMAX1(0.,1. - 0.4*(BIGT-THLF)) ENDIF IF (BIGT .LT. TL) THEN RTEFACT = 1.0 RTEFACV = 1.0 ELSE if ( i_old .eq. 1) then RTEFACT = 2.0 RTEFACV = 2.0 else RTEFACT = 1.+0.00*23.3/0.5 + 0.0 RTEFACV = 1.+0.00*1.67/0.5 + 0.0 endif ENDIF !bug IF (BIGT .LT. 2.*TSS) THEN IF (BIGT .LT. TL) THEN RTEFACG = 1.0 ELSE if ( i_old .eq. 1) then RTEFACG = 1.2 else RTEFACG = 1.+0.00*1.0 + 0.0 endif ENDIF C c JOE - cross wind and head wind are interpolated above, and c their values have not changed by this point. sanity c check debugging now included just as comments here. c !!!! UCRKP0 = UCRKP !!!! UCRKS0 = UCRKS !!!! UCRKP = FK(ZP,UCRZ,UCR,NZUCR,IZTOPU) !!!! UCRKS = FK(ZS,UCRZ,UCR,NZUCR,IZTOPU) c JOE - headwind !!!! UHWKP0 = UHWKP !!!! UHWKS0 = UHWKS !!!! UHWKP = FK(ZP,UHWZ,UHW,NZUHW,IZTOPHW) !!!! UHWKS = FK(ZS,UHWZ,UHW,NZUHW,IZTOPHW) !!!! if ( !!!! & (UCRKP.ne.UCRKP0) !!!! & .or. (UCRKS.ne.UCRKS0) !!!! & .or. (UHWKP.ne.UHWKP0) !!!! & .or. (UHWKS.ne.UHWKS0) !!!! & ) then !!!! write(6,*)" UCRKP,S,UHWKP,S changed. huh?" !!!! stop !!!! endif C DUKP = FK(ZP,UCRZ,DUCR,NZUCR,IZ2PDU) DUKS = FK(ZS,UCRZ,DUCR,NZUCR,IZ2PDU) DDUKP = FK(ZP,UCRZ,DDUCR,NZUCR,IZ2PDU) DDUKS = FK(ZS,UCRZ,DDUCR,NZUCR,IZ2PDU) c C c JOE - DBDT is not here. it is in RKDUMB2() CCC DBDT = (SEPTOT-BLAST)/H CCC BLAST = SEPTOT C C ************************************************************* ! hard-wired roll-up time - set Aeff(S), Aeff(B)), A_DB ! to zero if ROLLUP=1 (during roll-up time). if (time .lt. 20.0*BZ*ACSPEEDI) then ROLLUP = 1.0 else ROLLUP = 0.0 endif ! ROLLUP = 0.0 c JOE - a minus sign now appears for FACST. the equations c were also modified, so the net result is the same. c i've hard coded FACPO=1 and FACST=-1 so that these c terms are always present in the equations when not c in the rollup phase.. FACPO = 1. * (1.0-ROLLUP) FACST = -1. * (1.0-ROLLUP) C ************************************************************** C FACMAX = AMAX1(EPSTAR,0.08) HSFAC = HSEC(BETA1*(BIGT - TL - ALPH1)) c JOE - remove stratification effect on B23FAC B23FAC = BETA2 + BETA3*STRN**4 !!!! B23FAC = BETA2 HSFAC2 = HSEC(B23FAC*(BIGT - TSS) - ALPH2) C DBDT = DBDT1 BDBDT = BDBDT1 IF(ARGT.GT.ARGTB) THEN DBDT=(SEPTOT-SEPOLD)/(DBLE(ARGT)-DBLE(ARGTB)) BDBDT = DBDT * (SEPTOT+SEPOLD)*0.5 ENDIF C VAV = 0.5*(VP+VS) GAMAV = 0.5*(ABS(GAMP)+ABS(GAMS)) GAMAXP = AMAX1(ABS(GAMP),TWOPI*ABS(VP)*BZ) GAMAXS = AMAX1(ABS(GAMS),TWOPI*ABS(VS)*BZ) GAMAXAV = 0.5*(GAMAXP+GAMAXS) RRAT = RHOK/RHOC gamma_gas1 = 1./(gamma_gas-1.) ! JOE - in rollup, core temperature and density are set by ambient ! if (ROLLUP .EQ. 1.0) then ! RHOINIT = RHOK ! TDEGINIT = TDEGK ! endif RHOCORE = RHOINIT*(TDEGK/TDEGINIT)**gamma_gas1 C2RSB0SQ = 0.2 c JOE - C2RSB0SQ = C2(R/b0)^2 N2STAR2 = -(1.-RRAT)*GRAV*BSV*BSV/BZ * 3. N2STAR2 = N2STAR ! change N2STAR2's stratification dependence to the ! difference between the core and the oval buoyancy ! N2STAR2 = ABS(1.-RHOINIT/RHOC)*GRAV*BSV*BSV/BZ ! ! this should not be ABS(). using ABS() is producing a stable ! buoyancy effect for convectively unstable situations. test ! the sign against the nwra lab data. ! N2STAR2 = ABS(1.-RHOCORE/RHOC)*GRAV*BSV*BSV/BZ N2STAR2 = (1.-RHOCORE/RHOC)*GRAV*BSV*BSV/BZ ! N2STAR2 = N2STAR !! write(6,*)"N2STAR,N2STAR2,RHOC,RHOCORE = ", !! & N2STAR,N2STAR2,RHOC,RHOCORE ! get the EDR*=0, Lab-tuned model coefficients if ( i_old .eq. 1) then RTE0BASE = 0.0 RTE1BASE = 0.0 RTE2BASE = 0.0 else if (n2star2 .lt. 0.) then call rte_base(-STRN , RTE0BASE , RTE1BASE , RTE2BASE ) else call rte_base( STRN , RTE0BASE , RTE1BASE , RTE2BASE ) endif endif !! 1/(entrainment time) for momentum RTE1 = FLAG_ENTR_MOM + EPSTAR * FLAG_ENTR_MOM_EDR & + N2STAR2* FLAG_ENTR_MOM_N2 & + RTE1BASE RTE1 = RTE1 * RTEFACV !! RTE1 = RTE1 - RTEFACV !! 1/(entrainment time) for density RTE2 = FLAG_ENTR_DEN + EPSTAR * FLAG_ENTR_DEN_EDR & + N2STAR2* FLAG_ENTR_DEN_N2 & + RTE2BASE RTE2 = RTE2 * RTEFACT !! RTE2 = RTE2 - RTEFACT !! 1/(entrainment time) for particles RTE3 = FLAG_PARTICLE + EPSTAR * FLAG_ENTR_DEN_EDR & + N2STAR2* FLAG_ENTR_DEN_N2 & + RTE2BASE ! RTE3 = RTE2 RTE3 = RTE3 * RTEFACT !! RTE3 = RTE3 - RTEFACT c JOE - add EPSTAR dependence to RTE0 (0.272727 = 0.09/0.33) !! 1/(entrainment time) coefficient for circulation RTE0 = FLAG_ENTR_CIR + EPSTAR * FLAG_ENTR_CIR_EDR & + N2STAR2* FLAG_ENTR_CIR_N2 & + RTE0BASE RTE0 = RTE0 * RTEFACG !! RTE0 = RTE0 - RTEFACG RTE0 = RTE0 * ( 1.0 + 0.3 * TSSREPORT * (1.0+3.*EPSTAR) ) RTE1 = RTE1 * ( 1.0 - 0.3 * TSSREPORT * (1.0+3.*EPSTAR) ) !! write(6,*)" RTE0 = ",RTE0 c RTE0 = 1.0/(9.5*(gamav/gamz)**2) ! testing testing ! RTE0 = 0.0 ! RTE1 = 0.0 ! RTE2 = 0.0 ! RTE3 = 0.0 ! JOE - circulation evolution equations baroclinic term FCIRBARO = AEFF_B * (1.0-ROLLUP) FGRAV = FLAG_GRAV FVCDOTSC = FLAG_VCDOTSC FVBDOTSB = FLAG_VBDOTSB FLINK = FLAG_LINK c JOE - speed sign FUPP = AEFF_S FUPP = -AEFF_S * (1.0 + 0.0*EPSTAR) if (abs(FUPP).lt.0.0001) then FUPP = 0.0 endif SFUPP = 0.0 if (FUPP .lt. 0.) then SFUPP = -1.0 endif if (FUPP .gt. 0.) then SFUPP = 1.0 endif ADB = FLAG_DYNBUOY * (1.0-ROLLUP) ! factor for internal-structure term in the momentum equation CI = 1.0 CI = 0.1 CI = 0.2 CI = FLAG_CI CII = 1.0/CI - 1.0 FGE = FLAG_IGE_NGE !23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 ! ! how to improve this model? ! !?-1. T* is time to linking. there is a model for Gamma vs time that ! includes a rapid change at T*. what is the basis for this? ! (aircraft data supporting this?) ! ! 2. in the model T*'s impact is included as a time rate of change. ! rebounding may result because Gamma reverses (sign). is this the ! case? ! ! 3. once linking happens, then the nature of the prediction must ! change - topology of wake and quantification of the hazard ! changes fundamentally. do we have enough information to ! characterize it? ! (from lab data - is it sufficient?) ! ! 4. the lab experiments certainly demonstrate linking for the ! particle-entrainment measurements, yet they retain a fairly ! smooth descent rate. ! !?-5. T* is assigned a dependence on N (brunt-vaisala frequency). how ! is this dependence assigned? is there a better way to do it? ! T* dependence looks good for TASS sims. are these times ! consistent with lab results? ! ! -6. of the cases you tested already, what is T*? get it on the plots. ! ! -7. what does the literature say about entrainment of density versus ! momentum or particles? ! ! 8. does your current model (based on entrainment) include a sensible ! description of the ambient characteristics with which mixing is ! occurring? A: not sure. the time constants selected for density, ! particles, and momentum are rather disparate. why? it could ! result from differences in the s_e and s_d factors and the ! associated field gradients in the entrainment time t_e's for these ! different processes. e.g., for momentum, the lateral gradient ! depends on the dependence of the vertical velocity as a funciton ! of distance from the cell. so the entrained momentum will be ! greater than zero if the entrained fluid is also moving down as ! the cell is. the same goes for density/temperature. if we assume ! the wake flow is self-similar, then these factors can be fixed by ! the geometry of the wake. ! ! i could spend a lot of time setting the constants and getting them ! just right, then decide i need to change how the hazard circulation ! is handled, only to be back to setting the constants again. but it's ! not like setting the constants takes a lot of time, so quit whining ! and just get some results. ! ! c JOE - notes on what the variables mean: c c VP = VECY(1) speed, port vortex c ZP = VECY(2) height, port vortex c YP = VECY(3) horizontal position, port vortex c GAMP = VECY(4) hazard circulation, port vortex c VS = VECY(5) c ZS = VECY(6) c YS = VECY(7) c GAMS = VECY(8) c RHOC = VECY(9) c ptcls= VECY(10) c SEPSP = YS - YP starboard-port horizontal separation FNEW = 1.000 FNEW2= 0.000 if( debug2 ) then write(2,*) 'dimt,rhoc,rhok,rrat,gamav,vav = ', X dimt,rhoc,rhok,rrat,gamav,vav end if ! vertical velocities: c JOE - the 2nd term takes care of a near-ground effect that i am c not concerned with. leave it alone. CCC DVECY(2) = d(ZP)/dt c JOE - speed sign c DVECY(2) = VP*( COSANG*1.0 + 0.0 )*FFL + ZPCOR * FGE DVECY(2) = -VP*( COSANG*1.0 + 0.0 )*FFL + ZPCOR * FGE DVECY(2) = ( DVECY(2) - UHWKP*TAN3 ) * ( 1 - UHWKP*ACSPEEDI ) WP = DVECY(2) CCC DVECY(6) = d(ZS)/dt c JOE - speed sign c DVECY(6) = VS*( COSANG*1.0 + 0.0 )*FFL + ZSCOR * FGE DVECY(6) = -VS*( COSANG*1.0 + 0.0 )*FFL + ZSCOR * FGE DVECY(6) = ( DVECY(6) - UHWKS*TAN3 ) * ( 1 - UHWKS*ACSPEEDI ) WS = DVECY(6) c write(86,86) time, cosang, vp, vs, dvecy(2), dvecy(6), UHWKS*TAN3 c86 format(f8.3,1p,6e12.4) C c JOE - DVECY(4) = d(GAMP)/dt c JOE - comment out old terms c JOE - tie entrainment velocity to GAMAV, not GAMP or GAMS ! LINEAR=0, use entrainment model if (LINEAR.eq.0) then DVECY(4) = c 1st term = linking X -FAC21*FT*B23FAC*HSFAC2*HSFAC2 * FLINK c 2nd term = turbulent diffusion (entrainment should be ok) X -FAC22*FACMAX*GAMP * 0.0 c 3rd term = stratification X +FAC23*QB4 * 0.0 c X -RTE0*GAMP*GAMP/(TWOPI*BZ*BZ) * 1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMP*GAMAV /(TWOPI*BZ*BZ) *1.0 X -(RTE0-1*RTE2*(1.-RRAT))*GAMP*GAMAXP/(TWOPI*BZ*BZ) *1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMP*GAMP /(TWOPI*BZ*BZ) *1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMP*GAMZ /(TWOPI*BZ*BZ) *1.0 ! LINEAR=1, use linear model else ! what time is it? ! rollup phase: T*<1 if (time .le. BSV) then ! linear model: rollup DVECY(4) = -0.1020 else ! port vortex OGE phase: port vortex above ZIM (Zp>ZIM) if (ZP .gt. ZIM) then ! linear model: OGE DVECY(4) = -0.1020-0.687*EPSTAR & +0.048*sqrt(amax1(0.,N2STAR)) else if (time .lt. tminp+tlag) then ! linear model: NGE DVECY(4) = -0.0715-0.108*EPSTAR & +0.00178*sqrt(amax1(0.,N2STAR)) elseif (time .lt. tminp+tlag+tige1) then ! linear model: IGE1 DVECY(4) = -0.0591-0.915*EPSTAR & -0.119*sqrt(amax1(0.,N2STAR)) else ! linear model: IGE2 DVECY(4) = -0.0472-0.40*EPSTAR & +0.0173*sqrt(amax1(0.,N2STAR)) endif endif endif DVECY(4) = min( DVECY(4) * GAMZ*GAMZ/(TWOPI*BZ*BZ) , 0.0 ) endif ! add shear and shear-gradient effects DVECY(4) = DVECY(4) !JAW X +0.5*C2RSB0SQ*FACPO*(BZ/SEPTOT)*GAMS*DDUKP*BZ*COSANG X -FUPP *FACPO*BZ*BZ*VP*DDUKP*COSANG * 1. c JOE - this is the accidentally missing circulation term that goes with U" [1] !!!!!X -FUPP *FACPO*DUKP*RTE2*(1.-RRAT)*GAMAV /TWOPI *FNEW X -FUPP *FACPO*DUKP*RTE2*(1.-RRAT)*GAMAXP/TWOPI *FNEW !!!!!X -FUPP *FACPO*DUKP*RTE2*(1.-RRAT)*GAMP /TWOPI *FNEW !!!!!X -FUPP *FACPO*DUKP*RTE2*(1.-RRAT)*GAMZ /TWOPI *FNEW X +FGRAV*GRAV*(1.-RRAT)*FCIRBARO*BZ*COSANG c JOE - this is the accidentally missing circulation term that goes with rs [4] X -FACPO*(UCRKP+DUKP*SEPTOT*0.5*COSANG)*SEPTOT*0.5*COSANG !!!!!X *RTE1*(1.-RRAT)*GAMAV /(BZ*BZ) *FNEW2 X *RTE1*(1.-RRAT)*GAMAXP/(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMP /(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMZ /(BZ*BZ) *FNEW2 !JAW X + DCPDT_FIT ! store internal-structure momentum term for later use DVECY4 = DVECY(4) X - FACPO*BZ*BZ*VP*DDUKP*COSANG * 1. *( FUPP * CII ) !!!!!X - FACPO*DUKP*RTE2*(1.-RRAT)*GAMAV /TWOPI *FNEW*( FUPP * CII ) X - FACPO*DUKP*RTE2*(1.-RRAT)*GAMAXP/TWOPI *FNEW*( FUPP * CII ) !!!!!X - FACPO*DUKP*RTE2*(1.-RRAT)*GAMP /TWOPI *FNEW*( FUPP * CII ) !!!!!X - FACPO*DUKP*RTE2*(1.-RRAT)*GAMZ /TWOPI *FNEW*( FUPP * CII ) ! add headwind effect DVECY(4) = DVECY(4) * ( 1 - UHWKP*ACSPEEDI ) DVECY4 = DVECY4 * ( 1 - UHWKP*ACSPEEDI ) c JOE - DVECY(8) = d(GAMS)/dt c JOE - comment out old terms c JOE - tie entrainment velocity to GAMAV, not GAMP or GAMS ! LINEAR=0, use entrainment model if (LINEAR.eq.0) then DVECY(8) = c 1st term = linking X -FAC21*FT*B23FAC*HSFAC2*HSFAC2 * FLINK c 2nd term = turbulent diffusion (entrainment should be ok) X -FAC22*FACMAX*GAMS * 0.0 c 3rd term = stratification X +FAC23*QB8 * 0.0 c X -RTE0*GAMS*GAMS/(TWOPI*BZ*BZ) * 1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMS*GAMAV /(TWOPI*BZ*BZ) *1.0 X -(RTE0-1*RTE2*(1.-RRAT))*GAMS*GAMAXS/(TWOPI*BZ*BZ) *1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMS*GAMS /(TWOPI*BZ*BZ) *1.0 !!!!!X -(RTE0-1*RTE2*(1.-RRAT))*GAMS*GAMZ /(TWOPI*BZ*BZ) *1.0 ! LINEAR=1, use linear model else ! what time is it? ! rollup phase: T*<1 if (time .le. BSV) then ! linear model: rollup DVECY(8) = -0.1020 else ! starboard vortex OGE phase: port vortex above ZIM (ZS>ZIM) if (ZS .gt. ZIM) then ! linear model: OGE DVECY(8) = -0.1020-0.687*EPSTAR & +0.048*sqrt(amax1(0.,N2STAR)) else if (time .lt. tmins+tlag) then ! linear model: NGE DVECY(8) = -0.0715-0.108*EPSTAR & +0.00178*sqrt(amax1(0.,N2STAR)) elseif (time .lt. tmins+tlag+tige1) then ! linear model: IGE1 DVECY(8) = -0.0591-0.915*EPSTAR & -0.119*sqrt(amax1(0.,N2STAR)) else ! linear model: IGE2 DVECY(8) = -0.0472-0.40*EPSTAR & +0.0173*sqrt(amax1(0.,N2STAR)) endif endif endif DVECY(8) = min( DVECY(8) * GAMZ*GAMZ/(TWOPI*BZ*BZ) , 0.0 ) endif ! add shear and shear-gradient effects DVECY(8) = DVECY(8) !JAW X +0.5*C2RSB0SQ*FACST*(BZ/SEPTOT)*GAMP*DDUKS*BZ*COSANG X -FUPP *FACST*BZ*BZ*VS*DDUKS*COSANG * 1. c JOE - this is the accidentally missing circulation term that goes with U" [1] !!!!!X -FUPP *FACST*DUKS*RTE2*(1.-RRAT)*GAMAV /TWOPI *FNEW X -FUPP *FACST*DUKS*RTE2*(1.-RRAT)*GAMAXS/TWOPI *FNEW !!!!!X -FUPP *FACST*DUKS*RTE2*(1.-RRAT)*GAMS /TWOPI *FNEW !!!!!X -FUPP *FACST*DUKS*RTE2*(1.-RRAT)*GAMZ /TWOPI *FNEW X +FGRAV*GRAV*(1.-RRAT)*FCIRBARO*BZ*COSANG c JOE - this is the accidentally missing circulation term that goes with rs [4] X -FACST*(UCRKS+DUKS*SEPTOT*0.5*COSANG)*SEPTOT*0.5*COSANG !!!!!X *RTE1*(1.-RRAT)*GAMAV /(BZ*BZ) *FNEW2 X *RTE1*(1.-RRAT)*GAMAXS/(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMS /(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMZ /(BZ*BZ) *FNEW2 !JAW X + DCSDT_FIT ! store internal-structure momentum term for later use DVECY8 = DVECY(8) X - FACST*BZ*BZ*VS*DDUKS*COSANG * 1. *( FUPP * CII ) !!!!!X - FACST*DUKS*RTE2*(1.-RRAT)*GAMAV /TWOPI *FNEW*( FUPP * CII ) X - FACST*DUKS*RTE2*(1.-RRAT)*GAMAXS/TWOPI *FNEW*( FUPP * CII ) !!!!!X - FACST*DUKS*RTE2*(1.-RRAT)*GAMS /TWOPI *FNEW*( FUPP * CII ) !!!!!X - FACST*DUKS*RTE2*(1.-RRAT)*GAMZ /TWOPI *FNEW*( FUPP * CII ) ! add headwind effect DVECY(8) = DVECY(8) * ( 1 - UHWKS*ACSPEEDI ) DVECY8 = DVECY8 * ( 1 - UHWKS*ACSPEEDI ) c JOE - zeroed out turbulence and instability contribs to descent c VECY(1) = VP = speed of port vortex c 1st term = linking DVECY(1) = !!! X +FAC11*HSFAC*HSFAC * 0.0 c 2nd term = turb diffusion (entrainment should be ok) !!! X -FAC12*FACMAX*VP * 0.0 c JOE - you need to uncomment this line and comment c out the one below it to run Fred's model c 3rd term = buoyancy !!! X -FAC13*QB1 c JOE - i normalized with GAMS rather than GAMAV, and i added a cosang c factor to gravity since it only acts in the vertical, and this c code advances the speed. c JOE - tie entrainment velocity to GAMAV, not GAMP or GAMS !!! X -RTE1*VAV*RRAT - GRAV*(1.-RRAT) !!!!!X -VP *GAMAV /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 X -VP *GAMAXP/(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 !!!!!X -VP *GAMS /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 !!!!!X -VP *GAMZ /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 c JOE - this is the accidentally missing momentum term that goes with U [3] !!!!!X -UCRKP*SINANG*GAMAV /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW X -UCRKP*SINANG*GAMAXP/(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW !!!!!X -UCRKP*SINANG*GAMP /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW !!!!!X -UCRKP*SINANG*GAMZ /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW c JOE - speed sign c X -GRAV*(1.-RRAT)*COSANG * FGRAV X +GRAV*(1.-RRAT)*COSANG * FGRAV c JOE - this is the accidentally missing momentum term that goes with U'[2] X +DUKP*VP*COSANG*SINANG *1.0 *FNEW c this term = related to shear effects and varying separation !!! X -CSH*FACPO*SEPSP*VP*COSANG*DDUKP c JOEX +VP*( DVECY(8)/VECY(8) * FVCDOTSC c JOEX +ABS(GAMS)/(TWOPI*SEPTOT)*( DVECY(8)/VECY(8) * FVCDOTSC X +ABS(GAMS)/(TWOPI*SEPTOT)*( DVECY8 /VECY(8) * FVCDOTSC X - BDBDT/SEPTOTSQ * FVBDOTSB ) * CI !JAW X +VP*( DVECY(8)/VECY(8) - BDOTB_FIT/BB_FIT ) c JOE - this is the accidentally missing momentum term that goes with rs [5] X -(DUKP+DDUKP*SEPTOT*0.25*COSANG) X *SEPTOT*0.5*COSANG*SINANG !!!!!X *RTE1*(1.-RRAT)*GAMAV /(BZ*BZ) *FNEW2 X *RTE1*(1.-RRAT)*GAMAXP/(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMP /(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMZ /(BZ*BZ) *FNEW2 DVECY(1) = DVECY(1) * ( 1 - UHWKP*ACSPEEDI ) c JOE - zeroed out turbulence and instability contribs to descent c VECY(5) = VS = speed of starboard vortex c 1st term = linking DVECY(5) = !!! X +FAC11*HSFAC*HSFAC * 0.0 c 2nd term = turb diffusion (entrainment should be ok) !!! X -FAC12*FACMAX*VS * 0.0 c JOE - you need to uncomment this line and comment c out the one below it to run Fred's model c 3rd term = buoyancy !!! X -FAC13*QB5 c JOE - see cmments above for DVECY(1) c JOE - tie entrainment velocity to GAMAV, not GAMP or GAMS !!! X -RTE1*VAV*RRAT - GRAV*(1.-RRAT) !!!!!X -VS *GAMAV /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 X -VS *GAMAXS/(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 !!!!!X -VS *GAMS /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 !!!!!X -VS *GAMZ /(GAMZ*BSV)*(RTE1-1*RTE2*(1.-RRAT)) * 1.0 c JOE - this is the accidentally missing momentum term that goes with U [3] !!!!!X -UCRKS*SINANG*GAMAV /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW X -UCRKS*SINANG*GAMAXS/(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW !!!!!X -UCRKS*SINANG*GAMS /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW !!!!!X -UCRKS*SINANG*GAMZ /(GAMZ*BSV)*(RTE1-RTE2)*(1.-RRAT) *FNEW c JOE - speed sign c X -GRAV*(1.-RRAT)*COSANG * FGRAV X +GRAV*(1.-RRAT)*COSANG * FGRAV c JOE - this is the accidentally missing momentum term that goes with U'[2] X +DUKS*VS*COSANG*SINANG *1.0 *FNEW c this term = related to shear effects and varying separation !!! X -CSH*FACST*SEPSP*VS*COSANG*DDUKS c JOEX +VS*( DVECY(4)/VECY(4) * FVCDOTSC c JOEX +ABS(GAMS)/(TWOPI*SEPTOT)*( DVECY(4)/VECY(4) * FVCDOTSC X +ABS(GAMS)/(TWOPI*SEPTOT)*( DVECY4 /VECY(4) * FVCDOTSC X - BDBDT/SEPTOTSQ * FVBDOTSB ) * CI !JAW X +VS*( DVECY(4)/VECY(4) - BDOTB_FIT/BB_FIT ) c JOE - this is the accidentally missing momentum term that goes with rs [5] X -(DUKS+DDUKS*SEPTOT*0.25*COSANG) X *SEPTOT*0.5*COSANG*SINANG !!!!!X *RTE1*(1.-RRAT)*GAMAV /(BZ*BZ) *FNEW2 X *RTE1*(1.-RRAT)*GAMAXS/(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMS /(BZ*BZ) *FNEW2 !!!!!X *RTE1*(1.-RRAT)*GAMZ /(BZ*BZ) *FNEW2 DVECY(5) = DVECY(5) * ( 1 - UHWKS*ACSPEEDI ) c JOE - UCRKP is the horizontal drift velocity at the altitude of c the port vortex. the 2nd and 3rd terms are related to IGE. c the ADB term results from dynamical buoyancy !JW ADB = 3.0 * FDB *** ADB is now set to be equal to FLAG_DYNBUOY CCC DVECY(3) = UCRKP SAMEP = 1.0*UCRKP + ( YPCOR1 + YPCOR2 )*FGE + ADB*GAMP*BVFSQPSG IF(ZS.GE.ZP) DVECY(3) = SAMEP + ABSVP*SINANG * 1.0*FFL IF(ZS.LT.ZP) DVECY(3) = SAMEP - ABSVP*SINANG * 1.0*FFL DVECY(3) = DVECY(3) * ( 1 - UHWKP*ACSPEEDI ) CCC DVECY(7) = UCRKS SAMES = 1.0*UCRKS + ( YSCOR1 + YSCOR2 )*FGE - ADB*GAMS*BVFSQSSG IF(ZS.GE.ZP) DVECY(7) = SAMES + ABSVS*SINANG * 1.0*FFL IF(ZS.LT.ZP) DVECY(7) = SAMES - ABSVS*SINANG * 1.0*FFL DVECY(7) = DVECY(7) * ( 1 - UHWKS*ACSPEEDI ) C c JOE - didn't have the prefactor needed here ... doh! !!! DVECY(9) = RTE2*(RHOK-RHOC) c Tom - The mass equation needs to be modified with a source term that c accounts for adiabatic compression. w_avg = 0.5*( dvecy(2) + dvecy(6) ) DVECY(9) = GAMAXAV/(GAMZ*BSV)*RTE2*(RHOK-RHOC) * 1.0 !!!!! DVECY(9) = GAMAV /(GAMZ*BSV)*RTE2*(RHOK-RHOC) * 1.0 !!!!! DVECY(9) = GAMZ /(GAMZ*BSV)*RTE2*(RHOK-RHOC) * 1.0 DVECY(9) = DVECY(9) * ( 1 - UHWKA*ACSPEEDI ) & - rhoc*grav/(gamma_gas*R_gas*tdegk)*w_avg c JOE - added particle density variable and have it decay due to c detraiment using RTE3 DVECY(10)= GAMAXAV/(GAMZ*BSV)*RTE3*( -VECY(10) ) !!!!! DVECY(10)= GAMAV /(GAMZ*BSV)*RTE3*( -VECY(10) ) !!!!! DVECY(10)= GAMZ /(GAMZ*BSV)*RTE3*( -VECY(10) ) DVECY(10)= DVECY(10) * ( 1 - UHWKA*ACSPEEDI ) C GAMAVB = GAMAV ARGTB = ARGT SEPOLD = SEPTOT DBDT1 = DBDT BDBDT1 = BDBDT RETURN END C SUBROUTINE RK42(Y,DYDX,N,YOUT,TL,TSS,time) c JOE - changed NVAR to 10 so i can include particle density variable PARAMETER (NPTMAX=160, NVAR=2*NPTMAX) DIMENSION Y(*),DYDX(*),YOUT(*),YT(NVAR),DYT(NVAR),DYM(NVAR) c JOE - what time is it? in variable "time". REAL time C COMMON /HB/H,HH,H6,XH C DO 11 I=1,N YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE c JOE - added TSS as last arg to DERIVS2() call !!!!! CALL DERIVS2(HH,YT,DYT,TSS,time+HH) CALL DERIVS2(HH,YT,DYT,TL,TSS,time+HH) DO 12 I=1,N YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE c JOE - added TSS as last arg to DERIVS2() call !!!!! CALL DERIVS2(HH,YT,DYM,TSS,time+HH) CALL DERIVS2(HH,YT,DYM,TL,TSS,time+HH) DO 13 I=1,N YT(I)=Y(I)+H*DYM(I) DYM(I)=DYT(I)+DYM(I) 13 CONTINUE !!!!! CALL DERIVS2(H,YT,DYT,TSS,time+H) CALL DERIVS2(H,YT,DYT,TL,TSS,time+H) DO 14 I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE RETURN END C C C * * * * * * * * * * * * * * * QSF * * * * * * * * * * * * * * * * * C C SUBROUTINE QSF(H,YIN,ZOUT,NDIM) C C H -THE INCREMENT OF ARGUMENT VALUES. C YIN -THE INPUT VECTOR OF FUNCTION VALUES. C ZOUT -THE RESULTING VECTOR OF INTEGRAL VALUES. Z MAY BE C IDENTICAL WITH Y. C NDIM -THE DIMENSION OF VECTORS YIN AND ZOUT. C REAL H,YIN(*),ZOUT(*) INTEGER NDIM C PARAMETER (NG=513) C CC DOUBLE PRECISION Y(NG),Z(NG),HT,SUM1,AUX1,AUX2,SUM2,AUX,YIM2 REAL Y(NG),Z(NG) C IF(MOD(NDIM,2).EQ.0) THEN WRITE(6,91) NDIM 91 FORMAT(1X,'NDIM SHOULD BE ODD. NDIM = ',I5) STOP ENDIF C HT = (.3333333)*(H) C DO 100 IG=1,NDIM 100 Y(IG) = (YIN(IG)) C IF(NDIM-5)7,8,1 C C NDIM IS GREATER THAN 5, PREPARATIONS OF INTEGRATION LOOP 1 CONTINUE SUM1 = Y(2) + Y(2) SUM1 = SUM1 + SUM1 SUM1=HT*(Y(1)+SUM1+Y(3)) AUX1 = Y(4) + Y(4) AUX1 = AUX1 + AUX1 AUX1=SUM1+HT*(Y(3)+AUX1+Y(5)) AUX2=HT*(Y(1)+(3.875)*(Y(2)+Y(5)) X +(2.625)*(Y(3)+Y(4))+Y(6)) SUM2 = Y(5) + Y(5) SUM2 = SUM2 + SUM2 SUM2=AUX2-HT*(Y(4)+SUM2+Y(6)) Z(1) = (0.) AUX = Y(3) + Y(3) AUX = AUX + AUX Z(2)=SUM2-HT*(Y(2)+AUX+Y(4)) Z(3)=SUM1 Z(4)=SUM2 IF(NDIM-6)5,5,2 C C INTEGRATION LOOP 2 DO 4 I=7,NDIM,2 YIM2=Y(I-2) SUM1=AUX1 Z(I-2)=SUM1 SUM2=AUX2 AUX1 = Y(I-1) + Y(I-1) AUX1 = AUX1 + AUX1 AUX1=SUM1+HT*(YIM2+AUX1+Y(I)) IF(I-NDIM)3,6,6 3 CONTINUE AUX2 = Y(I) + Y(I) AUX2 = AUX2 + AUX2 AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1)) 4 Z(I-1)=SUM2 5 Z(NDIM-1)=AUX1 Z(NDIM)=AUX2 GO TO 12 6 Z(NDIM-1)=SUM2 Z(NDIM)=AUX1 GO TO 12 C END OF INTEGRATION LOOP C 7 IF(NDIM-3)12,11,8 C C NDIM IS EQUAL TO 4 OR 5 8 SUM2=(1.125)*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4)) SUM1 = Y(2) + Y(2) SUM1 = SUM1 + SUM1 SUM1=HT*(Y(1)+SUM1+Y(3)) Z(1)= (0.) AUX1 = Y(3) + Y(3) AUX1 = AUX1 + AUX1 Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4)) IF(NDIM-5)10,9,9 9 CONTINUE AUX1 = Y(4) + Y(4) AUX1 = AUX1 + AUX1 Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5)) 10 Z(3)=SUM1 Z(4)=SUM2 GO TO 12 C C NDIM IS EQUAL TO 3 11 SUM1=HT*((1.25)*Y(1)+Y(2)+Y(2)-(.25)*Y(3)) SUM2 = Y(2) + Y(2) SUM2 = SUM2 + SUM2 Z(3)=HT*(Y(1)+SUM2+Y(3)) Z(1)= (0.) Z(2)=SUM1 C 12 CONTINUE DO 200 IG=1,NDIM 200 ZOUT(IG) = (Z(IG)) RETURN END