!
      module mk_lats
      use hist,only: history 
      implicit none
!
! A latitude slice is a 2d contour at a selected latitude,
! with longitude on the x-axis, and height or pressure on 
! the y-axis. Lat slice plots can have optional extra bottom 
! axis in slt and extra right axis in ht or zp.
!
      type latslice
        real :: glat		! geog latitude of slice
        integer :: nx,ny	! number of x,y points
        real,pointer :: xx(:)=>NULL(),yy(:)=>NULL()	! x and y coords
        character(len=40) :: xlab,ylab	! x and y axis labels
        real,pointer :: data(:,:)=>NULL()	! the data array
        character(len=56) :: fname	! full field name
        character(len=8)  :: sname	! short field name
        character(len=16) :: funits	! field units
        character(len=16) :: ftype	! field type
        character(len=16) :: zptype	! midpoints or interfaces
        character(len=8)  :: difftype	! raw or percent
        integer :: ihtyax	! ht on yaxis if ihtyax==1, zp otherwise
        integer :: log10	! if > 0, then log10 of field is plotted
        character(len=240) :: tlabs(3),blabs(3)	! top and bottom labels
        logical :: known        ! is field known to the processor
      end type latslice
      real :: tlabchsz(3) = (/.022,.022,.022/),
     +        blabchsz(3) = (/.020,.020,.020/)
      contains
!-------------------------------------------------------------------
      subroutine mklats(f,fcntr,nf,h,hcntr)
!
! Make latitude slice plots and other output:
!
      use proc
      use fields,only: field,z_ifaces,z_midpts,
     |  zcntr_ifaces,zcntr_midpts 
      use input
      use plt
      use set_vert
      implicit none
!
! Args:
      integer,intent(in) :: nf
      type(field),intent(in) :: f(nf),fcntr(nf)
      type(history),intent(in) :: h,hcntr
!
! Locals:
      integer :: i,ixf,ilat,ier,ixlat,izpht0,izpht1,iht,nhtscale,ixz,
     +  iadvfr,k0zp,k1zp,nzprange,k,ix,logint
      type(latslice) :: latsli
      character(len=80) msgout
      real,pointer :: htscale(:)=>NULL() ! user requested y-axes
      real,pointer :: zprange(:)=>NULL() ! user requested y-axes
      real :: vp(4),vplat(4)=(/.10,.90,.30,.90/)! for 0 extra right axes
      real :: toffset=.05, boffset=.36
      real :: zsli(nlon,npress),zsli_cntr(nlon,npress) ! ht slices
      real,allocatable :: tmpsli1(:,:),tmpsli2(:,:)
      real :: flats_tmp
!
! hscale and fhscale are defined by sethscale when field is height-only:
! fhscale(nfhscale) will be height scale of height-only fields.
! hscale(nhscale)   will be height scale to use when plotting.
!
      real,pointer :: hscale(:)=>NULL()	! height scale to be used
      integer :: nhscale		! dimension for hscale
      real,pointer :: fhscale(:)=>NULL()! height scale of field
      integer :: nfhscale		! dimension for fhscale
!
! Externals:
      integer,external :: ixfind,ixfindc
!
! Begin exec:
      write(6,"(' ')")
      if (diffs) write(6,"('Difference fields: ')",advance='NO')
      write(6,"('Latitude slices at ut = ',f5.2,'  mtime = ',3i3)")
     +  h%ut,h%mtime
      if (multiadvfr > 0) nppf = 0
!
! Define x coords in longitude:
      latsli%nx = nlon
      if (.not.associated(latsli%xx)) then
        allocate(latsli%xx(latsli%nx),stat=ier)
        if (ier /= 0) 
     +    call allocerr(ier,"mklats allocating latsli%xx") 
      endif
      latsli%xx = gcmlon	! array op
      write(latsli%xlab,"('LONGITUDE (DEG)')")
!
! Locate index to heights field (needed for interp, etc):
!     ixz = ixfindc(f%fname8,nf,'Z       ')
!     if (ixz <= 0) then
!       write(6,"('>>> mklats: need heights in fields: field names',
!    +    ' = ',(8a8))") f%fname8
!       stop 'mklats z'
!     endif
!
! Field loop. Field must be requested, data pointer must be associated,
!             and field must not be height-independent:
      fields_loop: do ixf=1,nf
        if (.not.f(ixf)%requested.or.
     +      .not.associated(f(ixf)%data).or.
     +      f(ixf)%nlev==1) cycle fields_loop 
!
! Set up vertical scale(s) zp and/or ht:
        if (f(ixf)%zptype(1:3)=='MID') then
          call setvert(flat_zprange,gcmlev,npress,dlev,nzprange,k0zp,
     |      k1zp,zprange,flat_htscale,nhtscale,htscale,spval,1)
        else
          call setvert(flat_zprange,gcmilev,npress,dlev,nzprange,k0zp,
     |      k1zp,zprange,flat_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
!
        latsli%fname  = f(ixf)%fname56
        latsli%sname  = f(ixf)%fname8
        latsli%funits = f(ixf)%units
        latsli%ftype  = f(ixf)%type
        latsli%zptype = f(ixf)%zptype
        latsli%known  = f(ixf)%known
        latsli%difftype = f(ixf)%difftype
        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)
!
! Selected latitude loop:
        lat_loop: do ilat=1,mxslice 
          if (flats(ilat)==spval) cycle lat_loop
          latsli%glat = flats(ilat)
	  flats_tmp = flats(ilat)         ! switch to kind=4 if default
          ixlat = ixfind(gcmlat,nlat,flats_tmp,dlat)
!
! Save heights:
! 1/22/08 btf: use z_ifaces or z_midpts depending on zptype of
!   the current field (see sub mkgeopot in getflds.F): 
!
          if (trim(f(ixf)%zptype) == 'INTERFACES') then
            zsli(:,:) = z_ifaces(:,ixlat,:)
            if (diffs) zsli_cntr(:,:) = zcntr_ifaces(:,ixlat,:)
          else
            zsli(:,:) = z_midpts(:,ixlat,:)
            if (diffs) zsli_cntr(:,:) = zcntr_midpts(:,ixlat,:)
          endif
!
! 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
            latsli%ihtyax = iht
            if (associated(latsli%yy)) deallocate(latsli%yy)
            if (associated(latsli%data)) deallocate(latsli%data)
!
! - - - - - - - - - - - - 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 (ilat_yaxright==0) then
                vplat(2) = .90
              elseif (ilat_yaxright==1) then
                vplat(2) = .85
              else
                vplat(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 pressure:
              latsli%ny = nzprange
              allocate(latsli%yy(latsli%ny),stat=ier)
              if (ier /= 0) 
     +          call allocerr(ier,"mklats allocating latsli%yy") 
              latsli%yy = zprange	! array op
              latsli%ylab = ' '
              latsli%ylab = 'LN(P0/P) '//'('//trim(latsli%zptype)//')'
!
! Allocate and define latsli data for zp on y-axis:
              allocate(latsli%data(latsli%nx,latsli%ny),stat=ier)
              if (ier /= 0) 
     +          call allocerr(ier,"mklats allocating latsli%data") 
              latsli%data = f(ixf)%data(:,ixlat,k0zp:k1zp)
!
! Difference fields:
              if (diffs) then
                allocate(tmpsli1(latsli%nx,latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating tmpsli1")
                tmpsli1 = latsli%data ! perturbed from above
                allocate(tmpsli2(latsli%nx,latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating tmpsli2")
                tmpsli2 = fcntr(ixf)%data(:,ixlat,k0zp:k1zp)
                call mkdiffs(tmpsli1,tmpsli2,latsli%data,
     +            latsli%nx*latsli%ny,f(ixf)%difftype) 
                deallocate(tmpsli1) ; deallocate(tmpsli2)
              endif
!
! - - - - - - - - - - - - ht on y-axis: - - - - - - - - - - - - -
            else
              if (ilat_yaxright==0) then
                vplat(2) = .90
              else
                vplat(2) = .85
              endif
!
! Define y coords in height:
              write(latsli%ylab,"('HEIGHT (KM)')")
              if (trim(f(ixf)%vtype)/='HEIGHT') then
                latsli%ny = nhtscale
                allocate(latsli%yy(latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating latsli%yy") 
                latsli%yy = htscale	! array op
!
! Allocate and define latsli data for height on y-axis:
                allocate(latsli%data(latsli%nx,latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating latsli%data") 
!
! Interpolate to linear height scale:
                call cuthtint(f(ixf)%data(:,ixlat,:),zsli,nlon,npress,
     +            latsli%data,htscale,nhtscale,logint,spval,ier,1)
!
! Difference fields:
                if (diffs) then
                  allocate(tmpsli1(latsli%nx,latsli%ny),stat=ier)
                  if (ier /= 0) 
     +              call allocerr(ier,"mklats allocating tmpsli1")
                  tmpsli1 = latsli%data ! perturbed from above
                  allocate(tmpsli2(latsli%nx,latsli%ny),stat=ier)
                  if (ier /= 0) 
     +              call allocerr(ier,"mklats allocating tmpsli2")
!
! Interpolate control, then take diffs:
                  call cuthtint(fcntr(ixf)%data(:,ixlat,:),zsli_cntr,
     +              nlon,npress,tmpsli2,htscale,nhtscale,logint,spval,
     +              ier,1)
                  call mkdiffs(tmpsli1,tmpsli2,latsli%data,
     +              latsli%nx*latsli%ny,f(ixf)%difftype) 
                  deallocate(tmpsli1) ; deallocate(tmpsli2)
                endif
!
! Field is calculated in height only -- plot only that part of
! user's requested height scale (htscale) that is within height 
! range of the calculated field (fhscale). Also use resolution 
! of the height-only field, overriding requested resolution 
! (currently only oh-v and oh-b fields are height-only):
!
              else
                latsli%ny = nhscale
                allocate(latsli%yy(latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating latsli%yy") 
                latsli%yy = hscale	! array op
                allocate(latsli%data(latsli%nx,latsli%ny),stat=ier)
                if (ier /= 0) 
     +            call allocerr(ier,"mklats allocating latsli%data") 
                do k=1,nhscale
                  if (hscale(k) < fhscale(1).or.
     +                hscale(k) > fhscale(nfhscale)) then
                    latsli%data(:,k) = spval
                  else
                    ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev)
                    latsli%data(:,k) = f(ixf)%data(:,ixlat,ix)
                  endif
                enddo           
!
! Difference fields of height-only field:
                if (diffs) then
                  allocate(tmpsli1(latsli%nx,latsli%ny),stat=ier)
                  if (ier /= 0) 
     +              call allocerr(ier,"mklats allocating tmpsli1")
! Save perturbed from above:
                  tmpsli1 = latsli%data ! perturbed from above
                  allocate(tmpsli2(latsli%nx,latsli%ny),stat=ier)
                  if (ier /= 0) 
     +              call allocerr(ier,"mklats allocating tmpsli2")
! Get control:
                  do k=1,nhscale
                    if (hscale(k) < fhscale(1).or.
     +                  hscale(k) > fhscale(nfhscale)) then
                      tmpsli2(:,k) = spval
                    else
                      ix=ixfind(fhscale,nfhscale,hscale(k),f(ixf)%dlev)
                      tmpsli2(:,k) = fcntr(ixf)%data(:,ixlat,ix)
                    endif
                  enddo           
! Take diffs:
                  call mkdiffs(tmpsli1,tmpsli2,latsli%data,
     +              latsli%nx*latsli%ny,f(ixf)%difftype) 
                  deallocate(tmpsli1) ; deallocate(tmpsli2)
                endif  
              endif	! height-only field
            endif
! - - - - - - - - - - - - end zp or ht on y-axis: - - - - - - - - - 
!
! Take log10 if requested:
            if ((ilat_log10==1.and.trim(f(ixf)%type)=='DENSITY').or.
     +           ilat_log10==2) then
              call log10f(latsli%data,latsli%nx*latsli%ny,1.e-20,spval)
              latsli%log10 = 1
            else
              latsli%log10 = 0
            endif
            if (iplot > 0) then
!
! Set up viewport (for multiplt):
              call setmultivp(vplat,iadvfr,nppf+1,multiplt,ipltrowcol,
     +          vp)
! 
! Contour:
              if (flat_censlt /= spval) then
                call pltlat_rotate_xaxis(latsli,flat_censlt,h%ut)
              endif
              call pltlat(latsli,zsli,ilat_yaxright,vp,h%ut)
!
! Add top and bottom labels:
              call mklatlabs(latsli,h,hcntr,msgout)
              if (len_trim(clat_xaxes(1))==0.or.
     +            len_trim(clat_xaxes(2))==0) boffset=.25
              call wrlab6(latsli%tlabs,toffset,tlabchsz, 
     +                    latsli%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,
     +          'latitude slices',iframe)
!
! If not making plots, wrout will need zmin,zmax.
            else
              call fminmax(latsli%data,size(latsli%data),zmin,zmax)
            endif	! iplot > 0
!
! Write output data files:
            call wrout_lats(latsli,h,hcntr)
          enddo yax_loop
        enddo lat_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,' ','latitude slices',
     +    iframe)
!
! Release space (altho maybe latsli goes away anyway because
!   it's local?)
      if (associated(latsli%data)) deallocate(latsli%data)
      if (associated(latsli%xx))   deallocate(latsli%xx)
      if (associated(latsli%yy))   deallocate(latsli%yy)
      if (associated(htscale))     deallocate(htscale)
      if (associated(zprange))     deallocate(zprange)
      if (associated(fhscale))     deallocate(fhscale)
      if (associated(hscale))      deallocate(hscale)
      return
      end subroutine mklats
!-------------------------------------------------------------------
      subroutine mklatlabs(latsli,h,hcntr,msgout)
      use plt,only: zmin,zmax,ciu,scfac
      use proc,only: spval,diffs
      use input,only: ie5577,ie6300
      implicit none
!
! Construct 3 top labels (latsli%tlabs) and 3 bottom labels 
! (latsli%blabs) for lat slice:
!
! Args:
      type(latslice),intent(inout) :: latsli
      type(history),intent(in) :: h,hcntr
      character(len=*),intent(out) :: msgout
!
! Locals:
      integer :: lenlab
!
! Set up top and bottom text labels:
!
! Top (1st) top label is optionally provided by user:
      latsli%tlabs(1) = ' '
!
! Middle (2nd) top label is full field name + units:
! (prefaced with LOG10 if latsli%log10 > 0)
!
      if (latsli%known) then
        lenlab = len_trim(latsli%fname)+len_trim(latsli%funits)+3
        if (latsli%log10 > 0) lenlab = lenlab+6
        if (lenlab <= len(latsli%tlabs(2))) then
          if (latsli%log10 > 0) then
            latsli%tlabs(2) = 'LOG10 '//trim(latsli%fname)//' ('//
     +      trim(latsli%funits)//')'
          else
            latsli%tlabs(2) = trim(latsli%fname)//' ('//
     +      trim(latsli%funits)//')'
          endif
        else
          latsli%tlabs(2) = trim(latsli%fname)
        endif
      else ! unknown to proc
        lenlab = len_trim(latsli%sname)+2+len_trim(latsli%fname)+
     |           len_trim(latsli%funits)+3
        if (latsli%log10 > 0) lenlab = lenlab+6
        if (lenlab <= len(latsli%tlabs(2))) then
          if (latsli%log10 > 0) then
            if (trim(latsli%sname)/=trim(latsli%fname)) then
              latsli%tlabs(2) = 'LOG10 '//trim(latsli%sname)//': '//
     |          trim(latsli%fname)//' ('//trim(latsli%funits)//')'
            else
              latsli%tlabs(2) = 'LOG10 '//trim(latsli%fname)//
     |          ' ('//trim(latsli%funits)//')'
            endif
          else ! no log10
            if (trim(latsli%sname)/=trim(latsli%fname)) then
              latsli%tlabs(2) = trim(latsli%sname)//': '//
     |          trim(latsli%fname)//' ('//trim(latsli%funits)//')'
            else
              latsli%tlabs(2) = trim(latsli%fname)//
     |          ' ('//trim(latsli%funits)//')'
            endif
          endif
        else
          latsli%tlabs(2) = trim(latsli%fname)
        endif
      endif ! known or unknown
!
! Add DIFFS to label:
      if (diffs) then
        if (trim(latsli%difftype)=='RAW') then
          if (len_trim(latsli%tlabs(2))+7 <= len(latsli%tlabs(2))) then
            latsli%tlabs(2) = 'DIFFS: '//trim(latsli%tlabs(2))
          else
            latsli%tlabs(2) = 'DIFFS: '//trim(latsli%fname)
          endif
        else ! percent
          if (len_trim(latsli%tlabs(2))+15 <= len(latsli%tlabs(2))) then
            latsli%tlabs(2) = trim(latsli%difftype)//' DIFFS: '//
     |        trim(latsli%tlabs(2))
          else
            latsli%tlabs(2) = '% DIFFS: '//trim(latsli%fname)
          endif
        endif
      endif
!
! Bottom (3rd) top label is ut and grid info:
      write(latsli%tlabs(3),"('UT=',f5.2,' LAT=',f9.2,' (DEG)')")
     +  h%ut,latsli%glat
!
! 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(latsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=',
     +    1pe12.4)") zmin,zmax,ciu
      else
        write(latsli%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=',
     +    1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac
      endif
!
! Middle (2nd) bottom label is history info:
! Bottom (3rd) bottom label is also used if doing diffs.
!     if (.not.diffs) then
!       write(latsli%blabs(2),"(a,' ',a,' (DAY,HR,MIN=',
!    +    i3,',',i2,',',i2,')')") trim(h%version),trim(h%mssvol),
!    +    h%mtime 
!       latsli%blabs(3) = ' '
!       if (trim(latsli%ftype)=='EMISSION') then
!         call mkemislab(latsli%sname,ie5577,ie6300,latsli%blabs(3))
!         if (len_trim(latsli%blabs(3)) > 48) blabchsz(3) = .015
!       endif
!     else
!       if (7+len_trim(h%mssvol)    +12<=len(latsli%blabs(2)).and.
!    +      7+len_trim(hcntr%mssvol)+12<=len(latsli%blabs(3)))then
!         write(latsli%blabs(2),"('DIFFS: ',a,' (',i3,',',i2,',',i2,
!    +      ')')") trim(h%mssvol),h%mtime
!         write(latsli%blabs(3),"('MINUS ',a,' (',i3,',',i2,',',i2,
!    +      ')')") trim(hcntr%mssvol),hcntr%mtime
!       else
!         write(latsli%blabs(2),"('DIFFS: (',i3,',',i2,',',i2,
!    +      ') MINUS (',i3,',',i2,',',i2,')')")
!    +      h%mtime,hcntr%mtime
!         latsli%blabs(3) = ' '
!       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(latsli%blabs(2),"(a,' (DAY,HR,MIN=',
     |  i3,',',i2,',',i2,')')") trim(h%version),h%mtime
      if (.not.diffs) then
        write(latsli%blabs(3),"(a)") trim(h%histfile)
      else
        write(latsli%blabs(3),"(a,' MINUS ',a)") 
     |    trim(h%histfile),trim(hcntr%histfile)
      endif
!
! Return message to be printed to stdout:
      write(msgout,"(a,' glat=',f7.2,' min,max=',
     +  2e11.4,' iht=',i2)") latsli%sname,latsli%glat,
     +  zmin,zmax,latsli%ihtyax
!
      return
      end subroutine mklatlabs
!-------------------------------------------------------------------
      subroutine pltlat(latsli,hts,iyax,vp,ut)
      use plt
      use proc,only: gcmlev,gcmilev,npress,dlev,p0,spval
      use input,only: clat_xaxes,flat_censlt
!
! Make latitude slice 2d contour (lon on x-axis, ht or zp on y-axis)
! type(latslice) latsli has been fully defined.
!
! Args:
      type(latslice),intent(in) :: latsli
      real,intent(in) :: hts(:,:),vp(4),ut
      integer,intent(in) :: iyax
!
! Locals:
      integer :: iyaxright,nx,i
      integer,parameter :: idel=60 ! label lons every 60 deg
      integer,parameter :: nilab=360/idel+1
      real :: glonlabs(nilab)
      integer :: ilonlabs(nilab),idum(1)
      real :: dum(1)
!
! Set up conpack:
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',latsli%xx(1))
      call cpsetr('XCM',latsli%xx(latsli%nx))
      call cpsetr('YC1',latsli%yy(1))
      call cpsetr('YCN',latsli%yy(latsli%ny))
      call set(vp(1),vp(2),vp(3),vp(4),
     +  latsli%xx(1),latsli%xx(latsli%nx),
     +  latsli%yy(1),latsli%yy(latsli%ny),1)
!
! Contour:
      call contour(latsli%data,latsli%nx,latsli%nx,latsli%ny)
!
      glonlabs(1) = latsli%xx(1)
      do i=2,nilab
        glonlabs(i) = glonlabs(i-1)+float(idel)
      enddo
      do i=1,nilab
        if (glonlabs(i) >= 180.) glonlabs(i) = glonlabs(i)-360.
        if (glonlabs(i) < -180.) glonlabs(i) = glonlabs(i)+360.
        ilonlabs(i) = int(glonlabs(i))
      enddo
!     write(6,"('pltlat: glonlabs=',/,(8f9.2))") glonlabs
!     write(6,"('pltlat: ilonlabs=',/,(8i8))") ilonlabs
!
! Add perimeter and tic marks:
! clat_xaxes(2) is char*8 from user:
!   if clat_xaxes(1:2) == 'LON','SLT', then draw lon x-axis and 
!     an extra slt axis below the lon axis (this is the default)
!   if clat_xaxes(1:2) == 'LON',' ' , then draw only lon x-axis
!   if clat_xaxes(1:2) == 'SLT' ,' ' , then draw only slt x-axis
!   if clat_xaxes(1:2) == 'SLT' ,'LON' , then this is same as 'LON','LT'
!     (i.e., cannot plot slt on x-axis w/ extra axis in lon)
!
      nx = latsli%nx
      if (trim(clat_xaxes(1))=='LT')
     +  nx = 0	! tell labrect not to plot an x-axis if want slt only
!
      if (latsli%ihtyax == 0) then
        if (flat_censlt == spval) then
          call labrect(latsli%xx,nx,latsli%yy,latsli%ny,
     |      'LONGITUDE (DEG)','LN(P0/P) ('//trim(latsli%zptype)//')',
     |      pltchsize)
        else
          call labrect(latsli%xx,0,latsli%yy,latsli%ny,
     +      'LONGITUDE (DEG)','HEIGHT (KM)',pltchsize)
          call labaxes(nilab,ilonlabs,5,'LONGITUDE (DEG)',
     |      0,idum,2,' ',pltchsize)
        endif
      else
        if (flat_censlt == spval) then
          call labrect(latsli%xx,nx,latsli%yy,latsli%ny,
     +      'LONGITUDE (DEG)','HEIGHT (KM)',pltchsize)
        else
          call labrect(latsli%xx,0,latsli%yy,latsli%ny,
     +      'LONGITUDE (DEG)','HEIGHT (KM)',pltchsize)
          call labaxes(nilab,ilonlabs,5,'LONGITUDE (DEG)',
     |      0,idum,2,' ',pltchsize)
        endif
      endif
!
! X-axis in local time if requested:
      if (trim(clat_xaxes(1))=='LT'.or.
     +    trim(clat_xaxes(2))=='LT') then
        if (nx > 0) then	! lon + slt x-axes
          call sltxax(latsli%xx(1),latsli%xx(latsli%nx),ut,.18)
        else			! slt x-axis only
          call sltxax(latsli%xx(1),latsli%xx(latsli%nx),ut,0.)
        endif
      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
        if (latsli%zptype(1:3)=='MID') then
          call yaxzpht(gcmlev,hts,gcmlev,latsli%nx,npress,dlev,
     +      latsli%yy,latsli%ny,p0,latsli%ihtyax,iyaxright,vp,spval)
        else
          call yaxzpht(gcmilev,hts,gcmilev,latsli%nx,npress,dlev,
     +      latsli%yy,latsli%ny,p0,latsli%ihtyax,iyaxright,vp,spval)
        endif
      endif
      end subroutine pltlat
!-------------------------------------------------------------------
      subroutine wrout_lats(latsli,h,hcntr)
      use plt,only: zmin,zmax
      use proc
      use input
!
! Write ascii and/or xdr output data files for current lat slice frame:
!
! Args:
! (latslice must be intent(inout) because this routine may need
!  to call mklatlabs)
!
      type(latslice),intent(inout) :: latsli
      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:
!   latsli%tlabs(2) label is full field name + units
!   latsli%tlabs(3) label is ut and grid info
!   latsli%blabs(2) label is hist vol info
!
      if (iplot == 0) call mklatlabs(latsli,h,hcntr,msgout)
!
! Write to ascii data file:
!     subroutine wrdat(iwr,lu,f,nx,ny,xx,yy,xlab,ylab,fieldlab,
!    +  infolab,histlab,iframe,proclab,senddat)
!
      if (iwrdat > 0) then
        call wrdat(iwrdat,ludat,latsli%data,latsli%nx,latsli%ny,
     +    latsli%xx,latsli%yy,latsli%xlab,latsli%ylab,latsli%tlabs(2),
     +    latsli%tlabs(3),latsli%blabs(2),iframe_dat,'tgcmproc',senddat)
        if (iplot==0) then
          write(6,"('Data frame ',i4,': ',a)") iframe_dat,trim(msgout)
        endif
        iframe_dat = iframe_dat+1
      endif
!
! Write to xdr data file:
!     subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1,
!    +  lab2,lab3,lab4,mtime,iclose)
!
      if (iwrxdr > 0) then
        call wrxdr(flnm_xdr,latsli%data,latsli%nx,latsli%ny,
     +    latsli%xx,latsli%yy,latsli%xlab,latsli%ylab,
     +    latsli%tlabs(2),latsli%tlabs(3),latsli%blabs(2),
     +    latsli%blabs(1),h%mtime,0)
        flnm_xdr = flnm_xdr(1:len_trim(flnm_xdr)-1)
        if (iplot==0) then
          write(6,"('Xdr frame ',i4,': ',a)") iframe_xdr,trim(msgout)
        endif
        iframe_xdr = iframe_xdr+1
      endif
      return
      end subroutine wrout_lats
!-------------------------------------------------------------------
      subroutine pltlat_rotate_xaxis(latsli,censlt,ut)
!
! Args:
      type(latslice),intent(inout) :: latsli
      real,intent(in) :: censlt,ut
!
! Local:
      integer :: ilon0,ilon1,nx,icen,k
      real :: xx(latsli%nx),dum,glon0,glon1,dlon,slt0,slt1
      integer,save :: npts
!
! External:
      real,external :: fslt
!
! integer function ixfind8(z,nz,val,del)
      integer,external :: ixfind
!
      nx = latsli%nx
      icen = nx/2+1
      dlon = latsli%xx(2)-latsli%xx(1)   ! assume monotonically increasing
!
! Rotate x-axis:
!
      glon0 = latsli%xx(icen)
      slt0 = fslt(dum,ut,glon0,1)
      ilon0 = ixfind(latsli%xx,nx,glon0,dlon)
      glon1 = fslt(censlt,ut,dum,3)       ! lon at censlt
!
! Check if xx rotation is necessary (xx is not reset inside field loop): 
      if (latsli%xx(icen)==glon1.or.
     |   (abs(latsli%xx(icen))==180..and.abs(glon1)==180.)) then
        xx = latsli%xx
      else
        ilon1 = ixfind(latsli%xx,nx,glon1,dlon)
        if (ilon1==1) ilon1 = nx
        npts = ilon1-ilon0
        call rotate_periodic(latsli%xx,xx,nx,npts)
        latsli%xx = xx
      endif ! rotate x-axis
!
! Rotate the field (always need to do this):
      do k=1,latsli%ny
        call rotate_periodic(latsli%data(:,k),xx,nx,npts)
        latsli%data(:,k) = xx(:)
      enddo
      end subroutine pltlat_rotate_xaxis
!-------------------------------------------------------------------
      end module mk_lats