MODULE FIELD_LINE IMPLICIT NONE REAL :: theta, phi, height, R_zero, H_bottom, H_top, height_low, height_high,& theta_low, theta_high, besu = 31000., R_e=6375.0, R_i = 6475.0 END MODULE FIELD_LINE ! MODULE Conditions IMPLICIT NONE INTEGER :: daynr REAL :: f107, f107a, f107_12m, ap(7) END MODULE Conditions ! MODULE Constants IMPLICIT NONE ! PHYSICAL CONSTANTS REAL, PARAMETER :: q = 1.6021892E-19, & m_e = 9.1094E-31, & m_h = 1.6726E-27, & m_no_plus = 30.0*m_h, & m_o_plus = 16.0*m_h, & m_o2_plus = 32.0*m_h ! ! INTEGRATION PARAMETERS ! The following lines describe parameters for integration routines: ! - EPS_INT_1 is the fractional accuracy of computing integrals, ! - EPS_INT_2 is the accuracy for double integralsr,. ! - STEP_H_INIT is the initial step for height integration ! - STEP_H_MIN is the minimul allowed step for height integration ! - both are in km. Note: STEP_H_MIN can be set to zero. ! - STEP_T_INIT is the initial step for colatitude integration ! - STEP_T_MIN is the minimum allowed step for colatitude integration ! - both are in radians. ! Note: don't be greedy: the smaller the eps_int, the slower the code! INTEGER, PARAMETER :: nsteps_int_eq = 5 ! # of equally-sized steps for Equat. intgration REAL, PARAMETER :: eps_int_1 = 1.0E-2, & eps_int_2 = 1.0E-2, & step_h_init = 0.1, & step_h_min = 0.0001, & step_t_init = 0.0001, & step_t_min=0.0000001, & pi = 3.141592654 ! END MODULE Constants ! ! ! SUBROUTINE Integrate_theta (theta_low, theta_high, & int_sig_p, int_sig_h, int_pwn, & int_pwe, int_hwn, int_hwe) USE Constants IMPLICIT NONE REAL, INTENT (IN) :: theta_low, theta_high REAL, INTENT (OUT) :: int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe ! !-------------------------------------------------------------------------------- ! Subroutine to compute field-aligned integrals of pedersen ! and hall conductivities in the ionosphere ! The integral is along a (dipole) magnetic field line from ! colatitude theta_low to colatitude theta_high (both in rads). ! Stanislav, 9/25/97 ! This routine is used if the field line of integration is not too vertical. ! Integration is only in northern hemisphere, but the integrand ! is a sum of quantities in both hemispheres at conjugate points. ! ! Now, here is how integrals are computed: ! they are solutions of a system of 1st-order ODEs ! dy/dx = F(x) at x=b ! with the initial condition y(x=a)=0 ! where a and b are limits of integration, and y is the vector of ! unknown functions. y(1) is int_sig_p, y(2)=int_sig_h, etc (see below). !---------------------------------------------------------------------------- ! ! INTEGER :: i, nok, nbad INTEGER, PARAMETER :: num_int=6 REAL :: y_vector(NUM_int) ! EXTERNAL Derivs_theta, Rk_adapt_step ! ! ! Y_VECTOR contains 6 field-line integrated quantities as follows: ! Y_VECTOR(1) is int_sig_p, ! Y_VECTOR(2) is int_sig_h, ! Y_VECTOR(3) is int_pwn, ! Y_VECTOR(4) is int_pwe, ! Y_VECTOR(5) is int_hwn, ! Y_VECTOR(6) is int_hwe. ! none of the quantities above has sini factors (not yet). ! function derivs_theta returns a vector of RHS for ODEs (that is, ! a vector of integrands for our integrals; if you replace ! or modify the function, you need to preserve the calling ! sequence (number and type of arguments). ! y_vector(:)=0. CALL Ode_int (y_vector, num_int, theta_low, theta_high, eps_int_1, & step_t_init, step_t_min, nok, nbad, & Derivs_theta, Rk_adapt_step) int_sig_p = y_vector(1) int_sig_h = y_vector(2) int_pwn = y_vector(3) int_pwe = y_vector(4) int_hwn = y_vector(5) int_hwe = y_vector(6) RETURN END SUBROUTINE Integrate_theta ! ! ! ! SUBROUTINE Integrate_height (height_low, height_high, & int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe) USE Constants IMPLICIT none REAL, INTENT (IN) :: height_low, height_high REAL, INTENT (OUT) :: int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe ! !_____________________________________________________________________ ! ! subroutine to compute field-aligned integrals of pedersen ! and hall conductivities in the ionosphere ! Stanislav, 9/25/97 ! Use this subroutine when the field line is close to vertical, ! since integrating with respect to colatitude would lead to ! integrating large functions over small intervals. ! See the description of INTEGRATE_THETA too. !_____________________________________________________________________ ! ! INTEGER i, nok, nbad INTEGER, PARAMETER :: num_int=6 REAL :: y_vector(NUM_int) ! EXTERNAL Derivs_height, Rk_adapt_step ! ! y_vector(:)=0. CALL Ode_int (y_vector, num_int, height_low, height_high, eps_int_1, & step_h_init, step_h_min, nok, nbad, & Derivs_height, Rk_adapt_step) int_sig_p = y_vector(1) int_sig_h = y_vector(2) int_pwn = y_vector(3) int_pwe = y_vector(4) int_hwn = y_vector(5) int_hwe = y_vector(6) RETURN END SUBROUTINE Integrate_height ! ! ! ! SUBROUTINE Integrate_theta_eq (theta_low, theta_high, ss, sw) USE Constants IMPLICIT NONE REAL, INTENT (IN) :: theta_low, theta_high REAL, INTENT (OUT) :: ss, sw ! INTEGER i INTEGER, PARAMETER :: num_int=2 REAL :: y_vector_i(NUM_int), y_vector_e(num_int) ! EXTERNAL Derivs_theta_eq ! ! ! ! Y_VECTOR contains 2 field-line integrated quantities as follows: ! Y_VECTOR(1) is ss ! Y_VECTOR(2) is sw ! ! y_vector_i(:)=0. CALL Rkdumb (y_vector_i, y_vector_e, num_int, theta_low, theta_high, & nsteps_int_eq, Derivs_theta_eq) ss = y_vector_e(1) sw = y_vector_e(2) RETURN END SUBROUTINE Integrate_theta_eq ! ! ! SUBROUTINE Derivs_theta (n, theta_local, y, dydx) USE Field_line USE Conditions USE Constants IMPLICIT NONE ! ! Function that computes the integrands (dydx) when the integration ! variable is colatitude (theta). y is the value of integral; it is ! not needed, but is required to preserve the calling sequence. ! Stanislav, last mod. on 10/06/97 ! INTEGER, INTENT (IN) :: n REAL, INTENT (IN) :: theta_local, y(n) REAL, INTENT (OUT) :: dydx(n) ! ! REAL sigma_p,sigma_h,v_n,v_e, & sigma_p_conj,sigma_h_conj,v_n_conj,v_e_conj, & geom_factor ! IF (n /= 6) STOP 'N SHOULD BE 6' ! ! 1. The current point in the Northern hemisphere: ! For given theta and R_0 of the field line, compute: ! theta = theta_local height = R_zero*(SIN(theta))**2-R_e CALL Sigma_tensor (theta, phi, height, daynr, f107, f107a, f107_12m, ap, & sigma_p, sigma_h, v_n, v_e) ! ! ! 2. Do the same at the conjugate point (Southern hemisphere): ! theta = 4.0*ATAN(1.0)-theta_local CALL Sigma_tensor (theta, phi, height, daynr, f107, f107a, f107_12m, ap, & sigma_p_conj, sigma_h_conj, v_n_conj, v_e_conj) ! ! ! 3. Compute the geometrical factor ds/dtheta (same in both hemispheres): ! geom_factor = Dsdtheta (theta, R_zero) ! ! ! 4. RHS's: dydx(1)=(sigma_p+sigma_p_conj)*geom_factor dydx(2)=(sigma_h+sigma_h_conj)*geom_factor dydx(3)=(sigma_p*v_n+sigma_p_conj*v_n_conj)*geom_factor dydx(4)=(sigma_p*v_e+sigma_p_conj*v_e_conj)*geom_factor dydx(5)=(sigma_h*v_n+sigma_h_conj*v_n_conj)*geom_factor dydx(6)=(sigma_h*v_e+sigma_h_conj*v_e_conj)*geom_factor RETURN CONTAINS ! ! ! REAL FUNCTION Dsdtheta (theta_local, R_zero_local) IMPLICIT none REAL, INTENT (IN) :: theta_local, R_zero_local ! ! Derivative of S (field line length) with respect to colatitude: ! Factor of 1.E+3 here is since R_0 is in km, but we want ! conductances in SI units (that is, meters) ! Dsdtheta = (1.E+3)*R_zero_local*SIN(theta_local)* & SQRT(1.+3.*(COS(theta_local))**2) RETURN END FUNCTION Dsdtheta END SUBROUTINE Derivs_theta ! ! ! SUBROUTINE Derivs_height (n, height_local, y, dydx) USE FIELD_LINE USE Conditions USE Constants IMPLICIT none ! ! Function returns the derivative dS/dheight ! y is not needed, but required to preserve the calling sequence. ! INTEGER, INTENT (IN) :: n REAL, INTENT (IN) :: height_local, y(n) REAL, INTENT (OUT) :: dydx(n) ! ! REAL :: sigma_p,sigma_h,v_n,v_e, & sigma_p_conj,sigma_h_conj,v_n_conj,v_e_conj REAL :: geom_factor ! IF (n /= 6) STOP 'N IS NOT 6' ! ! 1. Current point of integration in northern hemisphere: ! here take theta= pi/2.0) then if ((theta_in-pi/2.0).le.0.000001) then print*,'setting derivs_eq to zero' dydx(1) = 0. dydx(2) = 0. return else write(6,'(2f15.10)')theta_in, pi/2.0 stop 'theta too large!!!' endif END IF ! R_zero = (R_e+h_bottom)/(SIN(theta_in)**2) height_low = h_bottom theta_low = theta_in IF ((R_zero-R_e) > H_top) THEN height_high = h_top theta_high = ASIN(SQRT((R_e+H_top)/R_zero)) ELSE height_high = R_zero-R_e theta_high = pi/2.0 END IF CALL Integrate_theta (theta_low, theta_high, & int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe) sin_i_local = 2.*COS(theta_in)/SQRT(1.+3.*(COS(theta_in))**2) B_local = BESU*(1.E-9)*SQRT(1.+3.*(COS(theta_in))**2)/ & (((R_e+h_bottom)/R_e)**3) dydx(1) = (int_sig_p+int_sig_h**2/int_sig_p)* & sin_i_local/SIN(theta_in) dydx(2) = ( -int_pwn+int_hwe & -(int_hwn+int_pwe)*int_sig_h/int_sig_p)*(R_i*1.E+3)*& sin_i_local*B_local ! y = 0.0 RETURN END SUBROUTINE Derivs_theta_eq ! ! SUBROUTINE Rkdumb (vstart,vend,nvar,x1,x2,nstep,derivs) ! ! Numerical Recipes subroutine to tabulate integrated functions. ! Computes the integral of vector of func(x) from x1 to x2 using initial ! values in vstart(1:nvar) and fourth-order Runge-Kutta for each step, ! the number of steps is nstep of equal size. The integrands are ! computed in the function derivs. The final values at x2 are ! returned in vend(1:nvar). ! ! Modified to return only the values at the end of the interval. ! To use all intermediate values, uncomment common block "path". ! IMPLICIT NONE ! NMAX is the max number of functions in the vector. INTEGER nstep,nvar,NMAX,NSTPMX PARAMETER (NMAX=50,NSTPMX=200) REAL x1,x2,vstart(nvar),vend(nvar),xx(NSTPMX),y(NMAX,NSTPMX) ! INTERFACE ! SUBROUTINE Derivs (n, x, y, dydx) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x ! REAL, INTENT (OUT) :: y(n), dydx(n) ! END SUBROUTINE Derivs ! END INTERFACE INTEGER i,k REAL h,x,dv(NMAX),v(NMAX) EXTERNAL Derivs v(1:nvar) = vstart(1:nvar) y(1:nvar,1) = v(1:nvar) xx(1)=x1 x=x1 h=(x2-x1)/nstep DO k = 1, nstep CALL Derivs (nvar, x, v, dv) CALL Rk4 (v, dv, nvar, x, h, v, Derivs) IF (x+h == x) STOP 'stepsize not significant in rkdumb' x = x + h xx(k+1)=x y(1:nvar,k+1)=v(1:nvar) END DO vend (1:nvar) = y(1:nvar,nstep+1) RETURN END SUBROUTINE Rkdumb ! ! !--------------------------------------------------------------------------------------- ! ! ! ! SUBROUTINE Ode_int (ystart, nvar, x1, x2, eps, h1, hmin, nok, nbad, Derivs, Rk_adapt_step) ! ! Runge-Kutta driver routines (after NR) with adaptive stepsize control. ! Integrate a set of ODEs (nvar equations) dy/dx = Derivs (y,x) from x1 to x2 ! with accuracy EPS. YSTART on input has y(x1), on output--y(x2). ! h1 is initial step, hmin is min allowed (>= 0). NOK and NBAD is number of ! good and bad steps taken. Derivs is function for RHS. Rkqs is subroutine ! to take one time step (now called Rk_adapt_step). ! IMPLICIT NONE INTEGER, INTENT (IN) :: nvar REAL, INTENT (IN) :: x1, x2, eps, h1, hmin INTEGER, INTENT (OUT) :: nok, nbad REAL, INTENT (IN OUT) :: ystart (nvar) ! INTERFACE ! SUBROUTINE Derivs (n,x,y,dxdy) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x, y(n) ! REAL, INTENT (OUT) :: dxdy(n) ! END SUBROUTINE Derivs ! END INTERFACE ! INTERFACE ! SUBROUTINE Rk_adapt_step (y, dydx, n, x, htry, eps, yscal, hdid, hnext, Derivs) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: htry, eps, yscal(n), dydx(n) ! REAL, INTENT (IN OUT) :: x, y(n) ! REAL, INTENT (OUT) :: hdid, hnext ! INTERFACE ! SUBROUTINE Derivs (n, x,y,dxdy) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x, y(n) ! REAL, INTENT (OUT) :: dxdy(n) ! END SUBROUTINE Derivs ! END INTERFACE ! END SUBROUTINE Rk_adapt_step ! END INTERFACE EXTERNAL Rk_adapt_step, Derivs ! INTEGER, PARAMETER :: MAXSTP=10000 REAL, PARAMETER :: tiny = 1.0E-30 INTEGER :: i, nstp REAL :: h, hdid, hnext, x, dydx(Nvar),y(Nvar), yscal(Nvar) ! x = x1 h = Sign (h1, x2-x1) nok = 0 nbad = 0 y(1:nvar) = ystart(1:nvar) ! DO nstp = 1, MAXSTP CALL Derivs (nvar, x, y, dydx) yscal (1:nvar) = ABS(y(1:nvar)) + ABS(h*dydx(1:nvar))+TINY IF ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x CALL Rk_adapt_step (y, dydx, nvar, x, h, eps, yscal, hdid, hnext, Derivs) IF (hdid == h) THEN nok = nok + 1 ELSE nbad = nbad + 1 END IF IF ((x-x2)*(x2-x1) >= 0.0) THEN ystart(1:nvar) = y(1:nvar) RETURN END IF IF (ABS(hnext) < hmin) STOP 'stepsize smaller than minimum in ode_int' h = hnext END DO STOP 'too many steps in ode_int' RETURN END SUBROUTINE Ode_int ! ! !******************************************************************************** SUBROUTINE Rk_adapt_step (y, dydx, n, x, htry, eps, yscal, hdid, hnext, Derivs) !******************************************************************************** ! ! Take one adaptive Runge-Kutta step to ensure accuracy by monitoring truncation ! error. Advance N ODEs dydx=Derivs(x,y) from point X by step H. YSCAL is the ! vectory to scale accuracy against. HTRY is initial step, HDID is the actual ! step taken, HNEXT is suggested next step. ! IMPLICIT NONE INTEGER, INTENT (IN) :: n REAL, INTENT (IN) :: htry, eps, yscal(n), dydx(n) REAL, INTENT (IN OUT) :: x, y(n) REAL, INTENT (OUT) :: hdid, hnext ! INTERFACE ! SUBROUTINE Derivs (n, x,y,dxdy) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x, y(n) ! REAL, INTENT (OUT) :: dxdy(n) ! END SUBROUTINE Derivs ! END INTERFACE EXTERNAL Derivs ! REAL :: ytemp(N), yerr(n), errmax, htemp, xnew, h REAL, PARAMETER :: SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4 ! h = htry 1 CALL Rk_cash_karp (y, dydx, n, x, h, ytemp, yerr, Derivs) errmax = MAXVAL (ABS(yerr(:)/yscal(:)))/ eps IF (errmax > 1.0) THEN htemp = SAFETY*h*(errmax**PSHRNK) h = SIGN(MAX(ABS(htemp),0.1*ABS(h)),h) xnew=x+h IF (xnew == x) STOP 'stepsize underflow in rk_cash_karp' GO TO 1 ELSE IF (errmax > ERRCON) THEN hnext=SAFETY*h*(errmax**PGROW) ELSE hnext=5.*h END IF hdid=h x=x+h y(:)=ytemp(:) RETURN END IF RETURN END SUBROUTINE Rk_adapt_step ! ! ! !************************************************************************************ SUBROUTINE Rk_cash_karp (y, dydx, n, x, h, yout, yerr, derivs) !************************************************************************************ ! ! Takes one Cash-Karp Runge-Kutta time step to advance system of N ! ODEs dy/dx = Derivs (y,x) over step H assuming Y and DYDX are ! known at point X. Result is in YOUT, and YERR has an estimate of ! the local truncation error. After NR. ! IMPLICIT NONE INTEGER, INTENT (IN) :: n REAL, INTENT (IN) :: h, x, dydx(n), y(n) REAL, INTENT (OUT) :: yerr(n), yout(n) ! INTERFACE ! SUBROUTINE Derivs (n, x,y,dydx) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x, y(n) ! REAL, INTENT (OUT) :: dydx(n) ! END SUBROUTINE Derivs ! END INTERFACE EXTERNAL Derivs ! REAL :: ak2(N),ak3(N),ak4(N),ak5(N),ak6(N), ytemp(N) REAL, PARAMETER :: & A2=0.2, A3=0.3, A4=0.6, A5=1.0, A6=0.875, & B21=0.2, & B31=3.0/40.0, B32=9./40., & B41=.3,B42=-.9,B43=1.2, & B51=-11./54.,B52=2.5, B53=-70./27.,B54=35./27., & B61=1631./55296.,B62=175./512., B63=575./13824.,B64=44275./110592.,B65=253./4096.,& C1=37./378., C3=250./621.,C4=125./594.,C6=512./1771., & DC1=C1-2825./27648., DC3=C3-18575./48384.,DC4=C4-13525./55296.,DC5=-277./14336., DC6=C6-.25 ! ytemp(:) = y(:) + B21*h*dydx(:) CALL Derivs (n, x+A2*h, ytemp, ak2) ytemp(:) = y(:) + h*(B31*dydx(:)+B32*ak2(:)) CALL Derivs (n, x+A3*h,ytemp,ak3) ytemp(:) = y(:) + h*(B41*dydx(:)+B42*ak2(:)+B43*ak3(:)) CALL Derivs (n, x+A4*h,ytemp,ak4) ytemp(:) = y(:) + h*(B51*dydx(:)+B52*ak2(:)+B53*ak3(:)+B54*ak4(:)) CALL Derivs (n, x+A5*h,ytemp,ak5) ytemp(:) = y(:) + h*(B61*dydx(:)+B62*ak2(:)+B63*ak3(:)+B64*ak4(:)+B65*ak5(:)) CALL Derivs (n, x+A6*h,ytemp,ak6) yout(:)=y(:) + h*(C1*dydx(:)+C3*ak3(:)+C4*ak4(:)+C6*ak6(:)) ! yerr(:)=h*(DC1*dydx(:)+DC3*ak3(:)+DC4*ak4(:)+DC5*ak5(:)+DC6*ak6(:)) RETURN END SUBROUTINE Rk_cash_karp ! !----------------------------------------------------------------------------------- ! ! !************************************************************** SUBROUTINE Rk4 (y,dydx,n,x,h,yout,derivs) !************************************************************** ! ! Numerical Recipes subroutine to carry out one Runge-Kutta step ! on a set of n differential equations. Given values for the ! variables y(1:n) and their derivatives dydx(1:n) known at x, ! it uses 4-order RK to advance the solution over interval h and ! return th incremented variables as yout(1:n). Subroutine derivs ! is supplied by the user to compute dydx at x. ! IMPLICIT NONE INTEGER, INTENT (IN) :: n REAL, INTENT (IN) :: h, x, dydx(n), y(n) REAL, INTENT (OUT) :: yout(n) ! INTERFACE ! SUBROUTINE Derivs (n, x, y, dydx) ! IMPLICIT NONE ! INTEGER, INTENT (IN) :: n ! REAL, INTENT (IN) :: x ! REAL, INTENT (OUT) :: y(n), dydx(n) ! END SUBROUTINE Derivs ! END INTERFACE EXTERNAL Derivs INTEGER, PARAMETER :: NMAX=50 INTEGER :: i REAL :: h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX) hh=h*0.5 h6=h/6.0 xh=x+hh yt(:)=y(:)+hh*dydx(:) CALL Derivs (n, xh,yt,dyt) yt(:) = y(:)+hh*dyt(:) CALL Derivs (n, xh,yt,dym) DO i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) END DO CALL Derivs(n, x+h,yt,dyt) yout(:)=y(:)+h6*(dydx(:)+dyt(:)+2.*dym(:)) RETURN END SUBROUTINE Rk4 ! ! ! ! ! PROGRAM RCM_conductances_IRI90 USE Rcm_mod_subs, ONLY : isize, jsize, colat, aloct, sini, & qtplam, qtped, qthall, ss, pwe, pwn, hwe, hwn, sw,& read_grid, write_qtcond USE FIELD_LINE USE Conditions USE Constants IMPLICIT NONE ! ! integer, parameter :: isize=155,jsize=99, lun=11 ! real, dimension (isize,jsize) :: colat, aloct, sini, qtplam, & ! qtped, qthall, pwe, pwn, hwe, hwn,& ! alpha, beta, vcorot, bir integer :: isize_pr, jsize_pr, jwrap_pr ! real :: dlam, dpsi, ri, re ! real, dimension (jsize) :: ss, sw character (len=80) :: form_string, form_length ! !________________________________________________________________________ ! ! Code to compute quiet-time conductances and winds for use with RCM. ! Computation of conductances and winds is based on IRI,MSIS, and HWM ! empirical models. ! Parts of code were taken from Pete Riley's programs written at ! Rice University. ! Stanislav Sazykin, USU, last modif on 9/24/97 !__________________________________________________________________________ ! ! INTEGER i,j, istat REAL int_sig_p, int_sig_h, int_pwn,int_pwe,int_hwn,int_hwe, & int_ss,int_sw ! WRITE (*,'(A)',advance='no') 'daynr=' READ (*,*) daynr ! WRITE (*,'(A)',advance='no') 'f107(PREV. DAY)=' READ (*,*) f107 ! WRITE (*,'(A)',advance='no') 'f107(81-DAY RUN.AVG)=' READ (*,*) f107a ! WRITE (*,'(A)',advance='no') 'f107(12-MONTHS RUN.AVG)=' READ (*,*) f107_12m ! ap = 0. WRITE (*,'(A)',advance='no') 'Ap(DAILY)=' READ (*,*) ap(1) ! CALL Read_grid () ! get RCM's grid ! ! H_top = 1000.0 H_bottom = 85.0 ! limits of integration ! DO i = isize, 1, -1 DO j = 1, jsize ! ! phi = aloct(i,j) R_zero = (R_e + h_bottom) / (SIN(colat(i,j))**2) theta = colat(i,j) height_low = h_bottom theta_low = theta IF ((R_zero - R_e) > H_top) THEN height_high = h_top theta_high = ASIN(SQRT((R_e+H_top)/R_zero)) ELSE height_high = R_zero-R_e theta_high = 0.5 * pi END IF ! IF (sini(i,j) > 0.77 ) THEN CALL Integrate_height (height_low, height_high,& int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe) ELSE CALL Integrate_theta (theta_low, theta_high,& int_sig_p, int_sig_h, & int_pwn, int_pwe, int_hwn, int_hwe) END IF ! WRITE (*,'(A3,I3.3,A3,I3.3,TR5,ES12.5,TR5,ES12.5)') 'I=',i,' J=',j,int_sig_p, int_sig_h ! qtplam (i,j) = int_sig_p / sini(i,j) qthall (i,j) = int_sig_h qtped (i,j) = int_sig_p * sini(i,j) ! pwn(I,J) = int_pwn/int_sig_p/sini(i,j) pwe(I,J) = int_pwe/int_sig_p hwn(I,J) = int_hwn/int_sig_h/sini(i,j) hwe(I,J) = int_hwe/int_sig_h END DO END DO ! ! Now equatorial quantities SS and SW: ! DO j = 1, jsize phi = aloct(isize,j) call integrate_theta_eq(colat(isize,j),pi/2.,int_ss,int_sw) ! CALL Integrate_theta_eq (colat(isize,j), pi/2.-0.0873, & ! int_ss,int_sw) WRITE (*,'(A,I3.3,2(A,F6.2))') 'Equatorial band, j=',j, ' SS=', int_ss, ' SW=',int_sw ss(j) = int_ss sw(j) = int_sw END DO ! CALL Write_qtcond () ! STOP END PROGRAM RCM_conductances_IRI90 ! ! ! SUBROUTINE Sigma_tensor (theta, phi, height, & daynr, f107, f107a, f107_12m, ap, & sig_p, sig_h, v_perp_n, v_e) IMPLICIT NONE INTEGER, INTENT (IN) :: daynr REAL, INTENT (IN) :: f107, f107a, f107_12m, ap(7), theta, phi, height REAL, INTENT (OUT) :: sig_p, sig_h, v_perp_n, v_e ! !******************************************************************* ! ! Compute Pedersen and Hall conductivities in the ionosphere ! using IRI,MSIS models, also two wind components: zonal (+east) ! and perp north (+ north and up) from HWM. ! Code modified, originally taken from P. Riley's IMP program. ! The coordinate system is geomagnetic (centered aligned dipole). ! LAST MODIFICATION: Stanislav, July 16,1997, 7/22/97, 8/15/97, 8/18/97 ! real b_tmp,gam1_e,gam1_no_plus,gam1_o2_plus,gam1_o_plus REAL omega_o2_plus,omega_o_plus,omega_e,omega_no_plus REAL IRI_TE,IRI_TI,IRI_O_PLUS,IRI_NO_PLUS,IRI_O2_PLUS,NE REAL MSIS_TN,MSIS_N2,MSIS_O2,MSIS_O REAL T_R ! REAL :: mlat_deg, mlong_deg, mlt_hours, glat_deg, glong_deg,stl, utsec, d(8), t(2) REAL :: doutf(11,50), oarr(35), wind(2) REAL :: msis_switches (25), hwm_switches (25) LOGICAL :: jf (12) ! REAL, PARAMETER :: pi = 3.141592654, & besu = 31000.0, & R_e = 6371.0, & Q = 1.6021892E-19, & M_E = 9.1094E-31, M_H = 1.6726E-27, & M_NO_PLUS=30.*M_H, & M_O_PLUS=16.*M_H, & M_O2_PLUS=32.*M_H ! ! msis_switches = 1.0 msis_switches (10:13) = 0.0 ! all UT-related effects off hwm_switches = 1.0 hwm_switches (10:13) = 0.0 ! all UT-related effects off ! mlat_deg = (pi/2.-theta)*180.0/pi ! IF (phi <= pi) THEN mlt_hours = phi*12.0/pi+12. ELSE mlt_hours = phi*12.0/pi-12. END IF mlong_deg = 87.0+90. stL = mlt_hours utsec = (mlt_hours - glong_deg/15.0)*3600.0 ! CALL ggm_iri (1, glong_deg, glat_deg, mlong_deg, mlat_deg) ! ! Call MSIS-90 to get neutral parameters: ! CALL Tselec (msis_switches) CALL Meter6 (.TRUE.) CALL Gtd6 (daynr, utsec, height, glat_deg, glong_deg, stL, & f107a, f107, ap, 48, d, t) msis_tn = t(2) msis_o = d(2) msis_n2 = d(3) msis_o2 = d(4) ! ! ! Call HWM93 to get two horizontal wind components: ! CALL Tselec (hwm_switches) CALL Gws5 (daynr, utsec, height, glat_deg, glong_deg, stL,& f107a, f107, ap, wind) v_e = wind (2) v_perp_n = wind (1) * & ! correct for SIN(I) 2.0*COS(theta)/SQRT(1.0+3.0*(COS(theta))**2) ! ! ! Call IRI-91 to get plasma parameters: ! Since it's impossible to turn off UT and LONG effects in IRI, ! let us choose the longitude of Jicamarca: ! jf (1:11) = .TRUE. jf (12) = .FALSE. IF (height < 80.0) then print*,height pause STOP 'TOO LOW FOR IRI, H<80 KM' END IF CALL Iris12 (jf, 1, mlat_deg, mlong_deg, -f107_12m, -daynr, & stL, height, height, 1.0, doutf, oarr) ne = doutf(1,1) iri_ti = doutf(3,1) iri_te = doutf(4,1) IF (height <= 100.0) THEN ! adjustment by Pete Riley from his code iri_o_plus = 0. iri_o2_plus = 0.5*ne iri_no_plus = 0.5*ne ELSE IRI_O_PLUS = DOUTF(5,1)*NE/100. IRI_O2_PLUS = DOUTF(8,1)*NE/100. IRI_NO_PLUS = DOUTF(9,1)*NE/100. END IF ! IF (height <= 120.0) THEN ! The temperature cross check/modifications: IRI_TE = MSIS_TN IRI_TI = MSIS_TN END IF T_R = (MSIS_TN + IRI_TI)/2.0 ! ! ! Calculate the electron and ion gyrofrequences ! b_tmp = (BESU*1.E-9)*SQRT(1.+3.*(COS(theta))**2)/ (((R_e+height)/R_e)**3) OMEGA_E = Q*B_TMP/M_E OMEGA_NO_PLUS = Q*B_TMP/M_NO_PLUS OMEGA_O2_PLUS = Q*B_TMP/M_O2_PLUS OMEGA_O_PLUS = Q*B_TMP/M_O_PLUS ! ! ! Calculation of the conductivities using the formulae in Rees(89) ! "Physics and chemistry of the upper atmosphere", Cambridge, pp.202-203. ! ! 1. The parallel conductivity is ! SIG_0 =(Q**2)*NE/ ! $ (M_E*NU_E(NU_EI(ME,IRI_TE), ! $ NU_EN(MSIS_N2,MSIS_O2,MSIS_O,IRI_TE)) ! $ + Q**2*( ! $ + IRI_NO_PLUS/(M_NO_PLUS* ! $ NU_NO_PLUS(MSIS_N2,MSIS_O2,MSIS_O)) ! $ + IRI_O_PLUS/(M_O_PLUS* ! $ NU_O_PLUS(MSIS_N2,MSIS_O2,MSIS_O,T_R)) ! $ ) ! ! 2. The Pedersen conductivity is ! gam1_e = NU_EN (msis_n2,msis_o2,msis_o,iri_te) * omega_e/ & (omega_e**2 + NU_EN(msis_n2,msis_o2,msis_o,iri_te)**2) gam1_no_plus = NU_NO_PLUS(msis_n2,msis_o2,msis_o)*omega_no_plus/ & (omega_no_plus**2+NU_NO_PLUS(msis_n2,msis_o2,msis_o)**2) gam1_o2_plus = NU_O2_PLUS(msis_n2,msis_o2,msis_o,t_r)* & omega_o2_plus/ & (omega_o2_plus**2+ & NU_O2_PLUS(msis_n2,msis_o2,msis_o,T_r)**2) gam1_o_plus = NU_O_PLUS(msis_n2,msis_o2,msis_o,T_r)* & omega_o_plus/ & (omega_o_plus**2+ & NU_O_PLUS(msis_n2,msis_o2,msis_o,T_r)**2) sig_p = (Q/B_TMP)*( & IRI_NO_PLUS*(gam1_e + gam1_no_plus)+ & IRI_O2_PLUS*(gam1_e + gam1_o2_plus)+ & IRI_O_PLUS *(gam1_e + gam1_o_plus ) & ) ! ! 3. The Hall conductivity is ! gam1_e = omega_e**2 / (omega_e**2+NU_EN(msis_n2,msis_o2,msis_o,iri_te)**2) gam1_no_plus = omega_no_plus**2 / & ((omega_no_plus**2)+NU_NO_PLUS(msis_n2,msis_o2,msis_o)**2) gam1_o2_plus = omega_o2_plus**2 / & (omega_o2_plus**2+NU_O2_PLUS(msis_n2,msis_o2,msis_o,t_r)**2) gam1_o_plus = omega_o_plus**2 / & (omega_o_plus**2+NU_O_PLUS(msis_n2,msis_o2,msis_o,t_r)**2) sig_h = (Q/B_TMP)*( & IRI_NO_PLUS*(gam1_e - gam1_no_plus)+ & IRI_O2_PLUS*(gam1_e - gam1_o2_plus)+ & IRI_O_PLUS* (gam1_e - gam1_o_plus) & ) ! RETURN ! ! CONTAINS ! ! REAL FUNCTION Nu_no_plus (den_n2, den_o2, den_o) IMPLICIT NONE REAL, INTENT (IN) :: den_n2, den_o2, den_o ! ! appendix of Schunk and Walker(73) ! Nu_no_plus = (4.34E-16)*den_n2 & + (4.28E-16)*den_o2 & + (2.44E-16)*den_o RETURN END FUNCTION Nu_no_plus ! ! ! REAL FUNCTION Nu_o2_plus (den_n2, den_o2, den_o, t_r) IMPLICIT NONE REAL, INTENT (IN) :: den_n2, den_o2, den_o, t_r ! ! appendix of Schunk and Walker(73) ! IF (T_R <= 800.) THEN Nu_o2_plus = (4.13E-16)*den_n2 & + (4.08E-16)*den_o2 & + (2.31E-16)*den_o ELSE Nu_o2_plus = (4.13E-16)*den_n2 & + (2.4E-19)*DEN_O2*SQRT(t_r)* & ((10.4-0.76*LOG10(t_r))**2) & + (2.31E-16)*den_o END IF RETURN END FUNCTION Nu_o2_plus ! ! ! REAL FUNCTION Nu_o_plus (den_n2, den_o2, den_o, t_r) IMPLICIT NONE REAL, INTENT (IN) :: den_n2, den_o2, den_o, t_r ! ! appendix of Schunk and Walker(73) ! Nu_o_plus = (6.82E-16)*den_n2 & + (6.66E-16)*den_o2 & + (4.0E-17)*den_o*SQRT(t_r) RETURN END FUNCTION Nu_o_plus ! ! ! REAL FUNCTION Nu_en (den_n2,den_o2,den_o,te) IMPLICIT NONE REAL, INTENT (IN) :: den_o2, den_n2, den_o,te ! ! The electron-neutral collision frequency is ! from Schunk and Nagy(1978) reported by Schlegel(84ish) ! Nu_en = (2.33E-17)*den_n2*(1.0-(1.21E-4)*te)*te & + (1.82E-16)*den_o2*(1.0+(3.6E-2)*SQRT(te))* & SQRT(te) & + (8.9E-17)*den_o*(1.0+(5.7E-4)*te)*SQRT(te) RETURN END FUNCTION Nu_en ! ! ! ! REAL FUNCTION NU_EI (NE,TE) ! IMPLICIT NONE ! REAL NE,TE ! the electron-ion collision frequency is from Hanson(1961) ! NU_EI = (1E-6)*(59.08 + 4.18*LOG10(TE**3/NE))* ! $ NE*TE**(-1.5) ! END ! ! REAL FUNCTION NU_E (NU_EI,NU_EN) ! IMPLICIT NONE ! REAL NU_EI,NU_EN ! The total electron collision frequency is ! NU_E = NU_EI + NU_EN ! END ! ! END SUBROUTINE Sigma_tensor