! module mk_utlat use hist,only: history implicit none ! ! Make ut vs latitude contours at selected vertical and longitude. ! Selected vertical can be zp or height, and selected longitude can ! be geographic longitude, local time, or zonal mean. ! type utlat_type integer,pointer :: mtimes(:,:)=>NULL() ! model times (ut on x-axis) real,pointer :: + xx(:)=>NULL(),yy(:)=>NULL(), ! x and y coords (x is ut, y is latitude) + xslt(:)=>NULL(), ! local times corresponding to ut's + data(:,:)=>NULL() ! data(nx,ny) to be contoured integer :: nx,ny ! number of x and y coordinates real :: zpht,glon,slt ! vert, lon, and slt of current data real :: pltlat0,pltlat1 ! min,max latitudes to plot character(len=3) :: lontype ! 'lon', 'slt', or 'zm ' character(len=5) :: vtype ! 'zp', 'ht', 'indep', or 'integ' character(len=8) :: sname ! field short name character(len=16) :: funits ! field units character(len=16) :: ftype ! field type character(len=40) :: xlab,ylab ! x and y axis labels character(len=56) :: fname ! full field name character(len=8) :: difftype ! raw or percent integer :: log10 ! if > 0, then log10 of field is plotted character(len=160) :: + 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) logical :: known ! is field known to the processor end type utlat_type real :: tlabchsz(3) = (/.02,.02,.02/), + blabchsz(3) = (/.018,.018,.018/) integer :: issltmap ! ! utlat is utlat_type to be plotted type(utlat_type),save :: utlat ! ! utlatdat(ntimes,nlat,nzphtlon,nf) is data from which utlat is set up. ! utlatdat(itime,:,:,:) is defined each time sub mkutlat is called. ! The last time mkutlat is called (itime==ntimes), utlat is defined ! from utlatdat, and contours are made. ! real,allocatable :: utlatdat(:,:,:,:) contains !------------------------------------------------------------------- subroutine mkutlat(f,fcntr,nf,h,hcntr,itime) ! ! Set up ut on x-axis, latitude on y-axis, at selected zpht,lon. ! (selected longitudes can include slt or zm) ! (contouring is done by pltutlat, in this module) ! This routine is called at every model time (index itime), and the ! data is saved in utlatdat(itime,:,:,:). ! Contours are made (from utlatdat) only if the current call is for the ! last time (itime==ntimes) ! use proc use fields,only: field,nohalt,oh_alt,z_ifaces,z_midpts, | zcntr_ifaces,zcntr_midpts use input use plt ! ! Args: integer,intent(in) :: nf,itime type(field),intent(in) :: f(nf),fcntr(nf) type(history),intent(in) :: h,hcntr ! ! Locals: integer :: izphtlon,nfld,ifld,ixf,nfemis, + i,j,k,n,ier,l,ixlon,izp,logint,iadvfr,isltax integer,save :: ixz,nzphtlon,nzphtlon_nointeg,nsel_lon real :: slt,dum,fmin,fmax,vp(4),boff,flat1(nlat),flat2(nlat) real,allocatable :: flon(:,:),zlon(:,:) real,allocatable :: flon_cntr(:,:),zlon_cntr(:,:) real :: vputlat(4) = (/.15,.90,.30,.90/), donelons(10) real :: toffset=.06, boffset=.35 real :: utlat_tmp ! character(len=80) msgout character(len=512) msgout logical :: isemis ! ! Externals: real,external :: fslt,fmean,quadrat integer,external :: ixfind,ixfindc,nunique_r logical,external :: isslt ! if (itime==ntimes) then if (.not.diffs) then write(6,"(/'Ut vs latitude contours at selected zpht/lon:')") else write(6,"(/'DIFFERENCE FIELDS: Ut vs latitude contours at ', + 'selected zpht/lon.')") endif if (multiadvfr > 0) nppf = 0 ! utlat%hvol1 = h%mssvol ! if (diffs) utlat%hvol1_cntr = hcntr%mssvol utlat%hvol1 = h%histfile if (diffs) utlat%hvol1_cntr = hcntr%histfile endif ! ! Allocate utlatdat if this is first call: ! Get number of valid zpht/lons and number of unique selected lons/zm: if (itime==1) then nzphtlon = 0 ! total number of slices for dimensioning utlatdat do i=1,mxslice if (utlat_zphtlon(1,i)/=spval.and.utlat_zphtlon(2,i)/=spval) | nzphtlon = nzphtlon+1 enddo nsel_lon = nunique_r(utlat_zphtlon(2,:),mxslice,spval) ! ! Add zphtlon slots for vertical integrations: ! (emission fields are integrated in height, once per selected longitude) nfemis = 0 ! number of requested emission fields do ixf=1,nf if (f(ixf)%requested.and.associated(f(ixf)%data).and. + (trim(f(ixf)%type)=='EMISSION'.or. + trim(f(ixf)%type)=='OH-BAND'.or. | trim(f(ixf)%fname8)=='XLO1')) nfemis = nfemis+1 enddo nzphtlon_nointeg = nzphtlon nzphtlon = nzphtlon+nsel_lon*nfemis if (nzphtlon==0) then write(6,"('>>> mkutlat: need selected zp/ht and lons', + ' (utlat_zphtlon). Turning off ipltutlat.')") ipltutlat = 0 return endif ! ! Get number of fields (also increment nzphtlon for integrations): nfld = 0 ! number of fields to allocate do ixf=1,nf if (f(ixf)%requested.and.associated(f(ixf)%data)) then nfld = nfld+1 endif enddo if (nfld==0) then write(6,"('>>> mkutlat: need fields: nfld==0.', + ' Turning ipltutlat off.')") ipltutlat = 0 return endif ! ! Do the allocation: allocate(utlatdat(ntimes,nlat,nzphtlon,nfld),stat=ier) if (ier /= 0) then write(6,"('mkutlat error allocating', + ' utlatdat(ntimes,nlat,nzphtlon,nfld): ntimes=', + i4,' nlat=',i3,' nzphtlon=',i3,' nfld=',i3)") + ntimes,nlat,nzphtlon,nfld call allocerr(ier,"mkutlat allocating utlatdat") else write(6,"('mkutlat: allocated utlatdat(ntimes,nlat,', + 'nzphtlon,nfld), where:',/,' ntimes=',i4,' nlat=',i2, + ' nzphtlon=',i2,' nfld=',i2)") ntimes,nlat,nzphtlon, + nfld endif ! ! Locate index to heights field (needed for interp, etc): ! (if diffs, pert hts were saved at ixz and control hts at ixzcntr) ! (note locals ixz and ixzcntr are saved across all calls) ! ixz = ixfindc(f%fname8,nf,'Z ') ! if (ixz <= 0) then ! write(6,"('>>> mkutlat: need heights in fields: field ', ! + 'names = ',(8a8))") f%fname8 ! stop 'mkutlat z' ! endif ! ! Set up utlat x-axis (utlat%xx(itime) and utlat%mtimes(:,itime) are ! defined at each call below): utlat%nx = ntimes allocate(utlat%xx(utlat%nx),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating utlat%xx") utlat%xlab = 'UT (HRS)' allocate(utlat%mtimes(3,ntimes),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating utlat%mtimes") ! utlat%hvol0 = h%mssvol ! if (diffs) utlat%hvol0_cntr = hcntr%mssvol utlat%hvol0 = h%histfile if (diffs) utlat%hvol0_cntr = hcntr%histfile ! ! Latitude min,max to plot on y-axis: ! utlat_latminmax(2) is from namelist read in input.F ! These are applied in the conpack calls in sub pltutlat. ! if (utlat_latminmax(1)==spval) then utlat%pltlat0 = gcmlat(1) else utlat%pltlat0 = utlat_latminmax(1) endif if (utlat_latminmax(2)==spval) then utlat%pltlat1 = gcmlat(nlat) else utlat%pltlat1 = utlat_latminmax(2) endif if (utlat%pltlat0 < gcmlat(1).or. | utlat%pltlat0 > gcmlat(nlat)) then write(6,"('>>> mkutlat: bad utlat%pltlat0=',f9.3)") | utlat%pltlat0 write(6,"('Please check namelist read utlat_latminmax(1)')") stop 'utlat%pltlat0' endif if (utlat%pltlat1 < gcmlat(1).or. | utlat%pltlat1 > gcmlat(nlat)) then write(6,"('>>> mkutlat: bad utlat%pltlat1=',f9.3)") | utlat%pltlat1 write(6,"('Please check namelist read utlat_latminmax(2)')") stop 'utlat%pltlat1' endif if (utlat%pltlat0 >= utlat%pltlat1) then write(6,"('>>> mkutlat: pltlat0 must be < pltlat1: ', | ' pltlat0=',f9.3,' pltlat1=',f9.3)") utlat%pltlat0, | utlat%pltlat1 write(6,"('Please check namelist read utlat_latminmax')") stop 'utlat%pltlat' endif write(6,"('mkutlat: utlat%pltlat0,1=',2f9.3)") utlat%pltlat0, | utlat%pltlat1 endif ! itime==1 ! ! Set up utlat y-axis if last time: if (itime==ntimes) then utlat%ny = nlat utlat%ylab = 'LATITUDE (DEG)' allocate(utlat%yy(utlat%ny),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating utlat%yy") utlat%yy = gcmlat endif ! ! Define current time in utlat: utlat%xx(itime) = float(h%mtime(1))*24.+h%ut utlat%mtimes(:,itime) = h%mtime ! use pert for now ! ! Field loop: ifld = 0 fields_loop: do ixf=1,nf if (.not.f(ixf)%requested.or..not.associated(f(ixf)%data)) + cycle fields_loop ifld = ifld+1 ! field index to utlatdat(:,:,:,ifld) isemis = .false. if (trim(f(ixf)%type)=='EMISSION'.or. + trim(f(ixf)%type)=='OH-BAND'.or. | trim(f(ixf)%fname8)=='XLO1') isemis = .true. ! ! If this is last time, define parts of utlat and prepare for contouring: if (itime==ntimes) then utlat%fname = f(ixf)%fname56 utlat%sname = f(ixf)%fname8 utlat%funits = f(ixf)%units utlat%ftype = f(ixf)%type utlat%difftype = f(ixf)%difftype utlat%known = f(ixf)%known ! ! Set field min, max, and interval to be used in contouring: ! (default f(i)%cmin,cmax,cint is 0.,0.,0., but may be provided ! by user via namelist input fmnmxint) ! (pltmin,pltmax,conint are in plt.f) pltmin = f(ixf)%cmin pltmax = f(ixf)%cmax conint = f(ixf)%cint scfac = f(ixf)%scalefac endif ! itime==ntimes in fields loop ! ! Selected zphtlon loop: izphtlon = 0 ! loc index to utlatdat(:,:,izphtlon,:) donelons = spval ! array op zphtlon_loop: do l=1,mxslice if (utlat_zphtlon(1,l)==spval.or.utlat_zphtlon(2,l)==spval) + cycle zphtlon_loop izphtlon = izphtlon+1 ! ! Zonal means: if (utlat_zphtlon(2,l) == zmflag) then ! zonal means utlat%lontype = 'zm ' utlat%glon = zmflag ! ! Regular longitude or local time: else if (isslt(utlat_zphtlon(2,l),slt)) then ! is local time utlat%lontype = 'slt' utlat%glon = fslt(slt,h%ut,dum,3) ixlon = ixfind(gcmlon,nlon,utlat%glon,dlon) ! write(6,"('mkutlat: slt=',f9.4,' utlat%glon=',f9.4, ! | ' ixlon=',i3,' gcmlon(ixlon)=',f9.4)") slt,utlat%glon, ! | ixlon,gcmlon(ixlon) ! utlat%glon = gcmlon(ixlon) utlat%slt = slt else ! is regular longitude utlat%lontype = 'lon' utlat%glon = utlat_zphtlon(2,l) ixlon = ixfind(gcmlon,nlon,utlat%glon,dlon) utlat%slt = fslt(dum,h%ut,utlat%glon,1) endif endif ! ! utlat%vtype can be zp, ht, indep, or integ (but is integ only ! after emission -- see goto 100) if (f(ixf)%nlev==1) then utlat%vtype = 'indep' elseif (utlat_zphtlon(1,l) <= f(ixf)%lev(npress)) then utlat%vtype = 'zp ' elseif (utlat_zphtlon(1,l) /= spval) then utlat%vtype = 'ht ' endif ! ! Save total heights zlon (if diffs, also save control heights): ! 1/23/08 btf: use Z consistent w/ zptype (midpoints or interfaces) ! if (allocated(zlon)) deallocate(zlon) allocate(zlon(nlat,npress),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating zlon") ! ! Field is on interfaces -> use Z on interfaces: if (f(ixf)%zptype == 'INTERFACES') then if (utlat%lontype /= 'zm ') then ! no diffs, no zm ! ! 5/31/07 btf: interpolate to utlat%glon: call interp_lonslice(nlon,nlat,npress,z_ifaces, | utlat%glon,gcmlon,zlon) else ! no diffs, zm do k=1,npress do j=1,nlat zlon(j,k) = fmean(z_ifaces(:,j,k),nlon,spval,0) enddo enddo endif if (diffs) then if (allocated(zlon_cntr)) deallocate(zlon_cntr) allocate(zlon_cntr(nlat,npress),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating zlon_cntr") if (utlat%lontype /= 'zm ') then ! diffs, no zm zlon_cntr(:,:) = zcntr_ifaces(ixlon,:,:) else ! diffs, zm do k=1,npress do j=1,nlat zlon_cntr(j,k) = fmean(zcntr_ifaces(:,j,k),nlon, + spval,0) enddo enddo endif endif ! ! Field is on midpoints -> use Z on midpoints: elseif (f(ixf)%zptype == 'MIDPOINTS') then if (utlat%lontype /= 'zm ') then ! no diffs, no zm ! ! 5/31/07 btf: interpolate to utlat%glon: call interp_lonslice(nlon,nlat,npress,z_midpts, | utlat%glon,gcmlon,zlon) else ! no diffs, zm do k=1,npress do j=1,nlat zlon(j,k) = fmean(z_midpts(:,j,k),nlon,spval,0) enddo enddo endif if (diffs) then if (allocated(zlon_cntr)) deallocate(zlon_cntr) allocate(zlon_cntr(nlat,npress),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating zlon_cntr") if (utlat%lontype /= 'zm ') then ! diffs, no zm zlon_cntr(:,:) = zcntr_midpts(ixlon,:,:) else ! diffs, zm do k=1,npress do j=1,nlat zlon_cntr(j,k) = fmean(zcntr_midpts(:,j,k),nlon, + spval,0) enddo enddo endif endif else ! write(6,"('>>> mkutlat: unknown zptype=',a)") f(ixf)%zptype endif ! ! Calculate field slice (if diffs, also save control field): if (allocated(flon)) deallocate(flon) allocate(flon(nlat,f(ixf)%nlev),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating flon") if (utlat%lontype /= 'zm ') then ! no zm ! ! 5/31/07 btf: interpolate to utlat%glon: call interp_lonslice(nlon,nlat,f(ixf)%nlev,f(ixf)%data, | utlat%glon,gcmlon,flon) else ! zm do k=1,f(ixf)%nlev do j=1,nlat flon(j,k) = fmean(f(ixf)%data(:,j,k),nlon,spval,0) enddo enddo endif if (diffs) then if (allocated(flon_cntr)) deallocate(flon_cntr) allocate(flon_cntr(nlat,f(ixf)%nlev),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating flon_cntr") if (utlat%lontype /= 'zm ') then ! no zm flon_cntr(:,:) = fcntr(ixf)%data(ixlon,:,:) else ! zm do k=1,f(ixf)%nlev do j=1,nlat flon_cntr(j,k) = + fmean(fcntr(ixf)%data(:,j,k),nlon,spval,0) enddo enddo endif endif 100 continue ! ! Define utlatdat(itime,:,izphtlon,ifld) (where 2nd dim is nlat): ! ! Ht-independent field: select case (utlat%vtype) case ('indep') if (.not.diffs) then utlatdat(itime,:,izphtlon,ifld) = flon(:,1) else call mkdiffs(flon(:,1),flon_cntr(:,1), + utlatdat(itime,:,izphtlon,ifld),nlat,f(ixf)%difftype) endif ! ! At selected zp: case ('zp ') if (trim(f(ixf)%vtype)=='HEIGHT') cycle zphtlon_loop utlat_tmp = utlat_zphtlon(1,l) ! switch to kind=4 if default izp = ixfind(f(ixf)%lev,npress,utlat_tmp,dlev) if (izp <= 0) then write(6,"('>>> utlat: bad selected zp', + ' (utlat_zphtlon(1,l)=',f10.2)") utlat_zphtlon(1,l) cycle zphtlon_loop endif utlat%zpht = f(ixf)%lev(izp) if (.not.diffs) then utlatdat(itime,:,izphtlon,ifld) = flon(:,izp) else call mkdiffs(flon(:,izp),flon_cntr(:,izp), + utlatdat(itime,:,izphtlon,ifld),nlat,f(ixf)%difftype) endif ! ! At selected ht: case ('ht ') ! ! Do not process if doing height-only field (e.g. oh-v/b) and ! selected height is outside vertical range of the field: if (trim(f(ixf)%vtype)=='HEIGHT'.and. + (utlat_zphtlon(1,l) < oh_alt(1) .or. + utlat_zphtlon(1,l) > oh_alt(nohalt))) then write(6,"('Note: height ',f10.2,' is outside ', + 'vertical resolution of field ',a)") + utlat_zphtlon(1,l),f(ixf)%fname8 cycle zphtlon_loop endif utlat%zpht = utlat_zphtlon(1,l) ! ! Height-only field (oh-v/b): if (trim(f(ixf)%vtype)=='HEIGHT') then k = ixfind(oh_alt,nohalt,utlat%zpht,1.) if (k <= 0) write(6,"('>>> mkutlat: cannot find', + ' utlat%zpht=',f10.2,' in oh_alt.')") utlat%zpht if (.not.diffs) then utlatdat(itime,:,izphtlon,ifld) = flon(:,k) else call mkdiffs(flon(:,k),flon_cntr(:,k), + utlatdat(itime,:,izphtlon,ifld),nlat, + f(ixf)%difftype) endif ! ! Interpolate to selected height: else logint = 0 if (trim(f(ixf)%type)=='DENSITY') logint = 1 if (.not.diffs) then ! dum = utlat_zphtlon(1,l) ! switch to kind=4 if default call cuthtint(flon,zlon,nlat,f(ixf)%nlev,flat1, + dum,1,logint,spval,ier,1) utlatdat(itime,:,izphtlon,ifld) = flat1(:) ! ! Diffs at selected height: else dum = utlat_zphtlon(1,l) ! switch to kind=4 if default call cuthtint(flon,zlon,nlat,f(ixf)%nlev,flat1, + dum,1,logint,spval,ier,0) call cuthtint(flon_cntr,zlon_cntr,nlat,f(ixf)%nlev, + flat2,dum,1,logint,spval,ier,0) call mkdiffs(flat1,flat2, + utlatdat(itime,:,izphtlon,ifld),nlat,f(ixf)%difftype) endif endif ! ! ht-integrated: case ('integ') do j=1,nlat if (trim(f(ixf)%vtype)/='HEIGHT') then utlatdat(itime,j,izphtlon,ifld) = + quadrat(flon(j,:),zlon(j,:),f(ixf)%nlev) else ! ht-only field (f(ixf)%nlev==nohalt) utlatdat(itime,j,izphtlon,ifld) = + quadrat(flon(j,:),oh_alt,f(ixf)%nlev) endif enddo ! ! Diffs of ht-integrations: if (diffs) then flat1 = utlatdat(itime,:,izphtlon,ifld) ! save pert do j=1,nlat if (trim(f(ixf)%vtype)/='HEIGHT') then flat2(j) = quadrat(flon_cntr(j,:),zlon(j,:), + f(ixf)%nlev) else ! ht-only field (f(ixf)%nlev==nohalt) flat2(j) = quadrat(flon_cntr(j,:),oh_alt,f(ixf)%nlev) endif enddo call mkdiffs(flat1,flat2,utlatdat(itime,:,izphtlon,ifld), + nlat,f(ixf)%difftype) endif case default write(6,"('>>> mkutlat: unknown utlat%vtype=',a)") + utlat%vtype end select ! ! Make plots and/or output data: if (itime==ntimes) then if (.not.associated(utlat%data)) then allocate(utlat%data(utlat%nx,utlat%ny),stat=ier) if (ier/=0) + call allocerr(ier,"mkutlat allocating utlat%data") endif utlat%data = utlatdat(:,:,izphtlon,ifld) ! ! Set ht-integrated units: if (utlat%vtype=='integ') then if (trim(f(ixf)%type)=='OH-BAND') then if (ibohv_watts <= 0) then utlat%data = utlat%data * 1.e-9 ! kilo-rayleighs utlat%funits = "KILO-RAYLEIGHS" else utlat%data = utlat%data / (4.*pi) ! watts/cm2-str utlat%funits = "WATTS/CM2-STR" endif else ! emission field if (trim(f(ixf)%fname8)=='XLO1') then ! O1 recombination rate utlat%funits = "RECOMB/CM2/S" else utlat%data = utlat%data * 1.e-6 ! rayleighs utlat%funits = "RAYLEIGHS" endif endif endif ! ! Take log10 if requested: if ((iutlat_log10==1.and.trim(f(ixf)%type)=='DENSITY').or. + iutlat_log10==2) then call log10f(utlat%data,utlat%nx*utlat%ny,1.e-20,spval) utlat%log10 = 1 else utlat%log10 = 0 endif if (iplot > 0) then call setmultivp(vputlat,iadvfr,nppf+1,multiplt, + ipltrowcol,vp) isltax = 0 if(utlat%lontype=='lon') isltax = 1 call pltutlat(utlat,isltax,vp) ! ! Add top and bottom labels: call mkutlatlabs(h,hcntr,msgout) boff = boffset if (isltax <= 0) boff=.20 call wrlab6(utlat%tlabs,toffset,tlabchsz, + utlat%blabs,boff ,blabchsz,vp,ilab_hq) ! ! Advance frame and write info to stdout: nppf = nppf+1 call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,iadvfr,nppf,msgout, + 'ut vs latitude',iframe) ! ! Make sltmaps if doing a local time plot: issltmap = 0 if (iutlat_sltmaps > 0.and.utlat%lontype=='slt') then issltmap = 1 call plt_sltmaps(ifld,izphtlon) endif ! ! If not making plots, wrout will need zmin,zmax. else call fminmax(utlat%data,size(utlat%data),zmin,zmax) endif ! iplot > 0 ! ! Write output data files: call wrout_utlat(h,hcntr) endif ! itime==ntimes ! ! Do height integration of emission field at current longitude/zm: ! (do this only once per selected longitude, after zpht's have been ! done at current lon) ! if (isemis .and. utlat%vtype /= 'integ' .and. | .not.any(donelons==utlat_zphtlon(2,l)).or. | (isemis.and.utlat%vtype /= 'integ'.and. | izphtlon==nzphtlon_nointeg.and.nsel_lon==1)) then utlat%vtype = 'integ' izphtlon = izphtlon+1 do i=1,mxslice if (donelons(i)==spval) then donelons(i) = utlat_zphtlon(2,l) exit endif enddo goto 100 endif ! ! Release slices: deallocate(flon) deallocate(zlon) if (allocated(flon_cntr)) deallocate(flon_cntr) if (allocated(zlon_cntr)) deallocate(zlon_cntr) enddo zphtlon_loop enddo fields_loop ! if (itime==ntimes) then ! ! Release allocated memory: if (allocated(utlatdat)) deallocate(utlatdat) if (associated(utlat%xx)) deallocate(utlat%xx) if (associated(utlat%yy)) deallocate(utlat%yy) if (associated(utlat%data)) deallocate(utlat%data) endif return contains !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine plt_sltmaps(ifld,izphtlon) ! internal to mkutlat use mk_maps,only: mapproj,inc ! ! Args: integer,intent(in) :: ifld,izphtlon ! ! Locals: integer :: i,ier,ixlon,ih real :: slt,fmin,fmax,dum,utime,utime1,deltalon type(mapproj) :: map real :: vp(4) real :: vpce(4) = (/.100,.960,.450,.880/) real :: vpst(4) = (/.1625,.8375,.195,.870/) real :: vpsv(4) = (/.150,.850,.190,.890/) real :: vpmo(4) = (/.060,.960,.250,.800/) integer,parameter :: nproj=4 ! glb,pol,sat,moll real :: vptmp(4,nproj) integer :: iproj(4) character(len=2) :: projs(nproj) ! ! These maps are made only when a full day, no more no less, is ! being processed: if (utlat%xx(utlat%nx)-utlat%xx(1) /= 24.) then write(6,"('>>> plt_sltmaps: maps are made only when one day', + ' is processed.')") write(6,"(' nx=',i2,' Increasing ut (including model day)', + ' = ',/(6f9.2))") utlat%nx,utlat%xx write(6,"('I will attempt to resolve this by setting the ', | 'last ut (',f9.2,') equal to the first ut + 24 (',f9.2,')')") | utlat%xx(utlat%nx),utlat%xx(1)+24. utlat%xx(utlat%nx) = utlat%xx(1)+24. endif ! ! Msg to stdout: ! write(6,"(/'Multi-ut maps at fixed solar local times: ', ! | ' ifld=',i3,' izphtlon=',i3)") ifld,izphtlon ! ! Calculate longitude x-axis array from mtimes and the fixed slt: ! (use nearest grid point longitudes) ! map%nx = ntimes allocate(map%xx(map%nx),stat=ier) if (ier/=0) + call allocerr(ier,"plt_sltmaps allocating map%xx") utime = utlat%mtimes(2,1)+utlat%mtimes(3,1)/60. map%xx(1) = fslt(slt,utime,dum,3) ixlon = ixfind(gcmlon,nlon,map%xx(1),dlon) map%xx(1) = gcmlon(ixlon) deltalon = 360./float(ntimes-1) ! ! Give ezmap increasing longitudes: do i=2,ntimes map%xx(i) = map%xx(i-1)+deltalon enddo map%xlab = 'LONGITUDE (DEG)' map%ylab = 'LATITUDE (DEG)' ! ! Full latitude on y-axis: map%ny = nlat allocate(map%yy(map%ny),stat=ier) if (ier/=0) + call allocerr(ier,"plt_sltmaps allocating map%yy") map%yy = gcmlat map%cenlat = 0. map%pltlat0 = map%yy(1) map%pltlat1 = map%yy(map%ny) ! ! Set data component of map structure: ! Assume increasing ut which means decreasing longitude. Because ! ezmap got increasing longitudes (see above), we need to reverse ! order of data in longitude: ! Note that unit changes may have been made to utlat%data by ! sub mkutlat. ! allocate(map%data(map%nx,map%ny),stat=ier) if (ier/=0) + call allocerr(ier,"plt_sltmaps allocating map%data") do i=1,map%nx map%data(i,:) = utlat%data(map%nx-i+1,:) enddo ! ! Set remaining map components (most of these are temporary): ! (map%u,v not yet allocated) map%icontinents = map_continents map%eradii = fmap_satview_eradii map%aveht = 0. ! dummy map%vectype = ' ' ! ! Make map projections and contours: ! if (iplot > 0) then iproj(1) = map_global ; projs(1) = 'CE' ; vptmp(:,1)=vpce iproj(2) = map_polar ; projs(2) = 'ST' ; vptmp(:,2)=vpst iproj(3) = map_satview ; projs(3) = 'SV' ; vptmp(:,3)=vpsv iproj(4) = map_mollweide ; projs(4) = 'MO' ; vptmp(:,4)=vpmo blabchsz = (/.017,.017,.017/) tlabchsz = (/.016,.016,.016/) proj_loop: do i=1,nproj if (iproj(i) <= 0) cycle proj_loop call setmultivp(vptmp(:,i),iadvfr,nppf+1,multiplt, + ipltrowcol,vp) map%proj = projs(i) boff = .22 call mkutlatlabs(h,hcntr,msgout) select case (i) ! ! Global (CE): case (1) if (map_global_cenlon /= ispval) then map%cenlon = float(map_global_cenlon) elseif (map_global_censlt /= ispval) then map%cenlon = fslt(float(map_global_censlt),h%ut, + dum,3) else map%cenlon = 0. endif msgout = trim(msgout)//' (global CE map)' call pltce(map,vp,spval) ! ! Polar (ST): case (2) do j=1,mxperimlat if (fmap_polar_perimlat(j) /= spval) then map%perimlat = fmap_polar_perimlat(j) call pltst(map,vp,h%ut,gcmlat,inc,nlat,spval) if (j < mxperimlat.and.fmap_polar_perimlat(j+1)/= | spval) then call mkutlatlabs(h,hcntr,msgout) msgout = trim(msgout)//' (polar maps)' call wrlab6(utlat%tlabs,toffset,tlabchsz, | utlat%blabs,boff ,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, | 'Multi-ut maps at fixed local time)',iframe) endif endif enddo ! Sat View (SV): case (3) if (fmap_satview_latlon(1)/=spval.and. + fmap_satview_latlon(2)/=spval) then map%cenlat = fmap_satview_latlon(1) map%cenlon = fmap_satview_latlon(2) elseif (fmap_satview_latslt(1)/=spval.and. + fmap_satview_latslt(2)/=spval) then map%cenlat = fmap_satview_latslt(1) map%cenlon = fslt(fmap_satview_latslt(2),h%ut,dum,3) else map%cenlat = 40. map%cenlon = 0. endif msgout = trim(msgout)//' (satview map)' call pltsv(map,vp,inc,nlat,spval) boff = .12 ! ! Mollweide (MO): case (4) if (fmap_mollweide_latlon(1)/=spval.and. + fmap_mollweide_latlon(2)/=spval) then map%cenlat = fmap_mollweide_latlon(1) map%cenlon = fmap_mollweide_latlon(2) elseif (fmap_mollweide_latslt(1)/=spval.and. + fmap_mollweide_latslt(2)/=spval) then map%cenlat = fmap_mollweide_latslt(1) map%cenlon = fslt(fmap_mollweide_latslt(2),h%ut, + dum,3) else map%cenlat = 0. map%cenlon = 0. endif msgout = trim(msgout)//' (mollweide map)' call pltmo(map,vp,spval) boff = .15 end select call wrlab6(utlat%tlabs,toffset,tlabchsz, + utlat%blabs,boff ,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, + 'Multi-ut maps at fixed local time)',iframe) enddo proj_loop if (multiplt > 0 .and. multiadvfr == 1 .and. iadvfr <= 0) + call advframe(iwk_cgm,igks_cgm,iwk_ps,igks_ps, + iwk_x11,igks_x11,multiplt,1,nppf,' ', + 'Multi-ut maps at fixed local time',iframe) endif ! iplot > 0 ! ! Write ascii output data files for sltmaps: if (iwrdat > 0) then ! ! Write actual longitudes (nearest gridpoint) corresponding to the ! fixed local time at each ut for x-axis values on the data file: do i=1,ntimes utime = utlat%mtimes(2,i)+utlat%mtimes(3,i)/60. map%xx(i) = fslt(slt,utime,dum,3) ixlon = ixfind(gcmlon,nlon,map%xx(i),dlon) map%xx(i) = gcmlon(ixlon) ! map%data(i,:) = utlatdat(i,:,izphtlon,ifld) map%data(i,:) = utlat%data(i,:) enddo call wrdat(iwrdat,ludat,map%data,map%nx,map%ny, | map%xx,map%yy,map%xlab,map%ylab,utlat%tlabs(2), | utlat%tlabs(3),utlat%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 ! ! Release local space: deallocate(map%data) deallocate(map%xx) deallocate(map%yy) end subroutine plt_sltmaps end subroutine mkutlat !------------------------------------------------------------------- subroutine mkutlatlabs(h,hcntr,msgout) use plt,only: zmin,zmax,ciu,scfac use proc,only: spval,diffs use input,only: ie5577,ie6300 ! ! Construct 3 top labels (utlat%tlabs) and 3 bottom labels ! (utlatv%blabs) for ut vs lat plots: ! ! Args: 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: utlat%tlabs(1) = ' ' ! ! Middle (2nd) top label is full field name + units: ! (prefaced with LOG10 if utlat%log10 > 0) ! (prefaced with DIFFS if doing diffs) ! if (utlat%known) then lenlab = len_trim(utlat%fname)+len_trim(utlat%funits)+3 if (utlat%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(utlat%tlabs(2))) then if (utlat%log10 > 0) then utlat%tlabs(2) = 'LOG10 '//trim(utlat%fname)//' ('// + trim(utlat%funits)//')' else utlat%tlabs(2) = trim(utlat%fname)//' ('// + trim(utlat%funits)//')' endif else utlat%tlabs(2) = trim(utlat%fname) endif else ! unknown to proc lenlab = len_trim(utlat%sname)+2+len_trim(utlat%fname)+ | len_trim(utlat%funits)+3 if (utlat%log10 > 0) lenlab = lenlab+6 if (lenlab <= len(utlat%tlabs(2))) then if (utlat%log10 > 0) then if (trim(utlat%sname)/=trim(utlat%fname)) then utlat%tlabs(2) = 'LOG10 '//trim(utlat%sname)//': '// | trim(utlat%fname)//' ('//trim(utlat%funits)//')' else utlat%tlabs(2) = 'LOG10 '//trim(utlat%fname)// | ' ('//trim(utlat%funits)//')' endif else ! no log10 if (trim(utlat%sname)/=trim(utlat%fname)) then utlat%tlabs(2) = trim(utlat%sname)//': '// | trim(utlat%fname)//' ('//trim(utlat%funits)//')' else utlat%tlabs(2) = trim(utlat%fname)// | ' ('//trim(utlat%funits)//')' endif endif else utlat%tlabs(2) = trim(utlat%fname) endif endif ! known or unknown ! ! Add DIFFS to label: if (diffs) then if (trim(utlat%difftype)=='RAW') then if (len_trim(utlat%tlabs(2))+7 <= len(utlat%tlabs(2))) then utlat%tlabs(2) = 'DIFFS: '//trim(utlat%tlabs(2)) else utlat%tlabs(2) = 'DIFFS: '//trim(utlat%fname) endif else ! percent if (len_trim(utlat%tlabs(2))+15 <= len(utlat%tlabs(2))) then utlat%tlabs(2) = trim(utlat%difftype)//' DIFFS: '// | trim(utlat%tlabs(2)) else utlat%tlabs(2) = '% DIFFS: '//trim(utlat%fname) endif endif endif ! ! Bottom (3rd) top label is grid info: if (utlat%lontype /= 'zm ') then if (utlat%lontype == 'lon') then ! longitude select case (utlat%vtype) case ('zp ') ! lon at zp write(utlat%tlabs(3),"(' ZP = ',f8.3,' LONGITUDE = ',f9.2)") + utlat%zpht,utlat%glon case ('ht ') ! lon at ht write(utlat%tlabs(3),"(' HEIGHT = ',f6.1,' LONGITUDE = ', + f9.2)") utlat%zpht,utlat%glon case ('indep') ! lon at ht-indep write(utlat%tlabs(3),"(' LONGITUDE = ',f9.2)") utlat%glon case ('integ') ! lon at ht-integ write(utlat%tlabs(3),"(' HEIGHT-INTEGRATED: LONGITUDE = ', + f9.2)") utlat%glon case default end select else ! local time select case (utlat%vtype) case ('zp ') ! slt at zp write(utlat%tlabs(3),"(' ZP = ',f8.3,' LOCAL TIME = ', + f5.2)") utlat%zpht,utlat%slt case ('ht ') ! slt at ht write(utlat%tlabs(3),"(' HEIGHT = ',f6.1,' LOCAL TIME = ', + f5.2)") utlat%zpht,utlat%slt case ('indep') ! slt at ht-indep write(utlat%tlabs(3),"(' LOCAL TIME = ',f5.2)") utlat%slt case ('integ') ! slt at ht-integ write(utlat%tlabs(3),"(' HEIGHT INTEGRATED: LOCAL TIME = ', + f5.2)") utlat%slt case default end select endif else ! zonal means select case (utlat%vtype) case ('zp ') ! zm at zp write(utlat%tlabs(3),"(' ZP = ',f8.3,' ZONAL MEANS')") + utlat%zpht case ('ht ') ! zm at ht write(utlat%tlabs(3),"(' HEIGHT = ',f6.1,' ZONAL MEANS')") + utlat%zpht case ('indep') ! zm at ht-indep write(utlat%tlabs(3),"(' ZONAL MEANS')") case ('integ') ! zm at ht-integ write(utlat%tlabs(3),"(' HEIGHT-INTEGRATED ZONAL MEANS')") case default end select endif ! ! Top (1st) bottom label is min,max,interval: ! (zmin,zmax, and ciu are in module plt and were defined by contour) ! (scfac is in module plt) if (scfac==1.) then write(utlat%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4)") zmin,zmax,ciu else write(utlat%blabs(1),"('MIN,MAX=',2(1pe12.4),' INTERVAL=', + 1pe12.4,' (X',1pe9.2,')')") zmin,zmax,ciu,scfac endif ! ! Middle (2nd) and bottom (3rd) bottom label are history info: ! if (.not.diffs) then ! if (len_trim(utlat%hvol0)+len_trim(utlat%hvol1)+15 <= ! + len(utlat%blabs(2))) then ! write(utlat%blabs(2),"('FIRST, LAST: ',a,', ',a)") ! + trim(utlat%hvol0), trim(utlat%hvol1) ! else ! write(utlat%blabs(2),"('FIRST: ',a)") trim(utlat%hvol0) ! endif ! utlat%blabs(3) = ' ' ! if (trim(utlat%ftype)=='EMISSION') then ! call mkemislab(utlat%sname,ie5577,ie6300,utlat%blabs(3)) ! if (len_trim(utlat%blabs(3)) > 48) blabchsz(3) = .015 ! endif ! else ! if (len_trim(utlat%hvol0)+len_trim(utlat%hvol1)+15 <= ! + len(utlat%blabs(2)).and. ! + len_trim(utlat%hvol0_cntr)+len_trim(utlat%hvol1_cntr)+15 <= ! + len(utlat%blabs(3))) then ! write(utlat%blabs(2),"('FIRST, LAST: ',a,', ',a,' MINUS ')") ! + trim(utlat%hvol0), trim(utlat%hvol1) ! write(utlat%blabs(3),"('FIRST, LAST: ',a,', ',a)") ! + trim(utlat%hvol0_cntr), trim(utlat%hvol1_cntr) ! else ! write(utlat%blabs(2),"('FIRST: ',a,' MINUS ')") ! + trim(utlat%hvol0) ! write(utlat%blabs(3),"('FIRST: ',a)") ! + trim(utlat%hvol0_cntr) ! endif ! endif ! ! 5/05 btf: Put first and last history files in blabs(2 and 3): if (.not.diffs) then if (trim(utlat%hvol0) /= trim(utlat%hvol1)) then write(utlat%blabs(2),"('FIRST: ',a)") trim(utlat%hvol0) write(utlat%blabs(3),"('LAST: ',a)") trim(utlat%hvol1) else write(utlat%blabs(2),"('MODEL ',a)") trim(h%version) write(utlat%blabs(3),"(a)") trim(utlat%hvol0) endif else ! diffs: show perturbed and control files: if (trim(utlat%hvol0) /= trim(utlat%hvol1)) then write(utlat%blabs(2),"('PERT FIRST: ',a,' LAST: ',a)") | trim(utlat%hvol0),trim(utlat%hvol1) write(utlat%blabs(3),"('CNTR FIRST: ',a,' LAST: ',a)") | trim(utlat%hvol0_cntr),trim(utlat%hvol1_cntr) else write(utlat%blabs(2),"('PERT FILE: ',a)") trim(utlat%hvol0) write(utlat%blabs(3),"('CNTR FILE: ',a)") | trim(utlat%hvol0_cntr) endif endif ! ! Return message to be printed to stdout: select case (utlat%lontype) case ('zm ') select case (utlat%vtype) case ('zp ') ! zm at zp write(msgout,"(a,' zp=',f6.1,' zonal means min,max=', + 2e11.4)") utlat%sname,utlat%zpht,zmin,zmax case ('ht ') ! zm at ht write(msgout,"(a,' ht=',f6.1,' zonal means min,max=', + 2e11.4)") utlat%sname,utlat%zpht,zmin,zmax case ('indep') ! zm at ht-indep write(msgout,"(a,' ht-indep zonal means min,max=', + 2e11.4)") utlat%sname,zmin,zmax case ('integ') ! zm at ht-integ write(msgout,"(a,' ht-integ zonal means min,max=', + 2e11.4)") utlat%sname,zmin,zmax case default end select case ('slt') select case (utlat%vtype) case ('zp ') ! slt at zp write(msgout,"(a,' zp=',f6.1,' slt=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%zpht,utlat%slt,zmin,zmax case ('ht ') ! slt at ht write(msgout,"(a,' ht=',f6.1,' slt=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%zpht,utlat%slt,zmin,zmax case ('indep') ! slt at ht-indep write(msgout,"(a,' ht-indep slt=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%slt,zmin,zmax case ('integ') ! slt at ht-integ write(msgout,"(a,' ht-integ slt=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%slt,zmin,zmax case default end select case ('lon') select case (utlat%vtype) case ('zp ') ! lon at zp write(msgout,"(a,' zp=',f6.1,' lon=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%zpht,utlat%glon,zmin,zmax case ('ht ') ! lon at ht write(msgout,"(a,' ht=',f6.1,' lon=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%zpht,utlat%glon,zmin,zmax case ('indep') ! lon at ht-indep write(msgout,"(a,' ht-indep lon=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%glon,zmin,zmax case ('integ') ! lon at ht-integ write(msgout,"(a,' ht-integ lon=',f7.1,' min,max=', + 2e11.4)") utlat%sname,utlat%glon,zmin,zmax case default end select case default end select end subroutine mkutlatlabs !------------------------------------------------------------------- subroutine wrout_utlat(h,hcntr) use plt,only: zmin,zmax use proc use input ! ! Args: type(history),intent(in) :: h,hcntr ! ! Locals: character(len=80) :: msgout ! if (iwrdat==0.and.iwrxdr==0) return ! ! Make labels if not done for plots: if (iplot == 0) call mkutlatlabs(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,utlat%data,utlat%nx,utlat%ny, + utlat%xx,utlat%yy,utlat%xlab,utlat%ylab,utlat%tlabs(2), + utlat%tlabs(3),utlat%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,utlat%data,utlat%nx,utlat%ny, + utlat%xx,utlat%yy,utlat%xlab,utlat%ylab, + utlat%tlabs(2),utlat%tlabs(3),utlat%blabs(2), + utlat%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 end subroutine wrout_utlat !------------------------------------------------------------------- subroutine pltutlat(utlat,isltax,vp) use plt use proc,only: dlat ! ! Args: type(utlat_type),intent(in) :: utlat real,intent(in) :: vp(4) integer,intent(inout) :: isltax ! ! Local: integer :: jlat0,jlat1,ny ! ! External: integer,external :: ixfind ! util.F ! ! Set y-axis to plot according to utlat%pltlat0,1: ! pltlat0,1 are from namelist read utlat_latminmax(2) ! utlat%yy, utlat%ny, and utlat%data are currently dimensioned ! and defined at the full gcmlat(nlat) range. ! jlat0 = ixfind(utlat%yy,utlat%ny,utlat%pltlat0,dlat) jlat1 = ixfind(utlat%yy,utlat%ny,utlat%pltlat1,dlat) if (jlat0<0.or.jlat0>utlat%ny.or.jlat1<0.or.jlat1>utlat%ny.or. | jlat0>=jlat1) then write(6,"('>>> pltutlat: bad jlat0,1=',2i4, | ' utlat%pltlat0,1=',2f9.3)") jlat0,jlat1,utlat%pltlat0, | utlat%pltlat1 write(6,"('>>> will resort to full latitude range..')") jlat0 = 1 jlat1 = utlat%ny endif ny = jlat1-jlat0+1 ! write(6,"('pltutlat: jlat0,1=',2i4,' pltlat0,1=',2f9.3, ! | ' ny=',i3)") jlat0,jlat1,utlat%pltlat0,utlat%pltlat1,ny ! ! Set up conpack: call cpseti('SET',0) call cpseti('MAP',0) call cpsetr('XC1',utlat%xx(1)) call cpsetr('XCM',utlat%xx(utlat%nx)) call cpsetr('YC1',utlat%yy(jlat0)) call cpsetr('YCN',utlat%yy(jlat1)) call set(vp(1),vp(2),vp(3),vp(4), + utlat%xx(1),utlat%xx(utlat%nx),utlat%yy(jlat0), + utlat%yy(jlat1),1) ! ! Contour: call contour(utlat%data(:,jlat0:jlat1),utlat%nx,utlat%nx,ny) ! ! Add axes and labels: call labutxy(utlat%mtimes,utlat%nx,utlat%yy(jlat0:jlat1),ny, + utlat%ylab,0.,isltax,utlat%glon) ! end subroutine pltutlat !----------------------------------------------------------------------- subroutine interp_lonslice(nlon,nlat,nlev,fglb,glon,gcmlon,flon) ! ! Interpolate to desired longitude glon. ! use proc,only: dlon ! ! Args: integer,intent(in) :: nlon,nlat,nlev real,intent(in) :: fglb(nlon,nlat,nlev),glon,gcmlon(nlon) real,intent(out) :: flon(nlat,nlev) ! ! Local: integer :: i,i0,i1,j,ier ! call bracket(glon,gcmlon,nlon,dlon,i0,i1,ier) do j=1,nlat flon(j,:) = fglb(i0,j,:)+(fglb(i1,j,:)-fglb(i0,j,:))* | (glon-gcmlon(i0))/(gcmlon(i1)-gcmlon(i0)) enddo end subroutine interp_lonslice !----------------------------------------------------------------------- end module mk_utlat