! module mk_maps use proc use hist,only: history use input implicit none ! ! Map projection structure: ! type mapproj character(len=2) :: proj ! projection (CE,ST,SV,MO) real :: zpht ! selected ht or zp real :: aveht ! ave hts if at zp (levtype=='zp') real :: perimlat ! perim lat for polar ST real :: cenlat,cenlon ! projection center real :: pltlat0,pltlat1 ! min,max latitudes to plot real :: eradii ! earth radii from satview proj character(len=5) :: levtype ! 'zp','ht','indep', or 'integ' integer :: icontinents ! flag for continental outlines integer :: nx,ny ! number of x,y points integer :: log10 ! if log10 > 0, log10 has been taken real,pointer :: xx(:)=>NULL() ! x coords (lon) real,pointer :: yy(:)=>NULL() ! y coords (lat) character(len=40) :: xlab,ylab ! x and y axis labels real,pointer :: data(:,:)=>NULL() ! the data array real,pointer :: u(:,:)=>NULL() ! vectors (s.a. vectype) real,pointer :: v(:,:)=>NULL() ! vectors (s.a. vectype) character(len=8) :: + sname, ! short field name + vectype ! 'UNVN','+UNVN','UIVI','+UIVI' character(len=56) :: fname ! full field name character(len=16) :: funits ! field units character(len=16) :: ftype ! field type character(len=240) :: tlabs(3),blabs(3) ! top and bottom labels character(len=8) :: difftype ! raw or percent logical :: known ! is field known to the processor end type mapproj ! ! inc is here so it can be used in mksltmaps: ! integer :: inc(nlat) ! real :: tlabchsz(3)=(/.02,.02,.02/), + blabchsz(3)=(/.018,.018,.018/) contains !----------------------------------------------------------------------- subroutine mkmaps(f,fcntr,nf,h,hcntr) ! ! Make map projections at current ut: ! use fields,only: field,oh_alt,printfield,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use plt implicit none ! ! Args: integer,intent(in) :: nf type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: ixf,izpht,ier,izp,ixz,nzpht,nzpht_sat,ih,iadvfr, + ixun,ixvn,ixui,ixvi,loginterp,k,nnans,ixhmf2,idim1,idim2,idim3 type(mapproj) :: map real :: dum=0.,fmin,fmax,hslice(nlon,nlat),hslice_cntr(nlon,nlat) real,allocatable :: fhts(:) ! heights of a height-only field character(len=80) :: msgout real :: vp(4) real :: vpce(4) = (/.100,.960,.450,.880/) real :: vpst(4) = (/.1625,.8375,.195,.870/) real :: vpsv(4) = (/.150,.850,.190,.890/) real :: vpmo(4) = (/.060,.960,.250,.800/) ! ! sattrack array (nlon, nlat, ntimes, nf, nzpht_sat) ! ntimes will be 25 for now (hrly for 1 day) real,allocatable,save :: satdat(:,:,:,:,:) real :: scratch(nlon,nlat,3) integer,save :: itime=0,iday=0 real,save :: utimes(25) integer :: i,j,kk real :: dntime,uptime,mapavg,frac integer,save :: precdays real,allocatable :: zglb(:,:,:), ! either z_ifaces or z_midpts | zglbcntr(:,:,:) ! ! Externals: integer,external :: ixfind,ixfind8,ixfindc,fseries real,external :: fslt,convlon,fglbm,fmean,gettime ! integer,parameter :: nlat5 = 36 integer :: inc5(nlat5) = | (/24, 8, 6, 4, 2, 2, 2, 2, 2, 2, | 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, | 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, | 2, 2, 4, 6, 8, 24/) ! ! Set inc for either 5.0 or 2.5 deg resolution: if (nlat==nlat5) then ! 5.0 deg res do j=1,nlat inc(j) = inc5(j) enddo elseif (nlat==nlat5*2) then ! 2.5 deg res do j=1,nlat5-1 inc(j*2-1) = inc5(j) ! 1,3,5,...,65,67,69 inc(j*2) = inc5(j) ! 2,4,6,...,66,68,70 enddo inc(nlat) = inc5(nlat5) inc(nlat-1) = inc5(nlat5) else write(6,"('mkmaps: unknown latitude resolution nlat=',i3)") | nlat stop 'mkmaps' endif ! if (isattrack <= 0.and.map_global <= 0.and.map_polar <= 0.and. | map_satview <= 0.and.map_mollweide <= 0) then write(6,"(/,'>>> mkmaps: returning because isattrack, ', | 'map_global, map_polar, map_satview, and map_mollweide ', | 'are all <= 0')") write(6,"(' To avoid this message, set ipltmaps=0')") return endif ! if(isattrack>0) then itime = itime + 1 utimes(itime) = h%mtime(2) + h%mtime(3)/60. if(itime==1) then iday = h%mtime(1) precdays = iday - ibaseday print *, ' SATTRACK: precessing orbit by ',precdays,' days.' nzpht_sat = fseries(fmap_zpht,mxzpht,-spval,spval) allocate(satdat(nlon,nlat,ntimes,nf,nzpht_sat),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating satdat") elseif(iday < h%mtime(1)) then utimes(itime) = utimes(itime) + 24. endif endif write(6,"(' ')") if (.not.diffs) then write(6,"('Map projections at ut = ',f5.2,' mtime = ',3i3)") + h%ut,h%mtime else write(6,"('Difference fields: ')",advance='NO') write(6,"('Map projections at ut = ',f5.2)") h%ut write(6,"(' perturbed mtime = ',3i4,' control mtime = ',3i4)") + h%mtime,hcntr%mtime endif if (multiadvfr > 0) nppf = 0 ! ! Locate index to heights field (needed for interp, etc): ixz = ixfindc(f%fname8,nf,'Z ') if (ixz <= 0) then write(6,"('>>> mkmaps: need heights in fields: field names', + ' = ',(8a8))") f%fname8 stop 'mkmaps z' endif ! ! Full lon grid on x-axis: map%nx = nlon if (.not.associated(map%xx)) then allocate(map%xx(map%nx),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%xx") endif map%xx = gcmlon write(map%xlab,"('LONGITUDE (DEG)')") ! ! Full lat grid on y-axis: ! map%ny = nlat ! if (.not.associated(map%yy)) then ! allocate(map%yy(map%ny),stat=ier) ! if (ier /= 0) ! + call allocerr(ier,"mkmaps allocating map%yy") ! endif ! map%yy = gcmlat ! ! Latitude min,max to plot on y-axis: ! fmap_latminmax(2) is from namelist read in input.F ! These are applied in the conpack calls in sub pltmaps. ! if (fmap_latminmax(1)==spval) then map%pltlat0 = gcmlat(1) else map%pltlat0 = fmap_latminmax(1) endif if (fmap_latminmax(2)==spval) then map%pltlat1 = gcmlat(nlat) else map%pltlat1 = fmap_latminmax(2) endif if (map%pltlat0 < gcmlat(1).or. | map%pltlat0 > gcmlat(nlat)) then write(6,"('>>> mkmaps: bad map%pltlat0=',f9.3)") | map%pltlat0 write(6,"('Please check namelist read fmap_latminmax(1)')") stop 'map%pltlat0' endif if (map%pltlat1 < gcmlat(1).or. | map%pltlat1 > gcmlat(nlat)) then write(6,"('>>> mkmaps: bad map%pltlat1=',f9.3)") | map%pltlat1 write(6,"('Please check namelist read fmap_latminmax(2)')") stop 'map%pltlat1' endif if (map%pltlat0 >= map%pltlat1) then write(6,"('>>> mkmaps: pltlat0 must be < pltlat1: ', | ' pltlat0=',f9.3,' pltlat1=',f9.3)") map%pltlat0, | map%pltlat1 write(6,"('Please check namelist read fmap_latminmax')") stop 'map%pltlat' endif ! write(6,"('mkmaps: map%pltlat0,1=',2f9.3)") map%pltlat0, ! | map%pltlat1 map%ny = nlat if (.not.associated(map%yy)) then allocate(map%yy(map%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%yy") endif map%yy = gcmlat write(map%ylab,"('LATITUDE (DEG)')") map%icontinents = map_continents ! ! Remove invalid fmap_zpht (below bottom zp boundary): ! (values above top zp boundary are assumed to be heights) ! do izpht=1,mxzpht ! if (fmap_zpht(izpht)/=spval.and.fmap_zpht(izpht)>> fmap_zpht requested hmf2 hts, but ', | 'could not find field HMF2, which is needed for', | ' interpolation..')") stop 'hmf2' endif write(6,"('Will interpolate to hmf2 hts: HMF2 global ', | 'min,max=',2f8.2)") minval(f(ixhmf2)%data(:,:,1)), | maxval(f(ixhmf2)%data(:,:,1)) endif enddo ! ! Allocate the map field data: if (.not.associated(map%data)) then allocate(map%data(map%nx,map%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%data") endif ! ! Allocate the map field data: if (.not.associated(map%data)) then allocate(map%data(map%nx,map%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%data") endif ! ! Field loop. Field must be requested, data pointer must be associated. fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data)) + cycle fields_loop map%fname = f(ixf)%fname56 map%sname = f(ixf)%fname8 map%ftype = f(ixf)%type map%funits = f(ixf)%units map%known = f(ixf)%known map%difftype = f(ixf)%difftype loginterp = 0 if (trim(f(ixf)%type)=='DENSITY') loginterp = 1 ! ! If ht-only field, construct array of heights: if (trim(f(ixf)%vtype)=='HEIGHT') then nzpht = f(ixf)%nlev allocate(fhts(nzpht),stat=ier) if (ier /= 0) + call allocerr(ier,"gethtfld allocating fhts") fhts(:) = f(ixf)%lev(:) else ! ! Set zglb according to current field's zptype: if (allocated(zglb)) deallocate(zglb) idim1 = size(f(ixf)%data(:,1,1)) idim2 = size(f(ixf)%data(1,:,1)) idim3 = size(f(ixf)%data(1,1,:)) allocate(zglb(idim1,idim2,idim3),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating zglb") if (trim(f(ixf)%zptype) == 'INTERFACES') then zglb = z_ifaces elseif (trim(f(ixf)%zptype) == 'MIDPOINTS') then zglb = z_midpts endif if (diffs) then if (allocated(zglbcntr)) deallocate(zglbcntr) idim1 = size(fcntr(ixf)%data(:,1,1)) idim2 = size(fcntr(ixf)%data(1,:,1)) idim3 = size(fcntr(ixf)%data(1,1,:)) allocate(zglbcntr(idim1,idim2,idim3),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating zglbcntr") if (trim(f(ixf)%zptype) == 'INTERFACES') then zglbcntr = zcntr_ifaces elseif (trim(f(ixf)%zptype) == 'MIDPOINTS') then zglbcntr = zcntr_midpts endif endif endif ! ! Set field min, max, and interval to be used in contouring: ! (default f(i)%cmin,cmax,cint is 0.,0.,0., but may be provided ! by user via namelist input fmnmxint) ! (pltmin,pltmax,conint are in plt.f) pltmin = f(ixf)%cmin pltmax = f(ixf)%cmax conint = f(ixf)%cint scfac = f(ixf)%scalefac ! ! Initialize vector low cutoff, high cutoff, scale, and arrow length: vlc = 0. ! vector low cutoff magnitude vhc = 0. ! vector high cutoff magnitude vsc = 0. ! vector scale arrow magnitude vrl = 0. ! vector realized length if (vmag_len>0.) vrl = vmag_len labelv = ivec_label ! ! Prepare for velocity vectors arrow plots: ! map%vectype = ' ' -> not plotting vectors ! map%vectype = 'UNVN', or '+UNVN' -> unvn vectors only, or vecs+contour ! map%vectype = 'UIVI', or '+UIVI' -> uivi vectors only, or vecs+contour ! map%vectype = ' ' if (f(ixf)%fname8=='UNVN '.or. + (f(ixf)%fname8=='TN '.and.map_tn_unvn > 0).or. + (f(ixf)%fname8=='Z '.and.map_ht_unvn > 0)) then ixun = ixfindc(f%fname8,nf,'UN ') ixvn = ixfindc(f%fname8,nf,'VN ') if (ixun>0.and.ixvn>0) then if (associated(f(ixun)%data).and.associated(f(ixvn)%data)) + then map%vectype = '+UNVN ' if (f(ixf)%fname8=='UNVN ') map%vectype = 'UNVN ' if (f(ixf)%fname8/='UNVN ') + map%sname = trim(f(ixf)%fname8)//trim(map%vectype) vlc = vn_scale(1) vhc = vn_scale(2) vsc = vn_scale(3) ! write(6,"('Vector low,high,scale for UN+VN=',3f10.2)") ! + vlc,vhc,vsc else write(6,"('>>> mkmaps: Need UN and VN to plot neutral ', + 'wind vectors on maps.', + /' (vectors will not be plotted).')") endif else write(6,"('>>> mkmaps: Need UN and VN to plot neutral ', + 'wind vectors on maps.', + /' (vectors will not be plotted).')") endif endif if (f(ixf)%fname8=='UIVI '.or. + (f(ixf)%fname8=='POTEN '.and.map_ep_uivi > 0)) then ixui = ixfindc(f%fname8,nf,'UI ') ixvi = ixfindc(f%fname8,nf,'VI ') if (ixui>0.and.ixvi>0) then if (associated(f(ixui)%data).and.associated(f(ixvi)%data)) + then map%vectype = '+UIVI ' if (f(ixf)%fname8=='UIVI ') map%vectype = 'UIVI ' if (f(ixf)%fname8/='UIVI ') map%sname = 'EP+UIVI ' vlc = vi_scale(1) vhc = vi_scale(2) vsc = vi_scale(3) ! write(6,"('Vector low,high,scale for UI+VI=',3f10.2)") ! + vlc,vhc,vsc else write(6,"('>>> mkmaps: Need UI and VI to plot ion drifts', + ' on maps.',/' (vectors will not be plotted).')") endif else write(6,"('>>> mkmaps: Need UI and VI to plot ion drifts', + ' on maps.',/' (vectors will not be plotted).')") endif endif if (len_trim(map%vectype) > 0) then if (associated(map%u)) deallocate(map%u) allocate(map%u(map%nx,map%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%u") if (associated(map%v)) deallocate(map%v) allocate(map%v(map%nx,map%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkmaps allocating map%v") endif 100 continue ! jump here from end of fields_loop for ht-integ nzpht = mxzpht if (map%levtype == 'integ') nzpht = 1 ! ht-integ (loop once only) ! ! Selected height/pressure loop: ! zpht_loop: do izpht=1,nzpht ! ! Determine levtype from field vtype and requested zpht ! (levtype integ is set at end of fields_loop). ! Field types (f(i)%vtype) may be 'ZP','HEIGHT', or 'HT-INDEP' ! (ZP types may be plotted at selected zp, or interpolated to ht) ! Possible values for map%levtype: ! map%levtype == 'zp ' -> at selected Zp's ! map%levtype == 'ht ' -> at selected Ht's (at hmf2 if fmap_zpht==hmf2flag) ! map%levtype == 'indep' -> ht-independent (nlev==1) (fof2,hmf2) ! map%levtype == 'integ' -> ht-integrated (emission or oh-band) ! if (map%levtype /= 'integ') then if (fmap_zpht(izpht)==spval) cycle zpht_loop if (trim(f(ixf)%vtype)=='HT-INDEP') then map%levtype = 'indep' izp = 1 elseif (trim(f(ixf)%vtype)=='ZP') then ! plot ZP var at hmf2 hts if (fmap_zpht(izpht) == hmf2flag) then map%levtype = 'ht ' map%zpht = hmf2flag elseif (fmap_zpht(izpht) > f(ixf)%lev(f(ixf)%nlev)+.001) | then map%levtype = 'ht ' map%zpht = fmap_zpht(izpht) else ! levtype is zp map%levtype = 'zp ' izp=ixfind8(f(ixf)%lev(:),f(ixf)%nlev,fmap_zpht(izpht), | f(ixf)%dlev) if (izp <= 0) then if (fmap_zpht(izpht) < f(ixf)%lev(1)) then write(6,"('>>> mkmaps: fmap_zpht ',f9.3,' is below', | ' bottom boundary for field ',a)") | fmap_zpht(izpht),f(ixf)%fname8 write(6,"('Will plot bottom boundary at zp ',f9.3)") | f(ixf)%lev(1) izp = 1 else write(6,"('>>> mkmaps: could not find zp near ', | f10.3,' for field ',a)") fmap_zpht(izpht), | f(ixf)%fname8 cycle zpht_loop endif endif map%zpht = f(ixf)%lev(izp) ! izp was determined above endif ! zp or ht elseif (trim(f(ixf)%vtype)=='HEIGHT') then map%levtype = 'ht ' if (fmap_zpht(izpht) >= f(ixf)%lev(1).and. | fmap_zpht(izpht) <= f(ixf)%lev(f(ixf)%nlev)) then map%zpht = fmap_zpht(izpht) else ! write(6,"('Note: Height-only field ',a,' not available', ! | ' at fmap_zpht=',f9.2)")f(ixf)%fname8,fmap_zpht(izpht) cycle zpht_loop endif endif ! vtype endif ! not integ ! ! At selected zp or is ht-independent: if (map%levtype=='zp '.or.map%levtype=='indep') then ! ! Define field map%data(nlon,nlat) from f(ixf)%data(nlon,nlat,npress): if (.not.diffs) then map%data(:,:) = f(ixf)%data(:,:,izp) else call mkdiffs(f(ixf)%data(:,:,izp), + fcntr(ixf)%data(:,:,izp),map%data,nlon*nlat, + f(ixf)%difftype) endif ! ! Calculate ave height for info label: ! ! 1/9/98: transferring heights to local array before calling fglbm ! is necessary to avoid bug in sgi (passing f(ixz)%data(:,:,izp) ! as 1st arg to fglbm results in operand range error w/ bad ixz ! subscript) ! 1/21/08: Use geopotential zglb set above according to field's ! zptype (interfaces or midpoints): ! hslice = zglb(:,:,izp) if (.not.diffs) then map%aveht = fglbm(hslice,nlon,nlat,gcmlat,dlat,dlon,spval) else hslice_cntr = zglbcntr(:,:,izp) map%aveht = 0.5* | (fglbm(hslice,nlon,nlat,gcmlat,dlat,dlon,spval)+ | fglbm(hslice_cntr,nlon,nlat,gcmlat,dlat,dlon,spval)) endif ! ! Get vectors at pressure surface if necessary: if (index(map%vectype,'UN') > 0) then ! get un and vn if (.not.diffs) then map%u(:,:) = f(ixun)%data(:,:,izp) map%v(:,:) = f(ixvn)%data(:,:,izp) else call mkdiffs(f(ixun)%data(:,:,izp), + fcntr(ixun)%data(:,:,izp),map%u,nlon*nlat, + f(ixun)%difftype) call mkdiffs(f(ixvn)%data(:,:,izp), + fcntr(ixvn)%data(:,:,izp),map%v,nlon*nlat, + f(ixvn)%difftype) endif elseif (index(map%vectype,'UI') > 0) then ! get ui and vi if (.not.diffs) then map%u(:,:) = f(ixui)%data(:,:,izp) map%v(:,:) = f(ixvi)%data(:,:,izp) else call mkdiffs(f(ixui)%data(:,:,izp), + fcntr(ixui)%data(:,:,izp),map%u,nlon*nlat, + f(ixui)%difftype) call mkdiffs(f(ixvi)%data(:,:,izp), + fcntr(ixvi)%data(:,:,izp),map%v,nlon*nlat, + f(ixvi)%difftype) endif endif ! ! Interpolate to constant height surface: ! (take diffs *after* interpolation) ! (if ht-only field, call gethtfld to transfer to map%data) ! elseif (map%levtype == 'ht ') then if (.not.diffs) then if (trim(f(ixf)%vtype)/='HEIGHT') then if (map%zpht /= hmf2flag) then call glbhtin(f(ixf)%data,zglb,nlon,nlat, + npress,map%data,map%zpht,1,loginterp,spval,ier,0) else ! ! Global interpolation to hmf2 heights: ! subroutine cuthtint(fin,fht,idim1,nzp,fout,hts,nhts,logint,spval,ier,iprnt) ! do j=1,nlat do i=1,nlon call cuthtint(f(ixf)%data(i,j,:), | zglb(i,j,:),1,npress, | map%data(i,j),f(ixhmf2)%data(i,j,1),1, | loginterp,spval,ier,0) enddo enddo endif else ! ht-only field call gethtfld(f(ixf),fhts,ier) ! local to mkmaps if (ier /= 0) cycle zpht_loop endif else if (trim(f(ixf)%vtype)/='HEIGHT') then call glbhtin(f(ixf)%data, ! interpolate perturbed + zglb,nlon,nlat,npress,hslice,map%zpht,1, + loginterp,spval,ier,0) call glbhtin(fcntr(ixf)%data, ! interpolate control + zglbcntr,nlon,nlat,npress,hslice_cntr, + map%zpht,1,loginterp,spval,ier,0) else ! diffs of ht-only field (interpolation not necessary) call gethtfld(f(ixf),fhts,ier) if (ier /= 0) cycle zpht_loop hslice = map%data call gethtfld(fcntr(ixf),fhts,ier) if (ier /= 0) cycle zpht_loop hslice_cntr = map%data endif call mkdiffs(hslice,hslice_cntr,map%data,nlon*nlat, + f(ixf)%difftype) endif ! ! Get vectors at height surface if necessary: if (index(map%vectype,'UN') > 0) then ! get un and vn if (.not.diffs) then call glbhtin(f(ixun)%data,zglb,nlon,nlat, + npress,map%u,map%zpht,1,0,spval,ier,0) call glbhtin(f(ixvn)%data,zglb,nlon,nlat, + npress,map%v,map%zpht,1,0,spval,ier,0) else ! un diffs at height: call glbhtin(f(ixun)%data,zglb,nlon,nlat,npress, + hslice,map%zpht,1,0,spval,ier,0) call glbhtin(fcntr(ixun)%data,zglbcntr,nlon, + nlat,npress,hslice_cntr,map%zpht,1,0,spval,ier,0) call mkdiffs(hslice,hslice_cntr,map%u,nlon*nlat, + f(ixun)%difftype) ! vn diffs at height: call glbhtin(f(ixvn)%data,zglb,nlon,nlat,npress, + hslice,map%zpht,1,0,spval,ier,0) call glbhtin(fcntr(ixvn)%data,zglbcntr,nlon, + nlat,npress,hslice_cntr,map%zpht,1,0,spval,ier,0) call mkdiffs(hslice,hslice_cntr,map%v,nlon*nlat, + f(ixvn)%difftype) endif elseif (index(map%vectype,'UI') > 0) then ! get ui and vi if (.not.diffs) then call glbhtin(f(ixui)%data,zglb,nlon,nlat, + npress,map%u,map%zpht,1,0,spval,ier,0) call glbhtin(f(ixvi)%data,zglb,nlon,nlat, + npress,map%v,map%zpht,1,0,spval,ier,0) else ! ui diffs at height: call glbhtin(f(ixui)%data,zglb,nlon,nlat,npress, + hslice,map%zpht,1,0,spval,ier,0) call glbhtin(fcntr(ixui)%data,zglbcntr,nlon, + nlat,npress,hslice_cntr,map%zpht,1,0,spval,ier,0) call mkdiffs(hslice,hslice_cntr,map%u,nlon*nlat, + f(ixui)%difftype) ! vi diffs at height: call glbhtin(f(ixvi)%data,zglb,nlon,nlat,npress, + hslice,map%zpht,1,0,spval,ier,0) call glbhtin(fcntr(ixvi)%data,zglbcntr,nlon, + nlat,npress,hslice_cntr,map%zpht,1,0,spval,ier,0) call mkdiffs(hslice,hslice_cntr,map%v,nlon*nlat, + f(ixvi)%difftype) endif endif map%aveht=0. ! ! Integrate in vertical (currently only emission fields, including OH): else ! ht-integrate if (.not.diffs) then if (trim(f(ixf)%vtype)/='HEIGHT') then call glb_integ(f(ixf)%data,zglb,map%data, + nlon,npress,nlat) ! ! 3/16/06 btf: There was a problem here where e5577 was INF at k==1,2 ! and lat >= 18, and elsewhere (194 total points). This was causing ! ncarg to hang in cpclam (called from contour in plt.F), even after ! replacing INF with spval with check_nans. I am leaving this call to ! check_nans with fatal stop in case this happens again. It was not ! obvious why mke5577 was making INF's, and also not obvious why cpclam ! was hanging even after replacing INF's with spval here..This was ! occurring only in a tiegcm-hist-lbc early development history at ! 80,0,0 ($TGCMDATA/TGCM.tiegcm_zbot-10_sres.p001.nc). This was with ! new lbc at zp=-10, so it was probably a lower boundary problem. ! call check_nans(f(ixf)%data,nlon,npress,nlat, | 'EMISS-INTEG',nnans,1,spval,1,1) ! last arg is ifatal ! else ! ht-only field (integ_htfld is local to mkmaps) call integ_htfld(f(ixf),fhts,nlon,nlat) endif else if (trim(f(ixf)%vtype)/='HEIGHT') then call glb_integ(f(ixf)%data,zglb,hslice,nlon, + npress,nlat) call glb_integ(fcntr(ixf)%data,zglbcntr, + hslice_cntr,nlon,npress,nlat) else ! ht-only field call integ_htfld(f(ixf),fhts,nlon,nlat) hslice = map%data call integ_htfld(fcntr(ixf),fhts,nlon,nlat) hslice_cntr = map%data endif call mkdiffs(hslice,hslice_cntr,map%data,nlon*nlat, + f(ixf)%difftype) endif if (trim(f(ixf)%type)=='EMISSION') then map%data = map%data * 1.e-6 ! rayleighs map%funits = "RAYLEIGHS" endif if (trim(f(ixf)%type)=='OH-BAND') then if (ibohv_watts <= 0) then map%data = map%data * 1.e-9 ! kilo-rayleighs map%funits = "KILO-RAYLEIGHS" else map%data = map%data / (4.*pi) ! watts/cm2-str map%funits = "WATTS/CM2-STR" endif endif map%aveht=0. endif !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Take log10 of field if requested: ! map_log10 is user input w/ default 0. If map_log10==1 take log ! of appropriate fields (mainly densities), if map_log10 > 1, ! take log10 of all fields. ! map%log10 = 0 if ((map_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + map_log10 > 1) then call log10f(map%data,map%nx*map%ny,1.e-20,spval) map%log10 = 1 endif !----------------------------------------------------------------------- ! ! Save data if doing sattrack run if(isattrack>0) then do i=1,nlon do j=1,nlat satdat(i,j,itime,ixf,izpht) = map%data(i,j) enddo ! j enddo ! i endif ! ! Make cyl eq global map: if (map_global > 0) then map%proj = 'CE' if (map_global_cenlon /= ispval) then map%cenlon = float(map_global_cenlon) elseif (map_global_censlt /= ispval) then map%cenlon = fslt(float(map_global_censlt),h%ut,dum,3) else map%cenlon = 0. endif if (iplot > 0) then ! ! Set up viewport (for multiplt): call setmultivp(vpce,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Contour: call pltce(map,vp,spval) ! ! Add local time axis at bottom: call sltxax(convlon(map%cenlon-180.,180), + convlon(map%cenlon+180.,180),h%ut,.25) ! ! Add top and bottom labels: call mkmaplabs(map,h,hcntr,msgout) call wrlab6(map%tlabs,.06,tlabchsz, + map%blabs,.55,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'global maps',iframe) endif endif ! map_global !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Make polar stereographic map: if (map_polar > 0) then map%proj = 'ST' ! ! Make map for each requested perimeter latitude: perimlat_loop: do ih=1,mxperimlat if (fmap_polar_perimlat(ih) == spval) cycle perimlat_loop map%perimlat = fmap_polar_perimlat(ih) if (iplot > 0) then ! ! Set up viewport (for multiplt): call setmultivp(vpst,iadvfr,nppf+1,multiplt,ipltrowcol, + vp) ! ! Make projection and contours: call pltst(map,vp,h%ut,gcmlat,inc,nlat,spval) ! ! Add top and bottom labels: call mkmaplabs(map,h,hcntr,msgout) call wrlab6(map%tlabs,.06,tlabchsz, + map%blabs,.10,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'polar maps',iframe) endif enddo perimlat_loop endif !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Make satellite view map: if (map_satview > 0) then map%proj = 'SV' ! ! Establish projection center lat,lon: if (fmap_satview_latlon(1)/=spval.and. + fmap_satview_latlon(2)/=spval) then map%cenlat = fmap_satview_latlon(1) map%cenlon = fmap_satview_latlon(2) elseif (fmap_satview_latslt(1)/=spval.and. + fmap_satview_latslt(2)/=spval) then map%cenlat = fmap_satview_latslt(1) map%cenlon = fslt(fmap_satview_latslt(2),h%ut,dum,3) else map%cenlat = 40. map%cenlon = 0. endif map%eradii = fmap_satview_eradii ! ! Set up viewport (for multiplt): if (iplot > 0) then call setmultivp(vpsv,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Make projection and contours: call pltsv(map,vp,inc,nlat,spval) ! ! Add top and bottom labels: call mkmaplabs(map,h,hcntr,msgout) call wrlab6(map%tlabs,.05,tlabchsz, + map%blabs,.10,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'satview maps',iframe) endif endif !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Make mollweide map: if (map_mollweide > 0) then map%proj = 'MO' ! ! Establish projection center lat,lon: if (fmap_mollweide_latlon(1)/=spval.and. + fmap_mollweide_latlon(2)/=spval) then map%cenlat = fmap_mollweide_latlon(1) map%cenlon = fmap_mollweide_latlon(2) elseif (fmap_mollweide_latslt(1)/=spval.and. + fmap_mollweide_latslt(2)/=spval) then map%cenlat = fmap_mollweide_latslt(1) map%cenlon = fslt(fmap_mollweide_latslt(2),h%ut,dum,3) else map%cenlat = 0. map%cenlon = 0. endif if (iplot > 0) then ! ! Set up viewport (for multiplt): call setmultivp(vpmo,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Make projection and contours: call pltmo(map,vp,spval) ! ! Add top and bottom labels: call mkmaplabs(map,h,hcntr,msgout) call wrlab6(map%tlabs,.05,tlabchsz, + map%blabs,.10,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'mollweide maps',iframe) endif endif !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Write output data files: if (iplot <= 0)call fminmax(map%data,size(map%data),zmin,zmax) if(isattrack==0.or.iwrsat==2) call wrout_glb(map,h,hcntr) if (map%levtype=='indep') exit zpht_loop ! ! make sattrack plots if(isattrack>0 .and. itime==ntimes) then ! ixf = 1 do i=1,nlon do j=1,nlat uptime = gettime(gcmlon(i),gcmlat(j),precdays,1) if(uptime < utimes(1)) uptime=uptime+24. if(uptime /= spval) then kk = ixfind(utimes,ntimes,uptime,1.) if(kk > 1) then frac = (uptime-utimes(kk-1))/(utimes(kk)-utimes(kk-1)) scratch(i,j,1) = satdat(i,j,kk-1,ixf,izpht) + frac* + (satdat(i,j,kk,ixf,izpht)-satdat(i,j,kk-1,ixf,izpht)) else frac = (uptime-utimes(kk))/(utimes(kk+1)-utimes(kk)) scratch(i,j,1) = satdat(i,j,kk,ixf,izpht) + frac* + (satdat(i,j,kk+1,ixf,izpht)-satdat(i,j,kk,ixf,izpht)) endif else scratch(i,j,1) = spval endif ! dntime = gettime(gcmlon(i),gcmlat(j),precdays,2) if(dntime < utimes(1)) dntime=dntime+24. if(dntime /= spval) then kk = ixfind(utimes,ntimes,dntime,1.) if(kk > 1) then frac = (dntime-utimes(kk-1))/(utimes(kk)-utimes(kk-1)) scratch(i,j,2) = satdat(i,j,kk-1,ixf,izpht) + frac* + (satdat(i,j,kk,ixf,izpht)-satdat(i,j,kk-1,ixf,izpht)) else frac = (dntime-utimes(kk))/(utimes(kk+1)-utimes(kk)) scratch(i,j,2) = satdat(i,j,kk,ixf,izpht) + frac* + (satdat(i,j,kk+1,ixf,izpht)-satdat(i,j,kk,ixf,izpht)) endif else scratch(i,j,2) = spval endif ! if(scratch(i,j,1)/=spval.and.scratch(i,j,2)/=spval) then scratch(i,j,3) = 0.5*(scratch(i,j,1)+scratch(i,j,2)) else scratch(i,j,3) = spval endif enddo enddo ! plot_loop: do isatplt=1,3 do kk=1,3 do j=1,nlat mapavg = fmean(scratch(:,j,kk),nlon,spval,1) do i=1,nlon if(idetrend/=0) then ! Detrend the data, if desired if(scratch(i,j,kk)/=spval) + scratch(i,j,kk) = scratch(i,j,kk) - mapavg endif if(kk==isatplt) map%data(i,j) = scratch(i,j,kk) enddo enddo enddo ! ! Make cyl eq global map: ! if (map_global > 0) then map%proj = 'CE' if (map_global_cenlon /= ispval) then map%cenlon = float(map_global_cenlon) elseif (map_global_censlt /= ispval) then map%cenlon = fslt(float(map_global_censlt),h%ut,dum,3) else map%cenlon = 0. endif if (iplot > 0) then ! ! Set up viewport (for multiplt): call setmultivp(vpce,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Contour: call pltce(map,vp,spval) ! ! Add top and bottom labels: call mkmaplabs(map,h,hcntr,msgout) call wrlab6(map%tlabs,.06,tlabchsz, + map%blabs,.55,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'global maps',iframe) endif ! endif ! map_global ! ! Make polar stereographic map: ! if (map_polar > 0) then !! map%proj = 'ST' ! ! Make map for each requested perimeter latitude: !! perimlat_loop2: do ih=1,mxperimlat !! if (fmap_polar_perimlat(ih) == spval) cycle perimlat_loop2 !! map%perimlat = fmap_polar_perimlat(ih) !! if (iplot > 0) then ! ! Set up viewport (for multiplt): !! call setmultivp(vpst,iadvfr,nppf+1,multiplt,ipltrowcol, !! + vp) ! ! Make projection and contours: !! call pltst(map,vp,h%ut,gcmlat,inc,nlat,spval) ! ! Add top and bottom labels: !! call mkmaplabs(map,h,hcntr,msgout) !! call wrlab6(map%tlabs,.06,tlabchsz, !! + map%blabs,.10,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: !! nppf = nppf+1 !! call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, !! + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, !! + 'polar maps',iframe) !! endif !! enddo perimlat_loop2 ! endif ! ! Write output data files: if (iplot <= 0)call fminmax(map%data,size(map%data),zmin,zmax) call wrout_glb(map,h,hcntr) enddo plot_loop endif enddo zpht_loop ! ! Do ht-integration (integht is user input w/ default 1): if (integht > 0 .and. (trim(f(ixf)%type)=='EMISSION'.or. + trim(f(ixf)%type)=='OH-BAND')) then if (map%levtype /= 'integ') then map%levtype = 'integ' goto 100 else ! just did ht-integ map%levtype = 'zp ' endif endif if (allocated(fhts)) deallocate(fhts) enddo fields_loop ! ! If doing multiplt and multiadvfr > 0, advance frame if page is only ! partially full (was advanced above if page is full) if (multiplt > 0 .and. multiadvfr == 1 .and. iadvfr <= 0.and. + iplot > 0) + call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,1,nppf,' ','maps',iframe) ! return contains !- - - - - - - - - - - - - - - internal to mkmaps - - - - - - - - - - - - subroutine gethtfld(ff,fhts,ier) ! ! Transfer a height-only field from ff (a field struct) to map%data at ! height surface map%zpht. ! This routine called only when the current field ff is height-only ! (i.e., OHBxx or OHVxx), and map%levtype=='ht', i.e., map%zpht ! is a height, not a pressure. ! ! Args: type(field),intent(in) :: ff real,intent(in) :: fhts(ff%nlev) integer,intent(out) :: ier ! ! Locals: integer :: kk ! ! Externals: integer,external :: ixfind ! ier = 0 ! ! Find field height closest to map%zpht: kk = ixfind(fhts,ff%nlev,map%zpht,ff%dlev) ! ! Desired height map%zpht is outside height range of the field: if (kk<=0) then write(6,"('Note: ',a,' not available at height ', + f10.2,/,6x,'This field is available only from ', + f10.2,' to ',f10.2,' km')") ff%fname8,map%zpht, + fhts(1),fhts(ff%nlev) ier = 1 ! ! Transfer to map%data: else map%zpht=fhts(kk) map%data = ff%data(:,:,kk) endif end subroutine gethtfld !- - - - - - - - - - - - - - - internal to mkmaps - - - - - - - - - - - - subroutine integ_htfld(ff,fhts,nlon,nlat) ! ! Args: type(field),intent(in) :: ff real,intent(in) :: fhts(ff%nlev) integer,intent(in) :: nlon,nlat ! ! Locals: integer :: i,j ! ! Externals: real,external :: quadrat ! do j=1,nlat do i=1,nlon map%data(i,j) = quadrat(ff%data(i,j,:),fhts,ff%nlev) enddo enddo end subroutine integ_htfld end subroutine mkmaps !----------------------------------------------------------------------- subroutine mkmaplabs(map,h,hcntr,msgout) use plt,only: zmin,zmax,ciu,vmin,vmax use input,only: mtimes ! ! Define map%tlabs and map%blabs: ! ! Args: type(mapproj),intent(inout) :: map type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Locals: character(len=16) :: units character(len=8 ) :: type integer :: i,l,len0,len1 ! ! Top (1st) top label is user optional: map%tlabs(1) = ' ' ! ! Middle (2nd) top label is full field name + units: ! map%tlabs(2) = ' ' map%tlabs(2) = trim(map%fname) if (map%known) then if (len_trim(map%fname)+len_trim(map%funits)+3 <= | len(map%tlabs(2))) | map%tlabs(2) = trim(map%fname)//' ('//trim(map%funits)//')' ! ! btf 6/23/04: if field is unknown, add sname to front of label. else if (trim(map%sname) /= trim(map%fname)) then if (len_trim(map%sname)+len_trim(map%fname)+ | len_trim(map%funits)+3+2 <= len(map%tlabs(2))) | map%tlabs(2) = trim(map%sname)//': '//trim(map%fname)// | ' ('//trim(map%funits)//')' else if (len_trim(map%fname)+len_trim(map%funits)+3 <= | len(map%tlabs(2))) | map%tlabs(2) = trim(map%fname)// | ' ('//trim(map%funits)//')' endif endif ! ! Add 'LOG10' at front of label if map%log10 > 0 if (map%log10 > 0) then if (len_trim(map%tlabs(2))+6 <= len(map%tlabs(2))) then map%tlabs(2) = 'LOG10 '//trim(map%tlabs(2)) else map%tlabs(2) = 'LOG10 '//trim(map%fname) endif endif ! ! Add 'DIFFS' at front of label if doing diffs: if (diffs) then if (trim(map%difftype)=='RAW') then if (len_trim(map%tlabs(2))+7 <= len(map%tlabs(2))) then map%tlabs(2) = 'DIFFS: '//trim(map%tlabs(2)) else map%tlabs(2) = 'DIFFS: '//trim(map%fname) endif else ! percent if (len_trim(map%tlabs(2))+15 <= len(map%tlabs(2))) then map%tlabs(2) = trim(map%difftype)//' DIFFS: '// | trim(map%tlabs(2)) else map%tlabs(2) = '% DIFFS: '//trim(map%fname) endif endif endif ! ! Bottom (3rd) top label is grid info: map%tlabs(3) = ' ' if (map%levtype == 'zp ') then ! at selected zp write(map%tlabs(3),"('UT=',f5.2,' ZP=',f7.3,' AVE HT=',f6.1)") + h%ut,map%zpht,map%aveht elseif (map%levtype == 'ht ') then ! at selected height if (map%zpht /= hmf2flag) then write(map%tlabs(3),"('UT=',f5.2,' HEIGHT=',f6.2,' (KM)')") + h%ut,map%zpht else write(map%tlabs(3),"('UT=',f5.2,' HMF2 HEIGHTS (KM)')") + h%ut endif elseif (map%levtype == 'indep') then ! ht-independent write(map%tlabs(3),"('UT=',f5.2)") h%ut else ! ht-integration write(map%tlabs(3),"('UT=',f5.2,' HT-INTEGRATED')") h%ut endif if (map%proj == 'ST') then ! perimeter latitude l = len_trim(map%tlabs(3)) write(map%tlabs(3)(l+1:l+16),"(' PERIM-LAT=',f5.1)") + map%perimlat endif ! ! if sattrack run, change "UT=" portion to up, down, or avg, label if(isattrack>0) then if(idetrend==0) then if(isatplt==1) map%tlabs(3)(1:8) = ' ASCND ' if(isatplt==2) map%tlabs(3)(1:8) = ' DSCND ' if(isatplt==3) map%tlabs(3)(1:8) = ' AVLEG ' else if(isatplt==1) map%tlabs(3)(1:8) = 'ASCND DT' if(isatplt==2) map%tlabs(3)(1:8) = 'DSCND DT' if(isatplt==3) map%tlabs(3)(1:8) = 'AVLEG DT' endif endif ! ! Top (1st) bottom label is min,max,interval: map%blabs(1) = ' ' if (trim(map%vectype) /= 'UNVN'.and.trim(map%vectype) /= 'UIVI') + then write(map%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(map%blabs(1),"('VECTOR MAGNITUDE MIN,MAX=',2(1pe12.4))") + vmin,vmax endif ! ! Middle (2nd) bottom label is history info: ! Bottom (3rd) bottom label is used for further info re emission fields, ! or for difference fields info. ! Only about 60 chars will show on the plot. If h%mssvol is longer ! than len1, it will be truncated to allow max of 56-char label. ! ! if (.not.diffs) then ! len0 = len_trim(h%version)+24 ! len1 = 60-len0 ! write(map%blabs(2),"(a,' ',a,' (DAY,HR,MIN=', ! | i3,',',i2,',',i2,')')") trim(h%version), ! | trim(h%mssvol(1:len1)),h%mtime ! map%blabs(3) = ' ' ! if (trim(map%ftype)=='EMISSION') then ! call mkemislab(map%sname,ie5577,ie6300,map%blabs(3)) ! if (len_trim(map%blabs(3)) > 48) blabchsz(3) = .015 ! endif ! else ! if (7+len_trim(h%mssvol) +12<=len(map%blabs(2)).and. ! + 7+len_trim(hcntr%mssvol)+12<=len(map%blabs(3)))then ! write(map%blabs(2),"('DIFFS: ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(h%mssvol),h%mtime ! write(map%blabs(3),"('MINUS ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(hcntr%mssvol),hcntr%mtime ! else ! write(map%blabs(2),"('DIFFS: (',i3,',',i2,',',i2, ! + ') MINUS (',i3,',',i2,',',i2,')')") ! + h%mtime,hcntr%mtime ! map%blabs(3) = ' ' ! endif ! endif write(map%blabs(2),"(a,' (DAY,HR,MIN=',i3,',',i2,',',i2,')')") | trim(h%version),h%mtime if (.not.diffs) then write(map%blabs(3),"(a)") trim(h%histfile) else write(map%blabs(3),"(a,' MINUS ',a)") | trim(h%histfile),trim(hcntr%histfile) endif ! ! if sattrack run, change blabs(2) to times if(isattrack>0) then write(map%blabs(2),"(a,' ',a, + ' (',i3.3,' ',i2.2,':',i2.2, + ' TO ',i3.3,' ',i2.2,':',i2.2,') ')") ! + trim(h%version),trim(h%mssvol), + trim(h%version),trim(h%histfile), + (mtimes(i,1),i=1,3),(mtimes(i,ntimes),i=1,3) endif ! ! Return message to be printed to stdout: if (map%proj == 'CE') type = 'global ' if (map%proj == 'ST') then if (map%perimlat <= 0.) type = "Spolar" if (map%perimlat > 0.) type = "Npolar" endif if (map%proj == 'MO') type = 'mollw ' if (map%proj == 'SV') type = 'satview ' if (map%levtype=='zp ') then write(msgout,"(a,' zp=',f8.3,' ',a,' min,max=', + 2e11.4)") map%sname,map%zpht,type,zmin,zmax elseif (map%levtype=='ht ') then if (map%zpht /= hmf2flag) then write(msgout,"(a,' ht=',f8.2,' ',a,' min,max=', + 2e11.4)") map%sname,map%zpht,type,zmin,zmax else write(msgout,"(a,' ht=hmf2 ',a,' min,max=', + 2e11.4)") map%sname,type,zmin,zmax endif elseif (map%levtype == 'indep') then write(msgout,"(a,' ht-indep',3x,' ',a,' min,max=', + 2e11.4)") map%sname,type,zmin,zmax else write(msgout,"(a,' ht-integ',3x,' ',a,' min,max=', + 2e11.4)") map%sname,type,zmin,zmax endif return end subroutine mkmaplabs !------------------------------------------------------------------- subroutine wrout_glb(map,h,hcntr) use plt,only: zmin,zmax ! ! Write ascii and/or xdr output data files for current glb slice frame: ! ! Args: ! (map must be intent(inout) because this routine may need ! to call mkmaplabs) ! type(mapproj),intent(inout) :: map type(history),intent(in) :: h,hcntr ! ! Locals: character(len=80) :: msgout ! if (iwrdat==0.and.iwrxdr==0) return ! ! Make labels if not done for plots: ! map%tlabs(2) label is full field name + units ! map%tlabs(3) label is ut and grid info ! map%blabs(2) label is hist vol info ! if (iplot == 0) call mkmaplabs(map,h,hcntr,msgout) ! ! Write to ascii data file: ! subroutine wrdat(iwr,lu,f,nx,ny,xx,yy,xlab,ylab,fieldlab, ! + infolab,histlab,iframe,proclab,senddat) ! if (iwrdat > 0) then call wrdat(iwrdat,ludat,map%data,map%nx,map%ny, + map%xx,map%yy,map%xlab,map%ylab,map%tlabs(2), + map%tlabs(3),map%blabs(2),iframe_dat,'tgcmproc',senddat) if (iplot==0) then write(6,"('Data frame ',i4,': ',a)") iframe_dat,trim(msgout) endif iframe_dat = iframe_dat+1 endif ! ! Write to xdr data file: if (iwrxdr > 0) then call wrxdr(flnm_xdr,map%data,map%nx,map%ny,map%xx,map%yy, + map%xlab,map%ylab,map%tlabs(2),map%tlabs(3),map%blabs(2), + map%blabs(1),h%mtime,0) if (iplot==0) then write(6,"('Xdr frame ',i4,': ',a)") iframe_xdr,trim(msgout) endif iframe_xdr = iframe_xdr+1 endif return end subroutine wrout_glb !------------------------------------------------------------------- end module mk_maps