Attachment 189171 Details for Bug 281331 – The source file that caused the segfault

The source file that caused the segfault

metfields_ACEv2.2MetO-gfortran.f (text/x-fortran), 68.44 KB, created by William Daffer on 2007-09-06 19:41:54 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: William Daffer

Created: 2007-09-06 19:41:54 UTC

Size: 68.44 KB

patch

obsolete

>      program pvACE   
>      
>      parameter(levamax=150,mlevmax=66,levtrop=30,levels=27,maxtm=2)
>c      include '/users/manney/epv/readlib/gridmax.f'
>      include 'gridmax.f'
>      include 'getopt.h'
>      include 'stdio.h'
>      integer pmax
>      parameter(pmax = 401, neqmax = 72)
>      real kap,faktor,deltaz,ka1,gamma,pliml,plimu,pvtrop
>      parameter(gamma=-0.002,pvtrop=3.5e-6)
>      parameter(kap=0.286)
>      parameter(faktor = -9.81/287.0)
>      parameter(deltaz = 2000.0)
>      parameter(ka1=kap-1.)
>      parameter(pliml=6500.,plimu=55000.)
>    
>      real*4 theta(mlevmax),phi1,lambda1,lambda2arg
>      real*4 alev(73),alevpv(73),alevin(73)
>      real*4 zint(levamax),xpos(levamax),ypos(levamax)
>
>      real*4 pvdata(maxdata,maxtm)
>      real*4 eqldata(maxdata,maxtm)
>      real*4 tmpdata(maxdata,maxtm)
>      real*4 hgtdata(maxdata,maxtm)
>      real*4 udata(maxdata,maxtm)
>      real*4 vdata(maxdata,maxtm)
>      real*4 wsdata(maxdata,maxtm)
>
>      real*4 eint(merm3,latm3),tcalc(maxdata,2)
>      real*4 pcalc(latm3,merm3,levels),ecalc(latm3,merm3,levels)
>      real*4 gcalc1(latm3,merm3),gcalc2(latm3,merm3)
>      real*4 hcalc1(latm3,merm3),hcalc2(latm3,merm3)
>      real*4 stabil(levels),grdwrk(latm3*merm3*mlevmax,2)
>      real*4 gdata(latm3*merm3*mlevmax,2),xing(1),ying(1)
>      real*4 xlat(latm3),xlatpv(latm3),xlon(merm3),xlonpv(merm3)
>      real*4 tgrd(mlevmax,1),zin(levtrop)
>      real*4 zint3(levamax,levels),zintp(1,levels),zoutp(1,1)
>      real*4 tmpin(levtrop),pnew(1)
>      real*4 pin(levtrop),pvtin(levtrop)
>      real*4 pvlev(pmax),pvgrde(latm3,mlevmax)
>      real*4 wseql(latm3,mlevmax),pvgrd2(latm3,mlevmax),pveql(latm3,mlevmax)
>      real*4 wsgrd(neqmax,1),pvcalc(latm3,merm3,mlevmax),xeqhem(neqmax)
>      real*4 pvin(latm3,merm3,mlevmax,1),eqcalc(latm3,merm3,mlevmax)
>      real*4 pvcalc2(pmax,merm3,mlevmax),eqcalc2(latm3,merm3,mlevmax)
>      real*4 pvhem(neqmax,1),pvghem(neqmax,1),pvg2hem(neqmax,1)
>      real*4 wshem(neqmax)
>
>      real*4 pvout(mlevmax),eqlout(mlevmax),spvout(mlevmax)
>      real*4 delhpv(mlevmax),losdpv(mlevmax),tmpout(mlevmax)
>      real*4 delhtmp(mlevmax),losdtmp(mlevmax)
>      real*4 hgtout(mlevmax),uout(mlevmax),vout(mlevmax)
>      real*4 eqledge(mlevmax),eqlinner(mlevmax),eqlouter(mlevmax)
>      real*4 junk(73)
>      real*4 dyntrop,tmptrop
>
>      real*4 pres(4),dtdzs(4),potemp(4)
>      integer*4 ir(latmax,mermax,levmax),ngoodS(mlevmax,maxtm),Tfit,pstart
>      integer*4 irp(latmax,mermax,levmax)
>      integer*4 ngoodN(mlevmax,maxtm),ngoodAllN(maxtm),ngoodAllS(maxtm)
>      logical norm,ftrop,thettrop,first,occdata,zinterp
>      character csfc,fldnam*11,cslin,avgmer,cpole
>      character*6 idt
>      character*16 keyword
>      character*30 value
>c------------------------------------------------------------------------
>      character*80 glc_name1, glc_name
>      character*80 ace_name1, ace_name
>      character*80 outfile1, outfile
>
>      character lonin*30,latin*29,datetime*48,nametag*34
>      character aceocc*7,year*4,doy*3,acever*3
>      real*4 press(levamax)         ! ACE Pressure profile (mb)
>      real*4 atemp(levamax)         ! ACE Tempereature profile (K)
>      real*4 zalt(levamax)  ! Altitude (km)
>
>
>      real*4 xsumsqr, ysumsqr !to check glc lat/lon (whd)      
>
>      integer i,j,k,l
>      logical found
>
>      integer index1, index2, index3
>      integer iostat
>      character*1 blank/' '/
>      character*80 rcsid;
>
>      data found/.false./
>      data lenpv/8368/
>      character*256 datatop,dmptop
>      character*4 dmpversion
>      data dmpversion/'v1.1'/
>      integer*4 ll
>
>
>
>      data pstart/37/
>c----------------------------------------------------------------------
>
>      data dud/1.e12/,duda/-999./,p1000/1000./
>      data pres/1013.25,226.32,54.748,8.6801/
>      data potemp/288.15,332.54,498.93,891.44/
>      data dtdzs/-6.5,0.0,1.0,5.6/
>      data rg,cp /287.,1004./
>      data pi/3.141592654/
>      data g/9.81/,rearth/6371./
>      character locflag*30
>      data locflag/'!LOCATION_FLAG             '/
>
>      character*30 date
>      character*55 dmp_version, dmp_date
>      data dmp_version/'!DMP_VERSION                       '/
>      data dmp_date   /'!DMP_DATE                                          '/
>c     '
>c      character*30 CTime(
>c      integer*8 Time8
>c      external CTime, Time
>
>      real*4 stdatm/1013.25/
>
>
>cc --------------------------------------------------
>cc    getopt stuff
>cc --------------------------------------------------
>
>      integer      getopt
>      character    flag
>      character*16 optstring/'vg:i:o:'/
>      logical      argset, firstcall/.true./
>
>      logical gotinfile/.false./,
>     &     gotoutfile/.false./,
>     &     gotglcfile/.false./,
>
>      logical infile_exists/.false./,
>     &        glcfile_exists/.false./
>      logical verbose/.false./
>
>
>      integer*4 slashloc(20)
>
>      character*20 units
>      
>
>      call CTime(TIME8(),date)
>C      date = CTime(Time())
>      dmp_date = dmp_date(1:23)//date
>      dmp_version = dmp_version(1:23)//'MetO '//dmpversion
>
>      i=1
>      status=1
>      do while (i.le.maxargs.and.status.eq.1) 
>        status=getopt(optstring,flag,argset,nargs,firstcall)
>        firstcall=.false.
>        if (argset) then
>          if (flag.eq.'i') then
>            read(optarg,'(a)') infile
>            gotinfile=.true.
>          elseif (flag.eq.'o') then
>            read(optarg,'(a)') outfile
>            gotoutfile=.true.
>          elseif (flag.eq.'g') then
>            read(optarg,'(a)') glcfile
>            gotglcfile=.true.
>          else 
>            print *,'unrecognized flag', flag, ' Aborting!'
>            call exit(1)
>          endif 
>        else 
>          if (flag.eq.'h') then
>            print *,'usage :'
>            print *,'metfields_ACEv2.2MetO [-hv] i infile -g glcfile -o outfile
>c           '
>            print *,''
>            print *,'-i in: name of input file (required)'
>            print *,'-o out: name of output file (required)'
>            print *,'-g glc: name of geolocation file (required)'
>            print *,'-h : print this message, then exit'
>            print *,'-v : verbose (not implemented yet ;-)'
>            call exit(0)
>          elseif (flag.eq.'v') then
>            verbose=.true.
>          endif 
>        endif 
>      enddo 
>
>      if (.not.(gotglcfile.and.gotinfile.and.gotoutfile)) then
>        print *,'Missing a required input'
>        
>        if (.not.gotglcfile) print *,'Missing name of glc file'
>        if (.not.gotinfile) print *,'Missing name of input .asc file'
>        if (.not.gotoutfile) print *,'Missing name of output MetO file'
>
>        call exit(1)
>      endif 
>
>      glc_name=glcfile
>      ace_name=infile
>      
>
>C      Check to see if the input files exist
>      inquire(10,file=glcfile(1:len_trim(glcfile)),exist=glcfile_exists)
>      inquire(10,file=infile(1:len_trim(infile)),exist=infile_exists)
>
>      if (.not.infile_exists) then 
>        print *'### Error: ',infile(1:len_trim(infile)),
>     &       ' doesn''t exist! -- Aborting'
>        call exit(1)
>      endif 
>
>
>c      print *,'dmp_version = ',dmp_version
>c      print *,'dmp_version = ',dmp_date
>
>      degtorad = pi/180.
>      
>      
>      print 1113,'glc_name = ',glc_name
>      print 1113,'ace_name = ',ace_name
>      print 1113,'outfile = ',outfile
>c     call exit(0)
>      
>cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>c     Open input ACE file (in this case, the .csv file
>c     what we call the GLC or 'geolocation' file
>cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>      
> 1113 format (2(a,1x))
>      
>      open(10,file=glc_name,status='old',iostat=ios)
>      if (ios.ne.0) then
>        print 1114,'### Error opening ', glc_name
>c        call exit(1)
>      else 
>        occdata = .false.
>        read(10,'(a34)')nametag
>        read(10,'(a48)')datetime
>        read(10,'(a29)')latin
>        read(10,'(a30)')lonin
>        do n=1,3
>          read(10,*,end=333)
>        enddo
>        occdata = .true.
> 333    continue
>        
>        read(datetime(32:33),'(i2)')istday
>        read(datetime(29:30),'(i2)')istmon
>        read(datetime(24:27),'(i4)')istyear
>        read(datetime(35:36),'(i2)')ihours
>        read(datetime(38:39),'(i2)')imins 
>        read(datetime(41:45),'(f5.2)')secs
>        rtime = float(ihours)*3600. + float(imins)*60. + secs
>        jstday = julday(istmon,istday,istyear)
>        read(latin(24:29),*)slat   
>        read(lonin(24:30),*)slon   
>        if(slon.lt.0.) slon = slon + 360.
>        levread=0
>        xsumsqr=0.0
>        ysumsqr=0.0
>        
>ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>c     
>c     Read the location data. Check for lat/lon identically == 0.  If
>c     this is a 'header only' file or the locations are == 0, fill the
>c     locations with the 'nominal' locations read from the header and
>c     set LOCATION_FLAG == 'FILL'
>c     
>ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>        
>        do n=1,levamax
>          read(10,*,end=222)zalt(n),press(n),atemp(n),ypos(n),xpos(n)
>          levread = levread+1
>          if(atemp(n).eq.0.)atemp(n) = -999.000
>          if(press(n).ne.-999.000 .and. press(n).ne.0.) then
>            press(n) = stdatm*press(n)
>          else 
>            press(n) = dud
>          endif
>          if (xpos(n).ne.-999.0) then 
>            if (xpos(n).lt.0) xpos(n) = xpos(n) + 360.
>            xsumsqr = xsumsqr + xpos(n)*xpos(n)
>          else 
>            xpos(n) = dud
>          endif 
>          if (ypos(n).ne.-999.0) then 
>            ysumsqr = ysumsqr + ypos(n)*ypos(n)
>          else 
>            ypos(n) = dud
>          endif 
>          print *,zalt(n),press(n),atemp(n),ypos(n),xpos(n)
>        enddo
> 222    print *,levread,' levels read from GLC file'
>        if (levread.eq.0) then
>          print*,'### Info: ', glc_name(1:len(glc_name)), 
>     &         ' is body-less file'
>        else 
>          levcalc = min(mlevmax,levread)
>        endif 
>        
>      endif                             ! Come from open of GLC file.
>      
>      
>c     
>c      print *,'after read/conversion:'
>      
>      
>      if (levread.eq.0.or.
>     &     (xsumsqr .eq. 0.0 .and. ysumsqr .eq. 0.0 )) then 
>        ll=len(glc_name)
>        locflag = locflag(1:23)//'FILL'
>        glc_ios=ios
>        if (glc_ios.ne.0) then 
>          print *,'%%% ', glc_name(1:ll),
>     &         ' is missing, reading .asc file for everything'
>        else 
>          print *,' %%% ', glc_name(1:ll), 
>     &         'is body-less or has <lat>=<lon>=0.0'
>          print *,' %%% Using scalar locations'
>          do n=1,mlevmax
>            xpos(n) = slon
>            ypos(n) = slat
>          enddo
>        endif 
>
>        
>c     
>c     Attempt to read the .asc file for the pressure/alt/temp
>c     
>
>        if (levread.eq.0.or.glc_ios.ne.0) then
>
>          close(10)
>          
>          open(10,file=ace_name,status='old',iostat=ios)
>          if (ios.ne.0) then
>            ll=len_trim(ace_name)
>            print 1114,'### Error opening ', ace_name(1:ll), 
>     &           ' Aborting!'
>            call exit(1)
>          endif 
>          
>          occdata = .false.
>          if (glc_ios.ne.0) then 
>c           GLC file wasn't there. Read the header of the .asc file to 
>c           get (scalar) GLC info           
>            nametag = '!OCCULTATION_NAME     '
>            latin   = '!SCALAR_LATITUDE      '
>            lonin   = '!SCALAR_LONGITUDE     '
>            datetime= '!SCALAR_DATE          '
>            
>            do n=1,12
>              read(10,'(a16,a)',end=444) keyword, value
>              if (keyword(1:8).eq.'latitude') then
>                read(value,'(F5.2)') slat
>                write(latin(23:),'(f7.2)') slat
>              elseif (keyword(1:9).eq.'longitude') then
>                read(value,'(F7.2)') slon
>                write(lonin(23:),'(f7.2)') slon
>                if(slon.lt.0.) slon = slon + 360.
>              elseif (keyword(1:4).eq.'date') then
>                ll=len_trim(value)
>                read(value,1234)istyear,istmon,istday,
>     &               ihours,imins,secs,ijunk
>                write(datetime(24:48),'(a)') value
>                rtime = float(ihours)*3600. + float(imins)*60. + secs
>                jstday = julday(istmon,istday,istyear)
>              elseif (keyword(1:4).eq.'name') then
>                write(nametag(24:34),'(a)') value
>              endif 
> 1234         format(i4,'-',i2,'-',i2,1x,2(i2,':'),f5.2,'+',i2)
>              
>            enddo
>            do n=1,mlevmax
>              xpos(n) = slon
>              ypos(n) = slat
>            enddo
>
>          else 
>c           GLC file was body-less, so we have info from the GLC file, skip header 
>c           and just read Z/P/T
>            do n=1,12
>              read(10,*,end=334)
>            enddo    
>          endif 
>          occdata = .true.
> 334      continue
>          print *,'%%% reading Z/T/P levels from .asc file'        
>          do n=1,levamax
>            read(10,*,end=223)zalt(n),atemp(n),junk(1),press(n)
>            levread = levread+1
>            if(atemp(n).eq.0.)atemp(n) = -999.000
>            if(press(n).ne.-999.000 .and. press(n).ne.0.) then
>              press(n) = stdatm*press(n)
>            else 
>              press(n) = dud
>            endif
>            print *,zalt(n),press(n),atemp(n)
>          enddo
> 223      print *,'%%%', levread,
>     &         ' levels read from .ASC file'
> 444      continue
>          if (levread.eq.0) then
>            print*,'### Error: ', ace_name(1:len(ace_name)), 
>     &           ' is body-less, can''t get location info! Aborting'
>            call exit(1)
>          else 
>            levcalc = min(mlevmax,levread)
>          endif 
>        endif 
>      else 
>        locflag = locflag(1:23)//'GLC'
>        print *,(xpos(i), i=1,levcalc)
>        print *,(ypos(i), i=1,levcalc)
>      endif 
>      print *,'levcalc = ',levcalc
>      
>      
> 1114 format(3(a,1x))
>      cpole = 'G'
>
>      print *,'occdata: ',occdata
>      
>      if(occdata)then
>        
>        if(rtime.lt.43200.)then
>          jday1 = jstday - 1
>          jday2 = jstday 
>        elseif(rtime.gt.43200.)then
>          jday1 = jstday
>          jday2 = jstday + 1
>        else
>          jday1 = jstday
>          jday2 = jstday
>        endif
>        index1 = 1
>        index2 = 2
>        if(rtime.lt.43200.)then
>          timesec = rtime + 86400.
>        else
>          timesec = rtime
>        endif
>        alpha = (timesec - 43200.)/86400.                
>        
>        print *,' jday1,jday2,index1,index2,alpha --------------------'
>        print *,jday1,jday2,index1,index2,alpha
>        print *,'istday,istmon,istyear,ihours,imins,secs,rtime,slon,slat ---------- '
>c     '
>        print *,istday,istmon,istyear,ihours,imins,secs,rtime,slon,slat
>        
>        jstart = jday1
>        jend = jday2
>        
>        num_ncep = jend - jstart + 1
>        
>        do n=1,num_ncep 
>          jday = jstart + n - 1
>          call caldat(jday,imon,iday,iyear)
>          write(idt(1:2),'(i2)')iday
>          write(idt(3:4),'(i2)')imon
>          if(iyear.lt.2000)then
>            write(idt(5:6),'(i2)')iyear - 1900
>          else
>            write(idt(5:6),'(i2)')iyear - 2000
>          endif
>          do ii=1,6
>            if(idt(ii:ii).eq.' ') idt(ii:ii) = '0'
>          enddo                                              
>          print *,'idt: ', idt
>          
>          csfc = 't'
>          itype = -7
>          call rukmo_read2(1,idt,itype,pvdata(1,n),alevin,merin,latsin,fldnam,ierr_r)
>          call indices(irp,latsin,merin,levmax)
>          do k=1,nint(alevin(73))
>            vmin = 1.e15
>            vmax = -1.e15
>            do i=1,merin 
>              do j=1,latsin
>                if(pvdata(irp(j,i,k),n).ne.dud) then
>                  pvdata(irp(j,i,k),n) = pvdata(irp(j,i,k),n)*1.e4
>c     vmin = min(vmin,pvdata(irp(j,i,k),n))
>c     vmax = max(vmax,pvdata(irp(j,i,k),n))
>                else
>                  pvdata(irp(j,i,k),n) = dud
>                endif
>              enddo
>            enddo
>c     print *,'PV input:',alevin(k),vmin,vmax
>          enddo
>          csfc = 't'
>          itype = 11
>          call rukmo_read2(1,idt,itype,eqldata(1,n),alevin,merin,latsin,fldnam,ierr_r)
>          if(n.eq.1)then
>            merpv = merin
>            latpv = latsin
>            levpv = nint(alevin(73))
>            do k=1,levpv
>              alevpv(k) = alevin(k)
>            enddo
>            
>c     If 2nd day and grid has changed, interpolate to 1st day's grid
>            
>          elseif(merin.ne.merpv .or. latsin.ne.latpv)then
>            print *,'*****NOTE: Lat/lon spacing has changed*****'
>            levin = nint(alevin(73))
>            zinterp = .true.
>c     print *,'calling gridinterp:',merin,latsin,merpv,latpv,' ',csfc,' ',zinterp
>c     print *,(alevin(k),k=1,levin)
>c     print *,(alevpv(k),k=1,levpv)
>            call gridinterp(merin,latsin,90.,0.,merpv,latpv,90.,
>     &           0.,levin,alevin,levpv,alevpv,zinterp,csfc,
>     &           pvdata(1,n),dud)
>            call gridinterp(merin,latsin,90.,0.,merpv,latpv,90.,
>     &           0.,levin,alevin,levpv,alevpv,zinterp,csfc,
>     &           eqldata(1,n),dud)
>            call indices(irp,latpv,merpv,levmax)
>          endif
>
>c     do k=1,levpv
>c     vmin = 1.e15
>c     vmax = -1.e15
>c     do i=1,merpv 
>c     do j=1,latpv
>c     if(pvdata(irp(j,i,k),n).ne.dud) then
>c     vmin = min(vmin,pvdata(irp(j,i,k),n))
>c     vmax = max(vmax,pvdata(irp(j,i,k),n))
>c     endif
>c     enddo
>c     enddo
>c     print *,'PV after gridinterp:',alevpv(k),vmin,vmax
>c     enddo
>          
>          csfc = 'p'
>          itype = 1
>          call rukmo_read2(1,idt,itype,udata(1,n),alevin,merin,latsin,
>     &         fldnam,ierr_r)
>          if (ierr_r.ne.0) then 
>            print *,'Bad return from rukmo_read for itype=',itype
>            call exit(1)
>          endif 
>          csfc = 'p'
>          itype = 2
>          call rukmo_read2(1,idt,itype,vdata(1,n),alevin,merin,latsin,
>     &         fldnam,ierr_r)
>          if (ierr_r.ne.0) then 
>            print *,'Bad return from rukmo_read for itype=',itype
>            call exit(1)
>          endif 
>          csfc = 'p'
>          itype = 5
>          call rukmo_read2(1,idt,itype,tmpdata(1,n),alevin,merin,latsin,
>     &         fldnam,ierr_r)
>          if (ierr_r.ne.0) then 
>            print *,'Bad return from rukmo_read for itype=',itype
>            call exit(1)
>          endif 
>          csfc = 'p'
>          itype = 6
>          call rukmo_read2(1,idt,itype,hgtdata(1,n),alevin,merin,latsin,
>     &         fldnam,ierr_r)
>          if (ierr_r.ne.0) then 
>            print *,'Bad return from rukmo_read for itype=',itype
>            call exit(1)
>          endif 
>          if(n.eq.1)then
>            merids = merin
>            lats = latsin
>            call indices(ir,lats,merids,levmax)
>            levp = nint(alevin(73))
>            do k=1,levp
>              alev(k) = alevin(k)
>            enddo
>c     If 2nd day and grid has changed, interpolate to 1st day's grid
>          elseif(merin.ne.merids .or. latsin.ne.lats)then
>            print *,'*****NOTE: Lat/lon spacing has changed*****'
>            levin = nint(alevin(73))
>            zinterp = .true.
>            call gridinterp(merin,latsin,90.,0.,merids,lats,90.,
>     &           0.,levin,alevin,levp,alev,zinterp,csfc,
>     &           udata(1,n),dud)
>            call gridinterp(merin,latsin,90.,0.,merids,lats,90.,
>     &           0.,levin,alevin,levp,alev,zinterp,csfc,
>     &           vdata(1,n),dud)
>            call gridinterp(merin,latsin,90.,0.,merids,lats,90.,
>     &           0.,levin,alevin,levp,alev,zinterp,csfc,
>     &           tmpdata(1,n),dud)
>            call gridinterp(merin,latsin,90.,0.,merids,lats,90.,
>     &           0.,levin,alevin,levp,alev,zinterp,csfc,
>     &           hgtdata(1,n),dud)
>          endif
>          do k=1,levp
>            vmin = 1.e15
>            vmax = -1.e15
>            do i=1,merids
>              do j=1,lats
>                if(udata(ir(j,i,k),n).ne.dud .and.
>     &               vdata(ir(j,i,k),n).ne.dud)then
>                  wsdata(ir(j,i,k),n) =
>     &                 sign( sqrt( udata(ir(j,i,k),n)*udata(ir(j,i,k),n) +
>     &                 vdata(ir(j,i,k),n)*vdata(ir(j,i,k),n) ),
>     &                 udata(ir(j,i,k),n) )
>                  vmin = min(vmin,wsdata(ir(j,i,k),n))
>                  vmax = max(vmax,wsdata(ir(j,i,k),n))
>                else
>                  wsdata(ir(j,i,k),n) = dud
>                endif
>              enddo
>            enddo
>c     print *,'windspeed:',alev(k),vmin,vmax
>          enddo
>c     Interpolate windspeed to PV grid
>          if(latpv.ne.lats .or. merpv.ne.merids)then
>            zinterp = .false.
>            call gridinterp(merids,lats,90.,0.,merpv,latpv,90.,0.,
>     &           levp,alev,levp,alev,zinterp,csfc,
>     &           wsdata(1,n),dud)
>          endif
>          
>          if(levp.lt.levels)then
>            do k=levp+1,levels
>              do i=1,merids
>                do j=1,lats
>                  udata(ir(j,i,k),n) = dud
>                  vdata(ir(j,i,k),n) = dud
>                  tmpdata(ir(j,i,k),n) = dud
>                  hgtdata(ir(j,i,k),n) = dud
>                enddo
>              enddo
>              do i=1,merpv
>                do j=1,latpv
>                  wsdata(irp(j,i,k),n) = dud
>                enddo
>              enddo
>            enddo
>          endif
>          
>          if(levpv.lt.levels)then
>            do k=levpv+1,levels
>              do i=1,merpv
>                do j=1,latpv
>                  pvdata(irp(j,i,k),n) = dud
>                  eqldata(irp(j,i,k),n) = dud
>                enddo
>              enddo
>            enddo
>          endif
>          
>          ngoodAllN(n) = 0
>          ngoodAllS(n) = 0
>          do k=1,levels
>            vmin = 1.e15
>            vmax = -1.e15
>            ngoodN(k,n) = 0
>            ngoodS(k,n) = 0
>            do i=1,merpv 
>              do j=1,latpv
>                if(j.gt.37 .and. pvdata(irp(j,i,k),n).ne.dud)then
>                  vmin = min(vmin,pvdata(irp(j,i,k),n))
>                  vmax = max(vmax,pvdata(irp(j,i,k),n))
>                  ngoodN(k,n) = ngoodN(k,n) + 1
>                elseif(pvdata(irp(j,i,k),n).ne.dud)then
>                  ngoodS(k,n) = ngoodS(k,n) + 1
>                  vmin = min(vmin,pvdata(irp(j,i,k),n))
>                  vmax = max(vmax,pvdata(irp(j,i,k),n))
>                endif
>              enddo
>            enddo
>            print *,'PV on theta:',k,alevpv(k),vmin,vmax
>            ngoodAllN(n) = ngoodAllN(n) + ngoodN(k,n)
>            ngoodAllS(n) = ngoodAllS(n) + ngoodS(k,n)
>          enddo
>
>
>        enddo                           !end of loop reading met data
>
>
>
>cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>        
>        dlat = 180./float(lats - 1)
>        dlatpv = 180./float(latpv - 1)
>        leq = lats/2 + 1
>        if(lats.gt.73)then
>          neqhem = 72
>          xeqlst = 36.
>          j35 = 169                     !indices of 36N/36S latitudes
>          jm35 = 73                
>          jsp = 2
>          j10 = 14                      !number of lat gridpoints in 10.5 degs
>        else
>          neqhem = 22
>          xeqlst = 35.
>          j35 = 51                      !indices of 35N/35S latitudes
>          jm35 = 23                
>          jsp = 2
>          j10 = 4                       !number of lat gridpoints in 10 degs
>        endif
>        
>        do j=1,pmax
>          pvlev(j) = -6.0 + 0.03*(j-1)
>        enddo
>        do j=1,neqhem
>          xeqhem(j) = xeqlst + (j-1)*dlatpv
>        enddo
>        do j=1,lats
>          xlat(j) = -90. + (j-1)*dlat
>        enddo
>        do j=1,lats
>          xlatpv(j) = -90. + (j-1)*dlatpv
>        enddo
>        dlon = 360./float(merids)
>        dlonpv = 360./float(merpv)
>        do i=1,merids
>          xlon(i) = (i-1)*dlon
>        enddo
>        do i=1,merids
>          xlonpv(i) = (i-1)*dlonpv
>        enddo
>c     print *,'xlat:',(xlat(n),n=1,lats)
>c     calculate standard atmosphere
>c     static stability (ref. Baldwin)
>        do l=1,levpv 
>          if(alevpv(l).lt.370.)then
>            index3 = 2
>            dtdzin = dtdzs(2)
>          elseif(alevpv(l).lt.440.)then
>            index3 = 2
>            frac = ( log(alevpv(l))-log(potemp(2)) )/
>     &           ( log(potemp(3))-log(potemp(2)) )
>            dtdzin = dtdzs(2)*(1.-frac) + dtdzs(3)*frac
>          else
>            index3 = 3
>            dtdzin = dtdzs(3)
>          endif
>          
>          pstth = pres(index3)*(alevpv(l)/potemp(index3))**(-g/(rg*
>     &         (dtdzin*0.001+g/cp)))
>          pvsign = -1./g
>          stabil(l) = (-(rg*alevpv(l))/(pstth*g))*(g/cp+0.001*dtdzin)
>        enddo
>c     print *,'stabil:',stabil
>c     Temperature interp
>        
>c     call minmax(tmpdata(1,index1),merids*lats*levp,vmin1,vmax1,dud,0)
>c     call minmax(tmpdata(1,index2),merids*lats*levp,vmin2,vmax2,dud,0)
>c     print *,'min, max T:',vmin1,vmin2,vmax1,vmax2
>        
>c     print *,'before T interp:'
>c     print *,(xpos(i),i=1,levcalc)
>c     print *,(ypos(i),i=1,levcalc)
>        do k=1,levp
>          rmin = 1.e15
>          rmax = -1.e15
>          do i=1,merids
>            do j=1,lats
>              if(tmpdata(ir(lats+1-j,i,k),index1).ne.dud .and.
>     &             tmpdata(ir(lats+1-j,i,k),index2).ne.dud)then
>                eint(i,j) = tmpdata(ir(lats+1-j,i,k),index1)*(1.-alpha) +
>     &               tmpdata(ir(lats+1-j,i,k),index2)*alpha
>              else
>                eint(i,j) = dud
>              endif
>              if(eint(i,j).ne.dud)then
>                rmin = min(eint(i,j),rmin)
>                rmax = max(eint(i,j),rmax)
>              endif
>            enddo
>          enddo
>          call grid_pnts(merm3,latm3,merids,lats,eint,90.,0.,
>     &         levcalc,zint,xpos,ypos,levamax,dud)
>          do np=1,levcalc
>            zint3(np,k) = zint(np)
>          enddo
>        enddo
>        
>        do np=1,levcalc
>          if(press(np).ne.dud)then
>            pnew(1) = press(np)
>            do k=1,levp
>              zintp(1,k) = zint3(np,k)
>            enddo
>            cslin = 'l'
>            call psplint_1d(1,levels,1,levp,alev,1,pnew,zintp,
>     &           zoutp,dud,cslin)
>            tmpout(np) = zoutp(1,1)
>c     theta from met analysis T
>            if(tmpout(np).ne.dud)then
>              theta(np) = tmpout(np)*(1000./press(np))**(2./7.)
>            else
>              theta(np) = dud
>            endif
>          else
>            tmpout(np) = dud
>            theta(np) = dud
>          endif
>        enddo
>        call minmax(tmpout,levcalc,vmin,vmax,dud,0)
>c     print *,'tmpout:',vmin,vmax
>        print *,'tmpout:',(tmpout(np),np=1,51)
>        call minmax(theta,levcalc,vmin,vmax,dud,0)
>        print *,'theta:',vmin,vmax
>c     print *,'theta:',(theta(k),k=1,levcalc)
>        levcalcpv = levcalc
>        do k=levcalc,1,-1
>          if(theta(k).eq.dud) levcalcpv = levcalcpv - 1
>        enddo
>        print *,'levcalc, levcalcpv:',levcalc, levcalcpv
>        
>c     tgrad trop height
>        do k=1,levtrop
>          tmpin(levtrop+1-k) = tmpout(k)
>          if(zalt(k).ne.-999.00)then
>            zin(levtrop+1-k) = zalt(k)
>          else
>            zin(levtrop+1-k) = dud
>          endif
>          if(press(k).ne.dud)then
>            pin(levtrop+1-k) = press(k)*100.
>          else
>            pin(levtrop+1-k) = dud
>          endif
>        enddo
>**************************************************************
>*     twmo routine inline
>**************************************************************
>        trp=dud                         ! negative means not valid
>        do j=levtrop,2,-1
>                                        ! dt/dz
>          pmk= .5 * (pin(j-1)**kap+pin(j)**kap)
>          pm = pmk**(1/kap)
>          a = (tmpin(j-1)-tmpin(j))/(pin(j-1)**kap-pin(j)**kap)
>          b = tmpin(j)-(a*pin(j)**kap)
>          tm = a * pmk + b
>          dtdp = a * kap * (pm**ka1)
>          dtdz = faktor*dtdp*pm/tm
>          
>                                        ! dt/dz valid?
>          if (j.eq.level)    go to 999  ! no, start level
>          if (dtdz.le.gamma) go to 999  ! no, dt/dz < -2 K/km
>          if (pm.gt.plimu)   go to 999  ! no, too low
>          
>                                        ! dtdz is valid, calculate tropopause pressure
>          if (dtdz0.lt.gamma) then
>            ag = (dtdz-dtdz0) / (pmk-pmk0)
>            bg = dtdz0 - (ag * pmk0)
>            ptph = exp(log((gamma-bg)/ag)/kap)
>c     print *,'first trp calc:',ptph,pin(j+1),pin(j),pin(j-1)
>          else
>            ptph = pm
>          endif
>          
>          if (ptph.lt.pliml) go to 999
>          if (ptph.gt.plimu) go to 999
>          
>                                        ! 2nd test: dtdz above 2 km must not exceed gamma
>          p2km = ptph + deltaz*(pm/tm)*faktor ! p at ptph + 2km
>          asum = 0.0                    ! dtdz above
>          icount = 0                    ! number of levels above
>          
>                                        ! test until apm < p2km
>          do jj=j,2,-1
>            
>            pmk2 = .5 * (pin(jj-1)**kap+pin(jj)**kap) ! p mean ^kappa
>            pm2 = pmk2**(1/kap)         ! p mean
>            if(pm2.gt.ptph) go to 110   ! doesn't happen
>            if(pm2.lt.p2km) go to 888   ! ptropo is valid
>            
>            a2 = (tmpin(jj-1)-tmpin(jj)) ! a
>            a2 = a2/(pin(jj-1)**kap-pin(jj)**kap)
>            b2 = tmpin(jj)-(a2*pin(jj)**kap) ! b
>            tm2 = a2 * pmk2 + b2        ! T mean
>            dtdp2 = a2 * kap * (pm2**(kap-1)) ! dt/dp
>            dtdz2 = faktor*dtdp2*pm2/tm2
>            asum = asum+dtdz2
>            icount = icount+1
>            aquer = asum/float(icount)  ! dt/dz mean
>            
>                                        ! discard ptropo ?
>            if (aquer.le.gamma) go to 999 ! dt/dz above < gamma
>            
> 110        continue
>          enddo                         ! test next level
>          
> 888      continue                      ! ptph is valid
>          trp = ptph
>          if(trp.lt.pin(j))then
>            ag = -log(pin(j))
>            bg = -log(pin(j-1))
>            zm = zin(j)
>            zm2 = zin(j-1)
>          else
>            ag = -log(pin(j+1))
>            bg = -log(pin(j))
>            zm = zin(j+1)
>            zm2 = zin(j)
>          endif
>          pmk = -log(trp)
>          frac = (pmk - bg)/(ag - bg)
>          trpz = zm*frac + zm2*(1.-frac)
>          print *,'found tropopause:',trpz
>          goto 111
>          
> 999      continue                      ! continue search at next higher level
>          pm0 = pm
>          pmk0 = pmk
>          dtdz0  = dtdz
>          
>        enddo
> 111    continue
>        if(trp.ne.dud)then
>          tmptrop = trpz
>        else
>          tmptrop = dud
>        endif
>**************************************************************
>*     end twmo
>**************************************************************
>        print *,'tmptrop:',tmptrop
>        
>c     Height interp
>        
>c     call minmax(hgtdata(1,index1),merids*lats*levp,vmin1,vmax1,dud,0)
>c     call minmax(hgtdata(1,index2),merids*lats*levp,vmin2,vmax2,dud,0)
>c     print *,'min, max Z:',vmin1,vmin2,vmax1,vmax2
>        
>c     print *,'before GPH interp:'
>c     print *,(xpos(i), i=1,levcalc)
>c     print *,(ypos(i), i=1,levcalc)
>        do k=1,levp
>          rmin = 1.e15
>          rmax = -1.e15
>          do i=1,merids
>            do j=1,lats
>              if(hgtdata(ir(lats+1-j,i,k),index1).ne.dud .and.
>     &             hgtdata(ir(lats+1-j,i,k),index2).ne.dud)then
>                eint(i,j) = hgtdata(ir(lats+1-j,i,k),index1)*(1.-alpha) +
>     &               hgtdata(ir(lats+1-j,i,k),index2)*alpha
>              else
>                eint(i,j) = dud
>              endif
>              if(eint(i,j).ne.dud)then
>                rmin = min(eint(i,j),rmin)
>                rmax = max(eint(i,j),rmax)
>              endif
>            enddo
>          enddo
>c     if(k.eq.13)print *,'eint:',k,rmin,rmax
>          call grid_pnts(merm3,latm3,merids,lats,eint,90.,0.,
>     &         levcalc,zint,xpos,ypos,levamax,dud)
>          do np=1,levcalc
>            zint3(np,k) = zint(np)
>          enddo
>        enddo
>        
>        do np=1,levcalc
>          if(press(np).ne.dud)then
>            pnew(1) = press(np)
>            do k=1,levp
>              zintp(1,k) = zint3(np,k)
>            enddo
>            cslin = 'l'
>            call psplint_1d(1,levels,1,levp,alev,1,pnew,zintp,
>     &           zoutp,dud,cslin)
>            hgtout(np) = zoutp(1,1)
>          else
>            hgtout(np) = dud
>          endif
>        enddo
>        call minmax(hgtout,levcalc,vmin,vmax,dud,0)
>        print *,'hgtout:',vmin,vmax
>        
>c     Zonal wind interp
>        
>c     call minmax(udata(1,index1),merids*lats*levp,vmin1,vmax1,dud,0)
>c     call minmax(udata(1,index2),merids*lats*levp,vmin2,vmax2,dud,0)
>c     print *,'min, max U:',vmin1,vmin2,vmax1,vmax2
>        
>c     print *,'before U interp:'
>c     print *,(xpos(i), i=1,levcalc)
>c     print *,(ypos(i), i=1,levcalc)
>        do k=1,levp
>          rmin = 1.e15
>          rmax = -1.e15
>          do i=1,merids
>            do j=1,lats
>              if(udata(ir(lats+1-j,i,k),index1).ne.dud .and.
>     &             udata(ir(lats+1-j,i,k),index2).ne.dud)then
>                eint(i,j) = udata(ir(lats+1-j,i,k),index1)*(1.-alpha) +
>     &               udata(ir(lats+1-j,i,k),index2)*alpha
>              else
>                eint(i,j) = dud
>              endif
>              if(eint(i,j).ne.dud)then
>                rmin = min(eint(i,j),rmin)
>                rmax = max(eint(i,j),rmax)
>              endif
>            enddo
>          enddo
>c     if(k.eq.13)print *,'eint:',k,rmin,rmax
>          call grid_pnts(merm3,latm3,merids,lats,eint,90.,0.,
>     &         levcalc,zint,xpos,ypos,levamax,dud)
>          do np=1,levcalc
>            zint3(np,k) = zint(np)
>          enddo
>        enddo
>        
>        do np=1,levcalc
>          if(press(np).ne.dud)then
>            pnew(1) = press(np)
>            do k=1,levp
>              zintp(1,k) = zint3(np,k)
>            enddo
>            cslin = 'l'
>            call psplint_1d(1,levels,1,levp,alev,1,pnew,zintp,
>     &           zoutp,dud,cslin)
>            uout(np) = zoutp(1,1)
>          else
>            uout(np) = dud
>          endif
>        enddo
>        call minmax(uout,levcalc,vmin,vmax,dud,0)
>        print *,'uout:',vmin,vmax
>        
>c     Meridional wind interp
>        
>c     call minmax(vdata(1,index1),merids*lats*levp,vmin1,vmax1,dud,0)
>c     call minmax(vdata(1,index2),merids*lats*levp,vmin2,vmax2,dud,0)
>c     print *,'min, max V:',vmin1,vmin2,vmax1,vmax2
>        
>c     print *,'before V interp:'
>c     print *,(xpos(i), i=1,levcalc)
>c     print *,(ypos(i), i=1,levcalc)
>        do k=1,levp
>          rmin = 1.e15
>          rmax = -1.e15
>          do i=1,merids
>            do j=1,lats
>              if(vdata(ir(lats+1-j,i,k),index1).ne.dud .and.
>     &             vdata(ir(lats+1-j,i,k),index2).ne.dud)then
>                eint(i,j) = vdata(ir(lats+1-j,i,k),index1)*(1.-alpha) +
>     &               vdata(ir(lats+1-j,i,k),index2)*alpha
>              else
>                eint(i,j) = dud
>              endif
>              if(eint(i,j).ne.dud)then
>                rmin = min(eint(i,j),rmin)
>                rmax = max(eint(i,j),rmax)
>              endif
>            enddo
>          enddo
>c     if(k.eq.13)print *,'eint:',k,rmin,rmax
>          call grid_pnts(merm3,latm3,merids,lats,eint,90.,0.,
>     &         levcalc,zint,xpos,ypos,levamax,dud)
>          do np=1,levcalc
>            zint3(np,k) = zint(np)
>          enddo
>        enddo
>        
>        do np=1,levcalc
>          if(press(np).ne.dud)then
>            pnew(1) = press(np)
>            do k=1,levp
>              zintp(1,k) = zint3(np,k)
>            enddo
>            cslin = 'l'
>            call psplint_1d(1,levels,1,levp,alev,1,pnew,zintp,
>     &           zoutp,dud,cslin)
>            vout(np) = zoutp(1,1)
>          else
>            vout(np) = dud
>          endif
>        enddo
>        call minmax(vout,levcalc,vmin,vmax,dud,0)
>        print *,'vout:',vmin,vmax
>        
>c     calcs involving PV only if some good data
>        
>        if( (ypos(k).lt.0 .and. ngoodallS(index1).gt.0 .and. 
>     &       ngoodallS(index2).gt.0) .or.
>     &       (ypos(k).gt.0 .and. ngoodallN(index1).gt.0 .and. 
>     &       ngoodallN(index2).gt.0) )then
>c     PV and sPV
>          
>c     call minmax(pvdata(1,index1),merpv*latpv*levpv,vmin1,vmax1,dud,0)
>c     call minmax(pvdata(1,index2),merpv*latpv*levpv,vmin2,vmax2,dud,0)
>c     print *,'min, max PV:',vmin1,vmin2,vmax1,vmax2
>          
>          do k=1,levpv
>            rmin = 1.e15
>            rmax = -1.e15
>c     <<<<<< check 1,2 instead of index1,index2 >>>>
>            do i=1,merpv 
>              do j=1,latpv
>                if(pvdata(irp(latpv+1-j,i,k),1).ne.dud .and.
>     &               pvdata(irp(latpv+1-j,i,k),2).ne.dud)then
>                  eint(i,j) = pvdata(irp(latpv+1-j,i,k),1)*(1.-alpha) +
>     &                 pvdata(irp(latpv+1-j,i,k),2)*alpha
>                else
>                  eint(i,j) = dud
>                endif
>                if(eint(i,j).ne.dud)then
>                  rmin = min(eint(i,j),rmin)
>                  rmax = max(eint(i,j),rmax)
>                endif
>              enddo
>            enddo
>            if(k.eq.13)print *,'eint:',k,rmin,rmax
>            call grid_pnts(merm3,latm3,merpv,latpv,eint,90.,0.,
>     &           levcalc,zint,xpos,ypos,levamax,dud)
>            do np=1,levcalc
>              zint3(np,k) = zint(np)
>            enddo
>          enddo
>          
>          call thint_pnts(levamax,levcalc,levels,levpv,
>     &         alevpv,zint3,theta,pvout,dud)
>          
>c     PV 
>c     PV
>          do np=1,levcalcpv
>            if(theta(np).ne.-999.00 .and. pvout(np).ne.dud)then
>              
>              if(theta(np).lt.370.)then
>                index3 = 2
>                dtdzin = dtdzs(2)
>              elseif(theta(np).lt.440.)then
>                index3 = 2
>                frac = ( log(theta(np))-log(potemp(2)) )/
>     &               ( log(potemp(3))-log(potemp(2)) )
>                dtdzin = dtdzs(2)*(1.-frac) + dtdzs(3)*frac
>              else
>                index3 = 3
>                dtdzin = dtdzs(3)
>              endif
>              pstth = pres(index3)*(theta(np)/potemp(index3))**(-g/(rg*
>     &             (dtdzin*0.001+g/cp)))
>              pvsign = -1./g
>              stabil1 = (-(rg*theta(np))/(pstth*g))*(g/cp+0.001*dtdzin)
>              
>              spvout(np) = pvout(np)*pvsign/stabil1*1.e2
>              
>            else
>              spvout(np) = dud
>            endif 
>          enddo
>          call minmax(pvout,levcalcpv,vmin,vmax,dud,0)
>          print *,'pvout:',vmin,vmax
>c     print *,'pvout:',(pvout(k),k=1,levcalcpv)
>          call minmax(spvout,levcalcpv,vmin,vmax,dud,0)
>          print *,'spvout:',vmin,vmax
>          
>c     dyn trop height
>          do k=1,levtrop
>            if(zalt(k).ne.-999.00)then
>              zin(levtrop+1-k) = zalt(k)
>            else
>              zin(levtrop+1-k) = dud
>            endif
>            if(press(k).ne.dud)then
>              pin(levtrop+1-k) = press(k)*100.
>            else
>              pin(levtrop+1-k) = dud
>            endif
>            if(pvout(k).ne.dud)then
>              pvtin(levtrop+1-k) = pvout(k)*1.e-4
>            else
>              pvtin(levtrop+1-k) = dud
>            endif
>            tmpin(levtrop+1-k) = tmpout(k)
>          enddo
>c     print *,zin
>c     print *,pin
>c     print *,pvtin
>c     print *,tmpin
>*****************************************************************
>*     troppv subroutine inline
>******************************************************************
>          trp=dud                       ! not valid
>          trpz = dud
>          ttheta = dud
>          thettrop = .false.
>          itstart = 0
>          do k=1,levtrop
>            if(pin(k).ne.dud .and. tmpin(k).ne.dud)then
>              ttheta = tmpin(k)*(100000./pin(k))**(2./7.)
>              if(ttheta.lt.380. .and. itstart.eq.0) itstart = k-1
>            endif
>          enddo
>c     print *,'starting PV trop search',itstart,
>c     &           tmpin(itstart)*(100000./pin(itstart))**(2./7.)
>          if(abs(pvtin(itstart)).gt.pvtrop)then
>            k = itstart
>            do while(k.lt.levtrop-1 .and. pin(k).le.plimu .and. 
>     &           abs(pvtin(k)).gt.pvtrop)
>              k = k + 1
>            enddo
>            if((k.lt.levtrop-1 .and. pin(k).le.plimu) .or. 
>     &           abs(pvtin(k)).lt.pvtrop)then
>              if(abs(pvtin(k)).lt.pvtrop)then
>                if(pvtin(k).ne.dud .and. pvtin(k-1).ne.dud .and.
>     &               pin(k).ne.dud .and. pin(k-1).ne.dud)then
>                  ag = -log(pin(k-1))
>                  bg = -log(pin(k))
>                  frac = (pvtrop - abs(pvtin(k)))/(abs(pvtin(k-1)) - abs(pvtin(k)))
>                  trp = exp(-(ag*frac + bg*(1.-frac)))
>                  trpt = tmpin(k-1)*frac + tmpin(k)*(1.-frac)
>                  trpz = zin(k-1)*frac + zin(k)*(1.-frac)
>                  ttheta = trpt*(100000./trp)**(2./7.)
>c     print *,'PV trp p, T, theta:',trp,trpt,trpz,ttheta
>                else
>                  trp = dud
>                endif
>              else
>                print *,'*******warning: pvtrop match failed*******'
>                trp = dud
>                thettrop = .false.
>              endif
>            else
>              trp = dud
>            endif
>          else
>            if(pin(itstart).ne.dud .and. tmpin(itstart).ne.dud)then
>              ttheta = tmpin(itstart)*(100000./pin(itstart))**(2./7.)
>              if(ttheta.gt.380)then
>                thettrop = .true.
>              else
>                thettrop = .false.
>                trp = dud
>              endif
>            else
>              trp = dud
>            endif
>          endif
>          if(ttheta.ge.380.) thettrop = .true.
>          
>          if(thettrop)then
>            k = 2
>            thetab = dud
>            thetaa = dud
>            if(tmpin(k).ne.dud) thetaa = tmpin(k)*(100000./pin(k))**(2./7.)
>            if(tmpin(k-1).ne.dud) thetab = tmpin(k-1)*(100000./pin(k-1))**(2./7.)
>            do while(k.lt.levtrop-1 .and. pin(k).le.plimu .and. thetaa.gt.380.)
>              k = k + 1
>              thetab = dud
>              thetaa = dud
>              if(tmpin(k).ne.dud) thetaa = tmpin(k)*(100000./pin(k))**(2./7.)
>              if(tmpin(k-1).ne.dud) thetab = tmpin(k-1)*(100000./pin(k-1))**(2./7.)
>            enddo
>            if((k.lt.levtrop-1 .and. pin(k).le.plimu) .or. thetaa.lt.380.)then
>              if(thetaa.lt.380.)then
>                if(thetaa.ne.dud .and. thetab.ne.dud .and.
>     &               pin(k).ne.dud .and. pin(k-1).ne.dud)then
>                  ag = -log(pin(k-1))
>                  bg = -log(pin(k))
>                  frac = (380. - thetaa)/(thetab - thetaa)
>                  trp = exp(-(ag*frac + bg*(1.-frac)))
>                  trpz = zin(k-1)*frac + zin(k)*(1.-frac)
>c     print *,'theta trop:',trp,trpz,frac,pin(k-1),pin(k)
>                  if(abs(ypos(1)).gt.25. .or. abs(ypos(1)).gt.25.) trpz = dud
>                else
>                  trpz = dud
>                endif
>              else
>                print *,'*******warning: theta trop match failed*******'
>                trpz = dud
>              endif
>            else
>              trpz = dud
>            endif
>          endif
>          
>          if (trp.gt.plimu .or. trp.lt.pliml) then
>            trpz = dud
>          endif
>          dyntrop = trpz
>          print *,'dyntrop:',dyntrop
>******************************************************************
>*     end troppv
>******************************************************************
>          
>c     EqL interpolation
>          
>          call minmax(eqldata(1,1),merpv*latpv*levpv,vmin1,vmax1,dud,0)
>          call minmax(eqldata(1,2),merpv*latpv*levpv,vmin2,vmax2,dud,0)
>c     print *,'min, max EqL:',vmin1,vmin2,vmax1,vmax2
>          
>c     print *,'before EqL interp:'
>          do k=1,levpv
>            rmin = 1.e15
>            rmax = -1.e15
>            do i=1,merpv
>              do j=1,latpv
>                if(eqldata(irp(latpv+1-j,i,k),1).ne.dud .and.
>     &               eqldata(irp(latpv+1-j,i,k),2).ne.dud)then
>                  eint(i,j) = eqldata(irp(latpv+1-j,i,k),1)*(1.-alpha) +
>     &                 eqldata(irp(latpv+1-j,i,k),2)*alpha
>                else
>                  eint(i,j) = dud
>                endif
>                if(eint(i,j).ne.dud)then
>                  rmin = min(eint(i,j),rmin)
>                  rmax = max(eint(i,j),rmax)
>                endif
>              enddo
>            enddo
>c     if(k.eq.13)print *,'eint:',k,rmin,rmax
>            call grid_pnts(merm3,latm3,merpv,latpv,eint,90.,0.,
>     &           levcalc,zint,xpos,ypos,levamax,dud)
>            do np=1,levcalc
>              zint3(np,k) = zint(np)
>            enddo
>          enddo
>          
>          call thint_pnts(levamax,levcalc,levels,levpv,
>     &         alevpv,zint3,theta,eqlout,dud)
>          
>          call minmax(eqlout,levcalc,vmin,vmax,dud,0)
>          print *,'eqlout:',vmin,vmax
>          
>c     Normalized horizontal PV gradient
>c     <<<<< 1,2 instead of index[12] >>>>
>          do k=1,levpv
>            do i=1,merpv
>              do j=1,latpv
>                if(pvdata(irp(j,i,k),1).ne.dud)then
>                  pcalc(j,i,k) = pvdata(irp(j,i,k),1)*pvsign/stabil(k)*1.e2
>                else
>                  pcalc(j,i,k) = dud
>                endif
>              enddo
>            enddo
>            call minmax(pcalc(1,1,k),merpv*latpv,vmin,vmax,dud,0)
>c     print *,'scaled PV1 for PVG calc:',alevpv(k),vmin,vmax
>          enddo
>          call thintl(latm3,merm3,levels,latpv,merpv,levpv,alevpv,
>     &         levcalc,theta,pcalc,eqcalc,dud)
>          do k=1,levcalcpv
>            do i=1,merpv
>              do j=1,latpv
>                grdwrk(irp(j,i,k),1) = eqcalc(j,i,k)
>              enddo
>            enddo
>            call minmax(eqcalc(1,1,k),merpv*latpv,vmin,vmax,dud,0)
>c     print *,'interp PV1 for PVG calc:',theta(k),vmin,vmax
>          enddo
>          do k=1,levpv
>            do i=1,merpv
>              do j=1,latpv
>                if(pvdata(irp(j,i,k),2).ne.dud)then
>                  pcalc(j,i,k) = pvdata(irp(j,i,k),2)*pvsign/stabil(k)*1.e2
>                else
>                  pcalc(j,i,k) = dud
>                endif
>              enddo
>            enddo
>            call minmax(pcalc(1,1,k),merpv*latpv,vmin,vmax,dud,0)
>c     print *,'scaled PV2 for PVG calc:',alevpv(k),vmin,vmax
>          enddo
>          call thintl(latm3,merm3,levels,latpv,merpv,levpv,alevpv,
>     &         levcalc,theta,pcalc,eqcalc,dud)
>          do k=1,levcalcpv
>            do i=1,merpv
>              do j=1,latpv
>                grdwrk(irp(j,i,k),2) = eqcalc(j,i,k)
>              enddo
>            enddo
>            call minmax(eqcalc(1,1,k),merpv*latpv,vmin,vmax,dud,0)
>c     print *,'interp PV2 for PVG calc:',theta(k),vmin,vmax
>          enddo
>          
>          do k=1,levcalcpv
>            do i=1,merpv
>              do j=1,latpv
>                hcalc1(j,i) = grdwrk(irp(j,i,k),1)
>                hcalc2(j,i) = grdwrk(irp(j,i,k),2)
>              enddo
>            enddo
>            call minmax(hcalc1,merpv*latpv,vmin1,vmax1,dud,0)
>            call minmax(hcalc2,merpv*latpv,vmin2,vmax2,dud,0)
>            merin = merpv
>            norm = .true.
>c     print *,'min, max PV1 for gradient:',theta(k),vmin1,vmax1
>            call hor_grad(merm3,merpv,latm3,latpv,xlatpv,xlonpv,
>     &           hcalc1,gcalc1,norm,dud)
>c     print *,'min, max PV2 for gradient:',theta(k),vmin2,vmax2
>            norm = .true.
>            call hor_grad(merm3,merpv,latm3,latpv,xlatpv,xlonpv,
>     &           hcalc2,gcalc2,norm,dud)
>            do i=1,merpv
>              do j=1,latpv
>                gdata(irp(j,i,k),1) = gcalc1(j,i)
>                gdata(irp(j,i,k),2) = gcalc2(j,i)
>              enddo
>            enddo
>          enddo
>          
>          call minmax(gdata(1,1),merpv*latpv*levcalc,vmin1,vmax1,dud,0)
>          call minmax(gdata(1,2),merpv*latpv*levcalc,vmin2,vmax2,dud,0)
>c     print *,'min, max nPVG:',vmin1,vmin2,vmax1,vmax2
>          
>          do k=1,levcalcpv
>            xing(1) = xpos(k)
>            ying(1) = ypos(k)
>            rmin = 1.e15
>            rmax = -1.e15
>            do i=1,merpv
>              do j=1,latpv
>                if(gdata(irp(latpv+1-j,i,k),1).ne.dud .and.
>     &               gdata(irp(latpv+1-j,i,k),2).ne.dud)then
>                  eint(i,j) = gdata(irp(latpv+1-j,i,k),1)*(1.-alpha) +
>     &                 gdata(irp(latpv+1-j,i,k),2)*alpha
>                  rmin = min(eint(i,j),rmin)
>                  rmax = max(eint(i,j),rmax)
>                else
>                  eint(i,j) = dud
>                endif
>              enddo
>            enddo
>c     print *,'PV grad for grid_pnts:',theta(k),xing,ying,rmin,rmax
>            call grid_pnts(merm3,latm3,merpv,latpv,eint,90.,0.,
>     &           1,zint,xing,ying,1,dud)
>            
>            if(zint(1).lt.dud .and. zint(1).ge.-888.)then !NaN should fail any loop?
>              delhpv(k) = zint(1)
>            else
>              delhpv(k) = dud
>            endif
>          enddo
>          print *,delhpv
>          
>          call minmax(delhpv,levcalcpv,vmin,vmax,dud,0)
>          print *,'delhpv:',vmin,vmax
>          
>c     vortex edge conditions
>          do k=1,levp
>            do i=1,merpv 
>              do j=1,latpv
>                if(wsdata(irp(j,i,k),index1).ne.dud .and. wsdata(irp(j,i,k),index2).ne.dud)then
>                  ecalc(j,i,k) = wsdata(irp(j,i,k),index1)*(1.-alpha) + wsdata(irp(j,i,k),index2)*alpha
>                else
>                  ecalc(j,i,k) = dud
>                endif
>              enddo
>            enddo
>            do i=1,merids
>              do j=1,lats
>                if(tmpdata(ir(j,i,k),index1).ne.dud .and. tmpdata(ir(j,i,k),index2).ne.dud)then
>                  pcalc(j,i,k) = tmpdata(ir(j,i,k),index1)*(p1000/alev(k))**(2./7.)*(1.-alpha) +
>     &                 tmpdata(ir(j,i,k),index2)*(p1000/alev(k))**(2./7.)*alpha
>                else
>                  pcalc(j,i,k) = dud
>                endif
>              enddo
>            enddo
>c     call minmax(ecalc(1,1,k),merpv*latpv,vmin,vmax,dud,0)
>c     print *,'windspeed1:',alev(k),vmin,vmax
>          enddo
>          if(merids.ne.merpv .or. lats.ne.latpv)then
>c     get T on PV grid for interpolation
>            zinterp = .false.
>            do k=1,levels
>              do i=1,merids
>                do j=1,lats
>                  hgtdata(ir(j,i,k),1) = pcalc(j,i,k) !using hgtdata for temp array
>                enddo
>              enddo
>            enddo 
>            zinterp = .false.
>            csfc = 'p'
>            call gridinterp(merids,lats,90.,0.,merpv,latpv,90.,0.,levels,alev(k),levels,
>     &           alev(k),zinterp,csfc,hgtdata(1,1),dud)
>            do k=1,levels
>              do i=1,merpv 
>                do j=1,latpv
>                  pcalc(j,i,k) = hgtdata(ir(j,i,k),1)
>                enddo
>              enddo
>            enddo 
>          endif
>          
>          do k=1,levcalcpv
>            if(k.eq.1)then
>              first = .true.
>            else
>              first = .false.
>            endif
>            call onthetp(latm3,merm3,levels,latpv,merpv,levp,alev,
>     &           theta(k),p1000,dud,pcalc,ecalc,gcalc1,
>     &           gcalc2,first,first,'l')
>            do i=1,merids
>              do j=1,lats
>                eqcalc(j,i,k) = gcalc1(j,i)
>              enddo
>            enddo
>c     call minmax(gcalc1,merpv*latpv,vmin,vmax,dud,0)
>c     print *,'interp windspeed1:',theta(k),vmin,vmax
>          enddo
>          
>          do k=1,levcalcpv
>            do i=1,merpv
>              do j=1,latpv
>                if(grdwrk(irp(j,i,k),1).ne.dud .and.
>     &               grdwrk(irp(j,i,k),2).ne.dud)then
>                  pvcalc(j,i,k) = grdwrk(irp(j,i,k),1)*(1.-alpha) +
>     &                 grdwrk(irp(j,i,k),2)*alpha
>                  pvin(j,i,k,1) = pvcalc(j,i,k)
>                else
>                  pvcalc(j,i,k) = dud
>                  pvin(j,i,k,1) = dud
>                endif
>              enddo
>            enddo
>          enddo
>          print *,'pvspace'
>          call pv_space(pmax,merm3,merm3,latm3,merpv,latpv,levcalc,
>     &         pmax,pvlev,pvcalc,eqcalc,pvcalc2,cpole,dud)
>          avgmer = 'm'
>          print *,'equiv_lat'
>          call equiv_lat(pmax,mlevmax,1,pmax,levcalc,1,pvlev,pvcalc2,
>     &         avgmer,merm3,latm3,merpv,latpv,dlatpv,pvin,
>     &         latm3,latpv,xlatpv,eqcalc2,cpole,dud)
>          do k=1,levcalcpv
>            do j=1,latpv
>              wseql(j,k) = 0.0
>              ng = 0
>              do i=1,merpv 
>                if(eqcalc2(j,i,k).ne.dud)then
>                  ng = ng + 1
>                  wseql(j,k) = wseql(j,k) + eqcalc2(j,i,k)
>                endif
>              enddo
>              if(ng.gt.0)then
>                wseql(j,k) = wseql(j,k)/float(ng)
>              else
>                wseql(j,k) = dud
>              endif
>            enddo
>c     call minmax(wseql(1,k),latpv,vmin,vmax,dud,0)
>c     print *,'wseql:',theta(k),vmin,vmax 
>          enddo
>c     call minmax(wseql,levcalc*latpv,vmin,vmax,dud,0)
>c     print *,'wseql:',vmin,vmax
>          do k=1,levcalcpv
>            do i=1,merpv
>              do j=1,latpv
>                eqcalc(j,i,k) = pvcalc(j,i,k)
>              enddo
>            enddo
>          enddo
>          call pv_space(pmax,merm3,merm3,latm3,merpv,latpv,levcalc,
>     &         pmax,pvlev,pvcalc,eqcalc,pvcalc2,cpole,dud)
>          avgmer = 'm'
>          call equiv_lat(pmax,mlevmax,1,pmax,levcalc,1,pvlev,pvcalc2,
>     &         avgmer,merm3,latm3,merpv,latpv,dlatpv,pvin,
>     &         latm3,latpv,xlatpv,eqcalc2,cpole,dud)
>          do k=1,levcalcpv
>            do j=1,latpv
>              ng = 0
>              pveql(j,k) = 0.0
>              do i=1,merpv 
>                if(eqcalc2(j,i,k).ne.dud)then
>                  ng = ng + 1
>                  pveql(j,k) = pveql(j,k) + eqcalc2(j,i,k)
>                endif
>              enddo
>              if(ng.gt.0)then
>                pveql(j,k) = pveql(j,k)/float(ng)
>              else
>                pveql(j,k) = dud
>              endif
>            enddo
>c     call minmax(pveql(1,k),latpv,vmin,vmax,dud,0)
>c     print *,'pveql:',theta(k),vmin,vmax 
>          enddo
>c     call minmax(pveql,levcalc*latpv,vmin,vmax,dud,0)
>c     print *,'pveql:',vmin,vmax
>          norm = .false.
>          if(tmptrop.ne.dud)then
>            if(dyntrop.ne.dud)then
>              tropmin = min(dyntrop,tmptrop)
>            else
>              tropmin = tmptrop
>            endif
>            kstart = max(5, (nint(tropmin)-2) ) !trop - 2 km
>          else
>            kstart = 5                  !or 5.5 km
>          endif
>          do k=1,kstart
>            eqledge(k) = 500.
>            eqlinner(k) = 500.
>            eqlouter(k) = 500.
>          enddo
>          eqmin1 = 1.e15
>          eqmax1 = -1.e15
>          eqmin2 = 1.e15
>          eqmax2 = -1.e15
>          eqmin3 = 1.e15
>          eqmax3 = -1.e15
>          do k=kstart,levcalcpv
>            if(theta(k).ge.345. .and. theta(k).ne.dud)then
>              if(ypos(k).lt.0.)then
>c     print *,'hemispheric arrays for SH',theta(k),jm35,jsp
>                do j=jm35,jsp,-1
>                  if(pveql(j,k).ne.dud)then
>                    pvhem(jm35+1-j,1) = -pveql(j,k)
>                  else
>                    pvhem(jm35+1-j,1) = dud
>                  endif
>                  wshem(jm35+1-j) = wseql(j,k)
>                enddo
>              else
>c     print *,'hemispheric arrays for NH',theta(k),j35,latpv,jsp
>                do j=j35,latpv-jsp+1
>                  pvhem(j-j35+1,1) = pveql(j,k)
>                  wshem(j-j35+1) = wseql(j,k)
>                enddo
>              endif
>              do j=neqhem-jsp-1,neqhem
>                if(pvhem(j,1).eq.dud .and. pvhem(j-1,1).ne.dud) 
>     &               pvhem(j,1) = pvhem(j-1,1)
>              enddo
>              norm = .true.
>c     print *,'calling latgrad for pv',theta(k),(pvhem(j,1),j=1,neqhem)
>              call lat_grad(1,neqhem,neqmax,xeqhem,pvhem,pvghem,norm,dud)
>              do j=1,neqhem
>                if(pvghem(j,1).ne.dud) then
>                  if(xeqhem(j).gt.80.)then
>                    filt = 1. - (xeqhem(j) - 80.)/10.
>                    pvghem(j,1) = pvghem(j,1)*filt
>                  endif
>                endif
>                if(wshem(j).ne.dud) then
>                  if(xeqhem(j).gt.80.)then
>                    filt = 1. - (xeqhem(j) - 80.)/10.
>                    wshem(j) = wshem(j)*filt
>                  endif
>                endif
>                if(pvghem(j,1).ne.dud .and. wshem(j).ne.dud)then
>                  wsgrd(j,1) = wshem(j)*pvghem(j,1)
>                else
>                  wsgrd(j,1) = dud
>                endif
>              enddo
>              norm = .false.
>c     use curvature of windspeed x pvgrad rather than pvgrad
>c     print *,'calling latgrad for wsgrd',theta(k),(wsgrd(j,1),j=1,neqhem)
>              call lat_grad(1,neqhem,neqmax,xeqhem,wsgrd,pvg2hem,norm,dud)
>              wmax = -1.e15
>              wsgmax = dud
>              wsmax = dud
>              gmax = -1.e15
>              wgmax = -1.e15 
>              pvgmax = dud
>              eqlmax = dud
>              do j=1,neqhem-1
>                if(wshem(j).ne.dud)then
>                  wmax = max(wmax,wshem(j))
>                  if(wmax.eq.wshem(j))then
>                    wsmax = xeqhem(j)
>                    jwmax = j
>                  endif
>                endif
>                if(wsgrd(j,1).ne.dud)then
>                  wgmax = max(wgmax,wsgrd(j,1))
>                  if(wgmax.eq.wsgrd(j,1))then
>                    wsgmax = xeqhem(j)
>                    jwgmax = j
>                  endif
>                endif
>                if(pvghem(j,1).ne.dud)then
>                  gmax = max(gmax,pvghem(j,1))
>                  if(gmax.eq.pvghem(j,1))then
>                    pvgmax = xeqhem(j)
>                    jgmax = j
>                  endif
>                endif
>              enddo
>c     print *,jwmax,wsmax,jgmax,pvgmax,jwgmax,wgmax
>              if( ((jwgmax.le.jgmax .and. jwgmax.ge.jwmax-1) .or. 
>     &             (jwgmax.ge.jgmax .and. jwgmax.le.jwmax+1) .or.
>     &             (jwgmax.le.jwmax .and. jwgmax.ge.jgmax-1) .or.
>     &             (jwgmax.ge.jwmax .and. jwgmax.le.jgmax+1)) .and. 
>     &             (abs(jwmax-jwgmax).le.j10 .or. abs(jgmax-jwgmax).le.j10) )then
>c     use windspeed x pvgrd as vortex edge criterion
>                eqlmax = wsgmax
>                jeqmax = jwgmax
>              else
>                print *,theta(k),'**wind and pvgrd maxima do not match**'
>                eqlmax = 500.
>              endif
>              if((wshem(jwmax) - wshem(1)) .lt. 4.0)then
>                print *,'*****windspeed max near boundary, ',theta(k),'****'
>                eqlmax = 500.
>              endif
>c     print *,'pvghem:',(pvghem(j,1),j=1,neqhem)
>c     print *,'pvg2hem:',(pvg2hem(j,1),j=1,neqhem)
>c     print *,'wshem:',(wshem(j),j=1,neqhem)
>c     print *,'wsgrd:',(wsgrd(j,1),j=1,neqhem)
>c     print *,'maxima',pvgmax,wsmax,wsgmax,jgmax,jwmax,jwgmax
>c     print *,' eqlmax:',theta(k),eqlmax,pvghem(jeqmax,1),wshem(jwmax),wshem(1)
>              if(abs(eqlmax).gt.80. .or. wshem(jwmax).lt.15.2 .or. 
>     &             pvghem(jgmax,1).lt.1.1)then 
>                print *,'****wind too weak, vortex too small, or PVgrad too weak****'
>                print *,theta(k),wshem(jwmax),eqlmax,pvghem(jgmax,1)
>                eqlmax = 500.
>              endif
>              if(eqlmax.ne.dud .and. eqlmax.ne.500.)then
>                j = jeqmax + 1
>                do while(pvg2hem(j,1).lt.pvg2hem(j-1,1) .and. j.le.neqhem)
>                  j = j + 1
>                enddo
>                if(j.ne.neqhem .or. pvg2hem(j,1).ge.pvg2hem(j-1,1))then
>                  eqlin = xeqhem(j-1)
>                else
>                  eqlin = 500.
>                endif
>                j = jeqmax - 1
>                do while(pvg2hem(j,1).gt.pvg2hem(j+1,1) .and. j.ge.1)
>                  j = j - 1
>                enddo
>                if(j.ne.1 .or. pvg2hem(j,1).le.pvg2hem(j+1,1))then
>                  eqlou = xeqhem(j+1)
>                else
>                  eqlou = 500.
>                endif
>              else
>                eqlin = eqlmax
>                eqlou = eqlmax
>              endif
>              if(eqlmax.eq.500.)then
>                eqledge(k) = 500.
>                eqlinner(k) = 500.
>                eqlouter(k) = 500.
>              elseif(eqlmax.eq.dud)then
>                eqledge(k) = dud
>                eqlinner(k) = dud 
>                eqlouter(k) = dud 
>              else
>                eqledge(k) = abs(eqlout(k)) - eqlmax
>                if(eqlin.ne.500.)then
>                  eqlinner(k) = abs(eqlout(k)) - eqlin
>                else
>                  eqlinner(k) = 500.
>                endif
>                if(eqlou.ne.500.)then
>                  eqlouter(k) = abs(eqlout(k)) - eqlou
>                else
>                  eqlouter(k) = 500.
>                endif
>              endif
>            else                        !don't define vortex below 345K or if theta undefined
>              if(theta(k).lt.345.)then
>                eqlouter(k) = 500.
>                eqledge(k) = 500.
>                eqlinner(k) = 500. 
>              else
>                eqlouter(k) = dud 
>                eqledge(k) = dud 
>                eqlinner(k) = dud  
>              endif
>            endif
>          enddo
>          do k=2,levcalc-1
>            if(eqledge(k).ne.500. .and. eqledge(k).ne.dud .and.
>     &           eqledge(k-1).eq.500. .and. eqledge(k+1).eq.500.)then
>              eqledge(k) = 500.
>              eqlinner(k) = 500.
>              eqlouter(k) = 500.
>              print *,'****vortex set undefined at isolated level, ',theta(k),'****'
>            endif
>            print *,'VE:',theta(k),eqlout(k),eqlouter(k),
>     &           eqledge(k),eqlinner(k)
>            if(eqledge(k).ne.500. .and. eqledge(k).ne.dud)then
>              eqmin1 = min(eqmin1,eqledge(k))
>              eqmax1 = max(eqmax1,eqledge(k))
>            endif
>            if(eqlinner(k).ne.500. .and. eqlinner(k).ne.dud)then
>              eqmin2 = min(eqmin2,eqlinner(k))
>              eqmax2 = max(eqmax2,eqlinner(k))
>            endif
>            if(eqlouter(k).ne.500. .and. eqlouter(k).ne.dud)then
>              eqmin3 = min(eqmin3,eqlouter(k))
>              eqmax3 = max(eqmax3,eqlouter(k))
>            endif
>          enddo
>          
>        else                            ! bad data for PV-related calcs
>          
>          do k=1,levcalc
>            pvout(k) = dud
>            spvout(k) = dud
>            eqlout(k) = dud
>            delhpv(k) = dud
>            eqledge(k) = dud
>            eqlouter(k) = dud
>            eqlinner(k) = dud
>          enddo
>          dyntrop = dud
>          
>        endif
>        
>        do k=levcalcpv+1,levcalc
>          pvout(k) = dud
>          spvout(k) = dud
>          eqlout(k) = dud
>          delhpv(k) = dud
>          eqledge(k) = dud
>          eqlouter(k) = dud
>          eqlinner(k) = dud
>        enddo
>        
>        call minmax(pvout,levcalc,vmin,vmax,dud,0)
>        print *,'pvout:',vmin,vmax
>        call minmax(spvout,levcalc,vmin,vmax,dud,0)
>        print *,'spvout:',vmin,vmax
>        print *,'dyntrop:',dyntrop
>        call minmax(eqlout,levcalc,vmin,vmax,dud,0)
>        print *,'eqlout:',vmin,vmax
>        call minmax(delhpv,levcalc,vmin,vmax,dud,0)
>        print *,'delhpv:',vmin,vmax
>        print *,'EqL outer, edge, inner:',eqmin3,eqmax3,eqmin1,
>     &       eqmax1,eqmin2,eqmax2
>        
>c     horizontal T gradient
>c     print *,'first p interp'
>        do k=1,levp
>          do i=1,merids
>            do j=1,lats
>              pcalc(j,i,k) = tmpdata(ir(j,i,k),index1)
>            enddo
>          enddo
>c     call minmax(pcalc(1,1,k),lats*merids,vmin,vmax,dud,0)
>c     print *,'T1 for grad:',alev(k),vmin,vmax
>        enddo
>        do k=1,levcalc
>          if(press(k).ne.dud)then
>            zint(k) = press(k)
>          else
>            zint(k) = dud
>          endif
>        enddo
>        cslin = 'l'
>c     print *,'calling psplint',(alev(k),k=1,levels)
>        call psplint(latm3,merm3,levels,lats,merids,levp,alev,
>     &       levcalc,zint,pcalc,eqcalc,dud,cslin)
>        do k=1,levcalc
>          do i=1,merids
>            do j=1,lats
>              grdwrk(ir(j,i,k),1) = eqcalc(j,i,k)
>            enddo
>          enddo
>c     call minmax(eqcalc(1,1,k),lats*merids,vmin,vmax,dud,0)
>c     print *,zint(k),vmin,vmax
>        enddo
>        print *,'second p interp'
>        do k=1,levp
>          do i=1,merids
>            do j=1,lats
>              pcalc(j,i,k) = tmpdata(ir(j,i,k),index2)
>            enddo
>          enddo
>c     call minmax(pcalc(1,1,k),lats*merids,vmin,vmax,dud,0)
>c     print *,'T2 for grad:',alev(k),vmin,vmax
>        enddo
>        call psplint(latm3,merm3,levels,lats,merids,levp,alev,
>     &       levcalc,zint,pcalc,eqcalc,dud,cslin)
>        do k=1,levcalc
>          do i=1,merids
>            do j=1,lats
>              grdwrk(ir(j,i,k),2) = eqcalc(j,i,k)
>            enddo
>          enddo
>c     call minmax(eqcalc(1,1,k),lats*merids,vmin,vmax,dud,0)
>c     print *,zint(k),vmin,vmax
>        enddo
>        
>c     print *,'k loop for hor_grad'
>        do k=1,levcalc
>          do i=1,merids
>            do j=1,lats
>              hcalc1(j,i) = grdwrk(ir(j,i,k),1)
>              hcalc2(j,i) = grdwrk(ir(j,i,k),2)
>            enddo
>          enddo
>          norm = .false.
>          call hor_grad(merm3,merids,latm3,lats,xlat,xlon,
>     &         hcalc1,gcalc1,norm,dud)
>          call hor_grad(merm3,merids,latm3,lats,xlat,xlon,
>     &         hcalc2,gcalc2,norm,dud)
>          do i=1,merids
>            do j=1,lats
>              gdata(ir(j,i,k),1) = gcalc1(j,i)
>              gdata(ir(j,i,k),2) = gcalc2(j,i)
>            enddo
>          enddo
>        enddo
>c     call minmax(gdata(1,1),merids*lats*levcalc,vmin1,vmax1,dud,0)
>c     call minmax(gdata(1,2),merids*lats*levcalc,vmin2,vmax2,dud,0)
>c     print *,'min, max TG:',vmin1,vmin2,vmax1,vmax2
>        
>        do k=1,levcalc
>          xing(1) = xpos(k)
>          ying(1) = ypos(k)
>          rmin = 1.e15
>          rmax = -1.e15
>          do i=1,merids
>            do j=1,lats
>              if(gdata(ir(lats+1-j,i,k),1).ne.dud .and.
>     &             gdata(ir(lats+1-j,i,k),2).ne.dud)then
>                eint(i,j) = gdata(ir(lats+1-j,i,k),1)*(1.-alpha) +
>     &               gdata(ir(lats+1-j,i,k),2)*alpha
>              else
>                eint(i,j) = dud
>              endif
>              if(eint(i,j).ne.dud)then
>                rmin = min(eint(i,j),rmin)
>                rmax = max(eint(i,j),rmax)
>              endif
>            enddo
>          enddo
>c     if(k.eq.13)print *,'eint:',k,rmin,rmax
>          call grid_pnts(merm3,latm3,merids,lats,eint,90.,0.,
>     &         1,zint,xing,ying,1,dud)
>          delhtmp(k) = zint(1)
>        enddo
>        
>        call minmax(delhtmp,levcalc,vmin,vmax,dud,0)
>        print *,'delhtmp:',vmin,vmax
>        
>        if(dyntrop.eq.dud) dyntrop = -999.000
>        if(tmptrop.eq.dud) tmptrop = -999.000
>        do k=1,levcalc
>          if(pvout(k).eq.dud .or. eqlout(k).eq.dud)then
>            pvout(k) = dud
>            eqledge(k) = dud
>            eqlinner(k) = dud
>            eqlouter(k) = dud
>            eqlout(k) = dud
>            delhpv(k) = dud
>          endif
>          if(xpos(k).eq.dud) xpos(k) = -999.000
>          if(ypos(k).eq.dud) ypos(k) = -999.000
>          if(tmpout(k).eq.dud) tmpout(k) = -999.000
>          if(press(k).eq.dud) press(k) = -999.000
>          if(tmpout(k).eq.dud) tmpout(k) = -999.000
>          if(theta(k).eq.dud) theta(k) = -999.000
>          if(delhtmp(k).eq.dud) delhtmp(k) = -999.000
>          if(hgtout(k).eq.dud) hgtout(k) = -999.000
>          if(uout(k).eq.dud) uout(k) = -999.000
>          if(vout(k).eq.dud) vout(k) = -999.000
>          if(pvout(k).eq.dud) pvout(k) = -999.000
>          if(spvout(k).eq.dud) spvout(k) = -999.000
>          if(eqlout(k).eq.dud) eqlout(k) = -999.000
>          if(delhpv(k).eq.dud) delhpv(k) = -999.000
>          if(eqledge(k).eq.dud) eqledge(k) = -999.000
>          if(eqlinner(k).eq.dud) eqlinner(k) = -999.000
>          if(eqlouter(k).eq.dud) eqlouter(k) = -999.000
>        enddo
>        
>        print *,'******Writing DMPs for Occ# ',aceocc,'******'
>
>cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>c     
>c     Open the outout DMP file
>c     
>cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
>
>
>        open(8,file=outfile,iostat=ios,status='replace')
>        if (ios.ne.0) then
>          print 1114,'### Error opening ', outfile, ' Aborting!'
>          call exit(1)
>        endif 
>
>
>        write(8,*)nametag
>        write(8,*)datetime
>        write(8,*)latin  
>        write(8,*)lonin  
>        write(8,*)locflag
>        write(8,*) dmp_version
>        write(8,*) dmp_date
>        write(8,*)
>
>        write(8,*)'Dyn-Trop WMO-Trop:'
>        write(8,1112)dyntrop,tmptrop
> 1112   format(2f8.2)
>        write(8,*)'z  ACE-T  ACE-p ACE-lon ACE-lat MetO-theta MetO-T delh-T GPH Zon-Wind Mer-Wind PV sPV EqL delh-PV EqL-edge EqL-inner EqL-outer'
>
>c '
>        do k=1,levcalc
>          write(8,1111)zalt(k),atemp(k),press(k),xpos(k),ypos(k),
>     &         theta(k),
>     &         tmpout(k),
>     &         delhtmp(k),
>     &         hgtout(k),
>     &         uout(k),
>     &         vout(k),
>     &         pvout(k),
>     &         spvout(k),
>     &         eqlout(k),
>     &         delhpv(k),
>     &         eqledge(k),
>     &         eqlinner(k),
>     &         eqlouter(k)
>        enddo
>        close(8)
>        
>c      else
>c        write(8,*)'NO ACE DATA...'
>
>      endif   
>      
> 1111 format(f6.1,17e12.4)
>      
>      rcsid = 
>     &     "$Id: metfields_ACEv2.2MetO.f,v 1.14 2007/01/22 20:50:35 whdaffer Exp $"
>      
>      end

Actions: View

Attachments on bug 281331: 189171