#include "dims.h" ! ! am_02/02: include the changes from Art Richmond done in tgcm12 ! SUBROUTINE INTGRLS(JM) implicit none C **** C **** PERFORMS APPROXIMATE FIELD LINE QUADRATURES AT ALL C **** GEOMAGNETIC GRID POINTS C **** #include "params.h" #include "consdyn.h" #include "coefm.h" #include "transmag.h" ! ! Args: integer,intent(in) :: jm ! ! Local: real :: z0 integer :: i,k,ii,kk real :: | SINDM(IMAXMP),COSDM(IMAXMP),RAM(IMAXMP),AAM(IMAXMP), | COSIAM(IMAXMP),CSTHDAM(IMAXMP),RTADRAM(IMAXMP), | RRM(IMAXMP,-2:KMAXP),SINIDM(IMAXMP,-2:KMAXP), | COSIDM(IMAXMP,-2:KMAXP),COSTHDM(IMAXMP,-2:KMAXP), | RTRAMRM(IMAXMP,-2:KMAXP) real :: TINT1(JMAXM),TINT2(JMAXM),TINT3(JMAXM) ! real :: sinlm,clm2,absinim,ra,sqomrra,sqrra,afac,htfac real :: rtrmn,rora,del,omdel,sig1,sig2,ue1,ue2 ! quantities needed for field-line integration real :: htfunc(imaxmp,-2:kmax),htfunc2(imaxmp,-2:kmax) ! ! Global: real :: SIGMA1E,SIGMA2E,ADOTVE common/equat/ SIGMA1E(IMAXMP,-2:KMAX),SIGMA2E(IMAXMP,-2:KMAX), 5 ADOTVE(IMAXMP,-2:KMAX,2) ! ! z0 = h0 ! evaluate 2-D fields first ! sinlm = sin(ylatm(jm)) ! sin(lam_m) clm2 = 1. - sinlm*sinlm ! cos^2(lam_m) absinim = abs(sinlm)/sqrt(1.-.75*clm2) ! | sin I_m | ra = R0/clm2 ! (R_E + H_0)/cos^2(lam_m) sqomrra = sqrt(1.-R0/ra) ! sqrt(1/ R_0/R_A) = sin(lam_m) sqrra = sqrt(R0/ra) ! sqrt(R_0/R_A) afac = 2.*sqrt(ra-R0) ! 2*sqrt(R_A-R_0) used for field line integration htfac = sqrt(ra-.75*R0) ! sqrt(R_A -3/4*R_0) ds = (A dh) / (2*sqrt(h_A-h)) ! DO 7 I = 1,IMAXM AAM(I) = afac/abs(SINIAM(I)) ! 2*sqrt( h_A - h_0 )/ |sin I_m | 7 CONTINUE ! w.r to reference height A(h_R) ! evaluate 3-D fields DO 8 K = -2,KMAXP ! height varying fields DO 8 I = 1,IMAXM ! ! RR = R0+Z-Z0 RRM(I,K) = R0+ZM(I,K)-Z0 ! radius of magnetic point ! ! RTRAMR = RA-R IF +IVE, ZERO OTHERWISE RTRAMRM(I,K) = AMAX1(0.,ra-RRM(I,K)) ! h_A - h (not negative) ! ! RTRAMR = SQRT(RA-R) RTRAMRM(I,K) = SQRT(RTRAMRM(I,K)) ! sqrt(h_A - h) 8 CONTINUE ! interpolation to the middle of horizontal slices DO 6 K = -2,KMAX ! calculate height varying factor and interpolate DO 6 I = 1,IMAXM ! to middle of horizontal slices RRM(I,K) = .5*(RRM(I,K)+RRM(I,K+1)) RTRAMRM(I,K) = RTRAMRM(I,K)-RTRAMRM(I,K+1) ! ds for line integration ! ! htfunc = factor by which to multiply AAM(I)*d(sqrt(ra-r)) = ds htfunc(I,K) = sqrt(ra-.75*RRM(I,K))/htfac htfunc2(I,K) = htfunc(I,K)**2 6 CONTINUE DO 9 I = 1,IMAXM ! clear temporary arrays ZIGM11(I,JM)= 0. ZIGM22(I,JM)= 0. ZIGM2(I,JM) = 0. ZIGMC(I,JM) = 0. RIM(I,JM,1) = 0. RIM(I,JM,2) = 0. 9 CONTINUE ! compute integrals ! DO 10 K = -2,KMAX ! loop over height-level DO 10 I = 1,IMAXM ! loop over longitude rora = amin1(1.,RRM(I,K)/ra) ! (R_E+h)/(R_E+h_A) < 1 -> h_A > h del = (sqomrra*sqrt(rora)-sqrra*sqrt(1.-rora))/ ! (lam_m - lam) / lam_m = 1 abs(YLATM(JM)) ! sqrt(1-r_0/r_A)sqrt(r/r_A) - sqrt(r_0/r_A)sqrt(1-r/r_A) omdel = 1. - del ! ! Interpolate conductivities and winds in latitude along field line, assuming ! linear variation between foot of field line and magnetic equator. ! (For field lines other than those near the magnetic equator, del is nearly ! zero, so that the interpolated values are essentially the values for the ! latitude of the foot of the field line; inaccuracy of the assumption of ! linear variation is thus unimportant for these field lines.) ! sig1 = omdel*SIGMA1M(I,K) + del*SIGMA1E(I,K) ! interpolation sig2 = omdel*SIGMA2M(I,K) + del*SIGMA2E(I,K) ue1 = omdel*ADOTVM(I,K,1) + del*ADOTVE(I,K,1) ue2 = omdel*ADOTVM(I,K,2) + del*ADOTVE(I,K,2) ! ! height varying factors: ds = AAM*htfunc ! d_1^2/D = 1/htfunc * ADOTAM(I,1) ! d_2^2/D = htfunc * ADOTAM(I,2) ! d_1*d_2/D = 1 * A1A2M(I) ! ZIGM11(I,JM) = ZIGM11(I,JM) + sig1*RTRAMRM(I,K) ! int (sigma_p*d_1^2/D) ds : d_1^2/D s. below loop 11 ZIGM22(I,JM) = ZIGM22(I,JM) + sig1*RTRAMRM(I,K)*htfunc2(I,K) ! int (sigma_p*d_2^2/D) ds : d_2^2/D s. below loop 11 ! ! Following uses Richmond convention for sign of ZIGMC (= -ZIGMC of Ridley) ZIGMC(I,JM) = ZIGMC(I,JM) + sig1*RTRAMRM(I,K)*htfunc(I,K) !int (sigma_p*d_1*d_2/D) ds ZIGM2(I,JM) = ZIGM2(I,JM) + sig2*RTRAMRM(I,K)*htfunc(I,K) !int (sigma_h) ds RIM(I,JM,1) = RIM(I,JM,1) + (sig1*ADOTAM(I,1)*ue2 +(sig2 - ! int [sigma_p*d_1^2/D u_e2 + (sigma_h - (sigma_p*d_1*d_2)/D) u_e1 ] ds 2 sig1*A1A2M(I))*htfunc(I,K)*ue1)*RTRAMRM(I,K) RIM(I,JM,2) = RIM(I,JM,2) + (sig1*ADOTAM(I,2)*htfunc2(I,K)* ! int [(sigma_h + sigma_p*d_1*d_2/D) u_e2 - sigma_p*d_2^2/D u_e1 ] ds 1 ue1 -(sig2 + sig1*A1A2M(I))*htfunc(I,K)*ue2)* 2 RTRAMRM(I,K) 10 CONTINUE ! complete calculation and place result in /coefm/ ! ZIGM's are in S, RIM's are in A/m ! multiply by 1/100 to convert from [cm] to [m] ! ! At this point, ZIGM11 is int[sig_p*d_1^2/D] ds, i.e. Sigma_(phi phi)/abs(sin Im) ! ZIGM22 is int[sig_p*d_2^2/D] ds, i.e. Sigma_(lam lam)*abs(sin Im) ! ZIGMC is int[sig_p*d_1*d_2/D] ds, i.e. Sigma_c ! ZIGM2 is int[sigma_h] ds, i.e. Sigma_h ! ! RIM1 is int[(sigma_h-sigma_p*d_1*d_2/D)u_e1 + sigma_p*d_1^2/D u_e2] *A(h_r)*B_e3 ds, i.e. K_(m phi)^D/abs(sin Im) ! RIM2 is int[(sigma_h+sigma_p*d_1*d_2/D)u_e2 - sigma_p*d_2^2/D u_e1] *A(h_r)*B_e3 ds, K_(m lam)^D ( minus in northern hemisphere ! Change sign of RIM(2) in S. hemisphere to be compatible with transf ! At this point, RIM2 is +-K_(m lam)^D ! ! ! DO 11 I = 1,IMAXM ZIGM11(I,JM) = 1.e-2*ZIGM11(I,JM)*AAM(I)*ADOTAM(I,1) ZIGM22(I,JM) = 1.e-2*ZIGM22(I,JM)*AAM(I)*ADOTAM(I,2) ZIGMC(I,JM) = 1.e-2*ZIGMC(I,JM)*AAM(I)*A1A2M(I) ZIGM2(I,JM) = 1.e-2*ZIGM2(I,JM)*AAM(I) RIM(I,JM,1) = 1.e-2*RIM(I,JM,1)*AAM(I)*BMODM(I) RIM(I,JM,2) = 1.e-2*RIM(I,JM,2)*AAM(I)*BMODM(I) 11 CONTINUE ! RETURN END