text:
!dis
!dis    Open Source License/Disclaimer, Forecast Systems Laboratory
!dis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
!dis
!dis    This software is distributed under the Open Source Definition,
!dis    which may be found at http://www.opensource.org/osd.html.
!dis
!dis    In particular, redistribution and use in source and binary forms,
!dis    with or without modification, are permitted provided that the
!dis    following conditions are met:
!dis
!dis    - Redistributions of source code must retain this notice, this
!dis    list of conditions and the following disclaimer.
!dis
!dis    - Redistributions in binary form must provide access to this
!dis    notice, this list of conditions and the following disclaimer, and
!dis    the underlying source code.
!dis
!dis    - All modifications to this software must be clearly documented,
!dis    and are solely the responsibility of the agent making the
!dis    modifications.
!dis
!dis    - If significant modifications or enhancements are made to this
!dis    software, the FSL Software Policy Manager
!dis    (softwaremgr@fsl.noaa.gov) should be notified.
!dis
!dis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
!dis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
!dis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
!dis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
!dis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
!dis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
!dis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
!dis
!dis

MODULE lapsprep_netcdf

! PURPOSE
! =======
! Module to contain the various output routines needed for lapsprep.
!
! SUBROUTINES CONTAINED
! =====================
! output_netcdf_format  - Used to support RAMS 4.x initializations
!
! REMARKS
! =======

!
! HISTORY
! =======
! 28 Nov 2000 -- Original -- Brent Shaw
! 14 Nov 2001 -- Modified -- John Snook

  USE setup
  USE laps_static
  USE date_pack

  PRIVATE
  PUBLIC output_netcdf_format
CONTAINS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE output_netcdf_format(pr,ht,tp,mr,uw,vw,ww,slp,spr,lwc,ice,rai,sno,pic,rh,snocov,tskin)

  ! Subroutine to output data in netcdf format.  The skin temp from
  ! LAPS (lsx/tgd) is used for SST.

  ! Note that the P array has a bogus 2001 mb value as the last entry
  ! to designate the surface for MM5 applications.  The surface values
  ! of the state variables are contained in the last layer of their
  ! respective 3D arrays.

  ! The u/v winds are grid-relative, and are rotated to make them true.

  ! P levels must be written as integer values

  ! RH must be written as a fraction

  ! 3d Variables are stored top down, and are flipped to bottom up.

  IMPLICIT NONE

  include 'netcdf.inc'

  REAL                   :: pr(z3+1)      !Pressure levels (mb)
  REAL                   :: uw(x,y,z3+1)  !U-component of wind wrt grid (m/s)
  REAL                   :: vw(x,y,z3+1)  !V-component of wind wrt grid (m/s)
  REAL                   :: ww(x,y,z3+1)  !W-component of wind (m/s)
  REAL                   :: tp(x,y,z3+1)  !Temperature (K)
  REAL                   :: ht(x,y,z3+1)  !Geopotential Height (m)
  REAL                   :: mr(x,y,z3+1)  !Mixing ratio (kg/kg)
  REAL                   :: slp(x,y)      !MSL Pressure (Pa)
  REAL                   :: spr(x,y)      !Surface Pressure (Pa)
  REAL                   :: lwc(x,y,z3)   !Cloud water mr (kg/kg)
  REAL                   :: ice(x,y,z3)   !Cloud ice mr (kg/kg)
  REAL                   :: rai(x,y,z3)   !Precip rain mr (kg/kg)
  REAL                   :: sno(x,y,z3)   !Precip snow mr (kg/kg)
  REAL                   :: pic(x,y,z3)   !Precip ice mr (kg/kg)
  REAL                   :: rh(x,y,z3+1)  !Relative Humidity (%)
  REAL                   :: snocov(x,y)   !Snow cover (fract)
  REAL                   :: tskin(x,y)    !Skin Temp (K)

  ! Local Variables

  real, allocatable              :: pr3(:)
  real, allocatable              :: ht3(:,:,:)
  real, allocatable              :: tp3(:,:,:)
  real, allocatable              :: mr3(:,:,:)
  real, allocatable              :: rh3(:,:,:)
  real, allocatable              :: uw3(:,:,:)
  real, allocatable              :: vw3(:,:,:)
  real, allocatable              :: ww3(:,:,:)
  real, allocatable              :: sht(:,:)
  real, allocatable              :: stp(:,:)
  real, allocatable              :: srh(:,:)
  real, allocatable              :: smr(:,:)
  real, allocatable              :: suw(:,:)
  real, allocatable              :: svw(:,:)
  real, allocatable              :: lwcf(:,:,:)
  real, allocatable              :: icef(:,:,:)
  real, allocatable              :: raif(:,:,:)
  real, allocatable              :: snof(:,:,:)
  real, allocatable              :: picf(:,:,:)
  real, allocatable              :: scvf(:,:)
  real, allocatable              :: tskf(:,:)

  INTEGER                        :: yyyyddd, valid_mm, valid_dd
  INTEGER                        :: icode,ncid,nid,idimid(3),start,count
  INTEGER                        :: i,j,k
  INTEGER*2                      :: short
  REAL, ALLOCATABLE              :: ut(:,:,:)
  REAL, ALLOCATABLE              :: vt(:,:,:)
  REAL                           :: lat0
  REAL*8                         :: reftime
  CHARACTER (LEN=256)            :: output_file_name
  CHARACTER (LEN=128)            :: agridtype
  CHARACTER (LEN=15)             :: date_string

  ! Separate surface data from upper air data.

  allocate(pr3(z3))
  pr3 = pr(1:z3)

  allocate(ht3(x,y,z3))
  ht3 = ht(1:x,1:y,1:z3)

  allocate(tp3(x,y,z3))
  tp3 = tp(1:x,1:y,1:z3)

  allocate(mr3(x,y,z3))
  mr3 = mr(1:x,1:y,1:z3)

  allocate(rh3(x,y,z3))
  rh3 = rh(1:x,1:y,1:z3)

  allocate(uw3(x,y,z3))
  uw3 = uw(1:x,1:y,1:z3)

  allocate(vw3(x,y,z3))
  vw3 = vw(1:x,1:y,1:z3)

  allocate(ww3(x,y,z3))
  ww3 = ww(1:x,1:y,1:z3)

  allocate(sht(x,y))
  sht = ht(1:x,1:y,z3+1)

  allocate(stp(x,y))
  stp = tp(1:x,1:y,z3+1)

  allocate(smr(x,y))
  smr = mr(1:x,1:y,z3+1)

  allocate(srh(x,y))
  srh = rh(1:x,1:y,z3+1)

  allocate(suw(x,y))
  suw = uw(1:x,1:y,z3+1)

  allocate(svw(x,y))
  svw = vw(1:x,1:y,z3+1)

  allocate(scvf(x,y))
  scvf = snocov(1:x,1:y)

  allocate(tskf(x,y))
  tskf = tskin(1:x,1:y)

  ! Flip 3d arrays.

  call flip_array(1,1,z3,pr3)
  call flip_array(x,y,z3,ht3)
  call flip_array(x,y,z3,tp3)
  call flip_array(x,y,z3,mr3)
  call flip_array(x,y,z3,rh3)
  call flip_array(x,y,z3,uw3)
  call flip_array(x,y,z3,vw3)

  if (hotstart) then
    allocate(lwcf(x,y,z3))
    lwcf = lwc
    call flip_array(x,y,z3,lwcf)

    allocate(icef(x,y,z3))
    icef = ice
    call flip_array(x,y,z3,icef)

    allocate(raif(x,y,z3))
    raif = rai
    call flip_array(x,y,z3,raif)

    allocate(snof(x,y,z3))
    snof = sno
    call flip_array(x,y,z3,snof)

    allocate(picf(x,y,z3))
    picf = pic
    call flip_array(x,y,z3,picf)
    call flip_array(x,y,z3,ww3)
  endif

  ! Build the output file name

  output_prefix = TRIM(laps_data_root)// '/lapsprd/lapsprep/cdf/LAPS'
  yyyyddd = valid_yyyy*1000 + valid_jjj
  CALL wrf_date_to_ymd(yyyyddd, valid_yyyy, valid_mm, valid_dd)
  WRITE(date_string,'(I4.4,"-",I2.2,"-",I2.2,"-",I2.2,I2.2)') &
          valid_yyyy, valid_mm, valid_dd, valid_hh, valid_min
  output_file_name = TRIM(output_prefix) // ':' // date_string

  !  Create LAPS model netcdf file.

  icode=nf_create(trim(output_file_name),nf_clobber,ncid)
  if (icode .ne. 0) then
     print *,'Could not open output file: ',trim(output_file_name)
     stop
  endif

  !  Define netcdf grid dimensions.

  icode=nf_def_dim(ncid,'nx',x,1)
  icode=nf_def_dim(ncid,'ny',y,2)
  icode=nf_def_dim(ncid,'nz',z3,3)
  icode=nf_def_dim(ncid,'len',128,4)
  icode=nf_def_dim(ncid,'pt',1,5)

  !  Define static and grid variables.

  icode=nf_def_var(ncid,'grid_type',nf_char,1,4,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',15,'grid projection')
  icode=nf_def_var(ncid,'Nx',nf_short,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',18,'number of x points')
  icode=nf_def_var(ncid,'Ny',nf_short,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',18,'number of y points')
  icode=nf_def_var(ncid,'Np',nf_short,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',25,'number of pressure levels')
  icode=nf_def_var(ncid,'Dx',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',16,'x grid increment')
  icode=nf_put_att_text(ncid,nid,'units',10,'kilometers')
  icode=nf_def_var(ncid,'Dy',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',16,'y grid increment')
  icode=nf_put_att_text(ncid,nid,'units',10,'kilometers')
  icode=nf_def_var(ncid,'Lat0',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',13,'pole latitude')
  icode=nf_put_att_text(ncid,nid,'units',13,'degrees north')
  icode=nf_def_var(ncid,'Lat1',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'first latitude')
  icode=nf_put_att_text(ncid,nid,'units',13,'degrees north')
  icode=nf_def_var(ncid,'Lat2',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',15,'second latitude')
  icode=nf_put_att_text(ncid,nid,'units',13,'degrees north')
  icode=nf_def_var(ncid,'Lon0',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'pole longitude')
  icode=nf_put_att_text(ncid,nid,'units',12,'degrees east')
  icode=nf_def_var(ncid,'SwLat',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',18,'Southwest latitude')
  icode=nf_put_att_text(ncid,nid,'units',13,'degrees north')
  icode=nf_def_var(ncid,'SwLon',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',19,'Southwest longitude')
  icode=nf_put_att_text(ncid,nid,'units',12,'degrees east')
  icode=nf_def_var(ncid,'NeLat',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',18,'Northeast latitude')
  icode=nf_put_att_text(ncid,nid,'units',13,'degrees north')
  icode=nf_def_var(ncid,'NeLon',nf_real,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',19,'Northeast longitude')
  icode=nf_put_att_text(ncid,nid,'units',12,'degrees east')
  icode=nf_def_var(ncid,'reftime',nf_double,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'reference time')
  icode=nf_put_att_text(ncid,nid,'units',35,'seconds since (1970-1-1 00:00:00.0)')
  icode=nf_def_var(ncid,'valtime',nf_double,1,5,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',10,'valid time')
  icode=nf_put_att_text(ncid,nid,'units',35,'seconds since (1970-1-1 00:00:00.0)')
  icode=nf_def_var(ncid,'level',nf_real,1,3,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',13,'level of data')
  icode=nf_put_att_text(ncid,nid,'units',2,'mb')
  icode=nf_def_var(ncid,'pr',nf_real,1,3,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',8,'pressure')
  icode=nf_put_att_text(ncid,nid,'units',7,'Pascals')

  !  Define surface variables.

  idimid(1)=1
  idimid(2)=2
  icode=nf_def_var(ncid,'sht',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'surface height')
  icode=nf_put_att_text(ncid,nid,'units',6,'meters')
  icode=nf_def_var(ncid,'spr',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',16,'surface pressure')
  icode=nf_put_att_text(ncid,nid,'units',7,'pascals')
  icode=nf_def_var(ncid,'slp',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',23,'mean sea-level pressure')
  icode=nf_put_att_text(ncid,nid,'units',7,'pascals')
  icode=nf_def_var(ncid,'stp',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',19,'surface temperature')
  icode=nf_put_att_text(ncid,nid,'units',6,'kelvin')
  icode=nf_def_var(ncid,'smr',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',20,'surface mixing ratio')
  icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
  icode=nf_def_var(ncid,'srh',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',25,'surface relative humidity')
  icode=nf_put_att_text(ncid,nid,'units',7,'percent')
  icode=nf_def_var(ncid,'suw',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'surface u-wind')
  icode=nf_put_att_text(ncid,nid,'units',14,'meters/seconds')
  icode=nf_def_var(ncid,'svw',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',14,'surface v-wind')
  icode=nf_put_att_text(ncid,nid,'units',14,'meters/seconds')
  icode=nf_def_var(ncid,'scv',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',10,'snow cover')
  icode=nf_put_att_text(ncid,nid,'units',8,'fraction')
  icode=nf_def_var(ncid,'tsk',nf_real,2,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',16,'skin temperature')
  icode=nf_put_att_text(ncid,nid,'units',6,'kelvin')

  !  Define upper-air variables.

  idimid(3)=3
  icode=nf_def_var(ncid,'ht',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',6,'height')
  icode=nf_put_att_text(ncid,nid,'units',6,'meters')
  icode=nf_def_var(ncid,'tp',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',11,'temperature')
  icode=nf_put_att_text(ncid,nid,'units',6,'kelvin')
  icode=nf_def_var(ncid,'mr',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',12,'mixing ratio')
  icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
  icode=nf_def_var(ncid,'rh',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',17,'relative humidity')
  icode=nf_put_att_text(ncid,nid,'units',7,'percent')
  icode=nf_def_var(ncid,'uw',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',19,'u-component of wind')
  icode=nf_put_att_text(ncid,nid,'units',13,'meters/second')
  icode=nf_def_var(ncid,'vw',nf_real,3,idimid,nid)
  icode=nf_put_att_text(ncid,nid,'long_name',19,'v-component of wind')
  icode=nf_put_att_text(ncid,nid,'units',13,'meters/second')

  if (hotstart) then
     icode=nf_def_var(ncid,'lwc',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',31  &
          ,'cloud liquid water mixing ratio')
     icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
     icode=nf_def_var(ncid,'ice',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',22  &
          ,'cloud ice mixing ratio')
     icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
     icode=nf_def_var(ncid,'sno',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',31  &
          ,'precipitating snow mixing ratio')
     icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
     icode=nf_def_var(ncid,'rai',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',31  &
          ,'precipitating rain mixing ratio')
     icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
     icode=nf_def_var(ncid,'pic',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',30  &
          ,'precipitating ice mixing ratio')
     icode=nf_put_att_text(ncid,nid,'units',5,'kg/kg')
     icode=nf_def_var(ncid,'ww',nf_real,3,idimid,nid)
     icode=nf_put_att_text(ncid,nid,'long_name',19  &
          ,'w-component of wind')
     icode=nf_put_att_text(ncid,nid,'units',13,'meters/second')
  endif

  icode=nf_enddef(ncid)

  !  Create reftime.

  call adate_to_i4time(laps_file_time,reftime)

  !  Fill static and grid variables.

  IF      ( grid_type(1:8)  .EQ. 'latlon'                   ) THEN
    agridtype='Lat-Lon'
  ELSE IF ((grid_type(1:24) .EQ. 'secant lambert conformal').OR. &
           (grid_type(1:28) .EQ. 'tangential lambert conformal')) THEN
    agridtype='Lambert-Conformal'
  ELSE IF ( grid_type(1:19) .EQ. 'polar stereographic'      ) THEN
    agridtype='Polar-Stereographic'
  ELSE
    PRINT '(A,A,A)','netcdf unsupported map projection: ', &
        TRIM(grid_type),'.  I quit.'
    STOP 'unsupported projection'
  END IF

  start=1
  count=len_trim(agridtype)
  icode=nf_inq_varid(ncid,'grid_type',nid)
  icode=nf_put_vara_text(ncid,nid,start,count,agridtype)
  short=x
  icode=nf_inq_varid(ncid,'Nx',nid)
  icode=nf_put_var_int2(ncid,nid,short)
  short=y
  icode=nf_inq_varid(ncid,'Ny',nid)
  icode=nf_put_var_int2(ncid,nid,short)
  short=z3
  icode=nf_inq_varid(ncid,'Np',nid)
  icode=nf_put_var_int2(ncid,nid,short)
  icode=nf_inq_varid(ncid,'Dx',nid)
  icode=nf_put_var_real(ncid,nid,dx*1000.)
  icode=nf_inq_varid(ncid,'Dy',nid)
  icode=nf_put_var_real(ncid,nid,dy*1000.)
  lat0=90.
  icode=nf_inq_varid(ncid,'Lat0',nid)
  icode=nf_put_var_real(ncid,nid,lat0)
  icode=nf_inq_varid(ncid,'Lat1',nid)
  icode=nf_put_var_real(ncid,nid,latin1)
  icode=nf_inq_varid(ncid,'Lat2',nid)
  icode=nf_put_var_real(ncid,nid,latin2)
  icode=nf_inq_varid(ncid,'Lon0',nid)
  icode=nf_put_var_real(ncid,nid,lov)
  icode=nf_inq_varid(ncid,'SwLat',nid)
  icode=nf_put_var_real(ncid,nid,la1)
  icode=nf_inq_varid(ncid,'SwLon',nid)
  icode=nf_put_var_real(ncid,nid,lo1)
  icode=nf_inq_varid(ncid,'NeLat',nid)
  icode=nf_put_var_real(ncid,nid,la2)
  icode=nf_inq_varid(ncid,'NeLon',nid)
  icode=nf_put_var_real(ncid,nid,lo2)
  icode=nf_inq_varid(ncid,'reftime',nid)
  icode=nf_put_var_double(ncid,nid,reftime)
  icode=nf_inq_varid(ncid,'valtime',nid)
  icode=nf_put_var_double(ncid,nid,reftime)
  icode=nf_inq_varid(ncid,'level',nid)
  icode=nf_put_var_real(ncid,nid,pr3)
  pr3=pr3*100.
  icode=nf_inq_varid(ncid,'pr',nid)
  icode=nf_put_var_real(ncid,nid,pr3)

  ! Rotate the surface u and v winds to true if polar grid projection.

  ALLOCATE (ut(x,y,1)) ! Array for true u-winds
  ALLOCATE (vt(x,y,1)) ! Array for true v-winds

  IF ((agridtype(1:3) .EQ. 'Pol').OR.(agridtype(1:3) .EQ. 'Lam')) THEN
    DO j=1,y
    DO i=1,x
      CALL uvgrid_to_uvtrue(suw(i,j),svw(i,j), &
                            ut(i,j,1),vt(i,j,1),&
                            lons(i,j))
    ENDDO
    ENDDO
  ELSE
    ut(:,:,1) = suw
    vt(:,:,1) = svw
  ENDIF

  !  Fill surface variables.

  icode=nf_inq_varid(ncid,'sht',nid)
  icode=nf_put_var_real(ncid,nid,sht)
  icode=nf_inq_varid(ncid,'spr',nid)
  icode=nf_put_var_real(ncid,nid,spr)
  icode=nf_inq_varid(ncid,'slp',nid)
  icode=nf_put_var_real(ncid,nid,slp)
  icode=nf_inq_varid(ncid,'stp',nid)
  icode=nf_put_var_real(ncid,nid,stp)
  icode=nf_inq_varid(ncid,'smr',nid)
  icode=nf_put_var_real(ncid,nid,smr)
  icode=nf_inq_varid(ncid,'srh',nid)
  icode=nf_put_var_real(ncid,nid,srh)
  icode=nf_inq_varid(ncid,'suw',nid)
  icode=nf_put_var_real(ncid,nid,ut)
  icode=nf_inq_varid(ncid,'svw',nid)
  icode=nf_put_var_real(ncid,nid,vt)
  icode=nf_inq_varid(ncid,'scv',nid)
  icode=nf_put_var_real(ncid,nid,scvf)
  icode=nf_inq_varid(ncid,'tsk',nid)
  icode=nf_put_var_real(ncid,nid,tskf)

  DEALLOCATE (ut)
  DEALLOCATE (vt)

  ! Rotate the u and v winds to true if polar grid projection.

  ALLOCATE (ut(x,y,z3)) ! Array for true u-winds
  ALLOCATE (vt(x,y,z3)) ! Array for true v-winds

  IF ((agridtype(1:3) .EQ. 'Pol').OR.(agridtype(1:3) .EQ. 'Lam')) THEN
    level_loop: DO k = 1,z3
    rotate_winds_j: DO j=1,y
    rotate_winds_i: DO i=1,x
      CALL uvgrid_to_uvtrue(uw3(i,j,k),vw3(i,j,k), &
                            ut(i,j,k),vt(i,j,k),&
                            lons(i,j))
    ENDDO rotate_winds_i
    ENDDO rotate_winds_j
    ENDDO level_loop
  ELSE
    ut = uw3
    vt = vw3
  ENDIF

  !  Fill upper air variables.

  icode=nf_inq_varid(ncid,'ht',nid)
  icode=nf_put_var_real(ncid,nid,ht3)
  icode=nf_inq_varid(ncid,'tp',nid)
  icode=nf_put_var_real(ncid,nid,tp3)
  icode=nf_inq_varid(ncid,'mr',nid)
  icode=nf_put_var_real(ncid,nid,mr3)
  icode=nf_inq_varid(ncid,'rh',nid)
  icode=nf_put_var_real(ncid,nid,rh3)
  icode=nf_inq_varid(ncid,'uw',nid)
  icode=nf_put_var_real(ncid,nid,ut)
  icode=nf_inq_varid(ncid,'vw',nid)
  icode=nf_put_var_real(ncid,nid,vt)

  if (hotstart) then
    icode=nf_inq_varid(ncid,'lwc',nid)
    icode=nf_put_var_real(ncid,nid,lwcf)
    icode=nf_inq_varid(ncid,'ice',nid)
    icode=nf_put_var_real(ncid,nid,icef)
    icode=nf_inq_varid(ncid,'rai',nid)
    icode=nf_put_var_real(ncid,nid,raif)
    icode=nf_inq_varid(ncid,'sno',nid)
    icode=nf_put_var_real(ncid,nid,snof)
    icode=nf_inq_varid(ncid,'pic',nid)
    icode=nf_put_var_real(ncid,nid,picf)
    icode=nf_inq_varid(ncid,'ww',nid)
    icode=nf_put_var_real(ncid,nid,ww3)
    DEALLOCATE(lwcf)
    deallocate(icef)
    deallocate(raif)
    deallocate(snof)
    deallocate(picf)
  endif

  DEALLOCATE (ut)
  DEALLOCATE (vt)
  deallocate (pr3)
  deallocate (ht3)
  deallocate (tp3)
  deallocate (mr3)
  deallocate (rh3)
  deallocate (uw3)
  deallocate (vw3)
  deallocate (ww3)
  deallocate (sht)
  deallocate (stp)
  deallocate (smr)
  deallocate (srh)
  deallocate (suw)
  deallocate (svw)
  deallocate (scvf)
  deallocate (tskf)

  ! Close the netcdf file.

  icode=nf_close(ncid)

  RETURN
  END SUBROUTINE output_netcdf_format

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE flip_array(nx,ny,nz,a)

  implicit none

  integer nx,ny,nz,i,j,k,kk

  real a(nx,ny,nz),tmp(nx,ny,nz/2)

  kk=nz+1
  do k=1,nz/2
     kk=kk-1
     do j=1,ny
     do i=1,nx
        tmp(i,j,k)=a(i,j,k)
        a(i,j,k)=a(i,j,kk)
     enddo
     enddo
  enddo

  kk=nz/2+1
  do k=(nz+1)/2+1,nz
     kk=kk-1
     do j=1,ny
     do i=1,nx
        a(i,j,k)=tmp(i,j,kk)
     enddo
     enddo
  enddo

  RETURN
  END SUBROUTINE flip_array

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE adate_to_i4time(adate,i4time)

  implicit none

  real*8 i4time

  integer*4 iyear,iday,ihour,imin,lp

  character*9 adate

  read(adate(1:2),'(i2)') iyear
  read(adate(3:5),'(i3)') iday
  read(adate(6:7),'(i2)') ihour
  read(adate(8:9),'(i2)') imin

  !  Valid for years 1960-2060.

  if (iyear .lt. 60) iyear = iyear + 100

  lp = (iyear + 3 - 60) / 4

  i4time = (iyear-60)  * 31536000  &
         + (iday-1+lp) * 86400     &
         + ihour       * 3600      &
         + imin        * 60

  RETURN
  END SUBROUTINE adate_to_i4time

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

END MODULE lapsprep_netcdf

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
MODULE mo_structured_datatypes
   !
   ! Written for UCLALES-SALSA (U-S).
   !
   ! Juha Tonttila, FMI, 2017
   !

   IMPLICIT NONE

   ! -------------------------------------------------------------
   ! Holds derived types for simple float arrays which are needed
   ! in order to use the unlimited polymorphic fields in the
   ! FieldArray and ArrayElement classes.

   ! The FloatArrayXX classes contain three elements. First, the pointer "d", which
   ! provides the access point to the data. "Alloc" is only for internal storage purposes
   ! and should not be directly accessed. The procedure pointer "onDemand" provides the
   ! optional possibility to link a specific subroutine for a specific instance of the
   ! FloatArrayXX class. Typical use example would be: for the instance containing bulk aerosol number
   ! concentration, one would link a subroutine calculating the bulk concentration from
   ! binned concentrations, and this subroutine would be invoked automatically, when writing
   ! the bulk concentration output.

   ! The constructors take 3 arguments:
   !   -  trgt: This is the TARGET data array for the FloatArray instance
   !   -  store: If TRUE, the field "alloc" will be allocated the memory required by "trgt"
   !             and the data is saved in "alloc" and the pointer "d" is associated with alloc.
   !             If False, the data in "trgt" is actually saved in an external TARGET array, and
   !             "d" is associated directly with "trgt".
   !             Thus, in the former case, the FieldArray instance is the storage of the data,
   !             and in the latter case, the FieldArray instance points to data stored somewhere
   !             else. In both casesm the data must be accessed by "d", i.e. myInstance%d(:,:,...)
   !   -  sub: OPTIONAL: the constructor can be given a subroutine name, which "onDemand" will be
   !           associated with. Repeating the example above, for FieldArray instance for bulk aerosol
   !           number concentration, a subroutine calculating the bulk number from binned number concentrations
   !           would be passed, and then the CALL myBulkNumber%onDemand(args) would calculate the bulk value.
   !           this ofcourse forces the external subroutines, especially their interfaces to be pretty generic
   !           and basically all the necessary data has to come from imports from other modules. For examples,
   !           see mo_derived_procedures.
   !
   !           The sub thing may be extended later to cover other than bulk values!

   TYPE FloatArray0d
      REAL, POINTER :: d => NULL()
      PROCEDURE(), NOPASS, POINTER :: onDemand => NULL()
      LOGICAL :: Initialized = .FALSE.
   END TYPE FloatArray0d
   INTERFACE FloatArray0d
      PROCEDURE :: FloatArray0d_constructor
   END INTERFACE

   TYPE FloatArray1d
      REAL, POINTER :: d(:) => NULL() ! This is used as the generic access name for the data
      PROCEDURE(), NOPASS, POINTER :: onDemand => NULL()
      LOGICAL :: Initialized = .FALSE.
   END TYPE FloatArray1d
   INTERFACE FloatArray1d
      PROCEDURE :: FloatArray1d_constructor
   END INTERFACE FloatArray1d

   TYPE FloatArray2d
      REAL, POINTER :: d(:,:) => NULL() ! This is used as the generic access name for the data
      PROCEDURE(), NOPASS, POINTER :: onDemand => NULL()
      LOGICAL :: Initialized = .FALSE.
   END TYPE FloatArray2d
   INTERFACE FloatArray2d
      PROCEDURE :: FloatArray2d_constructor
   END INTERFACE FloatArray2d

   TYPE FloatArray3d
      REAL, POINTER :: d(:,:,:) => NULL() ! This is used as the generic access name for the data
      PROCEDURE(), NOPASS, POINTER :: onDemand => NULL()
      LOGICAL :: Initialized = .FALSE.
   END TYPE FloatArray3d
   INTERFACE FloatArray3d
      PROCEDURE :: FloatArray3d_constructor
   END INTERFACE FloatArray3d

   TYPE FloatArray4d
      REAL, POINTER :: d(:,:,:,:) => NULL() ! This is used as the generic access name for the data
      PROCEDURE(), NOPASS, POINTER :: onDemand => NULL()
      LOGICAL :: Initialized = .FALSE.
   END TYPE FloatArray4d
   INTERFACE FloatArray4d
      PROCEDURE :: FloatArray4d_constructor
   END INTERFACE FloatArray4d

 CONTAINS

   ! ----------------------------------------------------
   FUNCTION FloatArray0d_constructor(trgt)
     IMPLICIT NONE
     TYPE(FloatArray0d), TARGET         :: FloatArray0d_constructor
     REAL, INTENT(in), TARGET, OPTIONAL :: trgt

     IF (PRESENT(trgt)) &
          FloatArray0d_constructor%d => trgt
     FloatArray0d_constructor%Initialized = .TRUE.

   END FUNCTION FloatArray0d_constructor

   !---------------------------------------------------------------------

   FUNCTION FloatArray1d_constructor(trgt)
     IMPLICIT NONE
     TYPE(FloatArray1d), TARGET         :: FloatArray1d_constructor
     REAL, INTENT(in), TARGET, OPTIONAL :: trgt(:)

     IF (PRESENT(trgt)) &
          FloatArray1d_constructor%d => trgt
     FloatArray1d_constructor%Initialized = .TRUE.

   END FUNCTION FloatArray1d_constructor

   ! ------------------------------------------------------------------

   FUNCTION FloatArray2d_constructor(trgt)
     IMPLICIT NONE
     TYPE(FloatArray2d), TARGET         :: FloatArray2d_constructor
     REAL, INTENT(in), TARGET, OPTIONAL :: trgt(:,:)

     IF (PRESENT(trgt)) &
          FloatArray2d_constructor%d => trgt
     FloatArray2d_constructor%Initialized = .TRUE.

   END FUNCTION FloatArray2d_constructor

   ! ----------------------------------------------------------------------------

   FUNCTION FloatArray3d_constructor(trgt)
     IMPLICIT NONE
     TYPE(FloatArray3d), TARGET         :: FloatArray3d_constructor
     REAL, INTENT(in), TARGET, OPTIONAL :: trgt(:,:,:)

     IF (PRESENT(trgt)) &
          FloatArray3d_constructor%d => trgt
     FloatArray3d_constructor%Initialized = .TRUE.

   END FUNCTION FloatArray3d_constructor

   ! ------------------------------------------------------------------------

   FUNCTION FloatArray4d_constructor(trgt)
     IMPLICIT NONE
     TYPE(FloatArray4d), TARGET         :: FloatArray4d_constructor
     REAL, INTENT(in), TARGET, OPTIONAL :: trgt(:,:,:,:)

     IF (PRESENT(trgt)) &
          FloatArray4d_constructor%d => trgt
     FloatArray4d_constructor%Initialized = .TRUE.

   END FUNCTION FloatArray4d_constructor

END MODULE mo_structured_datatypes

max_stars_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_issues_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_forks_repo_name: UCLALES-SALSA/UCLALES-SALSA


text: The East Asian Collection is a planting of trees and shrubs from temperate regions of the largest continent, Asia.  Three deodar cedars, each with a distinct weeping outline, stand guard by an evergreen mound of Ubame oak.   Young specimens of paper mulberry, Korean hackberry, and dawn redwood are growing vigorously here.  A new collection of Asian lilacs adds springtime color and fragrance to this group; additional plantings featuring the Asian ancestors of the cultivated hybrid tea rose are also planned. This garden is a popular place for picnics, with open lawns and lovely views of Lake Spafford. Cherry blossoms and fragrant daphne are standouts in winter and spring, while ginkgo and zelkova trees and ornamental grasses provide beautiful fall color.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C $Id: dmspsp.f,v 1.5 2008-07-27 00:17:08 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      REAL FUNCTION DMSPSP (SGNA,DEGS,MINS,SECS)
C
C FUNCTION TO CONVERT UNPACKED DMS TO PACKED DMS ANGLE
C SGNA : SIGN OF ANGLE
C DEGS : DEGREES PORTION OF ANGLE
C MINS : MINUTES PORTION OF ANGLE
C SECS : SECONDS PORTION OF ANGLE
C
      IMPLICIT REAL (A-Z)
      REAL SECS
      INTEGER DEGS,MINS
      CHARACTER*1 SGNA,NEG
      DATA CON1,CON2 /1000000.0E0,1000.0E0/
      DATA NEG /'-'/
C
      CON = DEGS * CON1 + MINS * CON2 + SECS
      IF (SGNA .EQ. NEG) CON = - CON
      DMSPSP = CON
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE getpresprofile
      USE vmec_main
      USE vsvd
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: i
      REAL(rprec) :: sigmin
C-----------------------------------------------

!
!       WRITE OVER THOMPSON DATA
!
      lpofr = .false.                   !!these data are at s-half nodes
      lpprof = .false.
      itse = 10
      pthommax = datathom(1)                  !!compute in final version
      sigMIN = 0.03*pthommax

      DO i = 1, itse
         rthom(i) = REAL(i - 1,rprec)/(itse - 1)
         datathom(i) = (pthommax - sigmin)*(1. - rthom(i))**2 + sigmin
         sigma_thom(i) = 0.2*pthommax
      END DO

      END SUBROUTINE getpresprofile


text:
	PROGRAM READ_CATALOG
C++************************************************************
C	PROGRAM READ_CATALOG
C To read the catalogue of shell galaxies
C
C I, IDENT : identification
C F, FIELD : number of the field in the ESO/SRC survey
C W, WEST  : distance of the object from the east edge of the plate
C N, NORTH : distance of the object from the south edge of the plate
C R, RA    : right ascension
C D, DEC   : declination
C
C E, 1, ENV   : environment (I: isolated, C: cluster, G: small group
C               (less than 10))
C C, 2, COMM  : comments on the shells and extenal features surrounding the
C	        galaxy.
C S, 3, SHAPE : description of the galaxy
C A, 4, AA    : other denomination (NGC or IC, or other)
C B, 5, BPHOT : photometry in circular apertures with 6 parameters
C               log(D25),log(aperture in 0.1'),log(D25)-log(apert.), V, B-V, U-B
C               each parameter is separated by a "tab" (Ctrl I)
C O, 6, OBS   : observations of this object.
C H, 7, HBIBLI: bibliography.
C
C PHOTOM(6,IDIM) : array of the photometric measurements
C
C JLP Version of 14-02-90
C--************************************************************
	PARAMETER (IDIM=200,NLINES=30)
	REAL*4 PHOTOM(6,IDIM)
	INTEGER*4 INDEX1,IDESCR
	CHARACTER*80 DESCR
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER TEST(IDIM)*3,OUTPUT(IDIM)*120,AAGAL*80
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

	OPEN(2,FILE='read_catalog.log',STATUS='unknown',
     1	ACCESS='SEQUENTIAL')
	PRINT 20
	WRITE(2,20)
20	FORMAT(' Program "READ_CATALOG",  Version of 14-02-90')

C Input the catalogue :
	CALL RD_INPUTCAT(NGALAXY)

C********************************************************************
C Menu :
C********************************************************************

88	PRINT 87
87	FORMAT(/,' Menu :',/,
     1	' 1: Galaxies in increasing right ascension',/,
     1	' 2: Galaxies in increasing field number',/,
     1	' 3: Extraction of a subset of this catalog',
     1	' (with a keyword)',/,
     1	' 4: Creation of a file with photometric measurements',/,
     1	'    in circular apertures',/,
     1	' 5: Creation of a correlation table (or TeX output)',/,
     1	' 6: Histogram',/,
     1	' 7: Diagram with two parameters',/,
     1	' 10: Exit',/,
     1	' Enter the number of the option you want:')
	READ(5,*) IOPT

	IF(IOPT.EQ.10)GO TO 8888

C------------------------------------------------------------------------
	IF(IOPT.EQ.1)THEN
	   NINDEX1=NGALAXY
	    DO I=1,NINDEX1
	      INDEX1(I)=I
	    END DO
C Writing the selection on unit 2
	   CALL RD_WRITE(NINDEX1)

C________________________________________________________
	ELSEIF(IOPT.EQ.2)THEN
	   DO I=1,NGALAXY
	     CHAR=FIELD(I)
	     TEST(I)=CHAR(1:3)
	   END DO
	  NINDEX1=NGALAXY
	  CALL CLASS_CHAR(TEST,INDEX1,NINDEX1)
C Writing the selection on unit 2
	  CALL RD_WRITE(NINDEX1)

C*******************************************************************
C Option 3 : Extracting a selection from the catalogue with a key-word
C*******************************************************************
	ELSEIF(IOPT.EQ.3)THEN
300	  PRINT 301
301	  FORMAT(' In which array do you want to search ?',/,
     1	' 1: Environment (E)',/,
     1	' 2: Comments (C)',/,
     1	' 3: Description of the galaxy (S)',/,
     1	' 4: Identification (NGC, IC, ARP or ESO number) (A)',/,
     1	' 5: Observations (O)',/,
     1	' 6: Bibliography (H)',/,
     1	' Enter the NUMBER you want :',$)
	  READ(5,*,ERR=300) ICHOICE
	  PRINT 302
302	  FORMAT(' Enter the keyword:')
	  READ(5,10) WORD
	  LWORD=INDEX(WORD,'  ')-1
	  PRINT *,' LENGTH :',LWORD

C Looking for the key-word:
	    PRINT *,' Option chosen: ',ICHOICE
	    PRINT *,' DESCR(1,1,3)',DESCR(1,1,3)
	    CALL SEARCH(WORD,LWORD,%DESCR(DESCR(1,1,ICHOICE)),
     1	IDESCR(1,ICHOICE),INDEX1,NINDEX1,NGALAXY)

C Output of the result :

	PRINT 967,WORD,LWORD,NINDEX1
	WRITE(2,967)WORD,LWORD,NINDEX1
967	FORMAT(/,' Output of a subset of the catalogue',/
     1	' KEY-WORD: ',A20,' (WITH THE FIRST',I4,' LETTERS)',/,
     1	5X,' NUMBER OF GALAXIES FOUND :',I5,/,/)

	  IF(NINDEX1.EQ.0)THEN
	    PRINT *,' The word is not in this field'
	  ELSE
C Writing the selection on unit 2
	    CALL RD_WRITE(NINDEX1)
	  ENDIF

C------------------------------------------------------------------
C Creation of a catalogue with photometric measurements :
C------------------------------------------------------------------
	ELSEIF(IOPT.EQ.4)THEN
	  PRINT *,' WHICH GALAXIES DO YOU WANT TO OUTPUT :'
	  PRINT *,' Identification (NGCxxx, ESOxxxx,...) OR "ALL" ?'
	  READ(5,10) WORD

	  IF(WORD(1:3).EQ.'ALL')THEN
	    ISTART=1
	    IEND=NGALAXY
	  ELSE
	    LWORD=8
	    CALL SEARCH1(WORD,LWORD,IDENT,INDEX1,NINDEX1,NGALAXY)
	      IF(NINDEX1.EQ.0)THEN
	        PRINT *,' There is not any galaxy with this name'
	        GO TO 88
	      ENDIF
	        IF(NINDEX1.EQ.1)PRINT *,'O.K.: one galaxy has been found'
	    PRINT *,DESCR(INDEX1(1),1,4)
	    ISTART=INDEX1(1)
	    IEND=ISTART

	  ENDIF

408	  PRINT *,' OUTPUT FILE ?'
	  READ(5,10) NAME
	  OPEN(4,FILE=NAME,STATUS='NEW',
     1	ACCESS='SEQUENTIAL',ERR=408)

	    DO 401 IGAL=ISTART,IEND
	       IF(IDESCR(IGAL,5).EQ.0)GO TO 401
	      WRITE(4,403)IDENT(IGAL)(1:8),DESCR(IGAL,1,4)
403	      FORMAT(/,3X,' IDENTIFICATION : ',A,
     1	6X,A)
	      WRITE(4,404)
404	      FORMAT(11X,'LOG(D25)',3X,'LOG(APER)',2X,
     1	'LOG(D/A)',7X,'V',9X,'B-V',8X,'U-B',/)

	      DO IK=1,IDESCR(IGAL,5)
	        BUFF=DESCR(IGAL,IK,5)
	        CALL READVALUES(BUFF,PHOTOM,IK)
	      END DO

	      DO IK=1,IDESCR(IGAL,5)
	        WRITE(4,405)(PHOTOM(KK,IK),KK=1,6)
	      END DO

401	    CONTINUE

405	  FORMAT(9X,6(F9.3,2X))
	  CLOSE(4)
	  GO TO 88

C*******************************************************************
C Option 5 : Creation of a correlation table
C*******************************************************************
	ELSEIF(IOPT.EQ.5)THEN
	  CALL RD_TABLE(NGALAXY)

C Return to the main menu :
	  GO TO 88
C*******************************************************************
C Option 6 : Histogram
C*******************************************************************
	ELSEIF(IOPT.EQ.6)THEN
	  CALL RD_HISTO(NGALAXY)

C Return to the main menu :
	  GO TO 88
C*******************************************************************
C Option 7 : Diagram with two parameters
C*******************************************************************
	ELSEIF(IOPT.EQ.7)THEN
	  CALL RD_DIAGRAM(NGALAXY)

C Return to the main menu :
	  GO TO 88

	ELSE
	  GO TO 88
	ENDIF

8888	PRINT *,' Output in "read_catalog.log"'

	CLOSE(2)
	STOP
	END

C******************************************************************
C	SUBROUTINE READVALUES(BUFF,PHOTOM,ILINE)
C To read the values of aperture photometry measurements
C and store the result in PHOTOM(6,IDIM)
C*******************************************************************
	SUBROUTINE READVALUES(BUFF,PHOTOM,ILINE)
	PARAMETER (IDIM=200)
	REAL*4 PHOTOM(6,IDIM),VALUE
	INTEGER*4 LIMIT(80)
	CHARACTER BUFF*80,WORD*80,TAB*1
	DATA TAB/'09'X/
10	FORMAT(A)

C LIMIT(IK) is the location of the IK th TAB.
	IK=1
	LIMIT(IK)=0

	DO I=1,80
	  IF(BUFF(I:I).EQ.TAB)THEN
	    IK=IK+1
	    LIMIT(IK)=I
	  ENDIF
	END DO

	NWORDS=MIN0(IK,6)
C Now the last limit:
	LIMIT(NWORDS+1)=80

C Sometimes  I put a * after the last value
C and the work is already done in that case
C Otherwise:
C	IF(IK.LT.NWORDS)THEN
C	DO I=LIMIT(NWORDS+1),80
C	  IF(BUFF(I:I).EQ.' ')THEN
C	    LIMIT(NWORDS+1)=I
C	    GO TO 200
C	  ENDIF
C	END DO
C	PRINT *,' I CANNOT FIND THE END OF THE LINE ...'
C200	CONTINUE
C	ENDIF

C Loop on all the words
	DO I=1,NWORDS

C Then we decode the words which lay between 2 limits
	  J1=LIMIT(I)+1
	  J2=LIMIT(I+1)-1
C	  ILENTH=J2-J1+1
C	  WORD(1:ILENTH)=BUFF(J1:J2)
C	  CALL ASCII_REAL(WORD,ILENTH,VALUE)
C	  PHOTOM(I,ILINE)=VALUE
	  READ(BUFF(J1:J2),*)PHOTOM(I,ILINE)
	END DO

	RETURN
	END
C-----------------------------------------------------------------------
C Subroutine RD_INPUTCAT to input the catalogue
C
C-----------------------------------------------------------------------
	SUBROUTINE RD_INPUTCAT(NGALAXY)
	PARAMETER (IDIM=200,NLINES=30)
	INTEGER*4 INDEX1,IDESCR
	CHARACTER*80 DESCR
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC
	CHARACTER CODE(7)*1
	DATA CODE/'E','C','S','A','B','O','H'/

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

	PRINT *,' NAME OF THE CATALOGUE ?'
	READ(5,10) NAME
	OPEN(1,FILE=NAME,STATUS='OLD',ACCESS='SEQUENTIAL')
	JMAX=IDIM*30

C********************************************************************
C Reading the catalogue :
C********************************************************************
	I=0
	DO J=1,JMAX

	  READ(1,10,END=9999)BUFF
C Ident
	  IF(BUFF(1:1).EQ.'I')THEN
	    I=I+1
	    IDENT(I)=BUFF(3:82)
C Field number
	  ELSEIF(BUFF(1:1).EQ.'F')THEN
	    FIELD(I)=BUFF(3:82)
C Declination
	  ELSEIF(BUFF(1:1).EQ.'D')THEN
	    DEC(I)=BUFF(3:82)
C Right Ascension
	  ELSEIF(BUFF(1:1).EQ.'R')THEN
	    RA(I)=BUFF(3:82)
C North position on the plate
	  ELSEIF(BUFF(1:1).EQ.'N')THEN
	    NORTH(I)=BUFF(3:82)
C West position
	  ELSEIF(BUFF(1:1).EQ.'W')THEN
	    WEST(I)=BUFF(3:82)
	  ENDIF

C Description:
	  DO K=1,7
	    IF(BUFF(1:1).EQ.CODE(K))THEN
	      IDESCR(I,K)=IDESCR(I,K)+1
	      IK=IDESCR(I,K)
	      DESCR(I,IK,K)=BUFF(3:82)
	    ENDIF
	  END DO

	END DO

9999	PRINT *,'END OF FILE'
	CLOSE(1)
	NGALAXY=I

	PRINT 80,NGALAXY
	WRITE(2,80)NGALAXY
80	FORMAT(/,2X,' TOTAL NUMBER OF GALAXIES IN THE CATALOGUE : ',
     1	I5,/)

	RETURN
	END
C--------------------------------------------------------------------------
C
C Printing all the objects which have their indices in the array "INDEX1"
C
C--------------------------------------------------------------------------
	SUBROUTINE RD_WRITE(NINDEX1)
	PARAMETER (IDIM=200,NLINES=30)
	INTEGER*4 INDEX1,IDESCR
	CHARACTER*80 DESCR
	CHARACTER*1 OUT(7)
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

	PRINT 1001
1001	FORMAT(' WHICH ARRAYS DO YOU WANT TO OUTPUT',/,
     1	' - ENVIRONMENT ? (N)',$)
	READ(5,10) OUT(1)
	PRINT 1002
1002	FORMAT(' - COMMENTS ?(N)',$)
	READ(5,10) OUT(2)
	PRINT 1003
1003	FORMAT(' - DESCRIPTION OF THE GALAXY ?(N)',$)
	READ(5,10) OUT(3)
	PRINT 1004
1004	FORMAT(' - NGC OR IC NUMBER ?(N)',$)
	READ(5,10) OUT(4)
	PRINT 1005
1005	FORMAT(' - PHOTOMETRIC MEASUREMENTS ?(N)',$)
	READ(5,10) OUT(5)
	PRINT 1006
1006	FORMAT(' - OBSERVATIONS ?(N)',$)
	READ(5,10) OUT(6)
	PRINT 1007
1007	FORMAT(' - BIBLIOGRAPHY ?(N)',$)
	READ(5,10) OUT(7)

C**************************************************************
C Now writing in the file...

	WRITE(2,1020)
1020	FORMAT(2X,' IDENT  ',3X,'FIELD',2X,' W ',3X,' N ',3X,
     1	'RIGHT ASCENSION',3X,'  DECLINATION  ')

	DO J=1,NINDEX1
	  I=INDEX1(J)
	  WRITE(2,1202) IDENT(I)(1:8),FIELD(I)(1:3),
     1	WEST(I)(1:3),NORTH(I)(1:3),RA(I)(1:16),DEC(I)(1:16)
1202	FORMAT(/,2X,A8,3X,A3,4X,A3,3X,A3,3X,A15,3X,A15)

	  DO K=1,7
	    IF((OUT(K).EQ.'Y'.OR.OUT(K).EQ.'y')
     1	.AND.IDESCR(I,K).NE.0.AND.K.NE.5)THEN
	       DO IK=1,IDESCR(I,K)
	         WRITE(2,1203) DESCR(I,IK,K)(1:68)
1203	         FORMAT(12X,A68)
	       END DO
	    ENDIF
	  END DO

	  IF((OUT(5).EQ.'Y'.OR.OUT(5).EQ.'y')
     1	.AND.IDESCR(I,5).NE.0)THEN
	    WRITE(2,1204)
1204	    FORMAT(/,13X,'LOG(D25)',X,'LOG(APER)',X,
     1	'LOG(D/A)',3X,'V',6X,'B-V',5X,'U-B',/)

	    DO IK=1,IDESCR(I,5)
	      WRITE(2,1205) DESCR(I,IK,5)
1205	      FORMAT(14X,A80)
	    END DO

            WRITE(2,1201)
1201 	    FORMAT(' ')
	  ENDIF

	END DO

	RETURN
	END
C*******************************************************************
C Creation of a correlation table
C*******************************************************************
	SUBROUTINE RD_TABLE(NGALAXY)
	PARAMETER (IDIM=200,NLINES=30)
	REAL*4 PHOTOM(6,IDIM),VALUES_OUT(IDIM)
	INTEGER*4 INDEX1,IDESCR
	LOGICAL OUT_COORD
	CHARACTER*80 DESCR
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER OUTPUT(IDIM)*120,AAGAL*80
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC
	CHARACTER COORD(2)*26

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

	  IBEGIN=1

C Clearing out OUTPUT :
	   DO IGAL=1,NGALAXY
	    OUTPUT(IGAL)='&&   '
	   END DO

C Output of a TEX file :
508	  PRINT *,' Output TeX file ?'
	  READ(5,10) NAME
	  OPEN(9,FILE=NAME,STATUS='NEW',
     1	ACCESS='SEQUENTIAL',ERR=508)

C Header of the TeX file:
	  WRITE(9,562) NGALAXY
562	     FORMAT('\magnification=\magstep1',/,
     1	'\hsize=16.5truecm',/,
     1	'\vsize=24.truecm',/,
     1	'\voffset=-0.5truecm',/,
     1	'\hoffset=-0.5truecm',/,/,
     1	' Number of galaxies: ',I5,/)

	  PRINT 555
555	  FORMAT(' Do you want also to output the coordinates? (N)')
	  READ(5,10),ANS
	  OUT_COORD=((ANS.EQ.'Y').OR.(ANS.EQ.'y'))

	  PRINT 511
511	  FORMAT(' Number of keywords ? ',/,
     1	'(Enter 0 if you want to output only names and coordinates)')
	  READ(5,*) NKWORDS
	  NCOL=NKWORDS+2
	  IF(OUT_COORD)NCOL=NCOL+2

C**********************************************************************
C Loop on all the keywords to generate OUTPUT
	DO IKWORD=1,NKWORDS

500	  PRINT 501
501	  FORMAT(' In which array do you want to search ?',/,
     1	' 1: Environment (E)',/,
     1	' 2: Comments (C)',/,
     1	' 3: Description of the galaxy (S)',/,
     1	' 4: Identification (NGC, IC, ARP or ESO number) (A)',/,
     1	' 5: Observations (O)',/,
     1	' 6: Bibliography (H)',/,
     1	' Enter the NUMBER you want :',$)
	  READ(5,*,ERR=500) ICHOICE
	  PRINT 502
502	  FORMAT(' Which key-word are you looking for ?')
	  READ(5,10) WORD

	  LWORD=INDEX(WORD,'  ')-1
	  PRINT *,' LENGTH :',LWORD

	  PRINT 517
517	  FORMAT('     Two possibilities :',/,
     1	' 1. Correlation table filled',
     1	' with crosses when the item is found',/,
     1	' 2. Table with the values of the given parameter',/,
     1	'     Enter your choice :')
	  READ(5,*) ICHO1

C Writing crosses when the item is found in OUTPUT :
	  IF(ICHO1.EQ.1)THEN
	     CALL RD_CORREL1(WORD,LWORD,%DESCR(DESCR(1,1,ICHOICE)),
     1	IDESCR(1,ICHOICE),OUTPUT,NUMBER,IBEGIN,NGALAXY)
	     IBEGIN=IBEGIN+5
	  ELSE
C Writing the values in OUTPUT :
	     CALL RD_CORREL2(WORD,LWORD,%DESCR(DESCR(1,1,ICHOICE)),
     1	IDESCR(1,ICHOICE),OUTPUT,
     1	VALUES_OUT,NUMBER,IBEGIN,NGALAXY)
	     IBEGIN=IBEGIN+10
	  ENDIF

C Writing the keyword in the output TEX file :
	  PRINT 567,WORD,NUMBER
	  WRITE(2,567)WORD,NUMBER
567	  FORMAT(/,' Keyword: ',A20,
     1	5X,' NUMBER OF OBJECTS SELECTED :',I5,/,/)
C Writing in the begining of the output file:
	  WRITE(9,565)WORD,NUMBER
565	  FORMAT(' Keyword: ',A20,
     1	2X,' Objects selected:',I5,/)

	ENDDO
C**********************************************************************
C Writing on the TeX file now:

C Beginning of the table:
	     WRITE(9,566)
566	     FORMAT('\vbox{\tabskip=0pt\offinterlineskip',/,
     1	'\hrule',/,'\halign to \hsize{\tabskip=\centering',/,
     1	<NCOL>('\vrule#&\strut#\hfil&',/),
     1	'\vrule#\tabskip=0pt\cr',/,
     1	'height2pt%',/,<NCOL>('&\omit&',/),'\cr',/,
     1	'& Ident. \hfil&',/,<NCOL-1>('&---&',/),'\cr',/,
     1	'height2pt%',/,<NCOL>('&\omit&',/),'\cr',/,
     1	'\noalign{\hrule}',/,
     1	'height2pt%',/,<NCOL>('&\omit&',/),'\cr',/)

	     IEND=MIN(IBEGIN,115)

	     DO IGAL=1,NGALAXY

C Looking for the identification in NGC numbers :
	       AAGAL=' '
	       IMAX=INDEX(DESCR(IGAL,1,4),';')-1
	        IF(IMAX.GE.1)THEN
		  AAGAL=DESCR(IGAL,1,4)(1:IMAX)
	        ELSE
	          AAGAL='$\ldots$'
	        ENDIF

C Writing the final "Carriage return" for each record:
	       WRITE(OUTPUT(IGAL)(IEND:120),568)
568	       FORMAT(' &\cr')

C Cutting the lines in less than 80 characters
	        IF(OUTPUT(1)(56:57).EQ.'&&')THEN
	           IEND1=55
	           IEND2=IEND+5-IEND1
	        ELSE
	           IEND1=50
	           IEND2=IEND+5-IEND1
	        ENDIF
C Writing on one line if less than 80 characters, or two lines otherwise.
	        IF(IEND2.LE.0)THEN
	          IF(OUT_COORD)THEN
		    CALL RDC_DECODE_COORD(RA(IGAL),DEC(IGAL),COORD)
	            WRITE(9,549) IDENT(IGAL),AAGAL,COORD(1),COORD(2),
     1	    OUTPUT(IGAL)
549	            FORMAT('&',A8,'&&',A11,'&&',A26,'&&',A26,/,A<IEND+5>)
	          ELSE
	            WRITE(9,569) IDENT(IGAL),AAGAL,OUTPUT(IGAL)
569	            FORMAT('&',A8,'&&',A11,A<IEND+5>)
	          ENDIF
	        ELSE
	          IF(OUT_COORD)THEN
		    CALL RDC_DECODE_COORD(RA(IGAL),DEC(IGAL),COORD)
	            WRITE(9,540) IDENT(IGAL),AAGAL,COORD(1),COORD(2),
     1	    OUTPUT(IGAL),OUTPUT(IGAL)(IEND1+1:IEND+5)
540	            FORMAT('&',A8,'&&',A11,'&&',A26,'&&',A26,/,
     1	    A<IEND1>,/,A)
	          ELSE
	            WRITE(9,570) IDENT(IGAL),AAGAL,
     1	    OUTPUT(IGAL),OUTPUT(IGAL)(IEND1+1:IEND+5)
570	            FORMAT('&',A8,'&&',A11,A<IEND1>,/,A)
	          ENDIF
	        ENDIF

	     IF(IGAL.EQ.NGALAXY)THEN
	       WRITE(9,572)
572	       FORMAT('height2pt%',/,<NCOL>('&\omit&',/),'\cr}',
     1	/,'\hrule}',/,'\vfill\eject',/,'\bye')
	     ELSEIF(MOD(IGAL,40).EQ.0)THEN
	        WRITE(9,571)
571	        FORMAT('height2pt%',/,<NCOL>('&\omit&',/),'\cr}',
     1	/,'\hrule}',/,'\vfill\eject')
	        WRITE(9,566)
	     ENDIF

	     END DO

	   CLOSE(9)

	RETURN
	END
C*******************************************************************
C Histogram
C*******************************************************************
	SUBROUTINE RD_HISTO(NGALAXY)
	PARAMETER (IDIM=200,NLINES=30)
	REAL*4 PHOTOM(6,IDIM)
	REAL*4 VALUES_IN(IDIM),XOUT(IDIM),YOUT(IDIM)
	INTEGER*4 INDEX1,IDESCR
	CHARACTER*80 DESCR
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER OUTPUT(IDIM)*120,AAGAL*80
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC
	CHARACTER NCHAR*4,CHAR1*30,CHAR2*30,TITLE*40,PLOTDEV*32
        CHARACTER PCOLOR*30

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

500	  PRINT 501
501	  FORMAT(' In which array do you want to search ?',/,
     1	' 1: Environment (E)',/,
     1	' 2: Comments (C)',/,
     1	' 3: Description of the galaxy (S)',/,
     1	' 4: Identification (NGC, IC, ARP or ESO number) (A)',/,
     1	' 5: Observations (O)',/,
     1	' 6: Bibliography (H)',/,
     1	' Enter the NUMBER you want :',$)
	  READ(5,*) ICHOICE
	  PRINT 502
502	  FORMAT(' WHICH KEY-WORD ARE YOU LOOKING FOR ?')
	  READ(5,10) WORD

	  LWORD=INDEX(WORD,'  ')-1
	  PRINT *,' LENGTH :',LWORD

C Writing the values in OUTPUT :
	IBEGIN=1
	CALL RD_CORREL2(WORD,LWORD,%DESCR(DESCR(1,1,ICHOICE)),
     1	IDESCR(1,ICHOICE),OUTPUT,
     1	VALUES_IN,NPTS_IN,IBEGIN,NGALAXY)

C Computing Min and Max :
	PRINT *,' ',NPTS_IN,' VALUES FOUND'
	XMIN=VALUES_IN(1)
	XMAX=VALUES_IN(1)
	DO I=2,NPTS_IN
	  XMIN=AMIN1(XMIN,VALUES_IN(I))
	  XMAX=AMAX1(XMAX,VALUES_IN(I))
	END DO
	PRINT *,' MIN AND MAX :',XMIN,XMAX

C Generating the histogram:
	PRINT *,' MIN, MAX, BIN WIDTH ?'
	READ(5,*) VMIN,VMAX,BIN_WIDTH
	CALL RD_HISTO2(VALUES_IN,NPTS_IN,XOUT,YOUT,NPTS_OUT,
     1	VMIN,VMAX,BIN_WIDTH)

C Output of the curve:
	NCHAR='L'
        PCOLOR='Default'
	KCUR=1
	CHAR2=' N'
	PRINT *,' X LABEL ?'
	READ(5,10) CHAR1
	PRINT *,' TITLE ?'
	READ(5,10) TITLE
	PRINT *,' PLOT DEVICE ?'
	READ(5,10) PLOTDEV
	CALL NEWPLOT(XOUT,YOUT,NPTS_OUT,IDIM,KCUR,CHAR1,CHAR2,
     1	TITLE,NCHAR,PCOLOR,PLOTDEV,' ',' ')

	RETURN
	END
C----------------------------------------------------------------------
C Subroutine RD_HISTO2
C Creates an histogram VALUES_OUT from a series of input values VALUES_IN
C----------------------------------------------------------------------
	SUBROUTINE RD_HISTO2(VALUES_IN,NPTS_IN,XOUT,YOUT,
     1	NPTS_OUT,VMIN,VMAX,BIN_WIDTH)
	REAL*4 VALUES_IN(*),XOUT(*),YOUT(*)

C Computing NPTS_OUT:
	NPTS_OUT=2*(NINT((VMAX-VMIN)/BIN_WIDTH)+1)
	  IF(NPTS_OUT.EQ.0.OR.NPTS_OUT.GT.200)THEN
	    PRINT *,' ERROR : BIN_WIDTH IS WRONG'
	    NPTS_OUT=1
	    RETURN
	  ENDIF

	DO J=1,NPTS_OUT/2
C	  XOUT(J)=VMIN+(FLOAT(J)+0.5)*BIN_WIDTH
	  XOUT(2*J-1)=VMIN+FLOAT(J-1)*BIN_WIDTH
	  XOUT(2*J)=VMIN+FLOAT(J)*BIN_WIDTH
	  YOUT(2*J-1)=0.
	  YOUT(2*J)=0.
	END DO

	DO I=1,NPTS_IN
	  WORK=AMIN1(VALUES_IN(I),VMAX)
	  WORK=AMAX1(VALUES_IN(I),VMIN)
	  J=NINT((WORK-VMIN)/BIN_WIDTH)+1
	  YOUT(2*J)=YOUT(2*J)+1.
	  YOUT(2*J-1)=YOUT(2*J-1)+1.
	END DO

	RETURN
	END
C*******************************************************************
C Diagram
C Draws a diagram with two parameters
C*******************************************************************
	SUBROUTINE RD_DIAGRAM(NGALAXY)
	PARAMETER (IDIM=200,NLINES=30)
	REAL*4 PHOTOM(6,IDIM),VOUT(IDIM)
	REAL*4 XPLOT(IDIM),YPLOT(IDIM)
	INTEGER*4 INDEX1,IDESCR
	CHARACTER*80 DESCR
	CHARACTER BUFF*90,CHAR*80,CHAR3*3,ANS*1,WORD*20,NAME*30
	CHARACTER OUTPUT(IDIM)*120,AAGAL*80
	CHARACTER*80 IDENT,FIELD,WEST,NORTH,RA,DEC
	CHARACTER NCHAR*4,CHAR1*30,CHAR2*30,TITLE*40,PLOTDEV*32
        CHARACTER PCOLOR*30

	COMMON /RD_CATAL1/ INDEX1(IDIM),IDENT(IDIM),FIELD(IDIM),
     1	WEST(IDIM),NORTH(IDIM),RA(IDIM),DEC(IDIM),
     1	DESCR(IDIM,NLINES,7),IDESCR(IDIM,7)

10	FORMAT(A)

509	  PRINT *,' OUTPUT MONGO FILE ?'
	  READ(5,10) NAME
	  OPEN(10,FILE=NAME,STATUS='NEW',
     1	ACCESS='SEQUENTIAL',ERR=509)

	  PRINT 489
489	  FORMAT(/,'    ENTER THE PARAMETER YOU WANT IN X :')
	DO I=1,2
	  IBEGIN=1+(I-1)*19
500	  PRINT 501
501	  FORMAT(/,' In which array do you want to search ?',/,
     1	' 1: ENVIRONMENT (E)',/,
     1	' 2: COMMENTS (C)',/,
     1	' 3: DESCRIPTION OF THE GALAXY (S)',/,
     1	' 4: NGC OR IC NUMBER (A)',/,
     1	' 6: OBSERVATIONS (O)',/,
     1	' 7: BIBLIOGRAPHY (H)',/,
     1	' 8: Right ascension ',/,
     1	' 9: Declination ',/,
     1	' Enter the NUMBER you want :')
	  READ(5,*,ERR=500) ICHOICE

C Coordinates:
	IF(ICHOICE.EQ.8)THEN
	     NOUT=NGALAXY
	     DO J=1,NOUT
		CALL DECODE_RA(RA(J),XRA)
	        WRITE(OUTPUT(J)(IBEGIN:),29) XRA
29	FORMAT('&& ',G15.7,2X)
	     END DO
	ELSEIF(ICHOICE.EQ.9)THEN
	     NOUT=NGALAXY
	     DO J=1,NOUT
		CALL DECODE_DEC(DEC(J),XDEC)
	        WRITE(OUTPUT(J)(IBEGIN:),29) XDEC
	     END DO
	ELSE
	  PRINT 502
502	  FORMAT(' WHICH KEY-WORD ARE YOU LOOKING FOR ?')
	  READ(5,10) WORD
	  LWORD=INDEX(WORD,'  ')-1
	  PRINT *,' LENGTH :',LWORD
C Writing the values in OUTPUT:
	  CALL RD_CORREL2(WORD,LWORD,%DESCR(DESCR(1,1,ICHOICE)),
     1	IDESCR(1,ICHOICE),OUTPUT,
     1	VOUT,NOUT,IBEGIN,NGALAXY)

	ENDIF

	PRINT *,' '
	PRINT *,' ',NOUT,' VALUES FOUND'

C Looking for the Y values :
	  PRINT *,' '
	  PRINT *,' ENTER THE PARAMETER YOU WANT IN Y :'
	  PRINT *,' '
	END DO

C Reading the points:
	NPTS=0
	DO I=1,NGALAXY
	   READ(OUTPUT(I)(3:19),*,ERR=11) X1
	   READ(OUTPUT(I)(22:40),*,ERR=11) Y1
	   NPTS=NPTS+1
	   XPLOT(NPTS)=X1
	   YPLOT(NPTS)=Y1
11	END DO

	PRINT *,' NUMBER OF POINTS FOR THE DIAGRAM :',NPTS

C Output in Mongo file:
	CALL WR_MONGO(XPLOT,YPLOT,NPTS)

C Output of the curve:
	KCUR=1
	PRINT *,' PLOT DEVICE ? (NONE if graph is not wanted)'
	READ(5,10) PLOTDEV
	IF((PLOTDEV(1:2).EQ.'NO').OR.(PLOTDEV(1:2).EQ.'no')) RETURN
	PRINT *,' X LABEL ?'
	READ(5,10) CHAR1
	PRINT *,' Y LABEL ?'
	READ(5,10) CHAR2
	PRINT *,' TITLE ?'
	READ(5,10) TITLE
	PRINT *,' SYMBOL (42?)'
	READ(5,10) NCHAR
        PCOLOR='Default'
	CALL NEWPLOT(XPLOT,YPLOT,NPTS,IDIM,KCUR,CHAR1,CHAR2,
     1	TITLE,NCHAR,PCOLOR,PLOTDEV,' ',' ')

	RETURN
	END
C----------------------------------------------------------------------
C Output in Mongo file:
	SUBROUTINE WR_MONGO(XPLOT,YPLOT,NPTS)
	REAL*4 XPLOT(NPTS),YPLOT(NPTS)
	INTEGER IFMT,NPTS

281	PRINT 280
280	FORMAT(' Format for the output  (MONGO FILE):',/,
     1	' 1: (2(X,G14.7))',/,' 2: (2(X,I10))',/,
     1	' 3: (X,I10,X,G14.7)',/,' 4: (X,G14.7,X,I10)',/,
     1	' Enter the NUMBER of the format you want :')
	READ(5,*,ERR=281) IFMT

	DO I=1,NPTS
	   X1=XPLOT(I)
	   Y1=YPLOT(I)
	IF(IFMT.EQ.1)THEN
	   PRINT 291,X1,Y1
	   WRITE(10,291) X1,Y1
291	FORMAT(2(X,G14.7))
	ELSEIF(IFMT.EQ.2)THEN
	   PRINT 292,INT(X1),INT(Y1)
	   WRITE(10,292) INT(X1),INT(Y1)
292	FORMAT(2(X,I10))
	ELSEIF(IFMT.EQ.3)THEN
	   PRINT 293,INT(X1),Y1
	   WRITE(10,293) INT(X1),Y1
293	FORMAT(X,I10,X,G14.7)
	ELSEIF(IFMT.EQ.4)THEN
	   PRINT 294,X1,INT(Y1)
	   WRITE(10,294) X1,INT(Y1)
294	FORMAT(X,G14.7,X,I10)
	ENDIF
	END DO
C Closing Mongo file:
	   CLOSE(10)
	RETURN
	END
C--------------------------------------------------------------
	SUBROUTINE DECODE_DEC(DEC,XDEC)
	CHARACTER DEC*80
	REAL*4 XDEC
	INTEGER IDEG,IMIN,ISEC
	I=INDEX(DEC,'DEG')
C	J=INDEX(DEC,''')
	J=I+6
	K=INDEX(DEC,'"')
	READ(DEC(1:I-1),*) IDEG
	READ(DEC(I+3:J-1),*) IMIN
	READ(DEC(J+1:K-1),*) ISEC
	IF(INDEX(DEC(1:4),'-').NE.0)THEN
	  XDEC=IDEG-FLOAT(IMIN)/60.-FLOAT(ISEC)/3600.
	ELSE
	  XDEC=IDEG+FLOAT(IMIN)/60.+FLOAT(ISEC)/3600.
	ENDIF
	RETURN
	END
C------------------------------------------------------------------
	SUBROUTINE DECODE_RA(RA,XRA)
	CHARACTER RA*80
	REAL*4 XRA
	INTEGER IHOU,IMIN,ISEC
	I=INDEX(RA,'H')
	J=INDEX(RA,'MN')
	K=INDEX(RA,'SEC')
	READ(RA(1:I-1),*) IHOU
	READ(RA(I+1:J-1),*) IMIN
	READ(RA(J+2:K-1),*) ISEC
	XRA=IHOU+FLOAT(IMIN)/60.+FLOAT(ISEC)/3600.
C Conversion in degrees to compare with ESO coordinates:
	XRA=XRA*15.
	RETURN
	END
C--------------------------------------------------------------------
C For TeX output:
	SUBROUTINE RDC_DECODE_COORD(RA1,DEC1,COORD)
	CHARACTER RA1*80,DEC1*80,COORD(2)*26
C RA:
	I=INDEX(RA1,'H')
	J=INDEX(RA1,'MN')
	K=INDEX(RA1,'SEC')
	READ(RA1(1:I-1),*) IHOU
	READ(RA1(I+1:J-1),*) IMIN
	READ(RA1(J+2:K-1),*) SEC
	WRITE(COORD(1),21)IHOU,IMIN,SEC
21	FORMAT(I2,'$^H$ ',I2,'$^{m}$ ',F4.1,'$^{s}$')

C DEC:
	I=INDEX(DEC1,'DEG')
C	J=INDEX(DEC1,''')
	J=I+6
	K=INDEX(DEC1,'"')
	READ(DEC1(1:I-1),*) IDEG
	READ(DEC1(I+3:J-1),*) IMIN
	READ(DEC1(J+1:K-1),*) SEC

	IDEG=JIABS(IDEG)
	IF(INDEX(DEC1(1:4),'-').NE.0)THEN
	  WRITE(COORD(2),22)IDEG,IMIN,SEC
22	  FORMAT('$-$',I3,'$^o$ ',I2,''' ',F4.1,'" ')
	ELSE
	  WRITE(COORD(2),23)IDEG,IMIN,SEC
23	  FORMAT('$+$',I3,'$^o$ ',I2,''' ',F4.1,'" ')
	ENDIF

	RETURN
	END
C--------------------------------------------------------------------
	include 'jlpsub:search_set.for'
	include 'jlpsub:sort_set.for'

ext: for


text:
Stevenson Bridge Road is located to the west of Davis, about halfway between Davis and Winters.  Stevenson Bridge Road starts at Sievers Road, just a few miles west of Pedrick Produce, and runs due north until it crosses Putah Creek at Stevenson Bridge and becomes County Road 95A, about half a mile south of Russell Blvd. The road is very close to the start and finish line of the Wednesday Night Time Trials.

As of 20070501, the Solano County Roads Department is in the early stages of a significant improvement on Stevenson Bridge Road from just south of the bridge down to Campbell Road.  The project will take about one full month and includes the following components:  (1) widening the shoulders on the east and west side of the road by about 3 feet; (2)  digging out the severely degraded sections of the roadway surface and re paving them; (3) applying a final chip seal treatment for the entire length of the repaired portion  about 2 or 3 miles.


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
        subroutine ref_fill_vert(ref_3d_io,ni,nj,nk
     1           ,l_low_fill,l_high_fill
     1           ,lat,lon,topo,heights_3d,rlat_radar,rlon_radar
     1           ,rheight_radar,istatus)

!       Called from read_radar_3dref, hence from (cloud/deriv/accum/lplot -
!       if NOWRAD data is unavailable). Operates on V0x, for the purpose
!       of mosaicing and writing out the VRZ file, but not VRC files.
!       This Routine Fills in the gaps in a 3D radar volume that has been
!       stored as a sparse array. A linear interpolation in the vertical
!       direction is used. An additional option 'l_low_fill' can be set to
!       .true.. This will fill in low level echoes that might have been
!       missed because the grid points are below the radar horizon or too
!       close to the terrain. This is done if radar echo is detected not
!       too far above such grid points.
!
!       If the input column is entirely missing, so is the output column.
!       If there are QC flags set anywhere in the column, then the entire
!       output column is set to 'ref_base'. This implies that we are in the
!       zone of radar coverage and that "no echoes" is a good enough
!       characterization of the column.

!       Steve Albers            1990
!       Steve Albers            1992          Declare l_low_fill,l_high_fill
!       Steve Albers            1994          Set msg data to ref_base
!       Steve Albers            1998          Read in ref_base, r_missing_data
!            "              Feb 1998          Allows output of missing data
!            "                  2014          More efficient version

        use mem_namelist, ONLY: grid_spacing_m

!       ni,nj,nk are input LAPS grid dimensions
!       rlat_radar,rlon_radar,rheight_radar are input radar coordinates

        real ref_3d_io(ni,nj,nk)               ! I/O   3D reflctvy grid
        real heights_3d(ni,nj,nk)              ! I
        real lat(ni,nj),lon(ni,nj),topo(ni,nj) ! I     2D grids

        real ref_3d(ni,nj,nk)                  ! Local 3D reflctvy grid
        integer isum_ref_2d(ni,nj)             ! Local array

        logical l_low_fill, l_high_fill, l_test, l_nonmissing(ni,nj)

        write(6,*)' ref_fill_vert: begin'

        rpd = 3.14159265/180.

        call get_r_missing_data(r_missing_data,istatus)
        if(istatus .ne. 1)then
            write(6,*)' ref_fill_vert: Error reading r_missing_data'
            return
        endif

        call get_ref_base(ref_base,istatus)
        if(istatus .ne. 1)then
            write(6,*)' ref_fill_vert: Error reading ref_base'
            return
        endif

        thresh_low_fill_elev = 1.6     ! Elevation of echoes used to extrapolate
                                       ! downward (degress)

        thresh_low_fill_ht = 3000.     ! Height of echoes used to extrapolate
                                       ! downward (meters AGL)

        call latlon_to_rlapsgrid(rlat_radar,
     &                           rlon_radar,
     &                           lat,lon,
     &                           ni,nj,
     &                           ri,rj,
     &                           jstatus)
        if(jstatus .lt. 0)then
            write(6,*)' ref_fill_vert: Error radar lat/lon unavailable'
            return
        else
            write(6,*)' Radar ri/rj is ',ri,rj
        endif

        call get_grid_spacing_actual_xy(rlat_radar,rlon_radar
     1                                 ,grid_spacing_actual_mx
     1                                 ,grid_spacing_actual_my
     1                                 ,istatus)
        if(istatus .ne. 1)then
            write(6,*)
     1          ' ref_fill_vert: Error from get_grid_spacing_actual_xy'
            return
        endif

!       Range limit is 460km plus a ~10% cushion
        rgrid_radarx = 500000. / grid_spacing_actual_mx ! rectangle around radar
        rgrid_radary = 500000. / grid_spacing_actual_my ! rectangle around radar
        ilow = max(nint(ri-rgrid_radarx),1)
        jlow = max(nint(rj-rgrid_radary),1)
        ihigh = min(nint(ri+rgrid_radarx),ni)
        jhigh = min(nint(rj+rgrid_radary),nj)

!       Set missing values to ref_base for internal & external processing
!       write(6,*)' ref_fill_vert: Setting r_missing_data/qc values to '
!    1             ,ref_base

        ref_3d(:,:,:) = ref_base
        do k = 1,nk
        do j = jlow,jhigh
        do i = ilow,ihigh
            if(     ref_3d_io(i,j,k) .eq. -101.                 ! QC flags
     1         .or. ref_3d_io(i,j,k) .eq. -102.           )then
                ref_3d_io(i,j,k) = ref_base
            endif

            if(     ref_3d_io(i,j,k) .ne. r_missing_data  )then
                ref_3d(i,j,k) = ref_3d_io(i,j,k)
            endif
        enddo
        enddo
        enddo

        I4_elapsed = ishow_timer()

        write(6,*)
     1    ' ref_fill_vert: Interpolating vertically through gaps'
     1     ,ni,ilow,ihigh,nj,jlow,jhigh

        n_low_fill = 0
        n_high_fill = 0

!       Column sum value assuming all levels have 'ref_base'
        isum_test = nint(ref_base) * nk

        isum_ref_2d(:,:) = 0
        l_nonmissing(:,:) = .false.

        do k = 1,nk
        do j = jlow,jhigh
        do i = ilow,ihigh
            isum_ref_2d(i,j) = isum_ref_2d(i,j)
     1                       + nint(ref_3d(i,j,k))
            if(ref_3d_io(i,j,k) .ne. r_missing_data)then
                l_nonmissing(i,j) = .true.
            endif
        enddo
        enddo
        enddo

        I4_elapsed = ishow_timer()

        if(l_low_fill)then
          if(lat(1,1) .eq. 0. .or. lon(1,1)  .eq. 0.
     1                        .or. topo(1,1) .eq. r_missing_data)then
            write(6,*)' Error in ref_fill_vert:'
     1               ,' lat/lon/topo has invalid value'
            istatus = 0
            return
          endif
        endif

        write(6,*)
     1    ' ref_fill_vert: start main fill loop'

        n_columns = 0

        do j = jlow,jhigh
c       write(6,*)' Doing Column ',j

        do i = ilow,ihigh

          if(isum_ref_2d(i,j) .ne. isum_test)then ! Test for presence of echo

            n_columns = n_columns + 1

            call latlon_to_radar(lat(i,j),lon(i,j),topo(i,j)
     1                  ,azimuth,slant_range,elev_topo
     1                  ,rlat_radar,rlon_radar,rheight_radar)

            if(l_high_fill)then ! Fill in between high level echoes
                k = 1

!               Test for presence of top
15              l_test = .false.
                ref_below = ref_3d(i,j,k)

                do while (k .lt. nk .and. (.not. l_test))
                    ref_above = ref_3d(i,j,k+1)

!                   Should we require 'ref_above' to be missing in io array?
!                   This way we'd only fill in levels not scanned by the radar
                    if(ref_below .gt. ref_base
!    1           .and. ref_3d_io(i,j,k+1) .eq. r_missing_data ! rf_io_above
     1           .and. ref_above .eq. ref_base)then
                        k_top = k
                        l_test = .true.
                    endif

                    k = k + 1
                    ref_below = ref_above

                enddo

                if(.not. l_test)goto100 ! No Top exists

!               Echo Top exists, Search higher for next echo bottom
                l_test = .false.
                ref_below = ref_3d(i,j,k)

                do while(k .lt. nk .and. (.not. l_test))
                    k = k + 1
                    ref_above = ref_3d(i,j,k)

!                   Should we require 'ref_below' to be missing in io array?
!                   This way we'd only fill in levels not scanned by the radar
                    if(ref_above .gt. ref_base
!    1           .and. ref_3d_io(i,j,k-1) .eq. r_missing_data ! rf_io_below
     1           .and. ref_below .eq. ref_base)then
                        k_bottom = k
                        l_test = .true.
                    endif

                    ref_below = ref_above

                enddo

                if(.not. l_test)goto100 ! No Higher Echo Bottom exists

                if(heights_3d(i,j,k_top) .eq. r_missing_data)then
                    write(6,*)' ref_fill_vert: height is missing'
                    istatus = 0
                    return
                endif

!               Fill in gap if it exists and is small enough
                call latlon_to_radar(lat(i,j),lon(i,j)
     1                  ,heights_3d(i,j,k_top)
     1                  ,azimuth,slant_range,elev_top_deg
     1                  ,rlat_radar,rlon_radar,rheight_radar)

!               Adjust for larger tilt gaps in upper tilts
                if(elev_top_deg .gt. 9.0)then
                    tilt_gap_rad = 6.0 * rpd
                elseif(elev_top_deg .gt. 5.5)then
                    tilt_gap_rad = 5.0 * rpd
                else
                    tilt_gap_rad = 2.0 * rpd
                endif

!               If gap_thresh is larger then more filling will be done
                gap_thresh = max(2000., slant_range * tilt_gap_rad)

                if(heights_3d(i,j,k_bottom) .lt.
     1             heights_3d(i,j,k_top) + gap_thresh)then

                    n_high_fill = n_high_fill + 1

c                   write(6,101)(nint(max(ref_3d(i,j,kwrt),ref_base)),kwrt=1,nk)

                    do k_bet = k_top+1,k_bottom-1
                        frac = float(k_bet-k_bottom)
     1                        /float(k_top-k_bottom)
                        ref_3d(i,j,k_bet) =
     1                              ref_3d(i,j,k_bottom)*(1.-frac)
     1                            + ref_3d(i,j,k_top)*(frac)
                    enddo ! k

                    if(n_high_fill .le. 8)then
                        write(6,*)' Filled gap between'
     1                            ,i,j,k_bottom,k_top
     1                            ,elev_top_deg,gap_thresh
                        write(6,101)
     1                  (nint(max(ref_3d(i,j,kwrt),ref_base))
     1                                            ,kwrt=1,nk)
101                     format(1x,100i4)
                    endif

                endif

                if(k .le. nk-2)then ! Look for another gap
                    goto15
                endif

            endif ! l_high_fill

100         if(l_low_fill)then ! Fill below low level echoes

!               Search for bottom
                k_bottom = 0
                l_test = .false.

                k = 2

                ref_below = ref_3d(i,j,1)

                do while((.not. l_test) .and. k .le. nk)
                    ref_above = ref_3d(i,j,k)

                    if(ref_above .gt. ref_base
     1           .and. ref_below .eq. ref_base)then
                        k_bottom = k
                        dbz_bottom = ref_above
                        l_test = .true.
                    endif

                    k = k + 1
                    ref_below = ref_above

                enddo

                istatus = 1

                k_topo = max(int(height_to_zcoord(topo(i,j),istatus)),1)

                topo_buffer = topo(i,j) + thresh_low_fill_ht
                k_topo_buffer =
     1              max(int(height_to_zcoord(topo_buffer,istatus)),1)

                if(istatus .ne. 1)then
                    write(6,*)' ERROR return in ref_fill_vert'
                    return
                endif

                if(k_bottom .gt. k_topo)then ! Echo base above ground

                    height_bottom = height_of_level(k_bottom)

                    call latlon_to_radar(lat(i,j),lon(i,j),height_bottom
     1                  ,azimuth,slant_range,elev_bottom
     1                  ,rlat_radar,rlon_radar,rheight_radar)

                    if(k_bottom .le. k_topo_buffer ! Echo base near ground
!    1         .or. elev_bottom .lt. (elev_topo + thresh_low_fill_elev) ! Echo below radar horizon
     1                 .or. elev_bottom .lt. thresh_low_fill_elev ! Echo base near radar horizon
     1                                                      )then ! Fill in
                        do k = k_topo,k_bottom ! -1
                            ref_3d(i,j,k) = ref_3d(i,j,k_bottom)
                        enddo ! k

                        n_low_fill = n_low_fill + 1

                        if(n_low_fill .le. 8)then
!                       if(nint(azimuth) .eq. 270)then
!                       if(j .eq. 258)then
                            write(6,211)i,j
     1                                 ,k_topo,k_topo_buffer,k_bottom
     1                                 ,elev_bottom
     1                                 ,slant_range,ref_3d(i,j,k_bottom)
211                         format(' low fill   ',5i5,f6.2,2f8.0)
                        endif

                   else
!                       if(j .eq. 258)then
                        if(.false.)then
!                       if(nint(azimuth) .eq. 270)then
                            write(6,212)i,j
     1                                 ,k_topo,k_topo_buffer,k_bottom
     1                                 ,elev_bottom
     1                                 ,slant_range,ref_3d(i,j,k_bottom)
212                         format(' no low fill',5i5,f6.2,2f8.0)
                        endif

                   endif ! Fill in from ground to bottom of echo

                endif ! Bottom of radar echo above ground

            endif ! l_low_fill

          endif ! echo present

!         Return valid ref values to io array
          if(l_nonmissing(i,j))then
              do k = 1,nk
                  ref_3d_io(i,j,k) = ref_3d(i,j,k)
              enddo ! k
          endif ! We have non-missing reflectivity values in the column

        enddo ! i
        enddo ! j

        write(6,*)' n_gridpts = ',ni*nj
        write(6,*)' n_columns = ',n_columns
        write(6,*)' n_low_fill = ',n_low_fill
        write(6,*)' n_high_fill = ',n_high_fill

        istatus = 1
        return

!       This section is under construction - new formulation
        i_low_fill = 0
        if(slant_range .lt. 24000.)then ! find nearest level containing data
            do kdelt = 0,nk
            do i_sign = -1,+1,2
                k = k_topo_buffer + i_sign*kdelt
                if(k .ge. k_topo .and. k .le. nk .and.
     1             ref_3d(i,j,k) .ne. ref_base            )then
                    k_bottom = k
                    i_low_fill = 1
                    go to 200
                endif
            enddo ! i_sign
            enddo ! kdelt

 200        continue

        else ! slant_range > 24000.
             ! Pick lowest altitude > 500m agl having ref > 0 dbz
             ! and elev_angle < 2.0 deg (0.5 or 1.5)

        endif

        return
        end


text:
C * * * * * * * * * * * * * * *
C --- DRIVER FOR RADAR5
C * * * * * * * * * * * * * * *
C       Use make command (relevant to Makefile)

        IMPLICIT REAL*8 (A-H,O-Z)

        INTEGER, PARAMETER :: DP=kind(1D0)
C --->  PARAMETERS FOR RADAR5 (FULL JACOBIAN) <---
        INTEGER, PARAMETER :: ND=1
        INTEGER, PARAMETER :: NRDENS=1
        INTEGER, PARAMETER :: NGRID=0
        INTEGER, PARAMETER :: NLAGS=1
        INTEGER, PARAMETER :: NJACL=1
        INTEGER, PARAMETER :: MXST=5000
        INTEGER, PARAMETER :: LWORK=30
        INTEGER, PARAMETER :: LIWORK=30
        REAL, dimension(2) :: TARRAY
        REAL(kind=DP), dimension(ND) :: Y
        REAL(kind=DP), dimension(NGRID+1) :: GRID
        REAL(kind=DP), dimension(LWORK) :: WORK
        REAL(kind=DP), dimension(1) :: RTOL
        REAL(kind=DP), dimension(1) :: ATOL
        INTEGER, dimension(LIWORK) :: IWORK
        INTEGER, dimension(NRDENS+1) :: IPAST
        INTEGER, dimension(22) :: ISTAT
        EXTERNAL  FCN,PHI,ARGLAG,JACLAG,QFUN,SOLOUT

C ------ FILE TO OPEN ----------
        OPEN(9,FILE='sol.out')
        OPEN(10,FILE='cont.out')
        REWIND 9
        REWIND 10

        WRITE(*,*) 'INSERT TOL '
        READ(*,*) TOL
C -----------------------------------------------------------------------
C
C --- DIMENSION OF THE SYSTEM
        N=ND
C --- COMPUTE THE STANDARD JACOBIAN ANALYTICALLY
        IJAC=0
C --- JACOBIAN IS A FULL MATRIX
        MLJAC=N
C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
        IMAS=0
        MLMAS=N
C --- OUTPUT ROUTINE IS USED DURING INTEGRATION
        IOUT=1
C --- INITIAL VALUES
        X=2.0D0
        Y(1)=1.0D0
C       Consistent with initial function
C --- ENDPOINT OF INTEGRATION
        XEND=5.5D0
C       XEND=2.D0
C --- REQUIRED (RELATIVE AND ABSOLUTE) TOLERANCE
        ITOL=0
        RTOL=TOL
        ATOL=RTOL
C --- INITIAL STEP SIZE
        H=1.0D-2
C --- DEFAULT VALUES FOR PARAMETERS
        DO 10 I=1,20
        IWORK(I)=0
  10    WORK(I)=0.0D0
C --- ERROR CONTROL
        IWORK(11)=1
C --- STEPSIZE RATIOS
C       WORK(8)=0.5D0
C       WORK(9)=2.0D0
C --- WORKSPACE FOR PAST
        IWORK(12)=MXST
C --- BOTH COMPONENTS USE RETARDED ARGUMENT
        IWORK(15)=NRDENS
        IPAST(1)=1
C --- CONTROL OF NEWTON ITERATION
        IWORK(3)=10
        IWORK(14)=1
C --- GRID
        IWORK(13)=NGRID

C _____________________________________________________________________
C --- CALL OF THE SUBROUTINE RADAR5
        CALL DTIME(TARRAY)
        CALL RADAR5(N,FCN,PHI,ARGLAG,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JFCN,IJAC,MLJAC,MUJAC,
     &                  JACLAG,NLAGS,NJACL,
     &                  IMAS,SOLOUT,IOUT,
     &                  WORK,IWORK,RPAR,IPAR,IDID,
     &                  GRID,IPAST,QFUN,MLMAS,MUMAS)
        CALL DTIME(TARRAY)

C ---
        WRITE (6,*) 'Time = ',TARRAY(1)
        YT=-2.D0*LOG(EXP(-2.D0)*(1.D0-XEND+LOG(4.D0*EXP(4.D0))))
C --- PRINT FINAL SOLUTION SOLUTION
        WRITE (6,90) X,Y(1),ABS(Y(1)-YT)
C --- PRINT STATISTICS
         DO J=13,20
            ISTAT(J)=IWORK(J)
         END DO
        WRITE(6,*)' ***** TOL=',RTOL,' ****'
        WRITE (6,91) (ISTAT(J),J=14,20)
        WRITE (6,92) ISTAT(13)
 90     FORMAT(1X,'X =',F8.2,'    Y =',2E23.15)
 91     FORMAT(' fcn=',I7,' jac=',I6,' step=',I6,
     &        ' accpt=',I6,' rejct=',I6,' dec=',I6,
     &        ' sol=',I7)
 92     FORMAT(' full Newt. its =',I7)
        WRITE(6,*) 'SOLUTION IS TABULATED IN FILES: sol.out & cont.out'

        CLOSE(9)
        CLOSE(10)

        STOP
        END
C
C
        SUBROUTINE SOLOUT (NR,XOLD,X,HSOL,Y,CONT,LRC,N,
     &                     RPAR,IPAR,IRTRN)
C ----- PRINTS THE DISCRETE OUTPUT AND THE CONTINUOUS OUTPUT
C       AT EQUIDISTANT OUTPUT-POINTS
        IMPLICIT REAL*8 (A-H,O-Z)
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), PARAMETER :: XSTEP=0.002D0
        REAL(kind=DP), dimension(N) :: Y
        REAL(kind=DP), dimension(LRC) :: CONT
        REAL(kind=DP), dimension(1) :: RPAR
C       EXTERNAL PHI
C       XOUT IS USED FOR THE DENSE OUTPUT
        COMMON /INTERN/XOUT

        WRITE (9,99) X,Y(1)

        IF (NR.EQ.1) THEN
           WRITE (10,99) X,Y(1)
           XOUT=XSTEP
        ELSE
 10        CONTINUE
           IF (X.GE.XOUT) THEN
              WRITE (10,99) XOUT,CONTR5(1,N,XOUT,CONT,X,HSOL)
              XOUT=XOUT+XSTEP
              GOTO 10
           END IF
        END IF
 99     FORMAT(1X,'X =',F12.8,'    Y =',2E18.10)
        RETURN
        END
C
        DOUBLE PRECISION FUNCTION ARGLAG(IL,X,Y,RPAR,IPAR,
     &                                   PHI,PAST,IPAST,NRDS)
        IMPLICIT REAL*8 (A-H,O-Z)
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), dimension(1) :: Y
        REAL(kind=DP), dimension(1) :: PAST
        INTEGER, dimension(1) :: IPAST
        INTEGER, dimension(1) :: IPAR
        REAL(kind=DP), dimension(1) :: RPAR

        ARGLAG=Y(1)

        RETURN
        END
C
        SUBROUTINE FCN(N,X,Y,F,ARGLAG,PHI,RPAR,IPAR,
     &                  PAST,IPAST,NRDS)
        IMPLICIT REAL*8 (A-H,K,O-Z)
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), dimension(N) :: Y
        REAL(kind=DP), dimension(N) :: F
        REAL(kind=DP), dimension(1) :: PAST
        INTEGER, dimension(1) :: IPAST
        REAL(kind=DP), dimension(1) :: RPAR
        REAL RN
        EXTERNAL PHI

        CALL LAGR5(1,X,Y,ARGLAG,PAST,THETA,IPOS,RPAR,IPAR,
     &                 PHI,IPAST,NRDS)
        Y1L1=YLAGR5(1,THETA,IPOS,PHI,RPAR,IPAR,
     &              PAST,IPAST,NRDS)

        F(1)= Y1L1
        RETURN
        END
C

        SUBROUTINE JACLAG(N,X,Y,DFYL,ARGLAG,PHI,IVE,IVC,IVL,
     &                    RPAR,IPAR,PAST,IPAST,NRDS)
C ----- JACOBIAN OF DELAY TERMS IN THE EQUATION
        IMPLICIT REAL*8 (A-H,O-Z)
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), dimension(N) :: Y
        REAL(kind=DP), dimension(1) :: DFYL
        REAL(kind=DP), dimension(1) :: PAST
        INTEGER, dimension(1) :: IPAST
        REAL(kind=DP), dimension(1) :: RPAR
        INTEGER, dimension(1) :: IVE,IVC,IVL
        EXTERNAL PHI

        IVL(1)=1
        IVE(1)=1
        IVC(1)=1
        DFYL(1)=1.D0

        RETURN
        END

C
CD      SUBROUTINE PHI(I,X,VAL,RPAR,IPAR)
        DOUBLE PRECISION FUNCTION PHI(I,X,RPAR,IPAR)
        IMPLICIT REAL*8 (A-H,O-Z)
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), dimension(1) :: RPAR

        PHI=0.5D0

        RETURN
        END

        SUBROUTINE QFUN(N,Q,LQ,RPAR,IPAR)
C --- MATRIX "M" FOR THE TEST PROBLEM
        INTEGER, PARAMETER :: DP=kind(1D0)
        REAL(kind=DP), dimension(1) :: RPAR
        REAL(kind=DP), dimension(LQ,N) :: Q

        Q(1,1)=1.D0
        RETURN
        END


text:
c  To get dgamma,  "send dgamma from fnlib".
c  To get d1mach, mail netlib
c       send d1mach from core
c
      subroutine gaussq(kind, n, alpha, beta, kpts, endpts, b, t, w)
c
c           this set of routines computes the nodes t(j) and weights
c        w(j) for gaussian-type quadrature rules with pre-assigned
c        nodes.  these are used when one wishes to approximate
c
c                 integral (from a to b)  f(x) w(x) dx
c
c                              n
c        by                   sum w  f(t )
c                             j=1  j    j
c
c        (note w(x) and w(j) have no connection with each other.)
c        here w(x) is one of six possible non-negative weight
c        functions (listed below), and f(x) is the
c        function to be integrated.  gaussian quadrature is particularly
c        useful on infinite intervals (with appropriate weight
c        functions), since then other techniques often fail.
c
c           associated with each weight function w(x) is a set of
c        orthogonal polynomials.  the nodes t(j) are just the zeroes
c        of the proper n-th degree polynomial.
c
c     input parameters (all real numbers are in double precision)
c
c        kind     an integer between 1 and 6 giving the type of
c                 quadrature rule:
c
c        kind = 1:  legendre quadrature, w(x) = 1 on (-1, 1)
c        kind = 2:  chebyshev quadrature of the first kind
c                   w(x) = 1/sqrt(1 - x*x) on (-1, +1)
c        kind = 3:  chebyshev quadrature of the second kind
c                   w(x) = sqrt(1 - x*x) on (-1, 1)
c        kind = 4:  hermite quadrature, w(x) = exp(-x*x) on
c                   (-infinity, +infinity)
c        kind = 5:  jacobi quadrature, w(x) = (1-x)**alpha * (1+x)**
c                   beta on (-1, 1), alpha, beta .gt. -1.
c                   note: kind=2 and 3 are a special case of this.
c        kind = 6:  generalized laguerre quadrature, w(x) = exp(-x)*
c                   x**alpha on (0, +infinity), alpha .gt. -1
c
c        n        the number of points used for the quadrature rule
c        alpha    real parameter used only for gauss-jacobi and gauss-
c                 laguerre quadrature (otherwise use 0.d0).
c        beta     real parameter used only for gauss-jacobi quadrature--
c                 (otherwise use 0.d0)
c        kpts     (integer) normally 0, unless the left or right end-
c                 point (or both) of the interval is required to be a
c                 node (this is called gauss-radau or gauss-lobatto
c                 quadrature).  then kpts is the number of fixed
c                 endpoints (1 or 2).
c        endpts   real array of length 2.  contains the values of
c                 any fixed endpoints, if kpts = 1 or 2.
c        b        real scratch array of length n
c
c     output parameters (both double precision arrays of length n)
c
c        t        will contain the desired nodes.
c        w        will contain the desired weights w(j).
c
c     underflow may sometimes occur, but is harmless.
c
c     references
c        1.  golub, g. h., and welsch, j. h., "calculation of gaussian
c            quadrature rules," mathematics of computation 23 (april,
c            1969), pp. 221-230.
c        2.  golub, g. h., "some modified matrix eigenvalue problems,"
c            siam review 15 (april, 1973), pp. 318-334 (section 7).
c        3.  stroud and secrest, gaussian quadrature formulas, prentice-
c            hall, englewood cliffs, n.j., 1966.
c
c        original version 20 jan 1975 from stanford
c        modified 21 dec 1983 by eric grosse
c          imtql2 => gausq2
c          hex constant => d1mach (from core library)
c          compute pi using datan
c          removed accuracy claims, description of method
c          added single precision version
c
      double precision b(n), t(n), w(n), endpts(2), muzero, t1,
     x gam, solve, dsqrt, alpha, beta
c

ccc      print *, kind, n, alpha, beta, kpts, endpts

      call class (kind, n, alpha, beta, b, t, muzero)
c
c           the matrix of coefficients is assumed to be symmetric.
c           the array t contains the diagonal elements, the array
c           b the off-diagonal elements.
c           make appropriate changes in the lower right 2 by 2
c           submatrix.
c
      if (kpts.eq.0)  go to 100
      if (kpts.eq.2)  go to  50
c
c           if kpts=1, only t(n) must be changed
c
      t(n) = solve(endpts(1), n, t, b)*b(n-1)**2 + endpts(1)
      go to 100
c
c           if kpts=2, t(n) and b(n-1) must be recomputed
c
   50 gam = solve(endpts(1), n, t, b)
      t1 = ((endpts(1) - endpts(2))/(solve(endpts(2), n, t, b) - gam))
      b(n-1) = dsqrt(t1)
      t(n) = endpts(1) + gam*t1
c
c           note that the indices of the elements of b run from 1 to n-1
c           and thus the value of b(n) is arbitrary.
c           now compute the eigenvalues of the symmetric tridiagonal
c           matrix, which has been modified as necessary.
c           the method used is a ql-type method with origin shifting
c
  100 w(1) = 1.0d0
      do 105 i = 2, n
  105    w(i) = 0.0d0
c
      call gausq2 (n, t, b, w, ierr)
      do 110 i = 1, n
  110    w(i) = muzero * w(i) * w(i)
c
      return
      end
c
c
c
      double precision function solve(shift, n, a, b)
c
c       this procedure performs elimination to solve for the
c       n-th component of the solution delta to the equation
c
c             (jn - shift*identity) * delta  = en,
c
c       where en is the vector of all zeroes except for 1 in
c       the n-th position.
c
c       the matrix jn is symmetric tridiagonal, with diagonal
c       elements a(i), off-diagonal elements b(i).  this equation
c       must be solved to obtain the appropriate changes in the lower
c       2 by 2 submatrix of coefficients for orthogonal polynomials.
c
c
      double precision shift, a(n), b(n), alpha
c
      alpha = a(1) - shift
      nm1 = n - 1
      do 10 i = 2, nm1
   10    alpha = a(i) - shift - b(i-1)**2/alpha
      solve = 1.0d0/alpha
      return
      end
c
c
c
      subroutine class(kind, n, alpha, beta, b, a, muzero)
c
c           this procedure supplies the coefficients a(j), b(j) of the
c        recurrence relation
c
c             b p (x) = (x - a ) p   (x) - b   p   (x)
c              j j            j   j-1       j-1 j-2
c
c        for the various classical (normalized) orthogonal polynomials,
c        and the zero-th moment
c
c             muzero = integral w(x) dx
c
c        of the given polynomial's weight function w(x).  since the
c        polynomials are orthonormalized, the tridiagonal matrix is
c        guaranteed to be symmetric.
c
c           the input parameter alpha is used only for laguerre and
c        jacobi polynomials, and the parameter beta is used only for
c        jacobi polynomials.  the laguerre and jacobi polynomials
c        require the gamma function.
c
      double precision a(n), b(n), muzero, alpha, beta
      double precision abi, a2b2, dgamma, pi, dsqrt, ab
c
      pi = 4.0d0 * datan(1.0d0)
      nm1 = n - 1
      go to (10, 20, 30, 40, 50, 60), kind
c
c              kind = 1:  legendre polynomials p(x)
c              on (-1, +1), w(x) = 1.
c
   10 muzero = 2.0d0
      do 11 i = 1, nm1
         a(i) = 0.0d0
         abi = i
   11    b(i) = abi/dsqrt(4*abi*abi - 1.0d0)
      a(n) = 0.0d0
      return
c
c              kind = 2:  chebyshev polynomials of the first kind t(x)
c              on (-1, +1), w(x) = 1 / sqrt(1 - x*x)
c
   20 muzero = pi
      do 21 i = 1, nm1
         a(i) = 0.0d0
   21    b(i) = 0.5d0
      b(1) = dsqrt(0.5d0)
      a(n) = 0.0d0
      return
c
c              kind = 3:  chebyshev polynomials of the second kind u(x)
c              on (-1, +1), w(x) = sqrt(1 - x*x)
c
   30 muzero = pi/2.0d0
      do 31 i = 1, nm1
         a(i) = 0.0d0
   31    b(i) = 0.5d0
      a(n) = 0.0d0
      return
c
c              kind = 4:  hermite polynomials h(x) on (-infinity,
c              +infinity), w(x) = exp(-x**2)
c
   40 muzero = dsqrt(pi)
      do 41 i = 1, nm1
         a(i) = 0.0d0
   41    b(i) = dsqrt(i/2.0d0)
      a(n) = 0.0d0
      return
c
c              kind = 5:  jacobi polynomials p(alpha, beta)(x) on
c              (-1, +1), w(x) = (1-x)**alpha + (1+x)**beta, alpha and
c              beta greater than -1
c
   50 ab = alpha + beta
      abi = 2.0d0 + ab
      muzero = 2.0d0 ** (ab + 1.0d0) * dgamma(alpha + 1.0d0) * dgamma(
     x beta + 1.0d0) / dgamma(abi)
      a(1) = (beta - alpha)/abi
      b(1) = dsqrt(4.0d0*(1.0d0 + alpha)*(1.0d0 + beta)/((abi + 1.0d0)*
     1  abi*abi))
      a2b2 = beta*beta - alpha*alpha
      do 51 i = 2, nm1
         abi = 2.0d0*i + ab
         a(i) = a2b2/((abi - 2.0d0)*abi)
   51    b(i) = dsqrt (4.0d0*i*(i + alpha)*(i + beta)*(i + ab)/
     1   ((abi*abi - 1)*abi*abi))
      abi = 2.0d0*n + ab
      a(n) = a2b2/((abi - 2.0d0)*abi)
      return
c
c              kind = 6:  laguerre polynomials l(alpha)(x) on
c              (0, +infinity), w(x) = exp(-x) * x**alpha, alpha greater
c              than -1.
c
   60 muzero = dgamma(alpha + 1.0d0)
      do 61 i = 1, nm1
         a(i) = 2.0d0*i - 1.0d0 + alpha
   61    b(i) = dsqrt(i*(i + alpha))
      a(n) = 2.0d0*n - 1 + alpha
      return
      end
c
c
      subroutine gausq2(n, d, e, z, ierr)
c
c     this subroutine is a translation of an algol procedure,
c     num. math. 12, 377-383(1968) by martin and wilkinson,
c     as modified in num. math. 15, 450(1970) by dubrulle.
c     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
c     this is a modified version of the 'eispack' routine imtql2.
c
c     this subroutine finds the eigenvalues and first components of the
c     eigenvectors of a symmetric tridiagonal matrix by the implicit ql
c     method.
c
c     on input:
c
c        n is the order of the matrix;
c
c        d contains the diagonal elements of the input matrix;
c
c        e contains the subdiagonal elements of the input matrix
c          in its first n-1 positions.  e(n) is arbitrary;
c
c        z contains the first row of the identity matrix.
c
c      on output:
c
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct but
c          unordered for indices 1, 2, ..., ierr-1;
c
c        e has been destroyed;
c
c        z contains the first components of the orthonormal eigenvectors
c          of the symmetric tridiagonal matrix.  if an error exit is
c          made, z contains the eigenvectors associated with the stored
c          eigenvalues;
c
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
c
c     ------------------------------------------------------------------
c
      integer i, j, k, l, m, n, ii, mml, ierr
      real*8 d(n), e(n), z(n), b, c, f, g, p, r, s, machep
      real*8 dsqrt, dabs, dsign, d1mach
c
      machep=d1mach(4)
c
      ierr = 0
      if (n .eq. 1) go to 1001
c
      e(n) = 0.0d0
      do 240 l = 1, n
         j = 0
c     :::::::::: look for small sub-diagonal element ::::::::::
  105    do 110 m = l, n
            if (m .eq. n) go to 120
            if (dabs(e(m)) .le. machep * (dabs(d(m)) + dabs(d(m+1))))
     x         go to 120
  110    continue
c
  120    p = d(l)
         if (m .eq. l) go to 240
         if (j .eq. 30) go to 1000
         j = j + 1
c     :::::::::: form shift ::::::::::
         g = (d(l+1) - p) / (2.0d0 * e(l))
         r = dsqrt(g*g+1.0d0)
         g = d(m) - p + e(l) / (g + dsign(r, g))
         s = 1.0d0
         c = 1.0d0
         p = 0.0d0
         mml = m - l
c
c     :::::::::: for i=m-1 step -1 until l do -- ::::::::::
         do 200 ii = 1, mml
            i = m - ii
            f = s * e(i)
            b = c * e(i)
            if (dabs(f) .lt. dabs(g)) go to 150
            c = g / f
            r = dsqrt(c*c+1.0d0)
            e(i+1) = f * r
            s = 1.0d0 / r
            c = c * s
            go to 160
  150       s = f / g
            r = dsqrt(s*s+1.0d0)
            e(i+1) = g * r
            c = 1.0d0 / r
            s = s * c
  160       g = d(i+1) - p
            r = (d(i) - g) * s + 2.0d0 * c * b
            p = s * r
            d(i+1) = g + p
            g = c * r - b
c     :::::::::: form first component of vector ::::::::::
            f = z(i+1)
            z(i+1) = s * z(i) + c * f
  200       z(i) = c * z(i) - s * f
c
         d(l) = d(l) - p
         e(l) = g
         e(m) = 0.0d0
         go to 105
  240 continue
c
c     :::::::::: order eigenvalues and eigenvectors ::::::::::
      do 300 ii = 2, n
         i = ii - 1
         k = i
         p = d(i)
c
         do 260 j = ii, n
            if (d(j) .ge. p) go to 260
            k = j
            p = d(j)
  260    continue
c
         if (k .eq. i) go to 300
         d(k) = d(i)
         d(i) = p
         p = z(i)
         z(i) = z(k)
         z(k) = p
  300 continue
c
      go to 1001
c     :::::::::: set error -- no convergence to an
c                eigenvalue after 30 iterations ::::::::::
 1000 ierr = l
 1001 return
c     :::::::::: last card of gausq2 ::::::::::
      end


text: I work at the Wildlife Health Center, part of School of Veterinary Medicine VetMed at UC Davis, doing their publications (print and electronic). I keep a blog, http://magpienest.org/feathersofhope/ Feathers of Hope, with my partner (we are Pica and Numenius there) which explores nature and place (and some baseball: go Giants!). We were founding members of the http://www.magpienest.org/scgibin/wiki.pl Ecotone Wiki.