! module mk_xylocs use hist,only: history use proc,only: nlon implicit none ! ! Make xy vertical profiles (line plots) at selected locations. ! Ht or zp on y-axis, field on x-axis. More than one curve (field) ! may be drawn on each plot (at same scale on x-axis), i.e., ncurve > 1. ! type xyvprof real :: glat,glon,slt ! geog lat,lon,slt of profile integer :: ny ! number of points real,pointer :: yy(:)=>NULL() ! y coords real,pointer :: data(:,:)=>NULL() ! field data (#curves in 2nd dim) integer :: ncurve ! number of curves (fields) on single plot integer :: ihtyax ! ht on yaxis if ihtyax==1, zp otherwise character(len=3) :: loctype ! 'lon', 'slt', 'zm ', or 'gm ' character(len=56) :: locname ! optional location name character(len=40) :: ylab ! y axis label character(len=56) :: fname ! full field name character(len=8) :: sname ! short field name character(len=8),pointer :: cnames(:)=>NULL() ! field names for each curve character(len=80) :: funits ! field units integer :: log10 ! if > 0, then log10 of field is plotted character(len=160) :: tlabs(3),blabs(3) ! top and bottom labels character(len=16) :: type ! scaler,vector,density,etc character(len=16) :: zptype ! midpoints or interfaces character(len=8) :: difftype ! raw or percent logical :: known ! true if field is known to the processor end type xyvprof real :: flon(nlon) contains !------------------------------------------------------------------- subroutine mkxyloc(f,fcntr,nf,h,hcntr) ! ! Make longitude slice plots and other output: ! use proc use fields,only: field,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use input use plt use set_vert use ampha use mk_auxflds implicit none ! ! Args: integer,intent(in) :: nf type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: type(xyvprof) :: xyprof integer :: nzprange,nhtscale,k0zp,k1zp,ixf,izpht0,izpht1, + iloc,iht,ixlon,ixlat,loginterp,ier,k,iadvfr,nohv,nohb real,pointer :: | htscale(:)=>NULL(), | zprange(:)=>NULL(),dum=>NULL() ! real :: vpxyloc(4)=(/.11,.85,.20,.90/) real :: zprof(npress),zprof_cntr(npress),fprof1(npress), + slt,vp(4) real,allocatable :: fdat1(:,:),fdat2(:,:) real :: fmin,fmax real(kind=8) :: glbminteg,glbminteg_cntr ! global mean integration real,parameter :: re=6371.e+5 ! earth radius (cm) real :: earth_area ! surface area of earth character(len=80) msgout ! ! Amphase variables: real,allocatable,save :: xyprofdat(:,:,:,:,:) ! z,nc,loc,fld,t integer,save :: itime = 0, mtime0(3) integer :: i,j,ll,m real,allocatable :: phase(:,:),amp(:,:),x(:) real :: w(200),zamp(npress) real :: vpxyloc(4)=(/.11,.85,.20,.90/) real :: tshift, cost, sint, at, bt character(len=8) :: amp_lab(4) = (/'DIURNAL ','SEMI-DI ', + 'TER-DIR ','QUAR-DI '/) ! ! Added for Bougher ascii file output: ! 11/6/13: Intel seg-faults if prof_all and aux_all are statically ! declared (mxloc,100,nf) here, so make them allocatable: ! real,allocatable :: prof_all(:,:,:) ! (mxloc,100,nf) (loc,zpht,fld) real,allocatable :: aux_all (:,:,:) ! (mxloc,100,nf) (loc,zpht,fld) real :: xylocations(2,mxloc) ! for wrdat_bf integer :: nflds=0 ! ! Externals: integer,external :: ixfind,ixfindc,ixohband logical,external :: isslt real,external :: fslt,fmean,fglbm,quadrat ! ! Begin exec: ! write(6,"(' ')") if (diffs) write(6,"('Difference fields: ')",advance='NO') write(6,"('XY location plots at ut = ',f4.1,' mtime = ',3i3)") + h%ut,h%mtime if (multiadvfr > 0) nppf = 0 ! earth_area = 4.*pi*re**2 ! surface area of earth for vertical integration ! nohv = 0 nohb = 0 xylocations(:,:) = spval ! for wrdat_bf if(iwrdat_bf > 0) then if (.not.allocated(prof_all)) then allocate(prof_all(mxloc,100,nf)) write(6,"('mkxyloc allocated prof_all: mxloc=',i5,' nf=',i4)") | mxloc,nf endif if (.not.allocated(aux_all)) then allocate(aux_all (mxloc,100,nf)) write(6,"('mkxyloc allocated aux_all: mxloc=',i5,' nf=',i4)") | mxloc,nf endif endif ! ! Field loop. Field must be requested, data pointer must be associated, ! and field must not be height-independent: ! (oh fields are plotted later with multi-curves by mkxyloc_oh) ! fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or. + .not.associated(f(ixf)%data).or. + f(ixf)%nlev==1) cycle fields_loop if (trim(f(ixf)%type)=='OH-VIB') then nohv = nohv+1 cycle fields_loop endif if (trim(f(ixf)%type)=='OH-BAND') then nohb = nohb+1 call setvert(xyloc_zprange,gcmlev,npress,dlev,nzprange,k0zp, | k1zp,zprange,xyloc_htscale,nhtscale,htscale,spval,1) cycle fields_loop endif ! ! Set up vertical scale(s) zp and/or ht: if (f(ixf)%zptype(1:3)=='MID') then ! midpoints call setvert(xyloc_zprange,gcmlev,npress,dlev,nzprange,k0zp, | k1zp,zprange,xyloc_htscale,nhtscale,htscale,spval,1) else call setvert(xyloc_zprange,gcmilev,npress,dlev,nzprange,k0zp, | k1zp,zprange,xyloc_htscale,nhtscale,htscale,spval,1) endif ! ! Set loop indices for zp and/or ht on y-axis: izpht0 = 0 izpht1 = 1 if (nhtscale == 0) izpht1 = 0 if (nzprange == 0) izpht0 = 1 ! nflds = nflds + 1 xyprof%ncurve = 1 if (associated(xyprof%cnames)) deallocate(xyprof%cnames) allocate(xyprof%cnames(xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprof%cnames") xyprof%cnames = ' ' xyprof%fname = f(ixf)%fname56 xyprof%sname = f(ixf)%fname8 xyprof%funits = f(ixf)%units if (f(ixf)%scalefac /= 1.) | write(xyprof%funits,"(a,' * ',1pe9.2)") trim(f(ixf)%units), | f(ixf)%scalefac xyprof%type = f(ixf)%type xyprof%zptype = f(ixf)%zptype xyprof%difftype = f(ixf)%difftype xyprof%known = f(ixf)%known ! ! Location loop: loc_loop: do iloc=1,mxloc if (xylocs(1,iloc)==spval.or.xylocs(2,iloc)==spval) + cycle loc_loop call setloc(ixlat,ixlon) ! ! Define locations for wrdat_bf (Bougher file format) xylocations(1,iloc) = xyprof%glat xylocations(2,iloc) = xyprof%glon if (isslt(xylocs(2,iloc),slt)) | xylocations(2,iloc) = xyprof%slt xyprof%locname = "" if (len_trim(xyloc_locname(iloc))>0) + xyprof%locname = ' ('//trim(xyloc_locname(iloc))//')' ! ! Save heights: ! 1/22/08 btf: Use z_ifaces or z_midpts, consistent w/ current fields's ! zptype (see sub mkgeopot in getflds.F) ! ! Field is plotted at interfaces -> use Z at interfaces (z_ifaces): if (trim(f(ixf)%zptype) == 'INTERFACES') then if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ')then zprof = z_ifaces(ixlon,ixlat,:) if (diffs) zprof_cntr = zcntr_ifaces(ixlon,ixlat,:) elseif (xyprof%loctype == 'zm ') then ! zonal means do k=1,npress zprof(k) = fmean(z_ifaces(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zprof_cntr(k)=fmean(z_ifaces(:,ixlat,k),nlon,spval,0) enddo endif else ! global means do k=1,npress zprof(k) = + fglbm(z_ifaces(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (diffs) then do k=1,npress zprof_cntr(k) = + fglbm(zcntr_ifaces(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif endif ! ! Field is plotted at midpoints -> use Z at midpoints (z_midpts): else ! MIDPOINTS if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ')then zprof = z_midpts(ixlon,ixlat,:) if (diffs) zprof_cntr = zcntr_midpts(ixlon,ixlat,:) elseif (xyprof%loctype == 'zm ') then ! zonal means do k=1,npress zprof(k) = fmean(z_midpts(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zprof_cntr(k)=fmean(z_midpts(:,ixlat,k),nlon,spval,0) enddo endif else ! global means do k=1,npress zprof(k) = + fglbm(z_midpts(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (diffs) then do k=1,npress zprof_cntr(k) = + fglbm(zcntr_midpts(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif endif endif ! ! Zp and/or ht on y-axis (0,0 for zp only, 1,1 for ht only, ! or 0,1 for zp and ht): yax_loop: do iht=izpht0,izpht1 xyprof%ihtyax = iht ! ! - - - - - - - - - - - - zp on y-axis: - - - - - - - - - - - - - if (iht == 0) then ! zp on y-axis ! ! Adjust right viewport according to number of extra right-hand y-axes: if (ixyloc_yaxright==0) then vpxyloc(2) = .90 elseif (ixyloc_yaxright==1) then vpxyloc(2) = .85 else vpxyloc(2) = .76 endif ! ! Define y coords in pressure: xyprof%ny = nzprange if (associated(xyprof%yy)) deallocate(xyprof%yy) allocate(xyprof%yy(xyprof%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprof%yy") xyprof%yy = zprange ! array op write(xyprof%ylab,"('LN(P0/P) (',a,')')") | trim(xyprof%zptype) ! ! Allocate and define xyprof data for zp on y-axis: if (associated(xyprof%data)) deallocate(xyprof%data) allocate(xyprof%data(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprof%data") ! ! Allocate temp space for diffs: if (diffs) then allocate(fdat1(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating fdat1") allocate(fdat2(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating fdat2") endif ! ! Lat,lon or lat,slt: ! (If diffs, save perturbed in fdat1 and control in fdat2) ! if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ') + then xyprof%data(:,1) = f(ixf)%data(ixlon,ixlat,k0zp:k1zp) if (diffs) then fdat1 = xyprof%data ! save perturbed fdat2(:,1) = fcntr(ixf)%data(ixlon,ixlat,k0zp:k1zp) endif ! ! Zonal means: elseif (xyprof%loctype == 'zm ') then do k=k0zp,k1zp flon(:) = f(ixf)%data(:,ixlat,k) xyprof%data(k-k0zp+1,1) = ! tmp only one curve + fmean(flon,nlon,spval,0) enddo if (diffs) then fdat1 = xyprof%data ! save perturbed do k=k0zp,k1zp flon(:) = fcntr(ixf)%data(:,ixlat,k) fdat2(k-k0zp+1,1) = ! tmp only one curve + fmean(flon,nlon,spval,0) enddo endif ! ! Global means: else do k=k0zp,k1zp ! zp levels are gcmlev(k0zp:k1zp) xyprof%data(k-k0zp+1,1) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo ! ! Integration of global means in zp: ! real function quadrat(f,hts,npts) is in util.F if (trim(f(ixf)%type)=='DENSITY') then glbminteg = quadrat(xyprof%data(:,1),zprof(k0zp:k1zp), | nzprange) ! write(6,"(/,'Field ',a,' glbmean vertical ', ! | 'integration (zprange) pert = ',e12.4,' total=', ! | e12.4)") trim(xyprof%sname),glbminteg, ! | glbminteg*earth_area endif if (diffs) then fdat1 = xyprof%data ! save perturbed do k=k0zp,k1zp fdat2(k-k0zp+1,1) = + fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo ! ! Diffs of vertical integration: if (trim(f(ixf)%type)=='DENSITY') then glbminteg_cntr = quadrat(fdat2(:,1),zprof(k0zp:k1zp) | ,nzprange) ! write(6,"('Field ',a,' glbmean vertical ', ! | 'integration (zprange) cntr = ',e12.4,' total=', ! | e12.4)") trim(xyprof%sname),glbminteg_cntr, ! | glbminteg_cntr*earth_area ! write(6,"('Field ',a,' glbmean vertical ', ! | 'integration (zprange) diff = ',e12.4,' total=', ! | e12.4)") trim(xyprof%sname), ! | glbminteg-glbminteg_cntr, ! | (glbminteg-glbminteg_cntr)*earth_area endif endif endif ! ! Difference fields in pressure: if (diffs) then call mkdiffs(fdat1,fdat2,xyprof%data, + xyprof%ny*xyprof%ncurve,f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif ! ! - - - - - - - - - - - - ht on y-axis: - - - - - - - - - - - - - else if (ixyloc_yaxright==0) then vpxyloc(2) = .90 else vpxyloc(2) = .85 endif ! ! Define y coords in height: xyprof%ny = nhtscale if (associated(xyprof%yy)) deallocate(xyprof%yy) allocate(xyprof%yy(xyprof%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprof%yy") xyprof%yy = htscale ! array op write(xyprof%ylab,"('HEIGHT (KM)')") ! ! Allocate and define xyprof data for zp on y-axis: if (associated(xyprof%data)) deallocate(xyprof%data) allocate(xyprof%data(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprof%data") loginterp = 0 if (trim(f(ixf)%type)=='DENSITY') loginterp = 1 ! ! Allocate temp space for diffs: if (diffs) then allocate(fdat1(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating fdat1") allocate(fdat2(xyprof%ny,xyprof%ncurve),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating fdat2") endif ! ! If doing zonal or global means, calculate means of the field ! (store in fprof1), and means of Z (store in zprof), then do ht interp ! of the means. ! If doing diffs, using ave pert+cntr heights in the interpolation ! (intloc) helps suppress jagged lines. (taking diffs before ! interpolation will also suppress jagged lines, but should ! always interpolate first) ! ! Lat,lon or lat,slt: if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ') + then if (.not.diffs) then call intloc(f(ixf)%data(ixlon,ixlat,:),zprof, + npress,htscale,nhtscale,loginterp,xyprof%data,1,1, + nhtscale,ier,spval,0) else ! diffs call intloc(f(ixf)%data(ixlon,ixlat,:),zprof, + npress,htscale,nhtscale,loginterp, + xyprof%data,1,1,nhtscale,ier,spval,0) fdat1 = xyprof%data call intloc(fcntr(ixf)%data(ixlon,ixlat,:), + zprof_cntr,npress,htscale,nhtscale, + loginterp,fdat2,1,1,nhtscale,ier,spval,0) endif ! ! Zonal means: elseif (xyprof%loctype == 'zm ') then do k=1,npress flon(:) = f(ixf)%data(:,ixlat,k) fprof1(k) = fmean(flon,nlon,spval,0) enddo if (.not.diffs) then call intloc(fprof1,zprof,npress,htscale, nhtscale, + loginterp,xyprof%data,1,1,nhtscale,ier,spval,0) else call intloc(fprof1,zprof,npress, + htscale, nhtscale,loginterp,xyprof%data,1,1, + nhtscale,ier,spval,0) fdat1 = xyprof%data do k=1,npress flon(:) = fcntr(ixf)%data(:,ixlat,k) fprof1(k) = fmean(flon,nlon,spval,0) enddo call intloc(fprof1,zprof_cntr,npress, + htscale,nhtscale,loginterp,fdat2,1,1,nhtscale, + ier,spval,0) endif ! ! Global means: else do k=1,npress fprof1(k) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (.not.diffs) then call intloc(fprof1,zprof,npress,htscale,nhtscale, + loginterp,xyprof%data,1,1,nhtscale,ier,spval,0) ! ! Calculate vertical integration of global means (height): if (trim(f(ixf)%type)=='DENSITY') then glbminteg = quadrat(xyprof%data(:,1),htscale, | nhtscale) write(6,"(/,'Field ',a,' glbmean vertical ', | 'integration (htscale)= ',e12.4,' total=',e12.4)") | trim(xyprof%sname),glbminteg,glbminteg*earth_area endif else ! diffs call intloc(fprof1,zprof,npress, + htscale,nhtscale,loginterp,xyprof%data,1,1, + nhtscale,ier,spval,0) fdat1 = xyprof%data do k=1,npress fprof1(k) = fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat, + gcmlat,dlat,dlon,spval) enddo call intloc(fprof1,zprof_cntr,npress, + htscale,nhtscale,loginterp,fdat2,1,1,nhtscale,ier, + spval,0) ! ! Diffs of vertical integration: if (trim(f(ixf)%type)=='DENSITY') then glbminteg = quadrat(fdat1(:,1),htscale,nhtscale) ! pert write(6,"(/,'Field ',a,' glbmean vertical ', | 'integration (htscale) pert = ',e12.4,' total=', | e12.4)") trim(xyprof%sname),glbminteg, | glbminteg*earth_area glbminteg_cntr = quadrat(fdat2,htscale,nhtscale) ! cntr write(6,"('Field ',a,' glbmean vertical ', | 'integration (htscale) cntr = ',e12.4,' total=', | e12.4)") trim(xyprof%sname),glbminteg_cntr, | glbminteg_cntr*earth_area write(6,"('Field ',a,' glbmean vertical ', | 'integration (htscale) diff = ',e12.4,' total=', | e12.4)") trim(xyprof%sname), | glbminteg-glbminteg_cntr, | (glbminteg-glbminteg_cntr)*earth_area endif endif endif ! lon, slt, zm, or gm ! ! Difference fields at heights: if (diffs) then call mkdiffs(fdat1,fdat2,xyprof%data, + xyprof%ny*xyprof%ncurve,f(ixf)%difftype) ! ! 5-pt smoother in height: if (ismooth > 0) ! input module (default == 1) | call smoother(xyprof%data,xyprof%ny,5,spval,1) ! deallocate(fdat1) ; deallocate(fdat2) endif endif ! zp or ht ! - - - - - - - - - - - - end zp or ht on y-axis: - - - - - - - - - ! ! Scale data if fscale is requested (scfac is in plt.f, s.a. sub pltmxy): scfac = f(ixf)%scalefac ! ! xyprof%log10 = 0 -> all fields linear ! xyprof%log10 = 1 -> log10 of densities ! xyprof%log10 = 2 -> log10 of all fields xyprof%log10 = 0 if (ixyloc_log10==1.and.trim(f(ixf)%type)=='DENSITY') + xyprof%log10 = ixyloc_log10 if (ixyloc_log10 > 1) xyprof%log10 = 1 ! if (all(xyprof%data==spval)) then write(6,"('>>> mkxyloc: skipping this plot because ', + 'all values are missing (spval).')") cycle yax_loop endif call setmultivp(vpxyloc,iadvfr,nppf+1,multiplt,ipltrowcol, + vp) ! ! Do not plot or write data if doing amphase plots if(amphase/=0) cycle yax_loop ! ! Set up viewport (for multiplt): if (iplot > 0) then ! ! Make labels and draw the xy line plot: call mkxyloclabs(xyprof,h,hcntr,msgout) call pltxyloc(xyprof,zprof,10*iht+ixyloc_yaxright,vp,h%ut) ! ! 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, + 'xy vert profs at locs',iframe) ! ! If not making plots, wrout will need zmin,zmax. else call fminmax(xyprof%data,size(xyprof%data),zmin,zmax) endif ! ! Write output data files: call wrout_xyloc(xyprof,h,hcntr) ! ! Save profile if writing a Bougher-formatted file if(iwrdat_bf > 0) prof_all(iloc,1:xyprof%ny,ixf) + = xyprof%data(1:xyprof%ny,1) enddo yax_loop ! ! Save data for amphase plot, if requested if (amphase > 0) then if (itime==0) then if (allocated(xyprofdat)) deallocate(xyprofdat) allocate(xyprofdat(xyprof%ny,xyprof%ncurve,mxloc, + nf,25),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc allocating xyprofdat") mtime0 = h%mtime endif itime = (h%mtime(1)-mtime0(1))*24 + + (h%mtime(2)-mtime0(2)) + 1 do i=1,xyprof%ny do j=1,xyprof%ncurve xyprofdat(i,j,iloc,ixf,itime) = xyprof%data(i,j) enddo ! j enddo ! i ! if (itime==25) then if (allocated(amp)) deallocate(amp) if (allocated(phase)) deallocate(phase) if (allocated(x)) deallocate(x) allocate(amp(xyprof%ny,amphase),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating amp") allocate(phase(xyprof%ny,amphase),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating phase") allocate(x(25),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating x") call rffti(25,w) do k=1,xyprof%ny x = xyprofdat(k,1,iloc,ixf,:) call rfftf(25,x,w) do ll=1,amphase m=ll+1 amp(k,ll)=sqrt(x(2*m)**2+x(2*m-1)**2)/ + float(25) tshift = float(ll)*pi*h%ut/12. cost = cos(tshift) sint = sin(tshift) at = x(2*m)*cost+x(2*m-1)*sint bt = x(2*m-1)*cost-x(2*m)*sint if(bt.eq.0.)go to 3 phase(k,ll) = 12./pi*atan2(at,bt) go to 1 3 continue if(at.eq.0.)go to 4 phase(k,ll) = 12./pi*2.*atan(at/abs(at)) go to 1 4 continue phase(k,ll) = 0. 1 continue ! ! Switch to local time, adjust for semi- or ter-diurnal phase(k,ll) = phase(k,ll)+xylocs(2,iloc)/15. if(phase(k,ll) .lt. 0) + phase(k,ll)=phase(k,ll)+24. if(ll == 2 .and. phase(k,ll) .gt. 12.) + phase(k,ll)=phase(k,ll)-12. if(ll == 3 .and. phase(k,ll) .gt. 16.) + phase(k,ll)=phase(k,ll)-8. if(ll == 3 .and. phase(k,ll) .gt. 8.) + phase(k,ll)=phase(k,ll)-8. enddo ! ll enddo ! k if (trim(f(ixf)%zptype) == 'INTERFACES') then zamp = z_ifaces(ixlon,ixlat,:) else zamp = z_midpts(ixlon,ixlat,:) endif ! do ll=1,amphase xyprof%data(:,1) = amp(:,ll) call mkxyloclabs(xyprof,h,hcntr,msgout) xyprof%tlabs(3)(1:8) = amp_lab(ll) call pltxyloc_amphase(xyprof,zamp, + 10*iht+ixyloc_yaxright,vp,h%ut,1) ! nppf = nppf+1 ! call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, ! + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, ! + 'xy amp profs at locs',iframe) iwr_label(4) = 0 iwr_label(5) = 0 iwr_label(6) = 0 call wrout_xyloc(xyprof,h,hcntr) ! xyprof%data(:,1) = phase(:,ll) call mkxyloclabs(xyprof,h,hcntr,msgout) xyprof%tlabs(3)(1:8) = amp_lab(ll) call pltxyloc_amphase(xyprof,zamp, + 10*iht+ixyloc_yaxright,vp,h%ut,2) nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'xy amp profs at locs',iframe) call wrout_xyloc(xyprof,h,hcntr) enddo ! ll ! endif endif ! enddo loc_loop enddo fields_loop ! ! Plot oh fields (multiple curves): if (nhtscale > 0) then if (nohv > 0) call mkxyloc_oh('OH-VIB ',nohv) if (nohb > 0) call mkxyloc_oh('OH-BAND ',nohb) endif ! ! 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,' ','xy vert profs at locs', + iframe) ! ! Output Bougher-formatted file, if desired if(iwrdat_bf > 0) then ! subroutine mkauxdev(prof_all,ny,nf,nflds,f,aux_all) ! real,intent(in) :: prof_all(mxloc,100,nf) ! real,intent(out) :: aux_all(mxloc,100,10) call mkauxdev(prof_all,xyprof%ny,nf,nflds,f,aux_all) ! ! subroutine wrdat_bf(lu,h,xylocations,prof_all,aux_all,nx,maxy, ! | ny,maxz,nz,yy,ut,proclab) ! real,intent(in) :: prof_all(nx,maxy,maxz) ! (loc,zpht,fld) ! allocate(prof_all(mxloc,100,nf)) ! above call wrdat_bf(ludat_bf,h,xylocations,prof_all,aux_all,mxloc, | 100,xyprof%ny,nf,nflds,xyprof%yy,h%ut,'tgcmproc_bf') endif ! ! Release local space: if (associated(xyprof%data)) deallocate(xyprof%data) if (associated(xyprof%yy)) deallocate(xyprof%yy) if (associated(htscale)) deallocate(htscale) if (associated(zprange)) deallocate(zprange) return contains !------------------------------------------------------------------- subroutine mkxyloc_oh(type,ntype) ! internal to mkxyloc ! ! Make multi-curve vertical profiles of oh-v or oh-b fields at selected ! locations (one frame per location, one curve per vib level or ! emission band). ! On input, type='OH-VIB' or 'OH-BAND', and ntype = number of fields ! of type that have been requested. ! ! Args: character(len=16),intent(in) :: type ! OH-VIB or OH-BAND integer,intent(in) :: ntype ! number of fields of type ! Locals: real,pointer :: hscale(:)=>NULL() ! height scale to be used integer :: nhscale ! dimension for hscale real,pointer :: fhscale(:)=>NULL()! height scale of field integer :: nfhscale ! dimension for fhscale integer :: ifoh,ix,len,i,hi,lo real,allocatable :: fprof(:),fprof1(:),fprof2(:) ! ! Get height scale for y-axis (assume all fields of type have same ! height range): floop: do ixf=1,nf if (f(ixf)%type==type) then call sethscale(f(ixf),fhscale,nfhscale,htscale,nhtscale, + hscale,nhscale) exit floop endif enddo floop xyprof%ny = nhscale xyprof%ihtyax = 1 if (associated(xyprof%yy)) deallocate(xyprof%yy) allocate(xyprof%yy(xyprof%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc_oh allocating xyprof%yy") xyprof%yy = hscale ! array op write(xyprof%ylab,"('HEIGHT')") ! ! fprof is temporary storage for a profile at the range of the field: ! (e.g., in the case of oh, 90-120 by 1 km) if (allocated(fprof)) deallocate(fprof) allocate(fprof(nfhscale),stat=ier) fprof = spval ! init if (diffs) then if (allocated(fprof1)) deallocate(fprof1) allocate(fprof1(nfhscale),stat=ier) fprof1 = spval ! init if (allocated(fprof2)) deallocate(fprof2) allocate(fprof2(nfhscale),stat=ier) fprof2 = spval ! init endif ! ! Set field names (this is for top title and x-lab -- individual ! curve names are set below (xyprof%cnames)): if (trim(type)=='OH-VIB') then xyprof%fname = 'OH VIBRATIONAL LEVELS' xyprof%sname = 'OH-VIBS ' endif if (trim(type)=='OH-BAND') then xyprof%fname = 'OH EMISSION BANDS' xyprof%sname = 'OH-BANDS' endif ! ! Location loop: loc_loop: do iloc=1,mxloc if (xylocs(1,iloc)==spval.or.xylocs(2,iloc)==spval) + cycle loc_loop ! write(6,"('mkxyloc before setloc: iloc=',i2)") iloc call setloc(ixlat,ixlon) xyprof%locname = "" if (len_trim(xyloc_locname(iloc))>0) + xyprof%locname = ' ('//trim(xyloc_locname(iloc))//')' if (associated(xyprof%data)) deallocate(xyprof%data) allocate(xyprof%data(nhscale,ntype),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc_oh allocating xyprof%data") ! write(6,"('mkxyloc allocated data')") xyprof%data = spval ! init if (associated(xyprof%cnames)) deallocate(xyprof%cnames) ! write(6,"('before cnames allocate')") allocate(xyprof%cnames(ntype),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyloc_oh allocating xyprof%cnames") ! write(6,"('after allocate xyprof%cnames: ntype=',i2)") ! | ntype ! ! Fields loop (oh-v and/or oh-b only): ifoh = 0 fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or. + .not.associated(f(ixf)%data).or. + f(ixf)%type/=type) cycle fields_loop ifoh = ifoh+1 xyprof%ncurve = ntype xyprof%funits = f(ixf)%units ! if (scfac /= 1.) ! | write(xyprof%funits,"(a,'X',1pe12.4)") trim(f(ixf)%units,scfac xyprof%type = f(ixf)%type ! ! Field names for oh-vib are: OHV-1, OHV-2,...,OHV-10 ! (label curves V-1, V-2, etc) if (trim(type)=='OH-VIB') then len = len_trim(f(ixf)%fname8) if (f(ixf)%fname8(len-1:len-1)=='-') then ! single digit level read(f(ixf)%fname8(len:len),"(i1)") i write(xyprof%cnames(ifoh),"('V-',i1)") i else read(f(ixf)%fname8(len-1:len),"(i2)") i write(xyprof%cnames(ifoh),"('V-',i2)") i endif endif ! ! Field names for oh-band are, eg: OHB-62, OHB-98, etc. ! (label curves 6-2, 9-8, etc) if (trim(type)=='OH-BAND') then i = ixohband(f(ixf)%fname8,hi,lo) write(xyprof%cnames(ifoh),"(i1,'-',i1)") hi,lo endif ! ! 1/22/08 btf: Use z_ifaces or z_midpts, consistent w/ current fields's ! zptype (see sub mkgeopot in getflds.F) ! ! Field is at interfaces -> save Z at interfaces: if (trim(f(ixf)%zptype) == 'INTERFACES') then if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ')then zprof(:) = z_ifaces(ixlon,ixlat,:) if (diffs) zprof_cntr(:) = zcntr_ifaces(ixlon,ixlat,:) elseif (xyprof%loctype == 'zm ') then do k=1,npress zprof(k) = fmean(z_ifaces(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zprof_cntr(k) = | fmean(zcntr_ifaces(:,ixlat,k),nlon,spval,0) enddo endif else ! global means do k=1,npress zprof(k) = fglbm(z_ifaces(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo if (diffs) then do k=1,npress zprof_cntr(k) = fglbm(zcntr_ifaces(:,:,k),nlon,nlat, + gcmlat,dlat,dlon,spval) enddo endif endif ! ! Field is at midpoints -> save Z at midpoints: else ! MIDPOINTS if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ')then zprof(:) = z_midpts(ixlon,ixlat,:) if (diffs) zprof_cntr(:) = zcntr_midpts(ixlon,ixlat,:) elseif (xyprof%loctype == 'zm ') then do k=1,npress zprof(k) = fmean(z_midpts(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zprof_cntr(k) = | fmean(zcntr_midpts(:,ixlat,k),nlon,spval,0) enddo endif else ! global means do k=1,npress zprof(k) = fglbm(z_midpts(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo if (diffs) then do k=1,npress zprof_cntr(k) = fglbm(zcntr_midpts(:,:,k),nlon,nlat, + gcmlat,dlat,dlon,spval) enddo endif endif endif ! ! Lat,lon or lat,slt: ! if (xyprof%loctype /= 'zm '.and.xyprof%loctype /= 'gm ') then fprof(:) = f(ixf)%data(ixlon,ixlat,:) if (diffs) then fprof1(:) = fcntr(ixf)%data(ixlon,ixlat,:) endif ! ! Zonal means: elseif (xyprof%loctype == 'zm ') then do k=1,nfhscale flon(:) = f(ixf)%data(:,ixlat,k) fprof(k) = fmean(flon,nlon,spval,0) enddo if (diffs) then do k=1,nfhscale flon(:) = fcntr(ixf)%data(:,ixlat,k) fprof1(k) = fmean(flon,nlon,spval,0) enddo endif ! ! Global means: else do k=1,nfhscale fprof(k) = fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo if (diffs) then do k=1,nfhscale fprof1(k) = fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat, + gcmlat,dlat,dlon,spval) enddo endif endif ! ! Difference fields: if (diffs) then call mkdiffs(fprof,fprof1,fprof2,npress,f(ixf)%difftype) fprof = fprof2 endif ! ! Define data within range of user requested height scale: do k=1,nhscale if (hscale(k) >= fhscale(1).and. + hscale(k) <= fhscale(nfhscale)) then ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev) xyprof%data(k,ifoh) = fprof(ix) endif enddo if (ifoh==ntype) exit fields_loop enddo fields_loop xyprof%log10 = ixyloc_log10 if (all(xyprof%data==spval)) then write(6,"('>>> mkxyloc_oh: skipping this plot because ', + 'all values are missing (spval).')") cycle loc_loop endif ! ! Set up viewport (for multiplt): if (iplot > 0 .and. amphase==0) then call setmultivp(vpxyloc,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Make labels and draw the xy line plot: call mkxyloclabs(xyprof,h,hcntr,msgout) call pltxyloc(xyprof,zprof,10+ixyloc_yaxright,vp,h%ut) ! ! 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, + 'xy vert profs at locs',iframe) ! ! If not making plots, wrout will need zmin,zmax. else call fminmax(xyprof%data,size(xyprof%data),zmin,zmax) endif ! ! Write output data files: call wrout_xyloc(xyprof,h,hcntr) enddo loc_loop if (associated(fhscale)) deallocate(fhscale) if (associated(hscale)) deallocate(hscale) if (allocated(fprof)) deallocate(fprof) if (allocated(fprof1)) deallocate(fprof1) if (allocated(fprof2)) deallocate(fprof2) return end subroutine mkxyloc_oh !------------------------------------------------------------------- subroutine setloc(ixlat,ixlon) ! internal to mkxyloc ! ! Set location of xyprof, and return ixlat,ixlon: ! (Called by mkxyloc and mkxyloc_oh) ! ! xylocs(1,iloc) is latitude, xylocs(2,iloc) may be lon, slt, or zm. ! (lats were validated and set to nearest grid point by chlats) ! (if regular lon, then nearest grid point was already set by chlons) ! (if xylocs(:,iloc)=='zm', (i.e., both lat,lon) then do global means) ! ! Args: integer,intent(out) :: ixlat,ixlon ! xyprof%loctype = ' ' ! ! Global means: if (xylocs(1,iloc)==zmflag.and.xylocs(2,iloc)==zmflag) then xyprof%loctype = 'gm ' ! ! Zonal means: elseif (xylocs(2,iloc) == zmflag) then ! zonal means xyprof%loctype = 'zm ' xyprof%glon = zmflag ! ! Regular longitude or local time: else if (isslt(xylocs(2,iloc),slt)) then ! is local time xyprof%loctype = 'slt' xyprof%glon = fslt(slt,h%ut,dum,3) ixlon = ixfind(gcmlon,nlon,xyprof%glon,dlon) xyprof%glon = gcmlon(ixlon) xyprof%slt = slt else ! is regular longitude xyprof%loctype = 'lon' xyprof%glon = xylocs(2,iloc) ixlon = ixfind(gcmlon,nlon,xyprof%glon,dlon) xyprof%slt = fslt(dum,h%ut,xyprof%glon,1) endif endif ! ! Set latitude if not global means: if (xyprof%loctype /= 'gm ') then xyprof%glat = xylocs(1,iloc) ixlat = ixfind(gcmlat,nlat,xyprof%glat,dlat) endif end subroutine setloc end subroutine mkxyloc !------------------------------------------------------------------- subroutine mkxyloclabs(xyprof,h,hcntr,msgout) use proc,only: diffs,spval use hist,only: hdr use plt,only: scfac ! ! Construct 3 top labels (xyprof%tlabs) and 3 bottom labels ! (xyprof%blabs) for xy vertical profiles: ! ! Args: type(xyvprof),intent(inout) :: xyprof type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Locals: integer :: lenlab,iday real :: fmin,fmax,sza ! ! Externals: real,external :: getsza ! ! Set up top and bottom text labels: ! ! Top (1st) top label is optionally provided by user: xyprof%tlabs(1) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if xyprof%log10 == 0) ! if (xyprof%known) then lenlab = len_trim(xyprof%fname)+len_trim(xyprof%funits)+3 if (xyprof%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(xyprof%tlabs(2))) then if (xyprof%log10 > 0) then xyprof%tlabs(2) = 'LOG10 '//trim(xyprof%fname)//' ('// + trim(xyprof%funits)//')' else xyprof%tlabs(2) = trim(xyprof%fname)//' ('// + trim(xyprof%funits)//')' endif else xyprof%tlabs(2) = trim(xyprof%fname) endif else ! unknown to proc lenlab = len_trim(xyprof%sname)+2+len_trim(xyprof%fname)+ | len_trim(xyprof%funits)+3 if (xyprof%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(xyprof%tlabs(2))) then if (xyprof%log10 > 0) then if (trim(xyprof%sname)/=trim(xyprof%fname)) then xyprof%tlabs(2) = 'LOG10 '//trim(xyprof%sname)//': '// | trim(xyprof%fname)//' ('//trim(xyprof%funits)//')' else xyprof%tlabs(2) = 'LOG10 '//trim(xyprof%fname)// | ' ('//trim(xyprof%funits)//')' endif else ! no log10 if (trim(xyprof%sname)/=trim(xyprof%fname)) then xyprof%tlabs(2) = trim(xyprof%sname)//': '// | trim(xyprof%fname)//' ('//trim(xyprof%funits)//')' else xyprof%tlabs(2) = trim(xyprof%fname)// | ' ('//trim(xyprof%funits)//')' endif endif else xyprof%tlabs(2) = trim(xyprof%fname) endif endif ! known or unknown ! ! Add DIFFS to label: if (diffs) then if (trim(xyprof%difftype)=='RAW') then if (len_trim(xyprof%tlabs(2))+7 <= len(xyprof%tlabs(2))) then xyprof%tlabs(2) = 'DIFFS: '//trim(xyprof%tlabs(2)) else xyprof%tlabs(2) = 'DIFFS: '//trim(xyprof%fname) endif else ! percent if (len_trim(xyprof%tlabs(2))+15 <= len(xyprof%tlabs(2))) then xyprof%tlabs(2) = trim(xyprof%difftype)//' DIFFS: '// | trim(xyprof%tlabs(2)) else xyprof%tlabs(2) = '% DIFFS: '//trim(xyprof%fname) endif endif endif ! ! Bottom (3rd) top label is ut and grid info plus optional location name: ! xyloc location may be lat,lon, lat,slt, zonal means, or global means: if (xyprof%loctype == 'gm ') then write(xyprof%tlabs(3),"('UT=',f5.2,' GLOBAL MEANS')") h%ut elseif (xyprof%loctype == 'zm ') then write(xyprof%tlabs(3),"('UT=',f5.2,' LAT=',f7.2, + ' ZONAL MEANS')") h%ut,xyprof%glat else ! lat,lon or lat,slt if (xyprof%loctype=='lon') then write(xyprof%tlabs(3),"('UT=',f5.2,' LAT=',f7.2, + ' LON=',f9.2,' (DEG)')") h%ut,xyprof%glat,xyprof%glon else write(xyprof%tlabs(3),"('UT=',f5.2,' LAT=',f7.2, + ' SLT=',f5.2,' (HRS))')") + h%ut,xyprof%glat,xyprof%slt endif endif if (len_trim(xyprof%locname)>0) then if (len_trim(xyprof%tlabs(3))+len_trim(xyprof%locname)<= + len(xyprof%tlabs(3))) then xyprof%tlabs(3) = trim(xyprof%tlabs(3))//trim(xyprof%locname) else write(6,"('Note: no room to add locname ',a)") + len_trim(xyprof%locname) endif endif ! ! Top (1st) bottom label is min,max,interval: ! call fminmax(xyprof%data,xyprof%ny*xyprof%ncurve,fmin,fmax) if (scfac == 1.) then write(xyprof%blabs(1),"('MIN,MAX=',2(1pe12.4))") + fmin,fmax else fmin = merge(fmin,fmin*scfac,fmin==spval) fmax = merge(fmax,fmax*scfac,fmax==spval) write(xyprof%blabs(1),"('MIN,MAX=',2(1pe12.4),' (X',1pe9.2, | ')')") fmin,fmax,scfac endif ! ! Middle (2nd) bottom label is history info: ! Bottom (3rd) bottom label is also used if doing diffs or ! excitation rates: ! ! if (.not.diffs) then ! write(xyprof%blabs(2),"(a,' ',a,' (DAY,HR,MIN=', ! + i3,',',i2,',',i2,')')") trim(h%version),trim(h%mssvol), ! + h%mtime ! xyprof%blabs(3) = ' ' ! if (trim(xyprof%type)=='EXCITED-STATE') then ! iday = h%iyd-h%iyd/1000*1000 ! sza = getsza(iday,xyprof%slt,xyprof%glat,xyprof%glon) ! write(xyprof%blabs(3),"('F10.7=',f6.2,' Solar Zenith', ! + ' Angle = ',f6.2,' (deg)')") hdr%f107d,sza ! endif ! else ! if (7+len_trim(h%mssvol) +12<=len(xyprof%blabs(2)).and. ! + 7+len_trim(hcntr%mssvol)+12<=len(xyprof%blabs(3)))then ! write(xyprof%blabs(2),"('DIFFS: ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(h%mssvol),h%mtime ! write(xyprof%blabs(3),"('MINUS ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(hcntr%mssvol),hcntr%mtime ! else ! write(xyprof%blabs(2),"('DIFFS: (',i3,',',i2,',',i2, ! + ') MINUS (',i3,',',i2,',',i2,')')") ! + h%mtime,hcntr%mtime ! xyprof%blabs(3) = ' ' ! endif ! endif ! ! 5/05: Put version name and model time in blabs(2), history disk ! file(s) in blabs(3) ! write(xyprof%blabs(2),"(a,' (DAY,HR,MIN=',i3,',',i2,',',i2,')')") | trim(h%version),h%mtime if (.not.diffs) then write(xyprof%blabs(3),"(a)") trim(h%histfile) else write(xyprof%blabs(3),"(a,' MINUS ',a)") | trim(h%histfile),trim(hcntr%histfile) endif ! ! Return message to be printed to stdout: if (xyprof%loctype == 'gm ') then ! global means write(msgout,"(a,' global means: min,max=', + 2e11.4,' iht=',i1)") xyprof%sname,fmin,fmax,xyprof%ihtyax elseif (xyprof%loctype == 'lon') then ! longitude write(msgout,"(a,' lat,lon=',f5.1,',',f6.1,' min,max=', + 2e11.4,' iht=',i1)") xyprof%sname,xyprof%glat,xyprof%glon, + fmin,fmax,xyprof%ihtyax elseif (xyprof%loctype == 'slt') then ! local time write(msgout,"(a,' lat,slt =',f7.2,',',f9.2,' min,max=', + 2e11.4,' iht=',i1)") xyprof%sname,xyprof%glat,xyprof%slt, + fmin,fmax,xyprof%ihtyax else ! zonal means write(msgout,"(a,' lat=',f7.2,' (ZM) min,max=', + 2e11.4,' iht=',i1)") xyprof%sname,xyprof%glat,fmin,fmax, + xyprof%ihtyax endif return end subroutine mkxyloclabs !------------------------------------------------------------------- subroutine pltxyloc(xyprof,hts,iyax,vp,ut) use plt,only: igks_cgm,igks_ps,igks_x11, + iwk_cgm, iwk_ps, iwk_x11, yaxzpht,scfac use proc,only: spval,gcmlev,gcmilev,npress,dlev,p0 use input,only: ilab_hq ! ! Make xy line plot with field on x-axis, vertical on y-axis ! at selected location and ut (called from mkxyloc): ! ! Args: type(xyvprof),intent(in) :: xyprof real,intent(in) :: hts(:),vp(4),ut integer,intent(in) :: iyax ! ! Locals: integer :: iplt_cgm,iplt_ps,iplt_x11,ier,multix,i real :: dum(1,1) real :: toffset=.04, boffset=.15 real :: tlabchsz(3)=(/.02,.02,.02/), + blabchsz(3)=(/.02,.02,.02/) ! ! Begin exec: multix = 0 if (associated(xyprof%cnames)) then do i=1,size(xyprof%cnames) if (len_trim(xyprof%cnames(i)) > 0) multix = multix+1 enddo endif ! ! Draw to cgm workstation: iplt_cgm = 0 call pltmxy(xyprof%data,xyprof%yy,xyprof%ny,xyprof%ncurve, + multix,xyprof%fname,xyprof%ylab,xyprof%cnames,0.,vp, + xyprof%log10,scfac,spval) if (xyprof%zptype(1:3)=='MID') then ! midpoints call yaxzpht(hts,dum,gcmlev,0,npress,dlev, + xyprof%yy,xyprof%ny,p0,xyprof%ihtyax,iyax,vp,spval) else ! interfaces call yaxzpht(hts,dum,gcmilev,0,npress,dlev, + xyprof%yy,xyprof%ny,p0,xyprof%ihtyax,iyax,vp,spval) endif call wrlab6(xyprof%tlabs,toffset,tlabchsz, + xyprof%blabs,boffset,blabchsz,vp,ilab_hq) return end subroutine pltxyloc !------------------------------------------------------------------- subroutine pltxyloc_amphase(xyprof,hts,iyax,vp,ut,ivpflag) use plt,only: igks_cgm,igks_ps,igks_x11, + iwk_cgm, iwk_ps, iwk_x11, yaxzpht use proc,only: spval,gcmlev,npress,dlev,p0 use input,only: ilab_hq ! ! Make xy line plot with field on x-axis, vertical on y-axis ! at selected location and ut (called from mkxyloc): ! ! Args: type(xyvprof),intent(in) :: xyprof real,intent(in) :: hts(:),ut real,intent(inout) :: vp(4) integer,intent(in) :: iyax,ivpflag ! ! Locals: integer :: iplt_cgm,iplt_ps,iplt_x11,ier,multix,i real :: dum(1,1) real :: toffset=.04, boffset=.15 real :: tlabchsz(3)=(/.02,.02,.02/), + blabchsz(3)=(/.02,.02,.02/) real :: vpsave(2) character(len=20) xname ! ! Begin exec: multix = 0 if (associated(xyprof%cnames)) then do i=1,size(xyprof%cnames) if (len_trim(xyprof%cnames(i)) > 0) multix = multix+1 enddo endif ! ! Draw to cgm workstation: iplt_cgm = 0 vpsave(1) = vp(1) vpsave(2) = vp(2) if(ivpflag == 1) then vp(1)=.05 vp(2)=.40 xname = 'AMPLITUDE' else vp(1)=.50 vp(2)=.85 xname = 'LT of MAX AMP' endif call pltmxy(xyprof%data,xyprof%yy,xyprof%ny,xyprof%ncurve, + multix,xname,xyprof%ylab,xyprof%cnames,0.,vp, + xyprof%log10,spval) if(ivpflag == 2) then call yaxzpht(hts,dum,gcmlev,0,npress,dlev, + xyprof%yy,xyprof%ny,p0,xyprof%ihtyax,iyax,vp,spval) endif vp(1) = vpsave(1) vp(2) = vpsave(2) if(ivpflag == 2) then call wrlab6(xyprof%tlabs,toffset,tlabchsz, + xyprof%blabs,boffset,blabchsz,vp,ilab_hq) endif return end subroutine pltxyloc_amphase !------------------------------------------------------------------- subroutine wrout_xyloc(xyprof,h,hcntr) use plt,only: zmin,zmax use proc use input ! ! Write ascii and/or xdr output data files for current xy loc frame: ! ! Args: ! (xyprof must be intent(inout) because this routine may need ! to call mkxyloclabs) ! type(xyvprof),intent(inout) :: xyprof type(history),intent(in) :: h,hcntr ! ! Locals: character(len=80) :: msgout real :: dum ! if (iwrdat==0.and.iwrxdr==0) return ! ! Make labels if not done for plots: ! xyprof%tlabs(2) label is full field name + units ! xyprof%tlabs(3) label is ut and grid info ! xyprof%blabs(2) label is hist vol info ! if (iplot == 0) call mkxyloclabs(xyprof,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,xyprof%data,1,xyprof%ny, + dum,xyprof%yy,xyprof%sname,xyprof%ylab,xyprof%tlabs(2), + xyprof%tlabs(3),xyprof%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 write(6,"('wrout_xyloc: xdr files NOT available')") endif return end subroutine wrout_xyloc end module mk_xylocs