program sunh2sp    
c converts sun hours monthly .nrm to sunshine percent (n/N)
      parameter (xmiss=-999.9,imiss=-9999)
      character*80 tdmfmt,tdmfmt2
      character*80 infl,outfl,bkfl
      character*80 ihd,jhd
      integer sunh(12)
      character nm*20,cn*13,id*9
      integer icode,ival
      integer dys(12)
      real day(12),dec(12),sunp(12)
c     
      DATA DYS/31,28.25,31,30,31,30,31,31,30,31,30,31/
      DATA DEC/-21.07,-13.4,-2.42,9.38,18.77,22.95,21.1,13.53,
     &      2.62,-9.08,-18.6,-23.35/                 
      pi=3.14159265
      rad=180.0/PI
c
      tdmfmt='(i7,i6,i7,i5,x,a20,x,a13,x,i4,x,i4,i7,a9,x,i2,x,i3)'
      tdmfmt2='(i4,12i5)'
c
      write(*,*)'Enter sunshine hours file name (.nrm)'
      read(*,'(a\)')infl
      open(1,file=infl,status='old')
      write(*,*)'Enter sunshine percent file name (.nrm) to write to'
      read(*,'(a\)')infl
      open(2,file=infl,status='unknown')
      
1     read(1,tdmfmt,end=99)iwmo,lat,lon,ielv,nm,cn,iy1,iy2,ibox,id,
     &      icode,ival
      write(2,tdmfmt)iwmo,lat,lon,ielv,nm,cn,iy1,iy2,ibox,id,icode,ival
      xlat=real(lat)/100
      do iy=1,1
       read(1,tdmfmt2)iyear,(sunh(im),im=1,12)
       DO im=1,12
        A1=-TAN(XLAT/RAD)*TAN(DEC(im)/RAD)
        IF(A1.GT.1.0) THEN 
         H=0
        ELSE
         IF(A1.LT.-1.0)THEN
          H=PI
         ELSE
          H=ACOS(A1)
         ENDIF
        ENDIF
        DAY(im)=(24*H)/PI
c
c Estimated from Table 3 in Doorenbus and Pruitt (1984)
c
        IF(DAY(im).NE.0) THEN
         DAY(im)=DAY(im)+0.1+(ABS(XLAT)*0.002)
         IF(DAY(im).GT.24)DAY(im)=24
         SUNp(im)=(100*((REAL(sunh(im))/10)/dys(im))/day(im))
         if(sunh(im).lt.0)sunp(im)=-999.9
        ELSE
	 sunp(im)=-999.9
        ENDIF
c						decide what to do when % is bad
        if(sunp(im).gt.100)sunp(im)=-999.9
        if(sunp(im).lt.  0)sunp(im)=-999.9
       enddo
       write(2,tdmfmt2)iyear,(nint(10*sunp(im)),im=1,12)
      enddo
      goto 1
99    end