c c------------------------------------------------------------------ c Begin file /home/sting/foster/tigcm/pltutlat.f c------------------------------------------------------------------ c subroutine pltutlat(ic,it) c c Contour tigcm fields with ut on x-axis, latitude on y-axis c include 'tgcmparam.h' include 'tgcmhdr.h' include 'input.h' include 'tigcmfld.h' include 'tiegcmfld.h' include 'selgrid.h' include 'tigcm.h' include 'tgcmlab.h' include 'cool.h' include 'color.h' dimension flat(jmx),tncol(kmx),htcol(kmx), + flat1(jmx),flat2(jmx) c dimension tmpsave(241,jmx) ! temp for save file data icp/1/ c c Axes and labels: c dimension numy(7),numx(mxtms) character*56 toplab,lab character*8 fmtx,fmty character*16 labx,laby data fmty/'(i3) '/ data nfmty/3/, mnry/3/ data laby/' LATITUDE (DEG) '/ data numy /-90,-60,-30,0,30,60,90/ data ny/7/ save icp c c Write latitude slices to scratch file for current ut: c lu = luscr1-1 do 120 i=1,nlon+nslt c c Find index to selected longitude, or local time: if (i.le.nlon) then if (slon(i).eq.r12flag) then ilon = ixslt(12.,ut(it),slon(i)) else ilon = islon(i) if (ilon.eq.izmflag) then write(6,"('pltutlat: zonal means not yet available')") goto 120 endif endif else ilon = ixslt(sslt(i-nlon),ut(it),dum) endif c c Selected pressure and/or height loop: do 121 k=1,npls+nhts kk = k-npls lu = lu+1 if (ic.gt.icp) rewind lu c c Field loop: do 125 ip=1,ntigcmf+ntiegcmf ipp = ip-ntigcmf c c (need everything if need cooling terms): if (ip.le.ntigcmf) then if (iptigcm(ip).le.0.and.iplcool.le.0) goto 125 if ((ip.eq.ixfof2.or.ip.eq.ixhmf2).and.k.gt.1) goto 125 else c c dynamo fields (only first 4 fields available): if (idyn.le.0.or.iptiegcm(ipp).le.0.or.ipp.ge.5) goto 125 endif if (ip.gt.ntigcmf) goto 126 c c NOTE: changes to this conditional should also be made to c conditional in do loops 215 and 220 below c if (ip.le.ngcmflds) then if (k.le.npls) then flat(:) = pnt(ilon,ispls(k),:,ip) else call tgcmint(pnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon(ilon),1,shts(kk),1,igcmlog(ip),3,flat,jmx,1,1, + jmx,gcmzp,ier,0,cpspval,0) endif c c fof2 and hmf2: elseif (ip.eq.ixfof2.and.k.le.1) then flat(:) = fof2(ilon,:) elseif (ip.eq.ixhmf2.and.k.le.1) then flat(:) = hmf2(ilon,:) c c Total density: elseif (ip.eq.ixrho) then if (k.le.npls) then flat(:) = td(ilon,ispls(k),:) else call tgcmint(td,pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon(ilon),1,shts(kk),1,igcmlog(ip),3,flat,jmx,1,1, + jmx,gcmzp,ier,0,cpspval,0) endif c c Atomic/mol ratio: elseif (ip.eq.ixrat) then if (k.le.npls) then flat(:) = pnt(ilon,ispls(k),:,ixo1) / + (pnt(ilon,ispls(k),:,ixo2) + + pnt(ilon,ispls(k),:,ixn2)) else call tgcmint(pnt(1,1,1,ixo1),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon(ilon),1,shts(kk),1,igcmlog(ip),3,flat,jmx,1,1, + jmx,gcmzp,ier,0,cpspval,0) call tgcmint(pnt(1,1,1,ixo2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon(ilon),1,shts(kk),1,igcmlog(ip),3,flat1,jmx,1,1, + jmx,gcmzp,ier,0,cpspval,0) call tgcmint(pnt(1,1,1,ixn2),pnt(1,1,1,ixz),gcmlat,jmx, + gcmlon(ilon),1,shts(kk),1,igcmlog(ip),3,flat2,jmx,1,1, + jmx,gcmzp,ier,0,cpspval,0) flat(:) = flat(:) / (flat1(:) + flat2(:)) endif endif c c Dynamo fields: 126 continue if (ip.gt.ntigcmf) then if (k.le.npls) then flat(:) = dynpnt(ilon,ispls(k),:,ipp) else call tgcmint(dynpnt(1,1,1,ipp),pnt(1,1,1,ixz),gcmlat, + jmx,gcmlon(ilon),1,shts(kk),1,idynlog(ipp),3,flat, + jmx,1,1,jmx,gcmzp,ier,0,cpspval,0) endif endif c c Write to scratch file: write(lu) flat c c End field loop: 125 continue c c End selected pressure and/or height loop: 121 continue c c End selected longitude and/or slt loop: 120 continue c c Integrate cooling terms and write to scratch file: if (iplcool.gt.0) then do 140 ip=1,ncoolt do 140 i=1,nlon+nslt c c Find index to selected longitude, or local time: if (i.le.nlon) then if (slon(i).eq.r12flag) then ilon = ixslt(12.,ut(it),slon(i)) else ilon = islon(i) endif else ilon = ixslt(sslt(i-nlon),ut(it),dum) endif lu = lu+1 if (ic.gt.icp) rewind lu do 135 j=1,jmx flat(j) = 0. tncol(:) = pnt(ilon,:,j,ixt) htcol(:) = pnt(ilon,:,j,ixz) do 130 k=2,kmx c c at k: tn = pnt(ilon,k,j,ixt) xo2 = pnt(ilon,k,j,ixo2) xo = pnt(ilon,k,j,ixo1) xn4s = pnt(ilon,k,j,ixn4s) xno = pnt(ilon,k,j,ixno) xn2d = pnt(ilon,k,j,ixn2d) xn2 = pnt(ilon,k,j,ixn2) xne = pnt(ilon,k,j,ixne) c c at k-1: tnm1 = pnt(ilon,k-1,j,ixt) xo2m1= pnt(ilon,k-1,j,ixo2) xom1 = pnt(ilon,k-1,j,ixo1) xn4sm1 = pnt(ilon,k-1,j,ixn4s) xnom1 = pnt(ilon,k-1,j,ixno) xn2dm1 = pnt(ilon,k-1,j,ixn2d) xn2m1 = pnt(ilon,k-1,j,ixn2) xnem1 = pnt(ilon,k-1,j,ixne) c c do integration: rz=1.66e-24*(16.*xo + 32.*xo2 + 28.*xn2) rzm=1.66e-24*(16.*xom1 + 32.*xo2m1 + 28.*xn2m1) coolk = cool(tn,xo2,xo,xn4s,xno,xn2d,xn2,xne,k, + tncol,htcol,ip) coolkm1 = cool(tnm1,xo2m1,xom1,xn4sm1,xnom1,xn2dm1, + xn2m1,xnem1,k-1,tncol,htcol,ip) rw = 0.5*(pnt(ilon,k,j,ixz) - pnt(ilon,k-1,j,ixz)) + *1.e+5 flat(j) = flat(j) + (coolk*rz + coolkm1*rzm) * rw 130 continue 135 continue write(lu) flat 140 continue endif c c If just wrote last ut, jump down to make plots: if (it.eq.ntms) goto 500 c c Otherwise exit, since not doing plots yet: goto 501 c c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> c c If doing last ut, make plots: c 500 continue c c Set up conpack: c call cpseti('SET',0) call cpseti('MAP - mapping flag',0) utmapl = 0.12 utmapr = 0.90 utmapb = 0.15 utmapt = 0.85 rilx = 0.5*(utmapl+utmapr) rily = utmapb - 0.35 c c Make room for color bar at bottom if doing color c (this will squash plot in y direction slightly compared to monochrome) if (iclrfill.gt.0) then utmapt = 0.92 utmapb = 0.30 rily = -0.240 endif call cpsetr('ILX - info label x-coord',rilx) call cpsetr('ILY - info label y-coord',rily) utinc1 = float(mhut(1)) utincn = utinc1 + float(ntms-1) xc = utinc1 xd = utincn call set(utmapl,utmapr,utmapb,utmapt,utinc1,utincn,-90.,+90.,1) yc = gcmlat(1) yd = gcmlat(jmx) c c Set up x-axis: c nx = 0 if (ndays.eq.1) then mnrx = 3 write(labx,"(' UT (day ',i5,') ')") iyd(1) do 45 i=1,ntms if (mod(mhut(i),6).eq.0) then nx = nx+1 numx(nx) = mhut(i) endif 45 continue elseif (ndays.le.3) then mnrx = 3 write(labx,"('days ',i5,'-',i5)") iyd(1),iyd(ndays) do 50 i=1,ntms if (mod(mhut(i),6).eq.0) then nx = nx+1 numx(nx) = mhut(i) endif 50 continue else mnrx = 4 do 55 i=1,ndays write(labx,"(' DAY OF ',i4,3x)") 1900+iyd(i)/1000 nx = nx+1 numx(nx) = iyd(i)-iyd(i)/1000*1000 55 continue nx = nx+1 iydlast = iyd(ndays)+1 numx(nx) = iydlast-iydlast/1000*1000 endif fmtx = '(i2) ' nfmtx = 2 if (numx(1).gt.99) then fmtx = '(i3) ' nfmtx = 3 endif c c Read fields back from scratch file: c lu = luscr1-1 do 200 i=1,nlon+nslt do 250 k=1,npls+nhts kk = k-npls write(6,"(' ')") if (i.le.nlon) then ! selected longitude if (k.le.npls) then ! selected pressure write(6,"('Pltutlat: lon=',f7.2,' zp=',f4.1)") + slon(i),spls(k) else ! selected height write(6,"('Pltutlat: lon=',f7.2,' ht=',f6.1)") + slon(i),shts(kk) endif else ! selected local time if (k.le.npls) then ! selected pressure write(6,"('Pltutlat: slt=',f7.2,' zp=',f4.1)") + sslt(i-nlon),spls(k) else ! selected height write(6,"('Pltutlat: slt=',f7.2,' ht=',f6.1)") + sslt(i-nlon),shts(kk) endif endif lu = lu + 1 rewind lu do 210 iit=1,ntms do 215 ip=1,ntigcmf+ntiegcmf ipp = ip-ntigcmf if (ip.le.ntigcmf) then if (iptigcm(ip).le.0.and.iplcool.le.0) goto 215 if ((ip.eq.ixfof2.or.ip.eq.ixhmf2).and.k.gt.1) goto 215 else if (idyn.le.0.or.iptiegcm(ipp).le.0) goto 215 endif c c Read from scratch file: read(lu) flat gcmuts(iit,:,ip) = flat(:) 215 continue 210 continue c c Contour it: c do 220 ip=1,ntigcmf+ntiegcmf ipp = ip-ntigcmf if (ip.le.ntigcmf) then if (iptigcm(ip).le.0.and.iplcool.le.0) goto 220 if ((ip.eq.ixfof2.or.ip.eq.ixhmf2).and.k.gt.1) goto 220 else if (idyn.le.0.or.iptiegcm(ipp).le.0) goto 220 endif c c fof2 and hmf2 are height independent: if ((ip.eq.ixfof2.or.ip.eq.ixhmf2).and.k.le.1) then if (i.le.nlon) then write(toplab,"(a,' LON=',f7.2,8x)") + runlab(1:lnblnk2(runlab)),slon(i) else write(toplab,"(a,' SLT=',f7.2,8x)") + runlab(1:lnblnk2(runlab)),sslt(i-nlon) endif else if (i.le.nlon) then if (k.le.npls) then write(toplab,"(a,' LON=',f7.2,' ZP=',f4.1)") + runlab(1:lnblnk2(runlab)),slon(i),spls(k) else write(toplab,"(a,' LON=',f7.2,' HT=',f5.1)") + runlab(1:lnblnk2(runlab)),slon(i),shts(kk) endif else if (k.le.npls) then write(toplab,"(a,' SLT=',f7.2,' ZP=',f4.1)") + runlab(1:lnblnk2(runlab)),sslt(i-nlon),spls(k) else write(toplab,"(a,' SLT=',f7.2,' HT=',f5.1)") + runlab(1:lnblnk2(runlab)),sslt(i-nlon),shts(kk) endif endif endif if (ip.le.ntigcmf) then cint = conints(ip) cmin = conmins(ip) cmax = conmaxs(ip) ipfm = ip else cint = cintdyn(ipp) cmin = cmindyn(ipp) cmax = cmaxdyn(ipp) ipfm = ip+ntigcmf endif c c Temp -- write data to save file: c c do ii=1,ntms c tmpsave(ii,:) = gcmuts(ii,:,ip) c enddo c write(90) tmpsave c write(6,"('pltutlat: xc xd=',2f9.3,' yc yd=',2f9.3, c + ' ntms=',i4,' ip=',i3)") xc,xd,yc,yd,ntms,ip c c End Temp c if (iclrfill.le.0) then call bwcon(gcmuts(1,1,ip),mxtms,ntms,jmx,xc,xd,yc,yd, + cint,cmin,cmax, + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) else lbar = 2 call clrcon(gcmuts(1,1,ip),mxtms,ntms,jmx,xc,xd,yc,yd, + cint,cmin,cmax, + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty, + lbar,utmapl,utmapr,.05,.13,ipfm) endif c call box(1) if (ip.le.ntigcmf) then call top2lab(labtigcm(ip),0.10,toplab,0.05,.017) else call top2lab(labtiegcm(ipp),0.10,toplab,0.05,.017) endif write(lab,"(8x,'TIGCM History = ',3a8,8x)") + (histvol(ii,ivf),ii=1,3) ybotlab = .05 if (iclrfill.gt.0) ybotlab = 0.20 c call botlab(lab,0.5,ybotlab) call frame iframe=iframe+1 if (ip.le.ntigcmf) then write(6,"('Pltutlat frame ',i3,': ',a)") iframe,labtigcm(ip) else write(6,"('Pltutlat frame ',i3,': ',a)") iframe, + labtiegcm(ipp) endif c c End field loop: 220 continue c c End selected pressure/height loop: 250 continue c c End selected lon/slt loop: 200 continue c c Height integrated cooling/heating terms: c (there are ncoolt terms: c no heat, no cool, no heat-cool, co2 cool, o3p cool, irt) c if (iplcool.gt.0) then do 305 ip=1,ncoolt do 305 i=1,nlon+nslt if (i.le.nlon) then write(toplab,"(a,' LON=',f7.2,' (Ht integrated)')") + runlab(1:lnblnk2(runlab)),slon(i) else write(toplab,"(a,' SLT=',f7.2,' (Ht integrated)')") + runlab(1:lnblnk2(runlab)),sslt(i-nlon) endif lu = lu + 1 rewind lu do 300 iit=1,ntms read(lu) flat cooluts(iit,:,ip) = flat(:) 300 continue call bwcon(cooluts(1,1,ip),mxtms,ntms,jmx,xc,xd,yc,yd, + 0.,1.,0., + 1,nx,numx,mnrx,labx,fmtx,nfmtx, + ny,numy,mnry,laby,fmty,nfmty) call top2lab(labcool(ip),0.10,toplab,0.05,.017) write(lab,"(8x,'TIGCM History = ',3a8,8x)") + (histvol(ii,ivf),ii=1,3) call botlab(lab,0.5,.05) call frame iframe=iframe+1 if (i.le.nlon) then write(6,"('Pltutlat frame ',i3,': ',a,' lon=',f6.1)") + iframe,labcool(ip)(1:43),slon(i) else write(6,"('Pltutlat frame ',i3,': ',a,' slt=',f6.1)") + iframe,labcool(ip)(1:43),sslt(i-nlon) endif 305 continue endif c c Doing plots, or current ut: 501 continue icp = ic return end