subroutine fcmp(mm,nlont,cp,sp,f) use ModSize implicit none save integer, intent(in) :: mm, nlont real, dimension(nlont), intent(in) :: cp, sp real, dimension(-mm:mm,nlont), intent(out) :: f real, parameter :: sq2 = 1.4142135623 integer :: m, lon do lon=1,nlont f(0,lon) = 1. f(-1,lon) = sq2*cp(lon) f(1,lon) = sq2*sp(lon) enddo do m=2,mmx do lon=1,nlont f(-m,lon) = cp(lon)*f(1-m,lon) - sp(lon)*f(m-1,lon) f(m,lon) = cp(lon)*f(m-1,lon) + sp(lon)*f(1-m,lon) enddo enddo return end subroutine fcmp subroutine solver(a,m,n,mdim,ipvt,z) ! columns m+1 to n contain right-hand-sides. mdim is first dimension ! of a in calling program. m.le.mdim and n.le.second dimension of ! a in calling program. on output columns m+1 to n contain solutions. ! n.gt.m ! ! save real, dimension(mdim) :: z real, dimension(mdim,n) :: a integer, dimension(mdim) :: ipvt call sgeco(a,mdim,m,ipvt,rcond,z) do j=m+1,n call sgesl(a,mdim,m,ipvt,a(1,j),0) enddo end subroutine solver SUBROUTINE MATSOLV(NUM,NData,DATFMA,WT,XWT,OBS,OBS2, & ER,ASUM,AO,A,A2,GV,GH,DCD, & ENERR,EEERR,VARIATION,ERRORS) ! ! Set up matrix and solve equation for fit of some parameter ! NUM = 1-4 (1=electric potential, 2=ln(conductance), ! 3=electron ave energy, 4=electron energy flux) use ModAMIE use ModIndices use ModConstants implicit none integer, intent(in) :: num, NData real, dimension(0:lonmx,0:ithmx) :: variation, errors real, dimension(mcoef) :: ao, asum, a, a2, ps, v, & colsme, culsm, colsmn, z integer, dimension(mcoef) :: ipvt real, dimension(MAX_DATA) :: wt, obs, obs2, xwt, er, gv, gh, dcd real, dimension(MAX_DATA, mcoef) :: datfma real, dimension(0:lonmx,0:ithmx) :: eeerr, enerr character (len=8), dimension(4) :: nam real*8, dimension(MAX_DATA) :: er_double real*8, dimension(mcoef,ncoef) :: ax_double real*8 :: x integer :: nump, i, ith, j, k, m, mm, imax, lon, ix real :: fac, xnorm, txnorm, sum, sum1, sum2, typvalc real :: dfint, sth, fint, xe, xn nam(1) = 'ELEC POT' nam(2) = 'LN(COND)' nam(3) = 'E- AVKEV' nam(4) = 'E- EFLUX' nump = 1 if (num == 2) nump = 2 fac = 1. if (num == 1) fac = 1.e+3 ! set UP MATRIX k=0 do m=-mmx,mmx ax_double(ibm(m):iem(m),1:imx+nump) = 0. if (NData /= 0) then do k=ibm(m),iem(m) do j=1,NData x = wt(j)*datfma(j,k) ax_double(k,imx+1) = ax_double(k,imx+1) + x*obs(j) do i=1,imx ax_double(k,i) = ax_double(k,i) + x*datfma(j,i) enddo enddo enddo if (nump /= 1) then do k=ibm(m),iem(m) do j=1,NData x = wt(j)*datfma(j,k) ax_double(k,imx+2) = ax_double(k,imx+2) + x*obs2(j) enddo enddo endif endif enddo ! find mean square data values for assumed spectral shape, normalize ! spectrum, and add diagonal contribution to ax. xnorm = 1. if (NData /= 0) then sum = 0. sum1 = 0. sum2 = 0. er_double = 0. !\ ! electric potential: !/ if (num == 1) then do j=1,NData do i=1,imx er_double(j) = er_double(j) + & dble(datfma(j,i))*dble(datfma(j,i))*dble(cu(i,i)) enddo enddo do i=1,imx-1 do k=i+1,imx do j=1,NData er_double(j) = er_double(j) + & 2.*dble(datfma(j,i))*dble(datfma(j,k))*dble(cu(k,i)) enddo enddo enddo !\ ! particles and conductance: !/ else do j=1,NData do i=1,imx er_double(j) = er_double(j) + & dble(datfma(j,i))*dble(datfma(j,i))*dble(ccu(i,i)) enddo enddo endif er = er_double do j=1,NData sum = sum + wt(j)*obs(j)*obs(j) sum2 = sum2 + wt(j)*er_double(j) if (wt(j).ne.0.) sum1 = sum1 + xwt(j) enddo if (sum .ne.0.) sum = sum*sum / (sum+sum1) if (sum2.ne.0.) xnorm = sum/sum2 !------------------------------------------------------------------------ ! dramatic change in the code: ! if iflxgnorm is 1 then we use the fitted xnorms based on the current ! AE. If not, use the calculated xnorm from above. These two are ! fundementally different because the fitted xnorms do not change much ! with time, while the fitted ones change significantly from time to ! time. ! A. Ridley Sept. 26, 1999 !------------------------------------------------------------------------ if (ReplaceXNORM) then if (num == 1) xnorm = 3.54196e09*ae(iteration_number) + 3.63769e12 if (num == 2) xnorm = -2.33317e-06*ae(iteration_number) + 0.0138756 if (num == 3) xnorm = 1.47461e-06*ae(iteration_number) + 0.00194865 if (num == 4) xnorm = 1.34341e-06*ae(iteration_number) + 0.0180242 endif endif typvalc = sqrt(xnorm) if (num == 1) typvalc = typvalc / ri imax = imx + nump if (Calculate_Errors) imax = 2*imx + nump if (imax > ncoef) then write (6,"(1x,'stopped because imax =',i4)") imax stop endif if (Calculate_Errors) then do j=imx+1+nump,imax do i=1,imx ax_double(i,i+imx+nump) = 1. ax_double(i,j) = 0. enddo enddo do i=1,imx ax_double(i,i+imx+nump) = 1. enddo endif !\ ! Potential Only !/ if (num .eq. 1) then do i=1,imx ax_double(i,i) = ax_double(i,i) + cuid(i)/xnorm enddo do i=1,imx-1 do k=i+1,imx ax_double(i,k) = ax_double(i,k) + cu(i,k)/xnorm ax_double(k,i) = ax_double(k,i) + cu(i,k)/xnorm enddo enddo !\ ! Particles and Conductances !/ else do i=1,imx ax_double(i,i) = ax_double(i,i) + ccuid(i)/xnorm do k=i+1,imx ax_double(i,k)=ax_double(i,k)+ccu(i,k)/xnorm ax_double(k,i)=ax_double(k,i)+ccu(i,k)/xnorm enddo enddo endif do j=1,ncoef do i=1,mcoef ax(i,j) = ax_double(i,j) enddo enddo call solver(ax,imx,imax,mcoef,ipvt,z) ! now a will hold solution do i=1,imx a(i) = ax(i,imx+1) enddo ! only for electric potential solution if (num .eq. 1) then do i=1,imx asum(i) = a(i) + ao(i) enddo endif ! only for conductance solution (i.e. a2 holds Hall conductance) if (num .eq. 2) then do i=1,imx a2(i) = ax(i,imx+2) enddo endif if (Calculate_Errors) then ! make lower-left off-diagonal terms of ax hold sum of conjugate ! elements, to cut in half the number of summing operations. ! first imx columns of ax will hold c(u(hat)-u). (or cc) ax(1,1) = ax(1,imx+nump+1) do i=2,imx do j=1,i-1 ax(i,j) = ax(i,j+imx+nump) + ax(j,i+imx+nump) enddo ax(i,i) = ax(i,i+imx+nump) enddo ! test whether c(u(hat)-u) is physical k = imx + 1 3042 format(1x,'i=',i4,'j=',i4,'ax(i,j)=',e10.3,'ax(i,i)=',& e10.3,'ax(j,j)=',e10.3) do i=1,imx do j=1,i-1 if (ax(i,j+k)*ax(i,j+k).gt.ax(i,i+k)*ax(j,j+k)) then write(6,3042) i,j,ax(i,j+k),ax(i,i+k),ax(j,j+k) endif enddo enddo ! test symmetry of c(u(hat)-u) ! do i=10,imx,10 ! do j=10,i,10 ! write(6,'(1x,2i4,2e10.3)') i,j,ax(i,j+k),ax(j,i+k) ! enddo ! enddo endif ! compute calculated data values. ! gv is difference between observed and computed components. if (NData /= 0) then do j=1,NData gv(j) = -obs(j) enddo do j=1,NData do i=2,imx colsmn(i-1) = 0. do k=1,i-1 colsmn(k) = datfma(j,i) * datfma(j,k) * ax(i,k) + colsmn(k) enddo enddo dcd(j) = datfma(j,imx) * datfma(j,imx) * ax(imx,imx) do i=1,imx-1 dcd(j) = dcd(j) + colsmn(i) + datfma(j,i)*datfma(j,i)*ax(i,i) enddo enddo do i=1,imx do j=1,NData gv(j) = gv(j) + datfma(j,i)*a(i) enddo enddo if (nump == 2) then do j=1,NData gh(j) = -obs2(j) enddo do i=1,imx do j=1,NData gh(j) = gh(j) + datfma(j,i)*a2(i) enddo enddo endif endif eeerr = -32767. enerr = -32767. if (Calculate_Errors) then ! compute expected rms variation in field and expected error in fit. ! variation holds variation and errors holds error. txnorm = 2.*xnorm ! initialize v and ps to zero v = 0.0 ps = 0.0 !\ ! electric potential solution using cu !/ if (num == 1) then do ith=0,ithtrns sth = st(ith) do lon=1,lonmx do m=-mmx,mmx mm = iabs(m) fint = amie_f(m,lon)/ri dfint = m*amie_f(-m,lon)/sth/ri do i=ibm(m),iem(m) ! ps temporarily holds etheta components ! v temporarily holds ephi components ix = i - ibm(m) v(i) = q(2*ix+ns(mm),ith)*dfint ps(i) = dq(2*ix+ns(mm),ith)*fint enddo if (mm.eq.1.and.ith.eq.0) then do i=ibm(m),iem(m) v(i) = ps(i) enddo endif enddo do i=2,imx colsme(i-1) = 0. colsmn(i-1) = 0. culsm(i-1) = 0. do j=1,i-1 xe = v(i)*v(j) colsme(j) = colsme(j) + ax(i,j)*xe xn = ps(i)*ps(j) colsmn(j) = colsmn(j) + ax(i,j)*xn culsm(j) = culsm(j) + txnorm*cu(i,j)*(xe + xn) enddo enddo xe = v(imx)*v(imx) eeerr(lon,ith) = ax(imx,imx)*xe xn = ps(imx)*ps(imx) enerr(lon,ith) = ax(imx,imx)*xn variation(lon,ith) = xnorm*cu(imx,imx)*(xe + xn) do i=1,imx-1 xe = v(i)*v(i) eeerr(lon,ith) = eeerr(lon,ith) + colsme(i) + ax(i,i)*xe xn = ps(i)*ps(i) enerr(lon,ith) = enerr(lon,ith) + colsmn(i) + ax(i,i)*xn variation(lon,ith) = variation(lon,ith) + & culsm(i) + cu(i,i)*(xe + xn)*xnorm enddo x = eeerr(lon,ith) + enerr(lon,ith) if (x >= 0.) errors(lon,ith) = sqrt(x)*fac if (eeerr(lon,ith) >= 0.) eeerr(lon,ith) = sqrt(eeerr(lon,ith)) if (enerr(lon,ith) >= 0.) enerr(lon,ith) = sqrt(enerr(lon,ith)) if (variation(lon,ith) >= 0.) & variation(lon,ith) = sqrt(variation(lon,ith))*fac enddo eeerr(0,ith) = eeerr(lonmx,ith) enerr(0,ith) = enerr(lonmx,ith) variation(0,ith) = variation(lonmx,ith) errors(0,ith) = errors(lonmx,ith) enddo else ! ln(conductance, eflux, or mean e) solution using ccu do ith=0,ithtrns sth = st(ith) do lon=1,lonmx do m=-mmx,mmx mm = iabs(m) do i=ibm(m),iem(m) ! ps temporarily basis functions ! qsclr is defined to ithmx, ! but is reasonable to ithtrns+~2 (4/94) ix = i - ibm(m) v(i) = 0. ps(i) = qsclr(i-ibm(m)+nss(mm),ith)*amie_f(m,lon) enddo enddo do i=2,imx colsme(i-1) = 0. colsmn(i-1) = 0. culsm(i-1) = 0. do j=1,i-1 xe = v(i)*v(j) colsme(j) = colsme(j) + ax(i,j)*xe xn = ps(i)*ps(j) colsmn(j) = colsmn(j) + ax(i,j)*xn culsm(j) = culsm(j) + txnorm*ccu(i,j)*xn enddo enddo xe = v(imx)*v(imx) eeerr(lon,ith) = ax(imx,imx)*xe xn = ps(imx)*ps(imx) enerr(lon,ith) = ax(imx,imx)*xn variation(lon,ith) = xnorm * ccu(imx,imx)*xn do i=1,imx-1 xe = v(i)*v(i) eeerr(lon,ith) = eeerr(lon,ith) + colsme(i) + ax(i,i)*xe xn = ps(i)*ps(i) enerr(lon,ith) = enerr(lon,ith) + colsmn(i) + ax(i,i)*xn variation(lon,ith) = variation(lon,ith)+ & culsm(i)+ccu(i,i)*xn*xnorm enddo x = eeerr(lon,ith) + enerr(lon,ith) if (x >= 0.) errors(lon,ith) = sqrt(x)*fac if (eeerr(lon,ith) >= 0.) eeerr(lon,ith) = sqrt(eeerr(lon,ith)) if (enerr(lon,ith) >= 0.) enerr(lon,ith) = sqrt(enerr(lon,ith)) if (variation(lon,ith) >= 0.) & variation(lon,ith) = sqrt(variation(lon,ith))*fac enddo eeerr(0,ith) = eeerr(lonmx,ith) enerr(0,ith) = enerr(lonmx,ith) variation(0,ith) = variation(lonmx,ith) errors(0,ith) = errors(lonmx,ith) enddo endif endif end subroutine matsolv