#include "dims.h" 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 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) ! ! Diagnostics for plotting: real,dimension(imaxmp,kmaxp+3) :: | zigm11_plt, zigmc_plt, zigm2_plt, zigm22_plt, rim1_plt, rim2_plt C **** C **** C **** EVALUATE 2-D FIELDS FIRST C **** DATA Z0/90.E5/ DO 7 I = 1,IMAXM RAM(I) = R0/COS(YLATM(JM))**2 COSIAM(I) = SQRT(1.-SINIAM(I)**2) AAM(I) = 4.*RAM(I)*(RAM(I)-R0)/((R0*SINIAM(I)/COSIAM(I)) 1 **2+4.*(RAM(I)-R0)**2) CSTHDAM(I) = SQRT((RAM(I)-AAM(I)*(RAM(I)-R0))/RAM(I)) RTADRAM(I) = SQRT(AAM(I)/RAM(I)) 7 CONTINUE ! write(6,"('fieldline_integrals: latm=',i3)") jm ! write(6,"('latm=',i3,' ram=' ,/,(6e12.4))") jm,ram ! write(6,"('latm=',i3,' siniam=' ,/,(6e12.4))") jm,siniam ! write(6,"('latm=',i3,' cosiam=' ,/,(6e12.4))") jm,cosiam ! write(6,"('latm=',i3,' aam=' ,/,(6e12.4))") jm,aam ! write(6,"('latm=',i3,' csthdam=',/,(6e12.4))") jm,csthdam ! write(6,"('latm=',i3,' rtadram=',/,(6e12.4))") jm,rtadram C **** C **** 3-D FIELDS C **** DO 8 K = -2,KMAXP DO 8 I = 1,IMAXM C **** C **** RR = R0+Z-Z0 C **** RRM(I,K) = R0+ZM(I,K)-Z0 C **** C **** RTRAMR = RA-R IF +IVE, ZERO OTHERWISE C **** RTRAMRM(I,K) = RAM(I)-RRM(I,K) if (rtramrm(i,k) < 0.) rtramrm(i,k) = 0. C **** C **** COSTHD = SQRT((RA-AA*(RA-R))/RA) C **** COSTHDM(I,K) = SQRT((RAM(I)-AAM(I)*RTRAMRM(I,K))/ 1 RAM(I)) C **** C **** RTRAMR = SQRT(RA-R) C **** RTRAMRM(I,K) = SQRT(RTRAMRM(I,K)) C **** C **** SINID = TAN(ID) C **** SINIDM(I,K) = SIGN(2./(RRM(I,K)*RTADRAM(I))* 1 COSTHDM(I,K)*RTRAMRM(I,K),SINIAM(I)) C **** C **** COSID = COS(ID) C **** COSIDM(I,K) = 1./SQRT(1.+SINIDM(I,K)**2) C **** C **** SINID = SIN(ID) C **** SINIDM(I,K) = SINIDM(I,K)*COSIDM(I,K) 8 CONTINUE C **** C **** INTERPOLATE TO MIDDLE OF HORIZONTAL SLICES C **** DO 6 K = -2,KMAX DO 60 I = 1,IMAXM RRM(I,K) = .5*(RRM(I,K)+RRM(I,K+1)) SINIDM(I,K) = .5*(SINIDM(I,K)+SINIDM(I,K+1)) COSIDM(I,K) = .5*(COSIDM(I,K)+COSIDM(I,K+1)) COSTHDM(I,K) = .5*(COSTHDM(I,K)+COSTHDM(I,K+1)) RTRAMRM(I,K) = RTRAMRM(I,K)-RTRAMRM(I,K+1) 60 CONTINUE ! write(6,"('integrals: latm=',i3,' k=',i3,' rtramrm(:,k)=',/, ! | (6e12.4))") jm,k,rtramrm(:,k) ! write(6,"('latm=',i3,' k=',i3,' rrm(:,k)=',/, ! | (6e12.4))") jm,k,rrm(:,k) ! ! Periodic point: ! rrm (imaxmp,k) = rrm (1,k) ! sinidm (imaxmp,k) = sinidm (1,k) ! cosidm (imaxmp,k) = cosidm (1,k) ! costhdm(imaxmp,k) = costhdm(1,k) ! rtramrm(imaxmp,k) = rtramrm(1,k) 6 CONTINUE ! call addfsech('RRM' ,' ',' ',rrm ,imaxmp,kmaxp+3,kmaxp+3,jm) ! call addfsech('SINIDM' ,' ',' ',sinidm ,imaxmp,kmaxp+3,kmaxp+3,jm) ! call addfsech('COSIDM' ,' ',' ',cosidm ,imaxmp,kmaxp+3,kmaxp+3,jm) ! call addfsech('COSTHDM',' ',' ',costhdm,imaxmp,kmaxp+3,kmaxp+3,jm) ! call addfsech('RTRAMRM',' ',' ',rtramrm,imaxmp,kmaxp+3,kmaxp+3,jm) C **** C **** COMPUTE INTEGRALS C **** C **** CLEAR TEMPORARY SPACE C **** DO 9 I = 1,IMAXM TINT1(I) = 0. TINT2(I) = 0. 9 CONTINUE C **** C **** QUADRATURES FOR A1**2/P C **** DO 10 K = -2,KMAX DO 100 I = 1,IMAXM C **** C **** TINT1 = QUAD(SIG1*R*D(SQRT(RA-R))/COS(THETAD)**2) C **** TINT1(I) = TINT1(I)+SIGMA1M(I,K)*RRM(I,K) 1 /(COSTHDM(I,K)**2) 2 *RTRAMRM(I,K) C **** C **** TINT2 = QUAD(SIG1*R*COS(ID)**2*D(SQRT(RA-R))/ C **** COS(THETAD)**2) C **** TINT2(I) = TINT2(I)+SIGMA1M(I,K)*RRM(I,K) 1 *(COSIDM(I,K)/COSTHDM(I,K))**2* 2 RTRAMRM(I,K) 100 CONTINUE ! write(6,"('integrals: k=',i3,' j=',i3,' sigma_pedm=',/, ! | (6e12.5))") k,jm,sigma1m(:,k) ! write(6,"('integrals: k=',i3,' j=',i3,' rrm=',/, ! | (6e12.5))") k,jm,rrm(:,k) ! write(6,"('integrals: k=',i3,' j=',i3,' costhdm=',/, ! | (6e12.5))") k,jm,costhdm(:,k) ! write(6,"('integrals: k=',i3,' j=',i3,' cosidm=',/, ! | (6e12.5))") k,jm,cosidm(:,k) ! write(6,"('integrals: k=',i3,' j=',i3,' rtramrm=',/, ! | (6e12.5))") k,jm,rtramrm(:,k) 10 CONTINUE C **** C **** COMPLETE CALCULATION AND PLACE RESULT IN ZIGM11 IN /COEFM/ C **** DO 11 I = 1,IMAXM ZIGM11(I,JM) = BMODM(I)*CSTHDAM(I)/COSIAM(I)*RTADRAM(I)* 1 (ADOTAM(I,1)*TINT1(I)+(AZM(I,1)/ 2 COSIAM(I))**2*TINT2(I))/PM(I) ! write(6,"('integrals: i=',i3,' j=',i3,' tint1=',e12.4,' tint2=', ! | e12.4,' zigm11=',e12.4)") i,jm,tint1(i),tint2(i), ! | zigm11(i,jm) ! write(6,"('integrals: i=',i3,' j=',i3,' bmodm=',e12.4, ! | ' csthdam=',e12.4,' cosiam=',e12.4,' rtadram=',e12.4, ! | ' adotam=',e12.4,' azm=',e12.4,' pm=',e12.4,' zigm11=', ! | e12.4)") i,jm,bmodm(i),csthdam(i),cosiam(i),rtadram(i), ! | adotam(i,1),azm(i,1),pm(i),zigm11(i,jm) ! write(6,"('integrals: i=',i3,' j=',i3,' bmodm=',e12.5, ! | ' csthdam=',e12.5,' cosiam=',e12.5,' rtadram=',e12.5, ! | ' adotam=',e12.5,' azm=',e12.5,' pm=',e12.5)") ! | i,jm,bmodm(i),csthdam(i),cosiam(i),rtadram(i), ! | adotam(i,1),azm(i,1),pm(i) 11 CONTINUE ! write(6,"(/,'integrals: latm=',i3,' zigm11(:,latm)=',/,(6e12.5))") ! | jm,zigm11(:,jm) ! write(6,"('integrals: latm=',i3,' bmodm=',/,(6e12.5))") jm,bmodm ! write(6,"('integrals: latm=',i3,' csthdam=',/,(6e12.5))") ! | jm,csthdam ! write(6,"('integrals: latm=',i3,' cosiam=',/,(6e12.5))") ! | jm,cosiam ! write(6,"('integrals: latm=',i3,' rtadram=',/,(6e12.5))") ! | jm,rtadram ! write(6,"('integrals: latm=',i3,' adotam(:,1)=',/,(6e12.5))") ! | jm,adotam(:,1) ! write(6,"('integrals: latm=',i3,' azm(:,1)=',/,(6e12.5))") ! | jm,azm(:,1) ! write(6,"('integrals: latm=',i3,' tint1=',/,(6e12.5))") ! | jm,tint1 ! write(6,"('integrals: latm=',i3,' tint2=',/,(6e12.5))") ! | jm,tint2 ! write(6,"('integrals: latm=',i3,' pm=',/,(6e12.5))") ! | jm,pm C **** C **** CLEAR TEMPORARY SPACE C **** DO 12 I = 1,IMAXM TINT1(I) = 0. TINT2(I) = 0. 12 CONTINUE C **** C **** QUADRATURES FOR A2**2/P C **** DO 13 K = -2,KMAX DO 13 I = 1,IMAXM C **** C **** TINT1 = QUAD(SIG1*R*TAN(ID)**2*D(SQRT(RA-R))) C **** TINT1(I) = TINT1(I)+SIGMA1M(I,K)*RRM(I,K) 1 *(SINIDM(I,K)/COSIDM(I,K))**2* 2 RTRAMRM(I,K) C **** C **** TINT2 = QUAD(SIG1*R*D(SQRT(RA-R))) C **** TINT2(I) = TINT2(I)+SIGMA1M(I,K)*RRM(I,K) 1 *RTRAMRM(I,K) 13 CONTINUE C **** C **** COMPLETE CALCULATION AND PLACE RESULT IN ZIGM22 IN /COEFM/ C **** DO 14 I = 1,IMAXM ZIGM22(I,JM) = BMODM(I)*COSIAM(I)/CSTHDAM(I)*RTADRAM(I)* 1 (ADOTAM(I,2)/SINIAM(I)**2*TINT1(I)+ 2 (AZM(I,2)/COSIAM(I))**2*TINT2(I))/PM(I) 14 CONTINUE C **** C **** CLEAR TEMPORARY SPACE C **** DO 15 I = 1,IMAXM TINT1(I) = 0. TINT2(I) = 0. 15 CONTINUE C **** C **** QUADRATURES FOR A1.A2/P C **** DO 16 K = -2,KMAX DO 16 I = 1,IMAXM C **** C **** TINT1 = QUAD(SIG1*R*TAN(ID)*D(SQRT(RA-R))/ C **** COS(THETAD)) C **** TINT1(I) = TINT1(I)+SIGMA1M(I,K)*RRM(I,K) 1 *SINIDM(I,K)/(COSIDM(I,K)*COSTHDM(I,K))* 3 RTRAMRM(I,K) C **** C **** TINT2 = QUAD(SIG1*R*COS(ID)*D(SQRT(RA-R))/ C **** COS(THETAD)) C **** TINT2(I) = TINT2(I)+SIGMA1M(I,K)*RRM(I,K) 1 *COSIDM(I,K)/COSTHDM(I,K) 2 *RTRAMRM(I,K) 16 CONTINUE C **** C **** COMPLETE CALCULATION AND PLACE RESULT IN ZIGMC IN /COEFM/ C **** DO 17 I = 1,IMAXM ZIGMC(I,JM) = BMODM(I)*RTADRAM(I)*(A1A2M(I) 1 /SINIAM(I)*TINT1(I)+AZM(I,1) 2 *AZM(I,2)/(COSIAM(I)**2)*TINT2(I))/PM(I) 17 CONTINUE C **** C **** CLEAR TEMPORARY SPACE C **** DO 18 I = 1,IMAXM TINT1(I) = 0. 18 CONTINUE C **** C **** QUADRATURE FOR B.(A2 X A1) = 1./B C **** DO 19 K = -2,KMAX DO 19 I = 1,IMAXM C **** C **** TINT1 = QUAD(SIG1*R/COS(ID)*D(SQRT(RA-R))/ C **** COS(THETAD)) C **** TINT1(I) = TINT1(I)+SIGMA2M(I,K)*RRM(I,K) 1 /(COSIDM(I,K)*COSTHDM(I,K))*RTRAMRM(I,K) 19 CONTINUE C **** C **** COMPLETE CALCULATION AND PLACE RESULT IN ZIGM2 IN /COEFM/ C **** DO 20 I = 1,IMAXM ZIGM2(I,JM) = RTADRAM(I)*TINT1(I) 20 CONTINUE C **** C **** NOW FOR RHS TERMS C **** C **** CLEAR TEMPORAY SPACE C **** DO 21 I = 1,IMAXM TINT1(I) = 0. TINT2(I) = 0. TINT3(I) = 0. 21 CONTINUE C **** C **** A1.(SIGMA1*(V X B)+SIGMA2*V) C **** A2.(SIGMA1*(V X B)+SIGMA2*V) C **** C **** TINT1 = INTEGRAL(A1(X)*Q(X)) C **** DO 22 K = -2,KMAX DO 220 I = 1,IMAXM TINT1(I) = TINT1(I)+(ADOTVM(I,K,1)*SIGMA2M(I,K)+ 1 AZM(I,1)*WM(I,K)*SIGMA2M(I,K)*COSIDM(I,K)/COSIAM(I)- 2 AXVM(I,K,1)*SINIDM(I,K)*SIGMA1M(I,K)+ 3 BXAM(I,1)*WM(I,K)*COSIDM(I,K)*SIGMA1M(I,K)+ 4 VXBM(I,K)*AZM(I,1)*SIGMA1M(I,K)*COSIDM(I,K)**2/COSIAM(I))* 5 RRM(I,K)**2/COSTHDM(I,K)*RTRAMRM(I,K) TINT2(I) = TINT2(I)+(ADOTVM(I,K,2)*SIGMA2M(I,K)*SINIDM(I,K)/ 1 (COSIDM(I,K)*SINIAM(I))+ 2 AZM(I,2)*WM(I,K)*SIGMA2M(I,K)/COSIAM(I)- 3 AXVM(I,K,2)*SINIDM(I,K)**2*SIGMA1M(I,K)/ 4 (SINIAM(I)*COSIDM(I,K))+ 5 BXAM(I,2)*WM(I,K)*SIGMA1M(I,K)*SINIDM(I,K)/SINIAM(I)+ 6 VXBM(I,K)*AZM(I,2)*SIGMA1M(I,K)*COSIDM(I,K)/COSIAM(I))* 7 RRM(I,K)**2*RTRAMRM(I,K) 220 CONTINUE 22 CONTINUE C **** C **** EVALUATE RHS TERMS RIM(1), RIM(2) C **** DO 23 I = 1,IMAXM RIM(I,JM,1) = BMODM(I)**2/R0*RTADRAM(I)/COSIAM(I)*TINT1(I) RIM(I,JM,2) = BMODM(I)**2/R0*RTADRAM(I)/CSTHDAM(I)*TINT2(I) 23 CONTINUE C **** C **** SCALE COEFFICIENTS AND RHS C **** DO 24 I = 1,IMAXM ZIGM11(I,JM) = ZIGM11(I,JM)*RCOS0S(JM)/(DT0DTS(JM)* 1 ABS(SINIAM(I))) ZIGMC(I,JM) = ZIGMC(I,JM)/SINIAM(I) ZIGM22(I,JM) = ZIGM22(I,JM)*DT0DTS(JM)/(RCOS0S(JM)* 1 ABS(SINIAM(I))) RIM(I,JM,1) = RIM(I,JM,1)*SIGN(1.,SINIAM(I))/DT0DTS(JM) RIM(I,JM,2) = RIM(I,JM,2)/RCOS0S(JM) 24 CONTINUE ! do i=1,imaxm ! zigm11_plt(i,:) = zigm11(i,jm) ! zigmc_plt (i,:) = zigmc (i,jm) ! zigm2_plt (i,:) = zigm2 (i,jm) ! zigm22_plt(i,:) = zigm22(i,jm) ! rim1_plt (i,:) = rim (i,jm,1) ! rim2_plt (i,:) = rim (i,jm,2) ! enddo ! i=1,imaxm ! Periodic point: ! zigm11_plt(imaxmp,:) = zigm11_plt(1,:) ! zigmc_plt (imaxmp,:) = zigmc_plt (1,:) ! zigm2_plt (imaxmp,:) = zigm2_plt (1,:) ! zigm22_plt(imaxmp,:) = zigm22_plt(1,:) ! rim1_plt (imaxmp,:) = rim1_plt (1,:) ! rim2_plt (imaxmp,:) = rim2_plt (1,:) ! call addfsech('ZIGM11',' ',' ',zigm11_plt,imaxmp,kmaxp+3,kmaxp,jm) ! call addfsech('ZIGMC' ,' ',' ',zigmc_plt ,imaxmp,kmaxp+3,kmaxp,jm) ! call addfsech('ZIGM2' ,' ',' ',zigm2_plt ,imaxmp,kmaxp+3,kmaxp,jm) ! call addfsech('ZIGM22',' ',' ',zigm22_plt,imaxmp,kmaxp+3,kmaxp,jm) ! call addfsech('RIM1' ,' ',' ',rim1_plt ,imaxmp,kmaxp+3,kmaxp,jm) ! call addfsech('RIM2' ,' ',' ',rim2_plt ,imaxmp,kmaxp+3,kmaxp,jm) RETURN END C