! Program Name: ! Author(s)/Contact(s): ! Abstract: ! History Log: ! ! Usage: ! Parameters: ! Input Files: ! ! Output Files: ! ! ! Condition codes: ! ! If appropriate, descriptive troubleshooting instructions or ! likely causes for failures could be mentioned here with the ! appropriate error code ! ! User controllable options: module module_ver implicit none type mapinfo_type integer :: ix integer :: jx character(len=40) :: projection real :: dx real :: dy real :: xlonc real :: reflat real :: reflon real :: refx real :: refy real :: truelat1 real :: truelat2 end type mapinfo_type type(mapinfo_type) :: map end module module_ver program ver use kwm_date_utilities use kwm_plot_utilities use module_ver implicit none integer, parameter :: maxstn = 200 integer, parameter :: maxtime = 2500 integer :: nstation character(len=4), dimension(maxstn) :: station real, dimension(maxstn) :: stnlat, stnlon, stnx, stny integer :: i, idt, indx ! character(len=10), parameter :: startdate = "2002051300" character(len=10), parameter :: startdate = "2002030100" character(len=10), parameter :: enddate = "2002063000" character(len=10) :: nowdate real, dimension(4, maxstn, maxtime) :: wc, sm character(len=256) :: wrfstatic ! wrfstatic filename call getarg(1, wrfstatic) print*, 'loc(station) = ', loc(station) ! Read the wrfstatic file to get map information call rdwrfstatic(wrfstatic) ! Find the OK Mesonet stations. call rdgeomeso(station, nstation, stnlat, stnlon, maxstn, 0) ! Find x/y do i = 1, nstation call lltoxy_generic(stnlat(i), stnlon(i), stnx(i), stny(i), & map%projection, 1.E-3*map%dx, map%reflat, map%reflon, map%refx, map%refy,& map%xlonc, map%truelat1, map%truelat2) enddo ! Loop over time call geth_newdate(nowdate, startdate, -1) DATELOOP: do while ( nowdate < enddate) call geth_newdate(nowdate, nowdate, 1) call geth_idts(nowdate, startdate, idt) indx = idt + 1 print*, 'nowdate = ', nowdate, indx ! Get station data from the OK Mesonet files do i = 1, nstation call rdsom(station(i), nowdate//"00", wc(1:4,i,indx)) enddo ! Get corresponding HRLDAS results from HRLDAS files call rdhrldas(station, stnlat, stnlon, nstation, nowdate, sm(1:4,1:nstation,indx)) enddo DATELOOP ! Make some plots call kwm_init_ncargks("sm.cgm") do i = 1, nstation call pltsm(station(i), wc(1:2,i,1:indx), sm(1:2,i,1:indx), 2, indx, startdate, & stnlat(i), stnlon(i), stnx(i), stny(i)) enddo call kwm_close_ncargks end program ver subroutine pltsm(stn, wc, sm, nlev, ndim, startdate, stnlat, stnlon, stnx, stny) use kwm_plot_utilities use kwm_date_utilities implicit none character(len=*), intent(in) :: stn integer, intent(in) :: nlev integer, intent(in) :: ndim real, dimension(nlev,ndim), intent(in) :: wc real, dimension(nlev,ndim), intent(in) :: sm character(len=*), intent(in) :: startdate real, intent(in) :: stnlat real, intent(in) :: stnlon real, intent(in) :: stnx real, intent(in) :: stny integer :: i,n character(len=256) :: text character(len=10) :: enddate call set(0., 1., 0., 1., 0., 1., 0., 1., 1) call pchiqu(0.14, 0.85, stn, 0.020, 0., -1.) write(text,'("Lat/Lon = ", F6.2, 1x, F7.2)') stnlat, stnlon call pchiqu(0.10, 0.81, trim(text), 0.012, 0., -1.) write(text,'("X/Y = ", F6.2, 1x, F7.2)') stnx, stny call pchiqu(0.10, 0.78, trim(text), 0.012, 0., -1.) text = "from "//startdate(1:4)//"-"//startdate(5:6)//"-"//startdate(7:8)//" "//startdate(9:10)//" Z" call pchiqu(0.44, 0.85, trim(text), 0.012, 0., -1.) call geth_newdate(enddate(1:10), startdate(1:10), ndim-1) text = "to "//enddate(1:4)//"-"//enddate(5:6)//"-"//enddate(7:8)//" "//enddate(9:10)//" Z" call pchiqu(0.44, 0.82, trim(text), 0.012, 0., -1.) call gasetc("YLF", '(F5.2)') call gasetc("XLF", '(I4)') do n = 1, nlev if (n == 1) then call set(.1, .9, .5, .75, 1., float(ndim), 0., 0.5, 1) else if (n == 2) then call set(.1, .9, .15, .4, 1., float(ndim), 0., 0.5, 1) endif call periml(0,0,5,1) call gsplci(color_index("blue")) do i = 1, ndim if (i == 1) then call frstpt(float(i), wc(n,i)) else call vector(float(i), wc(n,i)) endif enddo call plotif(0., 0., 2) call gsplci(color_index("red")) do i = 1, ndim if (i == 1) then call frstpt(float(i), sm(n,i)) else call vector(float(i), sm(n,i)) endif enddo call plotif(0., 0., 2) call gsplci(color_index("black")) enddo call frame() end subroutine pltsm subroutine rdgeomeso(station, nstation, stnlat, stnlon, maxstn, iverbose) implicit none integer, intent(in) :: maxstn character(len=4), dimension(maxstn), intent(out) :: station real, dimension(maxstn), intent(out) :: stnlat, stnlon integer, intent(out) :: nstation integer, intent(in) :: iverbose character(len=256) :: string integer :: ierr integer, parameter :: iunit = 10 station = " " stnlat = -999.0 stnlon = -999.0 nstation = 0 open(iunit, file="/scholar/kmanning/more_hrldas/OKMESO/geomeso.tbl", & status='old', form='formatted', action='read') PRELOOP : do read(iunit, '(A256)', iostat=ierr) string if (ierr /= 0) stop "PRELOOP" if (string(1:6) == "[STOP]") exit PRELOOP enddo PRELOOP RDLOOP : do read(iunit, '(A256)', iostat=ierr) string if (ierr /= 0) exit RDLOOP nstation = nstation + 1 read(string(6:9), '(A4)', iostat=ierr) station(nstation) if (ierr /= 0) stop "RDGEOMESO: problem extracting station name" read(string(64:81), '(F8.4,1x,F9.4)', iostat=ierr) stnlat(nstation), stnlon(nstation) if (ierr /= 0) stop "RDGEOMESO: problem extracting lat/lon" if (iverbose > 0) print*, station(nstation), stnlat(nstation), stnlon(nstation) enddo RDLOOP close(iunit) end subroutine rdgeomeso subroutine rdsom(stn, idate, wc) implicit none character(len=*), intent(in) :: stn, idate real, dimension(4), intent(out) :: wc integer :: ierr character(len=256) :: string, flnm character(len=12) :: rddate ! yyyymmddhhmm character(len=4) :: rdstn integer :: i real, dimension(4) :: st, tr, mp integer, dimension(4) :: stq, trq, mpq, wcq real :: tref integer, parameter :: iunit = 10 wc = -99999. ! write(flnm, & ! '("/wig/kmanning/HRLDAS_WRFDRIVER/VERIFICATION/OKMeso_Soil/",A8,".som")') & ! idate(1:8) write(flnm, & '("/wig/kmanning/HRLDAS_WRFDRIVER/VERIFICATION/OKMeso_SoilMoisture_SOM/",A8,".som")') & idate(1:8) open(iunit, file=trim(flnm), status='old', form='formatted', action='read') RDLOOP : do read(iunit, '(A256)', iostat=ierr) string if (ierr /= 0) exit RDLOOP read(string(6:17), '(A12)', iostat=ierr) rddate if (ierr /= 0) stop "Problem extracting date" if (rddate == idate) then read(string(1:4), '(A4)', iostat=ierr) rdstn if (ierr /= 0) stop "Problem extracting station name" if (rdstn == stn) then !KWM print*, trim(string) read(string(19:),*, iostat=ierr) ((st(i), stq(i), tr(i), trq(i), mp(i), mpq(i), wc(i), wcq(i)), i=1,4), & tref if (ierr /= 0) stop "Problem extracting data" !KWM print*, 'st = ', st !KWM print*, 'tr = ', tr !KWM print*, 'mp = ', mp !KWM print*, 'wc = ', wc !KWM print*, 'tref = ', tref exit RDLOOP endif endif enddo RDLOOP close(iunit) end subroutine rdsom subroutine rdwrfstatic(wrfstatic) use module_ver implicit none character(len=*), intent(in) :: wrfstatic #include integer :: ncid, iret, dimid real, dimension(16) :: corner_lat, corner_lon ! Open the NetCDF file, and get mapping information. write(*,'("wrfsi_static_flnm: ''", A, "''")') trim(wrfstatic) iret = nf_open(wrfstatic, NF_WRITE, ncid) if (iret /= 0) then write(*,'("Problem opening wrfsi_static file: ''", A, "''")') & trim(wrfstatic) stop endif iret = nf_inq_dimid(ncid, "west_east", dimid) if (iret /= 0) then stop "nf_inq_dimid: west_east" endif iret = nf_inq_dimlen(ncid, dimid, map%ix) if (iret /= 0) then stop "nf_inq_dimlen: west_east" endif iret = nf_inq_dimid(ncid, "south_north", dimid) if (iret /= 0) then stop "nf_inq_dimid: south_north" endif iret = nf_inq_dimlen(ncid, dimid, map%jx) if (iret /= 0) then stop "nf_inq_dimlen: south_north" endif iret = nf_inq_dimid(ncid, "land_cat", dimid) if (iret /= 0) then stop "nf_inq_dimid: land_cat" endif !KWM iret = nf_inq_dimlen(ncid, dimid, land_cat) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimlen: land_cat" !KWM endif !KWM !KWM iret = nf_inq_dimid(ncid, "soil_cat", dimid) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimid: soil_cat" !KWM endif !KWM !KWM iret = nf_inq_dimlen(ncid, dimid, soil_cat) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimlen: soil_cat" !KWM endif !KWM ! Various map attributes iret = nf_get_att_text(ncid, NF_GLOBAL, "map_projection", map%projection) if (iret /= 0) then stop "nf_get_att_float: projection" endif if (map%projection(1:17) == "LAMBERT CONFORMAL") map%projection = "LC" iret = nf_get_att_real(ncid, NF_GLOBAL, "TRUELAT1", map%truelat1) if (iret /= 0) then stop "nf_get_att_float: truelat1" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "TRUELAT2", map%truelat2) if (iret /= 0) then stop "nf_get_att_float: truelat2" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "DX", map%dx) if (iret /= 0) then stop "nf_get_att_float: dx" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "DY", map%dy) if (iret /= 0) then stop "nf_get_att_float: dy" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "STAND_LON", map%xlonc) if (iret /= 0) then stop "nf_get_att_float: xlonc" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "corner_lats", corner_lat) if (iret /= 0) then stop "nf_get_att_float: corner_lat" endif iret = nf_get_att_real(ncid, NF_GLOBAL, "corner_lons", corner_lon) if (iret /= 0) then stop "nf_get_att_float: corner_lon" endif map%reflat = corner_lat(1) map%reflon = corner_lon(1) map%refx = 1.0 map%refy = 1.0 iret = nf_close(ncid) end subroutine rdwrfstatic subroutine rdhrldas(stn, stnlat, stnlon, nstation, idate, sm) use module_ver implicit none integer, intent(in) :: nstation character(len=*), dimension(nstation), intent(in) :: stn real, dimension(nstation), intent(in) :: stnlat real, dimension(nstation), intent(in) :: stnlon character(len=*), intent(in) :: idate real, dimension(4,nstation), intent(out) :: sm character(len=256) :: flnm #include integer :: iret, ncid integer :: dimid, varid, idim, jdim real, allocatable, dimension(:,:) :: var real :: x, y, i, j integer :: ista, ilvl sm = -99999. write(flnm, & '("/scholar2/kmanning/VERY_TEMPORARY_HRLDAS/Run_3.2/",A10,".LDASOUT_DOMAIN1")') idate(1:10) iret = nf_open(flnm, NF_WRITE, ncid) if (iret /= 0) then write(*,'("Problem opening flnm: ''", A, "''")') & trim(flnm) stop endif !KWM iret = nf_inq_dimid(ncid, "idim", dimid) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimid: idim" !KWM endif !KWM !KWM iret = nf_inq_dimlen(ncid, dimid, idim) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimlen: idim" !KWM endif !KWM !KWM iret = nf_inq_dimid(ncid, "jdim", dimid) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimid: jdim" !KWM endif !KWM !KWM iret = nf_inq_dimlen(ncid, dimid, jdim) !KWM if (iret /= 0) then !KWM stop "nf_inq_dimlen: jdim" !KWM endif !KWM print*, 'ix, jx = ', map%ix, map%jx !KWM print*, 'map = ', map ! Read a soil-moisture field do ilvl = 1, 4 iret = nf_inq_varid(ncid, "SOIL_M_"//char(ilvl+ichar("0")), varid) if (iret /= 0) then stop "nf_inq_varid: SOIL_M_" endif !KWM print*, "SOIL_M_"//char(ilvl+ichar("0")), ' varid = ', varid allocate(var(map%ix, map%jx)) iret = nf_get_var_real(ncid, varid, var) if (iret /= 0) then stop "nf_get_var_real: SOIL_M_1" endif do ista = 1, nstation call lltoxy_generic(stnlat(ista), stnlon(ista), x, y, & map%projection, 1.E-3*map%dx, map%reflat, map%reflon, map%refx, map%refy,& map%xlonc, map%truelat1, map%truelat2) ! print*, 'lat, lon, x, y = ', stnlat(ista), stnlon(ista), x, y sm(ilvl, ista) = var(nint(x), nint(y)) enddo deallocate(var) enddo end subroutine rdhrldas