program yannos_ST_format800km_3_ellrot
c
c     The new format of the eigenfunction file being read has an additional
c     line with information on the portion of the mantle for which
c     eigenfunctions are stored. If the read eigenfunction set is not for 
c     the whole mantle, the code stops, because the complete set is needed
c     for the computation of the ellipticity and rotational splitting 
c     paramters. (Federica Apr 25, 2005)
c     Computation of the ellipticity and rotational splitting parameters
c     has been added. The value of these parameters is computed using the
c     C routine ellcoup_rot (in
c     /data/08/federica/SRC/KERNEL/ANISOZ/ellgravy_rot.c).
c     The computed parameters ellterm and rotterm have the same normalization
c     as the terms stored in the original XD eigenfuction files. Comments in
c     the routine ellcoup_rot and the readme
c     /data/08/federica/SRC/KERNEL/README_YANNOS give more details on the used
c     normalization. 
c     The values for ellterm and rotterm are stored for each mode in the 
c     eigenfunction header together with l,n,w,... For reading this reformatted
c     eigenfunction file into other codes use the modified routine
c     geteigys_ellrot.c   (Federica Feb 25, 2005)
c
c     Eigenfunctions for core and Stoneley computed inaccurately are set to 0. 
c     (Federica Feb 25, 2005)
c
c     further fix the n shift when l == 1 **** 2003/9/12 for 
c     both T and S modes
c     since there are cases that there is no shift in S, eg. for prem_228 
c
c     fix the n shift when l == 1 **** 2003/5/6
c
c     output top 800 km 10/2002 gung
c
c	!!! not sure about the sign of surface phi !!!
c
c
c     reformat the output bin eigfunction of yannos to be used by kernel
c     calculation and acal. 06/2002  gung
c
c     yannos output must be switched to UCB format (change yannos flag)
c  
c     input files 
c	1. lu_yannos.mdl : input model of yannos 
c	   (currently, only works for layer format. ifdeck=1) 
c	2. lu_yannos.asc ascii output of mode infromation from yannos 
c	3. lu_yannos.bin binary output of eigenfunction from yannos
c 
c
c
c     some convention in fortran
c     if (ACCESS=squential and FORM=unformatted ) ,each record is 
c     preceded and terminated with an integer*4 count, making each
c     record 8 characters longer than normal
c
c     --> it must be read the same way you wrote it
c
c     XD's structure of eigenfunction 
c     typedef struct { int n,type,l;
c     float w,q,cgp,a_vert,a_hori,phis;
c     float u,up,v,vp,ph,php; (array with the size of layers)
c     } eigen_st;
c
c
c     **   u,up,ph,php is not written for Toroidal mode
c     **   v,vp is 0 for Radial mode	
c
c	
c*****************************************************************
c     NORMALIZATION 	sql=1./sqrt(l*(l+1))
c     for radial mode, sql = 1
c
c     eig from (mineos, yannos)	eig used in kernel calculation
c     U			=	U   ( == U for radial mode)
c     Up		= 	Up  ( == Up for radial mode)
c     V			=	V*sql
c     Vp		=	Vp*sql
c     Phi		=	Phi
c     Phip		=	Phip
c     W			=	W*sql
c     Wp		=	Wp*sql
c	
c
c	a_vert = -U*w*w ( U is eig used in kernel calculation)
c       a_hori = -V*w*w ( V is eig used in kernel calculation, ie. normalized by sql already)      
c
c	 based on the option of iocean in Yann's minos model input
c	 a_vert,a_hori,and phis can be evaluated at the layer below
c	 seafloor or on the surface.  (nsfl)
c
c
c********************************************************************


      parameter ( nktmx=3000, mxlost=10000, ndimmx=nktmx*6+30,
     1     llmax=3000, ifinfo=1 )
      CHARACTER modetype(2),type

      common/mdl/r(nktmx),acon(nktmx),ccon(nktmx),lcon(nktmx)
     1     ,ncon(nktmx),fcon(nktmx),qshear(nktmx),qkappa(nktmx)  
     2     ,rho(nktmx),knot,nic,noc,nmoho,nsfl

      common/modef/nmax(2,0:llmax),nmin(2,0:llmax),istart(2,0:llmax),
     1     lmin(2),lmax(2),itotall(2),itotalmode(2),lostn(2,mxlost),
     2     lostl(2,mxlost)
     
      real*8 ww,qq,cg,BUF(ndimmx) !yannos output double precision
      real*4 ww1,qq1,cg1,buf1(ndimmx) !reformatted into real*4
      real*4 bufs(ndimmx),rm(nktmx) ! for top 800km
      
      real*4 ellterm,rotterm
      
      integer nvec(2),sizeofeig(2),nlost(2)
      integer nvecs(2),sizeofeigs(2),nshift_l1(2)
      integer nbceau,nbcou_lay
      
      data modetype/'T','S'/

      pi=3.1415926535
      radius=6371000.
      
      open(1,file='lu_yannos.mdl',form='formatted',status='old')
      open(4,FILE='formatted800km.eig',form='unformatted'
     1     ,ACCESS='sequential')

c---------------------------------------------------------------------
c     get knot from the input model of yannos
c
      call getmodel(1,ifinfo)
      close(1)
      
c-----------------------------------------------------------------
c     get useful mode information from ascii output information 
c     of yannos
c
      do jcom = 1,2
         type=modetype(jcom)
         open(2,file='lu_yannos.'//type//'.asc',
     1        form='formatted',status='old')
         if(ifinfo.eq.1)then
            if(jcom.eq.1) then
               print *, 'getting mode information for Toroidal modes'
            else if (jcom.eq.2) then
               print *, 'getting mode information for Spheroidal modes'
            endif   
         endif
         call getyannos_info(2,jcom,ifinfo,nshift_l1(jcom))
         close(2)
      
         open(2,file='lu_yannos.'//type//'.lost',
     1        form='formatted',status='old')
         call getyannos_lost(2,jcom,ifinfo,nlost(jcom))
         close(2)
         
         print *, 'lost modes in ', type, nlost(jcom)            

      enddo

c
c     writing header information 
c     sizeofeig based on the original structure of eigen_st (see geteig.c)
c     

      nvec(1)=knot*2
      nvec(2)=knot*6

c     Increased to include ellterm and rotterm
c      sizeofeig(1)=(nvec(1)+9)*4     
c      sizeofeig(2)=(nvec(2)+9)*4     
      sizeofeig(1)=(nvec(1)+11)*4     
      sizeofeig(2)=(nvec(2)+11)*4     

      do i=knot,1,-1
	if((radius-r(i)).ge.800000.) goto 21
      enddo
21    knot800=i
      knots=knot-knot800+1
      nsfls=nsfl-knot800+1 

      nvecs(1)=(knots)*2
      nvecs(2)=(knots)*6

c     Increased to include ellterm and rotterm
c      sizeofeigs(1)=(nvecs(1)+9)*4     
c      sizeofeigs(2)=(nvecs(2)+9)*4     
      sizeofeigs(1)=(nvecs(1)+11)*4     
      sizeofeigs(2)=(nvecs(2)+11)*4     


      write(4)(itotall(j),j=1,2),(itotalmode(j),j=1,2),
     1     knots,nsfls,(sizeofeigs(j),j=1,2),
     2     (r(i)/radius,i=knot800,knot),
     3     (lmin(j),j=1,2),(lmax(j),j=1,2),
     4     (nmin(1,l),l=0,lmax(1)),(nmin(2,l),l=0,lmax(2)),
     5     (nmax(1,l),l=0,lmax(1)),(nmax(2,l),l=0,lmax(2)),
     6     (istart(1,l),l=0,lmax(1)),(istart(2,l),l=0,lmax(2))
                     
      iheadersize=(8+knots+4+(lmax(1)+1+lmax(2)+1)*3)*4

      if(ifinfo.eq.1)then   
         print *, ' '
         print *, 'header size of outputeigenfunction file',iheadersize
         print *, 'size of eig structure for Toroidal mode',sizeofeig(1)
         print *, 'size of eig structure for Spheroidal mode',
     1        sizeofeig(2)
      endif
c
c     for S mode (L!=0)
c     buf(i) : U(r)              vertical
c     buf(i+knot) : dU(r)/dr 
c     buf(i+2*knot) : V(r)       horizontal 
c     buf(i+3*knot) : dV(r)/dr
c     buf(i+4*knot) : phi(r)
c     buf(i+5*knot) : dphi(r)/dr
c
c     for S mode and L=0, only u and du are ourput from yannos
c     buf(i) : U(r)              vertical
c     buf(i+knot) : dU(r)/dr 
c
c
c     for T mode 
c     buf(i) : W(r)
c     buf(i+knot) : dW(r)/dr 	
c 
      
      
      if(lmin(2).eq.0)then
         nrm = nmax(2,0)-nmin(2,0)+1 ! number of radial mode
      else
         nrm=0
      endif
      nrmcount=0

      do 1000 jcom = 1,2
         type=modetype(jcom)
         open(3,file='lu_yannos.'//type//'.bin',form='unformatted',
     1        access='sequential',status='old')

         read(3,end=901)nbceau,nbcou_lay,(rm(I),I=1,nbcou_lay)
         if ((rm(nbcou_lay)*radius).ne.r(knot)) then
	   print *, 'Inconsistencies in the model'
	   stop
	 endif
	 if ((nbcou_lay).ne.(knot)) then
	   print *, 'Inconsistencies in the model'
	   print *, knot, 'layers in the model'
	   print *, nbcou_lay, 'layers in eigenfunctions'
	   stop
	 endif
         
         kk=0
         nlst=1

 900     continue    
         if(jcom.eq.1) then
            read(3,end=901)nn,ll,ww,qq,cg,(buf(I),I=1,nvec(1))

         else 
            if(nrmcount.lt.nrm)then 
               read(3,end=901)nn,ll,ww,qq,cg,(buf(I),I=1,nvec(1))    
               nrmcount=nrmcount+1
            else
               read(3,end=901)nn,ll,ww,qq,cg,(buf(I),I=1,nvec(2))              
            endif
         endif               

         if(ll.eq.1.and.nshift_l1(jcom).ne.0)nn=nn+nshift_l1(jcom) ! nshift


        
         KK=KK+1
         ww1=ww            ! ouput real*4 
         qq1=1./qq
         cg1=cg         
         
         if(ll.eq.0)then
            sql=1
         else
            sql=1./sqrt(ll*(ll+1.))
         endif
c
c         write(*,*)kk, nn,ll,qq1,sql

	if(mod(kk,500).eq.0)write(*,'(i8,i6,a3,i6)')kk,nn,type,ll



         mvec=nvec(jcom)
         if(ll.eq.0.and.jcom.eq.2)then
            do jj=nvec(1)+1,nvec(2)
               buf(jj)=0.
            enddo
         endif
         

c     to be consistent with XD's eigenfunction convention

         if(jcom.eq.2.and.ll.ne.0)then               	  
            do i=1,knot*2
               buf1(i)=buf(i)                 
            enddo     
            i1=knot*2+1
            i2=knot*4

            do i=i2+1,mvec
               buf1(i)=buf(i)                 
            enddo 



         else

            i1=1
            i2=mvec
         endif         
       
         do i=i1,i2
            buf1(i)=buf(i)*sql
         enddo
         
c
c     retrive only the top 800km
c
         i1=knot800
         i2=knot
         j=0
         do i=i1,i2
            j=j+1
            bufs(j)=buf1(i) 
         enddo
         i1=i1+knot
         i2=i2+knot
         do i=i1,i2
            j=j+1
            bufs(j)=buf1(i) 
         enddo
         if(type.eq.'S') then
            do jj=1,4
               i1=i1+knot
               i2=i2+knot
               do i=i1,i2
                  j=j+1
                  bufs(j)=buf1(i) 
               enddo
            enddo
         endif
  
c     
c     assign mode type and renormalize surface eigenfunctions,
c     which are used in the source excitation calculation  
c
c
         
         ww2=ww1*ww1
         if(type.eq.'T')then            
            itype=2  
            avert=0.            
            phis=0.
            ahori=-buf1(nsfl)*ww2        
         else if(type.eq.'S') then
      
            avert=-buf1(nsfl)*ww2

            if(ll.eq.0)then
               ahori=0
               phis=0
               itype=1   
            else
               ahori=-buf1(knot*2+nsfl)*ww2
               phis=buf(knot*4+nsfl)
               itype=3
            endif   
         endif
         
	 ellterm=0.
	 rotterm=0.
	 	 
	 call ell(knot,nn,itype,ll,ww1,qq1,cg1,avert,ahori,phis,buf1,
     1        ellterm,ellrot)

c         write(*,'(i4,a2,i4,4f13.9)')nn,type,ll,ww1,qq1,ellterm,ellrot

c zero out the lost modes
c

         if(nlost(jcom).gt.0) then
               if(nn.eq.lostn(jcom,nlst).and.ll.eq.lostl(jcom,nlst)) 
     1	     then
                  do i3=1,nvecs(jcom)                  
                     bufs(i3)=0.
                  enddo
                  print *, ' get lost mode ', nlst, nn,type,ll
                  nlst=nlst+1                 
               endif
         endif
            

         write(4)nn,itype,ll,ww1,qq1,cg1,avert,ahori,phis,ellterm,
     1        ellrot,(bufs(I),I=1,nvecs(jcom))

         goto 900

 901  continue

      if ((nlost(jcom).ne.0).and.(nlst-1.ne.nlost(jcom))) then
          print *, ' ******************************************'
          print *, ' yannos_format : error 3 , inconsistent lostmodes'
          print *, ' ******************************************'
          stop
       endif

      if(kk.ne.itotalmode(jcom))then
         print *, ' total modes in lu_yannos.bin = ' ,kk
         print *, ' total modes in lu_yannos.asc = ' ,itotalmode(jcom)         
         print *, 'yannos_format : error 2'
         stop
      endif    
      close(3)
	
 1000  continue
       close(4)


      STOP	

      END


c------------------------------------
      subroutine getyannos_lost(luy,jcom,ifinfo,nlost)
c     get info of lost mode
c
      parameter ( llmax=3000 ,mxlost=10000 )
      common/modef/nmax(2,0:llmax),nmin(2,0:llmax),istart(2,0:llmax),
     1     lmin(2),lmax(2),itotall(2),itotalmode(2),lostn(2,mxlost),
     2     lostl(2,mxlost)

       character type

       i=1
 2     read(luy,*,end=3)lostn(jcom,i),type,lostl(jcom,i),wr,
     1      wmh,prd,grv,q,rayl

c       if(lostl(jcom,i).eq.1)lostn(jcom,i)=lostn(jcom,i)+1
       if(lostl(jcom,i).eq.1)lostn(jcom,i)=lostn(jcom,i)+nshift_l1


       if((type.eq.'T'.and.jcom.eq.2).or.
     1      (type.eq.'S'.and.jcom.eq.1))then
          print *, ' ******************************************'
          print*, 'getyannos_lost : error 1'
          print *, ' ******************************************'

          stop
       endif
       i=i+1       
       if(i.gt.mxlost) then
          print*, 'getyannos_lst: parameter mxlost too small'
          stop
       endif
       goto 2
 3     nlost=i-1
       return
       end

c------------------------------------
      subroutine getyannos_info(luy,jcom,ifinfo,nshift_l1)
c
c reading the ascii output of yannos to
c extract the useful header information for
c the reformat of it's output eigenfunction file
c
      parameter (llmax=3000,mxlost=10000)

      common/modef/nmax(2,0:llmax),nmin(2,0:llmax),istart(2,0:llmax),
     1     lmin(2),lmax(2),itotall(2),itotalmode(2),lostn(2,mxlost),
     2     lostl(2,mxlost)

      character c1*5,c2*2,type
   	      	
c skip the header	
1	read(luy,'(a5,a2)') c1,c2
	if ( (jcom.eq.2.and.c2(2:2).eq.'S') .or.
     1       (jcom.eq.1.and.c2(2:2).eq.'T')) then
	backspace(luy)
	goto 11
	endif

	goto 1
c now the mode information	
c
11	continue


        prdmn=9999.
        prdmx=-1.
        grvmn=9999.
        grvmx=-1.

	lmin(jcom)=9999
	lmax(jcom)=-1
	l1=-1
	nmode=0
        
        nshift_l1=0
        first=1

2	read(luy,*,end=3)n,type,l,wr,wmh,prd,grv,q,rayl

        if((type.eq.'T'.and.jcom.eq.2).or.
     1       (type.eq.'S'.and.jcom.eq.1))then
           print*, 'getyannos_info : error 1'
           stop
        endif

   	
	if(l.eq.1.and.first.eq.1) then
           first=0
  	   call l1_nshift(jcom,n,wr,nshift_l1)
        endif

        if(l.eq.1) n=n+nshift_l1   ! fixed n shift in Yannos 

	nmode=nmode+1
        prdmn=amin1(prdmn,prd)
        prdmx=amax1(prdmx,prd)
        grvmn=amin1(grvmn,grv)
        grvmx=amax1(grvmx,grv)

	if(nmode.eq.1)lmin(jcom)=l	
	if(l1.ne.-1.and.l.ne.l1)  nmax(jcom,l1)=n1
	  

	if(l.ne.l1) then
	  nmin(jcom,l)=n
c          write(*,*)l,'nmn',nmin(jcom,l)
	  istart(jcom,l)=nmode
	endif

	l1=l
	n1=n
	goto 2	
	
3	lmax(jcom)=l1
	nmax(jcom,l1)=n1

        if(ifinfo.eq.1)then
           print *, ' '
           print *, 'min L      = ',lmin(jcom)
           print *, 'max L      = ',lmax(jcom)
           
           print *, 'min period = ',prdmn
           print *, 'max period = ',prdmx
           print *, 'min gr vel = ',grvmn
           print *, 'max gr vel = ',grvmx        
        endif

	if(lmax(1).ge.llmax.or.lmax(2).ge.llmax) then
	   print *, ' Lmax too large, reset llmax'
	   stop
	endif


        if(jcom.eq.2)then
           do l=lmin(jcom),lmax(jcom)
              istart(jcom,l)=itotalmode(1)+istart(jcom,l)
           enddo
        endif

	nmode1=0
        if(ifinfo.eq.1) then
	write(*,101) 'L','Nmn','Nmx','#n','start'
	endif
	do l=lmin(jcom),lmax(jcom)
           if(ifinfo.eq.1.and.l.le.10) write(*,102)l,nmin(jcom,l),
     1	   nmax(jcom,l),nmax(jcom,l)-nmin(jcom,l)+1,
     2     istart(jcom,l)            
           nmode1=nmode1+(nmax(jcom,l)-nmin(jcom,l)+1)
	enddo

        if(ifinfo.eq.1)then
           print *, 'mode type   = ',type
           print *, 'total modes = ',nmode
           print *, ' '
        endif
        
        if(nmode.ne.nmode1) then
           print*, 'getyannos_info : error 2'
	   print*, 'nmode=',nmode,'  nmode1=',nmode1
           stop
        endif

        itotall(jcom)=lmax(jcom)-lmin(jcom)+1
        itotalmode(jcom)=nmode

        if(jcom.eq.2.and.ifinfo.eq.1)
     1       print *, 'total modes (T+S) = ',itotalmode(1)+itotalmode(2)


101	format(/'information for l<=10'//5a9)
102	format(5i9)	

	return
        end



c---------------------------------------------------------
        subroutine getmodel(lumd,ifinfo)

        parameter ( nktmx=3000 )

        CHARACTER*80 mdltit
        common/mdl/r(nktmx),acon(nktmx),ccon(nktmx),lcon(nktmx)
     1       ,ncon(nktmx),fcon(nktmx),qshear(nktmx),qkappa(nktmx)  
     2       ,rho(nktmx),knot,nic,noc,nmoho,nsfl


      read(lumd,*)mdltit
      read(lumd,*)ifanis,tref,ifdeck
      if(ifdeck.eq.0) then
         print *, 'yannos_format currently works for deck model only'
         stop
      endif
      read(lumd,*) knot,nic,noc,nsfl

      nsfl=knot-nsfl
      
	

      if(knot.ge.nktmx) then
         print *, ' knot too large, reset nktmx'
         stop
      endif

c      
c     the following information is not used for formating eigenfunction.
c     could be put as header for double check in application code, when
c     getprem is also used.
c     
       if(ifanis.eq.0) read(lumd,1001)(r(i),rho(i)
     1      ,ccon(i),lcon(i),qkappa(i),qshear(i),i=1,knot)
 1001 format(f8.0,5f9.2)

       if(ifanis.ne.0) read(lumd,1002)(r(i),rho(i)
     1      ,ccon(i),lcon(i),qkappa(i),qshear(i),acon(i)
     2      ,ncon(i),fcon(i),i=1,knot)
 1002 format(f8.0,8f9.2)
      
      rmoho = 6371000.-24370.
      dr=1500.
      r1=rmoho-dr
      r2=rmoho+dr
      do i=1,knot
         if(qshear(i).ne.0)qshear(i)=1./qshear(i)
         if(qkappa(i).ne.0)qkappa(i)=1./qkappa(i)
         if(r(i).ge.r1.and.r(i).le.r2.and.r(i).eq.r(i+1))nmoho=i
      enddo

      if (ifinfo.eq.1) then
      print *, 'knot  = ',knot
      print *, 'nic   = ',nic
      print *, 'noc   = ',noc

      if( nmoho.ne.0) then
         print *, 'nmoho = ',nmoho
      else
         print *, 'nmoho = no moho disc '
      endif
      print *, 'nsfl  = ',nsfl
      print *, ' '
      endif

      return
      end


c---------------------------------------------------------
      subroutine l1_nshift(jcom,n,wr,nshift_l1)

      pi=3.1415926535
  
      T11=808                   ! period for T11 in XD prem222
      T21=457                   ! period for T12 in XD prem222
      
      S11=19511  
      S21=2475 
      S31=1059

      print*, 'yannos: l=1, n=',n

      if(jcom.eq.1) then
         dT1=abs(2*pi/wr-T11)/T11
         dT2=abs(2*pi/wr-T21)/T21       
         
         print*, 'look for n shift in T modes for l=1'
         
         if(dT1.lt.0.1) then
            nshift_l1=1-n
            print*, 'XD T11 =',T11, '   yannos = ',2*pi/wr            
         else if (dT2.lt.0.1) then
            nshift_l1=2-n
            print*, 'XD T21 =',T21, '   yannos = ',2*pi/wr 
         else 
            print*, '*****************************'
            print*, 'Warning:getyannos_info - check this!'
            print*, 'T l=1,n=',n, ' period= ',  2*pi/wr
            print*, '*****************************'            
         endif

         print*, 'nshift for Tl=1 = ',nshift_l1	     

      else if(jcom.eq.2) then
              
         dT1=abs(2*pi/wr-S11)/S11
         dT2=abs(2*pi/wr-S21)/S21             
         dT3=abs(2*pi/wr-S31)/S31    
         print*, 'look for n shift in S modes for l=1'
               
         if(dT1.lt.0.1) then
            nshift_l1=1-n
            print*, 'XD S11 =',S11, '   yannos = ',2*pi/wr 

         else if (dT2.lt.0.1) then
            nshift_l1=2-n
            print*, 'XD S21 =',S21, '   yannos = ',2*pi/wr

         else if (dT3.lt.0.1) then
            nshift_l1=3-n       
            print*, 'XD S31 =',S31, '   yannos = ', 2*pi/wr
         else
            print*, '*****************************'
            print*, 'Warning:getyannos_info - check this!'
            print*, 'S l=1,n=',n, ' period= ',  2*pi/wr
            print*, '*****************************'
         endif
         print*, 'nshift for l=1 = ',nshift_l1       	
      endif

      return
      end