! module mk_satut ! ! Make ut vs zp/ht contours at selected locations: ! (ut on x-axis, zp or height on y-axis) ! 8/25/04 btf: ! Replaced hard-wired dimension 10000 with mxsatloc. ! Made satutdat and satuthtdat allocatable arrays, local to sub mksatut. ! Fixed cuthtint call for height interpolation. Height interpolation is ! for plots only -- the satut output nc file is at pressures only. ! use fields use hist use get_flds implicit none ! type satut 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=80) :: + 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 satut real :: tlabchsz(3) = (/.02,.02,.02/), + blabchsz(3) = (/.018,.018,.018/) integer,parameter :: mxsatloc = 10000 ! max # of satellite locations ! ! utv is satut type to be plotted type(satut),save :: utv ! ! utvdat(ntimes,npress,nloc,nf) is data from which utv is set up. ! utvdat(itime,:,:,:) is defined each time sub mksatut is called. ! The last time mksatut 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 mksatut(f,fcntr,fprev,nf,h,hcntr,flnm_satcdf,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 pltsatut, 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(field),intent(inout) :: fprev(nf) type(field),intent(in) :: fprev(nf) type(history),intent(in) :: h,hcntr character(len=120),intent(out) :: flnm_satcdf ! ! Locals: integer :: nloc,i,ixf,iloc,l,ifld,nfld,ier,iht,k,iyax,noh,ifoh, + izpht0,izpht1,nzprange,k0zp,k1zp,iadvfr,ixlat,ixlon,ix, + logint,isltax integer,save :: ixz,ixz_utvdat,nhtscale real,pointer,save :: ! user requested y-axes | htscale(:)=>NULL(), | zprange(:)=>NULL() real :: vp(4),dum,boff real :: vpsatut(4) = (/.12,.85,.29,.90/) real :: toffset=.06, boffset=.35 real,allocatable :: fdat1(:,:),fdat2(:,:) character(len=80) msgout ! ! Array to hold satellite attitude data at each location: real,save,dimension(mxsatloc,9) :: satlocs ! satellite attitude data ! ! Allocatable arrays to hold model data interpolated to satellite track: ! satutdat(mxsatloc,npress,nf) ! at pressures ! satuthtdat(mxsatloc,nhtscale,nf) ! at heights ! real,save,allocatable :: satutdat(:,:,:) real,save,allocatable :: satuthtdat(:,:,:) ! integer,save :: iprev,icur,nsat integer :: nlocs,istat,mins,mins_first,mins_last integer,save :: ncid,mins_prev,ipos real :: mtime_init,mtime_curr ! ! Externals: real,external :: fslt,fmean,fglbm 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(:) ! height scale to be used integer :: nhscale ! dimension for hscale real,pointer :: fhscale(:) ! height scale of field integer :: nfhscale ! dimension for fhscale ! ! write(6,"('Enter mksatut: itime=',i3)") itime ! 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 ! ! Read in satellite location data (lat/lon/times) into satlocs array mins = ((mtimes(1,itime)-1)*24. + + mtimes(2,itime))*60. + + mtimes(3,itime) mins_first = ((mtimes(1,1)-1)*24. + + mtimes(2,1))*60. + + mtimes(3,1) mins_last = ((mtimes(1,ntimes)-1)*24. + + mtimes(2,ntimes))*60. + + mtimes(3,ntimes) ! ! Run initial set up for first time if (itime==1) then icur = 0 ! satutdat = spval if(index(satfile, 'tidi') .gt. 0) then print *, ' assumed reading TIDI data...' call read_tidi_locs_nc(satlocs,nsat,mins_first,mins_last) else if(index(satfile, 'saber') .gt. 0) then print *, ' assumed reading SABER data...' ! call read_saber_locs(satlocs,nsat,mins_first,mins_last) !ascii call read_saber_locs_nc(satlocs,nsat,mins_first,mins_last) !netcdf else if(index(satfile, 'cedar') .gt. 0) then print *, ' assumed reading CEDAR data...' call read_cedar_locs_nc(satlocs,nsat,mins_first,mins_last) else if(index(satfile, 'champ') .gt. 0) then print *, ' assumed reading CHAMP data...' call read_champ_locs_nc(satlocs,nsat,mins_first,mins_last) else if(index(satfile, 'guvi') .gt. 0) then print *, ' assumed reading GUVI data...' call read_guvi_locs_nc(satlocs,nsat,mins_first,mins_last) else print *, ' not sure what type of satellite file in use.' print *, ' use tidi, saber, cedar, champ, guvi' print *, ' in sat path or fname.' stop endif if (nsat > mxsatloc) then write(6,"('>>> mksatut: nsat > mxsatloc: nsat=',i5, | ' mxsatloc=',i5)") nsat,mxsatloc write(6,"('To make this run, increase parameter mxsatloc', | ' (file mksatut.F) and rebuild the processor.')") stop 'mxsatloc' endif mins_prev = mins ixz = ixfindc(f%fname8,nf,'Z ') ! ! Set up vertical scale(s) zp and/or ht: call setvert(satut_zprange,gcmlev,npress,dlev,nzprange, | k0zp,k1zp,zprange,satut_htscale,nhtscale,htscale,spval,1) ! if (nzprange > 0) ! | write(6,"('mksatut after setvert: nzprange=',i3,' zprange=', ! | /,(8f9.2))") nzprange,zprange ! if (nhtscale > 0) ! | write(6,"('mksatut after setvert: nhtscale=',i3,' htscale=', ! | /,(8f9.2))") nhtscale,htscale ! ! Allocate arrays to hold model data interpolated to satellite track: ! allocate(satutdat(mxsatloc,npress,nf),stat=ier) if (ier /= 0) then write(6,"('mksatut error allocating', | ' satutdat(mxsatloc,npress,nf): mxsatloc=', | i5,' npress=',i3,' nf=',i3)") mxsatloc,npress,nf call allocerr(ier,"mksatut allocating satutdat") else write(6,"('mksatut: allocated satutdat(mxsatloc,npress,nf)', | ' where: mxsatloc=',i5,' npress=',i3,' nf=',i3)") | mxsatloc,npress,nf endif ! if (nhtscale > 0) then allocate(satuthtdat(mxsatloc,nhtscale,nf),stat=ier) if (ier /= 0) then write(6,"('mksatut error allocating', | ' satuthtdat(mxsatloc,nhtscale,nf): mxsatloc=', | i5,' nhtscale=',i3,' nf=',i3)") mxsatloc,nhtscale,nf call allocerr(ier,"mksatut allocating satuthtdat") else write(6,"('mksatut: allocated ', | 'satuthtdat(mxsatloc,nhtscale,nf) where: mxsatloc=', | i5,' nhtscale=',i3,' nf=',i3)") mxsatloc,nhtscale,nf endif endif write(6,"('mksatut: nsat=',i5)") nsat ! ! Init netcdf output file: call mksatcdf(satutdat(1,:,:),satlocs(1,:),1,ntimes,npress, | nf,flnm_satcdf,f,h,ncid,0,ier,satfile) endif ! itime==1 ! ! Find sat points within previous and current model times: iprev = icur+1 mtime_init = (mtimes(1,1)*24.*60. + + mtimes(2,1)*60. + + mtimes(3,1)) mtime_curr = ((mtimes(1,itime)-1)*24.*60. + + mtimes(2,itime)*60. + + mtimes(3,itime)) do i=iprev,nsat if (satlocs(i,1) .lt. mtime_curr) icur = i enddo ! i nlocs = icur - iprev + 1 ! write(6,"('mksatut: itime=',i5,' iprev=',i5,' icur=',i5, ! | ' nlocs=',i5,' nsat=',i5,' mtime_curr=',f12.2)") ! | itime,iprev,icur,nlocs,nsat,mtime_curr ! ! Interpolate model data to satlocs(iprev:icur) ! (when itime==1, iprev=1 and icur=0, so this loop is not executed when itime==1) do i=iprev,icur call interp_sat(f,nf,fprev,mins_prev,mins,npress, + satlocs(i,:),satutdat(i,:,:)) ! ! Interpolate to linear height scale: ! do ixf=1,nf if (nhtscale > 0) then logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 ! ! subroutine cuthtint(fin,fht,idim1,nzp,fout,hts,nhts,logint,spval,ier,iprnt) ! dimension fin(idim1,nzp),fht(idim1,nzp),fout(idim1,nhts),hts(nhts) ! call cuthtint(satutdat(i,:,ixf),satutdat(i,:,ixz),1, | npress,satuthtdat(i,:,ixf),htscale,nhtscale,logint, | spval,ier,1) ! write(6,"('Interpolated field ',i3,' to heights: ', ! | ' nhtscale=',i3,' satuthtdat(i,:,ixf) min,max=', ! | 2e12.4)") ixf,nhtscale,minval(satuthtdat(i,:,ixf)), ! | maxval(satuthtdat(i,:,ixf)) endif enddo ! ixf ! ! Update netcdf file for satellite interpolated data satlocs(i,1) = mod(satlocs(i,1),24.*60.) ! set to min of day call mksatcdf(satutdat(i,:,:),satlocs(i,:),i,ntimes,npress,nf, | flnm_satcdf,f,h,ncid,1,ier,satfile) if (ier /= 0) then write(6,"('>>> tgcmproc: error return from mksatcdf -- ', + 'will not make satellite netcdf file.')") sendsat = ' ' endif enddo ! i ! ! Save data for next input file mins_prev = mins ! ! If done with model files, then close netcdf file: if(itime==ntimes) then call mksatcdf(satutdat(icur,:,:),satlocs(icur,:),icur, + ntimes,npress,nf,flnm_satcdf,f,h,ncid,2,ier,satfile) istat = nf_close(ncid) if (istat /= NF_NOERR) then call handle_ncerr(istat, + 'Error return from nf_close to close netcdf file',ier) else write(6,"('Closed netcdf file ',a)") trim(flnm_satcdf) endif endif ! ! Allocate utvdat if this is first call: if (itime==1) then nloc = 0 ! number of valid locations do i=1,mxloc if (satut_locs(1,i)/=spval.and.satut_locs(2,i)/=spval) + nloc = nloc+1 enddo if (nloc==0) then write(6,"('>>> mksatut: need locations satut_locs ', + '(nloc==0). Turning off ipltsatut.')") ipltsatut = 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 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,"('>>> mksatut: need fields: nfld==noh==0.', + ' Turning ipltsatut off.')") ipltsatut = 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,"('mksatut error allocating', + ' utvdat(ntimes,npress,nloc,nfld): ntimes=', + i4,' npress=',i3,' nloc=',i3,' nfld=',i3)") + ntimes,npress,nloc,nfld call allocerr(ier,"mksatut allocating utvdat") endif if (diffs) then allocate(utvdat_cntr(ntimes,npress,nloc,nfld),stat=ier) if (ier /= 0) then write(6,"('mksatut error allocating', + ' utvdat_cntr(ntimes,npress,nloc,nfld): ntimes=', + i4,' npress=',i3,' nloc=',i3,' nfld=',i3)") + ntimes,npress,nloc,nfld call allocerr(ier,"mksatut 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,"('mksatut error allocating', + ' utvdat_oh(ntimes,nohalt,nloc,nfld): ntimes=', + i4,' nohalt=',i3,' nloc=',i3,' noh=',i3)") + ntimes,nohalt,nloc,noh call allocerr(ier,"mksatut allocating utvdat_oh") endif if (diffs) then allocate(utvdat_cntr_oh(ntimes,nohalt,nloc,noh),stat=ier) if (ier /= 0) then write(6,"('mksatut error allocating', + ' utvdat_cntr_oh(ntimes,nohalt,nloc,nfld): ntimes=', + i4,' nohalt=',i3,' nloc=',i3,' noh=',i3)") + ntimes,nohalt,nloc,noh call allocerr(ier,"mksatut allocating utvdat_cntr_oh") endif endif endif ! ! Locate index to heights field (needed for interp, etc): ! ixz = ixfindc(f%fname8,nf,'Z ') if (ixz <= 0) then write(6,"('>>> mksatut: need heights in fields: field ', + 'names = ',(8a8))") f%fname8 stop 'mksatut z' 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,"mksatut allocating utv%xx") utv%xlab = 'UT (HRS)' allocate(utv%mtimes(3,ntimes),stat=ier) if (ier/=0) + call allocerr(ier,"mksatut allocating utv%mtimes") utv%hvol0 = h%mssvol if (diffs) utv%hvol0_cntr = hcntr%mssvol 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 utv%hvol1 = h%mssvol if (diffs) utv%hvol1_cntr = hcntr%mssvol 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 ! ! Save heights in utvdat(:,:,:,ixz_utvdat) for use in cuthtint, etc below. iloc = 0 ! loc index to utvdat(:,:,iloc,:) loc_loop0: do l=1,mxloc if (satut_locs(1,l)==spval.or.satut_locs(2,l)==spval) + cycle loc_loop0 iloc = iloc+1 call setlocsat(h,satut_locs(1,l),satut_locs(2,l),ixlat,ixlon) if (utv%loctype /= 'zm '.and.utv%loctype /= 'gm ')then utvdat(itime,:,iloc,ixz_utvdat) = f(ixz)%data(ixlon,ixlat,:) if (diffs) utvdat_cntr(itime,:,iloc,ixz_utvdat) = + fcntr(ixz)%data(ixlon,ixlat,:) elseif (utv%loctype == 'zm ') then do k=1,npress utvdat(itime,k,iloc,ixz_utvdat) = + fmean(f(ixz)%data(:,ixlat,k),nlon,spval,0) enddo if (diffs) then do k=1,npress utvdat_cntr(itime,k,iloc,ixz_utvdat) = + fmean(fcntr(ixz)%data(:,ixlat,k),nlon,spval,0) enddo endif else do k=1,npress utvdat(itime,k,iloc,ixz_utvdat) = + fglbm(f(ixz)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo if (diffs) then do k=1,npress utvdat_cntr(itime,k,iloc,ixz_utvdat) = + fglbm(fcntr(ixz)%data(:,:,k),nlon,nlat,gcmlat,dlat, + dlon,spval) enddo endif endif enddo loc_loop0 ! ! Field loop: ifld = 0 ; ifoh = 0 fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or. + f(ixf)%nlev==1) cycle fields_loop if (trim(f(ixf)%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)=='MID') then ! midpoints call setvert(satut_zprange,gcmlev,npress,dlev,nzprange, + k0zp,k1zp,zprange,satut_htscale,nhtscale,htscale,spval,1) else ! interfaces call setvert(satut_zprange,gcmilev,npress,dlev,nzprange, + k0zp,k1zp,zprange,satut_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%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 (satut_locs(1,l)==spval.or.satut_locs(2,l)==spval) + cycle loc_loop iloc = iloc+1 call setlocsat(h,satut_locs(1,l),satut_locs(2,l),ixlat,ixlon) ! ! 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,:) 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 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 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 = satut_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 (isatut_yaxright > 0) vpsatut(2) = .77 ! ! Define y coords in pressure: utv%ny = nzprange allocate(utv%yy(utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mksatut allocating utv%yy") utv%yy = zprange ! array op write(utv%ylab,"('LN(P0/P)')") ! ! 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,"mksatut 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,"mksatut allocating fdat1") fdat1 = utvdat(:,k0zp:k1zp,iloc,ifld) allocate(fdat2(ntimes,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mksatut 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 vpsatut(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,"mksatut 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,"mksatut 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,"mksatut 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,"mksatut allocating fdat2") ! call cuthtint(utvdat_cntr(:,:,iloc,ifld), ! + utvdat_cntr(:,:,iloc,ixz_utvdat),ntimes,npress, ! + fdat2,htscale,nhtscale,logint,spval,ier,1) 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,"mksatut allocating utv%yy") utv%yy = hscale ! array op allocate(utv%data(utv%nx,utv%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mksatut 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,"mksatut 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 ((isatut_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + isatut_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(vpsatut,iadvfr,nppf+1,multiplt, + ipltrowcol,vp) ! ! Contour: isltax = 0 if(utv%loctype=='lon') isltax = 1 call pltsatut(utv,utvdat(:,:,iloc,ixz_utvdat), + isatut_yaxright,isltax,vp) ! ! 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_satut(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 end subroutine mksatut !------------------------------------------------------------------- subroutine setlocsat(h,reqlat,reqlon,ixlat,ixlon) ! ! Set location of utv, and return ixlat,ixlon: ! ! satut_locs(1,iloc) is latitude, satut_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 satut_locs(:,iloc)=='zm', (i.e., both lat,lon) then do global means) ! ! Args: real(kind=8),intent(in) :: reqlat,reqlon integer,intent(out) :: ixlat,ixlon type(history),intent(in) :: h ! ! Locals: real :: slt,dum ! ! External: logical,external :: isslt real,external :: fslt integer,external :: ixfind ! utv%loctype = ' ' ! ! Global means: if (reqlat==zmflag.and.reqlon==zmflag)then utv%loctype = 'gm ' ! ! Zonal means: elseif (reqlon == zmflag) then ! zonal means utv%loctype = 'zm ' utv%glon = zmflag ! ! Regular longitude or local time: else if (isslt(reqlon,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 = reqlon 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 = reqlat ixlat = ixfind(gcmlat,nlat,utv%glat,dlat) endif ! write(6,"('setlocsat: reqlat=',f9.2,' reqlon=',f9.2, ! | ' utv%glat=',f9.2,' utv%glon=',f9.2,' ixlat=',i3,' ixlon=', ! | i3)") reqlat,reqlon,utv%glat,utv%glon,ixlat,ixlon return end subroutine setlocsat !-------------------NETCDF------------------------------------------ subroutine read_tidi_locs_nc(satlocs,nsat,mins0,mins1) use input,only: satfile integer,intent(in) :: mins0,mins1 real,intent(out) :: satlocs(mxsatloc,9) ! satellite lat/lon/times integer,intent(out) :: nsat ! number of orbit points integer, external :: nextlu real :: orbnum, adflag character (len=1) :: adflag_ch integer :: lu,ios,k character(len=150) :: line character(len=1) :: char1 ! ! netcdf variables integer :: ncdfid, dimprofs, lenprofs,istat,trackid integer :: profid,dateid,timeid,lonid,latid,tpadid,tpdnid,wcsideid integer,dimension(mxsatloc) :: profnum,date,tpad,tpdn,wcside real,dimension(mxsatloc) :: longitude,latitude,time,lst,track lu=nextlu() ! open netcdf ios = NF_OPEN(satfile, NF_NOWRITE, ncdfid) if (ios /= NF_NOERR) then write(6,"('>>> WARNING read_tidi_locs_nc: error opening file ', + a)") satfile stop 'read_tidi_locs_nc' endif print *, ' <<< READ_TIDI_LOCS_NC >>> : opened satfile: ', + ncdfid, ' ', satfile ! get dim IDs and dim length istat=NF_INQ_DIMID(ncdfid,'profs',dimprofs) istat=NF_INQ_DIMLEN(ncdfid,dimprofs,lenprofs) ! get var IDs istat=NF_INQ_VARID(ncdfid,'profnum',profid) istat=NF_INQ_VARID(ncdfid,'date',dateid) istat=NF_INQ_VARID(ncdfid,'time',timeid) istat=NF_INQ_VARID(ncdfid,'lon',lonid) istat=NF_INQ_VARID(ncdfid,'lat',latid) istat=NF_INQ_VARID(ncdfid,'tpad',tpadid) istat=NF_INQ_VARID(ncdfid,'tpdn',tpdnid) istat=NF_INQ_VARID(ncdfid,'wcside',wcsideid) istat=NF_INQ_VARID(ncdfid,'track',trackid) ! get data istat=NF_GET_VARA_INT(ncdfid,profid,1,lenprofs,profnum) istat=NF_GET_VARA_INT(ncdfid,dateid,1,lenprofs,date) istat=NF_GET_VARA_REAL(ncdfid,timeid,1,lenprofs,time) istat=NF_GET_VARA_REAL(ncdfid,lonid,1,lenprofs,longitude) istat=NF_GET_VARA_REAL(ncdfid,latid,1,lenprofs,latitude) istat=NF_GET_VARA_INT(ncdfid,tpadid,1,lenprofs,tpad) istat=NF_GET_VARA_INT(ncdfid,tpdnid,1,lenprofs,tpdn) istat=NF_GET_VARA_INT(ncdfid,wcsideid,1,lenprofs,wcside) istat=NF_GET_VARA_REAL(ncdfid,trackid,1,lenprofs,track) ios=NF_CLOSE(ncdfid) do k=lenprofs+1,mxsatloc profnum(k)=-999 enddo nsat = lenprofs do k=1,nsat ! adflag=-tpad(k)+1 adflag=tpad(k) lst(k)=time(k)/60.+longitude(k)/15. if(lst(k) .lt. 0) lst(k)=lst(k)+24. if(lst(k) .ge. 24) lst(k)=lst(k)-24. time(k) = (mod(date(k),1000)-1)*24.+time(k)/60. if(longitude(k) .gt. 180.) longitude(k) = longitude(k) - 360. if (time(k)*60. .ge. mins0 .and. time(k)*60 .lt. mins1) then satlocs(k,1) = time(k)*60. satlocs(k,2) = longitude(k) satlocs(k,3) = latitude(k) satlocs(k,4) = date(k) satlocs(k,5) = adflag satlocs(k,6) = wcside(k) ! tidi warm(1)/cold(0) satlocs(k,7) = tpdn(k) satlocs(k,8) = lst(k) satlocs(k,9) = track(k) print *, ' <<< READ_TIDI_LOCS_NC: use point = >>> ', + k,satlocs(k,:) endif enddo ! k return end subroutine read_tidi_locs_nc subroutine read_saber_locs_nc(satlocs,nsat,mins0,mins1) ! internal to mksatut use input,only: satfile integer,intent(in) :: mins0,mins1 real,intent(out) :: satlocs(mxsatloc,9) ! satellite lat/lon/times integer,intent(out) :: nsat ! number of orbit points integer, external :: nextlu real :: orbnum, adflag character (len=1) :: adflag_ch integer :: lu,ios,k character(len=150) :: line character(len=1) :: char1 ! ! netcdf variables integer :: ncdfid, dimprofs, lenprofs,istat integer :: profid,dateid,timeid,lonid,latid,tpadid,tpdnid integer,dimension(mxsatloc) :: profnum,date,tpad,tpdn real,dimension(mxsatloc) :: longitude,latitude,time,lst lu=nextlu() ! open netcdf ios = NF_OPEN(satfile, NF_NOWRITE, ncdfid) if (ios /= NF_NOERR) then write(6,"('>>> WARNING read_saber_locs_nc: error opening file ', + a)") satfile stop 'read_saber_locs_nc' endif print *, ' <<< READ_SABER_LOCS_NC >>> : opened satfile: ', + ncdfid, ' ', satfile ! get dim IDs and dim length istat=NF_INQ_DIMID(ncdfid,'profs',dimprofs) istat=NF_INQ_DIMLEN(ncdfid,dimprofs,lenprofs) ! get var IDs istat=NF_INQ_VARID(ncdfid,'profnum',profid) istat=NF_INQ_VARID(ncdfid,'date',dateid) istat=NF_INQ_VARID(ncdfid,'time',timeid) istat=NF_INQ_VARID(ncdfid,'lon',lonid) istat=NF_INQ_VARID(ncdfid,'lat',latid) istat=NF_INQ_VARID(ncdfid,'tpad',tpadid) istat=NF_INQ_VARID(ncdfid,'tpdn',tpdnid) ! get data istat=NF_GET_VARA_INT(ncdfid,profid,1,lenprofs,profnum) istat=NF_GET_VARA_INT(ncdfid,dateid,1,lenprofs,date) istat=NF_GET_VARA_REAL(ncdfid,timeid,1,lenprofs,time) istat=NF_GET_VARA_REAL(ncdfid,lonid,1,lenprofs,longitude) istat=NF_GET_VARA_REAL(ncdfid,latid,1,lenprofs,latitude) istat=NF_GET_VARA_INT(ncdfid,tpadid,1,lenprofs,tpad) istat=NF_GET_VARA_INT(ncdfid,tpdnid,1,lenprofs,tpdn) ios=NF_CLOSE(ncdfid) lenprofs=lenprofs-1 do k=lenprofs+1,mxsatloc profnum(k)=-999 enddo nsat = lenprofs do k=1,nsat ! adflag=-tpad(k)+1 adflag=tpad(k) lst(k)=time(k)/60.+longitude(k)/15. if(lst(k) .lt. 0) lst(k)=lst(k)+24. if(lst(k) .ge. 24) lst(k)=lst(k)-24. time(k) = (mod(date(k),1000)-1)*24.+time(k)/60. if(longitude(k) .gt. 180.) longitude(k) = longitude(k) - 360. if (time(k)*60. .ge. mins0 .and. time(k)*60 .lt. mins1) then satlocs(k,1) = time(k)*60. satlocs(k,2) = longitude(k) satlocs(k,3) = latitude(k) satlocs(k,4) = date(k) satlocs(k,5) = adflag satlocs(k,6) = -999. ! no saber angles satlocs(k,7) = tpdn(k) satlocs(k,8) = lst(k) satlocs(k,9) = -999. ! no saber track angle print *, ' <<< READ_SABER_LOCS_NC: use point = >>> ', + k,satlocs(k,:),nsat endif enddo ! k return end subroutine read_saber_locs_nc !----------------------------------------------------------------------- subroutine read_champ_locs_nc(satlocs,nsat,mins0,mins1) ! internal to mksatut use input,only: satfile integer,intent(in) :: mins0,mins1 real,intent(out) :: satlocs(mxsatloc,9) ! satellite lat/lon/times integer,intent(out) :: nsat ! number of orbit points integer, external :: nextlu real :: orbnum, adflag character (len=1) :: adflag_ch integer :: lu,ios,k character(len=150) :: line character(len=1) :: char1 ! ! netcdf variables integer :: ncdfid, dimprofs, lenprofs,istat integer :: profid,dateid,timeid,lonid,latid,tpadid,tpdnid integer,dimension(mxsatloc) :: profnum,date,tpad,tpdn real,dimension(mxsatloc) :: longitude,latitude,time,lst lu=nextlu() ! open netcdf ios = NF_OPEN(satfile, NF_NOWRITE, ncdfid) if (ios /= NF_NOERR) then write(6,"('>>> WARNING read_champs_locs_nc: error opening file ', + a)") satfile stop 'read_champ_locs_nc' endif print *, ' <<< READ_CHAMP_LOCS_NC >>> : opened satfile: ', + ncdfid, ' ', satfile ! get dim IDs and dim length istat=NF_INQ_DIMID(ncdfid,'profs',dimprofs) istat=NF_INQ_DIMLEN(ncdfid,dimprofs,lenprofs) ! get var IDs istat=NF_INQ_VARID(ncdfid,'profnum',profid) istat=NF_INQ_VARID(ncdfid,'date',dateid) istat=NF_INQ_VARID(ncdfid,'time',timeid) istat=NF_INQ_VARID(ncdfid,'lon',lonid) istat=NF_INQ_VARID(ncdfid,'lat',latid) istat=NF_INQ_VARID(ncdfid,'tpad',tpadid) istat=NF_INQ_VARID(ncdfid,'tpdn',tpdnid) ! get data istat=NF_GET_VARA_INT(ncdfid,profid,1,lenprofs,profnum) istat=NF_GET_VARA_INT(ncdfid,dateid,1,lenprofs,date) istat=NF_GET_VARA_REAL(ncdfid,timeid,1,lenprofs,time) istat=NF_GET_VARA_REAL(ncdfid,lonid,1,lenprofs,longitude) istat=NF_GET_VARA_REAL(ncdfid,latid,1,lenprofs,latitude) istat=NF_GET_VARA_INT(ncdfid,tpadid,1,lenprofs,tpad) istat=NF_GET_VARA_INT(ncdfid,tpdnid,1,lenprofs,tpdn) ios=NF_CLOSE(ncdfid) lenprofs=lenprofs-1 do k=lenprofs+1,mxsatloc profnum(k)=-999 enddo nsat = lenprofs do k=1,nsat ! adflag=-tpad(k)+1 adflag=tpad(k) lst(k)=time(k)/60.+longitude(k)/15. if(lst(k) .lt. 0) lst(k)=lst(k)+24. if(lst(k) .ge. 24) lst(k)=lst(k)-24. time(k) = (mod(date(k),1000)-1)*24.+time(k)/60. if(longitude(k) .gt. 180.) longitude(k) = longitude(k) - 360. if (time(k)*60. .ge. mins0 .and. time(k)*60 .lt. mins1) then satlocs(k,1) = time(k)*60. satlocs(k,2) = longitude(k) satlocs(k,3) = latitude(k) satlocs(k,4) = date(k) satlocs(k,5) = adflag satlocs(k,6) = -999. ! no champ angles satlocs(k,7) = tpdn(k) satlocs(k,8) = lst(k) satlocs(k,9) = -999. ! no champ track angle print *, ' <<< READ_CHAMP_LOCS_NC: use point = >>> ', + k,satlocs(k,:),nsat endif enddo ! k return end subroutine read_champ_locs_nc !----------------------------------------------------------------------- subroutine read_guvi_locs_nc(satlocs,nsat,mins0,mins1) ! internal to mksatut use input,only: satfile integer,intent(in) :: mins0,mins1 real,intent(out) :: satlocs(mxsatloc,9) ! satellite lat/lon/times integer,intent(out) :: nsat ! number of orbit points integer, external :: nextlu real :: orbnum, adflag character (len=1) :: adflag_ch integer :: lu,ios,k character(len=150) :: line character(len=1) :: char1 ! ! netcdf variables integer :: ncdfid, dimprofs, lenprofs,istat integer :: profid,dateid,timeid,lonid,latid,tpadid,tpdnid integer,dimension(mxsatloc) :: profnum,date,tpad,tpdn real,dimension(mxsatloc) :: longitude,latitude,time,lst lu=nextlu() ! open netcdf ios = NF_OPEN(satfile, NF_NOWRITE, ncdfid) if (ios /= NF_NOERR) then write(6,"('>>> WARNING read_guvi_locs_nc: error opening file ', + a)") satfile stop 'read_guvi_locs_nc' endif print *, ' <<< READ_GUVI_LOCS_NC >>> : opened satfile: ', + ncdfid, ' ', satfile ! get dim IDs and dim length istat=NF_INQ_DIMID(ncdfid,'profs',dimprofs) istat=NF_INQ_DIMLEN(ncdfid,dimprofs,lenprofs) ! get var IDs istat=NF_INQ_VARID(ncdfid,'profnum',profid) istat=NF_INQ_VARID(ncdfid,'date',dateid) istat=NF_INQ_VARID(ncdfid,'time',timeid) istat=NF_INQ_VARID(ncdfid,'lon',lonid) istat=NF_INQ_VARID(ncdfid,'lat',latid) istat=NF_INQ_VARID(ncdfid,'tpad',tpadid) istat=NF_INQ_VARID(ncdfid,'tpdn',tpdnid) ! get data istat=NF_GET_VARA_INT(ncdfid,profid,1,lenprofs,profnum) istat=NF_GET_VARA_INT(ncdfid,dateid,1,lenprofs,date) istat=NF_GET_VARA_REAL(ncdfid,timeid,1,lenprofs,time) istat=NF_GET_VARA_REAL(ncdfid,lonid,1,lenprofs,longitude) istat=NF_GET_VARA_REAL(ncdfid,latid,1,lenprofs,latitude) istat=NF_GET_VARA_INT(ncdfid,tpadid,1,lenprofs,tpad) istat=NF_GET_VARA_INT(ncdfid,tpdnid,1,lenprofs,tpdn) ios=NF_CLOSE(ncdfid) lenprofs=lenprofs-1 do k=lenprofs+1,mxsatloc profnum(k)=-999 enddo nsat = lenprofs do k=1,nsat ! adflag=-tpad(k)+1 adflag=tpad(k) lst(k)=time(k)/60.+longitude(k)/15. if(lst(k) .lt. 0) lst(k)=lst(k)+24. if(lst(k) .ge. 24) lst(k)=lst(k)-24. time(k) = (mod(date(k),1000)-1)*24.+time(k)/60. if(longitude(k) .gt. 180.) longitude(k) = longitude(k) - 360. if (time(k)*60. .ge. mins0 .and. time(k)*60 .lt. mins1) then satlocs(k,1) = time(k)*60. satlocs(k,2) = longitude(k) satlocs(k,3) = latitude(k) satlocs(k,4) = date(k) satlocs(k,5) = adflag satlocs(k,6) = -999. ! no guvi angles satlocs(k,7) = tpdn(k) satlocs(k,8) = lst(k) satlocs(k,9) = -999. ! no guvi track angle print *, ' <<< READ_GUVI_LOCS_NC: use point = >>> ', + k,satlocs(k,:),nsat endif enddo ! k return end subroutine read_guvi_locs_nc !----------------------------------------------------------------------- subroutine read_cedar_locs_nc(satlocs,nsat,mins0,mins1) ! internal to mksatut use input,only: satfile integer,intent(in) :: mins0,mins1 real,intent(out) :: satlocs(mxsatloc,9) ! satellite lat/lon/times integer,intent(out) :: nsat ! number of orbit points integer, external :: nextlu real :: orbnum, adflag character (len=1) :: adflag_ch integer :: lu,ios,k character(len=150) :: line character(len=1) :: char1 ! ! netcdf variables integer :: ncdfid, dimprofs, lenprofs,istat integer :: profid,dateid,timeid,lonid,latid,tpadid,tpdnid integer,dimension(mxsatloc) :: profnum,date,tpad,tpdn real,dimension(mxsatloc) :: longitude,latitude,time,lst lu=nextlu() ! open netcdf ios = NF_OPEN(satfile, NF_NOWRITE, ncdfid) if (ios /= NF_NOERR) then write(6,"('>>> WARNING read_cedar_locs_nc: error opening file ', + a)") satfile stop 'read_cedar_locs_nc' endif print *, ' <<< READ_CEDAR_LOCS_NC >>> : opened satfile: ', + ncdfid, ' ', satfile ! get dim IDs and dim length istat=NF_INQ_DIMID(ncdfid,'profs',dimprofs) istat=NF_INQ_DIMLEN(ncdfid,dimprofs,lenprofs) ! get var IDs istat=NF_INQ_VARID(ncdfid,'profnum',profid) istat=NF_INQ_VARID(ncdfid,'date',dateid) istat=NF_INQ_VARID(ncdfid,'time',timeid) istat=NF_INQ_VARID(ncdfid,'lon',lonid) istat=NF_INQ_VARID(ncdfid,'lat',latid) istat=NF_INQ_VARID(ncdfid,'tpad',tpadid) istat=NF_INQ_VARID(ncdfid,'tpdn',tpdnid) ! get data istat=NF_GET_VARA_INT(ncdfid,profid,1,lenprofs,profnum) istat=NF_GET_VARA_INT(ncdfid,dateid,1,lenprofs,date) istat=NF_GET_VARA_REAL(ncdfid,timeid,1,lenprofs,time) istat=NF_GET_VARA_REAL(ncdfid,lonid,1,lenprofs,longitude) istat=NF_GET_VARA_REAL(ncdfid,latid,1,lenprofs,latitude) istat=NF_GET_VARA_INT(ncdfid,tpadid,1,lenprofs,tpad) istat=NF_GET_VARA_INT(ncdfid,tpdnid,1,lenprofs,tpdn) ios=NF_CLOSE(ncdfid) lenprofs=lenprofs-1 do k=lenprofs+1,mxsatloc profnum(k)=-999 enddo nsat = lenprofs do k=1,nsat ! adflag=-tpad(k)+1 adflag=tpad(k) lst(k)=time(k)/60.+longitude(k)/15. if(lst(k) .lt. 0) lst(k)=lst(k)+24. if(lst(k) .ge. 24) lst(k)=lst(k)-24. time(k) = (mod(date(k),1000)-1)*24.+time(k)/60. if(longitude(k) .gt. 180.) longitude(k) = longitude(k) - 360. if (time(k)*60. .ge. mins0 .and. time(k)*60 .lt. mins1) then satlocs(k,1) = time(k)*60. satlocs(k,2) = longitude(k) satlocs(k,3) = latitude(k) satlocs(k,4) = date(k) satlocs(k,5) = adflag satlocs(k,6) = -999. ! no cedar angles satlocs(k,7) = tpdn(k) satlocs(k,8) = lst(k) satlocs(k,9) = -999. ! no cedar track anlge print *, ' <<< READ_CEDAR_LOCS_NC: use point = >>> ', + k,satlocs(k,:),nsat endif enddo ! k return end subroutine read_cedar_locs_nc !------------------------------------------------------------------- subroutine interp_sat(f,nf,fprev,time0,time1,nlev,loc_sat, | data_int) integer,intent(in) :: nf,nlev type(field),intent(in) :: f(nf) type(field),intent(in) :: fprev(nf) real, intent(in) :: loc_sat(9) real, intent(out) :: data_int(nlev,nf) integer,intent(in) :: time0,time1 ! real :: frac_t,frac_lon,frac_lat real :: dlat1t0,dlat1t1,dlat2t0,dlat2t1,dt0,dt1 integer :: i,j,ixlon,ixlat integer,external :: ixfind,ixfindc ! ixlon = ixfind(gcmlon,nlon,loc_sat(2),dlon) if (gcmlon(ixlon) > loc_sat(2)) ixlon = ixlon-1 ixlat = ixfind(gcmlat,nlat,loc_sat(3),dlat) if (gcmlat(ixlat) > loc_sat(3)) ixlat = ixlat-1 frac_t = (loc_sat(1)-float(time0)) / float(time1-time0) frac_lon = (loc_sat(2)-gcmlon(ixlon)) / + (gcmlon(ixlon+1)-gcmlon(ixlon)) frac_lat = (loc_sat(3)-gcmlat(ixlat)) / + (gcmlat(ixlat+1)-gcmlat(ixlat)) if(frac_lon > 1.0 .or. frac_lon < 0.) + print *, ' << INTERP_SAT >> : frac_lon out of range ', + gcmlon(ixlon), loc_sat(2), gcmlon(ixlon+1) ! !# print *, ' << INTERP_SAT >> time0,time1,loc_sat: ', !# + time0, time1, loc_sat !# print *, ' << INTERP_SAT >> fracs: ', !# + frac_t, frac_lon, frac_lat do j=1,nf do i=1,nlev data_int(i,j) = spval if(j <= nf .and. i <= npress) then dlat1t0 = (1.-frac_lon) * + fprev(j)%data(ixlon ,ixlat ,i) + + frac_lon * + fprev(j)%data(ixlon+1,ixlat ,i) dlat1t1 = (1.-frac_lon) * + f(j)%data(ixlon ,ixlat ,i) + + frac_lon * + f(j)%data(ixlon+1,ixlat ,i) dlat2t0 = (1.-frac_lon) * + fprev(j)%data(ixlon ,ixlat+1,i) + + frac_lon * + fprev(j)%data(ixlon+1,ixlat+1,i) dlat2t1 = (1.-frac_lon) * + f(j)%data(ixlon ,ixlat+1,i) + + frac_lon * + f(j)%data(ixlon+1,ixlat+1,i) dt0 = (1.-frac_lat) * dlat1t0 + frac_lat * dlat2t0 dt1 = (1.-frac_lat) * dlat1t1 + frac_lat * dlat2t1 data_int(i,j) = (1.-frac_t) * dt0 + frac_t * dt1 !# !# dlat2t0 = fprev(j)%data(ixlon+1,ixlat+1,i) !# dlat2t1 = f(j)%data(ixlon+1,ixlat+1,i) !# dt0 = dlat2t0 !# dt1 = dlat2t1 !# data_int(i,j) = frac_t * dt0 + (1-frac_t) * dt1 !# if (i .eq. 5 .and. j .eq. 5) then !# print *, ' << INTERP_SAT >> lon interp: ', !# + fprev(j)%data(ixlon,ixlat,i), !# + dlat1t0, !# + fprev(j)%data(ixlon+1,ixlat,i), !# + f(j)%data(ixlon,ixlat,i), !# + dlat1t1, !# + f(j)%data(ixlon+1,ixlat,i), !# + fprev(j)%data(ixlon,ixlat+1,i), !# + dlat2t1, !# + fprev(j)%data(ixlon+1,ixlat+1,i), !# + f(j)%data(ixlon,ixlat+1,i), !# + dlat2t1, !# + f(j)%data(ixlon+1,ixlat+1,i) !# print *, ' << INTERP_SAT >> lat interp: ', !# + dlat1t0, dt0, dlat2t0, !# + dlat1t1, dt1, dlat2t1 !# print *, ' << INTERP_SAT >> time interp: ', !# + dt0,data_int(i,j),dt1 !# endif !# data_int(i,j) = f(j)%data(ixlon,ixlat,i) endif enddo ! j enddo ! i return end subroutine interp_sat !------------------------------------------------------------------- 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(satut),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 ! ! 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 pltsatut(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(satut) utv has been fully defined. ! ! Args: type(satut),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)=='MID') then call yaxzpht(gcmlev,hts,gcmlev,utv%nx,npress,dlev, + utv%yy,utv%ny,p0,utv%ihtyax,iyaxright,vp,spval) else call yaxzpht(gcmilev,hts,gcmilev,utv%nx,npress,dlev, + utv%yy,utv%ny,p0,utv%ihtyax,iyaxright,vp,spval) endif endif return end subroutine pltsatut !------------------------------------------------------------------- subroutine wrout_satut(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(satut),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) 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_satut end module mk_satut