! module mk_utvert use hist,only: history implicit none ! ! Make ut vs zp/ht contours at selected locations: ! (ut on x-axis, zp or height on y-axis) ! type utvert integer,pointer :: mtimes(:,:)=>NULL() ! model times (ut on x-axis) real,pointer :: + xx(:)=>NULL(),yy(:)=>NULL(), ! x and y coords (x is ut, y is zp or ht) + xslt(:)=>NULL(), ! local times corresponding to ut's + data(:,:)=>NULL() ! data(nx,ny) to be contoured integer :: nx,ny ! number of x and y coordinates integer :: ihtyax ! ht on yaxis if ihtyax==1, zp otherwise real :: glat,glon,slt ! lat, lon, and slt of current data character(len=3) :: loctype ! 'lon', 'slt', 'zm ', or 'gm ' character(len=56) :: fname ! full field name character(len=8) :: sname ! field short name character(len=16) :: funits ! field units character(len=16) :: ftype ! field type character(len=16) :: zptype ! midpoints or interfaces character(len=40) :: xlab,ylab ! x and y axis labels character(len=56) :: locname ! optional name of location integer :: log10 ! if > 0, then log10 of field is plotted character(len=240) :: + tlabs(3),blabs(3), ! top and bottom labels + hvol0,hvol1, ! 1st and last hist vols + hvol0_cntr,hvol1_cntr ! 1st and last hist vols (control) character(len=8) :: difftype ! raw or percent logical :: known ! true if field is known to the processor end type utvert real :: tlabchsz(3) = (/.02,.02,.02/), + blabchsz(3) = (/.018,.018,.018/) ! ! utv is utvert type to be plotted type(utvert),save :: utv ! ! utvdat(ntimes,npress,nloc,nf) is data from which utv is set up. ! utvdat(itime,:,:,:) is defined each time sub mkutvert is called. ! The last time mkutvert is called (itime==ntimes), utv is defined ! from utvdat, and contours are made. ! real,allocatable :: utvdat(:,:,:,:), utvdat_oh(:,:,:,:) real,allocatable :: utvdat_cntr(:,:,:,:), utvdat_cntr_oh(:,:,:,:) contains !------------------------------------------------------------------- subroutine mkutvert(f,fcntr,nf,h,hcntr,itime) ! ! Set up ut on x-axis, zp and/or ht on y-axis, at selected locations. ! (locations can include zonal and/or global means) ! (contouring is done by pltutvert, in this module) ! This routine is called at every model time (index itime), and the ! data is saved in utvdat(itime,:,:,:). ! Contours are made (from utvdat) only if the current call is for the ! last time (itime==ntimes) ! use proc use fields,only: field,nohalt use input use plt use set_vert ! ! Args: integer,intent(in) :: nf,itime type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: nloc,i,ixf,iloc,l,ifld,nfld,ier,iht,k,iyax,noh,ifoh, + izpht0,izpht1,nzprange,nhtscale,k0zp,k1zp,iadvfr,ixlat,ixlon,ix, + logint,isltax integer,save :: ixz_utvdat,ixne_utvdat real,pointer :: | htscale(:)=>NULL(), | zprange(:)=>NULL() ! user requested y-axes real :: vp(4),dum,boff real :: vputvert(4) = (/.12,.85,.29,.90/) real :: toffset=.06, boffset=.35 real,allocatable :: fdat1(:,:),fdat2(:,:) character(len=80) msgout integer :: ixne real,save :: fof2(100),hmf2(100),nmf2(100) real :: vl,vr,vb,vt,wl,wr,wb,wt integer :: lty real(kind=8) :: glbminteg,glbminteg_cntr real,parameter :: re = 6371.e5 ! earth radius (cm) ! ! Externals: real,external :: fslt,fmean,fglbm,quadrat integer,external :: ixfind,ixfindc logical,external :: isslt ! ! hscale and fhscale are defined by sethscale when field is height-only: ! fhscale(nfhscale) will be height scale of height-only fields. ! hscale(nhscale) will be height scale to use when plotting. ! 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 ! if (itime==ntimes) then if (.not.diffs) then write(6,"(/'Ut vs zp/ht contours at selected locations:')") else write(6,"(/'DIFFERENCE FIELDS: Ut vs zp/ht contours at ', + 'selected locations.')") endif endif ! ! Allocate utvdat if this is first call: if (itime==1) then nloc = 0 ! number of valid locations do i=1,mxloc if (utvert_locs(1,i)/=spval.and.utvert_locs(2,i)/=spval) + nloc = nloc+1 enddo if (nloc==0) then write(6,"('>>> mkutvert: need locations utvert_locs ', + '(nloc==0). Turning off ipltutvert.')") ipltutvert = 0 return endif nfld = 0 ! number of fields to allocate noh = 0 ! number of oh-v/b fields to allocate ixz_utvdat = 0 ! index to pert hts in utvdat ixne_utvdat = 0 ! index to NE in utvdat, for map_hmf2 do ixf=1,nf if (f(ixf)%requested.and.associated(f(ixf)%data).and. + f(ixf)%nlev/=1) then if (f(ixf)%vtype/='HEIGHT') then nfld = nfld+1 if (trim(f(ixf)%fname8)=='Z') ixz_utvdat = nfld else noh = noh+1 endif endif enddo if (nfld==0.and.noh==0) then write(6,"('>>> mkutvert: need fields: nfld==noh==0.', + ' Turning ipltutvert off.')") ipltutvert = 0 return endif ! ! if z not already requested, add extra field in utvdat to save ! total heights for later interpolation: if (ixz_utvdat==0) then nfld = nfld+1 ixz_utvdat = nfld endif ! ! Do the allocation: if (nfld > 0) then allocate(utvdat(ntimes,npress,nloc,nfld),stat=ier) if (ier /= 0) then write(6,"('mkutvert error allocating', + ' utvdat(ntimes,npress,nloc,nfld): ntimes=', + i4,' npress=',i3,' nloc=',i3,' nfld=',i3)") + ntimes,npress,nloc,nfld call allocerr(ier,"mkutvert allocating utvdat") endif if (diffs) then allocate(utvdat_cntr(ntimes,npress,nloc,nfld),stat=ier) if (ier /= 0) then write(6,"('mkutvert error allocating', + ' utvdat_cntr(ntimes,npress,nloc,nfld): ntimes=', + i4,' npress=',i3,' nloc=',i3,' nfld=',i3)") + ntimes,npress,nloc,nfld call allocerr(ier,"mkutvert allocating utvdat_cntr") endif endif endif if (noh > 0) then allocate(utvdat_oh(ntimes,nohalt,nloc,noh),stat=ier) if (ier /= 0) then write(6,"('mkutvert error allocating', + ' utvdat_oh(ntimes,nohalt,nloc,nfld): ntimes=', + i4,' nohalt=',i3,' nloc=',i3,' noh=',i3)") + ntimes,nohalt,nloc,noh call allocerr(ier,"mkutvert allocating utvdat_oh") endif if (diffs) then allocate(utvdat_cntr_oh(ntimes,nohalt,nloc,noh),stat=ier) if (ier /= 0) then write(6,"('mkutvert error allocating', + ' utvdat_cntr_oh(ntimes,nohalt,nloc,nfld): ntimes=', + i4,' nohalt=',i3,' nloc=',i3,' noh=',i3)") + ntimes,nohalt,nloc,noh call allocerr(ier,"mkutvert allocating utvdat_cntr_oh") endif endif endif ! ! Set up utv x-axis (utv%xx(itime) and utv%mtimes(:,itime) are ! defined at each call below): utv%nx = ntimes allocate(utv%xx(utv%nx),stat=ier) if (ier/=0) + call allocerr(ier,"mkutvert allocating utv%xx") utv%xlab = 'UT (HRS)' allocate(utv%mtimes(3,ntimes),stat=ier) if (ier/=0) + call allocerr(ier,"mkutvert allocating utv%mtimes") ! utv%hvol0 = h%mssvol ! if (diffs) utv%hvol0_cntr = hcntr%mssvol utv%hvol0 = h%histfile if (diffs) utv%hvol0_cntr = hcntr%histfile endif ! itime==1 ! ! Set up vertical scale and x-axis of utv if this is last time, i.e., ! will plot this time through the loops: if (itime==ntimes) then if (multiadvfr > 0) nppf = 0 endif ! ! Define current time in utv: utv%xx(itime) = float(h%mtime(1))*24.+h%ut utv%mtimes(:,itime) = h%mtime ! use pert for now ! ! Field loop: ifld = 0 ; ifoh = 0 fields_loop: do ixf=1,nf ! ! 10/10/07 btf: added vtype to this conditional to prevent attempt to ! plot height-independent fields fof2 and hmf2. For some unknown reason, ! nlev is > 1 here for fof2 and hmf2, however vtype is HT-INDEP. ! if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or. | f(ixf)%nlev==1.or.trim(f(ixf)%vtype)=='HT-INDEP') | cycle fields_loop if (trim(f(ixf)%vtype)/='HEIGHT') then ifld = ifld+1 ! field index to utvdat(:,:,:,ifld) else ifoh = ifoh+1 ! field index to utvdat_oh(:,:,:,ifoh) endif logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 ! ! If this is last time, define parts of utv and prepare for contouring: if (itime==ntimes) then ! ! Set up vertical scale(s) zp and/or ht: if (f(ixf)%zptype(1:3)=='INT') then ! interfaces call setvert(utvert_zprange,gcmilev,npress,dlev,nzprange, + k0zp,k1zp,zprange,utvert_htscale,nhtscale,htscale,spval,1) else ! midpoints call setvert(utvert_zprange,gcmlev,npress,dlev,nzprange, + k0zp,k1zp,zprange,utvert_htscale,nhtscale,htscale,spval,1) endif ! ! Set loop indices for zp and/or ht on y-axis (yax_loop): izpht0 = 0 izpht1 = 1 if (nhtscale == 0) izpht1 = 0 if (nzprange == 0) izpht0 = 1 utv%hvol1 = h%histfile if (diffs) utv%hvol1_cntr = hcntr%histfile ! utv%fname = f(ixf)%fname56 utv%sname = f(ixf)%fname8 utv%funits = f(ixf)%units utv%ftype = f(ixf)%type utv%zptype = f(ixf)%zptype utv%difftype = f(ixf)%difftype utv%known = f(ixf)%known ! ! 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 ! ! Set up height scale for fields that are plotted only in height: ! (i.e., oh-v or oh-b fields) if (trim(f(ixf)%vtype)=='HEIGHT'.and.nhtscale>0) + call sethscale(f(ixf),fhscale,nfhscale,htscale,nhtscale, + hscale,nhscale) endif ! itime==ntimes ! ! Location loop: iloc = 0 ! loc index to utvdat(:,:,iloc,:) loc_loop: do l=1,mxloc if (utvert_locs(1,l)==spval.or.utvert_locs(2,l)==spval) + cycle loc_loop iloc = iloc+1 call setloc(ixlat,ixlon) ! internal routine ! ! Save geopotential heights according to current field zptype. ! (sub setz is internal to mkutvert) ! call setz(utvdat(itime,:,iloc,ixz_utvdat), | utvdat_cntr(itime,:,iloc,ixz_utvdat),npress,ixlat,ixlon, | f(ixf)%zptype) ! ! Save fields in utvdat: if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ') then if (trim(f(ixf)%vtype)/='HEIGHT') then utvdat(itime,:,iloc,ifld) = f(ixf)%data(ixlon,ixlat,:) if (trim(f(ixf)%fname8)=='NE') ixne_utvdat = ifld else utvdat_oh(itime,:,iloc,ifoh) = f(ixf)%data(ixlon,ixlat,:) endif if (diffs) then if (trim(f(ixf)%vtype)/='HEIGHT') then utvdat_cntr(itime,:,iloc,ifld) = + fcntr(ixf)%data(ixlon,ixlat,:) else utvdat_cntr_oh(itime,:,iloc,ifoh) = + fcntr(ixf)%data(ixlon,ixlat,:) endif endif elseif (utv%loctype == 'zm ') then ! zonal means if (trim(f(ixf)%vtype)/='HEIGHT') then do k=1,npress utvdat(itime,k,iloc,ifld) = + fmean(f(ixf)%data(:,ixlat,k),nlon,spval,0) enddo else do k=1,nohalt utvdat_oh(itime,k,iloc,ifoh) = + fmean(f(ixf)%data(:,ixlat,k),nlon,spval,0) enddo endif if (diffs) then if (trim(f(ixf)%vtype)/='HEIGHT') then do k=1,npress utvdat_cntr(itime,k,iloc,ifld) = + fmean(fcntr(ixf)%data(:,ixlat,k),nlon,spval,0) enddo else do k=1,nohalt utvdat_cntr_oh(itime,k,iloc,ifoh) = + fmean(fcntr(ixf)%data(:,ixlat,k),nlon,spval,0) enddo endif endif else ! global means if (trim(f(ixf)%vtype)/='HEIGHT') then do k=1,npress utvdat(itime,k,iloc,ifld) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo ! ! 3/23/04: vertical integration of global means for Ray's dinos study: glbminteg =quadrat(utvdat(itime,:,iloc,ifld), | utvdat(itime,:,iloc,ixz_utvdat),npress) glbminteg = glbminteg*4.*pi*re**2 ! write(6,"('mkutvert: field ',a,' (pert) mtime=',3i4, ! | ' vertical integ of global means =',e12.4)") ! | f(ixf)%fname8,utv%mtimes(:,itime),glbminteg else do k=1,nohalt utvdat_oh(itime,k,iloc,ifoh) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif if (diffs) then if (trim(f(ixf)%vtype)/='HEIGHT') then do k=1,npress utvdat_cntr(itime,k,iloc,ifld) = + fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo glbminteg_cntr=quadrat(utvdat_cntr(itime,:,iloc,ifld), | utvdat_cntr(itime,:,iloc,ixz_utvdat),npress) glbminteg_cntr = glbminteg_cntr*4.*pi*re**2 ! write(6,"('mkutvert: field ',a,' (cntr) mtime=',3i4, ! | ' vertical integ of global means =',e12.4)") ! | f(ixf)%fname8,utv%mtimes(:,itime),glbminteg_cntr else do k=1,nohalt utvdat_cntr_oh(itime,k,iloc,ifoh) = + fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat,gcmlat, + dlat,dlon,spval) enddo endif endif endif if (itime==ntimes) utv%locname = utvert_locname(l) ! ! If not ready to plot, skip vertical axes loop: if (itime /= ntimes) cycle loc_loop ! ! Loop for vertical axes zp and/or ht: ! (This loop is executed only when itime==ntimes) yax_loop: do iht=izpht0,izpht1 utv%ihtyax = iht if (associated(utv%yy)) deallocate(utv%yy) if (associated(utv%data)) deallocate(utv%data) ! ! - - - - - - - - - - - - zp on y-axis: - - - - - - - - - - - - - - if (iht == 0) then ! zp on y-axis if (trim(f(ixf)%vtype)=='HEIGHT') cycle yax_loop if (iutvert_yaxright > 0) vputvert(2) = .77 ! ! Define y coords in pressure: utv%ny = nzprange allocate(utv%yy(utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%yy") utv%yy = zprange ! array op write(utv%ylab,"('LN(P0/P) (',a,')')") | trim(utv%zptype) ! ! Allocate and define utv data for zp on y-axis: allocate(utv%data(utv%nx,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%data") if (.not.diffs) then utv%data = utvdat(:,k0zp:k1zp,iloc,ifld) else allocate(fdat1(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating fdat1") fdat1 = utvdat(:,k0zp:k1zp,iloc,ifld) allocate(fdat2(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating fdat2") fdat2 = utvdat_cntr(:,k0zp:k1zp,iloc,ifld) call mkdiffs(fdat1,fdat2,utv%data,ntimes*utv%ny, + f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif ! ! - - - - - - - - - - - - ht on y-axis: - - - - - - - - - - - - - - else vputvert(2) = .86 ! ! Define y coords in height: write(utv%ylab,"('HEIGHT (KM)')") if (trim(f(ixf)%vtype)/='HEIGHT') then utv%ny = nhtscale allocate(utv%yy(utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%yy") utv%yy = htscale ! array op ! ! Allocate and define utv data for height on y-axis: allocate(utv%data(utv%nx,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%data") ! ! Interpolate to linear height scale: if (.not.diffs) then call cuthtint(utvdat(:,:,iloc,ifld), + utvdat(:,:,iloc,ixz_utvdat),ntimes,npress, + utv%data,htscale,nhtscale,logint,spval,ier,1) ! ! Save perturbed in fdat1, interpolate cntr and save in fdat2, ! and take diffs: else allocate(fdat1(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating fdat1") call cuthtint(utvdat(:,:,iloc,ifld), + 0.5*(utvdat(:,:,iloc,ixz_utvdat)+ + utvdat_cntr(:,:,iloc,ixz_utvdat)), + ntimes,npress,fdat1,htscale,nhtscale,logint,spval, + ier,1) allocate(fdat2(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating fdat2") call cuthtint(utvdat_cntr(:,:,iloc,ifld), + 0.5*(utvdat(:,:,iloc,ixz_utvdat)+ + utvdat_cntr(:,:,iloc,ixz_utvdat)), + ntimes,npress,fdat2,htscale,nhtscale,logint,spval, + ier,1) call mkdiffs(fdat1,fdat2,utv%data,ntimes*utv%ny, + f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif ! ! Field is calculated in height only -- plot only that part of ! user's requested height scale (htscale) that is within height ! range of the calculated field (fhscale). Also use resolution ! of the height-only field, overriding requested resolution ! (currently only oh-v and oh-b fields are height-only): ! else ! ht-only field (oh-v/b) utv%ny = nhscale allocate(utv%yy(utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%yy") utv%yy = hscale ! array op allocate(utv%data(utv%nx,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating utv%data") do k=1,nhscale if (hscale(k) < fhscale(1).or. + hscale(k) > fhscale(nfhscale)) then utv%data(:,k) = spval else ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev) utv%data(:,k) = utvdat_oh(:,ix,iloc,ifoh) endif enddo ! ! Difference fields of ht-only field: if (diffs) then allocate(fdat1(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkutvert allocating fdat1") fdat1 = utv%data do k=1,nhscale if (hscale(k) < fhscale(1).or. + hscale(k) > fhscale(nfhscale)) then fdat2(:,k) = spval else ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev) fdat2(:,k) = utvdat_cntr_oh(:,ix,iloc,ifoh) endif enddo call mkdiffs(fdat1,fdat2,utv%data,ntimes*utv%ny, + f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif endif ! height-only field endif ! - - - - - - - - - - - - end zp or ht on y-axis: - - - - - - - - - ! ! Take log10 if requested: if ((iutvert_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + iutvert_log10==2) then call log10f(utv%data,utv%nx*utv%ny,1.e-20,spval) utv%log10 = 1 else utv%log10 = 0 endif if (iplot > 0) then call setmultivp(vputvert,iadvfr,nppf+1,multiplt, + ipltrowcol,vp) ! ! Contour: isltax = 0 if(utv%loctype=='lon') isltax = 1 call pltutvert(utv,utvdat(:,:,iloc,ixz_utvdat), + iutvert_yaxright,isltax,vp) ! ! Add hmf2 line to plot, if requested if (map_hmf2 > 0) then ixne = ixfindc(h%fnames,h%nflds,'NE ') if (ixne==0) then write(6,"('WARNING mkderived: need NE and Z to ', + 'make HMF2')") cycle fields_loop endif call getset(vl,vr,vb,vt,wl,wr,wb,wt,lty) call set(vl,vr,vb,vt,1.,float(ntimes),wb,wt,lty) ! ! subroutine hnmf2(ht,ne,hmf2,nmf2,fof2,name,imx,kmx) ! real,intent(in),dimension(imx,kmx) :: ht,ne ! real,intent(out),dimension(imx) :: hmf2,nmf2,fof2 ! do i=1,ntimes ! call fof2int(utvdat(i,:,iloc,ixne_utvdat), ! + utvdat(i,:,iloc,ixz_utvdat), ! + npress,hmf2(i),fof2(i),0,1,1) call hnmf2(utvdat(i,:,iloc,ixz_utvdat), + utvdat(i,:,iloc,ixne_utvdat), + hmf2(i),nmf2(i),fof2(i),'HMF2',1,npress) call plch(float(i),hmf2(i),'*',.02,0.,0.) enddo endif ! ! Add top and bottom labels: call mkutvlabs(utv,h,hcntr,msgout) boff = boffset if (isltax <= 0) boff=.20 call wrlab6(utv%tlabs,toffset,tlabchsz, + utv%blabs,boff ,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, + 'ut vs vertical',iframe) ! ! If not making plots, wrout will need zmin,zmax. else call fminmax(utv%data,size(utv%data),zmin,zmax) endif ! iplot > 0 ! ! Write output data files: call wrout_utvert(utv,h,hcntr) ! ! End vertical axes loop: enddo yax_loop ! ! End location loop: enddo loc_loop ! ! End fields loop: 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 (itime==ntimes) then 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,' ','ut vs vertical', + iframe) ! ! Release allocated memory: if (allocated(utvdat)) deallocate(utvdat) if (allocated(utvdat_oh)) deallocate(utvdat_oh) if (associated(utv%data)) deallocate(utv%data) if (associated(utv%xx)) deallocate(utv%xx) if (associated(utv%yy)) deallocate(utv%yy) if (associated(utv%mtimes)) deallocate(utv%mtimes) if (associated(htscale)) deallocate(htscale) if (associated(zprange)) deallocate(zprange) if (associated(fhscale)) deallocate(fhscale) if (associated(hscale)) deallocate(hscale) endif return contains !------------------------------------------------------------------- subroutine setloc(ixlat,ixlon) ! internal to mkutvert ! ! Set location of utv, and return ixlat,ixlon: ! ! utvert_locs(1,iloc) is latitude, utvert_locs(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 utvert_locs(:,iloc)=='zm', (i.e., both lat,lon) then do global means) ! ! Args: integer,intent(out) :: ixlat,ixlon ! ! Locals: real :: slt ! utv%loctype = ' ' ! ! Global means: if (utvert_locs(1,l)==zmflag.and.utvert_locs(2,l)==zmflag)then utv%loctype = 'gm ' ! ! Zonal means: elseif (utvert_locs(2,l) == zmflag) then ! zonal means utv%loctype = 'zm ' utv%glon = zmflag ! ! Regular longitude or local time: else if (isslt(utvert_locs(2,l),slt)) then ! is local time utv%loctype = 'slt' utv%glon = fslt(slt,h%ut,dum,3) ixlon = ixfind(gcmlon,nlon,utv%glon,dlon) utv%glon = gcmlon(ixlon) utv%slt = slt else ! is regular longitude utv%loctype = 'lon' utv%glon = utvert_locs(2,l) ixlon = ixfind(gcmlon,nlon,utv%glon,dlon) utv%slt = fslt(dum,h%ut,utv%glon,1) endif endif ! ! Set latitude if not global means: if (utv%loctype /= 'gm ') then utv%glat = utvert_locs(1,l) ixlat = ixfind(gcmlat,nlat,utv%glat,dlat) endif return end subroutine setloc !----------------------------------------------------------------------- subroutine setz(z,zcntr,npress,ixlat,ixlon,zptype) ! ! Set geopotential heights z (utvdat(..ixz..)) according to ! zptype of current field (midpoints or interfaces): ! use fields,only: z_ifaces,z_midpts,zcntr_ifaces,zcntr_midpts integer,intent(in) :: npress,ixlat,ixlon real,dimension(npress),intent(out) :: z,zcntr character(len=*),intent(in) :: zptype ! ! Current field is at interfaces -> use Z at interfaces (z_ifaces): ! if (trim(zptype) == 'INTERFACES') then if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ')then z(:) = z_ifaces(ixlon,ixlat,:) if (diffs) zcntr(:) = zcntr_ifaces(ixlon,ixlat,:) elseif (utv%loctype == 'zm ') then do k=1,npress z(k) = fmean(z_ifaces(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zcntr(k) = fmean(zcntr_ifaces(:,ixlat,k),nlon,spval,0) enddo endif else do k=1,npress z(k) = fglbm(z_ifaces(:,:,k),nlon,nlat,gcmlat,dlat, | dlon,spval) enddo if (diffs) then do k=1,npress zcntr(k) = fglbm(zcntr_ifaces(:,:,k),nlon,nlat, | gcmlat,dlat,dlon,spval) enddo endif endif ! ! Current field is at midpoints -> use Z at midpoints (z_midpts): ! elseif (trim(zptype) == 'MIDPOINTS') then if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ')then z(:) = z_midpts(ixlon,ixlat,:) if (diffs) zcntr(:) = zcntr_midpts(ixlon,ixlat,:) elseif (utv%loctype == 'zm ') then do k=1,npress z(k) = fmean(z_midpts(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress zcntr(k) = fmean(zcntr_midpts(:,ixlat,k),nlon,spval,0) enddo endif else do k=1,npress z(k) = fglbm(z_midpts(:,:,k),nlon,nlat,gcmlat,dlat, | dlon,spval) enddo if (diffs) then do k=1,npress zcntr(k) = fglbm(zcntr_midpts(:,:,k),nlon,nlat, | gcmlat,dlat,dlon,spval) enddo endif endif else write(6,"('>>> mkutvert setz: unknown zptype=',a)") zptype endif end subroutine setz end subroutine mkutvert !------------------------------------------------------------------- subroutine mkutvlabs(utv,h,hcntr,msgout) use plt,only: zmin,zmax,ciu,scfac use proc,only: spval,diffs use input,only: ie5577,ie6300 implicit none ! ! Construct 3 top labels (utv%tlabs) and 3 bottom labels ! (utv%blabs) for lat slice: ! ! Args: type(utvert),intent(inout) :: utv type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Locals: integer :: lenlab ! ! Set up top and bottom text labels: ! ! Top (1st) top label is optionally provided by user: utv%tlabs(1) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if utv%log10 > 0) ! if (utv%known) then lenlab = len_trim(utv%fname)+len_trim(utv%funits)+3 if (utv%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(utv%tlabs(2))) then if (utv%log10 > 0) then utv%tlabs(2) = 'LOG10 '//trim(utv%fname)//' ('// + trim(utv%funits)//')' else utv%tlabs(2) = trim(utv%fname)//' ('// + trim(utv%funits)//')' endif else utv%tlabs(2) = trim(utv%fname) endif else ! unknown to proc lenlab = len_trim(utv%sname)+2+len_trim(utv%fname)+ | len_trim(utv%funits)+3 if (utv%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(utv%tlabs(2))) then if (utv%log10 > 0) then if (trim(utv%sname)/=trim(utv%fname)) then utv%tlabs(2) = 'LOG10 '//trim(utv%sname)//': '// | trim(utv%fname)//' ('//trim(utv%funits)//')' else utv%tlabs(2) = 'LOG10 '//trim(utv%fname)// | ' ('//trim(utv%funits)//')' endif else ! no log10 if (trim(utv%sname)/=trim(utv%fname)) then utv%tlabs(2) = trim(utv%sname)//': '// | trim(utv%fname)//' ('//trim(utv%funits)//')' else utv%tlabs(2) = trim(utv%fname)// | ' ('//trim(utv%funits)//')' endif endif else utv%tlabs(2) = trim(utv%fname) endif endif ! known or unknown ! ! Add DIFFS to label: if (diffs) then if (trim(utv%difftype)=='RAW') then if (len_trim(utv%tlabs(2))+7 <= len(utv%tlabs(2))) then utv%tlabs(2) = 'DIFFS: '//trim(utv%tlabs(2)) else utv%tlabs(2) = 'DIFFS: '//trim(utv%fname) endif else ! percent if (len_trim(utv%tlabs(2))+15 <= len(utv%tlabs(2))) then utv%tlabs(2) = trim(utv%difftype)//' DIFFS: '// | trim(utv%tlabs(2)) else utv%tlabs(2) = '% DIFFS: '//trim(utv%fname) endif endif endif ! ! Bottom (3rd) top label is ut and grid info: if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ') then if (utv%loctype == 'lon') then write(utv%tlabs(3),"(' LAT, LON=',f6.2,', ',f9.2)") + utv%glat,utv%glon else write(utv%tlabs(3),"(' LAT=',f6.2,' LOCAL TIME=',f5.2)") + utv%glat,utv%slt endif if (len_trim(utv%locname) > 0 .and. + len_trim(utv%tlabs(3))+len_trim(utv%locname)+3 <= + len(utv%tlabs(3))) + write(utv%tlabs(3)(len_trim(utv%tlabs(3))+1: + len_trim(utv%tlabs(3))+len_trim(utv%locname)+3), + "(' (',a,')')") + trim(utv%locname) elseif (utv%loctype == 'zm ') then write(utv%tlabs(3),"(' LAT=',f6.2,' (ZONAL MEANS)')") utv%glat else write(utv%tlabs(3),"('GLOBAL MEANS')") endif ! ! Top (1st) bottom label is min,max,interval: ! (zmin,zmax, and ciu are in module plt and were defined by contour) ! (scfac is in module plt) if (scfac==1.) then write(utv%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(utv%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac endif ! ! Middle (2nd) and bottom (3rd) bottom label are history info: ! if (.not.diffs) then ! if (len_trim(utv%hvol0)+len_trim(utv%hvol1)+15 <= ! + len(utv%blabs(2))) then ! write(utv%blabs(2),"('FIRST, LAST: ',a,', ',a)") ! + trim(utv%hvol0), trim(utv%hvol1) ! else ! write(utv%blabs(2),"('FIRST: ',a)") trim(utv%hvol0) ! endif ! utv%blabs(3) = ' ' ! if (trim(utv%ftype)=='EMISSION') then ! call mkemislab(utv%sname,ie5577,ie6300,utv%blabs(3)) ! if (len_trim(utv%blabs(3)) > 48) blabchsz(3) = .015 ! endif ! else ! if (len_trim(utv%hvol0)+len_trim(utv%hvol1)+15 <= ! + len(utv%blabs(2)).and. ! + len_trim(utv%hvol0_cntr)+len_trim(utv%hvol1_cntr)+15 <= ! + len(utv%blabs(3))) then ! write(utv%blabs(2),"('FIRST, LAST: ',a,', ',a,' MINUS ')") ! + trim(utv%hvol0), trim(utv%hvol1) ! write(utv%blabs(3),"('FIRST, LAST: ',a,', ',a)") ! + trim(utv%hvol0_cntr), trim(utv%hvol1_cntr) ! else ! write(utv%blabs(2),"('FIRST: ',a,' MINUS ')") trim(utv%hvol0) ! write(utv%blabs(3),"('FIRST: ',a)") trim(utv%hvol0_cntr) ! endif ! endif ! ! 5/05 btf: Put first and last history files in blabs(2 and 3): if (.not.diffs) then if (trim(utv%hvol0) /= trim(utv%hvol1)) then write(utv%blabs(2),"('FIRST: ',a)") trim(utv%hvol0) write(utv%blabs(3),"('LAST: ',a)") trim(utv%hvol1) else write(utv%blabs(2),"('MODEL ',a)") trim(h%version) write(utv%blabs(3),"(a)") trim(utv%hvol0) endif else ! diffs: show perturbed and control files: if (trim(utv%hvol0) /= trim(utv%hvol1)) then write(utv%blabs(2),"('PERT FIRST: ',a,' LAST: ',a)") | trim(utv%hvol0),trim(utv%hvol1) write(utv%blabs(3),"('CNTR FIRST: ',a,' LAST: ',a)") | trim(utv%hvol0_cntr),trim(utv%hvol1_cntr) else write(utv%blabs(2),"('PERT FILE: ',a)") trim(utv%hvol0) write(utv%blabs(3),"('CNTR FILE: ',a)") trim(utv%hvol0_cntr) endif endif ! ! Return message to be printed to stdout: if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ') then write(msgout,"(a,' lat,lon=',f5.1,',',f7.1,' min,max=', + 2e11.4,' iht=',i2)") utv%sname,utv%glat,utv%glon, + zmin,zmax,utv%ihtyax elseif (utv%loctype == 'zm ') then write(msgout,"(a,' lat=',f5.1,' zonal means min,max=', + 2e11.4,' iht=',i2)") utv%sname,utv%glat, + zmin,zmax,utv%ihtyax else write(msgout,"(a,' global means',9x,' min,max=', + 2e11.4,' iht=',i2)") utv%sname,zmin,zmax,utv%ihtyax endif return end subroutine mkutvlabs !------------------------------------------------------------------- subroutine pltutvert(utv,hts,iyax,isltax,vp) use plt use proc,only: gcmlev,gcmilev,npress,dlev,p0,spval ! ! Contour 2d ut vs zp/ht (ut on x-axis, ht or zp on y-axis) ! type(utvert) utv has been fully defined. ! ! Args: type(utvert),intent(in) :: utv real,intent(in) :: hts(:,:),vp(4) integer,intent(in) :: iyax integer,intent(inout) :: isltax ! ! Locals: integer :: iyaxright ! ! Set up conpack: call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',utv%xx(1)) call cpsetr('XCM',utv%xx(utv%nx)) call cpsetr('YC1',utv%yy(1)) call cpsetr('YCN',utv%yy(utv%ny)) call set(vp(1),vp(2),vp(3),vp(4), + utv%xx(1),utv%xx(utv%nx),utv%yy(1),utv%yy(utv%ny),1) ! ! Contour: call contour(utv%data,utv%nx,utv%nx,utv%ny) ! ! Add axes and labels: ! isltax = 0 ! if (utv%loctype=='lon') isltax = 1 ! labutxy may change isltax: call labutxy(utv%mtimes,utv%nx,utv%yy,utv%ny,utv%ylab,0., + isltax,utv%glon) ! ! Add extra right hand axis (in pressure if height is on yaxis, ! or in height if pressure is on yaxis): iyaxright = iyax-iyax/10*10 if (iyaxright > 0) then if (utv%zptype(1:3)=='INT') then call yaxzpht(gcmilev,hts,gcmilev,utv%nx,npress,dlev, + utv%yy,utv%ny,p0,utv%ihtyax,iyaxright,vp,spval) else call yaxzpht(gcmlev,hts,gcmlev,utv%nx,npress,dlev, + utv%yy,utv%ny,p0,utv%ihtyax,iyaxright,vp,spval) endif endif return end subroutine pltutvert !------------------------------------------------------------------- subroutine wrout_utvert(utv,h,hcntr) use plt,only: zmin,zmax use proc use input ! ! Write ascii and/or xdr output data files for current lat slice frame: ! ! Args: ! (utv must be intent(inout) because this routine may need ! to call mklatlabs) ! type(utvert),intent(inout) :: utv 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: if (iplot == 0) call mkutvlabs(utv,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,utv%data,utv%nx,utv%ny, + utv%xx,utv%yy,utv%xlab,utv%ylab,utv%tlabs(2), + utv%tlabs(3),utv%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: ! subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1, ! + lab2,lab3,lab4,mtime,iclose) ! if (iwrxdr > 0) then call wrxdr(flnm_xdr,utv%data,utv%nx,utv%ny, + utv%xx,utv%yy,utv%xlab,utv%ylab, + utv%tlabs(2),utv%tlabs(3),utv%blabs(2), + utv%blabs(1),h%mtime,0) flnm_xdr = flnm_xdr(1:len_trim(flnm_xdr)-1) 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_utvert end module mk_utvert