module type_def

  type :: dim_type
    sequence
    integer(4) :: x,y,z
  end type dim_type

  type :: reg_type
    sequence
    integer(4) 		:: val
    character(20) 	:: name
    character(80) 	:: fn
  end type reg_type

end module type_def

module data_def

  use type_def

  integer(4), public 		:: n_reg, ct_level, val_land, val_coast, key_area
  real(4), public		:: area_scale, key_scale
  character(250), public  	:: fn_par,fn_area,fn_mask,path1,path2
  character(8), public  	:: dat_min, dat_max

  type (dim_type), public 	:: dim
  type (reg_type), public 	:: reg(100)

end module data_def
!
subroutine get_params_from_file()
!================================

  use data_def

  character(250) 	:: buf,arg,val,msg
  character(1)		:: ch1,ch2

  open(1,file=fn_par,status='old',form ='formatted',iostat=ios,iomsg=msg);
  if (ios.ne.0) then; write (*,*) trim(msg); return; end if; 
!  write (*,'(a,i4,2a)') 'ios: ',ios,' msg: ',trim(buf); read(*,*);

  i = 0;

  do while (.true.)
    read (1,'(a)',iostat=ios) buf; if (ios.ne.0) exit;
    l1 = len_trim(buf); l2 = index(buf,'=');
!    write (*,*) trim(buf),l1,l2; read(*,*);

    if (buf(1:1).eq.';'.or.l1.lt.3.or.l2.eq.0.or.l2.eq.l1) cycle
    if (buf(1:1).eq.' ') then
      do while (.true.)
        buf(1:l1-1) = buf(2:l1); l1 = l1 - 1;
        if (buf(1:1).ne.' ') exit;
      end do
      l2 = index(buf,'='); if (l2.le.l1) cycle;
    end if
    arg = buf(1:l2-1); l3 = len_trim(arg); val = buf(l2+1:l1);

    if     (arg(1:l3).eq.'area_file') then
      read (val,*) key_area,key_scale,fn_area
    elseif (arg(1:l3).eq.'mask_file') then
      fn_mask = val
    elseif (arg(1:l3).eq.'path_input') then
      path1 = val
    elseif (arg(1:l3).eq.'path_output') then
      path2 = val
    elseif (arg(1:l3).eq.'land') then
      read (val,*) ch1,ch2; val_land = ichar(ch1); val_coast = ichar(ch2);
    elseif (arg(1:l3).eq.'ct_level') then
      read (val,*) ct_level
    elseif (arg(1:l3).eq.'area_scale') then
      read (val,*) area_scale
    elseif (arg(1:l3).eq.'dimensions') then
      read (val,*) dim%x,dim%y
    elseif (arg(1:l3).eq.'region') then
      n_reg = n_reg + 1;  
      read (val,*) ch1,reg(n_reg)%fn
      if (ch1.eq.'0') then
        reg(n_reg)%val = 0
      else
        reg(n_reg)%val = ichar(ch1)
      end if
    else
      write (*,*) 'param file: ',trim(fn_par),' param: ',trim(buf)
    end if
  end do

  close (1)

  return

end subroutine get_params_from_file

!=========================================================================
subroutine init_params()
!=========================================================================

  use data_def

  character(250) 	:: arg,val
  character(1) 		:: ch1

  ct_level = 15;
  val_land = ichar('-'); val_coast = ichar('+');
  fn_area = 'area_12.dat'; fn_mask = 'mask_12.dat'; 
  path1 ='./'; path2 ='./';
  area_scale = 0.001;
  dim%x = 304; dim%y = 448;
!  dim%x = 608; dim%y = 896;
!  dim%x = 316; dim%y = 332;
!  dim%x = 632; dim%y = 664;

  i = 0;

  do while (.true.)
    i = i + 1; if (i.gt.iargc()) exit;
    call getarg(i, arg)
    i = i + 1; if (i.gt.iargc()) exit;
    call getarg(i, val)

    if (arg(1:2).eq.'-i'.or.arg(1:2).eq.'/i') then
      fn_par = val; call get_params_from_file();
    elseif (arg(1:2).eq.'-a'.or.arg(1:2).eq.'/a') then          
      fn_area = val
    elseif (arg(1:2).eq.'-p'.or.arg(1:2).eq.'/p') then          
      path1 = val
    elseif (arg(1:2).eq.'-P'.or.arg(1:2).eq.'/P') then          
      path2 = val
    elseif (arg(1:2).eq.'-m'.or.arg(1:2).eq.'/m') then          
      fn_mask = val
    elseif (arg(1:2).eq.'-l'.or.arg(1:2).eq.'/l') then          
      read (val,*) ct_level
    elseif (arg(1:2).eq.'-s'.or.arg(1:2).eq.'/s') then          
      read (val,*) area_scale
    elseif (arg(1:2).eq.'-d'.or.arg(1:2).eq.'/d') then          
      read (val,*) dat_min, dat_max
    elseif (arg(1:2).eq.'-g'.or.arg(1:2).eq.'/g') then          
      read (val,*) dim%x,dim%y
    elseif (arg(1:2).eq.'-r'.or.arg(1:2).eq.'/r') then
      n_reg = n_reg + 1; read (val,*) ch1,reg(n_reg)%fn;
      if (ch1.eq.'0') then; reg(n_reg)%val = 0;
      else; reg(n_reg)%val = ichar(ch1);
      end if
    end if
  end do

  return

end subroutine init_params

program ssmi_reg
!================

  use data_def

  integer(4) 			:: ddJD
  integer(1), allocatable 	:: buf(:,:), buf1(:,:), land_mask(:,:), region_nsidc(:,:)
  integer(1) 			:: header(300)
  integer(4), allocatable 	:: grid(:,:), mask(:,:)

  real(4), allocatable 		:: area(:,:)
  real(4) 			:: ct_ratio       
                            
  double precision 		:: hist(0:255,0:255),sum 
  double precision 		:: extent(10,0:255),cover(10,0:255) 

  character(250) 		:: fn, msg

  integer(4) 			:: ct
  integer(4) 			:: JD, JD1, JD2, iDD,iMM,iYY

  integer(4) 			:: coastline = 253    
  integer(4) 			:: landmask  = 254  
  integer(4) 			:: polehole  = 251  
  integer(4) 			:: missingdata = 255
  logical			:: is_file  

  call init_params();
 
  allocate (buf(dim%x,dim%y),buf1(dim%x+2,dim%y),land_mask(dim%x,dim%y),region_nsidc(dim%x,dim%y))
  allocate (grid(dim%x,dim%y),mask(dim%x,dim%y),area(dim%x,dim%y))

  open(unit=1,file=fn_area,form='unformatted',access='stream',status='old',iostat=ios,iomsg=msg);
  if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if;
  if (key_area.eq.1) then 
    read (1,iostat=ios,iomsg=msg) area; 
    if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if;
  elseif (key_area.eq.2) then
    read (1,iostat=ios,iomsg=msg) mask; 
    if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if;
    area = key_scale * mask;
  end if 
  close (1);

  open(unit=1,file=fn_mask,form='unformatted',access='stream',status='old',iostat=ios,iomsg=msg)
  if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
  read (1,iostat=ios,iomsg=msg) buf1; 
  if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
  close (1);

  do j = 1,dim%y
    do i = 1,dim%x;
!      mask(i,j) = transfer(buf1(i,j),1)
!      mask(i,j) = iand(mask(i,j),X'000000ff')
      mask(i,j) = ibits(buf1(i,j), 0, 8);
    end do
  end do

  iu = 10
  do ireg = 1,n_reg
    fn = reg(ireg)%fn;
    inquire (file=trim(path2)//fn,exist=is_file); 
    if (is_file) then
      open(unit=iu+ireg,file=trim(path2)//fn,status='old',action='readwrite',iostat=ios,iomsg=msg)
      if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
      read (dat_min,'(i4,2i2)') iYY,iMM,iDD;  JD1 = ddJD(iDD,iMM,iYY); nn = 0;
      do while (.true.)
        read (iu+ireg,'(i2,1x,i2,1x,i4)',iostat=ios) iDD,iMM,iYY; if (ios.ne.0) exit
        nn = nn + 1; JD2 = ddJD(iDD,iMM,iYY);
        if (JD2.ge.JD1) exit
      end do
      if (nn.gt.0) backspace (iu+ireg)
    else
      open(unit=iu+ireg,file=trim(path2)//fn,status='new',action='write',iostat=ios,iomsg=msg)
      if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
    end if
  end do

  read (dat_min,'(i4,2i2)') iYY,iMM,iDD;  JD1 = ddJD(iDD,iMM,iYY);
  read (dat_max,'(i4,2i2)') iYY,iMM,iDD;  JD2 = ddJD(iDD,iMM,iYY);

  do JD = JD1, JD2
    call JDdd(JD,iDD,iMM,iYY);
    write (fn(1:12),'(i4.4,2i2.2,a)') iYY,iMM,iDD,'.dat'
    inquire (file=trim(path1)//fn(1:12),exist=is_file); if (.not.is_file) cycle;
    open(unit=1,file=trim(path1)//fn(1:12),form='unformatted',access='stream',status='old',iostat=ios,iomsg=msg)
    if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
    read (1,iostat=ios,iomsg=msg) buf; 
    if (ios.ne.0) then; write (*,*) trim(msg); goto 900; end if; 
    close (1);

    do j = 1,dim%y 
      do i = 1,dim%x
        grid(i,j) = ibits(buf(i,j), 0, 8);
!        grid(i,j) = transfer(buf(i,j),1); grid(i,j) = iand(grid(i,j),X'000000ff')
      end do 
    end do

    hist = 0D0; extent = 0D0; cover = 0D0;

    do j = 1,dim%y
      do i = 1,dim%x
        ct = grid(i,j); k = mask(i,j); if (ct.eq.polehole) ct = 100;
        if (k.ne.val_land.and.k.ne.val_coast) then
          hist(ct,k) = hist(ct,k) + area(i,j);
        end if
      end do 
    end do

    do ireg = 1,n_reg
      k = reg(ireg)%val
      if (k.ne.0) then
        do i = 0,255
          hist(i,0) = hist(i,0) + hist(i,k)
        end do
      end if
    end do

    do ireg = 1,n_reg
      k = reg(ireg)%val
      do i = ct_level,100
        extent(1,k)=extent(1,k)+hist(i,k); cover(1,k)=cover(1,k)+hist(i,k)*float(i)*0.01
      end do
      do i = 5,34
        extent(2,k)=extent(2,k)+hist(i,k); cover(2,k)=cover(2,k)+hist(i,k)*float(i)*0.01
      end do
      do i = 35,64
        extent(3,k)=extent(3,k)+hist(i,k); cover(3,k)=cover(3,k)+ hist(i,k)*float(i)*0.01
      end do
      do i = 65,84
        extent(4,k)=extent(4,k)+hist(i,k); cover(4,k)=cover(4,k)+hist(i,k)*float(i)*0.01
      end do
      do i = 85,100
        extent(5,k)=extent(5,k)+hist(i,k); cover(5,k)=cover(5,k)+hist(i,k)*float(i)*0.01
      end do
      extent(:,k) = extent(:,k) * area_scale; cover(:,k) = cover(:,k) * area_scale;
      if (extent(1,k).gt.0) then 
        ct_ratio = cover(1,k)/extent(1,k);
      else 
        ct_ratio = 0. 
      end if
      write (iu+ireg,'(2(i2.2,a1),i4.4,1x,i6,5f8.1,f6.3,5f8.1,f8.1)') iDD,'.',iMM,'.',iYY,JD, &
      (extent(j,k),j=1,5),ct_ratio,(cover(j,k),j=1,5),hist(missingdata,k)*area_scale
    end do  
  end do

900 continue

end program ssmi_reg

integer(4) function ddJD(iDD,iMM,iYY)
!------------------------------------ 

  integer(4) mm(12) /31,28,31,30,31,30,31,31,30,31,30,31/
  integer(4) iDD,iMM,iYY
  integer(4) days_in_february
  integer(4) days_in_year

  ddJD = 0
  do i = 1900,iYY-1
    if (days_in_february(i).eq.29) then
      days_in_year = 366
    else 
      days_in_year = 365
    end if
    ddJD = ddJD + days_in_year
  end do
  mm(2) = days_in_february(iYY)
  do i = 1,iMM-1
    ddJD = ddJD + mm(i)
  end do
  ddJD = ddJD + iDD

  return
end function ddJD

subroutine JDdd(JD,iDD,iMM,iYY)
!------------------------------- 

  integer(4) mm(12) /31,28,31,30,31,30,31,31,30,31,30,31/
  integer(4) iDD,iMM,iYY,JD,JDx
  integer(4) days_in_february
  integer(4) days_in_year

  JDx = JD
  do i = 1900,2999
    iYY = i
    if (days_in_february(i).eq.29) then
      days_in_year = 366
    else 
      days_in_year = 365
    end if
    JDx = JDx - days_in_year
    if (JDx.le.0) then
      JDx = JDx + days_in_year
      exit
    end if
  end do
  mm(2) = days_in_february(iYY)
  do i = 1,12
    iMM = i; JDx = JDx - mm(i)
    if (JDx.le.0) then
      iDD = JDx + mm(i); exit;
    end if
  end do

  return

end subroutine JDdd

integer function days_in_february(iYY)
!-------------------------------------- 

  integer(4) iYY

  if (iYY-(iYY/4)*4.eq.0) then
    if (iYY-(iYY/100)*100.ne.0) then
      days_in_february = 29
    else
      if (iYY-(iYY/400)*400.ne.0) then
        days_in_february = 28
      else
        days_in_february = 29
      end if
    end if
  else
    days_in_february = 28   
  end if

  return
end function days_in_february