! module mk_lons use hist,only: history implicit none ! ! A longitude slice is a 2d contour at a selected longitude, ! with latitude on the x-axis, and height or pressure on ! the y-axis. Lon slice plots can have optional extra ! right axis in ht or zp. ! User may specify selected longitudes, local times, zonal means, ! amplitudes and phases, or stream function. Lontype is set ! 'lon','slt','zm ','amp','pha', or 'str' accordingly. ! ! Longitude slice structure: type lonslice real :: glon ! geog longitude of slice real :: slt ! local time of slice character(len=3) :: lontype integer :: nx,ny ! number of x,y points real,pointer :: xx(:)=>NULL() ! x coords real,pointer :: yy(:)=>NULL() ! y coords real,pointer :: y(:,:)=>NULL() ! epflux vector (s.a. vectype) real,pointer :: z(:,:)=>NULL() ! epflux vector (s.a. vectype) character(len=40) :: xlab,ylab ! x and y axis labels real,pointer :: data(:,:)=>NULL() ! the data array character(len=80) :: fname ! full field name character(len=8) :: sname ! short field name character(len=32) :: vectype ! 'EPVYZ' or '+EPVYZ' ep vectors character(len=16) :: funits ! field units character(len=16) :: zptype ! zptype (MIDPOINTS or INTERFACES) character(len=16) :: ftype ! field type character(len=8) :: difftype ! raw or percent logical :: known ! true if field is known to the processor integer :: ihtyax ! ht on yaxis if ihtyax==1, zp otherwise integer :: irightax ! flag for extra right axis (ht or zp) integer :: log10 ! if > 0, then log10 of field is plotted character(len=240) :: tlabs(3),blabs(3) ! top and bottom labels integer :: nwave ! wave number if lontype='amp' or 'pha' end type lonslice ! real :: tlabchsz(3) = (/.022,.022,.022/), + blabchsz(3) = (/.020,.020,.020/) ! ! Timeav(nslon,nflds,0:1): array of lonslice structures for time-averaged data. ! (last dimension is for zp and ht on y axes) ! One structure for each selected lon/slt/zm and each field. ! Same for control case for differences. ! type(lonslice),allocatable :: timeav(:,:,:) type(lonslice),allocatable :: timeav_prev(:,:,:) ! for previous time type(lonslice),allocatable :: timeav_cntr(:,:,:) ! control case for diffs type(lonslice),allocatable :: timeav_cntr_prev(:,:,:) ! ! Local filename of lon slice netcdf file: character(len=16) :: flnm_cdf_lons = 'lonslices.nc' logical :: istimeave = .false. contains !------------------------------------------------------------------- subroutine mklons(f,fcntr,nf,h,hcntr,itime) ! ! Make longitude slice plots (latitude on x-axis, ht/zp on y-axis) at ! selected longitudes, local times, or zonal means. ! use proc use fields,only: field,oh_alt,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use input use plt use set_vert ! ! Args: integer,intent(in) :: nf,itime type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: i,ii,ixf,ilon,ier,ixlon,izpht0,izpht1,iht,nhtscale,ixz, + iadvfr,nlons,nslts,k0zp,k1zp,nzprange,ixlat0,ixlat1,j,k,ix, + nstrm,logint,jj,kk,nlats,ixepvy,ixepvz,nslon,it,nspv type(lonslice) :: lonsli character(len=80) :: msgout real,pointer :: htscale(:)=>NULL(),zprange(:)=>NULL() ! real :: vp(4),vplon(4)=(/.10,.87,.23,.88/) ! for 0 extra right axes real :: vp(4),vplon(4)=(/.10,.87,.26,.88/) ! for 0 extra right axes real :: toffset=.06, boffset=.15 real :: slons(mxslice),sslts(mxslice),dum,slt real :: zsli(nlat,npress),zsli_cntr(nlat,npress) real :: fmin,fmax,fislice(nlon),fislice_cntr(nlon) real,allocatable :: tmpsli1(:,:),tmpsli2(:,:),tmpsli3(:,:) ! ! hscale and fhscale are defined by sethscale when field is height-only: ! fhscale(nfhscale) will be height scale of height-only fields. ! hscale(nhscale) will be height scale to use when plotting. ! real,pointer :: hscale(:)=>NULL() ! height scale to be used integer :: nhscale ! dimension for hscale real,pointer :: fhscale(:)=>NULL()! height scale of ht-only field integer :: nfhscale ! dimension for fhscale ! ! Externals: integer,external :: ixfind,ixfindc real,external :: fslt,fmean logical,external :: isslt ! ! Begin exec: write(6,"(' ')") if (.not.diffs) then write(6,"('Longitude slices at ut = ',f5.2,' mtime = ',3i3)") + h%ut,h%mtime else write(6,"('Difference fields: ')",advance='NO') write(6,"('Longitude slices at ut = ',f5.2)") h%ut write(6,"(' perturbed mtime = ',3i4,' control mtime = ',3i4)") + h%mtime,hcntr%mtime endif if (multiadvfr > 0) nppf = 0 ! ! Determine number of selected lons and/or slts, and place in local ! arrays slons and sslts: ! (flons (selected lons) and fslts (selected slts) were provided by user) ! Also, convert any flons='LTxx' to fslts='xx' options (i.e. store in sslts) nlons = 0 nslts = 0 do i=1,mxslice if (flons(i) /= spval) then if (isslt(flons(i),slt)) then nslts = nslts+1 sslts(nslts) = slt else nlons = nlons+1 slons(nlons) = flons(i) endif endif if (fslts(i) /= spval) then nslts = nslts+1 sslts(nslts) = fslts(i) endif enddo if (nlons == 0 .and. nslts == 0) then write(6,"('mklons: NO plots made (no flons or fslts)')") return else ! write(6,"('Will make lon slices at ',i2,' longitudes ', ! + 'and ',i2,' local times.')") nlons,nslts endif ! ! Set up time-averaging structures (1st time only): nstrm = 0 if (istream > 0.and.any(f(:)%fname8=='VN ')) nstrm = 1 nslon = nlons+nslts+nstrm ! write(6,"('nlons=',i3,' nslts=',i3,' nstrm=',i3,' nslon=',i3)") ! | nlons,nslts,nstrm,nslon if (ilon_timeave > 0 .and. itime==1) then allocate(timeav(nslon,nf,0:1)) write(6,"('mklons: allocated timeav(nslon,nflds,0:1), where', | ' nslon=',i3,' nflds=',i3)") nslon,nf allocate(timeav_prev(nslon,nf,0:1)) write(6,"('mklons: allocated timeav_prev(nslon,nflds,0:1), ', | 'where nslon=',i3,' nflds=',i3)") nslon,nf ! ! Allocate control structures for diffs: if (diffs) then allocate(timeav_cntr(nslon,nf,0:1)) write(6,"('mklons: allocated timeav_cntr(nslon,nflds,0:1), ', | 'where nslon=',i3,' nflds=',i3)") nslon,nf allocate(timeav_cntr_prev(nslon,nf,0:1)) write(6,"('mklons: allocated timeav_cntr_prev', | '(nslon,nflds,0:1) where nslon=',i3,' nflds=',i3)") nslon,nf endif endif ! ! Define x coords in latitude: ! (User may optionally restrict latitude range with flon_xlatrange) ! (,flon_xlatrange is from user, ixlat0,ixlat1 are returned) ! subroutine setxlat(lonsli,xlatrange,ixlat0,ixlat1) ! call setxlat(lonsli,flon_xlatrange,ixlat0,ixlat1) nlats = ixlat1-ixlat0+1 if (ilon_timeave > 0 .and. itime==1) then do i=1,nf do ii=1,nslon call setxlat(timeav(ii,i,0),flon_xlatrange,ixlat0,ixlat1) call setxlat(timeav(ii,i,1),flon_xlatrange,ixlat0,ixlat1) enddo enddo endif ! ! Locate index to heights field (needed for interp, etc): ! ixz = ixfindc(f%fname8,nf,'Z ') if (ixz <= 0) then write(6,"('>>> mklons: need heights in fields: field names', + ' = ',(8a8))") f%fname8 stop 'mklons z' endif ! ! Field loop. Field must be requested, data pointer must be associated, ! and field must not be height-independent: ! (if stream function is requested, VN is a dependency which may or ! may not have been requested) ! fields_loop: do ixf=1,nf nstrm = 0 if (istream > 0.and.trim(f(ixf)%fname8)=='VN') then nstrm = 1 elseif (.not.f(ixf)%requested.or..not.associated(f(ixf)%data) + .or.f(ixf)%nlev==1) then cycle fields_loop endif ! ! Set up vertical scale(s) zp and/or ht: ! 3/16/06 btf: moved this to inside fields loop so it is called ! for zptype of current field. ! if (f(ixf)%zptype(1:3)=='INT') then ! interfaces call setvert(flon_zprange,gcmilev,npress,dlev,nzprange, | k0zp,k1zp,zprange,flon_htscale,nhtscale,htscale,spval,1) else ! midpoints call setvert(flon_zprange,gcmlev,npress,dlev,nzprange, | k0zp,k1zp,zprange,flon_htscale,nhtscale,htscale,spval,1) endif ! ! Set loop indices for zp and/or ht on y-axis: izpht0 = 0 izpht1 = 1 if (nhtscale == 0) izpht1 = 0 if (nzprange == 0) izpht0 = 1 ! ! Set some lonslice components from field components: lonsli%fname = f(ixf)%fname56 lonsli%sname = f(ixf)%fname8 lonsli%zptype = f(ixf)%zptype lonsli%known = f(ixf)%known lonsli%funits = f(ixf)%units lonsli%difftype = f(ixf)%difftype if (lonsli%sname=='EPVYZ '.and.vep_scale(3)<0.) + lonsli%funits = 'NORMALIZED ' lonsli%ftype = f(ixf)%type if (ilon_timeave > 0) then timeav(:,ixf,:)%fname = lonsli%fname timeav(:,ixf,:)%sname = lonsli%sname timeav(:,ixf,:)%funits = lonsli%funits timeav(:,ixf,:)%ftype = lonsli%ftype timeav(:,ixf,:)%zptype = lonsli%zptype timeav(:,ixf,:)%known = lonsli%known timeav(:,ixf,:)%difftype = lonsli%difftype timeav(:,ixf,:)%xlab = lonsli%xlab timeav(:,ixf,:)%ylab = lonsli%ylab endif logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 ! ! 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) ! ! Calculate stream function from VN: ! 5/8/98: stream function calculated on VN only (W is implicit in V): ! ! nstrm = 0 ! if (istream>0.and.trim(f(ixf)%fname8)=='VN') nstrm = 1 ! ! Initialize EPVY+EPVZ vectors: ! EPVYZ -> plot vector arrows of EPVY+EPVZ (lonsli%data, y, and z are ! set by mkepvec (lonsli%data used only to report min,max) ! EPVYZMAG -> contour EPVY+EVPVZ magnitudes (lonsli%data set by mkepvec) ! ilon_epvdiv_yz > 0 -> add EPVYZ vectors to EPVDIV contours ! ilon_epvyzmag_yz > 0 -> add EPVYZ vectors to EPVYZMAG contours ! ! Lonsli%data, when EPVYZ or EPVYZMAG must be recalculated by mkepvec ! so we get vector sum of means of y and z rather than means of the ! vector sums. ! lonsli%vectype = ' ' if (ilon_timeave > 0) then do i=1,nslon timeav(:,ixf,0)%vectype = ' ' timeav(:,ixf,1)%vectype = ' ' enddo endif if (f(ixf)%fname8=='EPVYZ '.or.f(ixf)%fname8=='EPVYZMAG'.or. + (f(ixf)%fname8=='EPVDIV '.and.ilon_epvdiv_yz > 0)) then ixepvy = ixfindc(f%fname8,nf,'EPVY ') ixepvz = ixfindc(f%fname8,nf,'EPVZ ') if (ixepvy<=0.or.ixepvz<=0) then write(6,"('>>> mklons: need EPVY and EPVZ to ', + 'make epflux vectors: vectors will not be plotted.')") else ! internal (5th arg==1 means init) call mkepvec(f(ixepvy),f(ixepvz),1,0,vepz_viewfac, + vep_scale(3)) endif endif if (len_trim(lonsli%vectype)>0) then blabchsz = .018 else blabchsz = .020 endif ! ! Selected longitude loop (lons are at nearest grid point): lon_loop: do ilon=1,nlons+nslts+nstrm if (nstrm>0.and.ilon <= nlons+nslts.and. + trim(f(ixf)%fname8)=='VN'.and..not.f(ixf)%requested) + cycle lon_loop lonsli%lontype = ' ' if (ilon <= nlons) then ! lon or zonal means if (slons(ilon) /= zmflag) then ! not zonal means ixlon = ixfind(gcmlon,nlon,slons(ilon),dlon) lonsli%glon = gcmlon(ixlon) lonsli%slt = fslt(dum,h%ut,lonsli%glon,1) lonsli%lontype = 'lon' else ! zonal means lonsli%lontype = 'zm ' lonsli%glon = zmflag endif elseif (ilon > nlons+nslts) then ! stream function lonsli%lontype = 'str' else ! local time lonsli%glon = fslt(sslts(ilon-nlons),h%ut,dum,3) ixlon = ixfind(gcmlon,nlon,lonsli%glon,dlon) lonsli%glon = gcmlon(ixlon) lonsli%slt = sslts(ilon-nlons) lonsli%lontype = 'slt' endif if (ilon_timeave > 0) then timeav(ilon,ixf,:)%lontype = lonsli%lontype timeav(ilon,ixf,:)%glon = lonsli%glon timeav(ilon,ixf,:)%slt = lonsli%slt endif ! ! Save heights in zsli: if (lonsli%lontype=='zm '.or.lonsli%lontype=='str') then do j=ixlat0,ixlat1 do k=1,npress if (trim(f(ixf)%zptype) == 'INTERFACES') then zsli(j,k) = fmean(z_ifaces(:,j,k),nlon,spval,0) if (diffs) then zsli_cntr(j,k) = fmean(zcntr_ifaces(:,j,k), | nlon,spval,0) endif else zsli(j,k) = fmean(z_midpts(:,j,k),nlon,spval,0) if (diffs) then zsli_cntr(j,k) = fmean(zcntr_midpts(:,j,k), | nlon,spval,0) endif endif enddo enddo else ! lon or slt if (trim(f(ixf)%zptype) == 'INTERFACES') then zsli(:,:) = z_ifaces(ixlon,:,:) if (diffs) zsli_cntr(:,:) = zcntr_ifaces(ixlon,:,:) else zsli(:,:) = z_midpts(ixlon,:,:) if (diffs) zsli_cntr(:,:) = zcntr_midpts(ixlon,:,:) endif endif ! ! Check for nans in zsli: ! subroutine check_nans(f,id1,id2,id3,name,n_total,ispval,spval, ! iprint,ifatal) ! Number of nans is returned in nnans_z. ! ! call check_nans(zsli(ixlat0:ixlat1,npress),ixlat1-ixlat0+1, ! | npress,1,'ZSLI',nnans_z,ispval,spval,1,0) ! ! Zp and/or ht on y-axis (0,0 for zp only, 1,1 for ht only, ! or 0,1 for zp and ht): yax_loop: do iht=izpht0,izpht1 lonsli%ihtyax = iht if (associated(lonsli%yy)) deallocate(lonsli%yy) if (associated(lonsli%data)) deallocate(lonsli%data) if (ilon_timeave > 0) timeav(ilon,ixf,iht)%ihtyax = iht ! ! - - - - - - - - - - - - zp on y-axis: - - - - - - - - - - - - - if (iht == 0) then ! zp on y-axis ! ! Adjust right viewport according to number of extra right-hand y-axes: if (ilon_yaxright==0) then vplon(2) = .90 elseif (ilon_yaxright==1) then vplon(2) = .85 else vplon(2) = .76 endif ! ! Skip fields that are plotted only in height (e.g. oh-v and oh-b): if (trim(f(ixf)%vtype)=='HEIGHT') cycle yax_loop ! ! Define y coords in log pressure (from setvert, called above): lonsli%ny = nzprange allocate(lonsli%yy(lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%yy") lonsli%yy = zprange ! array op write(lonsli%ylab,"('LN(P0/P) (',a,')')") | trim(lonsli%zptype) if (ilon_timeave > 0) then timeav(ilon,ixf,iht)%ny = nzprange allocate(timeav(ilon,ixf,iht)% | yy(timeav(ilon,ixf,iht)%ny),stat=ier) if (ier /= 0) call allocerr(ier, | "mklons allocating timeav(ilon,ixf,iht)%yy") timeav(ilon,ixf,iht)%yy = zprange ! array op write(timeav(ilon,ixf,iht)%ylab,"('LN(P0/P) (',a,')')") | trim(timeav(ilon,ixf,iht)%zptype) endif ! ! Allocate temporary lon slices: if (diffs) then ! allocate tmpsli1 and tmpsli2 allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocating tmpsli1 for perturbed case") allocate(tmpsli2(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocating tmpsli2 for control case") elseif (lonsli%lontype=='str') then allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocating tmpsli1 for stream function") endif ! ! Allocate and define lonsli data for zp on y-axis: allocate(lonsli%data(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%data") ! ! Allocate time-averaged data structures, first time only: if (ilon_timeave > 0.and.itime==1) then allocate(timeav(ilon,ixf,iht)%data(timeav(ilon,ixf, | iht)%nx,timeav(ilon,ixf,iht)%ny),stat=ier) if (ier /= 0) | call allocerr(ier,"mklons allocating timeav%data") ! write(6,"('mklons allocated timeav%data: ilon=',i3, ! | ' ixf=',i3,' iht=',i3)") ilon,ixf,iht ! ! Allocate time-averaged control case structures: if (diffs) then allocate(timeav_cntr(ilon,ixf,iht)%data | (timeav(ilon,ixf,iht)%nx,timeav(ilon,ixf,iht)%ny), | stat=ier) if (ier /= 0) call allocerr(ier, | "mklons allocating timeav_cntr%data") endif endif ! ! Zonal mean or stream: if (lonsli%lontype=='zm '.or.lonsli%lontype=='str') then do j=ixlat0,ixlat1 jj = j-ixlat0+1 do k=k0zp,k1zp kk = k-k0zp+1 ! ! Winterpark sgi crashes in fmean unless fislice is used in place of ! f(ixf)%data(:,j,k): ! lonsli%data(jj,kk) = fmean(f(ixf)%data(:,j,k), | nlon,spval,0) if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data(jj,kk)=lonsli%data(jj,kk) if (diffs) then tmpsli1(jj,kk) = fmean(f(ixf)%data(:,j,k), | nlon,spval,0) tmpsli2(jj,kk)= fmean(fcntr(ixf)%data(:,j,k), | nlon,spval,0) if (ilon_timeave > 0) timeav_cntr(ilon,ixf,iht)% | data(jj,kk)=tmpsli2(jj,kk) if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data(jj,kk) = | tmpsli2(jj,kk) endif enddo enddo ! ! Diffs of zm at zp: ! For time average differencing, we will take differences of averaged data ! (rather than average of diffs), so do not take diffs now. ! if (diffs) call mkdiffs(tmpsli1,tmpsli2,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) if (lonsli%lontype=='str') then if (.not.diffs) then call calcstrm(lonsli%data,lonsli%nx,lonsli%ny, + gcmlev(k0zp:k1zp),gcmlat(ixlat0:ixlat1),tmpsli1) lonsli%data = tmpsli1 if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data = lonsli%data else ! diffs of stream functions call calcstrm(tmpsli1,lonsli%nx,lonsli%ny, + gcmlev(k0zp:k1zp),gcmlat(ixlat0:ixlat1), + lonsli%data) call calcstrm(tmpsli2,lonsli%nx,lonsli%ny, + gcmlev(k0zp:k1zp),gcmlat(ixlat0:ixlat1),tmpsli1) tmpsli2 = lonsli%data call mkdiffs(tmpsli2,tmpsli1,lonsli%data, + lonsli%nx*lonsli%ny,'RAW') if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data = tmpsli1 endif endif else ! lon or slt if (.not.diffs) then lonsli%data= + f(ixf)%data(ixlon,ixlat0:ixlat1,k0zp:k1zp) if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data = lonsli%data else tmpsli1 = f(ixf)%data(ixlon,ixlat0:ixlat1,k0zp:k1zp) tmpsli2 = fcntr(ixf)%data(ixlon,ixlat0:ixlat1, + k0zp:k1zp) call mkdiffs(tmpsli1,tmpsli2,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data = tmpsli2 endif endif ! ! Get epflux vectors (epvy in lonsli%y, epvz in lonsli%z): ! (if we are doing EPVYZ or EPVYZMAG, lonsli%data will be reset ! by mkepvec so we get vector sum of means rather than means of ! vector sums) ! 7/9/01: time-averaging timeav: skipping this for now ! if (len_trim(lonsli%vectype) > 0.or. + f(ixf)%fname8=='EPVYZMAG') then if (.not.diffs) then call mkepvec(f(ixepvy),f(ixepvz),0,iht,vepz_viewfac, + vep_scale(3)) ! ! Diffs of epflux at zp: ! (tell mkepvec not to multiply by vepz_viewfac and not to normalize) else ! diffs call mkepvec(f(ixepvy),f(ixepvz),0,iht,1.,0.) if (f(ixf)%fname8/='EPVDIV ') then allocate(tmpsli3(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocate tmpsli3 for mkepvec diffs (iht=0)") endif ! ! Save perturbed from above: tmpsli1 = lonsli%y tmpsli2 = lonsli%z if (f(ixf)%fname8/='EPVDIV ') tmpsli3 = lonsli%data ! ! Get control into lonsli%y,z,data: ! (tell mkepvec not to multiply by vepz_viewfac and not to normalize) call mkepvec(fcntr(ixepvy),fcntr(ixepvz),0,iht,1.,0.) ! ! Take diffs of epvy,z,data at zp: call mkdiffs(tmpsli1,lonsli%y,lonsli%y, + lonsli%nx*lonsli%ny,f(ixf)%difftype) call mkdiffs(tmpsli2,lonsli%z,lonsli%z, + lonsli%nx*lonsli%ny,f(ixf)%difftype) call mkepvec(f(ixepvy),f(ixepvz),-1,iht,vepz_viewfac, + vep_scale(3)) ! scaling and normalization only if (f(ixf)%fname8/='EPVDIV ') + call mkdiffs(tmpsli3,lonsli%data,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) if (allocated(tmpsli3)) deallocate(tmpsli3) endif endif ! ! Release temporary lon slices: if (allocated(tmpsli1)) deallocate(tmpsli1) if (allocated(tmpsli2)) deallocate(tmpsli2) ! ! - - - - - - - - - - - - ht on y-axis: - - - - - - - - - - - - - else write(lonsli%ylab,"('HEIGHT (KM)')") if (ilon_timeave > 0) | write(timeav(ilon,ixf,iht)%ylab,"('HEIGHT (KM)')") if (ilon_yaxright==0) then vplon(2) = .90 else vplon(2) = .85 endif ! ! Allocate temporary lon slices: if (lonsli%lontype=='zm '.or.lonsli%lontype=='str') then allocate(tmpsli1(nlat,f(ixf)%nlev),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating tmpsli1 (iht=1)") allocate(tmpsli2(nlat,f(ixf)%nlev),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating tmpsli2 (iht=1)") endif if (diffs.and.lonsli%lontype=='str') then allocate(tmpsli3(nlat,f(ixf)%nlev),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating tmpsli3 (iht=1)") endif ! ! Store zonal means in tmpsli1: if (lonsli%lontype=='zm '.or.lonsli%lontype=='str') then do j=1,nlat do k=1,f(ixf)%nlev ! Sgi crashes here if fislice not used: ! fislice = f(ixf)%data(:,j,k) ! tmpsli1(j,k) = fmean(fislice,nlon,spval,0) tmpsli1(j,k) = fmean(f(ixf)%data(:,j,k), | nlon,spval,0) enddo enddo ! ! Store stream function from VN in tmpsli1: if (lonsli%lontype=='str') then call calcstrm(tmpsli1,nlat,f(ixf)%nlev,gcmlev,gcmlat, + tmpsli2) tmpsli1 = tmpsli2 endif ! ! Store control zonal means in tmpsli2: if (diffs) then do j=1,nlat do k=1,f(ixf)%nlev tmpsli2(j,k) = fmean(fcntr(ixf)%data(:,j,k), | nlon,spval,0) enddo enddo if (lonsli%lontype=='str') then call calcstrm(tmpsli2,nlat,f(ixf)%nlev,gcmlev, + gcmlat,tmpsli3) tmpsli2 = tmpsli3 endif endif endif ! ! Regular vtype field (vertical in pressure): if (trim(f(ixf)%vtype) /= 'HEIGHT') then ! ! Define y coords in height (from setvert, called above): lonsli%ny = nhtscale allocate(lonsli%yy(lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%yy") lonsli%yy = htscale ! array op ! ! Y-coords in height for time-averaging structure: if (ilon_timeave > 0) then timeav(ilon,ixf,iht)%ny = nhtscale allocate(timeav(ilon,ixf,iht)%yy(lonsli%ny),stat=ier) if (ier /= 0) | call allocerr(ier,"mklons allocating ', | timeav(ilon,ixf,iht)%yy") timeav(ilon,ixf,iht)%yy = htscale ! array op endif ! ! Allocate and define lonsli data for height on y-axis: allocate(lonsli%data(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%data") ! ! Allocate time-averaged structure data (first time only): if (ilon_timeave > 0 .and. itime==1) then allocate(timeav(ilon,ixf,iht)%data(timeav(ilon,ixf, | iht)%nx,timeav(ilon,ixf,iht)%ny),stat=ier) if (ier /= 0) | call allocerr(ier,"mklons allocating timeav%data") ! ! Allocate data for time-averaged control: if (diffs) then allocate(timeav_cntr(ilon,ixf,iht)%data | (timeav(ilon,ixf,iht)%nx,timeav(ilon,ixf,iht)%ny), | stat=ier) if (ier /= 0) call allocerr(ier, | "mklons allocating timeav_cntr%data") endif endif ! ! Need a 3rd temp slice for diffs: if (diffs) then if (allocated(tmpsli3)) deallocate(tmpsli3) allocate(tmpsli3(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating tmpsli3") endif ! ! If not doing zonal means, interpolate in height, defining lonsli%data. ! (if diffs, use pert tmpsli1 and cntr tmpsli2 defined above) ! If doing zonal means, use tmpsli and zsli calculated above. ! if (lonsli%lontype=='zm '.or.lonsli%lontype=='str') then ! ! Interpolate to htscale, store in lonsli%data (will be perturbed if diffs): call cuthtint(tmpsli1(ixlat0:ixlat1,:), + zsli(ixlat0:ixlat1,:),nlats,npress, + lonsli%data,htscale,nhtscale,logint,spval,ier,1) if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data = lonsli%data ! ! Difference fields of height-interpolated zm or stream function: if (diffs) then ! ! Save interpolated perturbed from above in tmpsli3: tmpsli3 = lonsli%data ! save interpolated pert if (allocated(tmpsli1)) deallocate(tmpsli1) allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating tmpsli1") ! ! Interpolate control to htscale, store in tmpsli1: call cuthtint(tmpsli2(ixlat0:ixlat1,:), + zsli_cntr(ixlat0:ixlat1,:),nlats,npress, + tmpsli1,htscale,nhtscale,logint,spval,ier,1) ! ! Take diffs: lonsli%data=tmpsli3-tmpsli1 call mkdiffs(tmpsli3,tmpsli1,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) ! ! 5-pt smoother in height: if (ismooth > 0) then ! input module (default == 1) do j=1,lonsli%nx ! scan across latitude x-axis call smoother(lonsli%data(j,:),lonsli%ny,5, | spval,0) enddo ! j=1,lonsli%nx endif if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data = tmpsli1 endif ! ! Make lon or slt slice (not zm or stream): else ! lon or slt (not zm or str) ! ! Interpolate to htscale, store in lonsli%data (is perturbed if diffs): call cuthtint(f(ixf)%data(ixlon,ixlat0:ixlat1,:), + zsli(ixlat0:ixlat1,:),nlats,npress,lonsli%data, + htscale,nhtscale,logint,spval,ier,1) if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data = lonsli%data ! ! Height-interpolated difference fields: if (diffs) then ! ! Allocate temporary lat vs ht slices: if (allocated(tmpsli1)) deallocate(tmpsli1) allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocating tmpsli1") if (allocated(tmpsli2)) deallocate(tmpsli2) allocate(tmpsli2(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, + "mklons allocating tmpsli2") ! ! Save interpolated perturbed from above in tmpsli1: tmpsli1 = lonsli%data ! save interpolated perturbed ! ! Interpolate control to htscale, store in tmpsli2: call cuthtint(fcntr(ixf)%data(ixlon,ixlat0:ixlat1,:) + ,zsli_cntr(ixlat0:ixlat1,:),nlats,npress,tmpsli2, + htscale,nhtscale,logint,spval,ier,1) ! ! Take difference fields lonsli%data = tmpsli1 - tmpsli2: call mkdiffs(tmpsli1,tmpsli2,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) ! ! 5-pt smoother in height: if (ismooth > 0) then ! input module (default == 1) do j=1,lonsli%nx ! smooth in height call smoother(lonsli%data(j,:),lonsli%ny,5, | spval,0) enddo ! j=1,lonsli%nx endif if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data = tmpsli2 ! endif ! ht-interpolated difference fields endif ! lon/slt or zm/str) ! ! Get epflux vectors (epvy in lonsli%y, epvz in lonsli%z): ! (if we are doing EPVYZ or EPVYZMAG, lonsli%data will be reset ! by mkepvec so we get vector sum of means rather than means of ! vector sums) ! 7/9/01: for timeav, skipping this for now. ! if (len_trim(lonsli%vectype) > 0.or. + f(ixf)%fname8=='EPVYZMAG') then if (.not.diffs) then call mkepvec(f(ixepvy),f(ixepvz),0,iht,vepz_viewfac, + vep_scale(3)) ! ! Diffs of epflux at ht: else if (allocated(tmpsli1)) deallocate(tmpsli1) allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier, + "mklons alloc tmpsli1 for mkepvec diffs (iht=1)") ! if (allocated(tmpsli2)) deallocate(tmpsli2) allocate(tmpsli2(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier, + "mklons alloc tmpsli2 for mkepvec diffs (iht=1)") ! if (f(ixf)%fname8/='EPVDIV ') then if (allocated(tmpsli3)) deallocate(tmpsli3) allocate(tmpsli3(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier, + "mklons alloc tmpsli3 for mkepvec diffs (iht=1)") endif ! ! Save perturbed from above: tmpsli1 = lonsli%y tmpsli2 = lonsli%z if (f(ixf)%fname8/='EPVDIV ') tmpsli3 = lonsli%data ! ! Get control into lonsli%y,z,data: ! (tell mkepvec not to multiply by vepz_viewfac and not to normalize) call mkepvec(fcntr(ixepvy),fcntr(ixepvz),0,iht, + 1.,0.) ! ! Take diffs of epvy,z,data at ht: call mkdiffs(tmpsli1,lonsli%y,lonsli%y, + lonsli%nx*lonsli%ny,f(ixf)%difftype) call mkdiffs(tmpsli2,lonsli%z,lonsli%z, + lonsli%nx*lonsli%ny,f(ixf)%difftype) call mkepvec(f(ixepvy),f(ixepvz),-1,iht, + vepz_viewfac,vep_scale(3)) ! scale and norm only if (f(ixf)%fname8/='EPVDIV ') + call mkdiffs(tmpsli3,lonsli%data,lonsli%data, + lonsli%nx*lonsli%ny,f(ixf)%difftype) if (allocated(tmpsli3)) deallocate(tmpsli3) endif 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) ! (note stream function not calculated for height-only fields) ! ! 7/9/01: for timeav, skipping this for now.. ! 8/7/07: adding timeave.. else ! vtype=='HEIGHT' (height-only field) lonsli%ny = nhscale allocate(lonsli%yy(lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%yy") lonsli%yy = hscale ! array op allocate(lonsli%data(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mklons allocating lonsli%data") if (ilon_timeave > 0) then timeav(ilon,ixf,iht)%ny = nhscale allocate(timeav(ilon,ixf,iht)%yy(nhscale),stat=ier) if (ier /= 0) | call allocerr(ier,"mklons allocating timeav%yy") timeav(ilon,ixf,iht)%yy = hscale ! write(6,"('timeave height-only field: ilon=',i3, ! | ' ixf=',i3,' iht=',i3,' timeav%ny=',i3, ! | ' timeav%yy=',/,(8f8.3))") ilon,ixf,iht, ! | timeav(ilon,ixf,iht)%ny,timeav(ilon,ixf,iht)%yy endif ! timeave if (ilon_timeave > 0.and.itime==1) then allocate(timeav(ilon,ixf,iht)%data(lonsli%nx, | lonsli%ny),stat=ier) if (ier /= 0) | call allocerr(ier,"mklons allocating timeav%data") write(6,"('mklons ht-only field: allocated ', | 'timeav%data')") if (diffs) then allocate(timeav_cntr(ilon,ixf,iht)%data(lonsli%nx, | lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier, | "mklons allocating timeav_cntr%data") write(6,"('mklons ht-only field: allocated ', | 'timeav_cntr%data')") endif endif ! timeave do k=1,nhscale if (hscale(k) < fhscale(1).or. + hscale(k) > fhscale(nfhscale)) then lonsli%data(:,k) = spval if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data(:,k) = spval else ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev) if (lonsli%lontype=='zm ') then ! zonal means if (.not.diffs) then lonsli%data(:,k) = tmpsli1(ixlat0:ixlat1,ix) else call mkdiffs(tmpsli1(ixlat0:ixlat1,k), + tmpsli2(ixlat0:ixlat1,k),lonsli%data(:,k), + nlats,f(ixf)%difftype) if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data(:,k) = | tmpsli2(ixlat0:ixlat1,k) endif ! diffs if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data(:,k) = tmpsli1(:,k) else ! lon or slt if (.not.diffs) then lonsli%data(:,k) = + f(ixf)%data(ixlon,ixlat0:ixlat1,ix) else call mkdiffs(f(ixf)%data(ixlon,ixlat0:ixlat1,ix), + fcntr(ixf)%data(ixlon,ixlat0:ixlat1,ix), + lonsli%data(:,k),nlats,f(ixf)%difftype) if (ilon_timeave > 0) | timeav_cntr(ilon,ixf,iht)%data(:,k) = | fcntr(ixf)%data(ixlon,ixlat0:ixlat1,ix) endif if (ilon_timeave > 0) | timeav(ilon,ixf,iht)%data(:,k) = | f(ixf)%data(ixlon,ixlat0:ixlat1,ix) endif ! zm or lon/slt endif ! hscale within range enddo ! k=1,nhscale endif ! height-only field ! ! Release temporary slices: if (allocated(tmpsli1)) deallocate(tmpsli1) if (allocated(tmpsli2)) deallocate(tmpsli2) if (allocated(tmpsli3)) deallocate(tmpsli3) endif ! zp or ht ! - - - - - - - - - - - - end zp or ht on y-axis: - - - - - - - - - - - ! ! Take log10 if requested: if ((ilon_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + ilon_log10==2) then call log10f(lonsli%data,lonsli%nx*lonsli%ny,1.e-20,spval) lonsli%log10 = 1 else lonsli%log10 = 0 endif ! ! Take running mean of time averaged data: ! (timeav_prev() are structures with data from previous time) ! if (ilon_timeave > 0) then if (itime==1) then do i=1,timeav(ilon,ixf,iht)%ny do ii=1,timeav(ilon,ixf,iht)%nx if (timeav(ilon,ixf,iht)%data(ii,i)/=spval) then timeav(ilon,ixf,iht)%data(ii,i) = | timeav(ilon,ixf,iht)%data(ii,i)/float(ntimes) endif enddo enddo allocate(timeav_prev(ilon,ixf,iht)%data( | timeav(ilon,ixf,iht)%nx,timeav(ilon,ixf,iht)%ny)) else ! itime > 1 do i=1,timeav(ilon,ixf,iht)%ny do ii=1,timeav(ilon,ixf,iht)%nx if (timeav(ilon,ixf,iht)%data(ii,i)/=spval) then timeav(ilon,ixf,iht)%data(ii,i) = | timeav_prev(ilon,ixf,iht)%data(ii,i)+ | (timeav(ilon,ixf,iht)%data(ii,i)/float(ntimes)) endif enddo enddo endif timeav_prev(ilon,ixf,iht)%data = timeav(ilon,ixf,iht)%data ! ! If diffs, take running mean of time-averaged control case. ! if (diffs) then if (itime==1) then ! no sum yet do i=1,timeav(ilon,ixf,iht)%ny do ii=1,timeav(ilon,ixf,iht)%nx if (timeav_cntr(ilon,ixf,iht)%data(ii,i)/=spval) | then timeav_cntr(ilon,ixf,iht)%data(ii,i) = | timeav_cntr(ilon,ixf,iht)%data(ii,i)/ | float(ntimes) endif enddo enddo ! if (associated(timeav_cntr_prev(ilon,ixf,iht)%data)) ! | deallocate(timeav_cntr_prev(ilon,ixf,iht)%data) allocate(timeav_cntr_prev(ilon,ixf,iht)%data( | timeav(ilon,ixf,iht)%nx,timeav(ilon,ixf,iht)%ny)) else ! itime > 1 -> sum with previous mtime do i=1,timeav(ilon,ixf,iht)%ny do ii=1,timeav(ilon,ixf,iht)%nx if (timeav_cntr(ilon,ixf,iht)%data(ii,i)/=spval) | then timeav_cntr(ilon,ixf,iht)%data(ii,i) = | timeav_cntr_prev(ilon,ixf,iht)%data(ii,i)+ | (timeav_cntr(ilon,ixf,iht)%data(ii,i)/ | float(ntimes)) endif enddo enddo endif timeav_cntr_prev(ilon,ixf,iht)%data = | timeav_cntr(ilon,ixf,iht)%data endif ! diffs ! ! If last time, take diffs of time-averaged data. ! (timeav = timeav-timeav_cntr) if (itime==ntimes.and.diffs) then call mkdiffs(timeav(ilon,ixf,iht)%data, | timeav_cntr(ilon,ixf,iht)%data, | timeav(ilon,ixf,iht)%data, | timeav(ilon,ixf,iht)%nx*timeav(ilon,ixf,iht)%ny, | f(ixf)%difftype) endif endif ! ilon_timeave > 0 ! ! Make plots: if (iplot > 0) then ! making plots ! ! Set up viewport (for multiplt): call setmultivp(vplon,iadvfr,nppf+1,multiplt,ipltrowcol, + vp) ! ! Contour: ! (if doing time averaged data, make these single-ut plots only if ! ilon_timeave_only > 0): ! if (ilon_timeave == 0 .or. ilon_timeave_only == 0) then if (lonsli%lontype /= 'str') then call pltlon(lonsli,zsli(ixlat0:ixlat1,:), | ixlat1-ixlat0+1,npress,ilon_yaxright,vp,h%ut) else ! ! subroutine pltstrm(f,fht,zp,xx,yy,nx,ny,nzp,dzp,glon,iyax,iht, ! + p0,vp,spval) ! call pltstrm(lonsli%data,zsli(ixlat0:ixlat1,:), + gcmlev,lonsli%xx,lonsli%yy,lonsli%nx,lonsli%ny, + npress,dlev,zmflag,ilon_yaxright,iht,p0,vp,spval) endif ! ! Add top and bottom labels: istimeave = .false. call mklonlabs(lonsli,h,hcntr,msgout) call wrlab6(lonsli%tlabs,toffset,tlabchsz, + lonsli%blabs,boffset,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, + 'longitude slices',iframe) endif ! ! Plot time averaged data (last mtime only): if (ilon_timeave > 0.and.itime==ntimes) then ! ! Take log10 of time averaged data if requested: if (allocated(tmpsli1)) deallocate(tmpsli1) if ((ilon_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + ilon_log10==2) then ! ! 3/5/09 btf: Save timeav data before taking log10, so it can be passed ! to sub timeav_glbm below. This avoids calculating global mean of log10 ! data in timeav_glbm (and avoids taking log again of the glb means in pltmxy) ! allocate(tmpsli1(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier, + "mklons alloc tmpsli1 for timeav data") tmpsli1(:,:) = timeav(ilon,ixf,iht)%data(:,:) call log10f(timeav(ilon,ixf,iht)%data, | lonsli%nx*lonsli%ny,1.e-20,spval) timeav(ilon,ixf,iht)%log10 = 1 else timeav(ilon,ixf,iht)%log10 = 0 endif ! ! Plot time averaged data: if (lonsli%lontype /= 'str') then call pltlon(timeav(ilon,ixf,iht),zsli(ixlat0:ixlat1,:) | ,ixlat1-ixlat0+1,npress,ilon_yaxright,vp,h%ut) else call pltstrm(timeav(ilon,ixf,iht)%data, | zsli(ixlat0:ixlat1,:),gcmlev,lonsli%xx,lonsli%yy, | timeav(ilon,ixf,iht)%nx,timeav(ilon,ixf,iht)%ny, | npress,dlev,zmflag,ilon_yaxright,iht,p0,vp,spval) endif istimeave = .true. call mklonlabs(timeav(ilon,ixf,iht),h,hcntr,msgout) call wrlab6(timeav(ilon,ixf,iht)%tlabs,toffset,tlabchsz, | timeav(ilon,ixf,iht)%blabs,boffset,blabchsz, | vp,ilab_hq) nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, | iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, | 'time averaged longitude slices',iframe) ! ! Make global mean line plots of time-averaged data: ! First restore the non-log data if log was taken: if (allocated(tmpsli1)) then timeav(ilon,ixf,iht)%data(:,:) = tmpsli1(:,:) deallocate(tmpsli1) endif call timeav_glbm(timeav(ilon,ixf,iht),lonsli%nx, | lonsli%ny,ixlat0,ixlat1,timeav(ilon,ixf,iht)%tlabs, | timeav(ilon,ixf,iht)%blabs,fmin,fmax,vp) write(msgout,"(a,' global means min,max=',2e11.4, | ' iht= time-ave')") lonsli%sname,fmin,fmax nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, | iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, | 'global mean time averaged lon slices',iframe) endif ! plot time averaged data ! ! If not making plots, wrout will need zmin,zmax. else call fminmax(lonsli%data,size(lonsli%data),zmin,zmax) endif ! iplot > 0 ! ! Write output data files: call wrout_lons(lonsli,h,hcntr) call wrcdf_lons(lonsli,h,hcntr) enddo yax_loop enddo lon_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 (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,' ','longitude slices', + iframe) ! ! Release space (including lonsli pointers): if (associated(lonsli%data)) deallocate(lonsli%data) if (associated(lonsli%xx)) deallocate(lonsli%xx) if (associated(lonsli%yy)) deallocate(lonsli%yy) if (associated(htscale)) deallocate(htscale) if (associated(zprange)) deallocate(zprange) if (associated(fhscale)) deallocate(fhscale) if (associated(hscale)) deallocate(hscale) if (allocated(tmpsli1)) deallocate(tmpsli1) if (allocated(tmpsli2)) deallocate(tmpsli2) if (allocated(tmpsli3)) deallocate(tmpsli3) return contains !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine mkepvec(fy,fz,init,iht,zviewfac,scale) ! ! This routine is internal to mklons. ! Prepare epflux vectors Y+Z, defining lonsli%y and lonsli%z. ! Field may be EPVYZ (vectors only), EPVYZMAG, or EPVDIV. If field is ! EPVYZ or EPVYZMAG, then also redefine lonsli%data as the vector ! sum of Y+Z. ! If init > 0, do only initialization (including setting lonsli%vectype). ! If init == 0, do not do initialization (define lonsli%y,z,data) ! If init < 0, do only scaling and normalization of lonsli%y,z. ! ! Args: type(field),intent(in) :: fy,fz integer,intent(in) :: init,iht real,intent(in) :: zviewfac,scale ! ! Locals: integer :: ier,j,k,jj,kk,len real :: epnorm,fmin,fmax 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 ! ! Init flag is set -> set lonsli%vectype, and insure presence of ! EPVY and EPVZ vectors. Also set vector cutoffs, etc. ! (Called from inside field loop but before lon and iht loops) ! if (init > 0) then ! initialize only if (.not.associated(fy%data).or. + .not.associated(fz%data)) then write(6,"('>>> mkepvec (mklons): need EPVY and EPVZ data', + 'to make epflux vectors: vectors will not be plotted.')") return endif ! ! vectype=='EPVYZ' -> plotting epflux vectors alone ! vectype=='+EPVYZ' -> plotting epflux vectors over epvdiv or epvyzmag ! vectype may remain blank after this, so vector magnitudes may be ! calculated below and put in lonsli%data for EPVYZ or EPVYZMAG. ! (see comments below) ! if (lonsli%sname=='EPVYZ ') then write(lonsli%vectype,"('EPVYZ')") elseif (lonsli%sname=='EPVDIV '.and.ilon_epvdiv_yz > 0) then write(lonsli%vectype,"('+EPVYZ')") elseif (lonsli%sname=='EPVYZMAG'.and.ilon_epvyzmag_yz > 0) then write(lonsli%vectype,"('+EPVYZ')") endif ! ! vlc,vhc,vsc are in module plt: vlc = vep_scale(1) vhc = vep_scale(2) vsc = vep_scale(3) vrl = 0. ! vector realized length if (vmag_len>0.) vrl = vmag_len return endif ! init if (init < 0) goto 100 ! do scaling and/or normalization only ! ! Init flag not set -> define lonsli%y(:,:) and lonsli%z(:,:) ! (Called from inside field, lon, and iht loops. Has been called ! previously with init==1) ! if (lonsli%lontype=='str') return if (associated(lonsli%y)) deallocate(lonsli%y) if (associated(lonsli%z)) deallocate(lonsli%z) allocate(lonsli%y(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier,"mkepvec allocating lonsli%y") allocate(lonsli%z(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) call allocerr(ier,"mkepvec allocating lonsli%z") ! ! iht==0 -> at zp, iht==1 -> at ht. ! ! Zp pressure surface: ! (use fislice to appease sgi compiler) ! if (iht==0) then ! at zp if (lonsli%lontype=='zm ') then do j=ixlat0,ixlat1 jj = j-ixlat0+1 do k=k0zp,k1zp kk = k-k0zp+1 lonsli%y(jj,kk)=fmean(fy%data(:,j,k),nlon,spval,0) lonsli%z(jj,kk)=fmean(fz%data(:,j,k),nlon,spval,0) enddo enddo else ! lon/slt lonsli%y = fy%data(ixlon,ixlat0:ixlat1,k0zp:k1zp) lonsli%z = fz%data(ixlon,ixlat0:ixlat1,k0zp:k1zp) endif ! zm or lon/slt ! ! Height surface: ! (epfluxes fy,fz, etc are at interface pressure coords, so can use npress ! or f(ixf)%nlev) ! else ! at ht if (allocated(tmpsli1)) deallocate(tmpsli1) allocate(tmpsli1(nlat,npress),stat=ier) if (ier /= 0) call allocerr(ier,"mkepv allocating tmpsli1") if (allocated(tmpsli2)) deallocate(tmpsli2) allocate(tmpsli2(nlat,npress),stat=ier) if (ier /= 0) call allocerr(ier,"mkepv allocating tmpsli2") ! ! Save zonal means in tmpsli1 and tmpsli2: if (lonsli%lontype=='zm ') then do j=1,nlat do k=1,npress ! fislice = fy%data(:,j,k) tmpsli1(j,k)=fmean(fy%data(:,j,k),nlon,spval,0) ! fislice = fz%data(:,j,k) tmpsli2(j,k)=fmean(fz%data(:,j,k),nlon,spval,0) enddo enddo ! note nlats==ixlat1-ixlat0+1 call cuthtint(tmpsli1(ixlat0:ixlat1,:),zsli(ixlat0:ixlat1,:), + nlats,npress,lonsli%y,htscale,nhtscale,logint,spval,ier,1) call cuthtint(tmpsli2(ixlat0:ixlat1,:),zsli(ixlat0:ixlat1,:), + nlats,npress,lonsli%z,htscale,nhtscale,logint,spval,ier,1) else ! lon/slt call cuthtint(fy%data(ixlon,ixlat0:ixlat1,:), + zsli(ixlat0:ixlat1,:),nlats,npress,lonsli%y,htscale, + nhtscale,logint,spval,ier,1) call cuthtint(fz%data(ixlon,ixlat0:ixlat1,:), + zsli(ixlat0:ixlat1,:),nlats,npress,lonsli%z,htscale, + nhtscale,logint,spval,ier,1) endif if (allocated(tmpsli1)) deallocate(tmpsli1) if (allocated(tmpsli2)) deallocate(tmpsli2) endif ! zp or ht ! ! Redefine lonsli%data and min,max for EPVYZ or EPVYZMAG: ! Getflds defined EPVYZ and EPVYZMAG as vector sum of EPVY+EPVZ. ! But since we may be plotting zonal means, lonsli%data should be ! redefined as vector sum of the means of EPVY and EPVZ (currently ! in lonsli%y and lonsli%z). This way we are assured to get the ! vector sum of the means zm(y)+zm(z) rather than means of the ! vector sum y+z. ! if (lonsli%sname /= 'EPVDIV ') then lonsli%data = vecsum(lonsli%y,lonsli%z, + size(lonsli%y,1),size(lonsli%y,2),spval) call fminmax(lonsli%data,size(lonsli%data),zmin,zmax) endif ! ! If not doing vectors, can skip the rest: if (len_trim(lonsli%vectype)==0) return ! ! vepz_viewfac is user provided factor to be multiplied by epvz ! so vector plots will not be dominated by the y component. ! Hanli recommends vepz_viewfac = pi*ar/(zmax-zmin)*aspect_ratio ! (pi*ar is half earth circumference, zmax/zmin is max/min altitude ! begin shown) (typically vepz_viewfac is in the range 50-100) ! ! (jump to here if init < 0) 100 continue if (zviewfac /= 1.) then lonsli%z = lonsli%z * vepz_viewfac if (index(lonsli%vectype,'Z*')==0) then len = len_trim(lonsli%vectype) write(lonsli%vectype(len+1:len+12),"(' (Z*',f7.1,')')") + vepz_viewfac endif endif ! ! If user has set vep_scale(3) < 0., then normalize the vectors: if (scale < 0.) then do k=1,lonsli%ny do j=1,lonsli%nx epnorm = sqrt(lonsli%y(j,k)**2+lonsli%z(j,k)**2) lonsli%y(j,k) = lonsli%y(j,k)/epnorm lonsli%z(j,k) = lonsli%z(j,k)/epnorm enddo enddo endif end subroutine mkepvec end subroutine mklons !------------------------------------------------------------------- subroutine mkampha(f,fcntr,nf,h,hcntr) ! ! Make amplitudes and phases for lon slices: ! use proc use fields,only: field use fields,only: field,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use input use plt use set_vert use ampha ! ! Args: integer,intent(in) :: nf type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: ixlat0,ixlat1,k0zp,k1zp,nzprange,nhtscale,iadvfr, + izpht0,izpht1,ixf,iap,iw,ixz,iht,logint,ier,j,k type(lonslice) :: lonsli character(len=80) :: msgout real,pointer :: htscale(:),zprange(:) ! ! amp(nlat,:,mxwave), phase(nlat,:,mxwave) are amplitudes and phases, ! where 2nd dimension will be npress or nhtscale ! fht(nlon,nlat,nhtscale) is global ht-interpolated field real,allocatable :: phase(:,:,:),amp(:,:,:),fht(:,:,:), | fht_cntr(:,:,:) real :: fzp(nlon,nlat,npress) ! ! 12/7/06 btf: Changed dimension of work array W from (200) to ! (2*nlon+nlon/2+15), as per rfft.F doc (this code was crashing ! with dres timegcm history, where nlon==145) ! real :: x(nlon),w(2*nlon+nlon/2+15),fmin,fmax,scale, | flonsli(nlat,npress),xleft,xright real :: fhi(mxwave) = (/24.,18.,16.,15./) real :: flo(mxwave) = (/0.,6.,8.,9./) real :: vp(4),vplon(4)=(/.13,.87,.23,.88/) ! for 0 extra right axes real :: toffset=.06,tlabchsz(3) = (/.02,.02,.02/) real :: boffset=.20,blabchsz(3) = (/.018,.018,.018/) real :: zsli(nlat,npress) ! slice of heights (for zm) ! ! Externals: real,external :: fmean integer,external :: ixfindc ! if (amphase <= 0) return write(6,"(/'Amplitudes and Phases at ut = ',f5.2,' mtime = ', + 3i3)") h%ut,h%mtime if (diffs) then write(6,"('Will calculate amplitudes and phases of ', | 'difference fields')") endif if (multiadvfr > 0) nppf = 0 lonsli%vectype = ' ' ! no vector plotting ! ! Define x coords in latitude: ! (User may optionally restrict latitude range with flon_xlatrange) ! (,flon_xlatrange is from user, ixlat0,ixlat1 are returned) call setxlat(lonsli,flon_xlatrange,ixlat0,ixlat1) xleft = lonsli%xx(1)-(lonsli%xx(2)-lonsli%xx(1))/2. xright = lonsli%xx(lonsli%nx)+(lonsli%xx(2)-lonsli%xx(1))/2. ! ! Adjust right viewport according to number of extra right-hand y-axes: if (ilon_yaxright==0) then vplon(2) = .90 elseif (ilon_yaxright==1) then vplon(2) = .85 else vplon(2) = .76 endif ! ! Locate index to heights field (needed for interp, etc): ! ixz = ixfindc(f%fname8,nf,'Z ') if (ixz <= 0) then write(6,"('>>> mkampha: need heights in fields: field names', + ' = ',(8a8))") f%fname8 stop 'mklons z' endif ! ! Field loop. Field must be requested, data pointer must be associated, ! and field must not be height-independent: ! (add height-only fields ohv/b later) fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data).or. + f(ixf)%nlev==1.or.trim(f(ixf)%vtype)=='HEIGHT') then ! ! 3/9/06 btf: Had to add this conditional to prevent loop from ! continuing after cycling when ixf==nf. I have *no* idea why ! it continued iterations with ixf > nf in this case, but ! exiting the loop rather than cycling when ixf >= nf as below ! solved the problem. ! if (ixf < nf) then cycle fields_loop else exit fields_loop endif endif ! ! Set up vertical scale(s) zp and/or ht: if (f(ixf)%zptype(1:3)=='MID') then ! midpoints ! ! Save zonal mean heights: do j=ixlat0,ixlat1 do k=1,npress zsli(j,k) = fmean(z_midpts(:,j,k),nlon,spval,0) enddo enddo call setvert(flon_zprange,gcmlev,npress,dlev,nzprange,k0zp, | k1zp,zprange,flon_htscale,nhtscale,htscale,spval,1) else ! interfaces do j=ixlat0,ixlat1 do k=1,npress zsli(j,k) = fmean(z_ifaces(:,j,k),nlon,spval,0) enddo enddo call setvert(flon_zprange,gcmilev,npress,dlev,nzprange,k0zp, | k1zp,zprange,flon_htscale,nhtscale,htscale,spval,1) endif ! ! Set loop indices for zp and/or ht on y-axis: izpht0 = 0 izpht1 = 1 if (nhtscale == 0) izpht1 = 0 if (nzprange == 0) izpht0 = 1 ! lonsli%fname = f(ixf)%fname56 lonsli%sname = f(ixf)%fname8 lonsli%funits = f(ixf)%units lonsli%ftype = f(ixf)%type lonsli%zptype = f(ixf)%zptype lonsli%known = f(ixf)%known lonsli%difftype = f(ixf)%difftype ! write(6,"('mkampha: ixf=',i3,' lonsli%fname=',a,' funits=', ! | a)") ixf,trim(lonsli%fname),trim(lonsli%funits) ! ! 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 logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 ! ! Zp and/or ht on y-axis (0,0 for zp only, 1,1 for ht only, ! or 0,1 for zp and ht): yax_loop: do iht=izpht0,izpht1 lonsli%ihtyax = iht if (associated(lonsli%yy)) deallocate(lonsli%yy) if (associated(lonsli%data)) deallocate(lonsli%data) ! ! Get amplitudes and phases for "amphase" wave numbers: ! (amphase is integer provided by user, in range 1-4) ! amp(nlat,:,mxwave), phase(nlat,:,mxwave) are amplitudes and phases, ! where 2nd dimension will be npress or nhtscale ! if (allocated(amp)) deallocate(amp) if (allocated(phase)) deallocate(phase) k = npress ! zp if (iht==1) k=nhtscale ! ht allocate(amp(nlat,k,mxwave),stat=ier) if (ier /= 0) + call allocerr(ier,"mkamphase allocating amp") allocate(phase(nlat,k,mxwave),stat=ier) if (ier /= 0) + call allocerr(ier,"mkamphase allocating phase") ! ! Zp on y-axis: if (iht==0) then ! zp on y-axis fzp = f(ixf)%data ! nlon,nlat,npress if (diffs) then ! ! Make difference fields fzp = fzp-fcntr(ixf)%data ! subroutine mkdiffs(pert,cntr,diffs,n,difftype) ! Then getampha will calculate amp and phase from the differences ! call mkdiffs(fzp,fcntr(ixf)%data,fzp,nlon*nlat*npress, | f(ixf)%difftype) endif call getampha(amphase,nlon,nlon,nlat,npress,fzp, + amp,phase,x,w,h%ut,spval) call fixphase(phase,nlat,npress,mxwave, + (f(ixf)%fname8=='N2 '),spval) ! ! Define y coords in log pressure (from setvert, called above): lonsli%ny = nzprange allocate(lonsli%yy(lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkampha allocating lonsli%yy") lonsli%yy = zprange ! array op write(lonsli%ylab,"('LN(P0/P) (',a,')')") | trim(lonsli%zptype) ! ! Height on y-axis: else ! ht on y-axis ! ! fht will be globally ht-interpolated field: if (.not.allocated(fht)) + allocate(fht(nlon,nlat,nhtscale),stat=ier) if (ier /= 0) call allocerr(ier,"mkamphase allocating fht") if (diffs) then if (.not.allocated(fht_cntr)) + allocate(fht_cntr(nlon,nlat,nhtscale),stat=ier) if (ier /= 0) call allocerr(ier,"mkamphase allocating', | ' fht_cntr") endif ! ! Height interpolation: if (trim(f(ixf)%zptype) == 'INTERFACES') then call glbhtin(f(ixf)%data,z_ifaces,nlon,nlat,npress, + fht,htscale,nhtscale,logint,spval,ier,0) if (diffs) call glbhtin(fcntr(ixf)%data,z_ifaces,nlon, | nlat,npress,fht_cntr,htscale,nhtscale,logint,spval, | ier,0) else call glbhtin(f(ixf)%data,z_midpts,nlon,nlat,npress, + fht,htscale,nhtscale,logint,spval,ier,0) if (diffs) call glbhtin(fcntr(ixf)%data,z_midpts,nlon, | nlat,npress,fht_cntr,htscale,nhtscale,logint,spval, | ier,0) endif ! ! Make difference fields fht = fht-fht_cntr (getampha will calculate ! amp and phases of the differences): if (diffs) then call mkdiffs(fht,fht_cntr,fht,nlon*nlat*nhtscale, | f(ixf)%difftype) endif ! ! Pass ht-interpolated field to getampha: ! call getampha(amphase,nlon,nlon,nlat,nhtscale,fht, + amp,phase,x,w,h%ut,spval) call fixphase(phase,nlat,nhtscale,mxwave, + (f(ixf)%fname8=='N2 '),spval) ! ! Define y coords in height (from setvert, called above): lonsli%ny = nhtscale allocate(lonsli%yy(lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkampha allocating lonsli%yy") lonsli%yy = htscale ! array op write(lonsli%ylab,"('HEIGHT (KM)')") endif ! zp or ht on y-axis ! ! Allocate and define lonsli data (lonsli%yy is zp or ht): allocate(lonsli%data(lonsli%nx,lonsli%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkampha allocating lonsli%data") ! ! Loop for wave numbers: wave_loop: do iw = 1,amphase lonsli%nwave = iw ! ! Loop for amplitudes, then phases: do iap = 1,2 if (iap==1) then ! amplitudes lonsli%lontype = 'amp' pltmin = f(ixf)%cmin pltmax = f(ixf)%cmax conint = f(ixf)%cint else ! phases lonsli%lontype = 'pha' pltmin = flo(iw) pltmax = fhi(iw) conint = (pltmax-pltmin)/12. scale = (pltmax-pltmin)/(2.*pi) endif ! ! Define lonsli%data (was allocated above): ! ! Zp on y-axis (use flonsli as temporary output): if (iht == 0) then if (iap == 1) then ! amplitudes call defamp(amp,nlat,npress,mxwave,flonsli,iw,spval) else ! phases call defphase(phase,nlat,npress,mxwave,flonsli,scale, + iw,spval) endif lonsli%data = flonsli(ixlat0:ixlat1,k0zp:k1zp) ! ! Height on y-axis (use fht(1,:,:) as temporary output): else if (iap == 1) then ! amplitudes call defamp(amp,nlat,nhtscale,mxwave,fht(1,:,:),iw, + spval) else ! phases call defphase(phase,nlat,nhtscale,mxwave,fht(1,:,:), + scale,iw,spval) endif lonsli%data = fht(1,ixlat0:ixlat1,:) endif ! ! Take log10 if requested (amplitudes only): if ((ilon_log10==1.and.trim(f(ixf)%type)=='DENSITY'.and. + iap==1).or.(ilon_log10==2.and.iap==1)) then call log10f(lonsli%data,lonsli%nx*lonsli%ny,1.e-20, + spval) lonsli%log10 = 1 else lonsli%log10 = 0 endif ! ! Make plots: if (iplot > 0) then ! ! Set up viewport (for multiplt): call setmultivp(vplon,iadvfr,nppf+1,multiplt,ipltrowcol, + vp) ! ! Contour: if (iap==1) then call pltlon(lonsli,zsli(ixlat0:ixlat1,:), + ixlat1-ixlat0+1,npress,ilon_yaxright,vp,h%ut) else ! ! Use old modified conrec to plot phases: ! subroutine pltconrec(plt,idx,nx,ny,xx,yy,xlab,ylab, ! + flo,fhi,fint,vp,wl,wr,wb,wt,spv,iyax,fht,zp,nzp,p0) ! call pltconrec(lonsli%data,lonsli%nx,lonsli%nx, + lonsli%ny,lonsli%xx,lonsli%yy,trim(lonsli%xlab), + trim(lonsli%ylab),lonsli%ihtyax,pltmin,pltmax, + conint,vp,xleft,xright,lonsli%yy(1), + lonsli%yy(lonsli%ny),spval,ilon_yaxright, + zsli(ixlat0:ixlat1,:),gcmlev,npress,p0) endif ! ! Add top and bottom labels: istimeave = .false. call mklonlabs(lonsli,h,hcntr,msgout) call wrlab6(lonsli%tlabs,toffset,tlabchsz, + lonsli%blabs,boffset,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, + 'amp-phase slices',iframe) ! ! If not making plots, wrout will need zmin,zmax. else ! no plots call fminmax(lonsli%data,size(lonsli%data),zmin,zmax) endif ! iplot > 0 ! ! Write output data files: call wrout_lons(lonsli,h,hcntr) call wrcdf_lons(lonsli,h,hcntr) enddo ! iap=1,2 for amp,phase enddo wave_loop enddo yax_loop enddo fields_loop ! ! Release space (altho maybe lonsli goes away anyway because ! it's local?) if (associated(lonsli%data)) deallocate(lonsli%data) if (associated(lonsli%xx)) deallocate(lonsli%xx) if (associated(lonsli%yy)) deallocate(lonsli%yy) if (associated(htscale)) deallocate(htscale) if (associated(zprange)) deallocate(zprange) if (allocated(fht)) deallocate(fht) if (allocated(amp)) deallocate(amp) if (allocated(phase)) deallocate(phase) if (allocated(fht_cntr)) deallocate(fht_cntr) end subroutine mkampha !------------------------------------------------------------------- subroutine setxlat(lonsli,xlatrange,ixlat0,ixlat1) use proc ! ! Set up latitude on x-axis of longitude slice lonsli. ! xlatrange(2) is user provided min,max. Return ixlat0,ixlat1 ! as indices to 1st and last latitude within gcmlat. ! This routine defines lonsli%nx and lonsli%xlag, and allocates ! and defines lonsli%xx(nx). ! ! Args: type(lonslice),intent(inout) :: lonsli real,intent(in) :: xlatrange(2) integer,intent(out) :: ixlat0,ixlat1 ! ! Locals: integer :: ier ! ! Externals: integer,external :: ixfind ! lonsli%nx = nlat ixlat0 = 1 ixlat1 = nlat if ((xlatrange(1)/=spval.and.xlatrange(2)/=spval).and. + (xlatrange(1)>gcmlat(1).or.xlatrange(2) 1) ixlat0 = ixlat0+1 if (ixlat1-ixlat0 > 0) then lonsli%nx = ixlat1-ixlat0+1 endif endif ! IBM croaks on this test: ! if (.not.associated(lonsli%xx)) then allocate(lonsli%xx(lonsli%nx),stat=ier) if (ier /= 0) call allocerr(ier,"setxlat allocating lonsli%xx") ! endif lonsli%xx = gcmlat(ixlat0:ixlat1) ! array op write(lonsli%xlab,"('LATITUDE (DEG)')") end subroutine setxlat !------------------------------------------------------------------- subroutine mklonlabs(lonsli,h,hcntr,msgout) use plt,only: zmin,zmax,ciu,scfac,pltmin,pltmax,conint use proc,only: spval,diffs use input,only: ie5577,ie6300,mtimes,ntimes,ilon_timeave use hist,only: hdr ! ! Construct 3 top labels (lonsli%tlabs) and 3 bottom labels ! (lonsli%blabs) for lon slice: ! ! Args: type(lonslice),intent(inout) :: lonsli type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Locals: integer :: lenlab, | iyear ! 4-digit year real :: dthrs character(len=16) :: char16 ! ! External: integer,external :: mtime_to_mins ! ! Set up top and bottom text labels: ! ! Top (1st) top label is optionally provided by user: lonsli%tlabs(:) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if lonsli%log10 > 0) ! if (lonsli%known) then lenlab = len_trim(lonsli%fname)+len_trim(lonsli%funits)+3 if (lonsli%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(lonsli%tlabs(2))) then if (lonsli%log10 > 0) then lonsli%tlabs(2) = 'LOG10 '//trim(lonsli%fname)//' ('// + trim(lonsli%funits)//')' else lonsli%tlabs(2) = trim(lonsli%fname)//' ('// + trim(lonsli%funits)//')' endif else lonsli%tlabs(2) = trim(lonsli%fname) endif ! write(6,"('mklonlabs: lonsli%fname=',a,' funits=',a, ! | /,' tlabs(2)=',a)") trim(lonsli%fname),trim(lonsli%funits), ! | trim(lonsli%tlabs(2)) else ! unknown to proc lenlab = len_trim(lonsli%sname)+2+len_trim(lonsli%fname)+ | len_trim(lonsli%funits)+3 if (lonsli%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(lonsli%tlabs(2))) then if (lonsli%log10 > 0) then if (trim(lonsli%sname)/=trim(lonsli%fname)) then lonsli%tlabs(2) = 'LOG10 '//trim(lonsli%sname)//': '// | trim(lonsli%fname)//' ('//trim(lonsli%funits)//')' else lonsli%tlabs(2) = 'LOG10 '//trim(lonsli%fname)// | ' ('//trim(lonsli%funits)//')' endif else ! no log10 if (trim(lonsli%sname)/=trim(lonsli%fname)) then lonsli%tlabs(2) = trim(lonsli%sname)//': '// | trim(lonsli%fname)//' ('//trim(lonsli%funits)//')' else lonsli%tlabs(2) = trim(lonsli%fname)// | ' ('//trim(lonsli%funits)//')' endif endif else lonsli%tlabs(2) = trim(lonsli%fname) endif endif ! known or unknown ! ! Add DIFFS to label: if (diffs) then if (trim(lonsli%difftype)=='RAW') then ! if (len_trim(lonsli%tlabs(2))+7 <= len(lonsli%tlabs(2))) then ! lonsli%tlabs(2) = 'DIFFS: '//trim(lonsli%tlabs(2)) ! else ! lonsli%tlabs(2) = 'DIFFS: '//trim(lonsli%fname) ! endif lonsli%tlabs(2) = 'DIFFS: '//trim(lonsli%fname) else ! percent ! if (len_trim(lonsli%tlabs(2))+15 <= len(lonsli%tlabs(2))) then ! lonsli%tlabs(2) = trim(lonsli%difftype)//' DIFFS: '// ! | trim(lonsli%tlabs(2)) ! else ! lonsli%tlabs(2) = '% DIFFS: '//trim(lonsli%fname) ! endif lonsli%tlabs(2) = '% DIFFS: '//trim(lonsli%fname) endif endif ! ! Bottom (3rd) top label is ut and grid info: ! 7/01: There should be a 4-digit h%iyear ! iyear = h%iyd ! yyddd (only 2-digit year on history?) iyear = iyear/1000 ! yy if (iyear > 10) then iyear = iyear+1900 else iyear = iyear+2000 endif if (.not.diffs) then write(char16,"('YEAR ',i4)") iyear else ! do not use iyear if doing diffs char16 = ' ' endif if (lonsli%lontype=='lon'.or.lonsli%lontype=='slt') then write(lonsli%tlabs(3),"(a,' UT=',f5.2,' LON=',f9.2, | ' (DEG) SLT=',f5.2,' (HRS)')") trim(char16),h%ut,lonsli%glon, | lonsli%slt elseif (lonsli%lontype == 'zm ') then write(lonsli%tlabs(3),"(a,' UT=',f5.2,' ZONAL MEANS')") | trim(char16),h%ut elseif (lonsli%lontype == 'amp') then write(lonsli%tlabs(3),"(a,' UT=',f5.2,' AMPLITUDE ', | '(WAVE ',i1,')')") trim(char16),h%ut,lonsli%nwave elseif (lonsli%lontype == 'pha') then write(lonsli%tlabs(3),"(a,'UT=',f5.2,' PHASE (WAVE ', | i1,')')") trim(char16),h%ut,lonsli%nwave elseif (lonsli%lontype == 'str') then write(lonsli%tlabs(3),"(a,'UT=',f5.2, | ' STREAM FUNCTION')") trim(char16),h%ut else write(6,"('>>> mklonlabs: unknown lontype=',a)") lonsli%lontype endif ! ! Top (1st) bottom label is min,max,interval: ! (zmin,zmax, and ciu are in module plt and were defined by contour, ! unless user specified fixed range and/or interval with namelist ! read fmnmxint (was transferred to pltmin,pltmax,conint, which ! are also in module plt) (scfac is in module plt) ! If doing vectors only (lonsli%vectype=='EPVYZ'), print only min,max ! if (pltmin >= pltmax) then ! user did not provide min,max,int if (len_trim(lonsli%vectype)==0.or.lonsli%vectype(1:1)=='+') + then ! contours if (scfac==1.) then write(lonsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(lonsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac endif else ! vectors only if (scfac==1.) then write(lonsli%blabs(1),"('MAGNITUDE MIN,MAX=',2(1pe12.4))") + zmin,zmax else write(lonsli%blabs(1),"('MAGNITUDE MIN,MAX=',2(1pe12.4), + ' (X',1pe9.2,')')") zmin,zmax,scfac endif endif else ! user provided min,max,int if (len_trim(lonsli%vectype)==0.or.lonsli%vectype(1:1)=='+') + then ! contours if (scfac==1.) then write(lonsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") pltmin,pltmax,conint else write(lonsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac endif else ! vectors only if (scfac==1.) then write(lonsli%blabs(1),"('MAGNITUDE MIN,MAX=',2(1pe12.4))") + pltmin,pltmax else write(lonsli%blabs(1),"('MAGNITUDE MIN,MAX=',2(1pe12.4), + ' (X',1pe9.2,')')") zmin,zmax,scfac endif endif endif ! ! Middle (2nd) bottom label is history info: ! Bottom (3rd) bottom label is also used if doing diffs. ! if (.not.diffs) then ! if (ilon_timeave==0) then ! write(lonsli%blabs(2),"(a,' ',a,' (DAY,HR,MIN=', ! + i3,',',i2,',',i2,')')") trim(h%version),trim(h%mssvol), ! + h%mtime ! else ! dthrs = float(mtime_to_mins(mtimes(:,2))- ! | mtime_to_mins(mtimes(:,1)))/60. ! write(lonsli%blabs(2),"('TIME-AVERAGED FROM ',3i3,' TO ', ! | 3i3,' (DT=',f4.1,' HOURS)')") ! | mtimes(:,1),mtimes(:,ntimes),dthrs ! endif ! lonsli%blabs(3) = ' ' ! if (trim(lonsli%ftype)=='EXCITED-STATE') ! + write(lonsli%blabs(3),"('F10.7=',f6.2)") hdr%f107d ! if (trim(lonsli%ftype)=='EMISSION') then ! call mkemislab(lonsli%sname,ie5577,ie6300,lonsli%blabs(3)) ! if (len_trim(lonsli%blabs(3)) > 48) blabchsz(3) = .015 ! endif ! else ! if (ilon_timeave==0) then ! if (7+len_trim(h%mssvol) +12<=len(lonsli%blabs(2)).and. ! + 7+len_trim(hcntr%mssvol)+12<=len(lonsli%blabs(3)))then ! write(lonsli%blabs(2),"('DIFFS: ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(h%mssvol),h%mtime ! write(lonsli%blabs(3),"('MINUS ',a,' (',i3,',',i2,',',i2, ! + ')')") trim(hcntr%mssvol),hcntr%mtime ! else ! write(lonsli%blabs(2),"('DIFFS: (',i3,',',i2,',',i2, ! + ') MINUS (',i3,',',i2,',',i2,')')") ! + h%mtime,hcntr%mtime ! lonsli%blabs(3) = ' ' ! endif ! else ! time-averaged diffs ! dthrs = float(mtime_to_mins(mtimes(:,2))- ! | mtime_to_mins(mtimes(:,1)))/60. ! write(lonsli%blabs(2),"('DIFFS OF TIME-AVERAGED DATA: ',3i3, ! | ' TO ',3i3)") mtimes(:,1),mtimes(:,ntimes) ! write(lonsli%blabs(3),"('(DT=',f4.1,' HOURS)')") dthrs ! endif ! endif ! ! 5/05: Put version name and model time in blabs(2), history disk ! file(s) in blabs(3) (give up on emission info label): ! write(lonsli%blabs(2),"(a,' (DAY,HR,MIN=', | i3,',',i2,',',i2,')')") trim(h%version),h%mtime ! if (ilon_timeave>0) then if (istimeave) then dthrs = float(mtime_to_mins(mtimes(:,2))- | mtime_to_mins(mtimes(:,1)))/60. ! write(lonsli%blabs(2),"('TIME-AVERAGED FROM ',3i3,' TO ', ! | 3i3,' (DT=',f4.1,' HOURS)')") ! | mtimes(:,1),mtimes(:,ntimes),dthrs write(lonsli%blabs(2),"('TIME-AVERAGE FROM ',3i3,' TO ', | 3i3,' BY ',f4.1,' HRS (',i4,' TIMES)')") | mtimes(:,1),mtimes(:,ntimes),dthrs,ntimes endif if (.not.diffs) then write(lonsli%blabs(3),"(a)") trim(h%histfile) else write(lonsli%blabs(3),"(a,' MINUS ',a)") | trim(h%histfile),trim(hcntr%histfile) endif ! ! Return message to be printed to stdout: if (lonsli%lontype == 'lon') then ! longitude write(msgout,"(a,' glon=',f7.2,' min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,lonsli%glon, + zmin,zmax,lonsli%ihtyax elseif (lonsli%lontype == 'slt') then ! local time write(msgout,"(a,' slt =',f7.2,' min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,lonsli%slt, + zmin,zmax,lonsli%ihtyax elseif (lonsli%lontype == 'amp') then ! amplitudes write(msgout,"(a,' amplitude (wave ',i1,') min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,lonsli%nwave, + zmin,zmax,lonsli%ihtyax elseif (lonsli%lontype == 'pha') then ! phases write(msgout,"(a,' phase (wave ',i1,') min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,lonsli%nwave, + zmin,zmax,lonsli%ihtyax elseif (lonsli%lontype == 'zm ') then ! zonal means write(msgout,"(a,' zonal means min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,zmin,zmax,lonsli%ihtyax elseif (lonsli%lontype == 'str') then ! stream function write(msgout,"(a,' stream func min,max=', + 2e11.4,' iht=',i2)") lonsli%sname,zmin,zmax,lonsli%ihtyax else write(6,"('>>> mklonlabs: unknown lontype=',a)") lonsli%lontype endif if (ilon_timeave > 0) | write(msgout(len_trim(msgout):len_trim(msgout)+8), | "('time-ave')") return end subroutine mklonlabs !------------------------------------------------------------------- subroutine pltlon(lonsli,hts,id1_hts,id2_hts,iyax,vp,ut) use plt use proc,only: gcmlev,gcmilev,npress,dlev,p0,spval ! ! Make longitude slice 2d contour (lat on x-axis, ht or zp on y-axis) ! type(lonslice) lonsli and type(history) h have been fully defined. ! ! Args: type(lonslice),intent(in) :: lonsli integer,intent(in) :: iyax,id1_hts,id2_hts real,intent(in) :: vp(4),ut,hts(id1_hts,id2_hts) ! ! Locals: integer :: iyaxright,iclr,k real :: xleft,xright,yboffset integer :: nnans_z ! number of nans ! ! Set up conpack: call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',lonsli%xx(1)) call cpsetr('XCM',lonsli%xx(lonsli%nx)) call cpsetr('YC1',lonsli%yy(1)) call cpsetr('YCN',lonsli%yy(lonsli%ny)) xleft = lonsli%xx(1)-(lonsli%xx(2)-lonsli%xx(1))/2. xright = lonsli%xx(lonsli%nx)+(lonsli%xx(2)-lonsli%xx(1))/2. call set(vp(1),vp(2),vp(3),vp(4),xleft,xright, + lonsli%yy(1),lonsli%yy(lonsli%ny),1) ! ! Contour: ! If lonsli%vectype== [blank] -> contour only (no vectors) ! If lonsli%vectype== '+EPVYZ' -> add vectors to contours ! If lonsli%vectype== 'EPVYZ' -> vectors only (no contour) ! (if vectors only, zmin,max was already set by mkepvec) ! if (len_trim(lonsli%vectype)==0.or.lonsli%vectype(1:1)=='+') then call contour(lonsli%data,lonsli%nx,lonsli%nx,lonsli%ny) else ! call fminmax(lonsli%data,size(lonsli%data),zmin,zmax) ! write(6,"('pltlon: zmin,max=',2e12.4)") zmin,zmax endif ! ! Add epflux vectors: if (len_trim(lonsli%vectype) > 0) then yboffset = .28 iclr = iwhite if (icolor>0.and.lonsli%vectype(1:1)=='+') iclr = iblack call pltvec(lonsli%y,lonsli%z,lonsli%nx,lonsli%ny, + lonsli%xx(1),lonsli%xx(lonsli%nx),lonsli%yy(1), + lonsli%yy(lonsli%ny),1,1,yboffset,iclr,labelv, + lonsli%vectype,spval,0) endif ! ! Add perimeter and tic marks: if (lonsli%ihtyax == 0) then call labrect(lonsli%xx,lonsli%nx,lonsli%yy,lonsli%ny, + 'LATITUDE (DEG)','LN(P0/P) ('//trim(lonsli%zptype)//')', | pltchsize) else call labrect(lonsli%xx,lonsli%nx,lonsli%yy,lonsli%ny, + 'LATITUDE (DEG)','HEIGHT (KM)',pltchsize) endif ! ! 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.gt.0) then ! ! Check for nans in zsli: ! subroutine check_nans(f,id1,id2,id3,name,n_total,ispval,spval, ! iprint,ifatal) ! If ispval /= 0, then replace any nans w/ spval. ! Number of nans is returned in nnans_z. ! Do not call yaxzpht if any nans are found. ! ! call check_nans(hts,lonsli%nx,npress,1,'ZSLICE',nnans_z, ! | 0,spval,0,0) ! if (nnans_z > 0) then ! write(6,"('>>> WARNING pltlons: Found ',i8,' NaNs', ! | ' in Z-slice(out of ',i8,')')") ! | nnans_z,lonsli%nx*npress ! else if (lonsli%zptype(1:3)=='MID') then ! midpoints call yaxzpht(gcmlev,hts,gcmlev,lonsli%nx,npress,dlev, | lonsli%yy,lonsli%ny,p0,lonsli%ihtyax,iyaxright,vp,spval) else ! interfaces call yaxzpht(gcmilev,hts,gcmilev,lonsli%nx,npress,dlev, | lonsli%yy,lonsli%ny,p0,lonsli%ihtyax,iyaxright,vp,spval) endif ! endif endif end subroutine pltlon !----------------------------------------------------------------------- subroutine wrout_lons(lonsli,h,hcntr) use plt,only: zmin,zmax use proc use input ! ! Write ascii and/or xdr output data files for current lon slice frame: ! ! Args: ! (lonslice must be intent(inout) because this routine may need ! to call mkloinlabs) ! type(lonslice),intent(inout) :: lonsli 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: ! lonsli%tlabs(2) label is full field name + units ! lonsli%tlabs(3) label is ut and grid info ! lonsli%blabs(2) label is hist vol info ! if (iplot == 0) then call mklonlabs(lonsli,h,hcntr,msgout) endif ! ! 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,lonsli%data,lonsli%nx,lonsli%ny, + lonsli%xx,lonsli%yy,lonsli%xlab,lonsli%ylab,lonsli%tlabs(2), + lonsli%tlabs(3),lonsli%blabs(2),iframe_dat,'tgcmproc',senddat) if (iplot==0) + write(6,"('Data frame ',i4,': ',a)") iframe_dat,trim(msgout) 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,lonsli%data,lonsli%nx,lonsli%ny, + lonsli%xx,lonsli%yy,lonsli%xlab,lonsli%ylab, + lonsli%tlabs(2),lonsli%tlabs(3),lonsli%blabs(2), + lonsli%blabs(1),h%mtime,0) flnm_xdr = flnm_xdr(1:len_trim(flnm_xdr)-1) if (iplot==0) + write(6,"('Xdr frame ',i4,': ',a)") iframe_xdr,trim(msgout) iframe_xdr = iframe_xdr+1 endif return end subroutine wrout_lons !----------------------------------------------------------------------- subroutine wrcdf_lons(lonsli,h,hcntr) ! ! Write netcdf file containing longitude slices. This routine is called ! for each lon slice frame (after sub wrout). Only one lontype (lons or ! slts) is written to the netcdf file. Also, only one y-axis type (zp ! or ht) is written. ! use proc,only: mxslice,spval,zmflag,mxtms,dlon,gcmlon,nlon use input,only: sendcdf_lons,flons,fslts use mkcdf_module,only: | netcdf_file, ! type def | nccreate,ncdefine,ncvalues ! subs include 'netcdf.inc' ! ! Args: type(lonslice),intent(inout) :: lonsli type(history),intent(in) :: h,hcntr ! ! Locals that are saved between calls: integer,save :: ncalls=0,ncid,ihtyax,nflds=1,nslice,mtime(3), | ntime=1 character(len=3),save :: lontype type(netcdf_file),save :: nclon real,save :: slices(mxslice),glons(mxslice,mxtms) ! ! Locals that are not saved between calls: integer :: i,istat,ids1(1),ids2(2),ids4(4),idv,ipos,count4(4), | start4(4),islice,ishape(2),ier,start2(2),count2(2),ixlon character(len=120) :: char120 character(len=80) :: char80 character(len=8) :: fname real,allocatable :: fdata(:,:) real :: dum logical :: newhist ! ! External: real,external :: fslt integer,external :: ixfind ! ! Return silently if not writing netcdf file: if (len_trim(sendcdf_lons) <= 0) return ! ncalls = ncalls+1 ! write(6,"('Enter wrcdf_lons: ncalls=',i3)") ncalls ! write(6,"('wrcdf_lons: lonsli%lontype=',a,' lonsli%glon=',f8.3)") ! | lonsli%lontype,lonsli%glon ! ! First call only: ! if (ncalls == 1) then ! ! Create new netcdf dataset: nclon%filename = flnm_cdf_lons ! local file name call nccreate(nclon%filename,nclon%ncid) ! returns nclon%ncid ncid = nclon%ncid ! ! Policy: only 1 type of vertical axis is written to the netcdf ! file. If both types are requested (both zp and ht), then only the ! type of the first call will be used. ! ihtyax = lonsli%ihtyax ! ht on y-axis if ihtyax==1, zp otherwise ! ! Likewise, only 1 of lons or slts will be written to the netcdf file. ! The type of the first call will be used. ! Determine type and number of slices, and the selected lons or slts: ! select case (lonsli%lontype) case('lon') lontype = 'lon' case('zm ') lontype = 'lon' case('slt') lontype = 'slt' case default write(6,"('>>> wrcdf_lons: unsupported lontype=',a)") | lonsli%lontype lontype = 'lon' end select nslice = 0 if (lontype=='lon') then do i=1,mxslice if (flons(i) /= spval) then nslice = nslice+1 slices(nslice) = flons(i) endif enddo else ! slt do i=1,mxslice if (fslts(i) /= spval) then nslice = nslice+1 slices(nslice) = fslts(i) endif enddo endif ! write(6,"('wrcdf_lons: lontype=',a,' nslice=',i3,' slices=', ! | /,(6f10.2))") lontype,nslice,slices(1:nslice) ! ! First model time: mtime(:) = h%mtime(:) ! ! Dimension var for number of slices (per time): istat = nf_def_dim(ncid,'slices',nslice,nclon%id_nslice) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: error defining nslice dimension.') ! ! Define coord var that will contain values of slices (lons or slts): ids1(1) = nclon%id_nslice istat = nf_def_var(ncid,"slices",NF_FLOAT,1,ids1, | nclon%idv_slice) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining slice coord variable') if (lontype=='lon') then write(char80,"('geographic longitude (-west, +east)')") istat = nf_put_att_text(ncid,nclon%idv_slice,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining long name of lon slice variable') istat = nf_put_att_text(ncid,nclon%idv_slice,"units", | 12,'degrees_east') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining units of lon slice variable') else write(char80,"('solar local time')") istat = nf_put_att_text(ncid,nclon%idv_slice,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining long name of slt slice variable') istat = nf_put_att_text(ncid,nclon%idv_slice,"units", | 5,'hours') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining units of slt slice variable') endif ! ! Define dimensions and coord vars: ! ! Time is the unlimited dimension: istat = nf_def_dim(ncid,"time",NF_UNLIMITED,nclon%id_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining time dimension') ! ! Time (total model time in minutes): ids1(1) = nclon%id_time istat = nf_def_var(ncid,"time",NF_INT,1,ids1,nclon%idv_time) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'ncdefine: Error defining time dimension variable') write(char80,"('total model time (day*24*60+hour*60+min)')") istat = nf_put_att_text(ncid,nclon%idv_time,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'ncdefine: Error defining long_name of time dimension variable') istat = nf_put_att_text(ncid,nclon%idv_time,"units",7,'minutes') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'ncdefine: Error defining units of time dimension variable') ! ! x-axis dimension is always latitude: istat = nf_def_dim(ncid,'lat',lonsli%nx,nclon%id_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: error defining lat (x) dimension') ! ! Latitude coord var: ids1(1) = nclon%id_lat istat = nf_def_var(ncid,"lat",NF_FLOAT,1,ids1,nclon%idv_lat) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining latitude dimension variable') write(char80,"('geographic latitude (-south, +north)')") istat = nf_put_att_text(ncid,nclon%idv_lat,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining long_name of latitude dimension variable') istat = nf_put_att_text(ncid,nclon%idv_lat,"units",13, | 'degrees_north') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining units of latitude dimension variable') ! ! y-axis dimension is either zp or ht: istat = nf_def_dim(ncid,'lev',lonsli%ny,nclon%id_lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: error defining ny dimension') ids1(1) = nclon%id_lev istat = nf_def_var(ncid,"lev",NF_FLOAT,1,ids1,nclon%idv_lev) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining levels dimension variable') if (ihtyax <= 0) then ! zp on y-axis write(char80,"('LOG PRESSURE')") istat = nf_put_att_text(ncid,nclon%idv_lev,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining long_name of zp lev coord variable') istat = nf_put_att_text(ncid,nclon%idv_lev,"units",8, | 'ln(p0/p)') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining units of zp lev coord variable') else ! ht on y-axis write(char80,"('HEIGHT')") istat = nf_put_att_text(ncid,nclon%idv_lev,"long_name", | len_trim(char80),trim(char80)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining long_name of ht lev coord variable') istat = nf_put_att_text(ncid,nclon%idv_lev,"units",2, | 'KM') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error defining units of ht lev coord variable') endif ! ! Define history vars and global attributes: call ncdefine(nclon%ncid,h,nclon) ! ! Define var giving longitudes corresponding to selected slts: ! (this must come after ncdefine because ncdefine sets nclon%id_time) if (lontype=='slt') then ids2(1) = nclon%id_nslice ids2(2) = nclon%id_time istat = nf_def_var(ncid,"longitudes",NF_FLOAT,2,ids2, | nclon%idv_lon) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining var of slt slice lons') write(char120,"('geographic longitude (-west, +east)', | ' corresponding to selected local time slices ', | '(nearest grid point)')") istat = nf_put_att_text(ncid,nclon%idv_lon,"long_name", | len_trim(char120),trim(char120)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining longname of slt slice lons') istat = nf_put_att_text(ncid,nclon%idv_lon,"units", | 12,'degrees_east') if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: Error defining units of slt slice lons') endif ! ! Global attribute giving slice type (lon or slt): istat = nf_put_att_text(ncid,NF_GLOBAL,"slice_type", | len_trim(lontype),trim(lontype)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: error defining slice_type global attribute.') ! ! Global attribute for zonal mean flag: if (lontype=='lon') then istat = nf_put_att_real(ncid,NF_GLOBAL,"zmflag",NF_FLOAT, | 1,zmflag) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'wrcdf_lons: error defining zmflag global attribute.') write(char80,"('Slice value of zmflag refers to zonal means ', | ' (longitudinal average)')") istat = nf_put_att_text(ncid,NF_GLOBAL,"zmflag_explain", | len_trim(char80),trim(char80)) endif nclon%fids(:) = 0 ! init ! ! Exit define mode: istat = nf_enddef(ncid) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error return from nf_enddef') ! ! Give values to coord vars lat, lev, slice: istat = nf_put_var_real(ncid,nclon%idv_lat,lonsli%xx) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error from nf_put_var_real for latitude coord var.') ! istat = nf_put_var_real(ncid,nclon%idv_lev,lonsli%yy) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error from nf_put_var_real for levels coord var.') ! istat = nf_put_var_real(ncid,nclon%idv_slice,slices(1:nslice)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error from nf_put_var_real for slices coord var.') ! ! Update values of history variables for new mtime: call ncvalues(ncid,h,nclon,ntime) ! ! End first call: endif ! ncalls==1 ! ! Return if we are not writing current slice type: if (ncalls > 1) then if ((lontype=='lon'.and.(lonsli%lontype /= 'lon' .and. | lonsli%lontype /= 'zm ')).or. | (lontype=='slt'.and.lonsli%lontype /= 'slt')) then ! write(6,"('wrcdf_lons: skipping this slice because ', ! | 'am writing lontype=',a,' but current type = ',a)") ! | lontype,lonsli%lontype return endif ! ! Re-open file for writing: istat = nf_open(nclon%filename,NF_WRITE,ncid) if (istat /= NF_NOERR) then write(char80,"('Error reopening file ',a)") | trim(nclon%filename) call handle_ncerr(istat,char80) endif ! ! Increment ntime if this is a new history: if (any(mtime(:) /= h%mtime(:))) then ntime = ntime+1 mtime = h%mtime ! ! Update values of history variables for new mtime: call ncvalues(ncid,h,nclon,ntime) newhist = .true. else newhist = .false. endif ! new history endif ! ncalls > 1 ! ! Find index to current slice: islice = 0 sliceloop: do i=1,nslice if ((lontype=='lon'.and.lonsli%glon==slices(i)).or. | (lontype=='slt'.and.lonsli%slt ==slices(i))) then islice = i ! ! If making slt slices, find corresponding longitude at current ut ! (nearest model grid point) for "longitudes" (glons) variable: if (lontype=='slt') then glons(islice,ntime) = fslt(slices(i),h%ut,dum,3) ixlon = ixfind(gcmlon,nlon,glons(islice,ntime),dlon) glons(islice,ntime) = gcmlon(ixlon) endif exit sliceloop endif enddo sliceloop if (islice==0) then if (lontype=='lon') then write(6,"('>>> WARNING wrcdf_lons: cannot find slice index', | ' to lonsli%glon=',e12.4,' slices=',/,(6e12.4))") | lonsli%glon,slices(1:nslice) else write(6,"('>>> WARNING wrcdf_lons: cannot find slice index', | ' to lonsli%slt=',e12.4,' slices=',/,(6e12.4))") | lonsli%slt,slices(1:nslice) endif endif ! ! Define current field, long name and units if it's not on the file yet: fname = lonsli%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,nclon%fids(nflds)) if (istat /= NF_NOERR) then ! error means field was not found istat = nf_redef(ncid) ! put in define mode ! ! Define 4-d field (lat,lev,time,nslice): ids4(1) = nclon%id_lat ids4(2) = nclon%id_lev ids4(3) = nclon%id_nslice ids4(4) = nclon%id_time ! unlimited dim must be in last place istat = nf_def_var(ncid,fname,NF_FLOAT,4,ids4,nclon%fids(nflds)) if (istat /= NF_NOERR) then write(char80,"('wrcdf_lons: Error defining 4d field var ',a)") | trim(fname) call handle_ncerr(istat,char80) else ! write(6,"('wrcdf_lons: defined field ',a,' on netcdf file.', ! | ' nflds =',i3,' id=',i3)") trim(fname),nflds, ! | nclon%fids(nflds) endif ! ! Long name attribute: istat = nf_put_att_text(ncid,nclon%fids(nflds),"long_name", | len_trim(lonsli%fname),trim(lonsli%fname)) if (istat /= NF_NOERR) then write(char80,"('Error defining long name of field ',a)") | trim(lonsli%fname) call handle_ncerr(istat,char80) endif ! ! Units attribute: istat = nf_put_att_text(ncid,nclon%fids(nflds),"units", | len_trim(lonsli%funits),trim(lonsli%funits)) if (istat /= NF_NOERR) then write(char80,"('Error defining units of field ',a)") | trim(lonsli%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 (lonsli%known) then ! known field if (trim(lonsli%zptype)=="MIDPOINTS") then istat = nf_put_att_text(ncid,nclon%fids(nflds),"zp_interp", | 3,'yes') elseif (trim(lonsli%zptype)=="INTERFACES") then istat = nf_put_att_text(ncid,nclon%fids(nflds),"zp_interp", | 2,'no') else write(6,"('>>> mkcdf: unknown lonsli%zptype=',a)") | trim(lonsli%zptype) istat = nf_put_att_text(ncid,nclon%fids(nflds),"zp_interp", | 7,'unknown') endif else ! unknown field istat = nf_put_att_text(ncid,nclon%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,lonsli%known,lonsli%zptype call handle_ncerr(istat,char120) endif ! ! Increment number of fields: nflds = nflds+1 ! ! Exit define mode: istat = nf_enddef(ncid) endif ! define current field ! ! Give values to longitudes of selected slts at current ut: if (lontype=='slt') then start2(1) = 1 start2(2) = ntime count2(1) = nslice count2(2) = 1 istat = nf_put_vara_real(ncid,nclon%idv_lon,start2,count2, | glons(1:nslice,ntime)) if (istat /= NF_NOERR) call handle_ncerr(istat, | 'Error from nf_put_vara_real for lons of slt slices.') endif ! ! Get id of current field: istat = nf_inq_varid(ncid,fname,idv) if (istat /= NF_NOERR) then write(6,"('>>> warning: could not find id of field ',a)") | fname call handle_ncerr(istat,'Error from nf_inq_varid') endif ! ! Remove log10 from data if necessary: ishape = shape(lonsli%data(:,:)) allocate(fdata(ishape(1),ishape(2)),stat=ier) if (ier /= 0) then write(6,"('>>> wrcdf_lons: Error allocating fdata with ishape=', | 2i4)") ishape call allocerr(ier,"wrcdf_lons allocate fdata with ishape") endif if (lonsli%log10 <= 0) then fdata(:,:) = lonsli%data(:,:) else fdata(:,:) = 10.**lonsli%data(:,:) endif ! ! Give values to current field at current slice (lat,lev,slice,time): start4(1:2) = 1 start4(3) = islice start4(4) = ntime count4(1) = lonsli%nx count4(2) = lonsli%ny count4(3:4) = 1 istat = nf_put_vara_real(ncid,idv,start4,count4,fdata) if (istat /= NF_NOERR) then write(char80,"('Error giving values to field ',a,' slice=', | f9.2,'(',a,') ntime=',i3,' mtime=',3i4)") trim(lonsli%sname), | slices(islice),lontype,ntime,mtime call handle_ncerr(istat,char80) endif if (allocated(fdata)) deallocate(fdata) ! ! Close file and report to stdout: istat = nf_close(ncid) if (newhist) write(6,"('Wrote to netcdf file ',a,' at mtime=',i3, | ':',i2,':',i2,' itime=',i5)") trim(flnm_cdf_lons),mtime,ntime end subroutine wrcdf_lons !----------------------------------------------------------------------- subroutine timeav_glbm(timeav,nx,ny,ixlat0,ixlat1,tlabs,blabs, | fmin,fmax,vp) use proc,only: gcmlat,dlon,dlat,spval use input,only: ntimes ! ! Plot global means of time-averaged lon slices. ! ! Args: integer,intent(in) :: nx,ny ! nlat,nlev integer,intent(in) :: ixlat0,ixlat1 type(lonslice),intent(in) :: timeav character(len=240),intent(in) :: tlabs(3),blabs(3) real,intent(in) :: vp(4) real,intent(out) :: fmin,fmax ! ! Local: integer :: k real :: glbm(ny) real,external :: fglbm ! util.F real :: tlabchsz(3) = (/.022,.022,.022/), + blabchsz(3) = (/.020,.020,.020/) real :: toffset=.06, boffset=.15 character(len=240) :: toplabs(3),botlabs(3) ! ! write(6,"('timeav_glbm: nx (lat) =',i3,' ny (lev) =',i3, ! | ' ixlat0,1=',2i4,' data min,max=',2e12.4)") nx,ny,ixlat0,ixlat1, ! | minval(timeav%data),maxval(timeav%data) ! ! Calculate global means of timeav%data ! real function fglbm(f,imx,jmx,gcmlat,dlat,dlon,spval) ! do k=1,ny glbm(k) = fglbm(timeav%data(:,k),1,nx,gcmlat(ixlat0:ixlat1), | dlat,dlon,spval) enddo ! write(6,"('timeav_glbm: glbm=',/,(6e12.4))") glbm ! ! subroutine pltmxy(x,y,npts,nxcurve,multix,xlab,ylab, ! + linelab,chsize,vp,ilog,scalefac,spval) ! call pltmxy(glbm,timeav%yy,ny,1,0,tlabs(2),timeav%ylab, | ' ',0.,vp,timeav%log10,1.,spval) ! ! Remake top and bottom labels for global mean of time-averaged lon slice: ! fmin = minval(glbm) fmax = maxval(glbm) toplabs(1) = tlabs(1) ! usually blank toplabs(2) = tlabs(2) ! diffs and field name toplabs(3) = 'Global Mean of Time-averaged Longitude Slice' write(botlabs(1),"('MIN,MAX=',e12.4,',',e12.4)") fmin,fmax write(botlabs(2),"('AVERAGED OVER ',i4,' TIMES')") ntimes botlabs(3) = ' ' select case (timeav%lontype) case('lon') write(botlabs(3),"('LON SLICE = ',f9.3,' (DEGREES)')") | timeav%glon case('zm ') write(botlabs(3),"('LON SLICE = ZONAL MEANS')") case('slt') write(botlabs(3),"('LOCAL TIME SLICE = ',f6.2,' (HOURS)')") | timeav%slt case default end select call wrlab6(toplabs,toffset,tlabchsz,botlabs,boffset,blabchsz,vp, | 0) end subroutine timeav_glbm !----------------------------------------------------------------------- end module mk_lons