#include "dims.h" ! am_02/02: calculate current desnity J_e1, K_(q,phi) & K_(q,lam) subroutine nosocrdens implicit none C **** C **** C **** #include "params.h" ! imx0 #include "consdyn.h" ! re,r0,ylatm #include "dynphi.h" ! zzm(imaxmp,jmaxm,-2:zkmxp) #include "transmag.h" ! sigma1m(imaxmp,-2:kmax), sigma2m(imaxmp,-2:kmax), ! adotvm(imaxmp,-2:kmax,2),bmodm(imaxmp), ! adotam(imaxmp,2),a1a2m(imaxmp) ! zkmx = 28; kmax = zkmx ! Global: common/ns3d/ sigma1m3d(imaxmp,jmaxm,-2:zkmxp), | sigma2m3d(imaxmp,jmaxm,-2:zkmxp), | bmodm3d(imaxmp,jmaxm), | adotv1m3d(imaxmp,jmaxm,-2:zkmxp), | adotv2m3d(imaxmp,jmaxm,-2:zkmxp), | a1a2m3d(imaxmp,jmaxm), | adotam3d(imaxmp,jmaxm), | sinim3d(imaxmp,jmaxm), | ed13d(imaxmp,jmaxm,-2:zkmxp),ed23d(imaxmp,jmaxm,-2:zkmxp) real :: sigma1m3d,sigma2m3d,bmodm3d,adotv1m3d,adotv2m3d, | a1a2m3d,adotam3d,sinim3d,ed13d,ed23d common/nscoeff/ nsc(imx0,jmaxm,10),nscrrt(imx0,jmaxm) ! imx0=imaxmp real :: nsc,nscrrt ! Local: real :: je13d(imaxmp,jmaxm,-2:zkmxp),kqphi(imaxmp,jmaxm,-2:zkmxp) real :: kqphi3d(imaxmp,jmaxm),dkqphi(imaxmp,jmaxm) real :: kqlam(imaxmp,jmaxm) real :: fac,lamm,sinlamm,sinim,dh,act,actpk,facq,facsin real :: coslamq2,sinlamq,r0m,ed1h,ed2h,facqm,facqp real :: fsumn(jmaxm),afsumn(jmaxm),epsn,difflm real :: fsums(jmaxm),afsums(jmaxm),epss integer :: i,j,jj,k,l,ip1f,ip2f,ip3f ! ! calculate current density component Je1 (Richmond: Ionospheric ! Electrodynamics using magnetic apex coordinates pp.203 (eq 5.7)) ! at 1/2 level ! ar half levels: sig_ped,sig_hall, d_1**2/D, d1*d2/D, ue1, ue2, be3 ! at full levels: ed1, ed2 ! ! je1/d = (sig_ped * d_1**2/D * (ed1 + ue2*be3) + ! (sig_ped* d1*d2/D - sig_hall) * (ed2 - ue1*be3) ! je13d = je1/d ! do j = 1,jmaxm do k = -2,zkmx do i = 1,imaxmp fac = sigma1m3d(i,j,k)*adotam3d(i,j) ed1h = 0.5*(ed13d(i,j,k)+ed13d(i,j,k+1)) ! ed1 at half level je13d(i,j,k) = fac*(ed1h+ | adotv2m3d(i,j,k)*bmodm3d(i,j)) fac = sigma1m3d(i,j,k)*a1a2m3d(i,j) - sigma2m3d(i,j,k) ed2h = 0.5*(ed23d(i,j,k)+ed23d(i,j,k+1)) ! ed2 at half level je13d(i,j,k) = je13d(i,j,k) + fac*(ed2h- | adotv1m3d(i,j,k)*bmodm3d(i,j)) enddo je13d(imaxmp,j,k) = je13d(1,j,k) enddo call addfsech('JE13D',' J_e1/D (magnetic)','A/m^2? ', | je13d(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j) enddo ! ! calculate K_(q,phi) (Richmond: Ionospheric ! Electrodynamics using magnetic apex coordinates pp.208 (eq 7.4)) ! K_(q,phi) = int_(h_l)^(h_u) [( [R_0/R]^0.5 * ! * je1/D * sin(lam_q)/sin(lam_m) sin(I_m)/sin(I)*D] dh ! with F = D*sin(lam_m)/sin(lam_q)*sin(I)/sin(I_m)*[R/R_0]^3 ! ! do j = 1,jmaxm sinlamq = sin(ylatm(j)) ! sin(lam_q) if(j.eq.jmaxmh) sinlamq = sin(ylatm(j-1)) coslamq2 = 1. - sinlamq*sinlamq ! cos^2(lam_q) do i = 1,imaxmp ! at equator sin lam_q/sin I is set to the average otherwise quotient = 0 ! check this later 010611 if(j.eq.jmaxmh) then facsin = sinlamq/sinim3d(i,j-1) ! sin(lam_q)/sin(I) k = -2 fac = r0/(re + zzm(i,j,k)) ! R_0 / R lamm = acos(sqrt(coslamq2*fac)) ! cos^2(lam_m) = R_0/R*cos^2(lam_q) lamm = sign(lamm,ylatm(j-1)) ! lam_m sinlamm = sin(lamm) ! sin(lam_m) sinim = sinlamm/sqrt(.25+.75*sinlamm**2) ! sin(I_m) = sin(lam_m)/sqrt(1/4+3/4*sin^2(lam_m)) facq = sinim/sinlamm*facsin ! sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) fac = fac**0.5 ! sqrt(R_0 / R) act = fac*facq ! sqrt(R_0 / R)*sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) else facsin = sinlamq/sinim3d(i,j) ! sin(lam_q)/sin(I) k = -2 fac = r0/(re + zzm(i,j,k)) ! R_0 / R lamm = acos(sqrt(coslamq2*fac)) ! cos^2(lam_m) = R_0/R*cos^2(lam_q) lamm = sign(lamm,ylatm(j)) ! lam_m sinlamm = sin(lamm) ! sin(lam_m) sinim = sinlamm/sqrt(.25+.75*sinlamm**2) ! sin(I_m) = sin(lam_m)/sqrt(1/4+3/4*sin^2(lam_m)) facq = sinim/sinlamm*facsin ! sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) fac = fac**0.5 ! sqrt(R_0 / R) act = fac*facq ! sqrt(R_0 / R)*sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) endif kqphi3d(i,j) = 0.0 do k = -2,zkmx ! height integration from k=-2 to zkmx kqphi(i,j,k) = 0.0 ! (R_e + h)/R_0 : R_0 = R_e + h0 ; h0 = 9.e6 cm fac = r0/(re + zzm(i,j,k+1)) ! R_0 / R lamm = acos(sqrt(coslamq2*fac)) ! cos^2(lam_m) = R_0/R*cos^2(lam_q) lamm = sign(lamm,ylatm(j)) ! lam_m sinlamm = sin(lamm) ! sin(lam_m) sinim = sinlamm/sqrt(.25+.75*sinlamm**2) ! sin(I_m) = sin(lam_m)/sqrt(1/4+3/4*sin^2(lam_m)) facq = sinim/sinlamm*facsin ! sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) dh = zzm(i,j,k+1) - zzm(i,j,k) ! dh dh = dh*1.e-2 ! convertion [cm] to [m] fac = fac**0.5 ! sqrt(R_0 / R) actpk = fac*facq ! sqrt(R_0 / R)*sin(I_m)/sin(lam_m)*sin(lam_q)/sin(I) ! kqphi3d(i,j) = kqphi3d(i,j) + 0.5*(actpk+act)* ! integration value at 1/2 level | je13d(i,j,k)*dh act = actpk enddo enddo enddo ! save on secondary history do j= 1,jmaxm do k = -2,zkmx kqphi(:,j,k) = kqphi3d(:,j) enddo call addfsech('KQPHI',' kqphi (magnetic)',' ', | kqphi(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j) enddo ! ! calculate K_(q,lam) (Richmond: Ionospheric ! Electrodynamics using magnetic apex coordinates pp.208 (eq 7.5)) ! K_(q,lam) = -1/cos(lam_q) int_(-pi/2)^(lam_q) [J_mr*R*cos(lam_q) + ! d(K_(q,phi))/d(phi_q) ] d(lam_q) ! ! d(K_(q,phi))/d(phi_q) fac = 0.5/dlonm ! 1/(2*d lon_m) do j = 1,jmaxm do i = 2,imaxmp-1 dkqphi(i,j) = (kqphi3d(i+1,j)-kqphi3d(i-1,j))*fac ! central difference enddo ! (kqphi3d(i+1/2)-kqphi3d(i-1/2))/2 dkqphi(1,j) = (kqphi3d(2,j)-kqphi3d(imaxmp-1,j))*fac dkqphi(imaxmp,j) = dkqphi(1,j) enddo r0m = r0*1.e-2 ! do i = 1,imaxmp fsums(1) = 0.0 afsums(1) = 0.0 kqlam(i,1)= -dkqphi(i,1) ! at south pole K_(q,lam) = -d(K_(q,phi))/d(phi) act = nscrrt(i,1)*r0m*cos(ylatm(1))+dkqphi(i,1) ! [J_mr*R*cos(lam_q)+d K_(qphi)]_(j=1) do j = 2,jmaxmh difflm = abs(ylatm(j)-ylatm(j-1)) ! d | (lam_q(j)-lam_q(j-1))| actpk = nscrrt(i,j)*r0m*cos(ylatm(j))+dkqphi(i,j) ! [J_mr*R*cos(lam_q)+d K_(qphi)]_(j) act = act +actpk ! [integrand]_(j-1)+[integrand]_(j) act = -act/2.0*difflm ! -[integrand]_(j-1/2)] d lam_q fsums(j) = fsums(j-1) + act ! int_(-pi/2)^(lam_q)[integrand] d lam_q afsums(j)= afsums(j-1) + abs(act) ! int_(-pi/2)^(lam_q) | [integrand] | d lam_q act = actpk enddo ! ! integrate from the north pole to equator j=jmaxm kqlam(i,jmaxm)= dkqphi(i,jmaxm) ! at north pole K_(q,lam) = d(K_(q,phi) fsumn(j) = 0.0 afsumn(j) = 0.0 act = nscrrt(i,j)*r0m*cos(ylatm(j))+dkqphi(i,j) ! [J_mr*R*cos(lam_q)+d K_(qphi)]_(jmaxm) do j = jmaxm-1,jmaxmh,-1 difflm = abs(ylatm(j)-ylatm(j+1)) ! d | (lam_q(j)-lam_q(j-1))| actpk = nscrrt(i,j)*r0m*cos(ylatm(j))+dkqphi(i,j) ! [J_mr*R*cos(lam_q)+d K_(qphi)]_(j) act = act +actpk ! [integrand]_(j+1)+[integrand]_(j) act = act/2.0*difflm ! [integrand]_(j+1/2)] d lam_q fsumn(j) = fsumn(j+1) + act ! int_(pi/2)^(lam_q)[integrand] d lam_q afsumn(j)= afsumn(j+1) + abs(act) ! int_(pi/2)^(lam_q) | [integrand] | d lam_q act = actpk enddo ! correction to equal integration from north to equator with ! integration from south to equator epsn = 0.5*(fsums(jmaxmh)-fsumn(jmaxmh))/afsumn(jmaxmh) ! half of error to south weighted by absolute value epss = 0.5*(fsums(jmaxmh)-fsumn(jmaxmh))/afsums(jmaxmh) ! half of error to north weighted by absolute value do j = 2,jmaxmh ! correct and copy into kqlam kqlam(i,j)= 0.0 kqlam(i,j)= fsums(j) - epss*afsums(j) ! kqlam_cor = kqlam - err/2/abs(kqlam)_south*abs(kqlam)_(lam_q) kqlam(i,j)= kqlam(i,j)/cos(ylatm(j)) ! kqlam_cor /cos(lam_q) enddo do j = jmaxm-1,jmaxmh,-1 kqlam(i,j)= 0.0 kqlam(i,j)= fsumn(j) + epsn*afsumn(j) ! kqlam_cor = kqlam + err/2/abs(kqlam)_north*abs(kqlam)_(lam_q) kqlam(i,j)= kqlam(i,j)/cos(ylatm(j)) ! kqlam_cor /cos(lam_q) enddo enddo ! end of i-loop do j= 1,jmaxm do i = 1,imaxmp do k = -2,zkmx kqphi(i,j,k) = kqlam(i,j) enddo enddo call addfsech('KQLAM',' kqlam(magnetic)',' ', | kqphi(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j) enddo ! return end