! module mk_xyut use hist,only: history use input,only: flnm_cdf_xyut implicit none ! ! Make xy ut vs field line plots at selected locations, and selected ! zp or ht (or ht-indep, ht-integ, etc). ! Ut on x-axis, field(s) on y-axis. More than one curve (field) ! may be drawn on each plot (at same scale on y-axis), i.e., ncurve > 1. ! (ncurve > 1 for oh-b and oh-v) ! type xyut_type integer,pointer :: | mtimes(:,:)=>NULL() ! model times (ut on x-axis) real :: glat,glon,slt ! geog lat,lon,slt of profile character(len=3) :: loctype ! 'lon', 'slt', 'zm ', or 'gm ' character(len=56) :: locname ! optional location name real :: zp,ht ! zp or ht of data character(len=5) :: levtype ! 'zp','ht','indep','integ','hmf2' integer :: iloc,ilev ! indices to data at current loc,lev integer :: nx ! number of points real,pointer :: xx(:)=>NULL() ! x coords (increasing ut) character(len=40) :: xlab ! x axis labels real,pointer :: data(:,:)=>NULL() ! field data (ntimes,ncurve) integer :: ncurve ! number of curves (fields) on single plot 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 character(len=16) :: ftype ! field type character(len=16) :: zptype ! zp type ("MIDPOINTS" or "INTERFACES") logical :: known ! true if field is known to the processor integer :: log10 ! if > 0, then log10 of field is plotted character(len=1024) :: + 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 integer :: nslev,nsloc ! number of selected levels and locs end type xyut_type real :: tlabchsz(3)=(/.02,.02,.02/), + blabchsz(3)=(/.02,.015,.015/) ! ! xyut is structure to be plotted after last time is ! saved in xyutdat: ! type(xyut_type),save :: xyut ! ! xyutdat(ntimes,nloc,nlev,nf) is data from which xyut is set up. ! (nloc is from user xyut_locs, and nlev is npress unless height-only ! fields such as oh-v and oh-b) ! xyutdat(itime,:,:,:) is defined each time sub mkxyut is called. ! The last time mkxyut is called (itime==ntimes), xyut is defined ! from xyutdat, and line plots are made. ! xyutdat_hmf2(ntimes,nloc) is hmf2, saved like xyutdat, for use ! when user requests xyut_zpht='hmf2', in which case fields are ! plotted at hmf2 heights. ! xyutdat_ohv and xyutdat_ohb are same as xyutdat, but for oh-v and ! oh-b. These must be separate, as number of levels allocated is ! nohvalt instead of npress (oh are height-only fields) ! xyutdat_emisdop and xyutdat_ohbdop are doppler fields. ! real,allocatable :: + xyutdat(:,:,:,:), ! (ntimes,nloc,nlev,nflds) + xyutdat_ohv(:,:,:,:), ! (ntimes,nloc,nohalt,nohv) + xyutdat_ohb(:,:,:,:), ! (ntimes,nloc,nohalt,nohb) + xyutdat_hmf2(:,:), ! (ntimes,nloc) + xyutdat_emisdop(:,:,:,:), ! (ntimes,nloc,nemis,nfdop) + xyutdat_ohbdop(:,:,:,:), ! (ntimes,nloc,nohb,nfdop) + udop(:),vdop(:) ! doppler winds (for 1d vectors) ! ! Similiar to xyutdat above, but control case for diffs: real,allocatable :: + xyutdat_cntr(:,:,:,:), ! (ntimes,nloc,nlev,nflds) + xyutdat_cntr_ohv(:,:,:,:), ! (ntimes,nloc,nohalt,nohv) + xyutdat_cntr_ohb(:,:,:,:), ! (ntimes,nloc,nohalt,nohb) + xyutdat_cntr_hmf2(:,:), ! (ntimes,nloc) + xyutdat_cntr_emisdop(:,:,:,:), ! (ntimes,nloc,nemis,nfdop) + xyutdat_cntr_ohbdop(:,:,:,:) ! (ntimes,nloc,nohb,nfdop) ! ! Local filename of lon slice netcdf file: ! character(len=16) :: flnm_cdf_xyut = 'xyut.nc' ! integer :: iadvfr contains !------------------------------------------------------------------- subroutine mkxyut(f,fcntr,nf,h,hcntr,itime) ! ! Set up ut on x-axis, field on y-axis at selected locations and ! vertical levels. ! (locations can include zonal and/or global means) ! (vertical level may be zp, ht, ht-independent, ht-integrated, or hmf2) ! This routine is called at every model time (index itime), and the ! data is saved in utxydat(itime,:,:,:). ! XY plots are made (from utxydat) only if the current call is for the ! last time (itime==ntimes) ! use proc use fields use input use plt,only: igks_cgm,igks_ps,igks_x11,iwk_cgm, iwk_ps, iwk_x11, + advframe,scfac ! ! Args: integer,intent(in) :: nf,itime type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: + i,ier,ixf,ifld,l,iloc,k,ilev,ixz,izp,ixlat,ixlon, + loginterp,it,idum,ixemis integer,save :: ixz_xyutdat,nslev,nsloc,ixhmf2,nohv,nohb real :: fmin,fmax,tmp real,allocatable :: futdat1(:),futdat2(:) character(len=80) msgout real(kind=8) :: glbminteg,glbminteg_cntr ! global mean vertical integration real :: earth_area ! surface area of earth real,parameter :: re = 6371.e5 ! earth radius (cm) ! ! Externals: real,external :: fmean,fglbm,quadrat integer,external :: ixfind,ixfindc character(len=16),external :: float_to_str interface function vecsum(u,v,id1,id2,spv) implicit none integer,intent(in) :: id1,id2 real,intent(in) :: u(id1,id2),v(id1,id2) real:: vecsum(id1,id2) ! result variable real,intent(in),optional :: spv end function vecsum end interface ! earth_area = 4.*pi*re**2 ! surface area of earth for vertical integration if (itime==ntimes) then if (.not.diffs) then write(6,"(/'Ut vs field line plots at selected locations:')") else write(6,"(/'DIFFERENCE FIELDS: Ut vs field line plots at ', + 'selected locations.')") endif if (multiadvfr > 0) nppf = 0 endif ! ! Allocate needed arrays in which to save data from each history: ! if (itime==1) then ! ! Allocate xyutdat_ohb and/or xyutdat_ohv (height-only fields): nohv = 0 ; nohb = 0 do ixf=1,nf if (trim(f(ixf)%type)=='OH-VIB') nohv = nohv+1 if (trim(f(ixf)%type)=='OH-BAND') nohb = nohb+1 enddo if (nohv>0) call alloc_xyutdat(f,nf,idum,'OH-VIB ', | xyut,ier) if (nohb>0) call alloc_xyutdat(f,nf,idum,'OH-BAND ', | xyut,ier) ! ! Allocate xyutdat(ntimes,nloc,npress,nf): ! (all fields except oh, which is height-only) ! (alloc_xyutdat will also allocate xyutdat_emisdop and/or xyutdat_ohbdop ! if necessary) (ixz_xyutdat is returned by alloc_xyutdat) ! call alloc_xyutdat(f,nf,ixz_xyutdat,'noheight ',xyut,ier) if (ier/=0) then if (ier<0) then write(6,"('>>> mkxyut: error allocating xyutdat -- ', + 'turning off ipltxyut.')") ipltxyut = 0 return else ! ier>0 means no fields were found if (nohv==0.and.nohb==0) then ipltxyut = 0 return else write(6,"('mkxyut: processing only oh fields..')") goto 200 endif endif endif ! ! If a selected zpht is "hmf2", then we need to plot fields at hmf2 ! heights. If this is requested, allocate xyutdat_hmf2(ntimes,nloc), ! where hmf2 will be saved. Sub alloc_xyutdat_hmf2 is internal, and ! will also save ixhmf2 as index to HMF2 if hmf2 heights are needed. ! If doing diffs, alloc_xyutdat_hmf2 will also allocate xyutdat_cntr_hmf2. ! If hmf2 is not needed, alloc_xyutdat_hmf2 leaves ixhmf2=0 ! call alloc_xyutdat_hmf2(f,nf,ixhmf2,ier) if (ier/=0) then write(6,"('>>> WARNING mkxyut: error allocating ', + 'xyutdat_hmf2')") endif ! ! Set up xyut x-axis (xyut%xx(itime) and xyut%mtimes(:,itime) are ! defined at each call below): xyut%nx = ntimes allocate(xyut%xx(xyut%nx),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut allocating xyut%xx") xyut%xlab = 'UT (HRS)' allocate(xyut%mtimes(3,ntimes),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut allocating xyut%mtimes") ! ! nslev = number of selected levels (not including ht-integ) nslev = 0 do k=1,mxzpht if (xyut_zpht(k)/=spval) nslev = nslev+1 enddo ! ! nsloc = number of selected locations: nsloc = 0 do k=1,mxloc if (xyut_locs(1,k)/=spval.and.xyut_locs(2,k)/=spval) | nsloc = nsloc+1 enddo endif ! itime==1 if (nslev <= 0) then write(6,"(/,'>>> mkxyut: nslev=',i3,': no xyut plots will ', | 'be made -- check user input xyut_zpht.')") nslev endif xyut%nslev = nslev if (nsloc <= 0) then write(6,"(/,'>>> mkxyut: nsloc=',i3,': no xyut plots will ', | 'be made -- check user input xyut_locs.')") nsloc endif xyut%nsloc = nsloc ! ! Define 1st and last history volumes: if (itime==1) then xyut%hvol0 = h%histfile if (diffs) xyut%hvol0_cntr = hcntr%histfile elseif (itime==ntimes) then xyut%hvol1 = h%histfile if (diffs) xyut%hvol1_cntr = hcntr%histfile endif ! ! Save heights in xyutdat(:,:,:,ixz_xyutdat) ! (if ixhmf2 > 0, then hmf2 heights are saved in xyutdat_hmf2) ! (if doing diffs, also save xyutdat_cntr and xyutdat_cntr_hmf2) ! ! call savehts_xyut(f,fcntr,nf,h,xyut,itime,ixz_xyutdat,ixhmf2) ! ! Define current time in xyut: xyut%xx(itime) = float(h%mtime(1))*24.+h%ut xyut%mtimes(:,itime) = h%mtime ! use pert for now ! ! Write history parameters to netcdf file (no xyut data yet): if (len_trim(sendcdf_xyut) > 0.or.len_trim(flnm_cdf_xyut) > 0) | call wrcdf_xyut(xyut,h,hcntr,itime,.true.) ! ! Field loop: ifld = 0 ; ixemis = 0 fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or. + trim(f(ixf)%vtype)=='HEIGHT') cycle fields_loop ifld = ifld+1 ! field index to xyutdat(:,:,:,ifld) if (any(f(ixf)%fname8==doppler_emis)) ixemis = ixemis+1 ! ! Allocate xyut%data, and define parts of xyut, if last time: if (itime==ntimes) then xyut%ncurve = 1 ! multiple curves are only for OH-V/B ! (see sub mkxyut_oh below) if (associated(xyut%data)) deallocate(xyut%data) allocate(xyut%data(ntimes,xyut%ncurve),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut allocating xyut%data") loginterp = 0 if (trim(f(ixf)%type)=='DENSITY') loginterp = 1 ! ! xyut%log10 = 0 -> all fields linear ! xyut%log10 = 1 -> log10 of densities ! xyut%log10 = 2 -> log10 of all fields xyut%log10 = 0 if (ixyut_log10==1.and.trim(f(ixf)%type)=='DENSITY') + xyut%log10 = ixyloc_log10 if (ixyut_log10 > 1) xyut%log10 = 1 endif ! ! Scale current field's data if fscale is requested ! (scfac is in plt.f, s.a. sub pltxmy): scfac = f(ixf)%scalefac ! ! Save heights in xyutdat(:,:,:,ixz_xyutdat) ! (if ixhmf2 > 0, then hmf2 heights are saved in xyutdat_hmf2) ! (if doing diffs, also save xyutdat_cntr and xyutdat_cntr_hmf2) ! call savehts_xyut(f,fcntr,nf,h,xyut,itime,ixz_xyutdat,ixhmf2, | f(ixf)%zptype) ! ! Location loop: iloc = 0 loc_loop: do l=1,mxloc if (xyut_locs(1,l)==spval.or.xyut_locs(2,l)==spval) + cycle loc_loop iloc = iloc+1 xyut%iloc = iloc ! ! Set loctype (lon,slt,zm,gm), and related indices or mean flags: call setloctype(xyut_locs(1,l),xyut_locs(2,l),h%ut, + xyut%glat,xyut%glon,xyut%slt,ixlat,ixlon,xyut%loctype) ! ! Define xyutdat(itime,iloc,:,ifld) at current loc: ! (note, use k=1,f(ixf)%nlev, instead of npress, for ht-indep fields) ! ! Save at regular longitude, latitude: if (xyut%loctype=='lon'.or.xyut%loctype=='slt') then do k=1,f(ixf)%nlev xyutdat(itime,iloc,k,ifld) = f(ixf)%data(ixlon,ixlat,k) enddo if (diffs) then do k=1,f(ixf)%nlev xyutdat_cntr(itime,iloc,k,ifld) = + fcntr(ixf)%data(ixlon,ixlat,k) enddo endif ! ! Save zonal means: elseif (xyut%loctype== 'zm ') then do k=1,f(ixf)%nlev xyutdat(itime,iloc,k,ifld) = + fmean(f(ixf)%data(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,f(ixf)%nlev xyutdat_cntr(itime,iloc,k,ifld) = + fmean(fcntr(ixf)%data(:,ixlat,k),nlon,spval,0) enddo endif ! ! Save global means: elseif (xyut%loctype=='gm ') then do k=1,f(ixf)%nlev xyutdat(itime,iloc,k,ifld) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat,dlon, + spval) enddo glbminteg = quadrat(xyutdat(itime,iloc,:,ifld), | xyutdat(itime,iloc,:,ixz_xyutdat), | size(xyutdat,3)) ! write(6,"('Field ',a,' mtime=',3i4, ! | ' glbmean vertical integ (pert) = ',e12.4, ! | ' total=',e12.4)") f(ixf)%fname8,xyut%mtimes(:,itime), ! | glbminteg,glbminteg*earth_area if (diffs) then do k=1,f(ixf)%nlev xyutdat_cntr(itime,iloc,k,ifld) = + fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo glbminteg_cntr = | quadrat(xyutdat_cntr(itime,iloc,:,ifld), | xyutdat_cntr(itime,iloc,:,ixz_xyutdat), | size(xyutdat_cntr,3)) ! write(6,"('Field ',a,' mtime=',3i4, ! | ' glbmean vertical integ (cntr) = ',e12.4, ! | ' total=',e12.4)") f(ixf)%fname8,xyut%mtimes(:,itime), ! | glbminteg_cntr,glbminteg_cntr*earth_area endif else write(6,"('>>> WARNING mkxyut: unknown loctype=',a)") + xyut%loctype endif ! ! Define doppler fields in xyutdat_emisdop: if (ixyut_doppler>0.and.ixemis>0.and. + trim(f(ixf)%type)=='EMISSION') then call savedop_xyut(f,nf,f(ixf)%fname8,f(ixf)%type, + xyutdat(itime,iloc,:,ifld),f(ixf)%nlev,ixlat,ixlon, + xyutdat_emisdop(itime,iloc,ixemis,:)) if (diffs) then call savedop_xyut(fcntr,nf,fcntr(ixf)%fname8, + fcntr(ixf)%type,xyutdat_cntr(itime,iloc,:,ifld), + fcntr(ixf)%nlev,ixlat,ixlon, + xyutdat_cntr_emisdop(itime,iloc,ixemis,:)) endif endif ! ! Set location name, if given by user: if (itime==ntimes) then xyut%locname = "" if (len_trim(xyut_locname(l))>0) + xyut%locname = ' ('//trim(xyut_locname(l))//')' ! ! Names are set here instead of in field loop above in case ! we have just finished doppler fields, and need to continue ! with emissions at new location (names need to be reset from ! doppler back to emissions): ! ! if (.not.diffs) then ! xyut%fname = f(ixf)%fname56 ! else if(trim(f(ixf)%difftype)=='RAW') then ! xyut%fname = 'DIFFS: '// f(ixf)%fname56 ! else ! xyut%fname = trim(f(ixf)%difftype)//' DIFFS: '// ! + f(ixf)%fname56 ! endif xyut%fname = f(ixf)%fname56 xyut%sname = f(ixf)%fname8 xyut%funits = f(ixf)%units if (f(ixf)%scalefac /= 1.) | write(xyut%funits,"(a,' * ',1pe9.2)") trim(f(ixf)%units), | f(ixf)%scalefac xyut%ftype = f(ixf)%type xyut%zptype = f(ixf)%zptype xyut%known = f(ixf)%known xyut%difftype = f(ixf)%difftype endif ! ! Enter selected levels loop only if this is last time: if (itime < ntimes) cycle loc_loop ! ! Levels loop (xyutdat now fully defined): ! Now define structure xyut%data(:,ncurve) and plot,etc. ! character(len=5) :: levtype ! 'zp','ht','indep','integ','hmf2' ! ilev = 0 lev_loop: do k=1,mxzpht if (((xyut_zpht(k)==spval).and. + .not.(nslev==0.and.trim(f(ixf)%type)=='EMISSION')) .or. + (ilev>0.and.trim(f(ixf)%vtype)=='HT-INDEP')) + cycle lev_loop ilev = ilev+1 xyut%ilev = ilev if (trim(f(ixf)%vtype)=='HT-INDEP') then xyut%levtype='indep' elseif (trim(float_to_str(xyut_zpht(k)))=='hmf2'.or. + trim(float_to_str(xyut_zpht(k)))=='HMF2'.or. | xyut_zpht(k)==hmf2flag) then xyut%levtype = 'hmf2 ' elseif (xyut_zpht(k)==spval) then ! + (trim(f(ixf)%type)=='EMISSION'.and.ilev==nslev)) then xyut%levtype = 'integ' elseif (xyut_zpht(k) > gcmlev(npress)) then xyut%levtype = 'ht ' else xyut%levtype = 'zp ' endif 100 continue select case (xyut%levtype) ! ! Define data at a selected zp: case ('zp ') izp = ixfind(gcmlev,npress,xyut_zpht(k),dlev) if (izp<=0) then write(6,"('>>> mkxyloc: bad zp=',e12.4)") xyut_zpht(k) cycle lev_loop endif xyut%zp = gcmlev(izp) if (.not.diffs) then xyut%data(:,1) = xyutdat(:,iloc,izp,ifld) else call mkdiffs(xyutdat(:,iloc,izp,ifld), + xyutdat_cntr(:,iloc,izp,ifld),xyut%data(:,1), + ntimes,f(ixf)%difftype) endif ! ! Define data at a selected ht: case ('ht ') xyut%ht = xyut_zpht(k) if (.not.diffs) then do it=1,ntimes call intloc(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),npress,xyut_zpht(k), + 1,loginterp,xyut%data(it,1),1,1,1,ier,spval,0) enddo else allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat2") do it=1,ntimes ! ! 4/98: 1st 2 calls use perturbed heights for perturbed interpolation, ! and control heights for control interpolation. The 2nd 2 calls use ! the average 0.5*(cntr+pert) heights for both interpolations. ! These 2 approaches yield quite different results, esp below ~100km. ! ! 7/18/11: Go back to using perturbed heights for perturbed interp, ! and control heights for control interp before taking diffs. ! This makes more sense when comparing full fields with difference fields. ! call intloc(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),npress,xyut_zpht(k), + 1,loginterp,futdat1(it),1,1,1,ier,spval,0) call intloc(xyutdat_cntr(it,iloc,:,ifld), + xyutdat_cntr(it,iloc,:,ixz_xyutdat),npress, + xyut_zpht(k),1,loginterp,futdat2(it),1,1,1,ier, + spval,0) ! call intloc(xyutdat(it,iloc,:,ifld), ! + 0.5*(xyutdat(it,iloc,:,ixz_xyutdat)+ ! + xyutdat_cntr(it,iloc,:,ixz_xyutdat)),npress, ! + xyut_zpht(k),1,loginterp,futdat1(it),1,1,1,ier, ! + spval,0) ! call intloc(xyutdat_cntr(it,iloc,:,ifld), ! + 0.5*(xyutdat(it,iloc,:,ixz_xyutdat)+ ! + xyutdat_cntr(it,iloc,:,ixz_xyutdat)),npress, ! + xyut_zpht(k),1,loginterp,futdat2(it),1,1,1,ier, ! + spval,0) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,1),ntimes, + f(ixf)%difftype) deallocate(futdat1) ; deallocate(futdat2) endif ! ! Define height-independent data: case ('indep') if (.not.diffs) then xyut%data(:,1) = xyutdat(:,iloc,1,ifld) else call mkdiffs(xyutdat(:,iloc,1,ifld), + xyutdat_cntr(:,iloc,1,ifld),xyut%data(:,1),ntimes, + f(ixf)%difftype) endif ! ! Interpolate to hmf2 heights: case ('hmf2 ') if (.not.diffs) then do it=1,ntimes call intloc(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),npress, + xyutdat_hmf2(it,iloc),1,loginterp, + xyut%data(it,1),1,1,1,ier,spval,0) enddo else allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat2") do it=1,ntimes call intloc(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),npress, + xyutdat_hmf2(it,iloc),1,loginterp, + futdat1(it),1,1,1,ier,spval,0) call intloc(xyutdat_cntr(it,iloc,:,ifld), + xyutdat_cntr(it,iloc,:,ixz_xyutdat),npress, + xyutdat_cntr_hmf2(it,iloc),1,loginterp, + futdat2(it),1,1,1,ier,spval,0) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,1),ntimes, + f(ixf)%difftype) deallocate(futdat1) ; deallocate(futdat2) endif ! ! Vertical integration (emission field): case ('integ') if (.not.diffs) then do it=1,ntimes xyut%data(it,1) = quadrat(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),size(xyutdat,3)) enddo else allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat2") do it=1,ntimes futdat1(it) = quadrat(xyutdat(it,iloc,:,ifld), + xyutdat(it,iloc,:,ixz_xyutdat),size(xyutdat,3)) enddo do it=1,ntimes futdat2(it) = quadrat(xyutdat_cntr(it,iloc,:,ifld), + xyutdat_cntr(it,iloc,:,ixz_xyutdat), + size(xyutdat_cntr,3)) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,1),ntimes, + f(ixf)%difftype) deallocate(futdat1) ; deallocate(futdat2) endif xyut%data = xyut%data * 1.e-6 ! rayleighs xyut%funits = "RAYLEIGHS" case default write(6,"('>>> WARNING: unknown xyut%levtype=',a)") + xyut%levtype end select ! xyut%levtype ! ! doxyut plots and writes output data: ! write(6,"('mkxyut: ',a,' (',a,') iloc=',i3,' ilev=',i3, ! | ' xyut%data=',/,(6e12.4))") f(ixf)%fname8, ! | trim(f(ixf)%zptype),iloc,ilev,xyut%data(:,1) call doxyut(xyut,h,hcntr,itime) ! ! Make doppler t,u,v from ht-integrated emission: ! (xyutdat_emisdop was defined by savedop_xyut) if (xyut%levtype == 'integ') then if (ixyut_doppler>0.and.any(f(ixf)%fname8==doppler_emis)) + then if (diffs) then allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut allocating futdat2") endif do i=1,nfdop do it=1,ntimes xyut%data(it,1) = xyutdat_emisdop(it,iloc,ixemis,i) enddo ! ! Diffs of doppler emission: if (diffs) then futdat1 = xyut%data(:,1) do it=1,ntimes futdat2(it) = + xyutdat_cntr_emisdop(it,iloc,ixemis,i) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,1),ntimes, + f(ixf)%difftype) write(xyut%fname,"(a,' DIFFS: ',a,' ',a)") + trim(f(ixf)%difftype),trim(f(ixf)%fname8), + trim(doppler_fnames(i)) else write(xyut%fname,"(a,' ',a)") + trim(f(ixf)%fname8),trim(doppler_fnames(i)) endif xyut%sname = doppler_snames(i) xyut%funits = doppler_units(i) call doxyut(xyut,h,hcntr,itime) ! ! Doppler u+v 1d vector plot: if (i==2) then ! load u allocate(udop(size(xyut%data,1)),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut allocating udop") udop = xyut%data(:,1) elseif (i==3) then ! load v and make plot allocate(vdop(size(xyut%data,1)),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut allocating vdop") vdop = xyut%data(:,1) xyut%data(:,1:1) = vecsum(udop,vdop,xyut%nx,1,spval) xyut%sname = 'DOP(U+V)' xyut%funits = doppler_units(i) if (.not.diffs) then write(xyut%fname,"(a,' DOPPLER U+V')") + trim(f(ixf)%fname8) else write(xyut%fname,"(a,' DIFFS: ',a, + ' DOPPLER U+V')") + trim(f(ixf)%difftype),trim(f(ixf)%fname8) endif call doxyut(xyut,h,hcntr,itime,udop,vdop) deallocate(udop) deallocate(vdop) endif enddo ! i=1,nfdop if (diffs) then deallocate(futdat1) deallocate(futdat2) endif endif ! do doppler exit lev_loop ! just did integ and/or doppler, so exit levs endif ! integ ! ! If emission was just plotted and this is last selected level, ! go back up and do ht-integrations (integht is user input): if (integht > 0.and.trim(f(ixf)%type)=='EMISSION'.and. + ilev==nslev) then xyut%levtype = 'integ' goto 100 endif enddo lev_loop enddo loc_loop enddo fields_loop ! ! Plot multi-curve oh-v and/or oh-b: 200 continue if (nohv>0) call mkxyut_oh('OH-VIB ',xyutdat_ohv, + xyutdat_cntr_ohv) if (nohb>0) call mkxyut_oh('OH-BAND ',xyutdat_ohb, + xyutdat_cntr_ohb) ! if (itime==ntimes) then if (multiplt > 0 .and. multiadvfr == 1 .and. iadvfr <= 0) + call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,1,nppf,' ', + 'xy vert profs at locs',iframe) ! ! Release allocated memory: ! ! Allocatable arrays: if (allocated(xyutdat)) deallocate(xyutdat) if (allocated(xyutdat_ohv)) deallocate(xyutdat_ohv) if (allocated(xyutdat_ohb)) deallocate(xyutdat_ohb) if (allocated(xyutdat_hmf2)) deallocate(xyutdat_hmf2) if (allocated(xyutdat_emisdop))deallocate(xyutdat_emisdop) if (allocated(xyutdat_ohbdop)) deallocate(xyutdat_ohbdop) ! if (allocated(xyutdat_cntr)) deallocate(xyutdat_cntr) if (allocated(xyutdat_cntr_ohv)) deallocate(xyutdat_cntr_ohv) if (allocated(xyutdat_cntr_ohb)) deallocate(xyutdat_cntr_ohb) if (allocated(xyutdat_cntr_hmf2)) deallocate(xyutdat_cntr_hmf2) if (allocated(xyutdat_cntr_emisdop)) + deallocate(xyutdat_cntr_emisdop) if (allocated(xyutdat_cntr_ohbdop)) + deallocate(xyutdat_cntr_ohbdop) ! ! Pointers: if (associated(xyut%data)) deallocate(xyut%data) if (associated(xyut%xx)) deallocate(xyut%xx) if (associated(xyut%mtimes)) deallocate(xyut%mtimes) if (associated(xyut%cnames)) deallocate(xyut%cnames) endif return contains !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine mkxyut_oh(type,xyutdat_oh,xyutdat_cntr_oh) ! ! Process oh fields (multiple curves per plot) (internal to mkxyut): ! use proc,only: pi use fields,only: field,nohalt,oh_alt use input,only: ibohv_watts ! ! Args: real,intent(inout) :: xyutdat_oh(:,:,:,:),xyutdat_cntr_oh(:,:,:,:) character(len=16),intent(in) :: type ! ! Locals: integer :: ixf,ifld,l,iloc,ixlat,ixlon,ntype,k,i,hi,lo,len integer,save :: nslev real :: fmin,fmax character(len=8) :: difftype ! ! Externals: integer,external :: ixohband ! ! Validate type: if (trim(type)/='OH-VIB'.and.trim(type)/='OH-BAND') then write(6,"('>>> mkxyut_oh: unknown type=',a)") type return endif if (trim(type)=='OH-VIB') then ntype = nohv else ntype = nohb endif ! ! First time: ! Set up xyut x-axis (xyut%xx(itime) and xyut%mtimes(:,itime) are ! defined at each call below): ! if (itime==1) then xyut%nx = ntimes if (associated(xyut%xx)) deallocate(xyut%xx) allocate(xyut%xx(xyut%nx),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut_oh allocating xyut%xx") xyut%xlab = 'UT (HRS)' if (associated(xyut%mtimes)) deallocate(xyut%mtimes) allocate(xyut%mtimes(3,ntimes),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut_oh allocating xyut%mtimes") xyut%hvol0 = h%histfile if (diffs) xyut%hvol0_cntr = hcntr%histfile ! ! Save number of selected levels: nslev = 0 levloop: do k=1,mxzpht if (xyut_zpht(k)==spval) cycle levloop if (trim(float_to_str(xyut_zpht(k)))=='hmf2'.or. + trim(float_to_str(xyut_zpht(k)))=='HMF2'.or. | xyut_zpht(k)==hmf2flag) cycle levloop if (xyut_zpht(k) <= gcmlev(npress)) cycle levloop nslev=nslev+1 enddo levloop ! ! Last time: ! Allocate xyut%data(ntimes,ncurve) and xyut%cnames(ncurve), where ! ncurve is number of fields of type (oh-v or oh-b) ! elseif (itime==ntimes) then xyut%ncurve = ntype if (associated(xyut%data)) deallocate(xyut%data) allocate(xyut%data(ntimes,xyut%ncurve),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut_oh allocating xyut%data") ! if (associated(xyut%cnames)) deallocate(xyut%cnames) allocate(xyut%cnames(ntype),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut_oh allocating xyut%cnames") write(6,"('mkxyut_oh allocated xyut%cnames where ntype=', | i3)") ntype loginterp = 0 xyut%hvol1 = h%histfile if (diffs) then xyut%hvol1_cntr = hcntr%histfile difftype = "PERCENT " if (idifpercent/=ispval.and.idifpercent<=0) + difftype = "RAW " endif endif ! itime==1 or itime==ntimes ! ! Define current time in xyut: xyut%xx(itime) = float(h%mtime(1))*24.+h%ut xyut%mtimes(:,itime) = h%mtime ! use pert for now ! ! Location loop: iloc = 0 loc_loop1: do l=1,mxloc if (xyut_locs(1,l)==spval.or.xyut_locs(2,l)==spval) + cycle loc_loop1 iloc = iloc+1 call setloctype(xyut_locs(1,l),xyut_locs(2,l),h%ut, + xyut%glat,xyut%glon,xyut%slt,ixlat,ixlon,xyut%loctype) ! ! Set location name, if given by user: if (itime==ntimes) then xyut%locname = "" if (len_trim(xyut_locname(l))>0) + xyut%locname = ' ('//trim(xyut_locname(l))//')' ! ! Set field names (this is for top title and x-lab -- individual ! curve names are set below (xyut%cnames)). ! (names are defined here instead of above in case doppler ! fields have just been made and names need to be reset for ! emission field at next location) ! if (trim(type)=='OH-VIB') then xyut%sname = 'OH-VIBS ' if (.not.diffs) then xyut%fname = 'OH VIBRATIONAL LEVELS' else xyut%fname = trim(difftype)//' DIFFS: '// + 'OH VIBRATIONAL LEVELS' endif endif if (trim(type)=='OH-BAND') then xyut%sname = 'OH-BANDS' if (.not.diffs) then xyut%fname = 'OH EMISSION BANDS' else xyut%fname = trim(difftype)//' DIFFS: '// + 'OH EMISSION BANDS' endif endif endif ! ! Fields loop: ifld = 0 fields_loop1: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or. + trim(f(ixf)%type)/=type) cycle fields_loop1 ifld = ifld+1 if (itime==ntimes) then xyut%funits = f(ixf)%units ! ! 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(xyut%cnames(ifld),"('V-',i1)") i else read(f(ixf)%fname8(len-1:len),"(i2)") i write(xyut%cnames(ifld),"('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(xyut%cnames(ifld),"(i1,'-',i1)") hi,lo endif endif ! itime==ntimes ! ! Define xyutdat_oh(itime,iloc,:,ifld) and xyutdat_cntr_oh if diffs: if (xyut%loctype=='lon'.or.xyut%loctype=='slt') then xyutdat_oh(itime,iloc,:,ifld) = f(ixf)%data(ixlon,ixlat,:) if (diffs) xyutdat_cntr_oh(itime,iloc,:,ifld) = + fcntr(ixf)%data(ixlon,ixlat,:) elseif (xyut%loctype== 'zm ') then do k=1,f(ixf)%nlev xyutdat_oh(itime,iloc,k,ifld) = + fmean(f(ixf)%data(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,f(ixf)%nlev xyutdat_cntr_oh(itime,iloc,k,ifld) = + fmean(fcntr(ixf)%data(:,ixlat,k),nlon,spval,0) enddo endif elseif (xyut%loctype== 'gm ') then do k=1,f(ixf)%nlev xyutdat_oh(itime,iloc,k,ifld) = + fglbm(f(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (diffs) then do k=1,f(ixf)%nlev xyutdat_cntr_oh(itime,iloc,k,ifld) = + fglbm(fcntr(ixf)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif else write(6,"('>>> WARNING mkxyut_oh: unknown loctype=',a)") + xyut%loctype endif ! ! Define xyutdat_ohbdop for each oh-band emission: if (ixyut_doppler>0.and.trim(type)=='OH-BAND') then call savedop_xyut(f,nf,f(ixf)%fname8,f(ixf)%type, + xyutdat_oh(itime,iloc,:,ifld),f(ixf)%nlev, + ixlat,ixlon,xyutdat_ohbdop(itime,iloc,ifld,:)) if (diffs) then call savedop_xyut(fcntr,nf,fcntr(ixf)%fname8, + fcntr(ixf)%type,xyutdat_cntr_oh(itime,iloc,:,ifld), + fcntr(ixf)%nlev,ixlat,ixlon, + xyutdat_cntr_ohbdop(itime,iloc,ifld,:)) endif endif ! ! End fields loop: enddo fields_loop1 if (itime < ntimes) cycle loc_loop1 ! ! Lev loop (enter only when itime==ntimes) ! (xyutdat_oh is now fully defined) ! (Process oh only for requested height surfaces) ! xyut%log10 = 0 ilev = 0 lev_loop1: do k=1,mxzpht ! ! Skip if not a selected height, unless there are no selected ! heights (nslev==0) and we are doing OH-B (in which case ! we need to do ht-integration of OH-B): ! (combining these 3 conditionals results in FPE when 'hmf2' ! is in xyut_zpht(k)) ! if (xyut_zpht(k)==spval.and. + .not.(nslev==0.and.type=='OH-BAND')) cycle lev_loop1 ! ! Oh fields not calculated at hmf2 heights: if (trim(float_to_str(xyut_zpht(k)))=='hmf2'.or. + trim(float_to_str(xyut_zpht(k)))=='HMF2'.or. | xyut_zpht(k)==hmf2flag) cycle lev_loop1 ! ! Oh fields calculated at heights only: if (xyut_zpht(k) <= gcmlev(npress)) cycle lev_loop1 ! ! Set levtype: if (xyut_zpht(k)==spval) then xyut%levtype = 'integ ' elseif (xyut_zpht(k) > gcmlev(npress)) then xyut%levtype = 'ht ' else write(6,"('>>> mkxyut_oh: cannot set xyut%levtype:', + ' k=',i3,' xyut_zpht(k)=',e12.4,' type= ',a)") + k,xyut_zpht(k),type cycle lev_loop1 endif ilev = ilev+1 ! ! Jump to here from below to do integrations: 100 continue if (diffs) then allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut_oh allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut_oh allocating futdat2") endif ! ! Interpolate to requested height (must be within oh range): if (xyut%levtype=='ht ') then if (xyut_zpht(k) < oh_alt(1).or. + xyut_zpht(k) > oh_alt(nohalt)) then write(6,"('>>> mkxyut_oh: selected height surface ', + f10.2,' is outside range of oh: oh bottom,top=', + 2f10.2)") xyut_zpht(k),oh_alt(1),oh_alt(nohalt) cycle lev_loop1 endif xyut%ht = xyut_zpht(k) do ifld = 1,ntype if (.not.diffs) then do it=1,ntimes call intloc(xyutdat_oh(it,iloc,:,ifld),oh_alt, + nohalt,xyut_zpht(k),1,loginterp,xyut%data(it,ifld), + 1,1,1,ier,spval,0) enddo else do it=1,ntimes call intloc(xyutdat_oh(it,iloc,:,ifld),oh_alt, + nohalt,xyut_zpht(k),1,loginterp,futdat1(it), + 1,1,1,ier,spval,0) call intloc(xyutdat_cntr_oh(it,iloc,:,ifld),oh_alt, + nohalt,xyut_zpht(k),1,loginterp,futdat2(it), + 1,1,1,ier,spval,0) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,ifld),ntimes, + difftype) endif enddo ! ifld=1,ntype ! ! Integrate in height (OH-BAND(s) only): else ! 'integ' do ifld=1,xyut%ncurve if (.not.diffs) then do it=1,ntimes xyut%data(it,ifld) = + quadrat(xyutdat_oh(it,iloc,:,ifld),oh_alt,nohalt) enddo else do it=1,ntimes futdat1(it) = quadrat(xyutdat_oh(it,iloc,:,ifld), + oh_alt,nohalt) futdat2(it) = quadrat(xyutdat_cntr_oh(it,iloc,:,ifld), + oh_alt,nohalt) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,ifld),ntimes, + difftype) endif enddo if (ibohv_watts <= 0) then xyut%data = xyut%data * 1.e-9 ! kilo-rayleighs xyut%funits = "KILO-RAYLEIGHS" else xyut%data = xyut%data / (4.*pi) ! watts/cm2-str xyut%funits = "WATTS/CM2-STR" endif endif ! ht or ht-integ if (diffs) then deallocate(futdat1) deallocate(futdat2) endif ! ! doxyut makes the plot and writes output data: call doxyut(xyut,h,hcntr,itime) ! ! Do doppler fields of integrated OH-B. ! Plot doppler t,u,v separately, each w/ multiple curves of oh-bands. ! (xyutdat_ohbdop was saved by savedop_xyut) ! if (xyut%levtype=='integ') then ! just plotted integ if (ixyut_doppler > 0) then if (diffs) then allocate(futdat1(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut_oh allocating futdat1") allocate(futdat2(ntimes),stat=ier) if (ier /= 0) + call allocerr(ier,"mkxyut_oh allocating futdat2") difftype = "RAW " if (idifpercent > 0) difftype = "PERCENT " endif do i=1,nfdop do ifld=1,xyut%ncurve do it=1,ntimes xyut%data(it,ifld) = xyutdat_ohbdop(it,iloc,ifld,i) enddo ! ! Diffs of integrated oh-b doppler fields: if (diffs) then futdat1 = xyut%data(:,ifld) do it=1,ntimes futdat2(it) = xyutdat_cntr_ohbdop(it,iloc,ifld,i) enddo call mkdiffs(futdat1,futdat2,xyut%data(:,ifld), + ntimes,difftype) write(xyut%fname,"(a,' DIFFS: OH EMISSION BANDS ', + a)") trim(difftype),trim(doppler_fnames(i)) else write(xyut%fname,"('OH EMISSION BANDS ',a)") + trim(doppler_fnames(i)) endif xyut%sname = doppler_snames(i) xyut%funits = doppler_units(i) enddo call doxyut(xyut,h,hcntr,itime) ! ! Doppler u+v 1d vector plot: if (i==2) then ! load u allocate(udop(size(xyut%data,1)),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut_oh allocating udop") udop = xyut%data(:,1) elseif (i==3) then ! load v and make plot allocate(vdop(size(xyut%data,1)),stat=ier) if (ier/=0) + call allocerr(ier,"mkxyut_oh allocating vdop") vdop = xyut%data(:,1) xyut%data(:,1:1) = vecsum(udop,vdop,xyut%nx,1,spval) xyut%funits = doppler_units(i) ! xyut%sname wrong (need f(ixf)%fname8) if (.not.diffs) then write(xyut%fname,"('OH EMISSION: DOPPLER U+V')") else write(xyut%fname,"(a,' DIFFS: OH EMISSION: ', + 'DOPPLER U+V')") trim(difftype) endif xyut%sname = 'DOP(U+V)' call doxyut(xyut,h,hcntr,itime,udop,vdop) deallocate(udop) deallocate(vdop) endif enddo ! i=1,nfdop if (diffs) then deallocate(futdat1) deallocate(futdat2) difftype = "PERCENT " if (idifpercent/=ispval.and.idifpercent<=0) + difftype = "RAW " endif endif ! ixyut_doppler exit lev_loop1 endif ! integ ! ! If OH-B was just plotted and this is last selected level, ! go back up and do ht-integration: if (integht > 0.and.trim(type)=='OH-BAND'.and.ilev==nslev)then xyut%levtype = 'integ' goto 100 endif enddo lev_loop1 enddo loc_loop1 end subroutine mkxyut_oh end subroutine mkxyut ! remaining routines are external to mkxyut !------------------------------------------------------------------- subroutine savedop_xyut(f,nf,fname,ftype,emis,kmx,ixlat,ixlon, + fdop) use proc,only: spval,dlat,dlon,nlat,nlon,gcmlat use fields,only: field,nfdop,oh_alt,nohalt ! ! Given emission emis(kmx) (OH-BAND or EMISSION field), return doppler ! t,u,v in fdop(nfdop). Emission has been calculated at lat,lon ! indices ixlat,ixlon. ! For OH-BAND emissions, t,u,v must be interpolated to OH heights ! (oh_alt(nohalt)) (use f(ix[tuv])%data(ixlon,ixlat,:)) ! (fname and ftype are short name and type of emission field) ! ! Args: integer,intent(in) :: nf type(field),intent(in) :: f(nf) character(len=*),intent(in) :: fname,ftype integer,intent(in) :: kmx,ixlat,ixlon real,intent(in) :: emis(kmx) real,intent(out) :: fdop(nfdop) ! ! Locals: real :: t(kmx),u(kmx),v(kmx),z(kmx),temp(kmx),sigma, + dopt,dopu,dopv,dattpv,attpz,dattpz integer :: ixt,ixu,ixv,ixz,k,ier,kbot ! ! Externals: integer,external :: ixfindc real,external :: fmean,fglbm ! ixt = ixfindc(f%fname8,nf,'TN ') ixu = ixfindc(f%fname8,nf,'UN ') ixv = ixfindc(f%fname8,nf,'VN ') ixz = ixfindc(f%fname8,nf,'Z ') if (ixt==0.or.ixu==0.or.ixv==0.or.ixz==0) then write(6,"('>>> savedop_xyut: need t,u,v,z to make doppler', + ' fields: ixt=',i2,' ixu=',i2,' ixv=',i2,' ixz=',i2)") + ixt,ixu,ixv,ixz endif ! ! Save at regular longitude, latitude: if (xyut%loctype=='lon'.or.xyut%loctype=='slt') then if (trim(ftype)=='EMISSION') then t(:) = f(ixt)%data(ixlon,ixlat,:) u(:) = f(ixu)%data(ixlon,ixlat,:) v(:) = f(ixv)%data(ixlon,ixlat,:) z(:) = f(ixz)%data(ixlon,ixlat,:) else call intloc(f(ixt)%data(ixlon,ixlat,:), + f(ixz)%data(ixlon,ixlat,:),size(f(ixt)%data,3),oh_alt, + nohalt,0,t,1,1,nohalt,ier,spval,0) call intloc(f(ixu)%data(ixlon,ixlat,:), + f(ixz)%data(ixlon,ixlat,:),size(f(ixu)%data,3),oh_alt, + nohalt,0,u,1,1,nohalt,ier,spval,0) call intloc(f(ixv)%data(ixlon,ixlat,:), + f(ixz)%data(ixlon,ixlat,:),size(f(ixv)%data,3),oh_alt, + nohalt,0,v,1,1,nohalt,ier,spval,0) z = oh_alt endif ! ! Save zonal means: elseif (xyut%loctype== 'zm ') then do k=1,kmx t(k) = fmean(f(ixt)%data(:,ixlat,k),nlon,spval,0) u(k) = fmean(f(ixu)%data(:,ixlat,k),nlon,spval,0) v(k) = fmean(f(ixv)%data(:,ixlat,k),nlon,spval,0) z(k) = fmean(f(ixz)%data(:,ixlat,k),nlon,spval,0) enddo if (trim(ftype)=='OH-BAND') then call intloc(t,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) t = temp call intloc(u,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) u = temp call intloc(v,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) v = temp call intloc(z,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) z = temp endif ! ! Save global means: elseif (xyut%loctype=='gm ') then do k=1,kmx t(k) = fglbm(f(ixt)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) u(k) = fglbm(f(ixu)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) v(k) = fglbm(f(ixv)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) z(k) = fglbm(f(ixz)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (trim(ftype)=='OH-BAND') then call intloc(t,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) t = temp call intloc(u,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) u = temp call intloc(v,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) v = temp call intloc(z,z,kmx,oh_alt,nohalt,0,temp,1,1,nohalt,ier, + spval,0) z = temp endif else write(6,"('>>> WARNING savedop_xyut: unknown loctype=',a)") + xyut%loctype endif ! ! If OH-BAND, then bottom of oh_alt (30km) is probably below bottom ! of tgcm, therefore intloc returned spval in t,u,v where bottom of ! tgcm z > 30km. Doppler does not handle spval, so call doppler at ! heights starting at bottom of tgcm (i.e., where spvals stop). ! If all spvals are not at bottom, print warning. ! Same for emis (because getoh returned spval for emis below model) ! if (trim(ftype)=='OH-BAND') then kbot = 0 tloop: do k=1,kmx if (t(k)/=spval) then kbot = k if (any(t(kbot:kmx)==spval)) + write(6,"('>>> savedop_xyut for OH-BAND: found ', + ' spval above non-spval t: t=',/(6e12.4))") t exit tloop endif enddo tloop uloop: do k=1,kmx if (u(k)/=spval) then kbot = max(k,kbot) if (any(u(kbot:kmx)==spval)) + write(6,"('>>> savedop_xyut for OH-BAND: found ', + ' spval above non-spval u: u=',/(6e12.4))") u exit uloop endif enddo uloop vloop: do k=1,kmx if (v(k)/=spval) then kbot = max(k,kbot) if (any(v(kbot:kmx)==spval)) + write(6,"('>>> savedop_xyut for OH-BAND: found ', + ' spval above non-spval v: v=',/(6e12.4))") v exit vloop endif enddo vloop eloop: do k=1,kmx if (emis(k)/=spval) then kbot = max(k,kbot) if (any(emis(kbot:kmx)==spval)) + write(6,"('>>> savedop_xyut for OH-BAND: found ', + ' spval above non-spval emis: emis=',/(6e12.4))") v exit eloop endif enddo eloop if (kbot==0) then write(6,"('>>> savedop_xyut: kbot=0 (t, u, or v all ', + 'spval')") endif else ! non-ohb emission (t,u,v should be ok) kbot = 1 endif ! ! Emission fields are: E6300 E5577 EO200 EOH83 ECO215u ENO53u ! (doppler fields not calculated for ECO215u and EN053u) ! if (trim(fname)=='E6300') then sigma = 15867.852 elseif (trim(fname)=='E5577') then sigma = 17924. elseif (trim(fname)=='EO200') then sigma = 13125. elseif (trim(fname)=='EOH83') then sigma = 3571. elseif (trim(ftype)=='OH-BAND') then sigma = 11900. else write(6,"('>>> savedop_xyut: do not know how to set sigma:', + ' fname=',a,' ftype=',a)") fname,ftype endif ! ! SUBROUTINE doppler(NP,HM,TM,ZM,VM,VER,sigma, ! + VELV,VELZ,ATTPV,DATTPV,ATTPZ,DATTPZ) ! call doppler(kmx-kbot+1,z(kbot:kmx),t(kbot:kmx),u(kbot:kmx), + v(kbot:kmx),emis(kbot:kmx),sigma,dopv,dopu,dopt,dattpv,attpz, + dattpz,spval) fdop(1) = dopt fdop(2) = dopu fdop(3) = dopv ! write(6,"('savedop returning: fdop=',3e12.4)") fdop end subroutine savedop_xyut !------------------------------------------------------------------- subroutine alloc_xyutdat(f,nf,ixz,type,xyuts,ier) ! ! Allocate xyutdat(ntimes,nloc,nlev,nf). ! ntimes number of times ! nloc from user-provided xyut_locs ! nlev npress if type=='noheight', otherwise is nohalt ! nf number of fields ! ixz is returned as index of heights field in xyutdat ! use proc use fields,only: field,nohalt,nfdop use input,only: xyut_locs,ixyut_doppler ! ! Args: integer,intent(in) :: nf integer,intent(out) :: ier,ixz type(field),intent(in) :: f(nf) character(len=16),intent(in) :: type type(xyut_type),intent(inout) :: xyuts ! ! Locals: integer :: nloc,nfld,nlev,i,istat,nemis ! ier = 0 ! ! Count number of locations: nloc = 0 ! number of valid locations do i=1,mxloc if (xyut_locs(1,i)/=spval.and.xyut_locs(2,i)/=spval) + nloc = nloc+1 enddo if (nloc==0) then write(6,"('>>> alloc_xyutdat: need locations xyut_locs ', + '(nloc==0).')") ier = 1 return endif xyuts%nsloc = nloc ! ! Count number of fields: ! nfld = 0 ! number of fields to allocate nemis = 0 ! number of emission fields ixz = 0 ! index to heights (output) if (trim(type)=='noheight') then do i=1,nf if (f(i)%requested.and.associated(f(i)%data).and. + trim(f(i)%vtype)/='HEIGHT') then nfld = nfld+1 if (trim(f(i)%type)=='EMISSION') nemis = nemis+1 ! ! 5/3 bugfix: ixz=nfld, NOT ixz=i ! if (trim(f(i)%fname8)=='Z') ixz = i if (trim(f(i)%fname8)=='Z') ixz = nfld endif enddo else do i=1,nf if (f(i)%requested.and.associated(f(i)%data).and. + f(i)%type==type) nfld = nfld+1 enddo endif if (nfld==0) then write(6,"('Alloc_xyutdat: no fields to allocate: type=',a)") + trim(type) ier = 1 return endif ! ! Add field for total heights: if (trim(type)=='noheight') then nlev = npress if (ixz==0) then nfld = nfld+1 ixz = nfld endif else nlev = nohalt endif ! ! Allocate xyutdat: if (trim(type)=='noheight') then allocate(xyutdat(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier,"alloc_xyutdat allocating xyutdat") else write(6,"('Allocated xyutdat(ntimes,nloc,nlev,nfld),', + ' where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif if (diffs) then allocate(xyutdat_cntr(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat_cntr(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier,"alloc_xyutdat allocating xyutdat_cntr") else write(6,"('Allocated xyutdat_cntr(ntimes,nloc,nlev,nfld),', + ' where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif endif ! ! Allocate xyutdat_emisdop: if (ixyut_doppler > 0.and.nemis > 0) then allocate(xyutdat_emisdop(ntimes,nloc,nemis,nfdop),stat=istat) if (istat /= 0) then call allocerr(ier, + "alloc_xyutdat allocating xyutdat_emisdop") else write(6,"('Allocated xyutdat_emisdop', + '(ntimes,nloc,nemis,nfdop) where:',/,' ntimes=',i4, + ' nloc=',i3,' nemis=',i3,' nfdop=',i3)") + ntimes,nloc,nemis,nfdop endif if (diffs) then allocate(xyutdat_cntr_emisdop(ntimes,nloc,nemis,nfdop), + stat=istat) if (istat /= 0) then call allocerr(ier, + "alloc_xyutdat allocating xyutdat_cntr_emisdop") else write(6,"('Allocated xyutdat_cntr_emisdop', + '(ntimes,nloc,nemis,nfdop) where:',/,' ntimes=',i4, + ' nloc=',i3,' nemis=',i3,' nfdop=',i3)") + ntimes,nloc,nemis,nfdop endif endif endif ! ! Allocate xyutdat_ohv: elseif (trim(type)=='OH-VIB') then allocate(xyutdat_ohv(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat_ohv(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier,"alloc_xyutdat allocating xyutdat_ohv") else write(6,"('Allocated xyutdat_ohv(ntimes,nloc,nlev,nfld),', + ' where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif if (diffs) then allocate(xyutdat_cntr_ohv(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat_cntr_ohv(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier, + "alloc_xyutdat allocating xyutdat_cntr_ohv") else write(6,"('Allocated xyutdat_cntr_ohv(ntimes,nloc,nlev,', + 'nfld) where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif endif ! ! Allocate xyutdat_ohb: elseif (trim(type)=='OH-BAND') then allocate(xyutdat_ohb(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat_ohb(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier,"alloc_xyutdat allocating xyutdat_ohb") else write(6,"('Allocated xyutdat_ohb(ntimes,nloc,nlev,nfld),', + ' where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif if (diffs) then allocate(xyutdat_cntr_ohb(ntimes,nloc,nlev,nfld),stat=istat) if (istat /= 0) then write(6,"('>>> Error allocating', + ' xyutdat_cntr_ohb(ntimes,nloc,nlev,nfld): ntimes=', + i4,' nloc=',i3,' nlev=',i3,' nfld=',i3,' type=',a)") + ntimes,nloc,nlev,nfld,trim(type) call allocerr(ier, + "alloc_xyutdat allocating xyutdat_cntr_ohb") else write(6,"('Allocated xyutdat_cntr_ohb(ntimes,nloc,nlev,', + 'nfld) where:',/,' ntimes=',i4,' nloc=',i2,' nlev=',i2, + ' nfld=',i2,' type=',a)") ntimes,nloc,nlev,nfld,trim(type) endif endif ! ! Allocate xyutdat_ohbdop: if (ixyut_doppler > 0) then allocate(xyutdat_ohbdop(ntimes,nloc,nfld,nfdop),stat=istat) if (istat /= 0) then call allocerr(ier, + "alloc_xyutdat allocating xyutdat_ohbdop") else write(6,"('Allocated xyutdat_ohbdop', + '(ntimes,nloc,nfld,nfdop) where:',/,' ntimes=',i4, + ' nloc=',i3,' nfld=',i3,' nfdop=',i3)") + ntimes,nloc,nfld,nfdop endif if (diffs) then allocate(xyutdat_cntr_ohbdop(ntimes,nloc,nfld,nfdop), + stat=istat) if (istat /= 0) then call allocerr(ier, + "alloc_xyutdat allocating xyutdat_cntr_ohbdop") else write(6,"('Allocated xyutdat_cntr_ohbdop', + '(ntimes,nloc,nfld,nfdop) where:',/,' ntimes=',i4, + ' nloc=',i3,' nfld=',i3,' nfdop=',i3)") + ntimes,nloc,nfld,nfdop endif endif endif endif end subroutine alloc_xyutdat !------------------------------------------------------------------- subroutine mkxyutlabs(xyut,h,hcntr,msgout) use proc,only: diffs,spval use input,only: ie5577,ie6300 use plt,only: scfac ! ! Construct 3 top labels (xyut%tlabs) and 3 bottom labels ! (ut%blabs) for xy vertical profiles: ! ! Args: type(xyut_type),intent(inout) :: xyut type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Locals: integer :: lenlab,ier real :: fmin,fmax real,allocatable :: dat(:,:) ! ! Set up top and bottom text labels: ! ! Top (1st) top label is optionally provided by user: xyut%tlabs(1) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if xyut%log10 > 0) ! if (xyut%known) then lenlab = len_trim(xyut%fname)+len_trim(xyut%funits)+3 if (xyut%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(xyut%tlabs(2))) then if (xyut%log10 > 0) then xyut%tlabs(2) = 'LOG10 '//trim(xyut%fname)//' ('// + trim(xyut%funits)//')' else xyut%tlabs(2) = trim(xyut%fname)//' ('// + trim(xyut%funits)//')' endif else xyut%tlabs(2) = trim(xyut%fname) endif else ! unknown to proc lenlab = len_trim(xyut%sname)+2+len_trim(xyut%fname)+ | len_trim(xyut%funits)+3 if (xyut%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(xyut%tlabs(2))) then if (xyut%log10 > 0) then if (trim(xyut%sname)/=trim(xyut%fname)) then xyut%tlabs(2) = 'LOG10 '//trim(xyut%sname)//': '// | trim(xyut%fname)//' ('//trim(xyut%funits)//')' else xyut%tlabs(2) = 'LOG10 '//trim(xyut%fname)// | ' ('//trim(xyut%funits)//')' endif else ! no log10 if (trim(xyut%sname)/=trim(xyut%fname)) then xyut%tlabs(2) = trim(xyut%sname)//': '// | trim(xyut%fname)//' ('//trim(xyut%funits)//')' else xyut%tlabs(2) = trim(xyut%fname)// | ' ('//trim(xyut%funits)//')' endif endif else xyut%tlabs(2) = trim(xyut%fname) endif endif ! known or unknown ! ! Add DIFFS to label: if (diffs) then if (trim(xyut%difftype)=='RAW') then if (len_trim(xyut%tlabs(2))+7 <= len(xyut%tlabs(2))) then xyut%tlabs(2) = 'DIFFS: '//trim(xyut%tlabs(2)) else xyut%tlabs(2) = 'DIFFS: '//trim(xyut%fname) endif else ! percent if (len_trim(xyut%tlabs(2))+15 <= len(xyut%tlabs(2))) then xyut%tlabs(2) = trim(xyut%difftype)//' DIFFS: '// | trim(xyut%tlabs(2)) else xyut%tlabs(2) = '% DIFFS: '//trim(xyut%fname) endif endif endif ! Bottom (3rd) top label is location and level: ! xyut location may be lat,lon, lat,slt, zonal means, or global means: select case (xyut%loctype) case ('lon') select case (xyut%levtype) case ('zp ') write(xyut%tlabs(3),"('LAT=',f5.1,' LON=',f6.1, + a,' ZP=',f6.2)") + xyut%glat,xyut%glon,trim(xyut%locname),xyut%zp case ('ht ') write(xyut%tlabs(3),"('LAT=',f5.1,' LON=',f6.1, + a,' HT=',f5.1,' (KM)')") + xyut%glat,xyut%glon,trim(xyut%locname),xyut%ht case ('indep') write(xyut%tlabs(3),"('LAT=',f5.1,' LON=',f6.1,a)") + xyut%glat,xyut%glon,trim(xyut%locname) case ('integ') write(xyut%tlabs(3),"('LAT=',f5.1,' LON=',f6.1,a, + ' (HT-INTEG)')") xyut%glat,xyut%glon, + trim(xyut%locname) case ('hmf2 ') write(xyut%tlabs(3),"('LAT=',f5.1,' LON=',f6.1,a, + ' (HMF2 HEIGHTS)')") + xyut%glat,xyut%glon,trim(xyut%locname) end select case ('slt') select case (xyut%levtype) case ('zp ') write(xyut%tlabs(3),"('LAT=',f5.1,' SLT=',f5.1,' HRS', + a,' ZP=',f6.2)") + xyut%glat,xyut%slt,trim(xyut%locname),xyut%zp case ('ht ') write(xyut%tlabs(3),"('LAT=',f5.1,' SLT=',f5.1,' HRS', + a,' HT=',f5.1,' (KM)')") + xyut%glat,xyut%slt,trim(xyut%locname),xyut%ht case ('indep') write(xyut%tlabs(3),"('LAT=',f5.1,' SLT=',f5.1,' HRS', + a)") xyut%glat,xyut%slt,trim(xyut%locname) case ('integ') write(xyut%tlabs(3),"('LAT=',f5.1,' SLT=',f5.1,' HRS', + a,' (HT-INTEG)')") xyut%glat,xyut%slt,trim(xyut%locname) case ('hmf2 ') write(xyut%tlabs(3),"('LAT=',f5.1,' SLT=',f5.1,' HRS', + a,' (HMF2 HEIGHTS)')") xyut%glat,xyut%slt, + trim(xyut%locname) end select case ('zm ') select case (xyut%levtype) case ('zp ') write(xyut%tlabs(3),"('LAT=',f5.1,' (ZONAL MEANS)', + a,' ZP=',f6.2)") xyut%glat,trim(xyut%locname),xyut%zp case ('ht ') write(xyut%tlabs(3),"('LAT=',f5.1,' (ZONAL MEANS)', + a,' HT=',f6.1)") xyut%glat,trim(xyut%locname),xyut%ht case ('indep') write(xyut%tlabs(3),"('LAT=',f5.1,' (ZONAL MEANS)',a)") + xyut%glat,trim(xyut%locname) case ('integ') write(xyut%tlabs(3),"('LAT=',f5.1,' (ZONAL MEANS)',a, + ' (HT-INTEG)')") xyut%glat,trim(xyut%locname) case ('hmf2 ') write(xyut%tlabs(3),"('LAT=',f5.1,' (ZONAL MEANS)',a, + ' (HMF2 HEIGHTS)')") xyut%glat,trim(xyut%locname) end select case ('gm ') select case (xyut%levtype) case ('zp ') write(xyut%tlabs(3),"('GLOBAL MEANS: ZP=',f6.2)") + xyut%zp case ('ht ') write(xyut%tlabs(3),"('GLOBAL MEANS: HT=',f5.1,' (KM)')") + xyut%ht case ('indep') write(xyut%tlabs(3),"('GLOBAL MEANS')") case ('integ') write(xyut%tlabs(3),"('GLOBAL MEANS (HT-INTEG)')") case ('hmf2 ') write(xyut%tlabs(3),"('GLOBAL MEANS (HMF2 HEIGHTS)')") end select end select ! ! Top (1st) bottom label is min,max,interval: allocate(dat(xyut%nx,xyut%ncurve),stat=ier) if (ier/=0) call allocerr(ier,'mkxyutlabs allocating dat') dat = xyut%data if (xyut%log10>0) call log10f(dat,xyut%nx*xyut%ncurve, + 1.e-20,spval) call fminmax(dat,xyut%nx*xyut%ncurve,fmin,fmax) if (scfac == 1.) then write(xyut%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(xyut%blabs(1),"('MIN,MAX=',2(1pe12.4),' (X',1pe9.2, | ')')") fmin,fmax,scfac endif deallocate(dat) ! ! Middle (2nd) and bottom (3rd) bottom label are history info: ! if (.not.diffs) then ! if (len_trim(xyut%hvol0)+len_trim(xyut%hvol1)+15 <= ! + len(xyut%blabs(2))) then ! write(xyut%blabs(2),"('FIRST, LAST: ',a,', ',a)") ! + trim(xyut%hvol0), trim(xyut%hvol1) ! else ! write(xyut%blabs(2),"('FIRST: ',a)") trim(xyut%hvol0) ! endif ! xyut%blabs(3) = ' ' ! if (trim(xyut%ftype)=='EMISSION') then ! call mkemislab(xyut%sname,ie5577,ie6300,xyut%blabs(3)) ! if (len_trim(xyut%blabs(3)) > 48) blabchsz(3) = .015 ! endif ! else ! if (len_trim(xyut%hvol0)+len_trim(xyut%hvol1)+22 <= ! + len(xyut%blabs(2)).and. ! + len_trim(xyut%hvol0_cntr)+len_trim(xyut%hvol1_cntr)+22 <= ! + len(xyut%blabs(3))) then ! write(xyut%blabs(2),"('FIRST, LAST: ',a,', ',a,' MINUS ')") ! + trim(xyut%hvol0), trim(xyut%hvol1) ! write(xyut%blabs(3),"('FIRST, LAST: ',a,', ',a)") ! + trim(xyut%hvol0_cntr), trim(xyut%hvol1_cntr) ! else ! write(xyut%blabs(2),"('FIRST: ',a,' MINUS ')") ! + trim(xyut%hvol0) ! write(xyut%blabs(3),"('FIRST: ',a)") trim(xyut%hvol0_cntr) ! endif ! endif ! ! 5/05 btf: Put first and last history files in blabs(2 and 3): if (.not.diffs) then if (trim(xyut%hvol0) /= trim(xyut%hvol1)) then write(xyut%blabs(2),"('FIRST: ',a)") trim(xyut%hvol0) write(xyut%blabs(3),"('LAST: ',a)") trim(xyut%hvol1) else write(xyut%blabs(2),"('MODEL ',a)") trim(h%version) write(xyut%blabs(3),"(a)") trim(xyut%hvol0) endif else ! diffs: show perturbed and control files: if (trim(xyut%hvol0) /= trim(xyut%hvol1)) then write(xyut%blabs(2),"('PERT FIRST: ',a,' LAST: ',a)") | trim(xyut%hvol0),trim(xyut%hvol1) write(xyut%blabs(3),"('CNTR FIRST: ',a,' LAST: ',a)") | trim(xyut%hvol0_cntr),trim(xyut%hvol1_cntr) else write(xyut%blabs(2),"('PERT FILE: ',a)") trim(xyut%hvol0) write(xyut%blabs(3),"('CNTR FILE: ',a)") trim(xyut%hvol0_cntr) endif endif ! ! Return message to be printed to stdout: ! ! At lon,lat: select case (xyut%loctype) case ('lon') select case (xyut%levtype) case ('zp ') write(msgout,"(a,' lat=',f5.1,' lon=',f6.1,' zp=',f5.1)") + xyut%sname,xyut%glat,xyut%glon,xyut%zp case ('ht ') write(msgout,"(a,' lat=',f5.1,' lon=',f6.1,' ht=',f6.1)") + xyut%sname,xyut%glat,xyut%glon,xyut%ht case ('indep') write(msgout,"(a,' lat=',f5.1,' lon=',f6.1,' ht-indep')") + xyut%sname,xyut%glat,xyut%glon case ('integ') write(msgout,"(a,' lat=',f5.1,' lon=',f6.1,' ht-integ')") + xyut%sname,xyut%glat,xyut%glon case ('hmf2 ') write(msgout,"(a,' lat=',f5.1,' lon=',f6.1,' hmf2 hts')") + xyut%sname,xyut%glat,xyut%glon end select ! ! At slt,lat: case ('slt') select case (xyut%levtype) case ('zp ') write(msgout,"(a,' lat=',f5.1,' slt=',f4.1,' zp=',f5.1)") + xyut%sname,xyut%glat,xyut%slt,xyut%zp case ('ht ') write(msgout,"(a,' lat=',f5.1,' slt=',f4.1,' ht=',f6.1)") + xyut%sname,xyut%glat,xyut%slt,xyut%ht case ('indep') write(msgout,"(a,' lat=',f5.1,' slt=',f4.1,' ht-indep')") + xyut%sname,xyut%glat,xyut%slt case ('integ') write(msgout,"(a,' lat=',f5.1,' slt=',f4.1,' ht-integ')") + xyut%sname,xyut%glat,xyut%slt case ('hmf2 ') write(msgout,"(a,' lat=',f5.1,' slt=',f4.1,' hmf2 hts')") + xyut%sname,xyut%glat,xyut%slt end select ! ! At zonal means: case ('zm ') select case (xyut%levtype) case ('zp ') write(msgout,"(a,' lat=',f5.1,' zonal means zp=',f5.1)") + xyut%sname,xyut%glat,xyut%zp case ('ht ') write(msgout,"(a,' lat=',f5.1,' zonal means ht=',f6.1)") + xyut%sname,xyut%glat,xyut%ht case ('indep') write(msgout,"(a,' lat=',f5.1,' zonal means ht-indep')") + xyut%sname,xyut%glat case ('integ') write(msgout,"(a,' lat=',f5.1,' zonal means ht-integ')") + xyut%sname,xyut%glat case ('hmf2 ') write(msgout,"(a,' lat=',f5.1,' zonal means hmf2 hts')") + xyut%sname,xyut%glat end select ! ! At global means: case ('gm ') select case (xyut%levtype) case ('zp ') write(msgout,"(a,' global means: zp=',f5.1)") + xyut%sname,xyut%zp case ('ht ') write(msgout,"(a,' global means: ht=',f6.1)") + xyut%sname,xyut%ht case ('indep') write(msgout,"(a,' global means: ht-indep')") xyut%sname case ('integ') write(msgout,"(a,' global means: ht-integ')") xyut%sname case ('hmf2 ') write(msgout,"(a,' global means: hmf2 hts')") xyut%sname end select end select ! ! Add min,max: write(msgout(43:72),"('min,max=',2e11.4)") fmin,fmax end subroutine mkxyutlabs !------------------------------------------------------------------- subroutine pltxyut(xyut,vp) ! ! Make line plot from xyut structure, ut on x-axis, field on y-axis: ! use plt,only: igks_cgm,igks_ps,igks_x11, + iwk_cgm, iwk_ps, iwk_x11, labutxy,scfac use proc,only: spval use input,only: ilab_hq ! ! Args: type (xyut_type),intent(inout) :: xyut real,intent(in) :: vp(4) ! ! Locals: integer :: iplt_cgm,iplt_ps,iplt_x11,ier real :: toffset=.04, boffset integer :: isltax,multiy character(len=32) :: ylab ! ylab = trim(xyut%sname)//' ('//trim(xyut%funits)//')' if (xyut%log10 == 1) then ylab = 'LOG10 '//trim(xyut%sname)//' ('// + trim(xyut%funits)//')' else ylab = trim(xyut%sname)//' ('// + trim(xyut%funits)//')' endif isltax = 0 ; boffset=.2 if (xyut%loctype=='lon') then isltax = 1 ; boffset = .35 endif multiy = 0 if (associated(xyut%cnames)) multiy = 1 ! ! subroutine pltxmy(x,y,npts,nycurve,multiy,xlab,ylab, ! + linelab,chsize,vp,ilog,scfac,spval) ! if (.not.associated(xyut%cnames)) then allocate(xyut%cnames(1),stat=ier) xyut%cnames(1) = ' ' endif if (xyut%ncurve<=1) multiy=0 call pltxmy(xyut%xx,xyut%data,xyut%nx,xyut%ncurve,multiy, + ' ',trim(ylab),xyut%cnames,0.,vp,xyut%log10,scfac,spval) call labutxy(xyut%mtimes,xyut%nx,xyut%data(:,1),0,'DUM',0., + isltax,xyut%glon) call wrlab6(xyut%tlabs,toffset,tlabchsz, + xyut%blabs,boffset,blabchsz,vp,ilab_hq) end subroutine pltxyut !------------------------------------------------------------------- subroutine wrout_xyut(xyut,h,hcntr,itime) use plt,only: zmin,zmax use proc use input ! ! Write ascii and/or xdr output data files for current xyut frame: ! ! Args: ! (xyut must be intent(inout) because this routine may need ! to call mkxyutlabs) ! type(xyut_type),intent(inout) :: xyut type(history),intent(in) :: h,hcntr integer,intent(in) :: itime ! ! Locals: character(len=80) :: msgout real :: dum ! if (iwrdat==0.and.iwrxdr==0.and.len_trim(sendcdf_xyut)==0.and. | len_trim(flnm_cdf_xyut)==0) return if (iplot == 0) call mkxyutlabs(xyut,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,xyut%data,xyut%nx,1, + xyut%xx,dum,xyut%xlab,xyut%sname,xyut%tlabs(2), + xyut%tlabs(3),xyut%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_xyut: xdr files NOT available for ', + 'xy vs ut')") endif ! ! Write netcdf file: if (len_trim(sendcdf_xyut) > 0.or.len_trim(flnm_cdf_xyut) > 0) | call wrcdf_xyut(xyut,h,hcntr,itime,.false.) return end subroutine wrout_xyut !------------------------------------------------------------------- subroutine wrcdf_xyut(xyut,h,hcntr,itime,histonly) ! ! Write netcdf data file. Current data and info are in xyut structure. ! This is called at every model time to update history variables, but ! data is written only after the last time. If histonly is true, only ! update history vars, otherwise (histonly=false) field data can is ! written. ! use proc,only: zmflag,gmflag,spval use mkcdf_module,only: | netcdf_file, ! type def | nccreate,ncdefine,ncvalues ! subs include 'netcdf.inc' ! ! Args: type(xyut_type),intent(inout) :: xyut type(history),intent(in) :: h,hcntr integer,intent(in) :: itime logical,intent(in) :: histonly ! ! Saved locals: integer,save :: ncid,ncalls=0,ntimes,nflds=1 integer,save :: id_nslev,id_nsloc,id_ncoords,id_char8, | id_char56,idv_slev,idv_sloc,idv_slev_type,idv_sloc_type, | idv_slocname character(len=120),save :: ncfile type(netcdf_file),save :: ncxyut real,parameter :: hmf2_flag=888., integ_flag=8888., | indep_flag=88888. ! ! Locals not saved: integer :: istat,ipos,ids1(1),ids2(2),ids3(3),icount3(3), | istart3(3),istart2(2),icount2(2),idv,ilev character(len=80) :: char80 character(len=120) :: char120 character(len=8) :: fname real :: rloc(2) real,allocatable :: spvals(:) ! ! Exec: ncalls = ncalls+1 ! write(6,"('Enter wrcdf_xyut: ncalls=',i3)") ncalls ncfile = flnm_cdf_xyut ntimes = xyut%nx if (ncalls == 1) then ! first call: set up new netcdf file ncxyut%filename = ncfile call nccreate(ncxyut%filename,ncxyut%ncid) ncid = ncxyut%ncid write(6,"('wrcdf_xyut: created new netcdf file ',a,' ncid=', | i3)") trim(ncfile),ncid ! ! Define dimensions: ! ! Time is the unlimited dimension: istat = nf_def_dim(ncid,"time",NF_UNLIMITED,ncxyut%id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining time dimension') ! ! Time (total model time in minutes): ids1(1) = ncxyut%id_time istat = nf_def_var(ncid,"time",NF_INT,1,ids1,ncxyut%idv_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_xyut: Error defining time dimension variable') write(char80,"('total model time (day*24*60+hour*60+min)')") istat = nf_put_att_text(ncid,ncxyut%idv_time,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_xyut: Error defining long_name of time dimension ' | //'variable') istat = nf_put_att_text(ncid,ncxyut%idv_time,"units",7, | 'minutes') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_xyut: Error defining units of time dimension variable') ! ! Char length of slev type (character position dimension): istat = nf_def_dim(ncid,"len_levtype",8,id_char8) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining len_levtype dimension') ! ! Char length of location name (character position dimension): istat = nf_def_dim(ncid,"len_locname",56,id_char56) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining len_locname dimension') ! ! Selected levels (dimension nslev and var slev): istat = nf_def_dim(ncid,"nslev",xyut%nslev,id_nslev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining nslev dimension') ids1(1) = id_nslev istat = nf_def_var(ncid,"slev",NF_FLOAT,1,ids1,idv_slev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining slev variable') write(char80, | "('Selected levels (type may be zp, ht, or hmf2)')") istat = nf_put_att_text(ncid,idv_slev,"long_name", | len_trim(char80),trim(char80)) ! ! hmf2_flag, integ_flag, indep_flag are attributes of slev var: istat = nf_put_att_real(ncid,idv_slev,'hmf2_flag',NF_FLOAT, | 1,hmf2_flag) ! ! Level types: ids2(1) = id_char8 ids2(2) = id_nslev istat = nf_def_var(ncid,"slev_type",NF_CHAR,2,ids2, | idv_slev_type) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining slev_type variable') write(char80,"('Selected level type (ht, zp, or hmf2)')") istat = nf_put_att_text(ncid,idv_slev_type,"long_name", | len_trim(char80),trim(char80)) ! ! Explanations of level types are attributes of slev_type: write(char80,"('Level type zp refers to log pressure', | ' (ln(p0/p)')") istat = nf_put_att_text(ncid,idv_slev_type,"zp", | len_trim(char80),trim(char80)) ! write(char80,"('Level type ht refers to height (km)')") istat = nf_put_att_text(ncid,idv_slev_type,"ht", | len_trim(char80),trim(char80)) ! write(char80,"('Level type hmf2 refers to height of ', | 'f2 layer (km)')") istat = nf_put_att_text(ncid,idv_slev_type,"hmf2", | len_trim(char80),trim(char80)) ! ! Selected locations: istat = nf_def_dim(ncid,"nsloc",xyut%nsloc,id_nsloc) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining nsloc dimension') istat = nf_def_dim(ncid,"ncoords_loc",2,id_ncoords) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining ncoords_loc dimension') ids2(1) = id_ncoords ids2(2) = id_nsloc istat = nf_def_var(ncid,"sloc",NF_FLOAT,2,ids2,idv_sloc) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining sloc variable') write(char80,"('Selected locations (type may be [lat,lon]', | ' [lat,slt], [lat,zm], or [gm,gm])')") istat = nf_put_att_text(ncid,idv_sloc,"long_name", | len_trim(char80),trim(char80)) ! ! zmflag and gmflag are attributes of sloc var: istat = nf_put_att_real(ncid,idv_sloc,'zm_flag',NF_FLOAT, | 1,zmflag) istat = nf_put_att_real(ncid,idv_sloc,'gm_flag',NF_FLOAT, | 1,gmflag) ! ! Location names: ids2(1) = id_char56 ids2(2) = id_nsloc istat = nf_def_var(ncid,"sloc_name",NF_CHAR,2,ids2,idv_slocname) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining slocname variable') write(char80,"('Selected location names')") istat = nf_put_att_text(ncid,idv_slocname,"long_name", | len_trim(char80),trim(char80)) ! ! Location types: ids2(1) = id_char8 ids2(2) = id_nsloc istat = nf_def_var(ncid,"sloc_type",NF_CHAR,2,ids2, | idv_sloc_type) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining sloc_type variable') write(char120,"('Selected location type (lat_lon,', | ' lat_zm, lat_slt, or gm_gm)')") istat = nf_put_att_text(ncid,idv_sloc_type,"long_name", | len_trim(char120),trim(char120)) ! ! Explanations of location types are attributes of sloc_type: write(char120,"('Location type lat_lon refers to ', | 'latitude (-90S to +90N), longitude (-180W to +180E)')") istat = nf_put_att_text(ncid,idv_sloc_type,"lat_lon", | len_trim(char120),trim(char120)) ! write(char120,"('Location type lat_zm refers to ', | 'latitude (-90S to +90N), zonal mean')") istat = nf_put_att_text(ncid,idv_sloc_type,"lat_zm", | len_trim(char120),trim(char120)) ! write(char120,"('Location type lat_slt refers to ', | 'latitude (-90S to +90N), solar local time (hrs)')") istat = nf_put_att_text(ncid,idv_sloc_type,"lat_slt", | len_trim(char120),trim(char120)) ! write(char120,"('Location type gm_gm refers to ', | 'global mean')") istat = nf_put_att_text(ncid,idv_sloc_type,"gm_gm", | len_trim(char120),trim(char120)) ! ! Define history vars and global attributes: call ncdefine(ncid,h,ncxyut) ! ! Exit define mode: istat = nf_enddef(ncid) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_enddef') ncxyut%fids(:) = 0 ! init else ! ncalls > 1 ! ! Re-open file for writing: istat = nf_open(ncxyut%filename,NF_WRITE,ncid) if (istat /= NF_NOERR) then write(char80,"('Error reopening file ',a)") | trim(ncxyut%filename) call handle_ncerr(istat,char80) endif endif ! ncalls==1 ! ! Update history variables for new history (do not write field data): if (histonly) then call ncvalues(ncid,h,ncxyut,itime) goto 100 endif ! ! Define current field, long name and units if it's not on the file yet: ! Fields will be dimensioned (nloc,nlev,ntimes) fname = xyut%sname ipos = index(fname,'+') if (ipos > 0) fname(ipos:ipos) = 'P' ipos = index(fname,'/') if (ipos > 0) fname(ipos:ipos) = '_' istat = nf_inq_varid(ncid,fname,ncxyut%fids(nflds)) if (istat /= NF_NOERR) then ! error means field was not found istat = nf_redef(ncid) ! put in define mode ids3(1) = id_nsloc ids3(2) = id_nslev ids3(3) = ncxyut%id_time ! unlimited dim must be in last place istat=nf_def_var(ncid,fname,NF_FLOAT,3,ids3,ncxyut%fids(nflds)) if (istat /= NF_NOERR) then write(char80,"('wrcdf_xyut: Error defining 3d field var ',a)") | trim(fname) call handle_ncerr(istat,char80) else write(6,"('wrcdf_xyut: defined field ',a,': ifld =',i3, | ' id=',i3,' nsloc=',i3,' nslev=',i3,' ntimes=',i4)") | trim(fname),nflds,ncxyut%fids(nflds),xyut%nsloc,xyut%nslev, | ntimes endif ! ! Long name attribute: istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"long_name", | len_trim(xyut%fname),trim(xyut%fname)) if (istat /= NF_NOERR) then write(char80,"('Error defining long name of field ',a)") | trim(xyut%fname) call handle_ncerr(istat,char80) endif ! ! Units attribute: istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"units", | len_trim(xyut%funits),trim(xyut%funits)) if (istat /= NF_NOERR) then write(char80,"('Error defining units of field ',a)") | trim(xyut%fname) call handle_ncerr(istat,char80) endif ! ! Field type: if (xyut%levtype/='indep'.and.xyut%levtype/='integ') then istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"field_type", | len_trim(xyut%ftype),trim(xyut%ftype)) elseif (xyut%levtype=='indep') then istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"field_type", | 8,'HT-INDEP') elseif (xyut%levtype=='integ') then istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"field_type", | 8,'HT-INTEG') endif if (istat /= NF_NOERR) then write(char80,"('Error defining field type ',a,' for field ', | a)") trim(xyut%ftype),trim(xyut%fname) call handle_ncerr(istat,char80) endif ! ! If field was calculated at midpoint levels (f%zptype=='MIDPOINTS'), then ! the processor interpolated the zp dimension to provide the field ! at interface levels (zp_interp=='yes'). Conversely if f%zptype='INTERFACES', ! then zp_interp=='no' (only fields NE,W,Z,POTEN). Note all "unknown" ! fields (secondary histories) are zp_interp=='no', i.e., the processor ! did not interpolate in zp. ! if (xyut%known) then ! known field if (trim(xyut%zptype)=="MIDPOINTS") then istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"zp_interp", | 3,'yes') elseif (trim(xyut%zptype)=="INTERFACES") then istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"zp_interp", | 2,'no') else write(6,"('>>> mkcdf: unknown xyut%zptype=',a)") | trim(xyut%zptype) istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"zp_interp", | 7,'unknown') endif else ! unknown field istat = nf_put_att_text(ncid,ncxyut%fids(nflds),"zp_interp", | 2,'no') endif if (istat /= NF_NOERR) then write(char120,"('Error return from nf_put_att_text', | ' of attribute zp_interp for field ',a,': known=',l1, | ' zptype=',a )") fname,xyut%known,xyut%zptype call handle_ncerr(istat,char120) endif ! ! Field index: nflds = nflds+1 ! ! Exit define mode: istat = nf_enddef(ncid) endif ! define current field ! write(6,"(/,'wrcdf_xyut: sname=',a,' ntimes=',i4)") ! | xyut%sname,ntimes ! write(6,"(' glat=',f9.2,' glon=',f9.2,' slt=',f9.2)") ! | xyut%glat,xyut%glon,xyut%slt ! write(6,"(' iloc=',i3,' loctype=',a,' locname=',a)") xyut%iloc, ! | xyut%loctype,trim(xyut%locname) ! write(6,"(' ilev=',i3,' zp=',f9.2,' ht=',f9.2,' levtype=',a)") ! | xyut%ilev,xyut%zp,xyut%ht,xyut%levtype ! write(6,"(' ncurve=',i3,' data=',/,(6e12.4))") ! | xyut%ncurve,xyut%data(:,1) ! ! Write current selected level: if (nflds-1==1.and.xyut%iloc==1) then select case (xyut%levtype) case ('zp ') istat = nf_put_var1_real(ncid,idv_slev,xyut%ilev,xyut%zp) case ('ht ') istat = nf_put_var1_real(ncid,idv_slev,xyut%ilev,xyut%ht) case ('hmf2 ') istat = nf_put_var1_real(ncid,idv_slev,xyut%ilev,hmf2_flag) case ('integ') istat = nf_put_var1_real(ncid,idv_slev,xyut%ilev,integ_flag) case ('indep') istat = nf_put_var1_real(ncid,idv_slev,xyut%ilev,indep_flag) case default write(6,"('>>> wrcdf_xyut: unknown levtype=',a)") | xyut%levtype end select if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_var1_real giving value ', | 'to slev var: levtype=',a,' ilev=',i3)") xyut%levtype, | xyut%ilev call handle_ncerr(istat,char120) endif ! ! Write 8-character level type (8,ilev): istart2(1) = 1 istart2(2) = xyut%ilev icount2(1) = len_trim(xyut%levtype)+1 ! extra for null terminate icount2(2) = 1 char80 = ' ' ; char80 = xyut%levtype ipos = len_trim(char80)+1 char80(ipos:ipos) = char(0) ! null terminate istat = nf_put_vara_text(ncid,idv_slev_type,istart2,icount2, | char80(1:8),8) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_text giving value ', | 'to slev_type var: levtype=',a,' ilev=',i3)") xyut%levtype, | xyut%ilev call handle_ncerr(istat,char120) endif endif ! nflds-1==1 and iloc==1 ! ! Write current selected location (ncoords==2,nsloc): if (nflds-1==1.and.xyut%ilev==1) then rloc(1) = xyut%glat istart2(1) = 1 icount2(1) = 2 ! 2 coords (e.g., [lat,lon]) icount2(2) = 1 ! 1 location (current) istart2(2) = xyut%iloc ! current location index select case (xyut%loctype) case ('lon') ! [lat,lon] ( latitude, longitude) rloc(2) = xyut%glon istat = nf_put_vara_real(ncid,idv_sloc,istart2,icount2,rloc) case ('slt') ! [lat,slt] (latitude, local time) rloc(2) = xyut%slt istat = nf_put_vara_real(ncid,idv_sloc,istart2,icount2,rloc) case ('zm ') ! [lat,zm] (zonal mean) rloc(2) = zmflag istat = nf_put_vara_real(ncid,idv_sloc,istart2,icount2,rloc) case ('gm ') ! [gm,gm] (global mean) rloc(1) = gmflag rloc(2) = gmflag istat = nf_put_vara_real(ncid,idv_sloc,istart2,icount2,rloc) case default write(6,"('>>> wrcdf_xyut: unknown loctype=',a)") | xyut%loctype end select if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_real giving value ', | 'to sloc var: loctype=',a,' iloc=',i3)") xyut%loctype, | xyut%iloc call handle_ncerr(istat,char120) endif ! ! Write 8-character location type (8,iloc): istart2(1) = 1 istart2(2) = xyut%iloc icount2(2) = 1 select case (xyut%loctype) case ('lon') ! lat_lon ( latitude, longitude) icount2(1) = 8 istat = nf_put_vara_text(ncid,idv_sloc_type,istart2,icount2, | 'lat_lon'//char(0),8) case ('slt') ! lat_slt (latitude, local time) icount2(1) = 8 istat = nf_put_vara_text(ncid,idv_sloc_type,istart2,icount2, | 'lat_slt'//char(0),8) case ('zm ') ! lat_zm (zonal mean) icount2(1) = 7 istat = nf_put_vara_text(ncid,idv_sloc_type,istart2,icount2, | 'lat_zm'//char(0),7) case ('gm ') ! gm_gm (global mean) icount2(1) = 6 istat = nf_put_vara_text(ncid,idv_sloc_type,istart2,icount2, | 'gm_gm'//char(0),6) case default write(6,"('>>> wrcdf_xyut: unknown loctype=',a)") | xyut%loctype end select if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_text giving value ', | 'to sloc_type var: loctype=',a,' iloc=',i3)") xyut%loctype, | xyut%iloc call handle_ncerr(istat,char120) endif ! ! Write 56-char location name (56,iloc): if (len_trim(xyut%locname) > 0) then istart2(1) = 1 istart2(2) = xyut%iloc icount2(1) = len_trim(xyut%locname)+1 ! extra for null terminate icount2(2) = 1 char80 = ' ' ; char80 = xyut%locname ipos = len_trim(char80)+1 char80(ipos:ipos) = char(0) ! null terminate istat = nf_put_vara_text(ncid,idv_slocname,istart2,icount2, | char80(1:56),56) if (istat /= NF_NOERR) then write(char120,"('Error from nf_put_vara_text giving value ', | 'to slocname var: locname=',a,' iloc=',i3)") | trim(xyut%locname),xyut%iloc call handle_ncerr(istat,char120) endif endif ! len_trim(xyut%locname) > 0 endif ! nflds-1==1.and.xyut%ilev==1 ! ! Write field data to netcdf file: ! Fields are dimensioned (xyut%nsloc,xyut%nslev,ntimes). idv = ncxyut%fids(nflds-1) istart3(1) = xyut%iloc istart3(2) = xyut%ilev istart3(3) = 1 icount3(1:2) = 1 icount3(3) = ntimes istat = nf_put_vara_real(ncid,idv,istart3,icount3, | xyut%data) if (istat /= NF_NOERR) then write(char80,"('Error giving values to field ',a,' iloc=',i3, | ' ilev=',i3,' ntimes=',i4)") xyut%sname,xyut%iloc,xyut%ilev, | ntimes call handle_ncerr(istat,char80) endif ! ! If field is height-integrated or height-independent, pad remaining selected ! levels with missing value: if ((xyut%levtype=='indep'.or.xyut%levtype=='integ').and. | xyut%nslev > 1.and.xyut%ilev==1) then allocate(spvals(ntimes),stat=istat) spvals(:) = spval do ilev=2,xyut%nslev istart3(1) = xyut%iloc ! start location istart3(2) = ilev ! start level istart3(3) = 1 ! start time icount3(1) = 1 ! 1 location icount3(2) = 1 ! 1 level icount3(3) = ntimes ! number of times istat = nf_put_vara_real(ncid,idv,istart3,icount3,spvals) if (istat /= NF_NOERR) then write(char80,"('Error padding missing vals for indep or ', | 'integ field ',a,' iloc=',i3,' ilev=',i3,' ntimes=',i4)") | xyut%sname,xyut%iloc,ilev,ntimes call handle_ncerr(istat,char80) endif enddo ! ilev=2,nslev endif ! pad integ or indep field ! 100 continue istat = nf_close(ncid) if (itime < ntimes) | write(6,"('Updated netcdf file ',a,' itime=',i4,' mtime=',3i4)") | trim(ncfile),itime,xyut%mtimes(:,itime) end subroutine wrcdf_xyut !----------------------------------------------------------------------- subroutine doxyut(xyut,h,hcntr,itime,u,v) use proc,only: spval,nppf,iframe use input,only: ixyut_log10,multiplt,multiadvfr,ipltrowcol, + ilab_hq,iplot use plt ! ! Args: type(xyut_type),intent(inout) :: xyut type(history),intent(in) :: h,hcntr real,optional,intent(in) :: u(size(xyut%data,1)), + v(size(xyut%data,1)) integer,intent(in) :: itime ! ! Locals: real :: vp(4),vpxyut(4),vscmag,dum character(len=80) msgout integer :: isltax real :: toffset=.04, boffset,yy(2),yboffset ! if (all(xyut%data==spval)) then write(6,"('>>> Skipping field ',a,' because all values are ', + 'missing (spval).')") xyut%sname return endif ! ! Set up viewport (for multiplt): if (iplot > 0) then vpxyut = (/.18,.93,.30,.90/) if (present(u).and.present(v)) vpxyut = (/.15,.90,.40,.90/) call setmultivp(vpxyut,iadvfr,nppf+1,multiplt,ipltrowcol,vp) ! ! Make labels and draw the xy line plot: call mkxyutlabs(xyut,h,hcntr,msgout) ! ! If u,v are present, make 1d vector plot from u+v: if (present(u).and.present(v)) then ! write(6,"('doxyut: u=',/,(6e12.4))") u ! write(6,"('doxyut: v=',/,(6e12.4))") v ! ! subroutine pltvec1d(u,v,xx,nx,vscalemag,vscmag,yboffset,vp, ! + spval) ! yboffset is distance below bottom of viewport to draw scale arrow, ! expressed as a fraction of the viewport height. ! yboffset = .55 ! fraction of viewport height call pltvec1d(u,v,xyut%xx,xyut%nx,0.,vscmag,yboffset,vp, + spval) yy(1) = -vscmag yy(2) = vscmag call labrect(xyut%xx,0,yy,2,' ',xyut%sname,.03) isltax = 0 ; boffset=.2 if (xyut%loctype=='lon') then isltax = 1 ; boffset = .4 endif call labutxy(xyut%mtimes,xyut%nx,xyut%data(:,1),0,'DUM',0., + isltax,xyut%glon) call wrlab6(xyut%tlabs,toffset,tlabchsz, + xyut%blabs,boffset,blabchsz,vp,ilab_hq) ! ! If u,v are not present, make line plot from xyut%data: else if (xyut%nx > 1) then call pltxyut(xyut,vp) else write(6,"(/,'>>> xyut: cannot make plots because ntimes=', | i4,' (ntimes must be > 1)')") xyut%nx endif endif ! ! 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(xyut%data,size(xyut%data),zmin,zmax) endif ! ! Write output data files: call wrout_xyut(xyut,h,hcntr,itime) end subroutine doxyut !----------------------------------------------------------------------- subroutine savehts_xyut(f,fcntr,nf,h,xyut,itime,ixz_xyutdat, | ixhmf2,zptype) use proc use fields use input ! ! Save heights in xyutdat(itime,:,:,ixz_xyutdat): ! (if diffs, save control heights in xyutdat_cntr(itime,:,:,ixz_xyutdat) ! Also save hmf2 heights in xyutdat_hmf2(ntimes,nloc) and ! xyutdat_cntr_hmf2(ntimes,nloc) if necessary. ! ! 2/98: this routine generates internal f90 compiler error w/ ! -g on winterpark sgi (is ok w/o -g). ! ! 1/23/08 btf: Use Z consistent w/ current fields zptype ! (midpoints or interfaces) ! ! Args: integer,intent(in) :: nf,itime,ixz_xyutdat,ixhmf2 type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h type(xyut_type),intent(inout) :: xyut character(len=*),intent(in) :: zptype ! ! Locals: integer :: ixz_cntr,iloc,l,k,ixlat,ixlon ! ! External: integer,external :: ixfind,ixfindc real,external :: fglbm,fmean ! iloc = 0 loc_loop0: do l=1,mxloc if (xyut_locs(1,l)==spval.or.xyut_locs(2,l)==spval) + cycle loc_loop0 iloc = iloc+1 call setloctype(xyut_locs(1,l),xyut_locs(2,l),h%ut, + xyut%glat,xyut%glon,xyut%slt,ixlat,ixlon,xyut%loctype) ! ! Save heights at lat,lon/slt: if (xyut%loctype=='lon'.or.xyut%loctype=='slt') then if (trim(zptype)=='INTERFACES') then xyutdat(itime,iloc,:,ixz_xyutdat) = z_ifaces(ixlon,ixlat,:) if (diffs) xyutdat_cntr(itime,iloc,:,ixz_xyutdat) = + zcntr_ifaces(ixlon,ixlat,:) else xyutdat(itime,iloc,:,ixz_xyutdat) = z_midpts(ixlon,ixlat,:) if (diffs) xyutdat_cntr(itime,iloc,:,ixz_xyutdat) = + zcntr_midpts(ixlon,ixlat,:) endif if (ixhmf2>0) then xyutdat_hmf2(itime,iloc) = f(ixhmf2)%data(ixlon,ixlat,1) if (diffs) xyutdat_cntr_hmf2(itime,iloc) = + fcntr(ixhmf2)%data(ixlon,ixlat,1) endif ! ! Save heights at zm: elseif (xyut%loctype== 'zm ') then if (trim(zptype)=='INTERFACES') then do k=1,npress xyutdat(itime,iloc,k,ixz_xyutdat) = + fmean(z_ifaces(:,ixlat,k),nlon,spval,1) if (diffs) + xyutdat_cntr(itime,iloc,k,ixz_xyutdat) = + fmean(zcntr_ifaces(:,ixlat,k),nlon,spval,0) enddo else do k=1,npress xyutdat(itime,iloc,k,ixz_xyutdat) = + fmean(z_midpts(:,ixlat,k),nlon,spval,1) if (diffs) xyutdat_cntr(itime,iloc,k,ixz_xyutdat) = + fmean(zcntr_midpts(:,ixlat,k),nlon,spval,0) enddo endif if (ixhmf2>0) then xyutdat_hmf2(itime,iloc) = + fmean(f(ixhmf2)%data(:,ixlat,1),nlon,spval,1) if (diffs) xyutdat_cntr_hmf2(itime,iloc) = + fmean(fcntr(ixhmf2)%data(:,ixlat,1),nlon,spval,0) endif ! ! Save heights at gm: elseif (xyut%loctype=='gm ') then if (trim(zptype)=='INTERFACES') then do k=1,npress xyutdat(itime,iloc,k,ixz_xyutdat) = + fglbm(z_ifaces(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) if (diffs) xyutdat_cntr(itime,iloc,k,ixz_xyutdat) = + fglbm(zcntr_ifaces(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo else do k=1,npress xyutdat(itime,iloc,k,ixz_xyutdat) = + fglbm(z_midpts(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) if (diffs) xyutdat_cntr(itime,iloc,k,ixz_xyutdat) = + fglbm(zcntr_midpts(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif if (ixhmf2>0) then xyutdat_hmf2(itime,iloc) = + fglbm(f(ixhmf2)%data(:,:,1),nlon,nlat,gcmlat,dlat, + dlon,spval) if (diffs) + xyutdat_cntr_hmf2(itime,iloc) = + fglbm(fcntr(ixhmf2)%data(:,:,1),nlon,nlat,gcmlat,dlat, + dlon,spval) endif endif ! write(6,"('savehts_xyut: zptype=',a,' itime=',i4,' iloc=',i3, ! | ' xyutdat(itime,iloc,:,ixz_xyutdat)=',/,(6f9.2))") ! | trim(zptype),itime,iloc,xyutdat(itime,iloc,:,ixz_xyutdat) enddo loc_loop0 ! save heights end subroutine savehts_xyut !----------------------------------------------------------------------- subroutine alloc_xyutdat_hmf2(f,nf,ixhmf2,ier) ! ! Allocate xyutdat_hmf2(ntimes,nloc) if necessary: ! use proc use fields use input ! ! Args: integer,intent(out) :: ixhmf2,ier integer,intent(in) :: nf type(field),intent(in) :: f(nf) ! ! Local: integer :: i,k,need_hmf2,nloc,istat ! ! External: integer,external :: ixfindc character(len=16),external :: float_to_str ! ! Include hmf2 heights if necessary: ier = 0 ; need_hmf2 = 0 do k=1,mxzpht if (xyut_zpht(k) /= spval.and. + (trim(float_to_str(xyut_zpht(k)))=='hmf2'.or. + trim(float_to_str(xyut_zpht(k)))=='HMF2'.or. | xyut_zpht(k)==hmf2flag)) + need_hmf2=k ! will need hmf2 heights enddo ! ! If hmf2 heights are needed, then ixhmf2 is index to HMF2 in f(nf). if (need_hmf2 > 0) then ixhmf2 = ixfindc(f%fname8,nf,'HMF2') if (ixhmf2 <= 0.or..not.associated(f(ixhmf2)%data)) then write(6,"('>>> WARNING alloc_xyutdat_hmf2: need hmf2', + ' heights, but cannot find HMF2 in fields (or hmf2 ', + 'data not associated).')") xyut_zpht(need_hmf2) = spval ixhmf2 = 0 ier = -1 return endif nloc = 0 ! number of valid locations do i=1,mxloc if (xyut_locs(1,i)/=spval.and.xyut_locs(2,i)/=spval) + nloc = nloc+1 enddo if (nloc==0) then write(6,"('>>> alloc_xyutdat_hmf2: need locations xyut_locs ', + '(nloc==0).')") ier = 1 return endif allocate(xyutdat_hmf2(ntimes,nloc),stat=istat) if (istat /= 0) then write(6,"('>>> alloc_xyutdat_hmf2: error allocating ', + 'xyutdat_hmf2(ntimes,nloc) where: ntimes=',i4, + ' nloc=',i3)") ntimes,nloc call allocerr(istat, + 'alloc_xyutdat_hmf2 allocating xyutdat_hmf2') else write(6,"('Allocated xyutdat_hmf2(ntimes,nloc), where:',/, + ' ntimes=',i4,' nloc=',i3)") ntimes,nloc endif if (diffs) then allocate(xyutdat_cntr_hmf2(ntimes,nloc),stat=istat) if (istat /= 0) then write(6,"('>>> alloc_xyutdat_hmf2: error allocating ', + 'xyutdat_cntr_hmf2(ntimes,nloc) where: ntimes=',i4, + ' nloc=',i3)") ntimes,nloc call allocerr(istat, + 'alloc_xyutdat_hmf2 allocating xyutdat_cntr_hmf2') else write(6,"('Allocated xyutdat_cntr_hmf2(ntimes,nloc),', + ' where:',/,' ntimes=',i4,' nloc=',i3)") ntimes,nloc endif endif endif end subroutine alloc_xyutdat_hmf2 end module mk_xyut