C23456
      program crear_NetCDF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C KOBI 04/2006:
C ============
C Este programa comvertira la informacion de una serie de
C tiempo en TXT al formato NetCDF. 
C 
C       crucen los dedos!!!1, :)
C
C Se compila de la siguiente manera :
C
C   f90 crear_NetCDF.f  /usr/local/netcdf/lib/libnetcdf.a
C
C KOBI 04/2013:
C ============
C Se hicieron modificaciones en el formato de lectura para
C que pueda ser compilado con gfortran
C
C KOBI 01/2014:
C ============
C Para compilar con gfortran la version 4 de Netcdf
C   gfortran crear_NetCDF.f -L/usr/lib -lnetcdf -lnetcdff
C 
C KOBI 14/01/2014:
C ===============
C    ifort crear_NetCDF.f /usr/local/netcdf-3.6.3-ifort10/lib64/libnetcdf.a
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      real data(10000,4)
      integer dimid,vid(4),iyi,iyf,imi,imf,idi,idf,ij,nreg,ndias
      character fecha*10,hora*8, estacion*20,sss*10,unidades*4
      character titulo*37, descripcion*35,variable*9
      real data_time1(10000),lat,lon,alt

      include '/usr/local/netcdf/include/netcdf.inc'

C KOBI 14/01/2014:
C =============== 
C     include '/usr/local/netcdf-3.6.3-ifort10/include/netcdf.inc'  

C ===========================================================
C Aqui defino algunos parametros como
C
C    undef= -999.0  lo usare cuando el valor no se conozca
C    lat  =  -5.    Latitud de la estacion 
C    lon  = -78.    Longitud de la estacion
C    alt  =   0.    Altura de la estacion 0. msnm
C    nreg =  72     Numero de datos en esta serie
C
C    variable="temp_aire"
C           Definimos el nombre la variable Temperatura del 
C           Aire que indetificara, y usara el programa
C           grafico para plotearlo
C
C    titulo ="Lat=-5. Lon=-78. Temperatura del Aire"
C           Entrara como un atributo global del archivo
C           y sera leido, al menos, por GrADS y FERRET 
C
C    descripcion="Temperatura del Aire en La Calurosa"
C           Descripcion de la variable temp_aire, es 
C           un atributo de la misma.
C
C    unidades="degC"
C           Unidades de la variable tempratura que significa
C           Grados Celsius, pueden colocar tambien C, Celsius,
C           etc
C
C
C ===========================================================

      data undef/-999.0/,lat/-5./,lon/-78./,alt/0./,nreg/72/
      data variable/"temp_aire"/
      data titulo/"Lat=-5. Lon=-78. Temperatura del Aire"/
      data descripcion/"Temperatura del Aire en La Calurosa"/
      data unidades/"degC"/
  

C =====================================
C    Inicializamos la data a cero
C =====================================

      do i=1,4
      do j=1,10000
         data(j,i)=0.
      enddo
      enddo

C =========================================
C    Leemos la data del archivo data_ta.txt
C =========================================

      open(1,file="data_ta.txt",form="formatted",status="old")
      read(1,*) sss

      do i=1,nreg
         read(1,'(A10,1X,A8,7X,f10.1)') fecha,hora,data(i,1)
         read(fecha(1:2),'(I2)') id
         read(fecha(4:5),'(I2)') im
         read(fecha(7:10),'(I4)') iy
         read(hora(1:2),'(I2)') ih
         call time_en_horas(iy,im,id,ih,data_time1(i))
         write(*,*) i,iy,im,id,ih,int(data_time1(i))
      enddo
      close(1)

C =====================================
C     empezamos a crear el archivo 
C     netcdf al cual le pondremos el 
C     nombre de   data_ta.nc
C =====================================

      status=nf_create('data_ta.nc',0,ncid)

C =====================================
C     definimos las dimensiones
C     (time,height,lat,lon)
C =====================================

      status=nf_def_dim(ncid,'time',nreg,dimid)
      write(*,*) dimid

      status=nf_def_dim(ncid,'height',1,dimid)
      write(*,*) dimid

      status=nf_def_dim(ncid,'lat',1,dimid)
      write(*,*) dimid

      status=nf_def_dim(ncid,'lon',1,dimid)
      write(*,*) dimid

C =====================================
C     definimos variables
C     (time,height,lat,lon,temp_aire)
C =====================================

      vid(1)=1
      vid(2)=0
      vid(3)=0
      vid(4)=0
      status=nf_def_var(ncid,'time',NF_FLOAT,1,vid,dimid)
      status=nf_put_att_text(ncid,dimid,'units',
     & 36,"hours since 2000-1-1 00:00:0.0 -5:00")
      status=nf_put_att_text(ncid,dimid,'long_name',4,"Time")
      status=nf_put_att_text(ncid,dimid,'delta_t',
     & 19,"0000-00-00 01:00:00")
      status=nf_put_att_text(ncid,dimid,'axis',1,'T')
      write(*,*) "ID Time: ",dimid

      vid(1)=2
      vid(2)=0
      vid(3)=0
      vid(4)=0
      status=nf_def_var(ncid,'height',NF_FLOAT,1,vid,dimid)
      status=nf_put_att_text(ncid,dimid,'units',2,"m ")
      status=nf_put_att_text(ncid,dimid,'axis',1,'Z')
      write(*,*) "ID DEPTH: ",dimid


      vid(1)=3
      vid(2)=0
      vid(3)=0
      vid(4)=0
      status=nf_def_var(ncid,'lat',NF_FLOAT,1,vid,dimid)
      status=nf_put_att_text(ncid,dimid,'units',13,"degrees_north")
      status=nf_put_att_text(ncid,dimid,'long_name',8,"Latitude")
      status=nf_put_att_text(ncid,dimid,'axis',1,'Y')
      write(*,*) "ID LAT: ",dimid

      vid(1)=4
      vid(2)=0
      vid(3)=0
      vid(4)=0
      status=nf_def_var(ncid,'lon',NF_FLOAT,1,vid,dimid)
      status=nf_put_att_text(ncid,dimid,'units',12,"degrees_east")
      status=nf_put_att_text(ncid,dimid,'long_name',9,"Longitude")
      status=nf_put_att_text(ncid,dimid,'axis',1,'X')
      write(*,*) "ID LON: ",dimid

      vid(1)=4
      vid(2)=3
      vid(3)=2
      vid(4)=1
      status=nf_def_var(ncid,variable,NF_FLOAT,4,vid,dimid)
      status=nf_put_att_text(ncid,dimid,'name',3,"AT ")
      status=nf_put_att_text(ncid,dimid,'long_name',14,"Temp. aire (C)")
      status=nf_put_att_text(ncid,dimid,'units',4,unidades)
      status=nf_put_att_real(ncid,dimid,'missing_value',NF_REAL,1,undef)
      write(*,*) "ID DATA: ",dimid,variable

C =====================================
C   Ingresamos los atributos globales
C =====================================

      status=nf_put_att_text(ncid,NF_GLOBAL,'Source',58,
     &"Instituto Geofisico del Peru ( http://www.met.igp.gob.pe )" )

      status=nf_put_att_text(ncid,NF_GLOBAL,'File_info',50,
     &"Contactar a Kobi Mosquera  kobi.mosquera@igp.gob.pe ")

      status=nf_put_att_text(ncid,NF_GLOBAL,'title',37,titulo)

      status=nf_put_att_text(ncid,NF_GLOBAL,'description',35,
     &descripcion)

      status=nf_put_att_real(ncid,NF_GLOBAL,'missing_value',
     &NF_REAL,1,undef)

      status=nf_put_att_text(ncid,NF_GLOBAL,'Convention',6,
     &"COARDS")

      status=nf_enddef(ncid)

C =====================================
C    Grabamos la data !!!!!
C =====================================

      status=nf_put_var_real(ncid,1,data_time1(1:nreg))
      status=nf_put_var_real(ncid,2,alt)
      status=nf_put_var_real(ncid,3,lat)
      status=nf_put_var_real(ncid,4,lon)
      status=nf_put_var_real(ncid,5,data(1:nreg,:))

      status=nf_close(ncid)
      end

C23456===============================================
      subroutine time_en_horas(iyf2,imf2,idf2,ihf2,horas)

      real horas
      character*10 sss
      integer iyi2,imi2,idi2,ihi2,ndias(12)
      integer iyf2,imf2,idf2,ihf2

      data iyi2/2000/,imi2/01/,idi2/01/,ihi2/00/
      data ndias/31,28,31,30,31,30,31,31,30,31,30,31/

      horas=0.
      do iy=iyi2,iyf2-1

      resto1=mod(iy,4)
      resto2=mod(iy,100)
      resto3=mod(iy,400)
      if (((resto1.eq.0).and.(resto2.ne.0)).or.(resto3.eq.0)) then
         ndias(2)=29
      else
         ndias(2)=28
      endif

      do im=1,12
         horas=horas+float(ndias(im))*24.
      enddo
      enddo

      do im=1,imf2-1
      resto1=mod(iyf2,4)
      resto2=mod(iyf2,100)
      resto3=mod(iyf2,400)
      if (((resto1.eq.0).and.(resto2.ne.0)).or.(resto3.eq.0)) then
         ndias(2)=29
      else
         ndias(2)=28
      endif
         horas=horas+float(ndias(im))*24.
      enddo

      horas=horas+float(idf2-1)*24.

      horas=horas+ihf2

      end