C+ C NAME: C ENLIL_Read.f C PURPOSE: C Read the external volumetric elements provided by the ENLIL 3D-MHD program. C These volumetric elements are read from files provided by ENLIL that are to be translated into C the heliographic coordinate system needed by the UCSD tomography program. To operate through the end of the year, C the program assumes that the starting hour of the tomography is at 3 UT. C The cadence is specified by nT and the values of XCEA. (BVJ 12/27/2018). C CATEGORY: C I/O C CALLING SEQUENCE: C call ENLIL_Read(cWild3DENLIL,INN,nLng,nLat,iMap,nT,nTmax,LLfstt,LLEndd,iYr,XCintDG,MJDRefB, C RRS,dRR,RADS,DDD1,TTT1,VVV3,BBB3,iNUM) C C INPUTS: C cWild3DENLIL character*80 input ENLIL file name C INN integer the location in the file name where ENLIL places the value of the year. C nLng integer Number of Longitude bins C nLat integer Number of Latitude bins C iMap integer Number of Map heights C nT integer Number of times read C nTmax integer Maximum number of times possible C LLEndd integer Ending number of 6-hour averages C iYr integer The year of the beginning input file C NinterD integer Number of intermediate spaces between tomography intervals C XCintDG(nTmax) real*8 Day of Year C MJDRefB real*8 The beginning MJD C RRS real Beginning Height of the IPS modeling C dRR real Outward step height of the modeling C C OUTPUTS: C RADS(iMap+1) real Heights input from the model (The first height is RADMS). C DDD1(nLng,nLat,iMap+1,nT) real Input 1 component density C TTT1(nLng,nLat,iMap+1,nT) real Input 1 component temperature C VVV3(nLng,nLat,iMap+1,nT,3) real Input 3 component velocity (r,t,n) C BBB3(nLng,nLat,iMap+1,nT,3) real Input 3 component magnetic field (r,t,n) C iNUM integer # of files read C C MODIFICATION HISTORY: C November-2015, Bernard Jackson (UCSD) C- subroutine ENLIL_Read(cWild3DENLIL,INN,nLng,nLat,iMap,nT,nTmax,LLfstt,LLEndd,iYr,XCintDG,MJDRefB, & RRS,dRR,RADS,DDD1,TTT1,VVV3,BBB3,iNUM,NTBAD) real DDD1(nLng,nLat,iMap+1,nTmax), & TTT1(nLng,nLat,iMap+1,nTmax), & VVV3(nLng,nLat,iMap+1,nTmax,3), & BBB3(nLng,nLat,iMap+1,nTmax,3), & RADS(iMap+1) real*8 XCintDG(nT),MJDRefB,MJDRef,JEpoch real NTBAD(nT) character cdate*8 character cfile*80 character ccfile*80 character cWild*80 character cheader*132 character cWild3DENLIL*80 character cFmtF4*9 /'(55F10.4)'/ character cMon*3 include 'dirspec.h' bad = badr4() badD = -999.9999 badV = -999.9999 badB = -999.9999 badT = -999.9999 print *, ' ' write (*,'(A,I3)')'Into external read for ENLIL outputs' Ifirst = 0 NNtot = 0 NNtotE = 0 iNUM = 0 do II=1,nT NNtotE = 0 if(II.ge.LLfstt.and.II.le.LLEndd) then ! Maybe needed someday. MJDref = XCintDG(II) - XCintDG(1) + MJDrefB call Julian(11,iYrr,Doy,MJDref,JEpoch) iDoy = int(Doy) call DATE_DOY(1,iYrr,cMon,iMon,iDay,iDoy) C MJDref = XCintDG(II) - XCintDG(1) + MJDrefB C Doy = sngl(XCintDG(1)) C iDoy = nint(Doy) C iH = nint((Doy - iDoy)*24.0) C call Julian(11,iYrr,Doy,MJDref,JEpoch) ! Where do the ENLIL data files begin in year and day of year? C iDoy = nint(Doy) C call DATE_DOY(1,iYrr,cMon,iMon,iDay,iDoy) mo = iMon iD = iDay iH = nint((Doy - iDoy)*24.0) C print *, 'ENLIL read 1', iYrr, DoY, iDoY, iH ! This is the time of Earth at the Carrington rotation value at Earth C print *, 'ENLIL read 2', iYrr, iMon, iDay, INN ! This is the year, month, and day at the beginning of the output sequence ! The first data begins at the next UT interval specified cfile = cWild3DENLIL write(cfile(INN:INN+9),'(I4,3I2.2)') iYrr,mo,iD,iH C write (*,'(A,A50)') 'cfile after additions =', cFile open (13, file=cFile,status='old',recl=1000,access='sequential',form='formatted',iostat=iRead3DMHD) C print*,'to here after file open attempt' C print *, ' ' if(iRead3DMHD.eq.0) write(*,'(A,I3,A,A40)') 'II =',II,' The file opened was ',cfile if(iRead3DMHD.ne.0) then close(13) write(*,'(A,I3,A,A40)') 'II =',II,' The file not found was ',cfile NTBAD(II) = 99 go to 1999 end if read (13,'(A)',iostat=iRead3DMHD) cheader ! read header - someday use parameters? if(iread3DMHD.ne.0) then close(13) print *, 'Something is wrong. The header lines could not be read.' go to 1999 end if Ifirst = Ifirst + 1 if(Ifirst.eq.1) write(*,'(A,I4)') 'Solar distances read; imap = ',iMap read (13,'(55F10.6)',iostat=iRead3DMHD) (RADS(N),N=1,iMAP) ! read solar distances if(iread3DMHD.ne.0) then close(13) print *, 'Something is wrong. The solar distances could not be read.' go to 1999 end if if(Ifirst.eq.1) write(*,'(31F10.6)') (RADS(N),N=1,iMAP) NBDDD1 = 0 NBTTT1 = 0 NBVVV1 = 0 NBVVV2 = 0 NBVVV3 = 0 NBBBB1 = 0 NBBBB2 = 0 NBBBB3 = 0 do IP=1,8 do N=1,iMAP FALLN = 1.0 FALLT = 1.0 FALLBR = 1.0 FALLBT = 1.0 FALLBN = 1.0 do JJ=1,nLat J=(nLat+1)-JJ ! BVJ 2019/03/20 ENLIL outputs begin at the UPPER left corner, unlike the tomography C print*, 'before all parameter read', IP,J,N, II if(IP.eq.1) read (13,cFmtF4,iostat=iRead3DMHD) (DDD1(I,J,N,II), I=1,nLng) ! read density If(IP.eq.2) read (13,cFmtF4,iostat=iRead3DMHD) (VVV3(I,J,N,II,1),I=1,nLng) ! read Vr C if(IP.eq.2.and.J.eq.5.and.N.eq. 1.and.II.eq.22) write (*,cFmtF4) (VVV3(I,J,N,II,1), I=1,nLng) C if(IP.eq.2.and.J.eq.5.and.N.eq.iMAP.and.II.eq.22) write (*,cFmtF4) (VVV3(I,J,N,II,1), I=1,nLng) If(IP.eq.3) read (13,cFmtF4,iostat=iRead3DMHD) (VVV3(I,J,N,II,2),I=1,nLng) ! read Vt If(IP.eq.4) read (13,cFmtF4,iostat=iRead3DMHD) (VVV3(I,J,N,II,3),I=1,nLng) ! read Vn If(IP.eq.5) read (13,cFmtF4,iostat=iRead3DMHD) (BBB3(I,J,N,II,1),I=1,nLng) ! read Br If(IP.eq.6) read (13,cFmtF4,iostat=iRead3DMHD) (BBB3(I,J,N,II,2),I=1,nLng) ! read Bt If(IP.eq.7) read (13,cFmtF4,iostat=iRead3DMHD) (BBB3(I,J,N,II,3),I=1,nLng) ! read Bn If(IP.eq.8) read (13,cFmtF4,iostat=iRead3DMHD) (TTT1(I,J,N,II), I=1,nLng) ! read temperature C print*, 'after all parameter read', IP,J,N, II do III=1, nlng if(IP.eq.1) then if(DDD1(III,J,N,II).gt.BadD) then DDD1(III,J,N,II) = DDD1(III,J,N,II)/FALLN C if(III.eq.20.and.J.eq.6.and.II.eq.24) write(*,'(A,4I4,F10.1,2F10.5)'), C & 'III, J, N, II, DDD1',III,J,N,II,DDD1(III,J,N,II),FALLN,RADS(N) else DDD1(III,J,N,II) = Bad NBDDD1 = NBDDD1 + 1 end if end if if(IP.eq.2) then if(VVV3(III,J,N,II,1).le.BadV) then VVV3(III,J,N,II,1) = Bad NBVVV1 = NBVVV1 + 1 end if end if if(IP.eq.3) then if(VVV3(III,J,N,II,2).le.BadV) then VVV3(III,J,N,II,2) = Bad NBVVV2 = NBVVV2 + 1 end if end if if(IP.eq.4) then if(VVV3(III,J,N,II,3).le.BadV) then VVV3(III,J,N,II,3) = Bad NBVVV3 = NBVVV3 + 1 end if end if if(IP.eq.5) then if(BBB3(III,J,N,II,1).gt.BadB) then BBB3(III,J,N,II,1) = BBB3(III,J,N,II,1)/FALLBR else BBB3(III,J,N,II,1) = Bad NBBBB1 = NBBBB1 + 1 end if end if if(IP.eq.6) then C print*, 'First 6 parameter read', IP,J,N, II, III, nLng, nLat, iMapm1, NF, IINC if(BBB3(III,J,N,II,2).gt.BadB) then BBB3(III,J,N,II,2) = BBB3(III,J,N,II,2)/FALLBT else BBB3(III,J,N,II,2) = Bad NBBBB2 = NBBBB2 + 1 C print*, 'Second 7 parameter read', IP,J,N, II, III, nLng, nLat, iMapm1, NF, IINC end if end if If(IP.eq.7) then if(BBB3(III,J,N,II,3).gt.BadB) then BBB3(III,J,N,II,3) = BBB3(III,J,N,II,3)/FALLBN else BBB3(III,J,N,II,3) = Bad NBBBB3 = NBBBB3 + 1 end if end if if(IP.eq.8) then if(TTT1(III,J,N,II).gt.BadT) then TTT1(III,J,N,II) = TTT1(III,J,N,II)/FALLT else TTT1(III,J,N,II) = Bad NBTTT1 = NBTTT1 + 1 end if end if end do end do end do end do iNUM = iNUM + 1 close(13) NNtot = NBDDD1+NBVVV1+NBVVV2+NBVVV3+NBBBB1+NBBBB2+NBBBB3+NBTTT1 C if(NNtot .eq. NNtotE) write(*,'(A,I7)') 'NNtot = ',NNtot if(NNtot .ne. 0) then if(NBDDD1.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBDDD1,' bad points in this density read' if(NBVVV1.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBVVV1,' bad points in this radial velocity read' if(NBVVV2.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBVVV2,' bad points in this tangential velocity read' if(NBVVV3.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBVVV3,' bad points in this normal velocity read' if(NBBBB1.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBBBB1,' bad points in this radial field read' if(NBBBB2.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBBBB2,' bad points in this tangential field read' if(NBBBB3.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBBBB3,' bad points in this normal field read' if(NBTTT1.ne.0) write(*,'(A,I3,A,I6,A)') 'II =',II,' There were',NBTTT1,' bad points in this temperature read' print *, ' ' end if C NNtotE = NNtot 1999 continue end if ! Maybe remove someday end do print *, ' ' if(iNUM.eq.0) write(*,'(A)') 'The external read for ENLIL input files was not successful.' if(iNUM.gt.0) write(*,'(A,I3,A,I3,A)') 'The external read for ENLIL files was successful.', & iNUM,' of',nT, ' files were read.' print *, ' ' return end