module util
type punto
integer:: i,j,a
end type punto

 end module util

program readv3

! This utility program is written in free-format Fortran 90.
!   It requires a Fortran 90 compiler to compile. On a DEC_Alpha
!   machine, type the following to compile:
!   
!   f90 -free -convert big_endian readv3.f
!

  use util

  implicit none
  integer, dimension(50,20) :: bhi
  real, dimension(20,20) :: bhr
  character(len=80), dimension(50,20) :: bhic
  character(len=80), dimension(20,20) :: bhrc
  character(len=120) :: flnm
  integer :: iunit = 99

  integer :: flag

  integer :: ndim
  real :: time, sample
  integer, dimension(4) :: start_index, end_index
  character (len= 4) :: staggering
  character (len= 4) :: ordering
  character (len=24) :: start_date
  character (len=24) :: current_date
  character (len= 9) :: name
  character (len=25) :: units
  character (len=46) :: description

  integer :: l

  real, allocatable, dimension(:,:,:,:) :: data, datam

  integer :: ierr, ier

  logical :: newtime = .TRUE.
  logical :: lmore

! ####### Definiciones ########

  type(punto),allocatable :: est(:)
  integer :: i1,i2,i3,i
  integer :: NumEst,NumVarLec,NumVarEsc,NumEnsembles,NumSaida
  integer :: R,G
  real    :: p00,ts0,tlp,ptop,sigmahsfc,phydro,z
  character(len=20) :: RetrChuviaChar
  logical :: imprime=.false.,imprime_cabeceira=.false.
  character(len=9) :: NameSrn,NameSrc,NameSrcrn

  character(len=9),allocatable :: NomeVarLec(:)
  character(len=50),allocatable :: NomeFichIn(:), NomeFichOu(:)
  real, allocatable :: VarEsc(:,:,:), VarLec(:,:,:)
  real, allocatable :: phydrosfc(:),tempest(:) 
  integer,allocatable :: NumSalvas(:)

  integer ::TmpInt,ContadorChuvia,RetrChuvia
  real, allocatable, dimension(:,:,:,:) :: PrecRc,PrecRn,Proba,PrecRcm,PrecRnm

  real, allocatable, dimension(:,:) :: VarTmp1,VarTmp2,VarTmp3,VarTmp4
  real, allocatable, dimension(:,:) :: ArrayTmp1,ArrayTmp2,ArrayTmp3


! #############################################
! ########  OPCIONES DEL PROGRAMA  ############
! #############################################

! ############ MEMBROS ENSEMBLE ###########

NumEnsembles = 10
allocate(NomeFichIn(NumEnsembles+1))

NomeFichIn(1)   = 'MMOUT_DOMAIN2-01_EX'
NomeFichIn(2)   = 'MMOUT_DOMAIN2-02_EX'
NomeFichIn(3)   = 'MMOUT_DOMAIN2-03_EX'
NomeFichIn(4)   = 'MMOUT_DOMAIN2-04_EX'
NomeFichIn(5)   = 'MMOUT_DOMAIN2-05_EX'
NomeFichIn(6)   = 'MMOUT_DOMAIN2-06_EX'
NomeFichIn(7)   = 'MMOUT_DOMAIN2-07_EX'
NomeFichIn(8)   = 'MMOUT_DOMAIN2-08_EX'
NomeFichIn(9)   = 'MMOUT_DOMAIN2-09_EX'
NomeFichIn(10)  = 'MMOUT_DOMAIN2-10_EX'


!############ FICHEIROS SAIDA ###########

NumSaida = 2

allocate(NomeFichOu(NumSaida))

NomeFichOu(1) = 'MMOUT_DOMAIN2-MS'
NomeFichOu(2) = 'MMOUT_DOMAIN2-S'

!############ ARRAY DE PRECIPITACION ######

RetrChuvia = 24  ! Indica el intervalo en intervalos temporales con los que 
                 ! se calcula el SPREAD de la precipitacion

!########################################

write(RetrChuviaChar,*) RetrChuvia

write(RetrChuviaChar,'(i2)') RetrChuvia

NameSrc = 'S'//trim(RetrChuviaChar)//'RC'     ! Nombres de las variables de SPREAD
NameSrn = 'S'//trim(RetrChuviaChar)//'RN'     ! de la precipitacion acumulada
NameSrcrn = 'S'//trim(RetrChuviaChar)//'RCRN' ! 

ContadorChuvia=0 ! Contador de dias de lluvia 

!########################################

  call arguments(flnm, lmore)

  do i=1,NumEnsembles
      iunit=i
      print*, 'flnm = ', trim(NomeFichIn(i))
      open(iunit, file=NomeFichIn(i), form='unformatted', status='old', action='read')
  enddo

  do i=1,NumSaida
      iunit=i+NumEnsembles
      open(iunit, file=NomeFichOu(i), form='unformatted', status='replace', action='write')
  enddo
  
  do i=1,NumEnsembles
  iunit=i
  read(iunit, iostat=ierr) flag
  enddo

  do i=1,NumSaida
      iunit=i+NumEnsembles
      if (ierr.eq.0) write(iunit) flag
  enddo

!################################
!###  BUCLE DE LECTURA ##########
!################################

  do while (ierr == 0)


!#######################
!###  FLAG == 0    #####
!#######################

     if (flag == 0) then

  do i=1,NumEnsembles
     iunit=i
     read(iunit,iostat=ier) bhi, bhr, bhic, bhrc
        if(ier/=0) then
           write(*,'("Error reading big header")')
           call abort()
        endif
     if(imprime_cabeceira)Then
        call printout_big_header(bhi, bhr, bhic, bhrc)
     endif
  enddo
  do i=1,NumSaida
        iunit=i+NumEnsembles
        write(iunit) bhi, bhr, bhic, bhrc
  enddo

!######## ARRAYS DE PRECIPITACION ######

  allocate(PrecRc(BHI(16,1),BHI(17,1),0:RetrChuvia,NumEnsembles))
  allocate(PrecRn(BHI(16,1),BHI(17,1),0:RetrChuvia,NumEnsembles))
  allocate(PrecRnm(BHI(16,1),BHI(17,1),0:RetrChuvia,1))
  allocate(PrecRcm(BHI(16,1),BHI(17,1),0:RetrChuvia,1))

  PrecRc = 0
  PrecRn = 0
  PrecRcm = 0
  PrecRnm = 0

!#####################################


!#######################
!###  FLAG == 1    #####
!#######################

     elseif (flag == 1) then

do i=1,NumEnsembles
        iunit=i
        READ (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering,&
             current_date, name, units, description

!#### IMPRIME A SUB CABECEIRA A MODO DE ESTUDIO - NADA MAIS- #####
!        print*, '############   ',iunit,'   ############'
!        sample=0
!        write(*,'(A8,1x,I1,4(1x,I3),1x,A,1x,A," : ", F20.8,1x,A)')&
!             name, ndim, end_index(1), end_index(2), end_index(3), end_index(4),&
!             staggering, ordering, sample, trim(units)
!#################################################################

        if(ier/=0) then
           write(*,'("Error reading subheader")')
           call abort()
        endif
!        print*, 'Lei la subcabecera de ', iunit 
enddo
    do i=1,NumSaida
            iunit=i+NumEnsembles
            WRITE (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering, &
                   current_date, name, units, description
    enddo

        if (lmore) then
           print*, 'ndim: ', ndim
           print*, 'start_index: ', start_index
           print*, 'end_index: ', end_index
           print*, 'time: ', time
           print*, 'staggering: #'//staggering//'#'
           print*, 'ordering: #'//ordering//'#'
           print*, 'date/time: #'//current_date//'#'
           print*, 'name: #'//name//'#'
           print*, 'units: #'//units//'#'
           print*, 'description: #'//description//'#'
        endif

        if (newtime) then
           print*, 'Estou lendo.....'
           write(*,'(/,A,2x, F15.5," Hours"/)') current_date, time/60.
           newtime = .FALSE.
        endif
!print*, 'cheguei 2'
        if (ndim == 1) then
           allocate(datam(end_index(1), 1, 1, NumSaida))
           allocate(data(end_index(1), 1, 1, NumEnsembles))
        elseif (ndim == 2) then
           allocate(datam(end_index(1), end_index(2), 1, NumSaida))
           allocate(data(end_index(1), end_index(2), 1, NumEnsembles))
        elseif (ndim == 3) then
           allocate(datam(end_index(1), end_index(2), end_index(3), NumSaida))
           allocate(data(end_index(1), end_index(2), end_index(3), NumEnsembles))
        endif
!print*, 'cheguei 3'
        datam=0
        data=0
!print*, 'cheguei 4'


do i=1,NumEnsembles
!print*, 'cheguei 5'

        iunit=i

        read(iunit) data(:,:,:,iunit)

   if(imprime) then
        if (ndim == 3) then
            sample = data( end_index(1)/2,end_index(2)/2,end_index(3)/2,1 )
        else if (ndim == 2) then
            sample = data( end_index(1)/2,end_index(2)/2,1,1)
        else if (ndim == 1) then
            sample = data( end_index(1)/2,1,1,1)
        end if

        write(*,'(A8,1x,I1,4(1x,I3),1x,A,1x,A," : ", F20.8,1x,A)')&
             name, ndim, end_index(1), end_index(2), end_index(3), end_index(4),&
             staggering, ordering, sample, trim(units)
   endif


enddo


!print*, 'cheguei 6'

! ### CALCULO DE LA MEDIA ##
    do  i=1,NumEnsembles
        datam(:,:,:,1) = datam(:,:,:,1) + data(:,:,:,i) 
    enddo
    datam(:,:,:,1) = datam(:,:,:,1)/NumEnsembles
! #########################


! ### CALCULO DEL SPREAD ##
    do  i=1,NumEnsembles
        datam(:,:,:,2) = datam(:,:,:,2) + (datam(:,:,:,1) - data(:,:,:,i))**2
    enddo
    datam(:,:,:,2) = sqrt( datam(:,:,:,2)/NumEnsembles )
! ##########################

     do  i=1,NumSaida
         iunit=i+NumEnsembles
         write(iunit) datam(:,:,:,i)
     enddo

!#### SALVA DE LA PRECIPITACION ###

if(name.eq.'RAIN CON ') then
  ContadorChuvia = ContadorChuvia + 1 
  PrecRc(:,:,0,:) = data(:,:,1,:)
  PrecRcm(:,:,0,1) = datam(:,:,1,1)

elseif(name.eq.'RAIN NON ') then
  ContadorChuvia = ContadorChuvia + 1
  PrecRn(:,:,0,:) = data(:,:,1,:)
  PrecRnm(:,:,0,1) = datam(:,:,1,1)

endif

!######################################################

! ### ESCRITURA DEL SPREAD EN EL FICHERO DE LA MEDIA ###

iunit = 1 + NumEnsembles
TmpInt = 1
i = 2

name = 'S'//name
 write(iunit) TmpInt 

 WRITE (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering, &
                   current_date, name, units, description


 write(iunit) datam(:,:,:,i)


!################################################

if(ContadorChuvia==2) then

!#####################################################
!##  ESTA BANDERA INDICA QUE SE HAN RECOGIDO LOS DATOS 
!##  DE LA RAINCON Y RAINNON. APROVECHO PARA CALCULAR
!##  Y ESCRIBIR EL SPREAD DE LAS VARIABLES
!##  DERIVADAS.
!########################################
!## SPREAD DE LA PREC ACUMULADA
!########################################

ContadorChuvia = 0

allocate(VarTmp1(end_index(1),end_index(2)))
allocate(VarTmp2(end_index(1),end_index(2)))
allocate(VarTmp3(end_index(1),end_index(2)))
allocate(VarTmp4(end_index(1),end_index(2)))
allocate(ArrayTmp1(end_index(1),end_index(2)))
allocate(ArrayTmp2(end_index(1),end_index(2)))
allocate(ArrayTmp3(end_index(1),end_index(2)))

VarTmp1 = 0
VarTmp2 = 0
VarTmp3 = 0
VarTmp4 = 0
ArrayTmp1 = 0
ArrayTmp2 = 0
ArrayTmp3 = 0

VarTmp3=PrecRcm(:,:,0,1)-PrecRcm(:,:,RetrChuvia,1)
VarTmp4=PrecRnm(:,:,0,1)-PrecRnm(:,:,RetrChuvia,1)
do i=1,NumEnsembles

  VarTmp1 = PrecRc(:,:,0,i)-PrecRc(:,:,RetrChuvia,i)
  VarTmp2 = PrecRn(:,:,0,i)-PrecRn(:,:,RetrChuvia,i)

  ArrayTmp1 = ArrayTmp1 + (VarTmp3-VarTmp1)**2
  ArrayTmp2 = ArrayTmp2 + (VarTmp4-VarTmp2)**2
  ArrayTmp3 = ArrayTmp3 + (VarTmp3 + VarTmp4 - VarTmp1 - VarTmp2)**2

!print*, PrecRc(25,25,0,i)+PrecRn(25,25,0,i),PrecRc(25,25,RetrChuvia,i)+PrecRn(25,25,RetrChuvia,i)

enddo

ArrayTmp1 = sqrt(ArrayTmp1/NumEnsembles)
ArrayTmp2 = sqrt(ArrayTmp2/NumEnsembles)
ArrayTmp3 = sqrt(ArrayTmp3/NumEnsembles)

!##########################################
!###### SALVA  EN EL FICHERO ##############
!##########################################


TmpInt = 1
iunit=NumEnsembles+1


! ### SRC ###

 write(iunit) TmpInt

 WRITE (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering, &
                   current_date, NameSrc, units, description

 write(iunit) ArrayTmp1

! ### SRN ###

name='S24RN'
 write(iunit) TmpInt

 WRITE (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering, &
                   current_date, NameSrn, units, description


 write(iunit) ArrayTmp2

! ### SRCRN ###

 write(iunit) TmpInt

 WRITE (iunit,iostat=ier) ndim, start_index, end_index, time, staggering, ordering, &
                   current_date, NameSrcrn, units, description

 write(iunit) ArrayTmp3

!### ACTUALIZO EL ARRAY DE LA PRECIPITACION  ###

do i=RetrChuvia,1,-1
  precRc(:,:,i,:) = precRc(:,:,i-1,:)
  precRcm(:,:,i,:) = precRcm(:,:,i-1,:)
  precRn(:,:,i,:) = precRn(:,:,i-1,:)
  precRnm(:,:,i,:) = precRnm(:,:,i-1,:)
enddo

deallocate(VarTmp1)
deallocate(VarTmp2)
deallocate(VarTmp3)
deallocate(VarTmp4)
deallocate(ArrayTmp1)
deallocate(ArrayTmp2)
deallocate(ArrayTmp3)

endif

!###################################################


        deallocate(data)
        deallocate(datam)

!####################################
!#######  FLAG == 2    ##############
!####################################


     elseif (flag == 2) then
        newtime = .TRUE.
     else
        stop
     endif
do i=1,NumEnsembles
     iunit=i
     read(iunit, iostat=ierr) flag
enddo

do  i=1,NumSaida
     iunit=i+NumEnsembles
     if (ierr.eq.0) write(iunit) flag
enddo
  enddo

  write(*,'(/,"Hit the end of file of unit ", I3)') iunit


end program readv3

   subroutine printout_big_header(bhi, bhr, bhic, bhrc)

     implicit none
     integer, dimension(50,20) :: bhi
     real, dimension(20,20) :: bhr
     character(len=80), dimension(50,20) :: bhic
     character(len=80), dimension(20,20) :: bhrc
     integer :: i, j, v3j
   
     write(*,'(/)')
     v3j = bhi(1,1)
     if (bhi(1,1) == 11) v3j = v3j+5
     do j = 1, v3j
      if (j < 8 .or. j>10) then
        if (j == 1) write(*, '("TERRAIN Portion of big header:")')
        if (j == 2) write(*, '(/,"REGRID Portion of big header:")')
        if (j == 3) write(*, '(/,"RAWINS Portion of big header:")')
        if (j == 4) write(*, '(/,"SFC RAWINS Portion of big header:")')
        if (j == 5) write(*, '(/,"INTERP Portion of big header:")')
        if (j == 11) write(*, '(/,"MM5 Portion of big header:")')
        if (j == 6) write(*, '(/,"MM5 Substrate Temp File big header:")')
        if (j == 7) write(*, '(/,"MM5 Boundary File big header:")')
        if (j == 8) write(*, '(/,"Interpolated MM5 Portion of big header:")')
        write(*,'(/,"***Integers:"/)')
        do i = 1, size(bhi,1)
           if (bhi(i,j) /= -999) then
              write(*,'("BHI(",I3,",",I2,"):",I8," : ",A)')&
                   i, j, bhi(i,j),trim(bhic(i,j))
           endif
        enddo
   
        write(*,'(/,"***Floats:"/)')
        do i = 1, size(bhr,1)
           if (bhr(i,j) /= -999.) then
              write(*,'("BHR(",I3,",",I2,"):",F9.2," : ",A)')&
                   i, j, bhr(i,j),trim(bhrc(i,j))
           endif
        enddo
        write(*,'(/)')
      endif
     enddo
   end subroutine printout_big_header
   
   subroutine arguments(v2file, lmore)
     implicit none
     character(len=*) :: v2file
     character(len=120) :: harg
     logical :: lmore
   
     integer :: ierr, i, numarg
     integer, external :: iargc
   
     numarg = iargc()
   
     i = 1
     lmore = .false.
   
     do while ( i < numarg) 
        call getarg(i, harg)
        print*, 'harg = ', trim(harg)
   
        if (harg == "-v") then
           i = i + 1
           lmore = .true.
        elseif (harg == "-h") then
           call help
        endif
   
     enddo
   
     call getarg(i,harg)
     v2file = harg
   end subroutine arguments
   
   subroutine help
     implicit none
     character(len=120) :: cmd
     call getarg(0, cmd)
   
     write(*,'(/,"Usage: ", A, " [-v] v2file ")') trim(cmd)
     write(*,'(8x, "-v     : Print extra info")')
     write(*,'(8x, "v3file : MM5v3 file name to read.")')
     write(*,'(8x, "-h     : print this help message and exit.",/)')
     stop
   end subroutine help