c subroutine mkmaps(mtime) c c Contour fields over mapped projections: c (projections are made and fields contoured by pltce, pltst, pltsv, c which are in lib plt.a) c include "flxproc.h" include "ccm.h" include "pltopt.h" real fglb0(imx,jmx), uglb(imx,jmx), vglb(imx,jmx), tlabchsz(3), + blabchsz(3) real vpce(4),vpst(4),vpsv(4),vp(4) integer inc(jmx),mtime(3) character*16 modstr character*4 modstr4 character*56 fldlabsave character*80 msgout,infolab,histlab character*96 tlabs(3),blabs(3) character*2 proj logical didgcmbot,isgcmbot data vpce /.10,.96,.45,.88/ data vpst /.155,.845,.180,.870/ data vpsv /.15,.85,.19,.89/ data tlabchsz/.02,.02,.02/, blabchsz/.015,.018,.018/ data ixuivi/0/,ixunvn/0/ ! not implemented c c inc(jmx) = increment in longitude at which to plot vectors around c each latitude (for polar stereographic and sat view projs only) c data inc / c -87.5 -82.5 -77.5 -72.5 -67.5 -62.5 -57.5 -52.5 -47.5 -42.5 + 24, 8, 6, 4, 2, 2, 2, 2, 2, 2, c -37.5 -32.5 -27.5 -22.5 -17.5 -12.5 -7.5 -2.5 2.5 7.5 + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, c 12.5 17.5 22.5 27.5 32.5 37.5 42.5 47.5 52.5 57.5 + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, c 62.5 67.5 72.5 77.5 82.5 87.5 + 2, 2, 4, 6, 8, 24/ c write(6,"(/'Map projections at ut = ',f4.1,' tgcm mtime = ',3i3, + ' ccmday = ',f10.4)") ut,(mtime(i),i=1,3),ccmday c c Write histlab for wrdat: c c if (iwrdat.gt.0) c + write(histlab,"(a,' ',a,' DAY,HR,MIN=',i3,',',i2,',',i2)") c + model(1:lenstr(model)), c + histvols(ivol)(1:lenstr(histvols(ivol))), c + (mtime(i),i=1,3) c c Field loop: c (non-zero values of ifproc are numbered in user defined order) c (common nfproc is number of fields to be plotted (max of ifproc values), c and nplt is current field in user defined order) c nplt = 1 101 continue do ip=1,mxfproc if (ifproc(ip).ne.nplt) goto 100 ivec = 0 iunvn = 0 iuivi = 0 if ((map_tn_unvn.gt.0.and.ip.eq.ixt).or. + (map_ht_unvn.gt.0.and.ip.eq.ixz).or.ip.eq.ixunvn) then iunvn = 1 ivec = 2 ! vectors on top of contours if (ip.eq.ixunvn) ivec = 1 ! vectors only endif if ((map_ep_uivi.gt.0.and.ip.eq.ixepot).or.ip.eq.ixuivi) then iuivi = 1 ivec = 2 if (ip.eq.ixuivi) ivec = 1 endif c c Set plot options (see pltopt.h): c pltmin = cmin(ip) pltmax = cmax(ip) conint = cint(ip) scfac = scalefac(ip) pltchsize = 0. vlc = 0. ! vector low cutoff magnitude vhc = 0. ! vector height cutoff magnitude vsc = 0. ! vector scale arrow magnitude if (iunvn.gt.0) then vlc = cmin(ixunvn) vhc = cmax(ixunvn) vsc = cint(ixunvn) elseif (iuivi.gt.0) then vlc = cmin(ixuivi) vhc = cmax(ixuivi) vsc = cint(ixuivi) endif labelv = ivec_label c c zp/ht loop: c nzpht = 0 didgcmbot = .false. isgcmbot = .false. do k=1,mxzpht if (fmap_zpht(k).eq.spval) goto 50 nzpht = nzpht+1 if (ihtindep(ip).gt.0.and.nzpht.gt.1) goto 50 c c Determine which model, get global ave height, and global field: c zpht = fmap_zpht(k) if (zpht.lt.cplzp(1)) zpht = cplzp(1) if (zpht.gt.cplzp(ncplzp)) zpht = cplzp(ncplzp) ! fix later for hts if (ifcpl(ip).gt.0.and.zpht.le.cplzp(nccmlev)) then modstr4 = 'CCM ' modstr = ccmres else modstr4 = 'TGCM' modstr = model isgcmbot = .false. if (zpht.le.cplzp(nccmlev)) then zpht = gcmzp(1) isgcmbot = .true. endif endif if (isgcmbot.and.didgcmbot) goto 50 call getfldglb(modstr4,fglb0,imx,jmx,mxfproc,zpht,logplt(ixz), + ixz,fhist,kmx,nfhist,gcmzp,ixfhist,fccm,nccmlev,nccmfld, + cplzp,ixfccm) avglbht = fglbm(fglb0,imx,jmx,gcmlat,dlat,dlon,spval) if (ip.ne.ixunvn.and.ip.ne.ixuivi) + call getfldglb(modstr4,fglb0,imx,jmx,mxfproc,zpht, + logplt(ip),ip,fhist,kmx,nfhist,gcmzp,ixfhist,fccm,nccmlev, + nccmfld,cplzp,ixfccm) c c Get velocities if necessary: c if (ivec.gt.0) then if (iunvn.gt.0) then ! get un and vn call getfldglb(modstr4,uglb,imx,jmx,mxfproc,zpht, + logplt(ixu),ixu,fhist,kmx,nfhist,gcmzp,ixfhist,fccm, + nccmlev,nccmfld,cplzp,ixfccm) call getfldglb(modstr4,vglb,imx,jmx,mxfproc,zpht, + logplt(ixv),ixv,fhist,kmx,nfhist,gcmzp,ixfhist,fccm, + nccmlev,nccmfld,cplzp,ixfccm) elseif (iuivi.gt.0) then ! get ui and vi call getfldglb(modstr4,uglb,imx,jmx,mxfproc,zpht, + logplt(ixui),ixui,fhist,kmx,nfhist,gcmzp,ixfhist,fccm, + nccmlev,nccmfld,cplzp,ixfccm) call getfldglb(modstr4,vglb,imx,jmx,mxfproc,zpht, + logplt(ixvi),ixvi,fhist,kmx,nfhist,gcmzp,ixfhist,fccm, + nccmlev,nccmfld,cplzp,ixfccm) endif if (ivec.eq.1) then zmin = vecmin(uglb,vglb,imx*jmx,spval) zmax = vecmax(uglb,vglb,imx*jmx,spval) endif endif c c Save global ascii data: c c if (iwrdat.gt.0) then c if (ihtindep(ip).gt.0) then ! check this first c write(infolab,"('UT=',f4.1)") ut c elseif (fmap_zpht(k).le.gcmzp(kmx)) then c write(infolab,"('UT=',f4.1,' ZP=',f5.1,' AVE HT=',f7.2, c + ' (KM)')") ut,fmap_zpht(k),avglbht c else c write(infolab,"('UT=',f4.1,' HT=',f6.1,' (KM)')") ut, c + fmap_zpht(k) c endif c iframe_dat = iframe_dat+1 c call wrdat(iwrdat,ludat,fglb0,imx,jmx,gcmlon,gcmlat, c + 'LONGITUDE (DEG)','LATITUDE (DEG)',flab(ip),infolab, c + histlab,iframe_dat,'flxproc',senddat) c endif c if (iwrxdr.gt.0) then iframe_xdr = iframe_xdr+1 call clearstr(msgout) if (mkplt.le.0) then call fminmax(fglb0,imx*jmx,zmin,zmax,spval) endif call mkmaplab(zpht,avglbht,0,ivec,zmin,zmax, + ciu,blabs,tlabs,msgout,'CE',mtime,modstr,ip) call wrxdr(trim(flnm_xdr),fglb0,imx,jmx,gcmlon,gcmlat, + 'LONGITUDE (DEG)','LATITUDE (DEG)', + tlabs(2),tlabs(1),blabs(2),blabs(1),mtime,0) if (mkplt.le.0) write(6,"('XDR frame ',i4,': ',a)") + iframe_xdr,msgout endif c if (mkplt.le.0) goto 50 c c Make global CE map: c if (map_global.gt.0) then proj = 'CE' call setmultivp(vpce,iadvfr,nppf+1,multiplt,ipltrowcol,vp) if (map_global_cenlon.ne.ispval) then cenlon = float(map_global_cenlon) else cenlon = fslt(float(map_global_censlt),ut,dum,3) endif call pltce(fglb0,uglb,vglb,gcmlon,gcmlat,imx,jmx,ivec, + cenlon,map_continents,ut,vp,ispval,spval) call mkmaplab(zpht,avglbht,0,ivec,zmin,zmax, + ciu,blabs,tlabs,msgout,proj,mtime,modstr,ip) toffset = .06 boffset = .50 call wrlab6(tlabs,toffset,tlabchsz, + blabs,boffset,blabchsz,vp,0) nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + multiplt,iadvfr,nppf,msgout,'global maps',iframe_plt) endif c c Make polar stereographic maps: c if (map_polar.gt.0) then proj = 'ST' do ih=1,mxperimlat if (fmap_polar_perimlat(ih).ne.spval) then call setmultivp(vpst,iadvfr,nppf+1,multiplt,ipltrowcol, + vp) call pltst(fglb0,uglb,vglb,gcmlon,gcmlat,imx,jmx,ivec, + inc,fmap_polar_perimlat(ih),map_continents,ut,vp, + spval) call mkmaplab(zpht,avglbht,0,ivec,zmin,zmax, + ciu,blabs,tlabs,msgout,proj,mtime,modstr,ip) toffset = .06 boffset = .065 call wrlab6(tlabs,toffset,tlabchsz, + blabs,boffset,blabchsz,vp,0) call clearstr(infolab) write(infolab,"('PERIMLAT=',f6.2)") + fmap_polar_perimlat(ih) vpw = vp(2)-vp(1) xpos = vp(1)+.10*vpw ypos = vp(3)-.05*vpw chsz = .018 call wrlab(infolab(1:lenstr(infolab)),xpos,ypos,chsz, + vp,0) nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + multiplt,iadvfr,nppf,msgout,'polar maps',iframe_plt) endif enddo endif c c Make satellite view maps: c if (map_satview.gt.0) then proj = 'SV' call setmultivp(vpsv,iadvfr,nppf+1,multiplt,ipltrowcol,vp) call pltsv(fglb0,uglb,vglb,gcmlon,gcmlat,imx,jmx,ivec, + inc,fmap_satview_eradii,fmap_satview_latlon, + fmap_satview_latslt,map_continents,ut,vp,spval) call mkmaplab(zpht,avglbht,0,ivec,zmin,zmax, + ciu,blabs,tlabs,msgout,proj,mtime,modstr,ip) toffset = .03 boffset = .01 call wrlab6(tlabs,toffset,tlabchsz, + blabs,boffset,blabchsz,vp,0) nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + multiplt,iadvfr,nppf,msgout,'satview maps',iframe_plt) endif if (isgcmbot) didgcmbot = .true. 50 continue enddo ! k=1,mxzpht c c Increment field count: c nplt = nplt+1 if (nplt.gt.nfproc) then ! done with all fields goto 102 else goto 101 ! go back for next field endif 100 continue enddo ! ip=1,mxfproc 102 continue ! nplt = 1,nfproc c c If doing multiplt and multiadvfr > 0, advance frame if page is only c partially full (was advanced above if page is full) c if (multiplt.gt.0.and.multiadvfr.eq.1.and.iadvfr.le.0) + call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + multiplt,1,nppf,' ','map projections',iframe_plt) return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine mkmaplab(zpht,avht,iht,ivec,zmin,zmax,ciu,blabs,tlabs, + msgout,proj,mtime,modstr,ip) c c Make labels for maps (any projection): c include "flxproc.h" include "ccm.h" include "tgcmhdr.h" integer mtime(3) character*(*) blabs(3),tlabs(3) character*80 msgout character*2 proj character*16 modstr character*4 modstr4 character*8 mapstr,datestr c lmodstr = lenstr(modstr) c c-------------------------- top labels -------------------------- c c Top (third) of 3 top labels is optional 80-char annotation: c call clearstr(tlabs(3)) if (lenstr(map_top_anno).gt.0) then write(tlabs(3),"(a)") map_top_anno endif c c Middle (second) of 3 top labels is field name: c call clearstr(tlabs(2)) write(tlabs(2),"(a,' FIELD ',a)") modstr(1:lmodstr),flab8(ip) if (ivec.eq.2) then len = lenstr(tlabs(2)) if (ip.eq.ixt) + write(tlabs(2)(len+1:len+12),"(' (UN+VN M/S)')") if (ip.eq.ixepot) + write(tlabs(2)(len+1:len+12),"(' (UI+VI M/S)')") endif c c Bottom (first) of 3 top labels is grid info: c Use model day (nday) to get date (date(2) frm header is unreliable) c Also use ifix(ut) with i2.2 so label doesn't jump in animxdr animations) c iyd = 96000 + nday call iyd2date(iyd,imo,ida,iyr) write(datestr,"(i2.2,'/',i2.2)") imo,ida c call clearstr(tlabs(1)) c write(tlabs(1),"(a,' UT=',f4.1,' ZP=',f5.1,' AVE HT=', c + f7.2,' (KM)')") datestr(1:lenstr(datestr)),ut,zpht,avht write(tlabs(1),"(a,' UT=',i2.2,' ZP=',f5.1,' AVE HT=', + f7.2,' (KM)')") datestr(1:lenstr(datestr)),ifix(ut),zpht,avht c c-------------------------- bottom labels -------------------------- c c Top (third) of 3 bottom labels is min,max,etc: c call clearstr(blabs(3)) if (ivec.ne.1) then if (scalefac(ip).eq.1.) then write(blabs(3),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(blabs(3),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe8.2,')')") zmin,zmax,ciu,scalefac(ip) endif else write(blabs(3),"('VECTOR MIN,MAX=',2(1pe12.4))") zmin,zmax endif c c Middle (second) of 3 bottom labels is histvol and mtime: c call clearstr(blabs(2)) if (modstr.eq.model) then ! is from tgcm write(blabs(2),"(a,' ',a,' (DAY,HR,MIN=', + i3.3,',',i2.2,',',i2.2,')')") modstr(1:lmodstr), + histvols(ivol)(1:lenstr(histvols(ivol))),(mtime(i),i=1,3) else ! is from ccm write(blabs(2),"(a,' ',a,' (DAY ',f7.3,')')") + modstr(1:lmodstr),ccmlsd(1:lenstr(ccmlsd)),ccmday endif c c Bottom (1st) of 3 bottom labels is unused: c call clearstr(blabs(1)) c c Make message to be sent to stdout by mkmaps: c call clearstr(msgout) if (modstr.eq.model) then modstr4 = 'tgcm' else modstr4 = 'ccm2' endif mapstr = ' ' if (proj.eq.'CE') mapstr = 'glb map ' if (proj.eq.'ST') mapstr = 'pol map ' if (proj.eq.'SV') mapstr = 'satv map' if (ihtindep(ip).gt.0) then ! ht-independent field write(msgout,"(a,' ',a,' ',a,' min,max=', + 1pe9.2,',',1pe9.2)") modstr4,flab8(ip),mapstr,zmin,zmax elseif (iht.eq.0) then ! at zp write(msgout,"(a,' ',a,' ',a,' Zp=',f5.1, + ' min,max=',2(1pe9.2)))") modstr4,flab8(ip),mapstr,zpht, + zmin,zmax c else ! at ht c if (proj.eq.'CE') then c write(msgout,"('Frame ',i3,': ',a,' global map: Ht=',f6.1, c + ': min,max=',1pe10.3,',',1pe10.3)") c + iframe_plt+1,flab8(ip),zpht,zmin,zmax c elseif (proj.eq.'ST') then c write(msgout,"('Frame ',i3,': ',a,' polar map: Ht=',f6.1, c + ': min,max=',1pe10.3,',',1pe10.3)") c + iframe_plt+1,flab8(ip),zpht,zmin,zmax c elseif (proj.eq.'SV') then c write(msgout,"('Frame ',i3,': ',a,' satv map: Ht=',f6.1, c + ': min,max=',1pe10.3,',',1pe10.3)") c + iframe_plt+1,flab8(ip),zpht,zmin,zmax c endif endif return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine getfldglb(modstr,fglbout,imx,jmx,mxfproc,zpht,log,ip, + fgcm,ngcmzp,nfgcm,gcmzp,ixfgcm, + fccm,nccmzp,nfccm,ccmzp,ixfccm) c c Call getfld for either tgcm or ccm depending on modstr: c character*4 modstr real fglbout(imx,jmx) real fgcm(imx,ngcmzp,jmx,nfgcm),gcmzp(ngcmzp) real fccm(imx,nccmzp,jmx,nfccm),ccmzp(nccmzp) integer ixfgcm(mxfproc),ixfccm(mxfproc) c if (modstr.eq.'TGCM') then call getfld(fgcm,imx,ngcmzp,jmx,nfgcm,fglbout,imx,jmx, + 'GLOBAL ',dum,dum,gcmzp,zpht,log,ixfgcm,ip, + ier) else call getfld(fccm,imx,nccmzp,jmx,nfccm,fglbout,imx,jmx, + 'GLOBAL ',dum,dum,ccmzp,zpht,log,ixfccm,ip, + ier) endif return end