!
! 16/01/2013 
!
!=========================================================================
module type_def
!=========================================================================
  type :: src_type
    sequence
    integer(1)	:: mix,aar,bsi,cis,dmi,jma,jic,nic,nav,kl1,kl2
  end type src_type

  type :: pt_type
    sequence
    integer(2)	:: dat
    integer(2)	:: num
    integer(1)	:: min
    integer(1)	:: max
    integer(1)	:: src
  end type pt_type

  type :: fn_ctl_type
    sequence
    integer(1)		:: src
    integer(4)		:: nx0,ny0
    integer(4)		:: ncol,nrow
    integer(4)		:: yy_start,yy_end
    character(250) 	:: name
  end type fn_ctl_type

end module type_def


!=========================================================================
module data_def
!=========================================================================

  use type_def

  type(src_type), public		:: src
  type(fn_ctl_type), public		:: fn_ctl(20)

  integer(4), public			:: fn_ctl_num
  integer(4), public			:: key_params
  integer(4), public			:: key_tsv,key_kl1
  integer(4), public			:: unit_dat, unit_ctl, unit_params, unit_src

  type(pt_type), allocatable, public	:: pt(:,:,:,:)

  integer(4), public			:: yy_start, yy_end, mm_start, mm_end 

  integer(2), public			:: nodata

  character(250), public 		:: fn_ctl_list, fn_params, fn_dat

  integer(4), public			:: src_counter(0:200,1901:2020)

end module data_def


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

  use data_def

  character(250) 	:: arg

  nodata = 122; 

  fn_dat = 'ct_blended'; 

  mm_start = 01; yy_start = 1901; 
  mm_end = 12; yy_end = 2013; 

  fn_ctl_list = 'ctl_ct.txt'; fn_ctl_num = 0; 
  fn_params='params.txt'; key_params = 0; key_tsv = 0; key_kl1 = 0;

  unit_dat = 1; unit_ctl = 2; unit_params = 3; unit_src = 4;

  src%aar = 40; src%cis = 50; src%nic = 60; src%mix = 5;    
  src%dmi = 10; src%jma = 20; src%bsi = 30; src%jic = 70; 
  src%nav = 80; src%kl1 = 90; src%kl2 = 100;

  i = 0;

  do while (.true.)
    i = i + 1; if (i.gt.iargc()) exit;
    call getarg(i, arg)
! -i
    if (arg(1:2).eq.'-i') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, fn_ctl_list)
! -o
    elseif (arg(1:2).eq.'-o') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, fn_dat)
! -y
    elseif (arg(1:2).eq.'-y') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, arg)
      read (arg,*) mm_start,yy_start,mm_end,yy_end
! -p
    elseif (arg(1:2).eq.'-p') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, arg)
      read (arg,*) key_params
! -tsv
    elseif (arg(1:4).eq.'-tsv') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, arg)
      read (arg,*) key_tsv
! -kelly
    elseif (arg(1:6).eq.'-kelly') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, arg)
      read (arg,*) key_kl1
! -9
    elseif (arg(1:2).eq.'-9') then
      i = i + 1; if (i.gt.iargc()) exit;
      call getarg(i, arg)
      read (arg,*) nodata
    end if 
  end do

  if (key_params.eq.1) then 
    open (unit_params,file=fn_params,action='write',status='unknown')
    write (unit_params,'(a,i5)')  ' nodata_val:',nodata
    write (unit_params,'(a,4i5)') ' mm_start, yy_start, mm_end, yy_end:', mm_start,yy_start,mm_end,yy_end
  end if

end subroutine init_params 

!=========================================================================
subroutine read_ctl_list()
!=========================================================================

  use data_def

  character(250) 	:: arg

  open (unit_ctl,file=fn_ctl_list,action='read',status='old')

  do while (.true.)

    read (unit_ctl,'(a)',iostat=ios) arg; if (ios.ne.0) exit;
    if (len_trim(arg).eq.0.or.arg(1:1).eq.';') cycle
    fn_ctl_num = fn_ctl_num + 1;
    read(arg,*)                		&
      fn_ctl(fn_ctl_num)%name, 		& 
      fn_ctl(fn_ctl_num)%nx0,  		& 
      fn_ctl(fn_ctl_num)%ny0,  		& 
      fn_ctl(fn_ctl_num)%ncol, 		& 
      fn_ctl(fn_ctl_num)%nrow, 		&
      fn_ctl(fn_ctl_num)%src,  		&
      fn_ctl(fn_ctl_num)%yy_start, 	& 
      fn_ctl(fn_ctl_num)%yy_end
!    write (*,*) fn(fn_num)%name, fn(fn_num)%ncol, fn(fn_num)%nrow, fn(fn_num)%src, fn(fn_num)%yy1, fn(fn_num)%yy2
  end do 
  close (unit_ctl)
  if (key_params.eq.1) then 
    do i = 1,fn_ctl_num 
      l = len_trim(fn_ctl(i)%name)
      write (unit_params,'(a,7i5)') 	& 
      fn_ctl(i)%name(1:l), 		& 
      fn_ctl(i)%nx0,  			& 
      fn_ctl(i)%ny0,  			& 
      fn_ctl(i)%ncol, 			& 
      fn_ctl(i)%nrow, 			&
      fn_ctl(i)%src,  			&
      fn_ctl(i)%yy_start, 		& 
      fn_ctl(i)%yy_end
    end do
  end if
end subroutine read_ctl_list 


!=========================================================================
program grd2bld
!=========================================================================

  use data_def

  call init_params()
  call read_ctl_list()

  allocate (pt(1440,180,12,yy_start:yy_end))

  pt%dat = 0; pt%src = nodata; pt%num = 0; pt%min = 110; pt%max = 0;

  src_counter = 0

  do i = 1, fn_ctl_num
    call read_grd(fn_ctl(i))
  end do

  call read_bld('gdsidb_blended.dat','gdsidb_blended.src')
  call write_src_stat()

  where (pt%num.eq.0)
    pt%dat = nodata; pt%src = nodata; pt%min = nodata; pt%max = nodata;
  end where

  where (pt%num.gt.0)
    pt%dat = pt%dat/pt%num;
  end where
  
  l = len_trim(fn_dat)
  call write_grd(fn_dat(1:l)//'_mea.dat',0)
  call write_grd(fn_dat(1:l)//'_num.dat',1)
  call write_grd(fn_dat(1:l)//'_max.dat',2)
  call write_grd(fn_dat(1:l)//'_min.dat',3)
  call write_grd(fn_dat(1:l)//'_src.dat',4)

  if (key_tsv.eq.1) then
    call write_tsv(fn_dat(1:l)//'_mea.tsv',0)
    call write_tsv(fn_dat(1:l)//'_num.tsv',1)
    call write_tsv(fn_dat(1:l)//'_max.tsv',2)
    call write_tsv(fn_dat(1:l)//'_min.tsv',3)
    call write_tsv(fn_dat(1:l)//'_src.tsv',4)
  end if
  
end program grd2bld

!=========================================================================
subroutine write_src_stat()
!=========================================================================

  use data_def
  character(9)	:: yy_tag

  if (key_params.eq.1) then 
    write (unit_params,'(a)')  ' Summary statistics by data source:'
    write (unit_params,'(a9,a1,9(a9,a1))')   & 
      '   Source',';', &
      '    AARI',';',   &
      '     CIS',';',   &
      '     NIC',';',   &
      '    BSIM',';',   &
      '     DMI',';',   &
      '     JMA',';',   &
      '     JIC',';',   &
      '    NAVO',';',   &
      '   Kelly',';' 
    write (unit_params,'(a9,a1,9(i8,a1))') & 
      '     Code',';',src%aar,';',src%cis,';',src%nic,';',src%bsi,';',src%dmi,';',src%jma,';',src%jic,';',src%nav,';',src%kl1,';'
    write (unit_params,'(i4,a1,i4,a1,9(i9,a1))')	&  
      yy_start,'-',yy_end,';', 		&
      sum(src_counter(src%aar,:)),';', 	&
      sum(src_counter(src%cis,:)),';', 	&
      sum(src_counter(src%nic,:)),';', 	&
      sum(src_counter(src%bsi,:)),';', 	&
      sum(src_counter(src%dmi,:)),';', 	&
      sum(src_counter(src%jma,:)),';', 	&
      sum(src_counter(src%jic,:)),';', 	&
      sum(src_counter(src%nav,:)),';', 	&
      sum(src_counter(src%kl1,:)),';' 
    do i = 1901,2011,10
      write (unit_params,'(i4,a1,i4,a1,9(i9,a1))')  &
      i,'-',i+9,';', 				&
      sum(src_counter(src%aar,i:i+9)),';', 	&
      sum(src_counter(src%cis,i:i+9)),';', 	&
      sum(src_counter(src%nic,i:i+9)),';', 	&
      sum(src_counter(src%bsi,i:i+9)),';', 	&
      sum(src_counter(src%dmi,i:i+9)),';', 	&
      sum(src_counter(src%jma,i:i+9)),';', 	&
      sum(src_counter(src%jic,i:i+9)),';', 	&
      sum(src_counter(src%nav,i:i+9)),';', 	&
      sum(src_counter(src%kl1,i:i+9)),';' 
    end do
    do i = 1901,2020
      write (unit_params,'(i9,a1,9(i9,a1))')  &
      i,';', 				&
      src_counter(src%aar,i),';', 	&
      src_counter(src%cis,i),';', 	&
      src_counter(src%nic,i),';', 	&
      src_counter(src%bsi,i),';', 	&
      src_counter(src%dmi,i),';', 	&
      src_counter(src%jma,i),';', 	&
      src_counter(src%jic,i),';', 	&
      src_counter(src%nav,i),';', 	&
      src_counter(src%kl1,i),';' 
    end do
  end if
  return 

end subroutine write_src_stat

!=========================================================================
subroutine write_grd(fn,key)
!=========================================================================

  use data_def

  integer(1) 		:: buf(1440,180)
  character(*)	 	:: fn

  integer(2) 		:: dummy2
  character(10)		:: date_tag

  dummy2 = Z'0a0d' 
  open (unit_dat,file=fn,action='write',FORM='unformatted',access='stream',STATUS='unknown')
  do iyy = yy_start, yy_end
    do imm = 1,12
      if (iyy.eq.yy_start.and.imm.lt.mm_start) cycle
      if (iyy.eq.yy_end.and.imm.gt.mm_end) exit
      write (date_tag,'(i2.2,a1,i2.2,a1,i4.4)') 15,'.',imm,'.',iyy
      select case (key)
      case (0)
        buf = pt(:,:,imm,iyy)%dat
      case (1)
        buf = pt(:,:,imm,iyy)%num
      case (2)
        buf = pt(:,:,imm,iyy)%max
      case (3)
        buf = pt(:,:,imm,iyy)%min
      case (4)
        buf = pt(:,:,imm,iyy)%src
      end select
      write (1) date_tag,dummy2
      do j = 1,180
        write (1) buf(:,j),dummy2
      end do
    end do
  end do

  return

end subroutine write_grd

!=========================================================================
subroutine read_grd(fn)
!=========================================================================

  use data_def

  type (fn_ctl_type) fn

  integer(2) 			:: dat0(1440,180), dat1(1440,180), dnum(1440,180)
  integer(1)			:: dmin(1440,180), dmax(1440,180), src1(1440,180)
  integer(1)			:: src0, src_mix
  integer(2)			:: dummy2
  character(10)			:: date_tag
  integer(1)			:: bs(100)

  integer(1), allocatable 	:: buf(:,:)
  integer(4)			:: counter
  bs = 8

  allocate (buf(fn%ncol+2,fn%nrow))

  open (unit_ctl,file=fn%name,action='read',form='unformatted',access='stream',status='old')
  l = len_trim(fn%name);
  counter = 0;
  src0 = fn%src; src_mix = src%mix;

  i1_1 = fn%nx0+1; i1_2 = fn%nx0+fn%ncol; l1_1 = 1; l1_2 = fn%ncol;
  j1 = fn%ny0+1; j2 = fn%ny0+fn%nrow-1;
  if (i1_2.le.1440) then
    key_split = 0;
  else
    key_split = 1; 
    l1_1 = 1;        l1_2 = fn%ncol - (i1_2 - 1440); i1_1 = i1_1; i1_2 = 1440; 
    l2_1 = l1_2 + 1; l2_2 = fn%ncol;                 i2_1 = 1;    i2_2 = (l2_2-l2_1)+1;
  end if

  if (j2.le.180) then
    ly = fn%nrow
  else
    ly = fn%nrow - (j2 - 180); j2 = 180; 
  end if

  write (*,*)

  do while (.true.)
    read (unit_ctl,iostat=ios) date_tag, dummy2, buf; if (ios.ne.0) exit;
    read (date_tag,'(i2,1x,i2,1x,i4)') idd,imm,iyy;

    if (iyy.lt.yy_start.or.iyy.gt.yy_end) cycle 
    if (iyy.lt.fn%yy_start.or.iyy.gt.fn%yy_end) cycle 

    write (*,'(a,1x,i4.4,a1,i2.2,a1,i2.2,100a)',advance="no") fn%name(1:l),iyy,'-',imm,'-',idd, bs(1:l+11);
    counter = counter + 1

    dat0 = nodata;
    if (key_split.eq.0) then
      dat0(i1_1:i1_2,j1:j2) = buf(l1_1:l1_2,1:ly); 
    else
      dat0(i1_1:i1_2,j1:j2) = buf(l1_1:l1_2,1:ly); 
      dat0(i2_1:i2_2,j1:j2) = buf(l2_1:l2_2,1:ly); 
    end if
    dat1 = pt(:,:,imm,iyy)%dat; dnum = pt(:,:,imm,iyy)%num;
    dmin = pt(:,:,imm,iyy)%min; dmax = pt(:,:,imm,iyy)%max;
    src1 = pt(:,:,imm,iyy)%src; 

    where (dat0.lt.0.or.dat0.gt.110)
      dat0 = nodata
    end where

    where (dat0.gt.100.and.dat0.le.110)
      dat0 = 101
    end where

    do j = 1,180
      do i = 1,1440
        if (dat0(i,j).ne.nodata) then
          src_counter(src0,iyy) = src_counter(src0,iyy) + 1
          if (dnum(i,j).eq.0) then
            src1(i,j) = src0
          elseif (src1(i,j).ne.src0) then
            src1(i,j) = src_mix
          end if
          dnum(i,j) = dnum(i,j) + 1
          dat1(i,j) = dat1(i,j) + dat0(i,j)
          if (dmin(i,j).gt.dat0(i,j)) dmin(i,j) = dat0(i,j) 
          if (dmax(i,j).lt.dat0(i,j)) dmax(i,j) = dat0(i,j) 
        end if
      end do
    end do

    pt(:,:,imm,iyy)%dat = dat1; pt(:,:,imm,iyy)%num = dnum; 
    pt(:,:,imm,iyy)%min = dmin; pt(:,:,imm,iyy)%max = dmax; 
    pt(:,:,imm,iyy)%src = src1;

  end do
  if (key_params.eq.1) then 
    write (unit_params,'(1x,a,a24,a,i6)')  ' Grids from ',fn%name(1:l),' :',counter
  end if

  close (unit_ctl);
  deallocate (buf)

  return
 
end subroutine read_grd

!====================================================
subroutine read_bld(fn1,fn2)
!====================================================

  use data_def

  character(*)			:: fn1, fn2

  integer(2) 			:: dat0(1440,180)
  integer(1)			:: src0(1440,180) 
  integer(2)			:: dat1(1440,180), dnum(1440,180)
  integer(1)			:: dmin(1440,180), dmax(1440,180), src1(1440,180)
  integer(1)			:: src_bld,src_mix
  integer(2)			:: dummy2
  character(10)			:: date_tag
  integer(1)			:: bs(100)

  integer(1), allocatable 	:: buf1(:,:), buf2(:,:)
  integer(4)			:: counter(2)
  bs = 8

  allocate (buf1(1440+2,181), buf2(1440+2,181))

  open (unit_ctl,file=fn1,action='read',FORM='unformatted',access='stream',STATUS='old')
  open (unit_src,file=fn2,action='read',FORM='unformatted',access='stream',STATUS='old')

  l = len_trim(fn1)
  write (*,*)

  do while (.true.)
    read (unit_ctl,iostat=ios) date_tag, dummy2, buf1; if (ios.ne.0) exit;
    read (date_tag,'(i2,1x,i2,1x,i4)') idd,imm,iyy; 
    read (unit_src,iostat=ios) date_tag, dummy2, buf2; if (ios.ne.0) exit;

    if (iyy.ge.1972) cycle

    write (*,'(a,1x,i4.4,a1,i2.2,a1,i2.2,100a)',advance="no") fn1(1:l),iyy,'-',imm,'-',idd, bs(1:l+11);

    dat0 = buf1(1:1440,1:180); src0 = buf2(1:1440,1:180); 

    dat1 = pt(:,:,imm,iyy)%dat; dnum = pt(:,:,imm,iyy)%num;
    dmin = pt(:,:,imm,iyy)%min; dmax = pt(:,:,imm,iyy)%max;
    src1 = pt(:,:,imm,iyy)%src; 

    where (dat0.lt.0.or.dat0.gt.110)
      dat0 = nodata
    end where

    where (dat0.gt.100.and.dat0.le.110)
      dat0 = 101
    end where

!  src%aar = 40; src%cis = 50;  src%nic = 60;  src%mix = 5;    
!  src%dmi = 10; src%jma = 20;  src%bsi = 30; src%jic = 70; 
!  src%nav = 80; src%kl1 = 90; src%kl2 = 100;

    where (buf2(1:1440,1:180).eq.10)
      src0 = src%dmi
    end where
    where (buf2(1:1440,1:180).eq.20)
      src0 = src%jma
    end where
    where (buf2(1:1440,1:180).eq.100)
      src0 = src%bsi
    end where
    where (buf2(1:1440,1:180).eq.120)
      src0 = src%aar
    end where
    where (buf2(1:1440,1:180).eq.110)
      src0 = src%cis
    end where
    where (buf2(1:1440,1:180).eq.90)
      src0 = src%nic
    end where
    where (buf2(1:1440,1:180).eq.50)
      src0 = src%jic
    end where
    where (buf2(1:1440,:).eq.30)
      src0 = src%nav
    end where
    where (buf2(1:1440,:).eq.40)
      src0 = src%kl1
    end where
    where (buf2(1:1440,:).eq.70)
      src0 = src%kl2 
    end where

    if (key_kl1.eq.1) then
      where ((buf2(1:1440,:).eq.60).or.(buf2(1:1440,:).eq.70).or.(buf2(1:1440,:).eq.80))
        src0 = nodata ! climatology and SSMR-SSM/I, temporal extension of Kelly grids 
      end where
    else
      where ((buf2(1:1440,:).eq.60).or.(buf2(1:1440,:).eq.70).or.(buf2(1:1440,:).eq.80).or.(buf2(1:1440,:).eq.40))
        src0 = nodata ! climatology and SSMR-SSM/I, temporal extension of/and Kelly grids 
      end where
    end if
!NEW  / OLD source:
! 10  /  10     - Danish Meteorological Institute
! 20  /  20     - Japan Meteorological Agency
! 30  / 100     - BSIM (1960 - 1979)
! 40  / 120     - AARI based on SIGRID (1933 - 1992)
! 41  /   -     - AARI based on SIGRID3 (1999 - 2012...)
! 50  / 110     - CIS based on SIGRID (1968/71 - 1998)
! 51  /   -     - CIS based on SIGRID3 (2006 - 2012...)
! 60  /  90     - NIC based on EASE-GRID (1972 - 2007) 
! 61  /  -      - NIC based on SIGRID3 (2003 - 2012...)
! 70  /  50     - Walsh and Johnson/Navy-NOAA Joint Ice Center
! 80  /  30     - Naval Oceanographic Office (NAVOCEANO)
! 90  /  40     - Kelly ice extent grids
! nodata / 60   - climatological values from basically AARI/BSIM/CIS/NIC observations 
! nodata / 70   - temporal extension of Kelly grids (see note below)
! nodata / 80   - Nimbus-7 SSMR Arctic Sea Ice Concentrations or DMSP SSM/I Sea Ice

    where (src0.eq.nodata)
      dat0 = nodata
    end where

    do j = 1,180
      do i = 1,1440
        if (dnum(i,j).gt.0) cycle
        if (dat0(i,j).ne.nodata) then
          src_bld = src0(i,j)
          src_counter(src_bld,iyy) = src_counter(src_bld,iyy) + 1
          src1(i,j) = src_bld
          dnum(i,j) = dnum(i,j) + 1
          dat1(i,j) = dat1(i,j) + dat0(i,j)
          if (dmin(i,j).gt.dat0(i,j)) dmin(i,j) = dat0(i,j) 
          if (dmax(i,j).lt.dat0(i,j)) dmax(i,j) = dat0(i,j) 
        end if
      end do
    end do

    pt(:,:,imm,iyy)%dat = dat1; pt(:,:,imm,iyy)%num = dnum; 
    pt(:,:,imm,iyy)%min = dmin; pt(:,:,imm,iyy)%max = dmax; 
    pt(:,:,imm,iyy)%src = src1;

  end do 
  
  close (unit_ctl);
  close (unit_src);

  return

end subroutine read_bld

!===========================
subroutine write_tsv(fn,key)
!===========================

  use data_def

  character(*)	 	:: fn
  character(117)	:: tag

  integer(4) 		:: x(1440,180)
  integer(2)		:: dummy2
  integer(1)		:: bs(100)
  character(7)          :: ch_lon(1440)
  character(7)     	:: ch_lat(180)
  character(5)		:: ch_x(0:122)
  character(10087)	:: buf_lon
  character(7207)	:: buf_x
  

  dummy2 = Z'0a0d'; bs = 8;

  l = len_trim(fn)

  tag = 'cpt:field=ct_ice, cpt:T=1901-01-15, cpt:nrow=1440, cpt:ncol=180, cpt:row=Y, cpt:col=X, cpt:units=%, cpt:missing=122.0'                

  xlat = 90. + 0.25
  do i = 1, 180
    xlat = xlat - 0.25; write (ch_lat(i),'(f7.2)') xlat;
  end do
  xlon = -0.25; buf_lon(1:7) = ' '; ix = 8; 
  do i = 1, 1440
    xlon = xlon + 0.25; write (buf_lon(ix:ix+6),'(f7.2)') xlon; ix = ix + 7;
  end do

  do i = 0, 101
    write (ch_x(i),'(f5.0)') float(i);
  end do
  ch_x(122) = ' 122.'

  open (unit_dat,file=fn,action='write',FORM='unformatted',access='stream',STATUS='unknown')

  write (unit_dat) 'xmlns:cpt=http://iri.columbia.edu/CPT/v10/',dummy2,'cpt:nfields=1',dummy2

  idd = 15

  do iyy = yy_start, yy_end
    do imm = 1,12

      if (iyy.eq.yy_start.and.imm.lt.mm_start) cycle
      if (iyy.eq.yy_end.and.imm.gt.mm_end) exit

      select case (key)
      case (0)
        x = pt(:,:,imm,iyy)%dat
      case (1)
        x = pt(:,:,imm,iyy)%num
      case (2)
        x = pt(:,:,imm,iyy)%max
      case (3)
        x = pt(:,:,imm,iyy)%min
      case (4)
        x = pt(:,:,imm,iyy)%src
      end select

      write (tag(25:34),'(i4.4,a1,i2.2,a1,i2.2)') iyy,'-',imm,'-',idd;
      write (unit_dat) tag,dummy2; 
      write (unit_dat) buf_lon,dummy2;
      do j = 1,180
        buf_x(1:7) = ch_lat(j); ix = 8;
        do i = 1,1440
          buf_x(ix:ix+4) = ch_x(x(i,j)); ix = ix + 5;
        end do
        write (unit_dat) buf_x,dummy2;
      end do

      write (*,'(a,1x,i4.4,a1,i2.2,a1,i2.2,100a1)',advance="no") fn(1:l),iyy,'-',imm,'-',idd, bs(1:l+11)

    end do
  end do
  close (unit_dat)

  return

end subroutine write_tsv