#include "dims.h" ! ! am_02/02: calculate current for both hemispheres subroutine nosocrrt use cons_module,only: pi implicit none ! **** ! **** calculate current for both hemisphere: ! **** [stencil*potential -RHS] = R**2 * J_mr / dt0dts / rcos0s ! **** #include "params.h" #include "dynphi.h" ! phim #include "transmag.h" #include "consdyn.h" ! Global: common/nscoeff/ nsc(imx0,jmaxm,10),nscrrt(imx0,jmaxm) ! imx0=imaxmp real :: nsc,nscrrt,tout(imaxmp,jmaxm,-2:zkmxp) ! Local: real :: vtmp,r0sq,fac,facmax,lat,lmin,lmax,pol integer :: j,jj,jjj,i,k,n,jmod,jmin,jmax ! External: real,external :: sddot ! in util.F ! nscrrt(:,:) = 0.0 ! ! for i=1 ! ! i = 1 j = 49 ! equator set J_mr = 0 nscrrt(i,j) = 0.0 ! do j = 2,jmaxmh-1 ! 2,48 no poles jj = jmaxmh+j-1 ! 50,96 jjj = jmaxmh-j+1 ! 48,2 ! contribution of stencil ! northern hemisphere ! nscrrt(i,jj) = nsc(i,jj,1)*phim(i+1,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,2)*phim(i+1,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,3)*phim(i ,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,4)*phim(imx0-1,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,5)*phim(imx0-1,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,6)*phim(imx0-1,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,7)*phim(i ,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,8)*phim(i+1,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,9)*phim(i ,jj) nscrrt(i,jj) = nscrrt(i,jj) - nsc(i,jj,10) ! ! southern hemisphere ! nscrrt(i,jjj) = nsc(i,jjj,1)*phim(i+1,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,2)*phim(i+1,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,3)*phim(i ,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,4)*phim(imx0-1,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,5)*phim(imx0-1,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,6)*phim(imx0-1,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,7)*phim(i ,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,8)*phim(i+1,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,9)*phim(i ,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) -nsc(i,jjj,10) ! enddo ! ! do i = 2,imx0-1 ! 2,80 ! j = 49 ! equator nscrrt(i,j) = 0.0 ! do j = 2,jmaxmh-1 ! 2,48 no poles jj = jmaxmh+j-1 ! 50,96 jjj = jmaxmh-j+1 ! 48,2 ! contribution of stencil ! northern hemisphere ! nscrrt(i,jj) = nsc(i,jj,1)*phim(i+1,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,2)*phim(i+1,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,3)*phim(i ,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,4)*phim(i-1,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,5)*phim(i-1,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,6)*phim(i-1,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,7)*phim(i ,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,8)*phim(i+1,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,9)*phim(i ,jj) nscrrt(i,jj) = nscrrt(i,jj) - nsc(i,jj,10) ! southern hemisphere ! nscrrt(i,jjj) = nsc(i,jjj,1)*phim(i+1,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,2)*phim(i+1,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,3)*phim(i ,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,4)*phim(i-1,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,5)*phim(i-1,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,6)*phim(i-1,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,7)*phim(i ,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,8)*phim(i+1,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) + nsc(i,jjj,9)*phim(i ,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) - nsc(i,jjj,10) enddo ! enddo ! ! for i=81 ! ! i = 81 j = 49 ! equator nscrrt(i,j) = 0.0 ! do j = 2,jmaxmh-1 ! 2,48 no poles jj = jmaxmh+j-1 ! 50,96 jjj = jmaxmh-j+1 ! 48,2 ! contribution of stencil ! northern hemisphere ! nscrrt(i,jj) = nsc(i,jj,1)*phim(2,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,2)*phim(2,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,3)*phim(i ,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,4)*phim(i-1,jj+1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,5)*phim(i-1,jj) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,6)*phim(i-1,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,7)*phim(i ,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,8)*phim(2,jj-1) nscrrt(i,jj) = nscrrt(i,jj) + nsc(i,jj,9)*phim(i ,jj) nscrrt(i,jj) = nscrrt(i,jj) - nsc(i,jj,10) ! southern hemisphere ! nscrrt(i,jjj) = nsc(i,jjj,1)*phim(2,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,2)*phim(2,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,3)*phim(i ,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,4)*phim(i-1,jjj-1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,5)*phim(i-1,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,6)*phim(i-1,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,7)*phim(i ,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,8)*phim(2,jjj+1) nscrrt(i,jjj) = nscrrt(i,jjj) +nsc(i,jjj,9)*phim(i ,jjj) nscrrt(i,jjj) = nscrrt(i,jjj) -nsc(i,jjj,10) ! enddo ! poles nscrrt(1,1) = sddot(imaxm,unit,nscrrt(1,2))/float(imaxm) nscrrt(1,jmaxm)= sddot(imaxm,unit,nscrrt(1,jmaxm-1))/float(imaxm) nscrrt(:,1) = nscrrt(1,1) ! extend in longitude nscrrt(:,jmaxm) = nscrrt(1,jmaxm) do j = 1,jmaxm lat = ylatm(j) if(lat.le.-pi/18.) then lmin = lat jmin = j elseif(lat.gt.pi/18.) then lmax = lat jmax = j goto 300 endif enddo 300 r0sq = r0*r0*1.e-4 ! in meter facmax = dt0dts(jmax)/r0sq*rcos0s(jmax) do j= 1,jmaxm fac = dt0dts(j)/r0sq*rcos0s(j) lat = ylatm(j) do i = 1,imaxmp ! linear interpolation of J_mr between -10 and 10 degrees if(j.gt.jmin.and.j.lt.jmax) then pol = nscrrt(i,jmin)-nscrrt(i,jmax)*facmax pol = pol / (lmin-lmax) nscrrt(i,j)= nscrrt(i,jmin) + pol*(lat-lmin) else nscrrt(i,j) = nscrrt(i,j)*fac endif enddo ! endo i-loop ! simple smoothing in longitudinal direction (weighted average) ! (nonsmooth values due to interpolation) do i = 1,imaxmp if (i.eq.1) then tout(i,j,1)= (nscrrt(i+1,j)+ | 3.*nscrrt(i,j)+ nscrrt(imaxmp-1,j))/5.0 elseif (i.eq.imaxmp) then tout(i,j,1)= (nscrrt(2,j)+ | 3.*nscrrt(i,j)+nscrrt(i-1,j))/5.0 else tout(i,j,1)= (nscrrt(i+1,j)+ | 3.*nscrrt(i,j)+nscrrt(i-1,j))/5.0 endif nscrrt(i,j) = tout(i,j,1) tout(i,j,:) = nscrrt(i,j) ! copy for secondary history field enddo ! endo i-loop call addfsech('NSCRRT','J_mr (magnetic)','A/m^2 ', | tout(:,j,:),imaxmp,kmaxp+3,kmaxp+3,j) enddo ! enddo j-loop return end