c c------------------------------------------------------------------ c subroutine mkuivi(phi,gcmht,gcmlat,velout,ivel,lu, + day,ut,un,vn,wn,xo2,xo,xn2,dipdec,ionvel,ivevu,vemag,vumag,spv) c include 'getgcm.h' c c Calculate ui,vi,wi from electric potential (input phi), and from c magnetic field coords bx,by,bz read from file assigned to unit lu c Magnetic field strengths suitable for assigning to lu are on mss c /FOSTER/magfield/magfield.cray c Modified 10/93 from the getgcm.a version to include option of adding c effects of neutral atmosphere (for tigcmloc): c c On input: c phi = global array (imx,kmx,jmx) of electric potential c gcmht = global array of model heights c gcmlat = array of model latitudes c velout = global array to receive ion drift velocity c ivel = 1,2,3 for ui,vi, or wi to be returned in velout c lu = fortran logical unit already attached to file containing magnetic c field strengths bx,by,bz (see above). A check is made to determine c if this unit is attached, and if it is not, a warning message c is printed to standard out. c day,ut = julian day and ut c ionvel = flag for whether or not to add neutral atmos effects: c ionvel = 1 -> ExB only c ionvel = 2 -> ExB + neutral effects c un,vn,wn,xo2,xo,xn2 (imx,kmx,jmx)= inputs necessary if ionvel = 2 c (may be dummies if ionvel = 1) c dipdec(imx,jmx,2) = dip and declination (radians) on global grid c (ignored if ionvel=1) c spv = special value c c On output: c if ivel=1, velout contains zonal drift velocities c if ivel=2, velout contains meridional drift velocities c if ivel=3, velout contains vertical drift velocities c (all velocities in m/s) c phi, gcmht, and gcmlat are untouched (as are neutral input params) c vemag and vumag are saved only if ionvel=2. They are Art's ion drifts c perpendicular to mag field where ve is +east and vu is +up c parameter(imxm1=imx-1,imxp3=imx+3,jmxp1=jmx+1,jmxp2=jmx+2, + jmxm1=jmx-1,imxp2=imx+2,dlon=5.,dlat=5.) parameter (pi=3.14159,dtr=pi/180.,rs=6.475165e+8) dimension phi(imx,kmx,jmx),gcmht(imx,kmx,jmx),un(imx,kmx,jmx), + vn(imx,kmx,jmx),wn(imx,kmx,jmx),xo2(imx,kmx,jmx),xo(imx,kmx,jmx), + xn2(imx,kmx,jmx),dipdec(imx,jmx,2),vemag(imx,kmx,jmx), + vumag(imx,kmx,jmx) dimension velout(imx,kmx,jmx) dimension wrk(imxp3,kmx,jmxp2),b(imx,jmx),gcmlat(jmx) dimension bx(imx,jmx),by(imx,jmx),bz(imx,jmx) dimension ex(kmx),ey(kmx),ez(kmx),uitmp(kmx),vitmp(kmx),witmp(kmx) logical isthere data scale/1.e6/ data ncalls/0/ data bgaus/.4/, xmaso/32./ save bx,by,bz,ncalls c ncalls = ncalls+1 c c Define work array with wrap around points: c do 100 k=1,kmx do 105 j=1,jmx do 105 i=1,imxm1 wrk(i+2,k,j+1) = phi(i,k,j) 105 continue c c Over the poles: c do 110 i=1,imxm1 wrk(i+2,k,1) = phi(1+mod(i+jmxm1,imxm1),k,1) wrk(i+2,k,jmxp2) = phi(1+mod(i+jmxm1,imxm1),k,jmx) 110 continue c c Periodic points: c do 115 j=1,jmxp2 do 115 i=1,2 wrk(i,k,j) = wrk(i+imxm1,k,j) wrk(i+imxp1,k,j) = wrk(i+2,k,j) 115 continue 100 continue c c Read magnetic coords: c if (ncalls.eq.1) then inquire(lu,exist=isthere) if (.not.isthere) then write(6,"('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>', + '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>')") write(6,"('mkuivi warning: unit ',i3,' is not attached', + /'It needs to be assigned to a file containing magnetic', + ' field strengths',/'Such a file exists on mss ', + '/FOSTER/magfield/magfield.cray')") lu write(6,"('>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>', + '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>')") endif rewind lu read(lu) bx,by,bz endif c c Grid loop: c do 200 j=1,jmx do 200 i=1,imx bxyzsum = bx(i,j)**2 + by(i,j)**2 + bz(i,j)**2 b(i,j) = bxyzsum**0.5 do 205 k=1,kmx c c Calculation of electric field: c ex(k) = -(wrk(i+3,k,j+1) - wrk(i+1,k,j+1)) / + (2.*dlon*dtr*(rs+gcmht(i,k,j)*1.e+5)* + cos(gcmlat(j)*dtr)) ey(k) = -(wrk(i+2,k,j+2) - wrk(i+2,k,j)) / + (2.*dlat*dtr*(rs+gcmht(i,k,j)*1.e+5)) 205 continue do 210 k=2,kmx-1 ez(k) = -((phi(i,k+1,j)-phi(i,k,j)) / + ((gcmht(i,k+1,j)-gcmht(i,k,j))*1.e+5) + + (phi(i,k,j)-phi(i,k-1,j)) / + ((gcmht(i,k,j)-gcmht(i,k-1,j))*1.e+5)) / 2. 210 continue ez(1) = ez(2)*2.-ez(3) ez(kmx) = ez(kmx-1)*2.-ez(kmx-2) c c Find the ExB drifts (m/s): c if (ionvel.lt.2) then if (ivel.eq.1) then do k=1,kmx velout(i,k,j) = (ey(k)*bz(i,j) - ez(k)*by(i,j)) * scale / + b(i,j)**2 enddo elseif (ivel.eq.2) then do k=1,kmx velout(i,k,j) = (ez(k)*bx(i,j) - ex(k)*bz(i,j)) * scale / + b(i,j)**2 enddo elseif (ivel.eq.3) then do k=1,kmx velout(i,k,j) = (ex(k)*by(i,j) - ey(k)*bx(i,j)) * scale / + b(i,j)**2 enddo endif c c Calculate effects of neutral atmosphere: c (taken partially from ~/timesloc/mkionvel.f and from communication with Art) c (note: Art may like to see line plots of ve and vu over time at selected c heights, like zphtline) c (11/10/93: moved ve,vu calculation into k loop where ui,vi,wi are c available) c else h = sqrt(bx(i,j)**2+by(i,j)**2) cosd = cos(dipdec(i,j,2)) sind = sin(dipdec(i,j,2)) si = sin(dipdec(i,j,1)) ci = cos(dipdec(i,j,1)) do k=1,kmx ui = (ey(k)*bz(i,j) - ez(k)*by(i,j)) * scale / b(i,j)**2 vi = (ez(k)*bx(i,j) - ex(k)*bz(i,j)) * scale / b(i,j)**2 wi = (ex(k)*by(i,j) - ey(k)*bx(i,j)) * scale / b(i,j)**2 ve = (ui*by(i,j)-vi*bx(i,j)) / h vu = -((ui*bx(i,j)+vi*by(i,j))/(h*b(i,j))*bz(i,j))+ + h/b(i,j)*wi if (ivevu.gt.0) then vemag(i,k,j) = ve vumag(i,k,j) = vu endif tun = un(i,k,j) tvn = vn(i,k,j) twn = wn(i,k,j) to1= xo(i,k,j) to2 = xo2(i,k,j) tn2 = xn2(i,k,j) if (to1.ne.spv.and.tn2.ne.spv.and.to2.ne.spv.and. + tun.ne.spv.and.tvn.ne.spv.and.twn.ne.spv) then vinp = 3.e-10 * (to1+tn2+to2) rhoi = vinp / (9.57946e+3 * bgaus / xmaso) coi = 1. / (1. + rhoi*rhoi) unmag = tun*cosd - tvn*sind vnmag = tun*sind + tvn*cosd uitmp(k) =coi*(rhoi*rhoi*unmag+rhoi*(-vnmag*si-twn*ci))+ + coi*(ve+rhoi*vu) vitmp(k) = coi*(rhoi*rhoi*vnmag+rhoi*unmag*si+ + vnmag*ci*ci-twn*si*ci)+coi*(vu*si-rhoi*ve*si) witmp(k) = coi*(rhoi*rhoi*twn+rhoi*unmag*ci-vnmag*si*ci+ + twn*si*si)+coi*(vu*ci-rhoi*ve*ci) else uitmp(k) = spv vitmp(k) = spv witmp(k) = spv endif enddo c c Add to output field: c if (ivel.eq.1) then do k=1,kmx if (uitmp(k).ne.spv.and.velout(i,k,j).ne.spv) then velout(i,k,j) = uitmp(k) else velout(i,k,j) = spv endif enddo endif if (ivel.eq.2) then do k=1,kmx if (vitmp(k).ne.spv.and.velout(i,k,j).ne.spv) then velout(i,k,j) = vitmp(k) else velout(i,k,j) = spv endif enddo endif if (ivel.eq.3) then do k=1,kmx if (witmp(k).ne.spv.and.velout(i,k,j).ne.spv) then velout(i,k,j) = witmp(k) else velout(i,k,j) = spv endif enddo endif endif ! ionvel= 1 or 2 200 continue return end