module mk_zmslt use hist,only: history implicit none type zmslt_type real,pointer :: | xx(:)=>NULL(),yy(:)=>NULL(), ! x and y coords (x is lat, y is zp or ht) | 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 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 integer :: log10 ! if > 0, then log10 of field is plotted integer :: mtime0(3),mtime1(3) real :: slt ! user-requested local time character(len=1024) :: + tlabs(3),blabs(3), ! top and bottom labels + hvol0,hvol1, ! 1st and last hist vols + hvol0_cntr,hvol1_cntr ! 1st and last hist vols (control) character(len=8) :: difftype ! raw or percent logical :: known ! true if field is known to the processor end type zmslt_type ! ! zmslt is zmslt_ type to be plotted type(zmslt_type),save :: zmslt ! ! zmsltdat(nlat,npress,nslt,nf) is data from which zmslt is set up. ! real,allocatable,save :: zmsltdat(:,:,:,:) ! (nlat,npress,nslt,nf) real,allocatable,save :: zmsltdat_cntr(:,:,:,:) ! (nlat,npress,nslt,nf) ! contains !------------------------------------------------------------------- subroutine mkzmslt(f,fcntr,nf,h,hcntr,itime) ! ! Make latitude vs zp/ht contours of averaged longitudes corresponding ! to selected local time (namelist zmslt) across a 1-day time series. ! use proc use fields,only: field,nohalt,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use input,only: zmsltinp=>zmslt,zmslt_zprange,zmslt_htscale, | ipltzmslt,izmslt_yaxright,multiplt,ipltrowcol,izmslt_log10, | iplot,ilab_hq,senddat 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 ! ! Local: integer :: j,k,ixf,noh,ier,ixlon,ifld,k0zp,k1zp,nzprange, | nhtscale,logint,izpht0,izpht1,iht,iadvfr,isltax,islt integer,save :: nfld,ixz_zmslt,nslt real :: ut,glon,fmin,fmax,vp(4) real,pointer :: | htscale(:)=>NULL(), | zprange(:)=>NULL() real :: vp_zmslt(4) = (/.12,.85,.29,.90/) real,allocatable :: fdat1(:,:),fdat2(:,:) real,save :: inv_ntimes character(len=1024) msgout real :: toffset=.06, boffset=.15 real :: tlabchsz(3) = (/.022,.022,.022/), + blabchsz(3) = (/.020,.020,.020/) ! ! External: real,external :: fslt integer,external :: ixfind,ixfindc ! write(6,"('mkzmslt: itime=',i4,' ntimes=',i4,' h%mtime=',3i4)") ! | itime,ntimes,h%mtime ! ! Allocate zmsltdat if this is first call: if (itime==1) then inv_ntimes = 1./real(ntimes) ixz_zmslt = 0 ! index to pert hts in zmsltdat nfld = 0 ! number of fields to allocate noh = 0 ! number of oh-v/b fields to allocate 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_zmslt = nfld else noh = noh+1 endif endif enddo if (nfld==0.and.noh==0) then write(6,"('>>> mkzmslt: need fields: nfld==noh==0.', + ' Turning ipltzmslt off.')") ipltzmslt = 0 return endif ! ! Number of requested local times: nslt = count(zmsltinp /= spval) if (nslt == 0) then write(6,"('>>> mkzmslt: no local times were requested,', | ' setting default zmslt=12.')") zmsltinp(1) = 12. nslt = 1 endif ! ! if z not already requested, add extra field in zmsltdat to save ! total heights for later interpolation: if (ixz_zmslt==0) then nfld = nfld+1 ixz_zmslt = nfld endif ! ! Do the allocation: if (nfld > 0) then allocate(zmsltdat(nlat,npress,nslt,nfld),stat=ier) if (ier /= 0) then write(6,"('mkzmslt error allocating', | ' zmsltdat(nlat,npress,nfld): nlat=',i3,' npress=',i3, | ' nslt=',i3,' nfld=',i3)") nlat,npress,nslt,nfld call allocerr(ier,"mkzmslt allocating zmsltdat") else write(6,"('mkzmslt: allocated zmsltdat: nlat=',i4, | ' npress=',i4,' nslt=',i3,' nfld=',i4)") | nlat,npress,nslt,nfld endif zmsltdat = 0. ! init if (diffs) then allocate(zmsltdat_cntr(nlat,npress,nslt,nfld),stat=ier) if (ier /= 0) then write(6,"('mkzmslt error allocating ', | 'mkzmslt_cntr(nlat,npress,nslt,nfld): nlat=',i4, | ' npress=',i4,' nslt=',i3,' nfld=',i4)") | nlat,npress,nslt,nfld call allocerr(ier,"mkzmslt allocating zmsltdat_cntr") else write(6,"('mkzmslt: allocated zmsltdat_cntr: nlat=',i4, | ' npress=',i4,' nslt=',i3,' nfld=',i4)") | nlat,npress,nslt,nfld endif endif ! diffs zmsltdat_cntr = 0. ! init zmslt%mtime0 = h%mtime endif ! nflds > 0 ! ! Set up zmslt x-axis: zmslt%nx = nlat allocate(zmslt%xx(zmslt%nx),stat=ier) if (ier/=0) + call allocerr(ier,"mkzmslt allocating zmslt%xx") zmslt%xx = gcmlat zmslt%xlab = 'Latitude' zmslt%hvol0 = h%histfile if (diffs) zmslt%hvol0_cntr = hcntr%histfile elseif (itime==ntimes) then zmslt%hvol1 = h%histfile if (diffs) zmslt%hvol1_cntr = hcntr%histfile zmslt%mtime1 = h%mtime endif ! itime==1 ! ! Field loop: ifld = 0 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)=='HT-INDEP') | cycle fields_loop if (trim(f(ixf)%vtype)/='HEIGHT') then ifld = ifld+1 ! field index to zmsltdat(:,:,:,ifld) else write(6,"('>>> mkzmslt will not plot height-only field',a)") | f(ixf)%fname8 cycle fields_loop endif logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 ! ! If this is last time, define parts of zmslt and prepare for contouring: if (itime==ntimes) then ! ! Set up vertical scale(s) zp and/or ht: if (f(ixf)%zptype(1:3)=='INT') then ! interfaces call setvert(zmslt_zprange,gcmilev,npress,dlev,nzprange, + k0zp,k1zp,zprange,zmslt_htscale,nhtscale,htscale,spval,1) else ! midpoints call setvert(zmslt_zprange,gcmlev,npress,dlev,nzprange, + k0zp,k1zp,zprange,zmslt_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 zmslt%hvol1 = h%histfile if (diffs) zmslt%hvol1_cntr = hcntr%histfile ! zmslt%fname = f(ixf)%fname56 zmslt%sname = f(ixf)%fname8 zmslt%funits = f(ixf)%units zmslt%ftype = f(ixf)%type zmslt%zptype = f(ixf)%zptype zmslt%difftype = f(ixf)%difftype zmslt%known = f(ixf)%known ! Set field min, max, and interval to be used in contouring: ! (pltmin,pltmax,conint,scfac are in plt.F) pltmin = f(ixf)%cmin pltmax = f(ixf)%cmax conint = f(ixf)%cint scfac = f(ixf)%scalefac endif ! itime==ntimes (still in fields loop) ! ! Local time loop: slt_loop: do islt = 1,nslt ! ! Calculate running average of fields(j,k) at longitudes corresponding ! to requested slt at this ut: ut = real(h%mtime(2)) glon = fslt(zmsltinp(islt),ut,glon,3) ! longitude for zmslt at current ut ixlon = ixfind(gcmlon,nlon,glon,dlon) zmslt%slt = zmsltinp(islt) ! write(6,"('slt_loop: islt=',i3,' nslt=',i3,' slt=',f6.2, ! | ' glon=',f9.2,' ixlon=',i4)") islt,nslt,zmslt%slt,glon,ixlon ! ! Define zmsltdat for current time: do k=1,npress do j=1,zmslt%nx if (f(ixf)%data(ixlon,j,k) /= spval) then zmsltdat(j,k,islt,ifld) = zmsltdat(j,k,islt,ifld)+ | f(ixf)%data(ixlon,j,k)*inv_ntimes else zmsltdat(j,k,islt,ifld) = spval endif enddo enddo ! call fminmax(zmsltdat(:,:,islt,ifld),zmslt%nx*npress, ! | fmin,fmax) ! write(6,"('Define zmsltdat for itime=',i4,' slt=',f6.2, ! | ' ixf=',i4,' ifld=',i4,' nfld=',i4,' field ',a,' min,max=', ! | 2e12.4)") itime,zmslt%slt,ixf,ifld,nfld, ! | f(ixf)%fname8,fmin,fmax if (diffs) then do k=1,npress do j=1,zmslt%nx if (fcntr(ixf)%data(ixlon,j,k) /= spval) then zmsltdat_cntr(j,k,islt,ifld) = | zmsltdat_cntr(j,k,islt,ifld)+ | fcntr(ixf)%data(ixlon,j,k)*inv_ntimes else zmsltdat_cntr(j,k,islt,ifld) = spval endif enddo enddo ! call fminmax(zmsltdat_cntr(:,:,islt,ifld), ! | zmslt%nx*npress,fmin,fmax) ! write(6,"('Define zmsltdat for itime=',i4,' slt=',f6.2, ! | ' ixf=',i4,' ifld=',i4,' nfld=',i4,' field ',a,' min,max=' ! | ,2e12.4)") itime,zmslt%slt,ixf,ifld,nfld, ! | f(ixf)%fname8,fmin,fmax endif ! diffs ! ! If not ready to plot, skip vertical axes loop: if (itime /= ntimes) cycle slt_loop ! ! Loop for vertical axes zp and/or ht: ! (This loop is executed only when itime==ntimes) yax_loop: do iht=izpht0,izpht1 zmslt%ihtyax = iht if (associated(zmslt%yy)) deallocate(zmslt%yy) if (associated(zmslt%data)) deallocate(zmslt%data) ! ! - - - - - - - - - - - - zp on y-axis: - - - - - - - - - - - - - - if (iht == 0) then ! zp on y-axis if (trim(f(ixf)%vtype)=='HEIGHT') cycle yax_loop if (izmslt_yaxright > 0) vp_zmslt(2) = .77 ! ! Define y coords in pressure: zmslt%ny = nzprange allocate(zmslt%yy(zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating zmslt%yy") zmslt%yy = zprange ! array op write(zmslt%ylab,"('LN(P0/P) (',a,')')") | trim(zmslt%zptype) ! ! Allocate and define zmslt data for zp on y-axis: allocate(zmslt%data(zmslt%nx,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating zmslt%data") if (.not.diffs) then zmslt%data = zmsltdat(:,k0zp:k1zp,islt,ifld) else allocate(fdat1(zmslt%nx,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating fdat1") fdat1 = zmsltdat(:,k0zp:k1zp,islt,ifld) allocate(fdat2(zmslt%nx,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating fdat2") fdat2 = zmsltdat_cntr(:,k0zp:k1zp,islt,ifld) call mkdiffs(fdat1,fdat2,zmslt%data,zmslt%nx*zmslt%ny, + f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif ! ! - - - - - - - - - - - - ht on y-axis: - - - - - - - - - - - - - - else vp_zmslt(2) = .86 ! ! Define y coords in height: write(zmslt%ylab,"('HEIGHT (KM)')") if (trim(f(ixf)%vtype)/='HEIGHT') then zmslt%ny = nhtscale allocate(zmslt%yy(zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating zmslt%yy") zmslt%yy = htscale ! array op ! ! Allocate and define zmslt data for height on y-axis: allocate(zmslt%data(zmslt%nx,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating zmslt%data") ! ! Interpolate to linear height scale: if (.not.diffs) then call cuthtint(zmsltdat(:,:,islt,ifld), + zmsltdat(:,:,islt,ixz_zmslt),nlat,npress, + zmslt%data,htscale,nhtscale,logint,spval,ier,1) else allocate(fdat1(nlat,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating fdat1") call cuthtint(zmsltdat(:,:,islt,ifld), + 0.5*(zmsltdat(:,:,islt,ixz_zmslt)+ + zmsltdat_cntr(:,:,islt,ixz_zmslt)), + nlat,npress,fdat1,htscale,nhtscale,logint,spval, + ier,1) allocate(fdat2(nlat,zmslt%ny),stat=ier) if (ier /= 0) + call allocerr(ier,"mkzmslt allocating fdat2") call cuthtint(zmsltdat_cntr(:,:,islt,ifld), + 0.5*(zmsltdat(:,:,islt,ixz_zmslt)+ + zmsltdat_cntr(:,:,islt,ixz_zmslt)), + nlat,npress,fdat2,htscale,nhtscale,logint,spval, + ier,1) call mkdiffs(fdat1,fdat2,zmslt%data,nlat*zmslt%ny, + f(ixf)%difftype) deallocate(fdat1) ; deallocate(fdat2) endif ! diffs endif ! not a height-only field endif ! zp or ht on y-axis ! - - - - - - - - - - - - end zp or ht on y-axis: - - - - - - - - - ! ! Take log10 if requested: if ((izmslt_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + izmslt_log10==2) then call log10f(zmslt%data,zmslt%nx*zmslt%ny,1.e-20,spval) zmslt%log10 = 1 else zmslt%log10 = 0 endif if (iplot > 0) then call setmultivp(vp_zmslt,iadvfr,nppf+1,multiplt, + ipltrowcol,vp) ! ! Contour: isltax = 0 call pltzmslt(zmslt,zmsltdat(:,:,islt,ixz_zmslt), + nlat,npress,izmslt_yaxright,vp,h%ut) ! ! Add top and bottom labels: call mkzmsltlabs(zmslt,h,hcntr,msgout) call wrlab6(zmslt%tlabs,toffset,tlabchsz, + zmslt%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, + 'lat vs vertical at zmslt',iframe) endif ! ! From lonslice: ! 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 ascii data file: if (iwrdat > 0) then call wrdat(iwrdat,ludat,zmslt%data,zmslt%nx,zmslt%ny, | zmslt%xx,zmslt%yy,zmslt%xlab,zmslt%ylab,zmslt%tlabs(2), | zmslt%tlabs(3),zmslt%blabs(2),iframe_dat,'tgcmproc', | senddat) write(6,"('Data frame ',i4,': ',a)") | iframe_dat,trim(msgout) iframe_dat = iframe_dat+1 endif enddo yax_loop enddo slt_loop enddo fields_loop end subroutine mkzmslt !------------------------------------------------------------------- subroutine pltzmslt(zmslt,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(zmslt_type) zmslt and type(history) h have been fully defined. ! ! Args: type(zmslt_type),intent(in) :: zmslt 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,fmin,fmax integer :: nnans_z ! number of nans ! ! Set up conpack: call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',zmslt%xx(1)) call cpsetr('XCM',zmslt%xx(zmslt%nx)) call cpsetr('YC1',zmslt%yy(1)) call cpsetr('YCN',zmslt%yy(zmslt%ny)) xleft = zmslt%xx(1)-(zmslt%xx(2)-zmslt%xx(1))/2. xright = zmslt%xx(zmslt%nx)+(zmslt%xx(2)-zmslt%xx(1))/2. call set(vp(1),vp(2),vp(3),vp(4),xleft,xright, + zmslt%yy(1),zmslt%yy(zmslt%ny),1) ! ! Contour: ! fmin = minval(zmslt%data) ; fmax = maxval(zmslt%data) ! write(6,"('pltzmslt call contour: zmslt%data min,max=',2e12.4)") ! | fmin,fmax call contour(zmslt%data,zmslt%nx,zmslt%nx,zmslt%ny) ! ! Add perimeter and tic marks: if (zmslt%ihtyax == 0) then call labrect(zmslt%xx,zmslt%nx,zmslt%yy,zmslt%ny, + 'LATITUDE (DEG)','LN(P0/P) ('//trim(zmslt%zptype)//')', | pltchsize) else call labrect(zmslt%xx,zmslt%nx,zmslt%yy,zmslt%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,zmslt%nx,npress,1,'ZSLICE',nnans_z, ! | 0,spval,0,0) ! if (nnans_z > 0) then ! write(6,"('>>> WARNING pltzmslt: Found ',i8,' NaNs', ! | ' in Z-slice(out of ',i8,')')") ! | nnans_z,zmslt%nx*npress ! else if (zmslt%zptype(1:3)=='MID') then ! midpoints call yaxzpht(gcmlev,hts,gcmlev,zmslt%nx,npress,dlev, | zmslt%yy,zmslt%ny,p0,zmslt%ihtyax,iyaxright,vp,spval) else ! interfaces call yaxzpht(gcmilev,hts,gcmilev,zmslt%nx,npress,dlev, | zmslt%yy,zmslt%ny,p0,zmslt%ihtyax,iyaxright,vp,spval) endif ! endif endif end subroutine pltzmslt !------------------------------------------------------------------- subroutine mkzmsltlabs(zmslt,h,hcntr,msgout) use input,only: izmslt_log10,ntimes,diffs use plt,only: zmin,zmax,ciu,scfac ! ! Args: type(zmslt_type),intent(inout) :: zmslt type(history),intent(in) :: h,hcntr character(len=*),intent(out) :: msgout ! ! Top (1st) top label is optionally provided by user: zmslt%tlabs(:) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if zmslt%log10 > 0) ! zmslt%tlabs(:) = ' ' if (zmslt%log10 > 0) then write(zmslt%tlabs(2),"('LOG10 ',a,' (',a,')')") | trim(zmslt%fname),trim(zmslt%funits) else write(zmslt%tlabs(2),"(a,' (',a,')')") | trim(zmslt%fname),trim(zmslt%funits) endif if (diffs) then if (trim(zmslt%difftype)=='RAW') then zmslt%tlabs(2) = 'DIFFS: '//trim(zmslt%tlabs(2)) else zmslt%tlabs(2) = '%DIFFS: '//trim(zmslt%tlabs(2)) endif endif write(zmslt%tlabs(3),"('Zonal Mean at constant SLT = ',f6.2)") | zmslt%slt ! ! 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(zmslt%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(zmslt%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac endif write(zmslt%blabs(2),"('Zonal Mean over ',i4,' times: ', | i3,',',i2,',',i2,' to ',i3,',',i2,',',i2)") | ntimes,zmslt%mtime0,zmslt%mtime1 if (trim(zmslt%hvol0)==trim(zmslt%hvol1)) then write(zmslt%blabs(3),"('File: ',a)") trim(zmslt%hvol0) else write(zmslt%blabs(3),"(a,' to ',a)") trim(zmslt%hvol0), | trim(zmslt%hvol1) endif if (diffs) then write(zmslt%blabs(3),"(a,' minus ',a)") trim(zmslt%hvol0), | trim(zmslt%hvol0_cntr) endif ! ! For stdout: if (zmslt%ylab(1:2)=='LN') then ! write(msgout,"('Zonal Mean at constant Local Time: Field ',a, ! | ' SLT=',f6.2,' Y-axis in Zp:',f7.2,' to ',f7.2, ! | ' (ln(p0/p))')") zmslt%sname,zmslt%slt,zmslt%yy(1), ! | zmslt%yy(zmslt%ny) ! add fmin,fmax to stdout msg: write(msgout,"('Zonal Mean at constant Local Time: Field ',a, | ' SLT=',f6.2,' Y-axis in Zp, min,max=',2e12.4)") | zmslt%sname,zmslt%slt,zmin,zmax else ! write(msgout,"('Zonal Mean at constant Local Time: Field ',a, ! | ' SLT=',f6.2,' Y-axis in Ht:',f9.2,' to ',f9.2,' (km)')") ! | zmslt%sname,zmslt%slt,zmslt%yy(1),zmslt%yy(zmslt%ny) write(msgout,"('Zonal Mean at constant Local Time: Field ',a, | ' SLT=',f6.2,' Y-axis in Ht, min,max=',2e12.4)") | zmslt%sname,zmslt%slt,zmin,zmax endif end subroutine mkzmsltlabs !------------------------------------------------------------------- end module mk_zmslt