C --- ------------------------------------------------------------------ C --- HHMS:TREC.07d: Read adjustable parameters from file Params.trec C Vmin, Vmax, Vxp, Efx, Ptotal C ------------------ C Vr = Vmin + (Vmax-Vmin)/(1. + (Ef/Efx)**Vxp) C Br = 10.*Bss*Bscale C Vp = Omega*1.5*Rs*sin(Th) C Bp = Br/Vr*(Vp - Omega*Rad(1)*sin(Th)) C Ro = FluxMass/Vr C Pfield = B^2/Pi8 C Tp = (Ptotal - Pfield)/(2*Ro*Rgas) C C --- INPUTS from t3din: (new structure of t3din) C GridFN C Fm0,Bsc0,SSmap0 C Fm1,Bsc1,SSmap1 C TITLE C Tshock, ShLat,ShLng,ShRad,ShSpd,ShTau (F10,3F4,2F5) C --> Shock Jumps are now computed internally C DTRUN, DTCHK1, DTCHK C AVC, CFL C C --- Owl (Offset West Longitude) = 30.0 C C --- Temporal-Roto-Extapo-Conversion (TREC) Model C Interpolate in Time between rotating Source Surface maps while C Extrapolating from Source Surface to 21.5 Rs and Converting C from source surface map parameters to MHD parameters. C --- ------------------------------------------------------------------ SUBROUTINE INITIALIZE( TCHECK,DTCHK,TSTOP,AVC,CFL) C --- Initialize the grid, source surface current sheet maps, C --- TREC empirical parameters, Shock input parameters, etc. PARAMETER ( H2S = 3600., S2H = 1./H2S ) PARAMETER ( D2R = 3.141593/180., KE = 12 ) INCLUDE 'igm_grid.ins' COMMON /BLOC15/ Tshock, ShTha,ShPhi,ShRad,ShSpd,ShTau, JCD,KCD COMMON /BLOC18/ RO0(KX,JX),VR0(KX,JX),VT0(KX,JX),VP0(KX,JX), & BR0(KX,JX),BT0(KX,JX),BP0(KX,JX),TP0(KX,JX) COMMON /BLOC21/ VRI(IX),VRU(IX),TPI(IX),TPU(IX) COMMON /BLOC22/ BPI(IX),BPU(IX),BTI(IX),BTU(IX),RC(IX) CHARACTER*80 TITLE, GridFN, SSmap0,SSmap1 REAL*8 TCHECK,DTCHK,TSTOP,Tshock REAL*8 Tmap0, Tmap1 COMMON /MAPS/ Bbc0(KX,JX),Vbc0(KX,JX), Bbc1(KX,JX),Vbc1(KX,JX) & , Ncr0,Long0,Tmap0,Fm0,Bsc0, Ncr1,Long1,Tmap1,Fm1,Bsc1 COMMON /GSTUFF/ TITLE,IRUN,ISEG,DRX,DTX,DPX COMMON /TRECGEAR/ Vmin,Vmax,Vxp,Efx, Ptotal, Owl C *** READ-IN INPUT DATA, AND ECHO IT OPEN(5,FILE='t3din',STATUS='OLD') READ (5,'(A)') GridFN READ (5,'(E7.2,F6.3,X,A)') Fm0,Bsc0,SSmap0 READ (5,'(E7.2,F6.3,X,A)') Fm1,Bsc1,SSmap1 READ (5,'(A)') TITLE READ (5,*) Tshock, ShLat,ShLng,ShRad,ShSpd,ShTau READ (5,*) DTRUN, DTCHK1, DTCHK READ (5,*) AVC, CFL CLOSE(5) WRITE (6,*) ' ' WRITE (6,*) 'INPUTS from t3din: ' WRITE (6,'(A,A)') ' INITIAL GRID fn: ', GridFN WRITE (6,'(A,E10.3,F7.3,X,A)') ' SSmap0: ', Fm0,Bsc0, SSmap0 WRITE (6,'(A,E10.3,F7.3,X,A)') ' SSmap1: ', Fm1,Bsc1, SSmap1 WRITE (6,'(A,A)') ' TITLE: ', TITLE WRITE (6,'(F12.0,5F8.2)') Tshock, ShLat,ShLng,ShRad,ShSpd,ShTau WRITE (6,*) DTRUN, DTCHK1, DTCHK WRITE (6,*) AVC, CFL WRITE (6,*) ' ' C --- Initialize grid by reading in a grid file CALL INITGRID(GridFN ) C --- Initial set up for subroutine SHOCKER C --- Convert disk position of shock center to grid position ShTha = (90. - ShLat)*D2R ShPhi = Phi(KE) + ShLng*D2R ShRad = ShRad*D2R ! convert shock input radius to radians ShTau = ShTau*3600. ! convert shock duration to seconds PRINT *, ' ' PRINT *, 'Initializing shocker' PRINT *, ' Shock Lat:', ShLat, ' ->', ShTha/D2R PRINT *, ' Shock Lng:', ShLng, ' ->', ShPhi/D2R C --- Find nearest grid point to shock center: -> KCD,JCD C --- Find Jlo,Jhi such that TH(Jhi) .GT. ShTha .GE. TH(Jlo) Jlo = 1 Jhi = JM 1 IF ( Jhi-Jlo .gt. 1 ) THEN J = (Jhi+Jlo)/2 IF ( Th(J) .gt. ShTha ) THEN Jhi = J ELSE Jlo = J ENDIF GO TO 1 ENDIF IF ( Th(Jhi) - ShTha .lt. 0.5*dTh ) THEN JCD = Jhi ELSE JCD = Jlo ENDIF C --- Find Khi,Klo such that Phi(Khi) .GT. ShPhi .GE. PHI(Klo) Klo = 1 Khi = KM 2 IF ( Khi-Klo .gt. 1 ) THEN K = (Khi+Klo)/2 IF ( Phi(K) .gt. ShPhi ) THEN Khi = K ELSE Klo = K ENDIF GO TO 2 ENDIF IF ( Phi(Khi) - ShPhi .lt. 0.5*dPh ) THEN KCD = Khi ELSE KCD = Klo ENDIF IF ( Phi(Khi) - ShPhi .lt. 0.5*dPh ) THEN KCD = Khi ELSE KCD = Klo ENDIF PRINT *, ' Set JCD,KCD:', JCD,KCD C --- TITLE & INITIATE 1-D MOVIES W/ STATE ALONG RADIAL OPEN(8,FILE='bqmovie' &,FORM='UNFORMATTED',STATUS='UNKNOWN',ACCESS='APPEND') WRITE (8) TITLE(1:80) WRITE (8) Time, IM K = KCD J = JCD DO 8 I=1,IM 8 WRITE (8) RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I), & BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I) C --- Initalize S.Surf Model based Lower Boundary Condition data C --- Owl is "Offset West Longitude", when a parcel reaches the inflow C grid boundary its point of orgin is 'Owl' degrees westward C Owl is used in two places to convert between Carrington & Time Owl = 30.0 OPEN(5,FILE='Params.trec',STATUS='OLD') WRITE (6,*) ' ' WRITE (6,*) 'INPUTS from Params.trec: ' READ (5,*) Vmin,Vmax,Vxp,Efx,Ptotal WRITE(6,*) Vmin,Vmax,Vxp,Efx,Ptotal CLOSE(5) CALL GETSSM(SSmap0,Bsc0, Ncr0,Long0,Tmap0,Bbc0,Vbc0) WRITE (6,101) ' ---> SSmap0:', Ncr0,Long0, Tmap0, Bsc0 CALL GETSSM(SSmap1,Bsc1, Ncr1,Long1,Tmap1,Bbc1,Vbc1) WRITE (6,101) ' ---> SSmap1:', Ncr1,Long1, Tmap1, Bsc1 WRITE (6,*) ' ' IF ( Time .LT. Tmap0 .OR. Time .GT. Tmap1 ) THEN PRINT *, 'Grid Time: ', Time PRINT *, ' outside interval of maps:' PRINT *, Tmap0, Tmap1 STOP 'Grid Time out of bounds.' ENDIF TCHECK = Time + DTCHK1*H2S DTCHK = DTCHK*H2S TSTOP = Time + DTRUN*H2S TSTOP = MIN(TSTOP,Tmap1) RETURN 120 STOP 'EOF IN INIT' 101 FORMAT (A,I5,I4,F12.1,F7.3) END C --- ------------------------------------------------------------------ SUBROUTINE INITGRID(GridFN ) C --- INITIALIZE THE GRID, THE ZERO SHELL (STEADY STATE) PARAMETER ( PI1 = 3.141593, PI4 = 4.*PI1, PI8 = 8.*PI1 ) INCLUDE 'igm_grid.ins' COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 COMMON /BLOC15/ Tshock, ShTha,ShPhi,ShRad,ShSpd,ShTau, JCD,KCD COMMON /BLOC18/ RO0(KX,JX),VR0(KX,JX),VT0(KX,JX),VP0(KX,JX), & BR0(KX,JX),BT0(KX,JX),BP0(KX,JX),TP0(KX,JX) COMMON /GSTUFF/ TITLE,IRUN,ISEG,DRX,DTX,DPX CHARACTER*80 TITLE, TITLEX, GridFN REAL*8 Tshock DATA PI/PI1/, XMU0/PI4/, XMU1/PI8/ DATA GAM/1.46/, R/8.32E-3/, GMS/1.3256E11/ ,RS/6.95E+5/ WRITE (6,'(/,A)') ' INITGRID: ' C --- INITALIZE GRID & ZERO SHELL WITH SOLUTION FROM FILE 'GridFN' OPEN(1,FILE=GridFN,STATUS='OLD',FORM='UNFORMATTED') READ (1) TITLEX(1:80) READ (1) ITIME, IrunX, IsegX, IM, JM, KM, DRX, DTX, DPX IF ((IM.GT.IX).OR.(JM.GT.JX).OR.(KM.GT.KX)) THEN PRINT *, 'WARNING: ARRAYS TOO BIG' STOP 'CHECK igm_grid.ins' ENDIF WRITE (6,'(X,A)') TITLEX WRITE (6,*) ITIME, IrunX, IsegX, IM, JM, KM, DRX, DTX, DPX Time = DBLE(ITIME) C --- SET UP GEOMETRICAL GRID AND PHYSICAL CONSTANTS DR = DRX*RS DTH = DTX*PI/180. DPH = DPX*PI/180. CALL GRID GAM1 = GAM - 1.0 GMM1 = GAM/GAM1 DMU0 = 1./XMU0 DMU1 = 1./XMU1 C --- READ IN ZERO SHELL (LOWER BOUNDARY ORIGINAL STATE) DO 4 J=1,JM READ (1) ( RO0(K,J),VR0(K,J),VT0(K,J),VP0(K,J), & BR0(K,J),BT0(K,J),BP0(K,J),TP0(K,J),K=1,KM ) 4 CONTINUE C --- READ IN FULL 3-D GRID DO 5 I=1,IM DO 5 J=1,JM READ (1) ( RO(K,J,I),VR(K,J,I),VT(K,J,I),VP(K,J,I), & BR(K,J,I),BT(K,J,I),BP(K,J,I),TP(K,J,I),K=1,KM ) 5 CONTINUE CLOSE(1) C --- SET COMPUTATIONAL DOMAIN JM1 = JM-1 KM1 = KM-1 IC1 = 2 IC2 = IM-1 JC1 = 2 JC2 = JM1 KC1 = 2 KC2 = KM1 RETURN 120 STOP 'EOF IN INIT' END C --- ------------------------------------------------------------------ SUBROUTINE GRID C --- GENERATE THE COMPUTATIONAL GRID MESH SYSTEM. INCLUDE 'igm_grid.ins' COMMON /BLOC7/ GMS,RS,GAM,GAM1,GMM1,R,PI,XMU0,XMU1,DMU0,DMU1 TH(1) = 17.5*PI/180. PHI(1) = 2.5*PI/180. RAD(1) = 21.5*RS DO 10 I=2,IM RAD(I)=RAD(I-1)+DR 10 CONTINUE DO 20 J=2,JM TH(J)=TH(J-1)+DTH 20 CONTINUE DO 30 K=2,KM PHI(K)=PHI(K-1)+DPH 30 CONTINUE RETURN END C --- ------------------------------------------------------------------ SUBROUTINE GETSSM (MapFile,Bscale, Ncr,Long,Tyme,Bbc,Vbc) C --- Read in and store Source Surface Map file for SSMLBC INCLUDE 'igm_grid.ins' COMMON /TRECGEAR/ Vmin,Vmax,Vxp,Efx, Ptotal, Owl CHARACTER*80 MapFile REAL*8 Tyme DIMENSION Bbc(KX,JX), Vbc(KX,JX) DIMENSION Blbc(KX,JX), Vlbc(KX,JX) DIMENSION bssp(72,37), fexp(72,37) Cx DIMENSION theta(72,37), fii(72,37), bphot(72,36) WRITE (6,*) ' ' WRITE (6,'(A)') ' GETSSM:' WRITE (6,'(A,A)') ' mapfile: ', MapFile OPEN (4,file=MapFile) READ (4,100) Ncr, Long READ (4,105) ((bssp(i,j),i=1,72),j=1,37) READ (4,103) ((fexp(i,j),i=1,72),j=1,37) Cx READ (4,104) ((theta(i,j),i=1,72),j=1,37) Cx READ (4,104) ((fii(i,j),i=1,72),j=1,37) Cx READ (4,105) ((bphot(i,j),i=1,72),j=1,36) CLOSE (4) C --- Translate to lbc conditions DO 10 J=1,JM Js = 34 - J DO 10 K=1,KM Blbc(K,J) = 10.*bssp(K,Js)*Bscale Vlbc(K,J) = Vmin + (Vmax-Vmin)/(1. + (fexp(K,Js)/Efx)**Vxp) 10 CONTINUE WRITE (6,*) ' translation done.' C --- Convert Map Coord's Ncr and Long into Posix time at grid boundary C --- CarrF is Carrington Nbr with fractional part C --- Owl is "Offset West Long." point of orgin is Owl deg. W of CM CarrF = Ncr - Long/360. Tyme = ((CarrF - Owl/360.)*27.2753D0 - 42420.5D0)*86400.D0 WRITE (6,*) ' Tyme = ', Tyme C --- Shift map cut to Long = 0 K0 = (360 - Long)/5 fp = ( MOD( 360.-REAL(Long), 5.) )/5. fm = 1. - fp WRITE (6,*) ' K0, Km, Jm:', K0, Km, Jm WRITE (6,*) ' fm, fp:', fm, fp DO 20 J=1,JM DO 20 K=1,KM K1 = mod(K-1+K0,KX) + 1 K2 = mod(K+K0,KX) + 1 Bbc(K,J) = fm*Blbc(K1,J) + fp*Blbc(K2,J) Vbc(K,J) = fm*Vlbc(K1,J) + fp*Vlbc(K2,J) 20 CONTINUE 100 FORMAT (2i5) 103 FORMAT (F10.3,11F11.3) 104 FORMAT (F10.4,11F11.4) 105 FORMAT (F10.5,11F11.5) END C --- ------------------------------------------------------------------ SUBROUTINE SSMLBC C --- Source Surface Map Lower Boundary Condition C --- Interpolate in time between maps, then rotate INCLUDE 'igm_grid.ins' COMMON/BLOC18/Ro0(KX,JX),Vr0(KX,JX),Vt0(KX,JX),Vp0(KX,JX), & Br0(KX,JX),Bt0(KX,JX),Bp0(KX,JX),Tp0(KX,JX) REAL*8 Tmap0, Tmap1, CarrF COMMON /MAPS/ Bbc0(KX,JX),Vbc0(KX,JX), Bbc1(KX,JX),Vbc1(KX,JX) & , Ncr0,Long0,Tmap0,Fm0,Bsc0, Ncr1,Long1,Tmap1,Fm1,Bsc1 COMMON /TRECGEAR/ Vmin,Vmax,Vxp,Efx, Ptotal, Owl REAL*4 Bbct(KX,JX), Vbct(KX,JX) PARAMETER ( Pi = 3.141593 ) PARAMETER ( Omega = 2.*Pi/(27.2753*86400.) ) PARAMETER ( Pi8 = 8.*Pi, Rgas = 8.32E-3 ) DATA Rs / 6.95E5/ C --- Interpolate in time between maps f1 = (Time-Tmap0)/(Tmap1-Tmap0) f0 = 1.0 - f1 DO 10 J=1,JM DO 10 K=1,KM Vbct(K,J) = f0*Vbc0(K,J) + f1*Vbc1(K,J) Bbct(K,J) = f0*Bbc0(K,J) + f1*Bbc1(K,J) FluxMass = f0*Fm0 + f1*Fm1 10 CONTINUE C --- Convert grid Posix time into Map Carrington Number and Longitude C --- CarrF is Carrington Nbr with fractional part C --- Owl is "Offset West Long." point of orgin is Owl deg. W of CM CarrF = (Time/86400.D0 + 42420.5D0)/27.2753D0 + Owl/360. rLong = ( 1. - (CarrF - AINT(CarrF)))*360. K0 = rLong/5 fp = mod(rLong,5.)/5. fm = 1. - fp C --- Rotate map and interpolate within map to grid points on lower boundary DO 20 J=1,JM DO 20 K=1,KM K1 = mod(K-1+K0,KX) + 1 K2 = mod(K+K0,KX) + 1 Vr(K,J,1) = fm*Vbct(K1,J) + fp*Vbct(K2,J) Br(K,J,1) = fm*Bbct(K1,J) + fp*Bbct(K2,J) Vt(K,J,1) = 0.0 Vp(K,J,1) = Omega*1.5*Rs*sin(TH(J)) Bt(K,J,1) = Br(K,J,1)*Vt(K,J,1)/Vr(K,J,1) Bp(K,J,1) = Br(K,J,1)/Vr(K,J,1) & *(Vp(K,J,1) - Omega*RAD(1)*sin(TH(J))) Ro(K,J,1) = FluxMass/Vr(K,J,1) Pfield = ( Br(K,J,1)**2 + Bp(K,J,1)**2 + Bt(K,J,1)**2 )/Pi8 Tp(K,J,1) = (Ptotal - Pfield)/(2.*Ro(K,J,1)*Rgas) 20 CONTINUE 101 FORMAT (A,I5,I4,F12.1) END C --- ----------------------------------------------------------------- SUBROUTINE BoAng(Time, Bang) C --- Compute Solar "B-angle" from POSIX Time PARAMETER ( R2D = 180.D0/3.1415926 ) REAL*8 Time, Hulien, eNN, eLL, GG, Sun, Omega, BoSin, BoCos REAL Bang Hulien = 2440587.5D0 + Time/86400.D0 eNN = Hulien - 2451545.0D0 eLL = 280.472D0 + 0.9856474D0*eNN GG = 357.528D0 + 0.9856003D0*eNN Sun = eLL + 1.915*SIN(GG/R2D) + 0.020*SIN(2.*GG/R2D) Omega = 73. + 4./6. + 50.25*(Hulien - 2396756.75D0)/365.25/3600. BoSin = SIN((Sun - Omega)/R2D)*SIN(7.25/R2D) BoCos = SQRT( 1.D0 - BoSin**2 ) Bang = ATAN2(BoSin,BoCos)*R2D END