C2GNUc
C
C An example of an external analysis function for VIS-5D version 4.0
C
C The parts in uppercase should be the same for all analysis functions,
C while the parts in lowercase are specific to this example (except
C comments).


C All VIS-5D analysis functions must be named USERFUNC and take the
C following arguments:
C
C   Output arguments:
C     OUTGRID - This is the array which the function computes
C     OUTNL - Number of levels in computed OUTGRID
C     OUTLOWLEV - Lowest level in computed OUTGRID
C
C   Input arguments:
C     INGRID -  Array of all 3-D grids for each available parameter
C     NR,NC - Number of rows and columns in all 3-D grids
C     NL - number of levels in *each* variable's 3-D grid.
C          An array of NL values is needed because each variable's
C          grid can contain a different number of levels.  However,
C          it's often the case that the number of grid levels is
C          the same for all grids.
C     LOWLEV - lowest level in *each* variable's 3-D grid.
C     MAXNL - maximum number of grid levels (max value in NL() array)
C     NVARS - number of physical variables
C     NAMES - names of each parameter
C     DATE - date of this grid in days since January 1, 1900
C            (use the IYYDDD function to convert to YYDDD format)
C     TIME - time of this grid in seconds since midnight
C            (use the IHMS function to convert to HHMMSS format)
C     PROJECTION - which map projection is being used.  See the README
C            file's sections on using the map projection and the
C            v5dCreate function to learn how this parameter and PROJ_ARGS
C            are used.
C     PROJ_ARGS - projection-specific parameters  (See README file)
C     VERTICAL - which vertical coordinate system is being used.  Again,
C            see the README file's documentation on v5dCreate to see how
C            this parameter and VERT_ARGS are used.
C     VERT_ARGS - vertical coordinate system specific parameters.
C
C Here is the row, column, level arrangement of the OUTGRID (and INGRID)
C arrays:
C     OUTGRID(1, 1, 1 ) = North-west-bottom corner
C     OUTGRID(NR,1, 1 ) = South-west-bottom corner
C     OUTGRID(1, NC,1 ) = North-east-bottom corner
C     OUTGRID(NR,NC,1 ) = South-east-bottom corner
C     OUTGRID(1, 1, NL) = North-west-top corner
C     OUTGRID(NR,1, NL) = South-west-top corner
C     OUTGRID(1, NC,NL) = North-east-top corner
C     OUTGRID(NR,NC,NL) = South-east-top corner
C
C The function should return 0 if it completed successfully or C a
C non-zero value if there was an error.
C
C
      INTEGER FUNCTION USERFUNC( OUTGRID, OUTNL, OUTLOWLEV,
     *                           INGRID, NR, NC, NL, LOWLEV, MAXNL,
     *                           NVARS, NAMES,
     *                           DATE, TIME,
     *                           PROJECTION, PROJ_ARGS,
     *                           VERTICAL, VERT_ARGS )
C     ARGUMENTS:
	implicit none
      INTEGER NVARS
      INTEGER NR, NC, NL(NVARS), LOWLEV(NVARS), MAXNL
      REAL OUTGRID(NR,NC,MAXNL)
      INTEGER OUTNL, OUTLOWLEV
      REAL INGRID(NR,NC,MAXNL,NVARS),zpo(nr,100),zp(nr,100)
      CHARACTER*8 NAMES(NVARS)
      INTEGER DATE, TIME
      INTEGER PROJECTION
      REAL PROJ_ARGS(*)
      INTEGER VERTICAL,IYYDDD,IHMS
      REAL VERT_ARGS(*),zpmin,zpmax,p(maxnl),pf(100),p1,p2
	character*20 name
	character*6 chartime
	character*6 chardate
	character*4 charlev


C     LOCAL VARS:
      integer iv, ir, ic, il, ifi,iln,numl,lev,leng,ilf,il1,il2,
     $        ilf_filt,ir_filt,inames,numnames,nfilt
      
C DAYS WIHIN THE 12 MONTHS (FOR NORMAL YEARS)
       Integer days_month(12),iyy,imm,idd,iyymmdd,julian
       Data days_month / 31,28,31,30,31,30,31,31,30,31,30,31 /          
          
C     Specify number of levels in OUTGRID
      OUTNL = MAXNL
      OUTLOWLEV = 0
C     look for variables
	do iln=1,maxnl
       p(iln)=1012.5*exp(vert_args(iln)/-7.2)
       enddo
	open(17,file='2GNUc.txt')
C Llegim numero de variables que volem passar a gnuplot
        read(17,*) numnames	
	Do 444 inames=1,numnames
	read(17,*)	
	
C Llegim el nom de la variable que volem passar a gnuplot
	read(17,'(A)')name
	
C Llegim nfilt: Rows and Columns on either side of center to apply matrix smooth 
C (average method with weight of matrix center = 2):	
	read(17,*) nfilt
		
C Cercam la longitud de name
	if (name(2:2).eq.'"') then
	leng=1	 
	else if (name(3:3).eq.'"') then
	leng=2	 
	else if (name(4:4).eq.'"') then
	leng=3	 
	else if (name(5:5).eq.'"') then
	leng=4	 
	else if (name(6:6).eq.'"') then
	leng=5	 
	else if (name(7:7).eq.'"') then
	leng=6	 
	else if (name(8:8).eq.'"') then
	leng=7	 
	else if (name(9:9).eq.'"') then
	leng=8	 
	else if (name(10:10).eq.'"') then
	leng=9	 
	endif	 
	do iv=1,nvars
	if (names(iv).eq.name(2:leng)) then
	ifi=iv
	endif
	enddo
	print*,'Variable: ',names(ifi),'nfilt=',nfilt
	read(17,'(I2)')numl
        do 333 iln=1,numl
        read(17,'(I4)')ic
	print*,'Column:',ic
	
	do 3000 ilf=1,100
	 pf(ilf)=p(1)+(ilf-1)*(p(maxnl)-p(1))/99.	
         do il=1,maxnl-1
          If (p(il).Ge.pf(ilf).And.p(il+1).Le.pf(ilf)) Then
	   il1=il 
	   p1=p(il)
	   il2=il+1
	   p2=p(il+1)
	   Goto 20
	  Endif
20        Continue	  
         enddo 
	 do ir=1,nr
	  zpo(ir,ilf)=ingrid(nr-ir+1,ic,il1,ifi)+(pf(ilf)-p1)*
     $    (ingrid(nr-ir+1,ic,il2,ifi)-ingrid(nr-ir+1,ic,il1,ifi))/
     $    (p2-p1)	  
          zp(ir,ilf)=0
	 enddo  
3000    continue	  

        zpmin=10000000
	zpmax=-10000000	
	do 5000 ilf=1+nfilt,100-nfilt
	do 5000 ir=1+nfilt,nr-nfilt
	  zp(ir,ilf)=zpo(ir,ilf) 
	  do ilf_filt=ilf-nfilt,ilf+nfilt
	  do ir_filt=ir-nfilt,ir+nfilt	  	   	
	    zp(ir,ilf)=zp(ir,ilf)+zpo(ir_filt,ilf_filt)
          enddo
	  enddo
	  zp(ir,ilf)=zp(ir,ilf)/((2.*nfilt+1)*(2.*nfilt+1)+1.)
	  if (zp(ir,ilf).lt.zpmin) zpmin=zp(ir,ilf)
	  if (zp(ir,ilf).gt.zpmax) zpmax=zp(ir,ilf)	 	
5000    continue
	

C Cercam la longitud de ic	 
	 if (ic.ge.1000) write (charlev,'(I4)')ic
	 if (ic.lt.1000.and.ic.ge.100) write (charlev,'(A,I3)')'0',ic
	 if (ic.lt.100.and.ic.ge.10) write (charlev,'(A,I2)')'00',ic
	 if (ic.lt.10.and.ic.ge.1) write (charlev,'(A,I1)')'000',ic
	
	 iyy=Mod(IYYDDD(date),100000)/1000
	 imm=0
	 idd=0
	 If (Mod(iyy-1964,4).Eq.0) days_month(2)=29
	 julian=Mod(IYYDDD(date),100000)-iyy*1000
	 Do imm=1,12
	  If (julian.Gt.days_month(imm)) Then 
	   julian=julian-days_month(imm)
	  Else
	   idd=julian
	   goto 3333
	  Endif
	 Enddo
3333     Continue
         iyymmdd=iyy*10000+imm*100+idd	
	 
	write(chardate,'(I6)') iyymmdd
	if (chardate(1:1).eq.' ') chardate(1:1)='0'
	if (chardate(2:2).eq.' ') chardate(2:2)='0'	
	write(chartime,'(I6)') IHMS(time)
	if (chartime(1:1).eq.' ') chartime(1:1)='0'
	if (chartime(2:2).eq.' ') chartime(2:2)='0'
	if (chartime(3:3).eq.' ') chartime(3:3)='0'
	if (chartime(4:4).eq.' ') chartime(4:4)='0'		
	print*,'gnufiles/GNUc/'//chardate//chartime(1:4)//'-GNU'
     $	//name(2:leng)//'-'//'c'//charlev//'.grd'
     	open(24,file='gnufiles/GNUc/'//chardate//chartime(1:4)//'-GNU'
     $  //name(2:leng)//'-'//'c'//charlev//'.grd')

 	do ir=3,nc-3
	WRITE (24,'(300E16.6)')(zp(ir,ic),ic=3,96)  
        enddo       
              
	close(24)

333     Continue	

444     Continue	
		
	close(17)
	USERFUNC = 0
      RETURN
      END