subroutine read_inputs use ModAMIE use ModFiles use ModData use ModConductance use ModIndices implicit none include "mpif90.h" integer :: ierr, i, n, m, noff logical :: done, exist logical :: starttime_read = .false. logical :: timelist_read = .false. logical :: dt_read = .false. logical :: ntimes_read = .false. logical :: mhd_solar_wind = .false. logical :: old_indices = .false. logical :: ngdc_indices = .false. logical :: noaa_indices = .false. character (len=100) :: line integer, parameter :: unitin=20 integer, parameter :: maxline=500 character (len=100), dimension(maxline) :: text ! Total number of lines, and current line number integer :: iline, nline character (LEN=*), parameter :: inputfile="AMIE.in" integer :: year, month, day, hour, minute, second integer, dimension(7) :: itime real*8 :: time integer :: iConductanceSolutions call set_defaults if(me_world==0)then nline=0 inquire(file=inputfile,EXIST=exist) if(.not.exist)call stop_mpi(& inputfile//" cannot be found by read_inputs") open(unitin,file=inputfile,status="old") do read(unitin,'(a)',END=100,ERR=100) line if(index(line,'#END')<1)then ! Store line into buffer nline=nline+1 if(nline>maxline)& call stop_mpi("Too many lines of input in read_inputs") text(nline)=line cycle else exit end if 100 continue close (unitin) end do if(nline==0)call stop_mpi(& "No lines of input read by read_inputs") end if ! Broadcast the number of lines and the text itself to all processors call MPI_Bcast(nline,1,MPI_Integer,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Number of lines could not be broadcast by read_mpi") call MPI_Bcast(text,len(text(1))*nline,MPI_Character,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Text could not be broadcast by read_mpi") done = .false. iline = 1 do while (.not.done) line = text(iline) if (line(1:1) == "#") then if (DebugLevel > 2) & write(*,*) " ===> Working on Command : ",line, me_world if (index(line,"#STARTTIME") > 0) then call read_in_int(year) call read_in_int(month) call read_in_int(day) call read_in_int(hour) call read_in_int(minute) call read_in_int(second) itime(1) = year itime(2) = month itime(3) = day itime(4) = hour itime(5) = minute itime(6) = second itime(7) = 0 call time_int_to_real(itime, starttime) starttime_read = .true. endif if (index(line,"#TIMELIST") > 0) then if (.not. ntimes_read) then write(*,*) "#NTIMES needs to be before #TIMELIST!!" call stop_mpi("Please reorganize your AMIE.in file.") endif do i=1,ntimes_total call read_in_time(time) timelist(i) = time call time_real_to_int(time, itime) enddo starttime_read = .true. timelist_read = .true. endif if (index(line,"#DT") > 0) then call read_in_int(dt) call read_in_int(dt_window) window_dt = dt_window dt_read = .true. endif if (index(line,"#NTIMES") > 0) then call read_in_int(ntimes_total) ntimes_read = .true. endif if (index(line,"#NORTH") > 0) then call read_in_logical(north) if (.not.north) ihsoln = -1 endif if (index(line,"#OUTPUT") > 0) then call read_in_string(output_file) call read_in_logical(IsOutputBinary) call read_in_logical(OutputData) call read_in_logical(OutputMagnetometers) endif if (index(line,"#MAGNETOMETER") > 0) then use_magnetometers = .true. call read_in_string(magnetometer_file) call read_in_string(psi_file) call read_in_logical(use_zed) call read_in_logical(use_ahn) call read_in_real(weight_on_x) call read_in_real(weight_on_y) call read_in_real(weight_on_z) endif if (index(line,"#EFIELD") > 0) then if (index(line,"WEIGHT") > 0) then do i= 1, nefield_files call read_in_11weights(efield_weightings(i,:)) enddo else call read_in_int(nefield_files) call read_in_string(efield_dir) do i= 1, nefield_files call read_in_flag(efield_flag(i)) call read_in_string(efield_file(i)) enddo endif endif if (index(line,"#SOLVECONDUCTANCE") > 0) then call read_in_logical(UseConductanceSolve) endif if (index(line,"#CONDUCTANCE") > 0) then call read_in_int(ncond_files) call read_in_string(cond_dir) do i= 1, ncond_files call read_in_string(cond_file(i)) enddo endif if (index(line,"#CURRENT") > 0) then call read_in_int(ncurrent_files) call read_in_string(current_dir) do i= 1, ncurrent_files call read_in_string(current_file(i)) ! call read_in_flag2(current_flag(i)) iline = iline + 1 read (text(iline),*) & current_flag(i),current_weightings(i,1:11) enddo endif if (index(line,"#DEBUG") > 0) then call read_in_int(DebugLevel) if(me_world==0 .and. DebugLevel > 0)write(*,*)& 'Read_inputs read and broadcast ',nline,' lines of text' endif if (index(line,"#BACKGROUND") > 0) then call read_in_string(models_dir) call read_in_string(efield_model) call read_in_string(cond_model) call read_in_string(solar_model) call merge_str(models_dir, apex_file) call merge_str(models_dir, qset_file) if (index(cond_model,'IHP') > 0) cond_model = 'ihp' if (index(cond_model,'PEM') > 0) cond_model = 'pem' if (index(efield_model,'weimer01') > 0) efield_model = 'weimer01' if (index(efield_model,'Weimer01') > 0) efield_model = 'weimer01' if (index(efield_model,'WEIMER01') > 0) efield_model = 'weimer01' if (index(efield_model,'weimer') > 0 .and. & index(efield_model,'01') < 0) efield_model = 'weimer96' if (index(efield_model,'Weimer') > 0 .and. & index(efield_model,'01') < 0) efield_model = 'weimer96' if (index(efield_model,'WEIMER') > 0 .and. & index(efield_model,'01') < 0) efield_model = 'weimer96' if (index(efield_model,'weimer96') > 0) efield_model = 'weimer96' if (index(efield_model,'Weimer96') > 0) efield_model = 'weimer96' if (index(efield_model,'WEIMER96') > 0) efield_model = 'weimer96' if (index(efield_model,'SAMIE') > 0) efield_model = 'samie' endif ! if (index(line,"#INDICES") > 0) then ! call read_in_string(indices_file) ! call read_in_real(time_delay) ! endif if (index(line,"#MHD_INDICES") > 0) then call read_in_int(nMHDFiles) do i=1,nMHDFiles call read_in_string(cMHDFile(i)) write(*,*) "MHD FILE : ", i, cMHDFile(i) enddo mhd_solar_wind = .true. endif if (index(line,"#NOAA_INDICES") > 0) then call read_in_string(cNOAAFile) noaa_indices = .true. endif if (index(line,"#OLD_INDICES") > 0) then call read_in_string(cOldFile) call read_in_real(time_delay) old_indices = .true. endif if (index(line,"#NGDC_INDICES") > 0) then call read_in_string(cNGDCFile) call read_in_real(time_delay) ngdc_indices = .true. endif if (index(line,"#FLUX_LIMITS") > 0) then do iConductanceSolutions = 1, nConductanceSolutions call read_in_real(ConductanceBackground(iConductanceSolutions)) enddo call read_in_real(bkgcer(1)) call read_in_real(bkgcer(2)) call read_in_real(bkgcer(3)) call read_in_real(bkgcer(4)) call read_in_real(flxmax) call read_in_real(dbymax) call read_in_real(dbzmax) endif endif iline = iline + 1 if (iline >= nline) done = .true. enddo if (DebugLevel > 2) & write(*,*) " ===> Done with reads : ",me_world if ((starttime_read .and. dt_read .and. ntimes_read).or. & (ntimes_read .and. timelist_read)) then m = ntimes_total / numprocs n = 0 if (m * numprocs /= ntimes_total) then n = ntimes_total - m * numprocs if (me_world < n) m = m + 1 endif ntimes = m currenttime = starttime noff = 1 do i = 0, me_world-1 m = ntimes_total / numprocs if (i < n) m = m + 1 currenttime = currenttime + m * dt noff = noff + m enddo if (timelist_read) then do i=noff,noff+ntimes_total-1 timelist(i-noff+1) = timelist(i) enddo currenttime = timelist(1) endif else call stop_mpi("Need to specify start time, dt, and ntimes or timelist") endif call time_real_to_int(currenttime, itime) if (DebugLevel > 1) & write(*,*) " ==> ntimes, me_world, currenttime ", ntimes, me_world, itime if (DebugLevel > 2) & write(*,*) " ===> Before read indices : ",me_world if (DebugLevel > 2) & write(*,*) " ===> old_indices : ",old_indices if (old_indices) call read_old_Indices if (DebugLevel > 2) & write(*,*) " ===> ngdc_indices : ",ngdc_indices if (ngdc_indices) call read_NGDC_Indices if (DebugLevel > 2) & write(*,*) " ===> ngdc_indices : ",noaa_indices if (noaa_indices) call read_NOAA_Indices if (DebugLevel > 2) & write(*,*) " ===> mhd_solar_wind : ",mhd_solar_wind if (mhd_solar_wind) call read_MHD_Indices where(Indices(:,hpi_) > 0) Indices(:,hpi_norm_) = 2.09 * ALOG(Indices(:,hpi_)) * 1.0475 endwhere where(Indices(:,hpi_norm_) > 10.0) Indices(:,hpi_norm_) = 10.0 endwhere where(Indices(:,hpi_norm_) < 1.0) Indices(:,hpi_norm_) = 1.0 endwhere !\ ! This is a hack! If kp is not read in, then use HPI data!! !/ where(Indices(:,kp_) < 0.0) Indices(:,kp_) = Indices(:,hpi_norm_) endwhere if (Indices(1,f107_) < 50.0) then write(*,*) "You must give at least on value of F10.7!!!" call stop_mpi("Please set this value in the AMIE.in file!") endif if (DebugLevel > 2) & write(*,*) " ===> Indices read : ",me_world contains subroutine read_in_int(variable) integer, intent(out) :: variable iline = iline + 1 read(text(iline),*) variable end subroutine read_in_int subroutine read_in_logical(variable) logical, intent(out) :: variable iline = iline + 1 read(text(iline),*) variable end subroutine read_in_logical subroutine read_in_string(variable) character (len=100), intent(out) :: variable iline = iline + 1 variable = text(iline) end subroutine read_in_string subroutine read_in_flag(variable) character (len=1), intent(out) :: variable iline = iline + 1 variable = text(iline) end subroutine read_in_flag subroutine read_in_flag2(variable) character (len=2), intent(out) :: variable iline = iline + 1 variable = text(iline) end subroutine read_in_flag2 subroutine read_in_11weights(variable) real, dimension(11) :: variable integer :: ierror iline = iline + 1 variable(1) = -2.0 read(text(iline), *, iostat = ierror) variable if (ierror /= 0) then if (variable(1) /= -1.0) then if (me_world == 0) then write(*,*) "It appears that you have an improper number of" write(*,*) "Weight lines. If you want the default weight values" write(*,*) "for the instrument, then you MUST have a -1.0 in the" write(*,*) "line. If you don't, I can't tell if you accidently" write(*,*) "skipped this instrument, or whether you wanted the" write(*,*) "default values..." endif call mpi_barrier(mpi_comm_world, ira) call stop_mpi("Error in read_in_11weights") else variable(htwindow_) = 900.0 variable(timeerror1_) = 1.0 variable(timeerror2_) = 1.0 variable(errorper1_) = 0.0 variable(errorper2_) = 0.0 variable(valueper1_) = 0.0 variable(valueper2_) = 0.0 variable(minerror1_) = 5.0 variable(minerror2_) = 5.0 variable(xwt1_) = 1.0 variable(xwt2_) = 1.0 endif endif end subroutine read_in_11weights subroutine read_in_real(variable) real :: variable iline = iline + 1 read(text(iline),*) variable end subroutine read_in_real subroutine read_in_time(variable) real*8 :: variable iline = iline + 1 read(text(iline),*) itime(1:6) itime(7) = 0 call time_int_to_real(itime, variable) end subroutine read_in_time subroutine set_defaults ConductanceBackground(pedersen_) = 0.2 ConductanceBackground(hall_) = 0.2 ConductanceBackground(eflux_) = 0.01 ConductanceBackground(avee_) = 1.0 bkgc(3) = 0.01 bkgc(4) = 3.0 bkgcer(1) = 0.2 bkgcer(2) = 0.2 bkgcer(3) = 1.0 bkgcer(4) = 5.0 UseConductanceSolve = .false. LinearDropOffConductance = .false. LinearDropOffTimeConductance = 600.0 LinearDropOffAveE = 0.5 LinearDropOffEFlux = 0.2 MaxeFoldingDeltaT = 1800.0 ExponentialDropOffeFolding = 600.0 LinearDropOffTimeEfield = 900.0 LinearDropOffEfield = 1.0 LinearDropOffTimeCurrent = 900.0 LinearDropOffCurrent = 1.0 flxmax = 85.0 dbymax = 50.0 dbzmax = 50.0 timelist = -1.0e32 Efield_Weightings(:,htwindow_) = 900.0 Efield_Weightings(:,timeerror1_) = 1.0 Efield_Weightings(:,timeerror2_) = 1.0 Efield_Weightings(:,errorper1_) = 0.0 Efield_Weightings(:,errorper2_) = 0.0 Efield_Weightings(:,valueper1_) = 0.0 Efield_Weightings(:,valueper2_) = 0.0 Efield_Weightings(:,minerror1_) = 5.0 Efield_Weightings(:,minerror2_) = 5.0 Efield_Weightings(:,xwt1_) = 1.0 Efield_Weightings(:,xwt2_) = 1.0 end subroutine set_defaults end subroutine read_inputs !============================================================================== subroutine read_MHD_Indices use ModAMIE use ModIndices use ModFiles use ModConstants implicit none include "mpif90.h" integer :: ierror, i, input_coor_system logical :: done ! One line of input character (len=100) :: line integer, dimension(1:7) :: Tmp_Time integer, dimension(Max_Upstream_Npts,7) :: Upstream_integer_time real*8 :: time_now real :: dt_min integer :: iyr, iday, ihour, imin, isec, idir, j, iFile real :: XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE, TimeDelayMHD integer, external :: jday Input_Coor_System = GSM if (me_world == 0) then Upstream_Npts = 0 do iFile = 1, nMHDFiles write(*,*) "MHD File : ", iFile, cMHDFile(iFile) open(lun_indices, file=cMHDFile(iFile), status="old", iostat = ierror) if (ierror /= 0) then call stop_mpi("Unable to open file "//cMHDFile//". Stopping.") endif done = .false. Input_Coor_System = GSM Propagation_Plane_XY = 0.0 Propagation_Plane_XZ = 0.0 write(6,*) "=> Reading Upstream Solar Wind Conditions File." do while (.not.done) read(lun_indices,'(a)', iostat = ierror ) line if (ierror.ne.0) done = .true. if(index(line,"#COOR")>0)then read(lun_indices,'(a)') line if(index(line,'GSE')>0) Input_Coor_System = GSE if(index(line,'GSM')>0) Input_Coor_System = GSM if(index(line,'GEO')>0) Input_Coor_System = GEO endif if(index(line,"#TIMEDELAY")>0)then read(lun_indices,*) TimeDelayMHD endif if(index(line,"#PLANE")>0)then read(lun_indices,*) Propagation_Plane_XY read(lun_indices,*) Propagation_Plane_XZ Propagation_Plane_XY = Propagation_Plane_XY * pi / 180.0 Propagation_Plane_XZ = Propagation_Plane_XZ * pi / 180.0 endif if(index(line,"#POSITION")>0)then read(lun_indices,*) Satellite_Y_Pos read(lun_indices,*) Satellite_Z_Pos endif if(index(line,"#START")>0)then do while (.not.done) Upstream_Npts = Upstream_Npts + 1 read(lun_indices,*,iostat=ierror) & (Tmp_Time(i),i=1,7), & (Upstream_Data(Upstream_Npts,i),i=1,8) if (ierror.ne.0.0) then done = .true. Upstream_Npts = Upstream_Npts - 1 else if (Upstream_Npts >= Max_Upstream_Npts) then done = .true. write(6,*) "=> This upstream boundary condition file" write(6,*) " contains too many lines! (Max lines =",& Max_Upstream_Npts else Upstream_integer_time(Upstream_Npts,:) = Tmp_Time(:) if (Upstream_integer_time(Upstream_Npts,1) >= 65 .and. & Upstream_integer_time(Upstream_Npts,1) < 100) & Upstream_integer_time(Upstream_Npts,1) = & Upstream_integer_time(Upstream_Npts,1) + 1900 if (Upstream_integer_time(Upstream_Npts,1) < 65) & Upstream_integer_time(Upstream_Npts,1) = & Upstream_integer_time(Upstream_Npts,1) + 2000 if (Input_Coor_System /= GSM) then if (Input_Coor_System == GSE) then iyr = Upstream_integer_time(Upstream_Npts,1) iday = & jday(Upstream_integer_time(Upstream_Npts,1),& Upstream_integer_time(Upstream_Npts,2), & Upstream_integer_time(Upstream_Npts,3)) ihour = Upstream_integer_time(Upstream_Npts,4) imin = Upstream_integer_time(Upstream_Npts,5) isec = Upstream_integer_time(Upstream_Npts,6) call RECALC(iyr,iday,ihour,imin,isec) idir = -1 xgse = Upstream_Data(Upstream_Npts,1) ygse = Upstream_Data(Upstream_Npts,2) zgse = Upstream_Data(Upstream_Npts,3) call GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,idir) Upstream_Data(Upstream_Npts,1) = xgsm Upstream_Data(Upstream_Npts,2) = ygsm Upstream_Data(Upstream_Npts,3) = zgsm xgse = Upstream_Data(Upstream_Npts,4) ygse = Upstream_Data(Upstream_Npts,5) zgse = Upstream_Data(Upstream_Npts,6) call GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,idir) Upstream_Data(Upstream_Npts,4) = xgsm Upstream_Data(Upstream_Npts,5) = ygsm Upstream_Data(Upstream_Npts,6) = zgsm else write(*,*) & "I don't know how to transform between the" write(*,*) & "input coordinate system and the code system!" stop endif endif endif endif enddo write(*,*) "pts read : ",Upstream_Npts,Upstream_integer_time(1,:) endif enddo close(lun_indices) enddo endif call MPI_Bcast(Upstream_Npts,1,MPI_Integer,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Upstream_Npts could not be broadcast by read_upstream_input_file") call MPI_Bcast(Upstream_integer_time,Max_Upstream_Npts*7,MPI_Integer, & 0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Upstream_integer_Time could not be broadcast by read_upstream_input_file") call MPI_Bcast(TimeDelayMHD,1,MPI_Real,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "TimeDelayMHD could not be broadcast by read_upstream_input_file") dt_min = 366.0*24.0*3600.0 write(*,*) "time : ",Upstream_integer_Time(1,:) do i=1,Upstream_Npts call time_int_to_real(Upstream_integer_Time(i,:), Upstream_Time(i)) Upstream_Time(i) = Upstream_Time(i) + TimeDelayMHD if (dt_min > abs(Upstream_Time(i) - currenttime)) then dt_min = abs(Upstream_Time(i) - currenttime) endif end do if ((dt_min > 24.0*3600.0).and.(me_world == 0)) then write(*,*) "******************************************************************************" write(*,*) "* *" write(*,*) "* Warning! Warning! Warning! Warning! Warning! Warning! Warning! Warning! *" write(*,*) "* *" write(*,*) "* Time dependent solar wind file disagrees with the starting time of the *" write(*,*) "* simulation by more than 24 hours. This could cause results which you *" write(*,*) "* may not enjoy. *" write(*,*) "* *" write(*,*) "* I am assuming that you are smarter than I am. Continuing with simulation. *" write(*,*) "* *" write(*,*) "******************************************************************************" endif call MPI_Bcast(Upstream_Data,Max_Upstream_Npts*8,MPI_Real, & 0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Upstream_Data could not be broadcast by read_upstream_input_file") call MPI_Bcast(Input_Coor_System,1,MPI_Integer,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Input_Coor_System could not be broadcast by read_upstream_input_file") call MPI_Bcast(Propagation_Plane_XY,1,MPI_Real,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Propagation_Plane_XY could not be broadcast by read_upstream_input_file") call MPI_Bcast(Propagation_Plane_XZ,1,MPI_Real,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Propagation_Plane_XZ could not be broadcast by read_upstream_input_file") call MPI_Bcast(Satellite_Y_Pos,1,MPI_Real,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Satellite_Y_Pos could not be broadcast by read_upstream_input_file") call MPI_Bcast(Satellite_Z_Pos,1,MPI_Real,0,MPI_COMM_WORLD,ira) if(ira>0)call stop_mpi(& "Satellite_Z_Pos could not be broadcast by read_upstream_input_file") if (Upstream_Npts == 1) then IMF_Bx = Upstream_Data(1,1) IMF_By = Upstream_Data(1,2) IMF_Bz = Upstream_Data(1,3) SW_v = -Upstream_Data(1,4) SW_n = Upstream_Data(1,7) else j = 1 do i=1,ntimes time_now = currenttime + (i-1) * dt do while ((j < Upstream_Npts).and. & (Upstream_Time(j) < time_now)) j = j + 1 enddo if ((j == Upstream_Npts .and. & Upstream_Time(j) <= time_now).or.(j==1)) then IMF_Bx(i) = Upstream_Data(j,1) IMF_By(i) = Upstream_Data(j,2) IMF_Bz(i) = Upstream_Data(j,3) SW_v(i) = -Upstream_Data(j,4) SW_n(i) = Upstream_Data(j,7) else dt_min = 1.0 - (Upstream_Time(j) - time_now) / & (Upstream_Time(j) - Upstream_Time(j-1) + 1.0e-6) IMF_Bx(i) = dt_min * Upstream_Data(j,1) + & (1.0 - dt_min) * Upstream_Data(j-1,1) IMF_By(i) = dt_min * Upstream_Data(j,2) + & (1.0 - dt_min) * Upstream_Data(j-1,2) IMF_Bz(i) = dt_min * Upstream_Data(j,3) + & (1.0 - dt_min) * Upstream_Data(j-1,3) SW_v(i) = dt_min * Upstream_Data(j,4) + & (1.0 - dt_min) * Upstream_Data(j-1,4) SW_n(i) = dt_min * Upstream_Data(j,7) + & (1.0 - dt_min) * Upstream_Data(j-1,7) SW_v(i) = -1.0 * SW_v(i) endif Indices(i,imf_bx_) = IMF_Bx(i) Indices(i,imf_by_) = IMF_By(i) Indices(i,imf_bz_) = IMF_Bz(i) Indices(i,sw_v_) = SW_v(i) j = j - 2 if (j < 1) j = 1 enddo endif end subroutine read_MHD_Indices !============================================================================== subroutine read_NGDC_Indices use ModAMIE use ModIndices use ModFiles implicit none integer :: ierror, year, month, i, j, npts, npts_hpi, k integer :: idir, input_coor_system, iday logical :: done, found ! One line of input character (len=100) :: line real, dimension(6,MAX_NTIMES) :: tmp real*8, dimension(MAX_NTIMES) :: ut_new, ut_keep, ut_tmp real, dimension(MAX_NTIMES) :: data_new, data_keep, data_tmp real :: xgse, ygse, zgse, xgsm, ygsm, zgsm real*8 :: time_now integer, dimension(7) :: itime integer, external :: jday open(lun_indices, file=cNGDCFile, status="old", iostat = ierror) if (ierror.ne.0) then call stop_mpi("Unable to open file "//cNGDCFile//". Stopping.") endif done = .false. npts_hpi = 0 Input_Coor_System = GSM do while (.not.done) read(lun_indices,'(a)', iostat = ierror ) line ! if (index(line, "#Element") > 0) then ! write(*,*) "New Line -->",line ! endif if (ierror.ne.0) done = .true. if(index(line,'#Element: dst')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, dst_) endif if(index(line,'#Element: kp')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, kp_) endif if(index(line,'#Element: bx')>0)then read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#Table: IMFMin')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, imf_bx_) endif endif if(index(line,'#Element: by')>0)then read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#Table: IMFMin')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, imf_by_) endif endif if(index(line,'#Element: bz')>0)then read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#Table: IMFMin')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, imf_bz_) endif endif if(index(line,'#Element: flow')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, sw_v_) endif if(index(line,'#Element: swVelX')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, sw_v_) endif if(index(line,'#Element: observed')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, f107_) endif if(index(line,'#Element: flux')>0)then call read_values call Insert_into_Indices_Array_Large(tmp, f107_) endif ! This has to be before the DMSP stuff because of the ! extra N and S on the end of the Element line. if (north) then if(index(line,'#Element: powerN')>0)then read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#Table: Hpinoaa')>0)then call read_values call merge_hpi_data endif endif else if(index(line,'#Element: powerS')>0)then read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#Table: Hpinoaa')>0)then call read_values call merge_hpi_data endif endif endif if(index(line,'#Element: power')>0)then read(lun_indices,'(a)', iostat = ierror ) line ! if(index(line,'#Table: Hpidmsp')>0)then call read_values call merge_hpi_data ! endif endif enddo if (npts_hpi > 0) then tmp = 0.0 do i=1,npts_hpi call time_real_to_int(ut_keep(i), itime) tmp(1:5,i) = itime(1:5) tmp(6,i) = data_keep(i) enddo call Insert_into_Indices_Array_Large(tmp, hpi_) endif if (Input_Coor_System == GSE) then idir = -1 do i=1,ntimes time_now = currenttime + (i-1) * dt call time_real_to_int(time_now, itime) iday = jday(itime(1),itime(2),itime(3)) call RECALC(itime(1),iday,itime(4),itime(5),itime(6)) xgse = Indices(i,imf_bx_) ygse = Indices(i,imf_by_) zgse = Indices(i,imf_bz_) call GSMGSE(XGSM,YGSM,ZGSM,XGSE,YGSE,ZGSE,idir) Indices(i,imf_bx_) = xgsm Indices(i,imf_by_) = ygsm Indices(i,imf_bz_) = zgsm enddo endif contains subroutine read_values logical :: done_inner real :: missing, missingPlusTol, missingMinusTol done_inner = .false. tmp = 0.0 missing = -1.0e32 do while (.not.done_inner) read(lun_indices,'(a)', iostat = ierror ) line if(index(line,'#yyyy-mm-dd')>0) done_inner = .true. if (ierror.ne.0) done = .true. if (index(line,'#Missing value:')>0) read(line,'(15x,f7.1)') missing if (index(line,'GSE')>0) Input_Coor_System = GSE enddo if (DebugLevel > 2) & write(*,*) "Missing value : ",missing missingPlusTol = missing + abs(missing)*0.01 missingMinusTol = missing - abs(missing)*0.01 if (ierror.eq.0) then done_inner = .false. i = 1 do while (.not.done_inner) read(lun_indices,'(f4.0,4(1x,f2.0),f7.1)', iostat = ierror ) tmp(1:6,i) if (ierror /= 0) then done_inner = .true. else if (tmp(6,i) < missingMinusTol .or. & tmp(6,i) > missingPlusTol) then i = i + 1 else tmp(1:6,i) = 0.0 endif endif enddo npts = i-1 end if end subroutine read_values subroutine merge_hpi_data itime(6) = 0 do i=1,npts itime(1:5) = tmp(1:5,i) call time_int_to_real(itime, ut_new(i)) data_new(i) = tmp(6,i) enddo if (npts_hpi == 0) then npts_hpi = npts ut_keep(1:npts) = ut_new(1:npts) data_keep(1:npts) = data_new(1:npts) else ut_tmp = ut_keep data_tmp = data_keep ut_keep = 0.0 data_keep = 0.0 j = 1 i = 1 k = 1 do while (i <= npts .or. j <= npts_hpi) if (i > npts) then ut_keep(k) = ut_tmp(j) data_keep(k) = data_tmp(j) k = k + 1 j = j + 1 else if (j > npts_hpi) then ut_keep(k) = ut_new(i) data_keep(k) = data_new(i) k = k + 1 i = i + 1 else if (ut_tmp(j) < ut_new(i)) then ut_keep(k) = ut_tmp(j) data_keep(k) = data_tmp(j) k = k + 1 j = j + 1 else ut_keep(k) = ut_new(i) data_keep(k) = data_new(i) k = k + 1 i = i + 1 endif endif endif enddo npts_hpi = npts_hpi + npts endif end subroutine merge_hpi_data end subroutine read_NGDC_Indices !============================================================================== subroutine read_NOAA_Indices use ModAMIE use ModIndices use ModFiles implicit none integer :: ierror, year, day, hour, minute, i integer :: idir, input_coor_system, iday, hpii logical :: done, found real :: tmp_hpi, normal ! One line of input character (len=100) :: line character (len=10) :: shortline real, dimension(6,MAX_NTIMES) :: tmp real*8, dimension(MAX_NTIMES) :: ut_new, ut_keep, ut_tmp real, dimension(MAX_NTIMES) :: data_new, data_keep, data_tmp real*8 :: time_now integer, dimension(7) :: itime open(lun_indices, file=cNOAAFile, status="old", iostat = ierror) if (ierror.ne.0) then call stop_mpi("Unable to open file "//cNGDCFile//". Stopping.") endif done = .false. do while (.not.done) read(lun_indices,'(a)', iostat = ierror ) line if (index(line, "Normalizing factor") > 0) done = .true. enddo read(lun_indices,*) year i = 1 done = .false. ierror = 0 tmp = 0.0 do while (.not.done) read(lun_indices,'(a10,i3,i2,i2,f8.1,i3,f8.3)', iostat = ierror ) & shortline, day, hour, minute, tmp_hpi, hpii, normal if (ierror.ne.0) then done = .true. else itime(1) = year itime(2) = 1 itime(3) = day itime(4) = hour itime(5) = minute itime(6) = 0 itime(7) = 0 call time_int_to_real(itime, time_now) if (time_now >= starttime-24.0*3600.0 .and. & time_now <= starttime + dt * (ntimes-1) + 24.0*3600.0) then tmp(1,i) = year tmp(2,i) = 1 tmp(3,i) = day tmp(4,i) = hour tmp(5,i) = minute tmp(6,i) = tmp_hpi if (index(shortline, "(N)") > 0 .and. north) i = i + 1 if (index(shortline, "(S)") > 0 .and. .not.north) i = i + 1 endif endif enddo write(*,*) "Done reading", i call Insert_into_Indices_Array_Large(tmp, hpi_) end subroutine read_NOAA_Indices !============================================================================== subroutine read_old_Indices use ModAMIE use ModIndices use ModFiles implicit none integer :: ierror, year, month, i, j logical :: done, found real, dimension(4,MAX_NTIMES) :: & power, aureb0, aehrly, f107d, bximf, byimf, bzimf, rkp, dst, swvel namelist /Data/ & power, aureb0, aehrly, f107d, bximf, byimf, bzimf, rkp, dst, swvel integer :: n_power, n_aureb0, n_aehrly, n_f107d integer :: n_bximf, n_byimf, n_bzimf, n_rkp, n_dst, n_swvel ! One line of input character (len=100) :: line open(lun_indices, file=cOldFile, status="old", iostat = ierror) if (ierror.ne.0) then call stop_mpi("Unable to open file "//cOldFile//". Stopping.") endif done = .false. found = .false. do while (.not.done) read(lun_indices,'(a)', iostat = ierror ) line if (ierror.ne.0) then done = .true. else if (index(line,'&Data') > 0) then backspace(lun_indices) found = .true. done = .true. endif endif enddo if (.not.found) then write(*,*) 'Error - Could not find the &Data' write(*,*) 'in the amiein file. Please check this file: ',cOldFile close(lun_indices) stop endif if (DebugLevel > 2 .and. me_world == 0) then write(*,*) "---------------------------------------------------------" write(*,*) "Warning:" write(*,*) "" write(*,*) "There is a small difference between the old namelist" write(*,*) "reads and the new namelist reads. If the code crashes" write(*,*) "with the following error : " write(*,*) "Invalid input for real editing" write(*,*) "Program terminated by fatal I/O error" write(*,*) "file: read_inputs.f90, line xxxx" write(*,*) "Abort" write(*,*) "" write(*,*) "The ending of the namelist should end with:" write(*,*) "/" write(*,*) "&" write(*,*) "Instead of:" write(*,*) "" write(*,*) "&" write(*,*) "---------------------------------------------------------" endif read(lun_indices,Data) close(lun_indices) call Insert_into_Indices_Array(power, hpi_) call Insert_into_Indices_Array(dst, dst_) call Insert_into_Indices_Array(rkp, kp_) call Insert_into_Indices_Array(bximf, imf_bx_) call Insert_into_Indices_Array(byimf, imf_by_) call Insert_into_Indices_Array(bzimf, imf_bz_) call Insert_into_Indices_Array(swvel, sw_v_) call Insert_into_Indices_Array(aehrly, ae_) call Insert_into_Indices_Array(f107d, f107_) ! Ignore aureb0 end subroutine read_old_Indices subroutine Insert_into_Indices_Array(array, label_) use ModAMIE use ModIndices real, dimension(4,MAX_NTIMES), intent(in) :: array real, dimension(MAX_NTIMES) :: index_value integer, intent(in) :: label_ real*8, dimension(MAX_NTIMES) :: index_time real*8 :: time_now, dt_min integer, dimension(7) :: itime !\ ! We have to assume that the dates coming in are in the current ! year and month !/ call time_real_to_int(currenttime, itime) itime(3:7) = 0 npts = Number_of_Points(array) if (npts == -1) then Indices(:,label_) = -1.0e32 elseif (npts == 0) then Indices(:,label_) = array(1,1) else do i = 1, npts itime(3) = array(1,i) itime(4) = array(2,i) itime(5) = array(3,i) itime(6) = 0 itime(7) = 0 call time_int_to_real(itime, index_time(i)) index_value(i) = array(4,i) enddo j = 1 do i=1,ntimes if (timelist(i) == -1.0e32) then time_now = currenttime + (i-1) * dt else time_now = timelist(i) endif do while ((j < npts).and. & (index_time(j) < time_now)) j = j + 1 enddo if ((j == npts .and. & index_time(j) <= time_now).or.(j==1)) then Indices(i,label_) = index_value(j) else dt_min = 1.0 - (index_time(j) - time_now) / & (index_time(j) - index_time(j-1) + 1.0e-6) Indices(i,label_) = dt_min * index_value(j) + & (1.0 - dt_min) * index_value(j-1) endif j = j - 2 if (j < 1) j = 1 enddo endif contains integer function Number_of_Points(array) real, dimension(4,MAX_NTIMES), intent(in) :: array integer :: n n = 0 if (maxval(abs(array)) == 0) then Number_of_Points = -1 else if (maxval(abs(array(:,2:MAX_NTIMES))) == 0 .and. & array(1,1) /= 0.0 .and. array(4,1) == 0.0) then Number_of_Points = 0 else Number_of_Points = 1 do while (maxval(abs(array(:,Number_of_Points))) > 0) Number_of_Points = Number_of_Points + 1 enddo Number_of_Points = Number_of_Points - 1 endif endif end function Number_of_Points end subroutine Insert_into_Indices_Array subroutine Insert_into_Indices_Array_Large(array, label_) use ModAMIE use ModIndices real, dimension(6,MAX_NTIMES), intent(in) :: array real, dimension(MAX_NTIMES) :: index_value integer, intent(in) :: label_ real*8, dimension(MAX_NTIMES) :: index_time real*8 :: time_now, dt_min integer, dimension(7) :: itime npts = Number_of_Points(array) if (npts == -1) then Indices(:,label_) = -1.0e32 elseif (npts == 0) then Indices(:,label_) = array(6,1) else do i = 1, npts itime(1) = array(1,i) itime(2) = array(2,i) itime(3) = array(3,i) itime(4) = array(4,i) itime(5) = array(5,i) itime(6) = 0 itime(7) = 0 call time_int_to_real(itime, index_time(i)) index_value(i) = array(6,i) enddo j = 1 do i=1,ntimes if (timelist(i) == -1.0e32) then time_now = currenttime + (i-1) * dt else time_now = timelist(i) endif do while ((j < npts).and. & (index_time(j) < time_now)) j = j + 1 enddo if ((j == npts .and. & index_time(j) <= time_now).or.(j==1)) then Indices(i,label_) = index_value(j) else dt_min = 1.0 - (index_time(j) - time_now) / & (index_time(j) - index_time(j-1) + 1.0e-6) Indices(i,label_) = dt_min * index_value(j) + & (1.0 - dt_min) * index_value(j-1) endif j = j - 2 if (j < 1) j = 1 enddo endif contains integer function Number_of_Points(array) real, dimension(6,MAX_NTIMES), intent(in) :: array integer :: n n = 0 if (maxval(abs(array)) == 0) then Number_of_Points = -1 else if (maxval(abs(array(:,2:MAX_NTIMES))) == 0 .and. & array(1,1) /= 0.0 .and. array(4,1) == 0.0) then Number_of_Points = 0 else Number_of_Points = 1 do while (maxval(abs(array(:,Number_of_Points))) > 0) Number_of_Points = Number_of_Points + 1 enddo Number_of_Points = Number_of_Points - 1 endif endif end function Number_of_Points end subroutine Insert_into_Indices_Array_Large