c f95 -w -maxcontin=50 read_gfs.f

! This file contains the description of the mrf station data derived from
! the MRF-AVN runs intended for ARM single model comparison and other uses
! of the column models.
! record 1 ... 12 integers

      character*100 input
      integer   hour       !(e.g.00)                                       
      integer   month      !(e.g.11)
      integer   day        !(e.g.7)
      integer   year       !(e.g.2000)
      integer   nsfc       !# of surface variables (e.g. 25)
      integer   nflx       !# of flux variables (e.g.28)
      integer   nvar       !# of variables for each sounding (e.g.7)
      integer   levss      !# of vertical levels for each sounding (e.g.42)
      integer   npoint     !# of station points (e.g.50)
      integer   starting   !forecast hour (e.g.0)
      integer   ending     !forecast hour (e.g.48)
      integer   forecast   !output step (e.g.3)
      integer   yyyymmddhh

! record 2 ... vertical sounding levels (real*4)
      parameter (levs=64)       
      real*4    sigi(levs+1), sigl(levs)
!	where sigi denotes the interface and sigl the model levels
!	All vertical level variables start from the near surface level

!the following are looped over # of station points and forecast time
!record 3 ... surface variables (real*4)
      real*4	latitude  !of station (degree)
      real*4	longitude !of station (degree)
      real*4	zsfc      !(model surface height for the station) (m)
      real*4	psfc      !(model surface pressure for the station) (mPa)
      real*4    dpsdt     !(surface pressure tendency) (Pa/sec)
      real*4	tsfc      !(model surface temperature) (K)
      real*4	soilm1    !(first soil layer volumetric water content) ()
      real*4	soilm2    !(second soil layer volumetric water content) ()
      real*4	soilm3    !(third soil layer volumetric water content) ()
      real*4	soilm4    !(fourth soil layer volumetric water content) ()
      real*4	snow      !(water equivalent snow depth) (m)
      real*4	soilt1    !(first soil layer temperature) (K)
      real*4	soilt2    !(second soil layer temperature) (K)
      real*4	soilt3    !(third soil layer temperature) (K)
      real*4	soilt4    !(fourth soil layer temperature) (K)
      real*4	tg3       !(deep soil layer temperature, fixed for time) (K)
      real*4	z0        !(model surface roughness length) (cm)
      real*4	cv        !(model convective cloud fraction) ()
      real*4	cvb       !(model convective cloud base pressure) (?)
      real*4	cvt       !(model convective cloud top pressure) (?)
      real*4    albedo1   !(one of four albedo used in the model) ()
      real*4    albedo2   !(one of four albedo used in the model) ()
      real*4    albedo3   !(one of four albedo used in the model) ()
      real*4    albedo4   !(one of four albedo used in the model) ()
      real*4	slimsk    !(sea(0) land(1) and ice(2) flag)
      real*4	vegfrac   !(model vegetation fraction) ()
      real*4	f10m      !(model 10-meter to sigma level1 similarity profile ratio) ()
      real*4	canopy    !(canopy water content) (m)
      real*4 	vegtype   !(vegetation type 1-13) ()
      real*4 	soiltype  !(soil type 1-9) ()
      real*4 	vegfrac1  !(vegfrac type 1 used in radiation)
      real*4 	vegfrac2  !(vegfrac type 2 used in radiation)
      real*4    uustar     !()
      real*4    ffmm       !()
      real*4    ffhh       !()
      real*4    hice       !sea ice thickness)
      real*4    fice       !sea ice fraction)
      real*4    tprcp      !precipitation) (m)
      real*4    srflag     !snow-rain flag (1=snow, 0=rain)
      real*4    snwdpth    !actual snow depth over land (mm)
      real*4    slc1       !first layer total soil water
      real*4    slc2       !second layer total soil water
      real*4    slc3       !third layer total soil water
      real*4    slc4       !fourth layer total soil water
      real*4    shdmin     !minimum vegetation fraction
      real*4    shdmax     !maximum vegetation fraction
      real*4    slope      !slope type (integer 1-9)
      real*4    snoalb     !maximum snow albedo
      real*4    orog       !orography (m)
      real*4    tice       !orography (m)


!N.B. record 4 should be skipped if the number nflx in the header record is zero
!as this would be the case for GDAS files.
!record 4 ... flux type variables (real*4)
       real*4		lsmask   !(land(0) sea(1) mask)
       real*4		ustress   !(stress in zonal direction)   !(N/m**2)
       real*4		vstress   !(stress in meridional direction)   !(N/m**2)
       real*4		shf   !(sensible heat flux)   !(W/m**2)
       real*4		lhf   !(latent heat flux)   !(W/m**2)
       real*4		dlws   !(downward long wave flux at surface)   !(W/M**2)
       real*4		ulws   !(upward long wave flux at surface)   !(W/M**2)
       real*4		ulwt   !(upward long wave flux at toa)   !(W/M**2)
       real*4		uswt   !(upward short wave flux at toa)   !(W/M**2)
       real*4		usws   !(upward short wave flux at surface)   !(W/M**2)
       real*4		dsws   !(downward short wave flux at surface)   !(W/M**2)
       real*4		cldh   !(high cloud fraction)   !()
       real*4		cldm   !(mid cloud fraction)   !()
       real*4		cldl   !(low cloud fraction)   !()
       real*4		precip   !(precipitation accumulated between outputs)   !(kg/m**2)
       real*4		cnvprc   !(conv precipitation accumulated between outputs)(kg/m**2)
       real*4		gflx   !(ground heat flux)   !(W/m**2)
       real*4		u10   !(model derived 10-meter zonal wind)   !(m/s)
       real*4		v10   !(model derived 10-meter meridional wind)   !(m/s)
       real*4		t2   !(model derived 2-meter temperature)   !(K)
       real*4		q2   !(model derived 2-meter specific humidity)   !(g/g)
       real*4		psfc1  !(surface pressure again)   !(hPa)
       real*4		tmax   !(maximum temp between outputs)   !(K)
       real*4		tmin   !(minimum temp between outputs)   !(K)
       real*4		ugrw   !(zonal stress due to gravity wave drag)   !(W/m**2)
       real*4		vgrw   !(meridional stress due to gravity wave drag)   !(W/m**2)
       real*4		hpbl   !(model diagnosed planetary boundary layer depth)   !(m)
       real*4		albedo   !(effective model albedo..derived)   !()
       real*4		cldpbl   !(pbl cloud fraction)   !()

! All records that follows are real*4
! record 5 ... levs of model zonal wind velocity (m/s)
! record 6 ... levs of model meridional wind velocity (m/s)
! record 7 ... levs of model temperature (K)
! record 8 ... levs of model specific humidity (g/g)
! record 9 ... levs of model pressure (mPa)
       real*4 u_wind(levs)
       real*4 v_wind(levs)
       real*4 tl(levs)
       real*4 ql(levs)
       real*4 pl(levs)

! if(nvar.gt.5)
! record 10 ... levs of model derived omega (mPa/sec)
!             This is an inadvertant unit selection due to making two unit
!             conversions. mPa is milli-Pascal.
       real*4 omega(levs)

! record 11 ... levs of model derived dtdt (advection) (K/sec)
! record 12 ... levs of model derived dqdt (advection) (g/g/sec)
       real*4 dtdt(levs)
       real*4 dqdt(levs)

! if(nvar.gt.8)
! record 13 ... levs of model cloud water/ice concentration (g/g)
! record 14 ... levs of model derived cloud water tendency (advection) (g/g/sec)
! record 15 ... levs of model cloud fraction 
       real*4 wc(levs)
       real*4 dwc(levs)
       real*4 cld(levs)
!----------------------------------------------------------------------

      print*,"input file name?"
      read(5,*) input 
      open(10,file=input,form="unformatted",
     >        status="old")



      read(10) hour ,!(e.g.00)                                       
     >    month     ,!(e.g.11)
     >    day       ,!(e.g.7)
     >    year      ,!(e.g.2000)
     >    nsfc      ,!# of surface variables (e.g. 25)
     >    nflx      ,!# of flux variables (e.g.28)
     >    nvar      ,!# of variables for each sounding (e.g.7)
     >    levss     ,!# of vertical levels for each sounding (e.g.42)
     >    npoint    ,!# of station points (e.g.50)
     >    starting  ,!forecast hour (e.g.0)
     >    ending    ,!forecast hour (e.g.48)
     >    forecast   !output step (e.g.3)

      write(11,*) "h m d y nsfc nflx nvar lev np st end int"
      write(11,'(i8,i16, 20i8)') hour ,
     >    month     ,
     >    day       ,
     >    year      ,
     >    nsfc      ,
     >    nflx      ,
     >    nvar      ,
     >    levss     ,
     >    npoint    ,
     >    starting  ,
     >    ending    ,
     >    forecast   

      read(10) sigi,sigl
      write(11,*) "sigi=",sigi
      write(11,*) "sigl=",sigl

      ntime=(ending-starting)/forecast+1
      do 100 i=1,ntime
      do 100 j=1,npoint

      read(10)	latitude  !of station (degree)
     > ,	longitude !of station (degree)
     > ,	zsfc      !(model surface height for the station) (m)
     > ,	psfc      !(model surface pressure for the station) (mPa)
     > ,        dpsdt     !(surface pressure tendency) (Pa/sec)
     > ,	tsfc      !(model surface temperature) (K)
     > ,	soilm1    !(first soil layer volumetric water content) ()
     > ,	soilm2    !(second soil layer volumetric water content) ()
     > ,	soilm3    !(third soil layer volumetric water content) ()
     > ,	soilm4    !(fourth soil layer volumetric water content) ()
     > ,	snow      !(water equivalent snow depth) (m)
     > ,	soilt1    !(first soil layer temperature) (K)
     > ,	soilt2    !(second soil layer temperature) (K)
     > ,	soilt3    !(third soil layer temperature) (K)
     > ,	soilt4    !(fourth soil layer temperature) (K)
     > ,	tg3       !(deep soil layer temperature, fixed for time) (K)
     > ,	z0        !(model surface roughness length) (cm)
     > ,	cv        !(model convective cloud fraction) ()
     > ,	cvb       !(model convective cloud base pressure) (?)
     > ,	cvt       !(model convective cloud top pressure) (?)
     > ,        albedo1   !(one of four albedo used in the model) ()
     > ,        albedo2   !(one of four albedo used in the model) ()
     > ,        albedo3   !(one of four albedo used in the model) ()
     > ,        albedo4   !(one of four albedo used in the model) ()
     > ,	slimsk    !(sea(0) land(1) and ice(2) flag)
     > ,	vegfrac   !(model vegetation fraction) ()
     > ,	f10m      !(model 10-meter to sigma level1 similarity profile ratio) ()
     > ,	canopy    !(canopy water content) (m)
     > , 	vegtype   !(vegetation type 1-13) ()
     > , 	soiltype  !(soil type 1-9) ()
     > , 	vegfrac1  !(vegfrac type 1 used in radiation)
     > , 	vegfrac2  !(vegfrac type 2 used in radiation)
     > ,        uustar     !()
     > ,        ffmm       !()
     > ,        ffhh       !()
     > ,        hice       !sea ice thickness)
     > ,        fice       !sea ice fraction)
     > ,        tprcp      !precipitation) (m)
     > ,        srflag     !snow-rain flag (1=snow, 0=rain)
     > ,        snwdpth    !actual snow depth over land (mm)
     > ,        slc1       !first layer total soil water
     > ,        slc2       !second layer total soil water
     > ,        slc3       !third layer total soil water
     > ,        slc4       !fourth layer total soil water
     > ,        shdmin     !minimum vegetation fraction
     > ,        shdmax     !maximum vegetation fraction
     > ,        slope      !slope type (integer 1-9)
     > ,        snoalb     !maximum snow albedo
     > ,        orog       !orography (m)
     > ,        tice       !sea ice skin temperature


      write(11,*) "timeID=",i," stationID=",j   
      write(21,'("latitude ,longitude: ", 10f10.3)') latitude ,longitude 
      write(11,'("latitude ,longitude: ", 10f10.3)') latitude ,longitude 
      write(11,'(10e13.5)')	zsfc      
     > ,	psfc      
     > ,        dpsdt     
     > ,	tsfc      
     > ,	soilm1    
     > ,	soilm2    
     > ,	soilm3    
     > ,	soilm4    
     > ,	snow      
     > ,	soilt1    
     > ,	soilt2    
     > ,	soilt3    
     > ,	soilt4    
     > ,	tg3    
     > ,	z0        
     > ,	cv        
     > ,	cvb       
     > ,	cvt       
     > ,        albedo1   
     > ,        albedo2   
     > ,        albedo3   
     > ,        albedo4   
     > ,	slimsk    
     > ,	vegfrac   
     > ,	f10m      
     > ,	canopy    
     > , 	vegtype   
     > , 	soiltype
     > , 	vegfrac1 
     > , 	vegfrac2  
     > ,        uustar     
     > ,        ffmm      
     > ,        ffhh     
     > ,        hice       
     > ,        fice      
     > ,        tprcp    
     > ,        srflag     
     > ,        snwdpth   
     > ,        slc1     
     > ,        slc2    
     > ,        slc3   
     > ,        slc4       
     > ,        shdmin    
     > ,        shdmax   
     > ,        slope   
     > ,        snoalb     
     > ,        orog      
     > ,        tice    


!N.B. record 4 should be skipped if the number nflx in the header record is zero
!as this would be the case for GDAS files.
!record 4 ... flux type variables (real*4)
      if(nflx.ne.0) then
      read(10) 		lsmask   !(land(0) sea(1) mask)
     >  ,		ustress   !(stress in zonal direction)   !(N/m**2)
     >  ,		vstress   !(stress in meridional direction)   !(N/m**2)
     >  ,		shf   !(sensible heat flux)   !(W/m**2)
     >  ,		lhf   !(latent heat flux)   !(W/m**2)
     >  ,		dlws   !(downward long wave flux at surface)   !(W/M**2)
     >  ,		ulws   !(upward long wave flux at surface)   !(W/M**2)
     >  ,		ulwt   !(upward long wave flux at toa)   !(W/M**2)
     >  ,		uswt   !(upward short wave flux at toa)   !(W/M**2)
     >  ,		usws   !(upward short wave flux at surface)   !(W/M**2)
     >  ,		dsws   !(downward short wave flux at surface)   !(W/M**2)
     >  ,		cldh   !(high cloud fraction)   !()
     >  ,		cldm   !(mid cloud fraction)   !()
     >  ,		cldl   !(low cloud fraction)   !()
     >  ,		precip   !(precipitation accumulated between outputs)   !(kg/m**2)
     >  ,		cnvprc   !(conv precipitation accumulated between outputs)(kg/m**2)
     >  ,		gflx   !(ground heat flux)   !(W/m**2)
     >  ,		u10   !(model derived 10-meter zonal wind)   !(m/s)
     >  ,		v10   !(model derived 10-meter meridional wind)   !(m/s)
     >  ,		t2   !(model derived 2-meter temperature)   !(K)
     >  ,		q2   !(model derived 2-meter specific humidity)   !(g/g)
     >  ,		psfc1  !(surface pressure again)   !(hPa)
     >  ,		tmax   !(maximum temp between outputs)   !(K)
     >  ,		tmin   !(minimum temp between outputs)   !(K)
     >  ,		ugrw   !(zonal stress due to gravity wave drag)   !(W/m**2)
     >  ,		vgrw   !(meridional stress due to gravity wave drag)   !(W/m**2)
     >  ,		hpbl   !(model diagnosed planetary boundary layer depth)   !(m)
     >  ,		albedo   !(effective model albedo..derived)   !()
     >  ,		cldpbl   !(pbl cloud fraction)   !()


      write(12,*) "timeID=",i," stationID=",j   
      write(12,*)	"lsmask ", lsmask   
     >  ,"ustress ",		ustress   
     >  ,"vstress ",		vstress  
     >  ,"shf     ",		shf   
     >  ,"lhf     ",		lhf  
     >  ,"dlws    ",		dlws  
     >  ,"ulws    ",	ulws 
     >  ,"ulwt    ",		ulwt   
     >  ,"uswt    ",		uswt  
     >  ,		usws 
     >  ,		dsws
     >  ,		cldh   
     >  ,		cldm  
     >  ,		cldl 
     >  ,		precip  
     >  ,		cnvprc 
     >  ,		gflx  
     >  ,		u10  
     >  ,		v10 
     >  ,		t2 
     >  ,		q2
     >  ,		psfc1  
     >  ,		tmax  
     >  ,		tmin 
     >  ,		ugrw
     >  ,		vgrw   
     >  ,		hpbl  
     >  ,		albedo   
     >  ,		cldpbl  

      endif

       read(10) u_wind
       read(10) v_wind
       read(10) tl
       read(10) ql
       read(10) pl

      write(13,*) "timeID=",i," stationID=",j   
       write(13,'("u_wind", 10e13.5)') u_wind
       write(13,'("v_wind", 10e13.5)') v_wind
       write(13,'("tl    ", 10e13.5)') tl
       write(13,'("ql    ", 10e13.5)') ql
       write(13,'("pl    ", 10e13.5)') pl

      if(nvar.gt.5) then
       read(10) omega
       read(10) dtdt
       read(10) dqdt

      write(14,*) "timeID=",i," stationID=",j   
       write(14,'("omega ", 10e13.5)') omega
       write(14,'("dtdt  ", 10e13.5)') dtdt
       write(14,'("dqdt  ", 10e13.5)') dqdt

      endif

! if(nvar.gt.8)
! record 13 ... levs of model cloud water/ice concentration (g/g)
! record 14 ... levs of model derived cloud water tendency (advection) (g/g/sec)
! record 15 ... levs of model cloud fraction 
      if(nvar.gt.8) then
       read(10) wc
       read(10) dwc
       read(10) cld

      write(15,*) "timeID=",i," stationID=",j   
       write(15,'("wc  ", 10e13.5)') wc
       write(15,'("dwc ", 10e13.5)') dwc
       write(15,'("cld ", 10e13.5)') cld
      endif

 100  continue
      end