PROGRAM Setup_bf USE Rcm_mod_subs IMPLICIT NONE ! ! INTEGER (iprec) :: n, nbf, i, j, i_msm, j_msm, istat, ncol=80, Itype_tracer REAL (rprec) :: bbi, bbj, fstoff_tmp REAL (rprec), PARAMETER :: v_flag = 1.0E+30 ! LOGICAL :: ST integer :: iyear,iday,ihour,iminute,isec, imin_j_temp(jsize) REAL (rprec) :: p_dyn, g_2, g_3,ps, alpha_ratio,fkp REAL (rprec) :: parmod(10) = -1.0E+30 REAL (rprec) :: X_mgnp, y_mgnp, z_mgnp INTEGER (iprec) :: id EXTERNAL :: T01_s, Dip, T96_01 ! INTEGER :: ierr, numproc, myid, ierrorcode INTEGER :: nblk = 10 INTEGER, DIMENSION(:), ALLOCATABLE :: i1_blk, i2_blk, j1_blk, j2_blk REAL(rprec), ALLOCATABLE :: vm_blk (:,:), xmin_blk(:,:), & ymin_blk(:,:), bmin_blk(:,:) INTERFACE FUNCTION Pdynamic_sw (sw_den, sw_vel, ratio_den) ! sw_den is solar wind proton density in cm-3 ! sw_vel is solar wind bulk velocity in km/s ! ratio_den is ratio of alpha particles to proton number density ! Pdynamic_sw is rho*v**2 in nPa ! IF ratio_den is not known, set it to 0 or <0 and the default ! 4% value will be used. IMPLICIT NONE REAL, INTENT (IN) :: sw_den, sw_vel, ratio_den REAL :: Pdynamic_sw END FUNCTION Pdynamic_sw END INTERFACE INCLUDE "mpif.h" ! ! CALL MPI_Init (ierr) CALL MPI_Comm_Size (MPI_COMM_WORLD, numproc, ierr) CALL MPI_Comm_rank (MPI_COMM_WORLD, myid, ierr) ! nblk = numproc ALLOCATE (i1_blk(nblk), i2_blk(nblk), j1_blk(nblk), j2_blk(nblk)) i1_blk (:) = 1 i2_blk (:) = isize DO n = 1, nblk j1_blk(n) =1+(jsize/numproc)*(n-1) IF (n /= nblk) THEN j2_blk(n) =j1_blk(n) + jsize/numproc -1 ELSE j2_blk(n) = jsize END IF END DO ! ! ! WRITE (*,'(//A)') 'CHOOSE TYPE OF MAGNETIC FIELD MODEL TO USE:' ! WRITE (*,'(A)') ' 1--HILMER-VOIGT' ! WRITE (*,'(A)') ' 2--PURE DIPOLE' ! WRITE (*,'(A)') ' 3--T01-STORM' ! WRITE (*,'(A)') ' 4--T96' ! WRITE (*,'(A)') ' 5--T01' ! WRITE (*,'(A)',ADVANCE='NO') ' ENTER BFIELD_ID: ' ! READ (*,*) I_bfield_id ! ! Itype_tracer = 1 IF (Itype_tracer == 1) THEN IF (myid == 0) PRINT*,'USING TRACE ROUTINE FROM GEOPACK-2003' ELSE IF (Itype_tracer == 2) THEN IF (myid == 0) PRINT*,'USING TRACER ROUTINE FROM FRICUS' ELSE IF (myid == 0) PRINT*,'INVALID VALUE FOR ITYPE_TRACE' CALL MPI_Abort (MPI_COMM_WORLD, ierrorcode, ierr) END IF ! IF (myid == 0) THEN OPEN (UNIT = LUN_2, FILE = 'setup_bf.list', STATUS = 'REPLACE',& ACTION = 'WRITE', FORM = 'FORMATTED', ACCESS = 'SEQUENTIAL') END IF ! CALL Read_grid () ! CALL Rcm_read_params () ! ! ! Read number of bfields from 'definebf.dat': ! OPEN (LUN, FILE = 'inputs.dat', STATUS = 'OLD') READ (LUN,*) I_bfield_id READ (LUN,*) nbf = 0 DO READ (LUN,'(A)', END = 20) nbf = nbf + 1 END DO 20 CLOSE (LUN) ! ! ALLOCATE (ibtime (nbf)) ALLOCATE ( vm_blk(isize,jsize)) ALLOCATE (xmin_blk(isize,jsize)) ALLOCATE (ymin_blk(isize,jsize)) ALLOCATE (bmin_blk(isize,jsize)) ! ! ! ! Calculate the BF-matrices: ! OPEN (LUN_3, FILE = 'inputs.dat', STATUS = 'OLD') READ (LUN_3,*) READ (LUN_3,*) ! DO n = 1, nbf ! READ (LUN_3,*) ibtime(n),sw_n,sw_v,alpha_ratio,vdrop,fdst,fmeb,fkp,imf_bz,g_2,g_3 ! p_dyn = 1.67E-06*(sw_n*(sw_v)**2) ! sw ram pressure in nPa p_dyn = Pdynamic_sw (sw_n, sw_v, alpha_ratio) imf_by = 0.0 IF (myid == 0) print*,n, ibtime(n),fdst, imf_bz, g_2, g_3 parmod(1) = p_dyn parmod(2) = fdst parmod(3) = imf_by parmod(4) = imf_bz parmod(5) = g_2 parmod(6) = g_3 ps = 0.0 iyear = 1970 !2001 iday = 89 ihour = 0 iminute = 0 isec = 0 ! ! Need to find what the magnetopause standoff distance is in T01_s. ! Note that in T01_S, it is different from T01 and T96. All ! coefficients are in T01_s code, we deduce fstoff here as: ! Fstoff = 10.218556 / xappa, where xappa = (p_dyn/2)**0.15954 ! IF (i_bfield_id == 1) THEN ! T01_s fstoff = 10.218556 /((0.5*p_dyn)**0.15954) ELSE IF (i_bfield_id == 4) THEN ! T96 fstoff = 11.08 / ((0.5*p_dyn)**0.14) ELSE IF (i_bfield_id == 0) THEN ! DIPOLE fstoff = 12.0 ELSE PRINT*,'NOT IMPLEMENTED FSTOFF YET' CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) END IF ! label%intg(6) = ibtime(n) label%real(11) = p_dyn label%real(12) = 0.0 !fmeb label%real(13) = fstoff label%real(14) = fdst label%real(15) = 1.0 ! fclpse label%real(18) = imf_by label%real(19) = imf_bz label%real(09) = g_2 label%real(10) = g_3 label%intg (01)= n label%intg (20)= nbf label%intg (19)= I_bfield_id ! t01_S B_fields ! ! For each set of inputs, calculate matrices on grid: ! must zero arrays to use mpi_reduce (sum) later! ! vm_blk = 0.0 bmin_blk = 0.0 xmin_blk = 0.0 ymin_blk = 0.0 ! ! IF (I_BFIELD_id == 4 .AND. ibnd_type == 6) THEN ! CALL Get_boundary (boundary, bndloc) !that also gets imin_j ! i1_blk (:) = imin_j(1) - 2 ! i2_blk (:) = isize ! ! vm_blk (i1_blk(myid+1):i2_blk(myid+1),j1_blk(myid+1):j2_blk(myid+1)) = -1.0 ! xmin_blk (i1_blk(myid+1):i2_blk(myid+1),j1_blk(myid+1):j2_blk(myid+1)) = v_flag ! ymin_blk (i1_blk(myid+1):i2_blk(myid+1),j1_blk(myid+1):j2_blk(myid+1)) = v_flag ! bmin_blk (i1_blk(myid+1):i2_blk(myid+1),j1_blk(myid+1):j2_blk(myid+1)) = v_flag ! END IF ! IF (i_bfield_id == 1) THEN CALL Btrace_new(iyear,iday,ihour,iminute,isec, & Ri, Re, colat, aloct, ps, v_flag, & parmod, isize, jsize, SIZE(parmod), & vm_blk, bmin_blk, xmin_blk, ymin_blk, T01_s, Dip, n, & i1_blk(myid+1),i2_blk(myid+1),j1_blk(myid+1),j2_blk(myid+1), & Itype_tracer) ELSE IF (i_bfield_id == 4) THEN CALL Btrace_new(iyear,iday,ihour,iminute,isec, & Ri, Re, colat, aloct, ps, v_flag, & parmod, isize, jsize, SIZE(parmod), & vm_blk, bmin_blk, xmin_blk, ymin_blk, T96_01, Dip, n, & i1_blk(myid+1),i2_blk(myid+1),j1_blk(myid+1),j2_blk(myid+1), & Itype_tracer) ELSE IF (i_bfield_id == 0) THEN ! DIPOLE, ANALYTICAL EXPRESSIONS: DO i=i1_blk(myid+1),i2_blk(myid+1) DO j=j1_blk(myid+1),j2_blk(myid+1) xmin_blk (i,j) = Dipole_Bfield (colat(i,j), aloct(i,j), 1) ymin_blk (i,j) = Dipole_Bfield (colat(i,j), aloct(i,j), 2) bmin_blk (i,j) = Dipole_Bfield (colat(i,j), aloct(i,j), 3) vm_blk (i,j) = Dipole_Bfield (colat(i,j), aloct(i,j), 4) END DO END DO ELSE STOP 'I_BFIELD_ID NOT IMPLEMENTED' END IF ! CALL MPI_Barrier (MPI_COMM_WORLD, ierr) CALL MPI_Reduce ( vm_blk(:,:), vm, isize*jsize,MPI_REAL,& MPI_SUM,0,MPI_COMM_WORLD,ierr) CALL MPI_Reduce (xmin_blk(:,:), xmin, isize*jsize,MPI_REAL,& MPI_SUM,0,MPI_COMM_WORLD,ierr) CALL MPI_Reduce (ymin_blk(:,:), ymin, isize*jsize,MPI_REAL,& MPI_SUM,0,MPI_COMM_WORLD,ierr) CALL MPI_Reduce (bmin_blk(:,:), bmin, isize*jsize,MPI_REAL,& MPI_SUM,0,MPI_COMM_WORLD,ierr) ! IF (myid == 0) THEN ! ! Find the "actual" magnetopause standoff distance: j = 3 fstoff_tmp = 0.0 do_find_mgnp: DO i = isize, 2, -1 IF ((xmin(i-1,j) < xmin(i,j)) .OR. (xmin(i-1,j) > 0.9*v_flag)) THEN fstoff_tmp = xmin (i,j) EXIT do_find_mgnp END IF END DO do_find_mgnp IF (fstoff_tmp == 0 .AND. I_bfield_id > 0) then print*,'fstoff_tmp problem:', fstoff, fstoff_tmp pause 'fstoff_tmp problem' print*,xmin(:,j) pause fstoff_tmp = maxval(xmin(:,j)) end if IF (myid == 0) print*,' ',fstoff, fstoff_tmp ! ! Mark i=1 grid points as open field lines no matter what: ! xmin (1,:) = v_flag ymin (1,:) = v_flag vm (1,:) = -1.0 ! ! Find the outermost closed field line along each J-line, ! and set all points poleward of it to the values of that ! "last closed field line" (thus, get rid of V_flag values) ! (don't want them in Get_boundary in the RCM!) ! do j = 1, jsize imin_j_temp(j) = -1 do_i: do i = isize, 2, -1 if ((abs(xmin(i,j)) < v_flag .AND. & abs(xmin(i-1,j)) >= v_flag*0.9) ) then imin_j_temp (j) = i exit do_i end if end do do_i if (imin_j_temp(j) < 0) then print*,xmin(:,j) print*,'nbf',n STOP 'negative imin_j_temp' end if xmin(1:imin_j_temp(j)-1,j) = xmin(imin_j_temp(j),j) ymin(1:imin_j_temp(j)-1,j) = ymin(imin_j_temp(j),j) ! vm (1:imin_j_temp(j)-1,j) = vm (imin_j_temp(j),j) end do ! ! Place rcm boundary on the eq. mapping: and check for open field lines inside: ! CALL Get_boundary (boundary, bndloc) !that also gets imin_j DO j = 1, jsize DO i = imin_j(j), isize IF (xmin(i,j) > v_flag .OR. ymin(i,j) > v_flag & .OR. vm(i,j) > v_flag) THEN PRINT*,'OPEN FIELD LINE INSIDE MODELING REGION' PRINT*,I,J,XMIN(I,J), YMIN(I,J), NBF stop END IF END DO END DO ! END IF ! ! ! ! CALL Circle (vm) ! CALL Circle (bmin) ! CALL Circle (xmin) ! CALL Circle (ymin) ! IF (myid == 0) THEN CALL Outp (vm, 1,isize,1,1,jsize,1,10.0_rprec,label%intg,'VM ', LUN_2, ncol) CALL Outp (xmin,1,isize,1,1,jsize,1,10.0_rprec,label%intg,'XMIN', LUN_2, ncol) CALL Outp (ymin,1,isize,1,1,jsize,1,10.0_rprec,label%intg,'YMIN', LUN_2, ncol) CALL Outp (bmin,1,isize,1,1,jsize,1,1.E+2_rprec,label%intg,'BMIN',LUN_2, ncol) ! IF (n == 1) THEN ST = .TRUE. ELSE ST = .FALSE. END IF CALL Write_array ('rcmbmin_inp', n, label, bmin, SETUP = ST, ASCI = asci_flag) CALL Write_array ('rcmvm_inp', n, label, vm, SETUP = ST, ASCI = asci_flag) CALL Write_array ('rcmxmin_inp', n, label, xmin, SETUP = ST, ASCI = asci_flag) CALL Write_array ('rcmymin_inp', n, label, ymin, SETUP = ST, ASCI = asci_flag) END IF ! END DO CLOSE (UNIT = LUN_3) CLOSE (UNIT = LUN_2 ) ! DEALLOCATE (vm_blk, xmin_blk, ymin_blk, bmin_blk) DEALLOCATE (i1_blk, i2_blk, j1_blk, j2_blk) CALL MPI_Finalize (ierr) STOP END PROGRAM Setup_bf ! ! ! SUBROUTINE Btrace_new (iyear, iday, ihour, imin, isec, & R_i, R_e, colat, aloct, pss, V_flag, & parmod, imax, jmax, ipdim, & vm, bmin, xmin, ymin, & exname, inname, n_bf, i1,i2, j1, j2, & itype_trace) USE Rcm_mod_subs, ONLY : rprec, iprec IMPLICIT NONE INTEGER, INTENT (IN) :: iyear, iday, ihour, imin, isec, & imax, jmax, ipdim, n_bf, i1, i2, j1,j2,& itype_trace REAL(rprec), INTENT (IN) :: R_i, R_e, pss, parmod(ipdim), V_flag REAL(rprec), DIMENSION (imax,jmax), INTENT (IN) :: colat, aloct REAL(rprec), DIMENSION (Imax,jmax), INTENT (OUT):: vm, bmin, xmin, ymin EXTERNAL :: exname, inname ! ! LOGICAL :: L_mgnp COMMON /stanislav_1/ L_mgnp ! REAL :: aa(10), sps, cps, bb(3), ps, cc(19) COMMON /geopack1/ aa,sps,cps,bb,ps,cc ! INTEGER :: iopt, i, j, m, n_open, n, ierrorcode, ierr, iopen REAL :: xx(10000), yy(10000), zz(10000), btot(10000), ftv(10000), pi, & dir, r0, gphi, gthet, grr, glat, xf, yf, zf, & bx, by, bz, bbx, bby, bbz, sum, h, xgsm, ygsm, zgsm, tmp,& xmin_p, ymin_p, zmin_p ! real x_min,x_max,y_min,y_max,z_max,z_min, rlim common /limits/ x_min,x_max,y_min,y_max,z_max,z_min ! INCLUDE "mpif.h" ! ! ! * * field line tracing from ionospheric grid specified through colat ! and aloct to compute flux tube volumes. ! ! ! * * paramerers: colat, aloct = ionospheric grip point coordinates ! ! parmod = input vector ! parmod(1) = solar wind pressure p_dyn (nPa) ! parmod(2) = Dst (nT) ! parmod(3) = imf_by (nT) ! parmod(4) = imf_bz (nT) ! parmod(5) = g_2 index ! parmod(6) = g_3 index ! (see Tsyganenko, Singer, and Kasper [2003] for def. of g_2 and g_3) ! ! ps = dipole tilt angle ! ! IF (i1 < 1 ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (i1 > imax ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (i2 < 1 ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (i2 > imax ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (j1 < 1 ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (j1 > jmax ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (j2 < 1 ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) IF (j2 > jmax ) CALL MPI_abort (MPI_COMM_WORLD, ierrorcode, ierr) ! pi = acos(-1.000000) ! vm (i1:i2,j1:j2) = -1.0 bmin(i1:i2,j1:j2) = v_flag xmin(i1:i2,j1:j2) = v_flag ymin(i1:i2,j1:j2) = v_flag ! CALL Recalc (iyear,iday,ihour,imin,isec) ps = pss; sps = SIN(ps); cps = COS(ps) ! rlim = 40.0 ! limit the tracing region within 40 Re ! this is for geopack x_min = -40.0; x_max = +15.0; y_min = -20.0; y_max = +20.0; z_min=-15.0;z_max=15.0 IF (itype_trace == 1) THEN ! we trace from northern hemispher dir = 1.0 ! this is tsyganenko's convention ELSE dir = -1.0 ! this is toffoletto's convention END IF r0 = R_i/R_e ! landing point on Earth's surface iopt = 0 ! dummy unused parameter ! do_j_loop: DO j = j1, j2, 1 ! n_open = 0 ! do_i_loop: DO i = i2, i1, -1 ! print*,'calling btrace with i,j=',i,j ! gphi = aloct(i,j) gthet = colat(i,j) glat = pi/2. - gthet grr = R_i / R_e ! xgsm = grr*sin(gthet)*cos(gphi) ygsm = grr*sin(gthet)*sin(gphi) zgsm = grr*cos(gthet) ! L_mgnp = .TRUE. ! Flag if a point is inside magnetosphere ! false means either outside magnetopause ! or in the "boundary layer" ! IF (itype_trace == 1) THEN ! use geopack tracer: ! CALL Trace (xgsm, ygsm,zgsm, dir, rlim, r0, iopt, parmod, & exname, inname, xf, yf, zf, xx, yy, zz, m) ! IF (m>998 .OR. .not.L_mgnp .OR. & SQRT(xf**2+yf**2+zf**2)>1.6) THEN !OPEN FIELD n_open = n_open + 1 IF (n_open > 4) EXIT do_i_loop ! ELSE ! closed field line, proceed with integration ! DO n=1,m ! Get B-field values along a field line: CALL Exname (iopt,parmod,ps,xx(n),yy(n),zz(n),bx,by,bz) CALL Inname (xx(n), yy(n), zz(n), bbx, bby, bbz) btot(n) = SQRT ((bx+bbx)**2 + (by+bby)**2 + (bz+bbz)**2) END DO ! if (i==imax .and. j== 3) then ! do n =1 ,m ! print*,n,xx(n),yy(n),zz(n),btot(n) ! end do ! end if ! sum = 0.0 DO n = 1, m-1 ! Integrate flux-tube volume: h = SQRT((xx(n+1)-xx(n))**2 + & (yy(n+1)-yy(n))**2 + & (zz(n+1)-zz(n))**2) sum = sum + 2.*h/(btot(n+1)+btot(n)) END DO vm(i,j) = sum**(-2.0/3.0) ! ! Find equatorial crossings (note: cannot use min B, because ! B minimizes off the z=0 plane): ! do_find_eq_crossing: DO n = 1, m-1 IF (zz(n) > 0.0 .AND. zz(n+1) < 0.0) THEN xmin (i,j) = (xx(n)+xx(n+1))/2.0 ymin (i,j) = (yy(n)+yy(n+1))/2.0 bmin (i,j) = (btot(n)+btot(n+1))/2.0 EXIT do_find_eq_crossing END IF END DO do_find_eq_crossing ! END IF ! ELSE IF (itype_trace == 2) THEN ! use mapd: ! CALL Mapd (xgsm, ygsm, zgsm, xf, yf, zf, & xmin_p, ymin_p, zmin_p, dir, iopen, & xx, yy, zz, btot, ftv, SIZE(xx), m, ierr, & parmod, exname, inname, grr) IF (iopen == -1) THEN vm (i,j) = ftv(m)**(-2.0/3.0) xmin(i,j)= xmin_p ymin(i,j)= ymin_p END IF ELSE STOP 'ITYPE_TRACE ILLEGAL' END IF ! END DO do_i_loop END DO do_j_loop ! RETURN END SUBROUTINE Btrace_new ! ! ! ! !----------------------------------------------------- ! SUBROUTINE Tracf (xa,ya,za,fa,xn,yn,zn,fn,er,ep,h,ierr,& parmod, exname, inname) ! !================================================================ ! purpose: ! field line tracing program, based on the Voigt first order tracer ! keeps track of changes in flux tube volume to adjust stepsizes ! Take one step along a field line, and compute contribution to ! flux-tube volume. ! ! inputs: ! xa,ya,za: starting location ! fa: starting flux tube volume value ! er: error tolerance for tracer ! ep: direction (along =1, against -1) along field line ! h: stepsize (new value will be outputted) ! outputs: ! xf,yf,zf: ending location ! vf: end flux tube volume value ! h: new stepsize ! ierr: error flag ! !================================================================ IMPLICIT NONE REAL, INTENT (IN) :: xa, ya, za, fa, er, ep REAL, INTENT (OUT):: xn, yn, zn, fn REAL, INTENT (IN) :: parmod(10) INTEGER, INTENT (OUT) :: ierr REAL, INTENT (IN OUT) :: h EXTERNAL :: exname, inname ! real dx1,dy1,dz1,bf1 real dxt,dyt,dzt,bft real xn1,yn1,zn1,fn1 real xn2,yn2,zn2,fn2 real xn3,yn3,zn3,fn3 real dx2,dy2,dz2,bf2,df2 real rtes1,rtes2,rtes3,rtesf,rtes real delx,dely,delz,delf real h1,h2 real, PARAMETER :: ftv_fact=200.0 ! REAL a(15),psi,aa(10),ds3,bb(8) common /geopack1/ a,psi,aa,ds3,bb ! ds3 = 1.0 CALL Rhand_b (xa, ya, za, dx1, dy1, dz1, bf1, 0, parmod, exname, inname) ! 1 continue h1 = h * ep h2 = 0.5 * h1 xn1 = xa + h1*dx1 yn1 = ya + h1*dy1 zn1 = za + h1*dz1 CALL Rhand_b (xn1, yn1, zn1, dxt, dyt, dzt, bft, 0, parmod, exname, inname) fn1 = fa + 2.0*abs(h1/(bf1+bft)) xn2 = xa + h2*dx1 yn2 = ya + h2*dy1 zn2 = za + h2*dz1 CALL Rhand_b (xn2, yn2, zn2, dxt, dyt, dzt, bft, 0, parmod, exname, inname) fn2 = fa + 2.0*abs(h2/(bf1+bft)) CALL Rhand_b (xn2, yn2, zn2, dx2, dy2, dz2, bf2, 0, parmod, exname, inname) xn3 = xn2 + h2*dx2 yn3 = yn2 + h2*dy2 zn3 = zn2 + h2*dz2 CALL Rhand_b (xn3, yn3, zn3, dxt, dyt, dzt, bft, 0, parmod, exname, inname) fn3 = fn2 + 2.0*abs(h2/(bf2+bft)) delx = xn1 - xn3 dely = yn1 - yn3 delz = zn1 - zn3 delf = fn1 - fn3 rtes1 = delx*delx rtes2 = dely*dely rtes3 = delz*delz ! decrease tolerance of ftv by ftv_fact rtesf = delf*delf*ftv_fact rtes = sqrt(rtes1 + rtes2 + rtes3 + rtesf) ! rtes = rtes/sqrt(xa**2+ya**2+za**2) if (rtes-er) 4,4,2 4 continue xn = xn3 - delx yn = yn3 - dely zn = zn3 - delz fn = fn3 - delf ! bb = bf2 if(rtes-(0.25*er)) 7,6,6 7 h = 2.0*h 6 go to 5 2 continue h = 0.5*h go to 1 5 continue RETURN END SUBROUTINE Tracf ! ! ! SUBROUTINE Mapd (xs, ys, zs, xf, yf, zf, xe, ye, ze, dir, iopen, & xx, yy, zz, bbb, ftv, imax, ival, ierr, & parmod, exname, inname, grr) ! ! Trace a magnetic field line. Based on Voigt-Toffoletto code. ! ! inputs: ! xs,ys,zs = starting location ! xf,yf,zf = ending location, in the ionosphere (closed) or ! some defined boundary (open) ! xe,ye,ze = equatorial crossing point (if there is one) ! dir = direction for mapping (either -1/+1) wrt field line ! imax = max size of 3 1D arrays that have field line coordinates ! ! outputs: ! ! iopen = flag if point is connected to inosphere (-1) or bndy (0/1) ! xx,yy,zz = coordinates of a single field line (size imax) ! bbb = magnitude of field along field line track ! ftv = flux tube volume as a function of field line location ! ival = ending index where field line tracing stopped ! ierr = error flag (=0 ok) !============================================================== ! IMPLICIT NONE INTEGER, INTENT (IN) :: imax REAL, INTENT (IN) :: parmod(10) INTEGER, INTENT (OUT):: ival, ierr, iopen REAL, INTENT (IN) :: xs, ys, zs, dir, grr REAL, INTENT (OUT):: xf, yf, zf, xe, ye, ze REAL, INTENT (OUT):: xx(imax), yy(imax), zz(imax), bbb(imax),ftv(imax) EXTERNAL :: exname, inname real dx1,dy1,dz1,bb,h,ep,er real, PARAMETER :: er1=1.0E-4,er2=1.0E-4 real dr,bex,r,rrr real xn,yn,zn,xm,ym,zm,fn,fm,ff real bbx,bby,bbz,bbt real rn,rm real bbmin logical :: L_mgnp common /stanislav_1/ L_mgnp ! ! REAL a(15),psi,aa(10),ds3,bb_b(8) common /geopack1/ a,psi,aa,ds3,bb_b ! ! integer icr1,i ! common /bstep/ dx1,dy1,dz1,bb,h,ep,er ! common /trac_error/ er1,er2 ! real xmin,xmax,ymin,ymax,zmax,zmin common /limits/ xmin,xmax,ymin,ymax,zmax,zmin ! data bex /10.0/ data bbmin /10.0/ ! ierr = 0 ival = 0 xe = 0.0 ye = 0.0 ze = 990.0 icr1 = 1 er = er1 iopen = -1.0 h = 0.5 ! used to be 0.1 ep = dir ! reset ftv do i=1,imax ftv(i) = 0.0 end do ! xm = xs ym = ys zm = zs fm = 0.0 ! xx(1) = xm yy(1) = ym zz(1) = zm ftv(1) = 0.0 ! ds3 = 1.0 CALL Rhand_b (xs, ys, zs, bbx, bby, bbz, bbt, 0, parmod, exname, inname) bbb(1) = bbt h = 0.1 ! DO i = 2,imax ! CALL Tracf (xm,ym,zm,fm,xn,yn,zn,fn,er,ep,h,ierr, & parmod, exname, inname) CALL Rhand_b (xn, yn, zn, bbx, bby, bbz,bbt, 0, parmod, exname, inname) r = SQRT(xn*xn + yn*yn + zn*zn) ! ! safety switch ! dr = SQRT((xm-xn)**2+(ym-yn)**2+(zm-zn)**2) ! IF (dr < 1.0e-3 .AND. r <= 2.0) RETURN ! bb = bbt bbb(i) = bb ! ! set limits on h IF (h > 0.1) h=0.1 ! ! check to see if z=0 is crossed; look for crossing point for bfield ! IF (zm*zn <= 0.0) THEN IF (zm-zn.ne.0.0) THEN xe = (xn*zm-xm*zn)/(zm-zn) ye = (yn*zm-ym*zn)/(zm-zn) ze = 0.0 ELSE xe = xn ye = yn ze = zn END IF END IF ! .......need to be careful about what error limits exist and where... ! ! condition open field line !-----------note er is now input IF (r < 4.0) THEN er = er2 ELSE er = er1 END IF ! for ftv calc IF (bbt < bbmin) THEN er = er2 END IF ! ! condition closed field line: IF (r < grr .AND. i > 10) THEN iopen = - 1 ! reset to r = grr by interpolation ! rn=sqrt(xn**2+yn**2+zn**2) rm=sqrt(xm**2+ym**2+zm**2) xf = (xm*(grr-rn)+xn*(rm-grr))/(rm-rn) yf = (ym*(grr-rn)+yn*(rm-grr))/(rm-rn) zf = (zm*(grr-rn)+zn*(rm-grr))/(rm-rn) ff = (fm*(grr-rn)+fn*(rm-grr))/(rm-rn) xn = xf yn = yf zn = zf fn = ff xx(i) = xf yy(i) = yf zz(i) = zf ftv(i) = ff RETURN END IF ! ! condition don't know ! if((.not.L_mgnp).or.(xn.le.xmin).or. & (yn.le.ymin).or. & (zn.le.zmin).or. & (xn.ge.xmax).or. & (yn.ge.ymax).or. & (zn.ge.zmax).or. & (bb.lt.0.1)) then iopen = 0 xf = xn yf = yn zf = zn ff = fn xx(i) = xf yy(i) = yf zz(i) = zf ftv(i) = ff return END IF xx(i) = xn yy(i) = yn zz(i) = zn ftv(i) = fn ival = i xm = xn ym = yn zm = zn fm = fn END DO write(6,*)' i greater than ',imax write(6,*)'xn,yn,zn,fn',xn,yn,zn,fn STOP RETURN END SUBROUTINE Mapd ! ! ! !==================================================================================== ! subroutine rhand_b (x,y,z,r1,r2,r3,btot, iopt,parmod,exname,inname) ! ! Identical to the GEOPACK routine RHAND, but also returns total B-field. ! ! calculates the components of the right hand side vector in the geomagnetic field ! line equation (a subsidiary subroutine for the subroutine step) ! ! last modification: march 31, 2003 ! ! author: n. a. tsyganenko ! dimension parmod(10) ! ! exname and inname are names of subroutines for the external and internal ! parts of the total field ! common /geopack1/ a(15),psi,aa(10),ds3,bb(8) call exname (iopt,parmod,psi,x,y,z,bxgsm,bygsm,bzgsm) call inname (x,y,z,hxgsm,hygsm,hzgsm) bx=bxgsm+hxgsm by=bygsm+hygsm bz=bzgsm+hzgsm b=ds3/sqrt(bx**2+by**2+bz**2) r1=bx*b r2=by*b r3=bz*b btot = SQRT (bx**2+by**2+bz**2) return end ! !===================================================================================