c subroutine mkutkmx(fields,htutkmx,utkmx,glat,glon,it,l,idop) include 'timesloc.h' include 'tgcmhdr.h' c c call mkutkmx(fields,htutkmx,utkmx,glat,glon,it,l,1) c real fkmx10(kmx,10),hts(kmx),fields(imx,kmx,jmx,nfget), + utkmx(ntms,kmx,nfplt),htutkmx(ntms,kmx,nloc) logical didions pointer(pfhts,fhts(1)) c didions = .false. do ip=1,nftot if (ifplt(ip).gt.0) then if (ip.le.nfhist) then call defutkmx(fields(1,1,1,ifget(ip)), + utkmx(1,1,ifplt(ip),l),glat,glon,it) else c c Desired field is not on history, and must be calculated: c c Ion drifts: c (use fkmx10(kmx,10) for t,u,v,w,o2,o,n2) c (igrf (from ~/lib/igrf) is used to determine dec and dip at loc) c (note igrf extrapolates beyond 1985 -- avoid this) c Zonal or global means of ion drifts with efield not available c (this would require finding average rdip,rdec by multiple calls c to igrf, and calling efield over all lons or entire globe, so c is not worth it at this time) c if (ip.ge.ixui.and.ip.le.ixwi) then if (.not.didions) then if (ionvel.eq.2.and.(glat.eq.zmflag.or. + glon.eq.zmflag)) then didions = .true. goto 100 endif ! means with efield not available altgrf = 300. yr = float(date(1)) if (yr.gt.1985.) yr = 1985. call igrf(glat,glon,altgrf,yr,gmlat,gmlon,rdip,rdec) rdip = rdip*dtr rdec = rdec*dtr call defkmx(fields(1,1,1,ifget(ixu)),fkmx10, + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,2), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixw)),fkmx10(1,3), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,4), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,5), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,6), + glat,glon,it) c c subroutine mkionvel(fout,kmx,glat,glon,gmlat,gmlon,rdip,rdec, c + day,ut,u,v,w,xo2,xo,xn2,iefield,spv) c call mkionvel(fkmx10(1,8),kmx,glat,glon,gmlat, + gmlon,rdip,rdec,float(date(2)),ut(it),fkmx10, + fkmx10(1,2),fkmx10(1,3),fkmx10(1,4),fkmx10(1,5), + fkmx10(1,6),ionvel-1,cpspval) utkmx(it,:,ifplt(ixui),l) = fkmx10(:,8) utkmx(it,:,ifplt(ixvi),l) = fkmx10(:,9) utkmx(it,:,ifplt(ixwi),l) = fkmx10(:,10) c c If doing linear height scale on y-axis, we need to interpolate c the 6 needed neutral fields first, then call mkionvel to calculate c drifts from the interpolated fields; this defines the array c htions(ntms,nhtscale,3,nloc) (phtions pointer in common, allocated c in getinp), which is plotted from pltloc: c if (nhtscale.gt.0) then call alloc(pfhts,nhtscale*10) hts(:) = htutkmx(it,:,l) do i=1,6 call intloc(fkmx10(1,i),hts,kmx, + htscale,nhtscale,0,fhts(nhtscale*(i-1)+1), + 1,1,nhtscale,ier,cpspval,0) enddo call mkionvel(fhts(nhtscale*7+1),nhtscale,glat,glon, + gmlat,gmlon,rdip,rdec,float(date(2)),ut(it), + fhts,fhts(nhtscale+1),fhts(nhtscale*2+1), + fhts(nhtscale*3+1),fhts(nhtscale*4+1), + fhts(nhtscale*5+1),ionvel-1,cpspval) c c pointers: htions(ntms,nhtscale,3,nloc) and fhts(nhtscale,10), c with ui,vi,wi in fhts(:,8-10) c do i=1,3 do k=1,nhtscale ix0 = (l-1)*3*nhtscale*ntms+(i-1)*nhtscale* + ntms+(k-1)*ntms+it ix1 = (i+6)*nhtscale+k htions(ix0) = fhts(ix1) enddo enddo call hpdeallc(pfhts,ier,1) endif didions = .true. endif c c E5577 (greenline volume emission, see mke5577): c (note defkmx returns loc, global, or zonal means as needed) c elseif (ip.eq.ixe5577) then call defkmx(fields(1,1,1,ifget(ixt)),fkmx10, + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,2), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,3), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,4), + glat,glon,it) if (idop5577.gt.0) then ! put un,vn in fkmx10 5,6 call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,5), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,6), + glat,glon,it) sigma = 17924. ! for doppler from greenline endif call mke5577(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),kmx,fkmx10(1,10),cpspval) utkmx(it,:,ifplt(ixe5577),l) = fkmx10(:,10) if (nhtscale.gt.0) then call alloc(pfhts,nhtscale*7) hts(:) = htutkmx(it,:,l) do i=1,4 ! interpolate tn,o,o2,n2 to htscale lg = 0 if (i.gt.1) lg = 1 call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale, + lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier, + cpspval,0) enddo call mke5577(fhts,fhts(nhtscale+1), + fhts(nhtscale*2+1),fhts(nhtscale*3+1), + nhtscale,fhts(nhtscale*6+1),cpspval) c c pointers: hte5577(ntms,nhtscale,nloc), fhts(nhtscale,7), with c e5577 currently in fhts(:,7) c do k=1,nhtscale ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it hte5577(ix0) = fhts(nhtscale*6+k) enddo c c If doing doppler (i.e., doing e5577, nhtscale > 0), c then define 4 doppler fields: tn,un,vn,brightness in dop5577(ntms,nloc,4) c (dop5577 is pdop pointee, pdop5577 is in common) c (doppler and quadrat routine in ~/lib/doppler.f) c SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma, c + VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ) c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER.. c (t,u,v currently in fkmx10, slots 1,5,6; put ver in fkmx10(:,10): c c interpolate u,v -> fhts 5,6 if (idop5577.gt.0.and.idop.gt.0) then call mke5577(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),kmx,fkmx10(1,10),cpspval) hts(:) = htutkmx(it,:,l) call doppler(kmx,hts,fkmx10,fkmx10(1,5), + fkmx10(1,6),fkmx10(1,10),sigma,dopvn,dopun,doptn, + dattpv,attpz,dattpz) hts(:) = htutkmx(it,:,l)*1.e5 dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6 do ipp=1,4 ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it if (ipp.eq.1) dop5577(ix) = doptn if (ipp.eq.2) dop5577(ix) = dopun if (ipp.eq.3) dop5577(ix) = dopvn if (ipp.eq.4) dop5577(ix) = dopb enddo endif call hpdeallc(pfhts,ier,1) endif ! nhtscale > 0 c c E6300 (redline volume emission, see mke6300): c (note defkmx returns loc, global, or zonal means as needed) c elseif (ip.eq.ixe6300) then call defkmx(fields(1,1,1,ifget(ixt)),fkmx10, + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixte)),fkmx10(1,2), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,3), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,4), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,5), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo2p)),fkmx10(1,6), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixne)),fkmx10(1,7), + glat,glon,it) if (idop6300.gt.0) then ! put un,vn in fkmx10 8,9 call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,8), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,9), + glat,glon,it) sigma = 15867.852 ! for doppler from redline endif c pass in t,te,o2,o,n2,o2+, get out redline ver in slot 10 of fkmx10: c subroutine mke6300(tn,te,xo2,xo,xn2,xo2p,xne,id,fout,spval) call mke6300(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),fkmx10(1,5),fkmx10(1,6),fkmx10(1,7),kmx, + fkmx10(1,10),cpspval) utkmx(it,:,ifplt(ixe6300),l) = fkmx10(:,10) if (nhtscale.gt.0) then call alloc(pfhts,nhtscale*10) hts(:) = htutkmx(it,:,l) do i=1,7 ! interpolate tn,te,o,o2,n2,o2+,ne to htscale lg = 0 if (i.gt.2) lg = 1 call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale, + lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier, + cpspval,0) enddo c ht interpolated ver will be in fhts(:,10): call mke6300(fhts,fhts(nhtscale+1), + fhts(nhtscale*2+1),fhts(nhtscale*3+1), + fhts(nhtscale*4+1),fhts(nhtscale*5+1), + fhts(nhtscale*6+1), + nhtscale,fhts(nhtscale*9+1),cpspval) c c pointers: hte6300(ntms,nhtscale,nloc), fhts(nhtscale,10), with c e6300 currently in fhts(:,10) c do k=1,nhtscale ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it hte6300(ix0) = fhts(nhtscale*9+k) enddo c c If doing doppler (i.e., doing e6300, nhtscale > 0) c then define 4 doppler fields: tn,un,vn,brightness in dop6300(ntms,nloc,4) c (dop6300 is pdop pointee, pdop6300 is in common) c (doppler and quadrat routine in ~/lib/doppler.f) c SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma, c + VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ) c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER.. c (t,u,v currently in fkmx10, slots 1,8,9; put ver in fkmx10(:,10): c if (idop6300.gt.0.and.idop.gt.0) then call mke6300(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),fkmx10(1,5),fkmx10(1,6),fkmx10(1,7), + kmx,fkmx10(1,10),cpspval) hts(:) = htutkmx(it,:,l) call doppler(kmx,hts,fkmx10,fkmx10(1,8), + fkmx10(1,9),fkmx10(1,10),sigma,dopvn,dopun,doptn, + dattpv,attpz,dattpz) hts(:) = htutkmx(it,:,l)*1.e5 dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6 do ipp=1,4 ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it if (ipp.eq.1) dop6300(ix) = doptn if (ipp.eq.2) dop6300(ix) = dopun if (ipp.eq.3) dop6300(ix) = dopvn if (ipp.eq.4) dop6300(ix) = dopb enddo endif call hpdeallc(pfhts,ier,1) endif ! nhtscale > 0 c c Eo200 (O2 (0-0) atmos band emission, see ~/lib/mkeo200): c (note defkmx returns loc, global, or zonal means as needed) c elseif (ip.eq.ixeo200) then call defkmx(fields(1,1,1,ifget(ixt)),fkmx10, + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo2)),fkmx10(1,2), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixo1)),fkmx10(1,3), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixn2)),fkmx10(1,4), + glat,glon,it) if (idopo200.gt.0) then ! put un,vn in fkmx10 5,6 call defkmx(fields(1,1,1,ifget(ixu)),fkmx10(1,5), + glat,glon,it) call defkmx(fields(1,1,1,ifget(ixv)),fkmx10(1,6), + glat,glon,it) sigma = 17924. ! for doppler from greenline endif call mkeo200(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),kmx,fkmx10(1,10),cpspval) utkmx(it,:,ifplt(ixeo200),l) = fkmx10(:,10) if (nhtscale.gt.0) then call alloc(pfhts,nhtscale*7) hts(:) = htutkmx(it,:,l) do i=1,4 ! interpolate tn,o,o2,n2 to htscale lg = 0 if (i.gt.1) lg = 1 call intloc(fkmx10(1,i),hts,kmx,htscale,nhtscale, + lg,fhts(nhtscale*(i-1)+1),1,1,nhtscale,ier, + cpspval,0) enddo call mkeo200(fhts,fhts(nhtscale+1), + fhts(nhtscale*2+1),fhts(nhtscale*3+1), + nhtscale,fhts(nhtscale*6+1),cpspval) c c pointers: hteo200(ntms,nhtscale,nloc), fhts(nhtscale,7), with c eo200 currently in fhts(:,7) c do k=1,nhtscale ix0 = (l-1)*nhtscale*ntms + (k-1)*ntms + it hteo200(ix0) = fhts(nhtscale*6+k) enddo c c If doing doppler (i.e., doing eo200, nhtscale > 0), c then define 4 doppler fields: tn,un,vn,brightness in dopo200(ntms,nloc,4) c (dopo200 is pdop pointee, pdopo200 is in common) c (doppler and quadrat routine in ~/lib/doppler.f) c SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma, c + VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ) c doppler takes ht,t,u,v, and ver as input (..HM,TM,ZM,VM,VER.. c (t,u,v currently in fkmx10, slots 1,5,6; put ver in fkmx10(:,10): c c interpolate u,v -> fhts 5,6 if (idopo200.gt.0.and.idop.gt.0) then call mkeo200(fkmx10,fkmx10(1,2),fkmx10(1,3), + fkmx10(1,4),kmx,fkmx10(1,10),cpspval) hts(:) = htutkmx(it,:,l) call doppler(kmx,hts,fkmx10,fkmx10(1,5), + fkmx10(1,6),fkmx10(1,10),sigma,dopvn,dopun,doptn, + dattpv,attpz,dattpz) hts(:) = htutkmx(it,:,l)*1.e5 dopb = quadrat(kmx,hts,fkmx10(1,10))*1.e-6 do ipp=1,4 ix = (ipp-1)*nloc*ntms+(l-1)*ntms+it if (ipp.eq.1) dopo200(ix) = doptn if (ipp.eq.2) dopo200(ix) = dopun if (ipp.eq.3) dopo200(ix) = dopvn if (ipp.eq.4) dopo200(ix) = dopb enddo endif call hpdeallc(pfhts,ier,1) endif ! nhtscale > 0 else write(6,"('mkfields: unknown ip=',i3)") ip endif endif endif 100 continue enddo ! fields loop return end