c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/dynsatv.f c------------------------------------------------------------------ c subroutine dynsatv(utc,id) c c Make sat view projection contours: c include 'tgcmparam.h' include 'input.h' include 'tigcm.h' include 'tgcmhdr.h' include 'tiegcmfld.h' include 'tigcmfld.h' include 'selgrid.h' include 'tgcmlab.h' include 'cool.h' include 'plt.h' include 'color.h' dimension plotp(imx,jmx) character*56 toplab,lab,labgcm1 character*24 latlab character*18 radlab save rmintot,rmaxtot data rmintot,rmaxtot/1.e36,-1.e36/ c c Below are essentially dummies for bwcon and clrcon calls: parameter (nx=7,ny=7) dimension numx(nx),numy(ny) data numx/-180,-120,-60,0,60,120,180/ data numy/-90,-60,-30,0,30,60,90/ character*8 fmtx,fmty character*16 labx,laby data fmtx/'(i4) '/, fmty/'(i3) '/ data nfmtx/4/, nfmty/3/, mnrx/3/, mnry/3/ data laby/' LATITUDE (DEG) '/ data labx/' LONGITUDE (DEG)'/ c c Default distance from center of earth in earth radii: data sa /5.0/, dgrid/15./ c c Temp: data maps/1/, iccont/1/ c c Set up EZMAP: c proj = 'SV ' iel = 1 pmapl = .05 pmapr = .80 pmapb = .13 pmapt = .88 xmid = 0.5*(pmapl+pmapr) xc = -180. xd = +180. yc = -90. yd = +90. call mappos(pmapl,pmapr,pmapb,pmapt) call mapsti('EL',iel) call mapstc('OU','CO') call mapsti('DO',1) call mapstr('SA',sa) call mapset('MA',0.,0.,0.,0.) call mapsti('G1',1) call mapsti('G2',2) call mapsti('VS',0) call mapstr('GR',dgrid) call mapint c c Set up conpack: c (parameters common to all frames by this routine) c rilx = 0.5 rily = pmapb-0.20 call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) call cpsetr('ORV - out of range value',1.e12) call cpseti('MAP - mapping flag',1) call cpseti('SET',0) c c Loop through desired lat/lon points for center of projection: c do 100 l=1,nloc c c input longitude of r12flag means we want longitude at local noon: c projlat = rloc(1,l) projlon = rloc(2,l) if (rloc(2,l).eq.r12flag) then ixloc(2,l) = ixslt(12.,utc,projlon) endif c c Tell ezmap lat/lon and projection: c call maproj(proj,projlat,projlon,0.) call mapint call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c c Field loop: write(6,"(' ')") write(6,"('Satellite View projections (dynamo fields):')") do 200 ip=1,ntiegcmf if (iptiegcm(ip).le.0) goto 200 c c Selected pressure and heights loop: c do 110 k=1,npls+nhts kk = k-npls if (k.le.npls) then write(toplab,"('UT = ',f6.2,' ZP = ',f4.1,' DAY = ', + i5,23x)") utc,spls(k),iyd(id) else write(toplab,"('UT = ',f6.2,' HT = ',f6.1,' DAY = ', + i5,21x)") utc,shts(kk),iyd(id) endif c c At constant pressure surface: if (k.le.npls) then plotp(:,:) = dynpnt(:,ispls(k),:,ip) c c Interpolate to constant height surface: else call tgcmint(dynpnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ip),2,plotp,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) endif c c Take logs if necessary: if (idynlog(ip).gt.0.and.log10pl.gt.0) then do 125 j=1,jmx do 125 i=1,imx if (plotp(i,j).eq.cpspval.or.plotp(i,j).le.1.e-30) + then plotp(i,j) = cpspval else plotp(i,j) = alog10(plotp(i,j)) endif 125 continue write(lab,"('LOG10 ',a50)") labtiegcm(ip)(1:50) else lab = labtiegcm(ip) endif c c Find min/max at this projection: c (Used to fix contour levels for movies) rmin = 1.e36 rmax = -1.e36 do 130 j=1,jmx do 130 i=1,imx call maptrn(gcmlat(j),gcmlon(i),xpos,ypos) if (xpos.eq.1.e12.or.ypos.eq.1.e12) goto 130 if (plotp(i,j).ne.cpspval) then if (plotp(i,j).lt.rmin) rmin = plotp(i,j) if (plotp(i,j).gt.rmax) rmax = plotp(i,j) endif 130 continue rmintot = amin1(rmintot,rmin) rmaxtot = amax1(rmaxtot,rmax) c c Contour it, with axes and axis labels: c if (iclrfill.le.0) then call bwcon(plotp,imx,imx,jmx,xc,xd,yc,yd, + cintdyn(ip),cmindyn(ip),cmaxdyn(ip), + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) else lbar = 1 call clrcon(plotp,imx,imx,jmx,xc,xd,yc,yd, + cintdyn(ip),cmindyn(ip),cmaxdyn(ip), + 0,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,0.83,0.98,pmapb,pmapt,ip) endif c call box(1) c c Change color indices for continental outlines and lat/lon grid: c if (maps.gt.0) then call gqplci(ierr,isplci) call gqpmci(ierr,ispmci) if (iclrfill.le.0) then call gsplci(iwhite) call gspmci(iwhite) else call gsplci(iccont) call gspmci(iccont) endif c c Add continental outlines and lat/lon grid and labels c call mapsti('PE',1) call mapsti('LA',1) call maplot call mapgrd call maplbl call gsplci(isplci) call gspmci(ispmci) endif c c Mark location (center of projection): c szfs = .013 call maptrn(projlat,projlon,poslat,poslon) call plchhq(poslat,poslon,':KRL:G',szfs,0.,0.) c c Label location in upper right: c write(latlab,"(' at (',f5.1,',',f6.1,')',6x)") + projlat,projlon call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) call set(0.,1.,0.,1.,0.,1.,0.,1.,1) ypos = vt + 0.135*(vt-vb) call plchhq(0.70,ypos,':KRL:G',szfs,0.,0.) call plchhq(0.96,ypos,latlab(1:18),szfs,0.,1.) c c Label earth radii: c write(radlab,"('(',f4.1,' earth radii)')") sa ypos = ypos-0.045*(vt-vb) call plchhq(0.96,ypos,radlab,szfs,0.,1.) call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg) c c Mark longitude at local noon c rlon12 = 15.*(12.-utc) if (rlon12.ge.180.) rlon12 = rlon12-360. if (rlon12.lt.-180.) rlon12 = rlon12+360. c c For black and white, use stars near each pole, and at equator: c if (iclrfill.le.0) then rlat12 = 90.-1.3*dgrid call maptrn(rlat12,rlon12,xnoon,ynoon) szfs = .025 if (xnoon.ne.1.e12.and.ynoon.ne.1.e12) + call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.) rlat12 = -90.+1.3*dgrid call maptrn(rlat12,rlon12,xnoon,ynoon) if (xnoon.ne.1.e12.and.ynoon.ne.1.e12) + call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.) call maptrn(0.,rlon12,xnoon,ynoon) if (xnoon.ne.1.e12.and.ynoon.ne.1.e12) + call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.) else c c For color, draw wide red line along local noon meridian: c (default line width (LW) is 1000) c call setusv('LW',1900) call gsplci(ired) call gspmci(ired) c do 105 j=-90,88,2 c call maptrn(float(j),rlon12,xp,yp) c call maptrn(float(j+2),rlon12,xp1,yp1) c if (xp.ne.1.e12.and.yp.ne.1.e12.and. c + xp1.ne.1.e12.and.yp1.ne.1.e12) c + call line(xp,yp,xp1,yp1) c105 continue call setusv('LW',1000) call gsplci(iwhite) call gspmci(iwhite) endif c c Add top labels: c if (idynlog(ip).gt.0.and.log10pl.gt.0) then write(labgcm1,"('LOG10 ',a50)") labtiegcm(ip)(1:50) call top2lab(labgcm1,0.135,toplab,0.09,.017) else call top2lab(labtiegcm(ip),0.135,toplab,0.09,.017) endif write(lab,"(12x,a,12x)") runlab call botlab(lab,xmid,.05) c c Add neutral wind vectors if doing temp: c if (ip.eq.ixpot.and.ipltuv.eq.1) then if (k.le.npls) then udat(:,:) = ui(:,ispls(k),:) vdat(:,:) = vi(:,ispls(k),:) else c c Do height interpolation of ui,vi: c call tgcmint(ui,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ixui),2,udat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) call tgcmint(vi,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon,imx,shts(kk),1,idynlog(ixvi),2,vdat,imx,jmx, + 1,kmx,gcmzp,ier,1,cpspval,0) endif uvmax = 0. m = 1 ispvec = 0 spvec = 1.e36 iskip = 0 call pltvec(udat,vdat,uvmax,jmx,iskip,gcmlat(1), + gcmlat(jmx),m,ispvec,spvec) endif c c Finish up: call frame iframe=iframe+1 if (k.eq.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) then c write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, c + ' lat,lon=',2f8.2,/1x,' min,max = ',2e12.4, c + ' mintot maxtot = ',2e10.4)") c + iframe,labtiegcm_short(ip),utc,projlat,projlon, c + rmin,rmax,rmintot,rmaxtot write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, + ' lat,lon=',2f8.2)") + iframe,labtiegcm_short(ip),utc,projlat,projlon else if (k.le.npls) then c write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, c + ' ZP=',f4.1,' lat,lon=',2f8.2,/1x,' min,max = ', c + 2e12.4,' mintot maxtot=',2e12.4)") c + iframe,labtiegcm_short(ip),utc,spls(k),projlat,projlon, c + rmin,rmax,rmintot,rmaxtot write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, + ' ZP=',f4.1,' lat,lon=',2f8.2)") + iframe,labtiegcm_short(ip),utc,spls(k),projlat,projlon else c write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, c + ' HT=',f6.1,' lat,lon=',2f8.2,/1x,' min,max = ',2e12.4, c + ' mintot maxtot = ',2e12.4)") c + iframe,labtiegcm_short(ip),utc,shts(kk),projlat, c + projlon,rmin,rmax,rmintot,rmaxtot write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1, + ' HT=',f6.1,' lat,lon=',2f8.2)") + iframe,labtiegcm_short(ip),utc,shts(kk),projlat, + projlon endif endif c c End pressure/height loop: 110 continue c c End field loop: 200 continue c c End loc loop: 100 continue c return end