!0000000011111111112222222222333333333344444444445555555555666666666677777777778
!2345678901234567890123456789012345678901234567890123456789012345678901234567890
!*******************************************************************************

!***************** program to read extracted NOAA TDR SSMIS data ***************
!******* (data format is the same as the TDRs but without file header) *********
!*******************************************************************************

!******************************* scan declarations *****************************
!*********** these arrays are to store the individual 'scene' data *************
	integer*1 ishdr(36),iemp(60),imag(4320),ienv(1800)
	integer*1 ilas(1440),iuas(480),iaux(1456)

!**************************** scan header declarations *************************
	integer*2 jsdy,ishr,ismn,iscann
	integer*4 isyr
	real*4 ztime(4100)

!********************** imager scener declarations (imag) **********************
	real*4 zv91(180,4100),zh91(180,4100),zh150(180,4100)
	real*4 zh183c(180,4100),zh183b(180,4100),zh183a(180,4100)
	real*4 zlata(180,4100),zlona(180,4100)
	real*4 zlatb(180,4100),zlonb(180,4100)
	integer*2 iscim(180,4100)
	integer*1 isfim(180,4100),irnfl(180,4100)
	real*4 coeff_v91(180),coeff_h91(180),coeff_h150(180)
	real*4 coeff_h183c(180),coeff_h183b(180),coeff_h183a(180)

!******************** environmental scene declarations (ienv) ******************
	real*4 zv19(90,4100),zh19(90,4100),zv22(90,4100)
	real*4 zv37(90,4100),zh37(90,4100)
	real*4 zlatc(90,4100),zlonc(90,4100)
	real*4 zlatd(90,4100),zlond(90,4100)
	integer*1 iscie(90,4100)
	integer*1 isfie(90,4100)
	real*4 coeff_v19(90),coeff_h19(90),coeff_v22(90)
	real*4 coeff_v37(90),coeff_h37(90)

!***************** lower atmosphere scene declarations (ilas) ******************
	real*4 zch01(60,4100),zch02(60,4100)
	real*4 zch03(60,4100),zch04(60,4100)
	real*4 zch05(60,4100),zch06(60,4100)
	real*4 zch07(60,4100),zch24(60,4100)
	real*4 zlate(60,4100),zlone(60,4100)
	integer*2 iscil(60,4100)
	integer*2 isfil(60,4100)
	real*4 coeff_ch01(60),coeff_ch02(60)
	real*4 coeff_ch03(60),coeff_ch04(60)
	real*4 coeff_ch05(60),coeff_ch06(60)
	real*4 coeff_ch07(60),coeff_ch24(60)

!***************** upper atmosphere scene declarations (iuas) ******************
	real*4 zch19(30,4100),zch20(30,4100)
	real*4 zch21(30,4100),zch22(30,4100)
	real*4 zch23(30,4100)
	real*4 zlatf(30,4100),zlonf(30,4100)
	integer*2 isciu(30,4100)
	real*4 coeff_ch19(30),coeff_ch20(30)
	real*4 coeff_ch21(30),coeff_ch22(30)
	real*4 coeff_ch23(30)

!***************************** general declarations ****************************
	integer*2 iscan,isdy,ismo
	character*40 fileip


!*******************************************************************************
!*************************** start of main program *****************************
!*******************************************************************************

!********************* get the data filename as an argument ********************
	call getarg(1,fileip)

!******************* read in the scan correction coefficients ******************
	if(fileip(7:9).eq.'f16')isat=1
	if(fileip(7:9).eq.'f17')isat=2
	if(fileip(7:9).eq.'f18')isat=3
	call read_coeff(isat,coeff_v91,coeff_h91,coeff_h150,
     +coeff_h183a,coeff_h183b,coeff_h183c,
     +coeff_v19,coeff_h19,coeff_v22,coeff_v37,coeff_h37)

!****************************** read in the data ******************************
	open(unit=1,file=fileip,access='direct',status='old',
     +form='unformatted',recl=9592)
	do iscan=1,9999
	read(1,rec=iscan,err=999)ishdr,iemp,imag,ienv,ilas,iuas,iaux
	call proc_ihdr(ishdr,iscan,isyr,jsdy,ishr,ismn,iscann,ztime)
	ichr=int(ztime(iscan)/3600.0)
	icmn=int(ztime(iscan)/60.0)-float(ichr)*60
	zcsc=ztime(iscan)-float(ichr)*3600-float(icmn)*60
	call conv_jday(isyr,jsdy,ismo,isdy)
!	call proc_iemp(iemp,iscan)                 ! ignore the ephemeris data for now....
	call proc_imag(imag,iscan,zv91,zh91,zh150,zh183a,zh183b,zh183c,
     +zlata,zlona,zlatb,zlonb,iscim,isfim,irnfl,coeff_v91,coeff_h91,
     +coeff_h150,coeff_h183a,coeff_h183b,coeff_h183c)
	call proc_ienv(ienv,iscan,zv19,zh19,zv22,zv37,zh37,
     +zlatc,zlonc,zlatd,zlond,iscie,isfie,
     +coeff_v19,coeff_h19,coeff_v22,coeff_v37,coeff_h37)
	call proc_ilas(ilas,iscan,zch01,zch02,zch03,zch04,zch05,
     +zch06,zch07,zch24,zlate,zlone,iscil,isfil)
	call proc_iuas(iuas,iscan,zch19,zch20,zch21,zch22,zch23,
     +zlatf,zlonf,isciu)
!	call proc_iaux(iaux,iscan)                 ! ignore the aux. data for now.......
	enddo
999	close(unit=1)


!*******************************************************************************
!*************************     ENTER YOUR CODE HERE.... ************************
!*******************************************************************************



	end


!*******************************************************************************
!********************************  SUBROUTINES  ********************************
!*******************************************************************************

!*******************************************************************************
!*********************************** scan header *******************************
!*******************************************************************************
	subroutine proc_ihdr(ishdr,iscan,isyr,jsdy,ishr,ismn,iscann,ztime)
	integer*1 ishdr(36)
	integer*2 iscan,jsdy,ishr,ismn,iscann
	integer*4 itime,isyr
	real*4 ztime(4100)

	call conv_b4si4(ishdr(1),ishdr(2),ishdr(3),ishdr(4),isyr)
	call conv_b2si2(ishdr(5),ishdr(6),jsdy)
	ishr=ishdr(7)
	ismn=ishdr(8)
	call conv_b2si2(ishdr(11),ishdr(12),iscann)
	call conv_b4si4(ishdr(13),ishdr(14),ishdr(15),ishdr(16),itime)
	ztime(iscan)=float(itime)/1000.0

	return
	end

!*******************************************************************************
!******************************** Ephemeris scene ******************************
!*******************************************************************************
	subroutine proc_iemp(iemp,iscan)
	integer*1 iemp(60)
	integer*2 iscan

	return
	end

!*******************************************************************************
!********************************** imager scene *******************************
!*******************************************************************************
	subroutine proc_imag(imag,iscan,
     +zv91,zh91,zh150,zh183a,zh183b,zh183c,
     +zlata,zlona,zlatb,zlonb,iscim,isfim,irnfl,
     +coeff_v91,coeff_h91,coeff_h150,coeff_h183a,
     +coeff_h183b,coeff_h183c)
	integer*1 imag(4320)
	integer*2 iscan
	integer*2 ilata,ilona,ih150,ih183c,ih183b,ih183a
	integer*2 ilatb,ilonb,iv91,ih91
	real*4 zv91(180,4100),zh91(180,4100)
	real*4 zh150(180,4100),zh183c(180,4100)
	real*4 zh183b(180,4100),zh183a(180,4100)
	real*4 zlata(180,4100),zlona(180,4100)
	real*4 zlatb(180,4100),zlonb(180,4100)
	integer*2 iscim(180,4100)
	integer*1 isfim(180,4100),irnfl(180,4100)
	real*4 coeff_v91(180),coeff_h91(180)
	real*4 coeff_h150(180),coeff_h183c(180)
	real*4 coeff_h183b(180),coeff_h183a(180)

!***** currently commented out since coefficients are read in from file ********
!	do i=1,180                   ! assign dummy values to the coefficients
!	coeff_v91(i)=1.0  
!	coeff_h91(i)=1.0
!	coeff_h150(i)=1.0
!	coeff_h183c(i)=1.0
!	coeff_h183b(i)=1.0
!	coeff_h183a(i)=1.0
!	enddo

	do ipos=1,180                ! loop through the 180 scan positions
	ioff=24*(ipos-1)             ! each group is 24 bytes

	call conv_b2si2(imag(1+ioff),imag(2+ioff),ilata)
	call conv_b2si2(imag(3+ioff),imag(4+ioff),ilona)
	zlata(ipos,iscan)=float(ilata)*0.01
	zlona(ipos,iscan)=float(ilona)*0.01

	call conv_b2si2(imag(5+ioff),imag(6+ioff),iscim(ipos,iscan))
	isfim(ipos,iscan)=imag(7+ioff)
	irnfl(ipos,iscan)=imag(8+ioff)

	call conv_b2si2(imag(9+ioff),imag(10+ioff),ih150)
	call conv_b2si2(imag(11+ioff),imag(12+ioff),ih183c)
	call conv_b2si2(imag(13+ioff),imag(14+ioff),ih183b)
	call conv_b2si2(imag(15+ioff),imag(16+ioff),ih183a)
	zh150(ipos,iscan)=(float(ih150)*0.01+273.15)*coeff_h150(ipos)
	zh183c(ipos,iscan)=(float(ih183c)*0.01+273.15)*coeff_h183c(ipos)
	zh183b(ipos,iscan)=(float(ih183b)*0.01+273.15)*coeff_h183b(ipos)
	zh183a(ipos,iscan)=(float(ih183a)*0.01+273.15)*coeff_h183a(ipos)

	call conv_b2si2(imag(17+ioff),imag(18+ioff),ilatb)
	call conv_b2si2(imag(19+ioff),imag(20+ioff),ilonb)
	zlatb(ipos,iscan)=float(ilatb)*0.01
	zlonb(ipos,iscan)=float(ilonb)*0.01

	call conv_b2si2(imag(21+ioff),imag(22+ioff),iv91)
	call conv_b2si2(imag(23+ioff),imag(24+ioff),ih91)
	zv91(ipos,iscan)=(float(iv91)*0.01+273.15)*coeff_v91(ipos)
	zh91(ipos,iscan)=(float(ih91)*0.01+273.15)*coeff_h91(ipos)

	enddo

	return
	end


!*******************************************************************************
!****************************** environmental scene ****************************
!*******************************************************************************
	subroutine proc_ienv(ienv,iscan,zv19,zh19,zv22,zv37,zh37,
     +zlatc,zlonc,zlatd,zlond,iscie,isfie,
     +coeff_v19,coeff_h19,coeff_v22,coeff_v37,coeff_h37)
	integer*1 ienv(1800)
	integer*2 iscan
	integer*2 ilatc,ilonc,iv19,ih19,iv22
	integer*2 ilatd,ilond,iv37,ih37
	real*4 zv19(90,4100),zh19(90,4100)
	real*4 zv22(90,4100)
	real*4 zv37(90,4100),zh37(90,4100)
	real*4 zlatc(90,4100),zlonc(90,4100)
	real*4 zlatd(90,4100),zlond(90,4100)
	integer*1 iscie(90,4100)
	integer*1 isfie(90,4100)
	real*4 coeff_v19(90),coeff_h19(90)
	real*4 coeff_v22(90)
	real*4 coeff_v37(90),coeff_h37(90)

!***** currently commented out since coefficients are read in from file ********
!	do i=1,90                   ! assign dummy values to the coefficients
!	coeff_v19(i)=1.0
!	coeff_h19(i)=1.0
!	coeff_v22(i)=1.0
!	coeff_v37(i)=1.0
!	coeff_h37(i)=1.0
!	enddo

	do ipos=1,90                 ! loop through the 90 scan positions
	ioff=20*(ipos-1)             ! each group is 20 bytes

	call conv_b2si2(ienv(1+ioff),ienv(2+ioff),ilatc)
	call conv_b2si2(ienv(3+ioff),ienv(4+ioff),ilonc)
	zlatc(ipos,iscan)=float(ilatc)*0.01
	zlonc(ipos,iscan)=float(ilonc)*0.01

	iscie(ipos,iscan)=ienv(5+ioff)
	isfie(ipos,iscan)=ienv(6+ioff)

	call conv_b2si2(ienv(9+ioff),ienv(10+ioff),iv19)
	call conv_b2si2(ienv(7+ioff),ienv(8+ioff),ih19)
	call conv_b2si2(ienv(11+ioff),ienv(12+ioff),iv22)
	zv19(ipos,iscan)=(float(iv19)*0.01+273.15)*coeff_h19(ipos)
	zh19(ipos,iscan)=(float(ih19)*0.01+273.15)*coeff_v19(ipos)
	zv22(ipos,iscan)=(float(iv22)*0.01+273.15)*coeff_v22(ipos)

	call conv_b2si2(ienv(13+ioff),ienv(14+ioff),ilatd)
	call conv_b2si2(ienv(15+ioff),ienv(16+ioff),ilond)
	zlatd(ipos,iscan)=float(ilatd)*0.01
	zlond(ipos,iscan)=float(ilond)*0.01
	call conv_b2si2(ienv(19+ioff),ienv(20+ioff),iv37)
	call conv_b2si2(ienv(17+ioff),ienv(18+ioff),ih37)
	zv37(ipos,iscan)=(float(iv37)*0.01+273.15)*coeff_v37(ipos)
	zh37(ipos,iscan)=(float(ih37)*0.01+273.15)*coeff_h37(ipos)

	enddo

	return
	end


!*******************************************************************************
!************************************ LAS scene ********************************
!*******************************************************************************
	subroutine proc_ilas(ilas,iscan,zch01,zch02,zch03,zch04,zch05,
     +zch06,zch07,zch24,zlate,zlone,iscil,isfil)
	integer*1 ilas(1440)
	integer*2 iscan
	integer*2 ilate,ilone
	integer*2 ich01,ich02,ich03,ich04,ich05,ich06,ich07,ich24
	real*4 zch01(60,4100),zch02(60,4100)
	real*4 zch03(60,4100),zch04(60,4100)
	real*4 zch05(60,4100),zch06(60,4100)
	real*4 zch07(60,4100),zch24(60,4100)
	real*4 zlate(60,4100),zlone(60,4100)
	integer*2 iscil(60,4100)
	integer*2 isfil(60,4100)
	real*4 coeff_ch01(60),coeff_ch02(60)
	real*4 coeff_ch03(60),coeff_ch04(60)
	real*4 coeff_ch05(60),coeff_ch06(60)
	real*4 coeff_ch07(60),coeff_ch24(60)

!***** currently commented out since coefficients are read in from file ********
!	do i=1,60                   ! assign dummy values to the coefficients
!	coeff_ch01(i)=1.0
!	coeff_ch02(i)=1.0
!	coeff_ch03(i)=1.0
!	coeff_ch04(i)=1.0
!	coeff_ch05(i)=1.0
!	coeff_ch06(i)=1.0
!	coeff_ch07(i)=1.0
!	coeff_ch24(i)=1.0
!	enddo

	do ipos=1,60                 ! loop through the 90 scan positions
	ioff=24*(ipos-1)             ! each group is 20 bytes

	call conv_b2si2(ilas(1+ioff),ilas(2+ioff),ilate)
	call conv_b2si2(ilas(3+ioff),ilas(4+ioff),ilone)
	zlate(ipos,iscan)=float(ilate)*0.01
	zlone(ipos,iscan)=float(ilone)*0.01

	call conv_b2si2(ilas(5+ioff),+ilas(6+ioff),iscil(ipos,iscan))
	call conv_b2si2(ilas(7+ioff),ilas(8+ioff),isfil(ipos,iscan))

	call conv_b2si2(ilas(9+ioff),ilas(10+ioff),ich01)
	call conv_b2si2(ilas(11+ioff),ilas(12+ioff),ich02)
	call conv_b2si2(ilas(13+ioff),ilas(14+ioff),ich03)
	call conv_b2si2(ilas(15+ioff),ilas(16+ioff),ich04)
	call conv_b2si2(ilas(17+ioff),ilas(18+ioff),ich05)
	call conv_b2si2(ilas(19+ioff),ilas(20+ioff),ich06)
	call conv_b2si2(ilas(21+ioff),ilas(22+ioff),ich07)
	call conv_b2si2(ilas(23+ioff),+ilas(24+ioff),ich24)

	zch01(ipos,iscan)=(float(ich01)*0.01+273.15)*coeff_ch01(ipos)
	zch02(ipos,iscan)=(float(ich02)*0.01+273.15)*coeff_ch02(ipos)
	zch03(ipos,iscan)=(float(ich03)*0.01+273.15)*coeff_ch03(ipos)
	zch04(ipos,iscan)=(float(ich04)*0.01+273.15)*coeff_ch04(ipos)
	zch05(ipos,iscan)=(float(ich05)*0.01+273.15)*coeff_ch05(ipos)
	zch06(ipos,iscan)=(float(ich06)*0.01+273.15)*coeff_ch06(ipos)
	zch07(ipos,iscan)=(float(ich07)*0.01+273.15)*coeff_ch07(ipos)
	zch24(ipos,iscan)=(float(ich24)*0.01+273.15)*coeff_ch24(ipos)

	enddo

	return
	end


!*******************************************************************************
!*********************************** UAS scene *********************************
!*******************************************************************************
	subroutine proc_iuas(iuas,iscan,zch19,zch20,zch21,zch22,zch23,
     +zlatf,zlonf,isciu)
	integer*1 iuas(480)
	integer*2 iscan,ilatf,ilonf
	integer*2 ich19,ich20,ich21,ich22,ich23
	real*4 zch19(30,4100),zch20(30,4100)
	real*4 zch21(30,4100),zch22(30,4100)
	real*4 zch23(30,4100)
	real*4 zlatf(30,4100),zlonf(30,4100)
	integer*2 isciu(30,4100)
	real*4 coeff_ch19(30),coeff_ch20(30)
	real*4 coeff_ch21(30),coeff_ch22(30)
	real*4 coeff_ch23(30)

!***** currently commented out since coefficients are read in from file ********
!	do i=1,30                   ! assign dummy values to the coefficients
!	coeff_ch19(i)=1.0
!	coeff_ch20(i)=1.0
!	coeff_ch21(i)=1.0
!	coeff_ch22(i)=1.0
!	coeff_ch23(i)=1.0
!	enddo

	do ipos=1,30                 ! loop through the 90 scan positions
	ioff=16*(ipos-1)             ! each group is 20 bytes

	call conv_b2si2(iuas(1+ioff),iuas(2+ioff),ilatf)
	call conv_b2si2(iuas(3+ioff),iuas(4+ioff),ilonf)
	zlatf(ipos,iscan)=float(ilatf)*0.01
	zlonf(ipos,iscan)=float(ilonf)*0.01

	call conv_b2si2(iuas(5+ioff),iuas(6+ioff),isciu(ipos,iscan))

	call conv_b2si2(iuas(9+ioff),iuas(10+ioff),ich19)
	call conv_b2si2(iuas(11+ioff),iuas(12+ioff),ich20)
	call conv_b2si2(iuas(13+ioff),iuas(14+ioff),ich21)
	call conv_b2si2(iuas(15+ioff),iuas(16+ioff),ich22)
	call conv_b2si2(iuas(17+ioff),iuas(18+ioff),ich23)

	zch19(ipos,iscan)=(float(ich19)*0.01+273.15)*coeff_ch19(ipos)
	zch20(ipos,iscan)=(float(ich20)*0.01+273.15)*coeff_ch20(ipos)
	zch21(ipos,iscan)=(float(ich21)*0.01+273.15)*coeff_ch21(ipos)
	zch22(ipos,iscan)=(float(ich22)*0.01+273.15)*coeff_ch22(ipos)
	zch23(ipos,iscan)=(float(ich23)*0.01+273.15)*coeff_ch23(ipos)

	enddo

	return
	end


!*******************************************************************************
!***************** section to read in the scan coefficients ********************
!*******************************************************************************
	subroutine read_coeff(isat,coeff_v91,coeff_h91,coeff_h150,
     +coeff_h183a,coeff_h183b,coeff_h183c,
     +coeff_v19,coeff_h19,coeff_v22,coeff_v37,coeff_h37)
	real*4 coeff_h150(180),coeff_h183a(180)
	real*4 coeff_h183b(180),coeff_h183c(180)
	real*4 coeff_v91(180),coeff_h91(180)
	real*4 coeff_v19(90),coeff_h19(90),coeff_v22(90)
	real*4 coeff_v37(90),coeff_h37(90)
	integer*2 isat
	character*80 fileip,junk

! isat defines the current satellite: F16= F17= F18=
! environmental channels are 12-16: h19, v19,v22,h37,v37
! imager scene channels are 8-11 and 17/18: h150, h183, h183, h183, v91, h91

! read in the environmental scene corrections first:
	if(isat.eq.1)fileip="ssmis_f16_env_scan_correction_coefficients.txt"
	if(isat.eq.2)fileip="ssmis_f17_env_scan_correction_coefficients.txt"
	if(isat.eq.3)fileip="ssmis_f18_env_scan_correction_coefficients.txt"
	open(unit=2,file=fileip,status='old')
	read(2,*)junk   ! reads in the test at the start of the file
	read(2,*)junk   ! reads in the test at the start of the file
	read(2,*)junk   ! reads in the test at the start of the file
	do i=1,90
	read(2,*)npos,coeff_h19(i),coeff_v19(i),coeff_v22(i),
     +coeff_h37(i),coeff_v37(i)
	enddo
	close(unit=1)

! read in the imager scene corrections first:
	if(isat.eq.1)fileip="ssmis_f16_img_scan_correction_coefficients.txt"
	if(isat.eq.2)fileip="ssmis_f17_img_scan_correction_coefficients.txt"
	if(isat.eq.3)fileip="ssmis_f18_img_scan_correction_coefficients.txt"
	open(unit=2,file=fileip,status='old')
	read(2,*)junk   ! reads in the test at the start of the file
	read(2,*)junk   ! reads in the test at the start of the file
	read(2,*)junk   ! reads in the test at the start of the file
	do i=1,180
	read(2,*)npos,coeff_h150(i),coeff_h183c(i),coeff_h183b(i),
     +coeff_h183a(i),coeff_v91(i),coeff_h91(i)
	enddo
	close(unit=1)

	return
	end	

!*******************************************************************************
!*********** subroutine to convert bytes to integers and byte swap *************
!*******************************************************************************

!******************************* integer*2 byte swap ***************************
	subroutine byteswapi2(k)
	integer*1 ii(2),jj(2)
	integer*2 i,j,k
	equivalence (i,ii)
	equivalence (j,jj)

	i=k
	jj(1)=ii(2)
	jj(2)=ii(1)
	k=j

	return
	end

!******************************* integer*4 byte swap ***************************
	subroutine byteswapi4(k)
	integer*1 ii(4),jj(4)
	integer*4 i,j,k
	equivalence (i,ii)
	equivalence (j,jj)

	i=k
	jj(1)=ii(4)
	jj(2)=ii(3)
	jj(3)=ii(2)
	jj(4)=ii(1)
	k=j

	return
	end


!******************************* real*4 byte swap ******************************
	subroutine byteswapr4(r)
	integer*1 ii(4),jj(4)
	real*4 r,s,t
	equivalence (s,ii)
	equivalence (t,jj)

	s=r
	jj(1)=ii(4)
	jj(2)=ii(3)
	jj(3)=ii(2)
	jj(4)=ii(1)
	r=t

	return
	end

!********************** 2 bytes to signed integer*2 ****************************
	subroutine conv_b2si2(i1,i2,ii)
	integer*1 i1,i2
	integer*2 ii,ii1,ii2
	ii1=i1
	ii2=i2	
	if(ii2.lt.0)ii2=ii2+256
	ii=ii1*256+ii2
	return
	end
	
!********************** 4 bytes to signed integer*4 ****************************
	subroutine conv_b4si4(i1,i2,i3,i4,ii)
	integer*1 i1,i2,i3,i4
	integer*4 ii,ii1,ii2,ii3,ii4
	ii1=i1
	ii2=i2
	ii3=i3
	ii4=i4
	if(ii2.lt.0)ii2=ii2+256
	if(ii3.lt.0)ii3=ii3+256
	if(ii4.lt.0)ii4=ii4+256
	ii=ii1*256*256*256+ii2*256*256+ii3*256+ii4
	return
	end

!*******************************************************************************
!************* subroutine to generate month/day from julian day ****************
	subroutine conv_jday(isyr,jsdy,ismo,isdy)
	integer*2 jsdy,ismo,isdy
	integer*4 isyr,ndays(12)

	data ndays /31,28,31,30,31,30,31,31,30,31,30,31/

	if(mod(isyr,4).eq.0)then
	ndays(2)=29
	endif

	ndy=jsdy
	do i=1,12
	if(ndy.ge.1.and.ndy.le.ndays(i))then
	isdy=ndy
	ismo=i
	endif
	ndy=ndy-ndays(i)
	enddo

	return
	end