hwrf-baseline-20120103-1354/util/src/plotfmt.F

   1 program plotfmt
   2
   3   use read_met_module
   4
   5   implicit none
   6 !
   7 ! Utility program to plot up files created by pregrid / SI / ungrib.
   8 ! Uses NCAR graphics routines.  If you don't have NCAR Graphics, you're
   9 ! out of luck.
  10 !
  11    INTEGER :: istatus
  12    integer :: idum, ilev
  13
  14    CHARACTER ( LEN =132 )            :: flnm
  15
  16    TYPE (met_data)                   :: fg_data
  17
  18 !
  19 !   Set up the graceful stop (Sun, SGI, DEC).
  20 !
  21    integer, external :: graceful_stop
  22 #if (defined(_DOUBLEUNDERSCORE) && defined(MACOS)) || defined(NO_SIGNAL)
  23    ! we do not do any signaling
  24 #else
  25    integer, external :: signal
  26 #endif
  27    integer :: iii
  28
  29 #if (defined(_DOUBLEUNDERSCORE) && defined(MACOS)) || defined(NO_SIGNAL)
  30   ! still more no signaling
  31 #else
  32   iii = signal(2, graceful_stop, -1)
  33 #endif
  34
  35   call getarg(1,flnm)
  36
  37    IF ( flnm(1:1) == ' ' ) THEN
  38       print *,'USAGE: plotfmt.exe <filename>'
  39       print *,'       where <filename> is the name of an intermediate-format file'
  40       STOP
  41    END IF
  42
  43   call gopks(6,idum)
  44   call gopwk(1,55,1)
  45   call gopwk(2,56,3)
  46   call gacwk(1)
  47   call gacwk(2)
  48   call pcseti('FN', 21)
  49   call pcsetc('FC', '~')
  50
  51   call gscr(1,0, 1.000, 1.000, 1.000)
  52   call gscr(1,1, 0.000, 0.000, 0.000)
  53   call gscr(1,2, 0.900, 0.600, 0.600)
  54
  55    CALL read_met_init(TRIM(flnm), .true., '0000-00-00_00', istatus)
  56
  57    IF ( istatus == 0 ) THEN
  58
  59       CALL  read_next_met_field(fg_data, istatus)
  60
  61       DO WHILE (istatus == 0)
  62
  63          ilev = nint(fg_data%xlvl)
  64
  65          if (fg_data%iproj == PROJ_LATLON) then
  66             call plt2d(fg_data%slab, fg_data%nx, fg_data%ny, fg_data%iproj, &
  67                        fg_data%startlat, fg_data%startlon, fg_data%deltalon, &
  68                        fg_data%deltalat, fg_data%xlonc, fg_data%truelat1, fg_data%truelat2, &
  69                        fg_data%field, ilev, fg_data%units, fg_data%version, fg_data%desc, &
  70                        fg_data%map_source, TRIM(flnm))
  71          else
  72             call plt2d(fg_data%slab, fg_data%nx, fg_data%ny, fg_data%iproj, &
  73                        fg_data%startlat, fg_data%startlon, fg_data%dx, fg_data%dy, fg_data%xlonc, &
  74                        fg_data%truelat1, fg_data%truelat2, fg_data%field, ilev, fg_data%units, &
  75                        fg_data%version, fg_data%desc, fg_data%map_source, TRIM(flnm))
  76          end if
  77
  78
  79          IF (ASSOCIATED(fg_data%slab)) DEALLOCATE(fg_data%slab)
  80
  81          CALL  read_next_met_field(fg_data, istatus)
  82       END DO
  83
  84       CALL read_met_close()
  85
  86    ELSE
  87
  88       print *, 'File = ',TRIM(flnm)
  89       print *, 'Problem with input file, I can''t open it'
  90       STOP
  91
  92    END IF
  93
  94   call stopit
  95
  96 end program plotfmt
  97
  98 subroutine plt2d(tcr2d, iz, jz, llflag, &
  99      lat1, lon1, dx, dy, lov, truelat1, truelat2, &
 100      field, ilev, units, ifv, Desc, source, flnm)
 101
 102   use misc_definitions_module
 103
 104   implicit none
 105
 106   integer :: llflag
 107   integer :: iz, jz, ifv
 108   real, dimension(iz,jz) :: tcr2d(iz,jz)
 109   real :: lat1, lon1, lov, truelat1, truelat2
 110   real :: dx, dy
 111   character(len=*) :: field
 112   character(len=*) ::  units
 113   character(len=*) :: Desc
 114   character(len=*) :: source
 115   character(len=30) :: hunit
 116   character(len=32) :: tmp32
 117   character (len=*) :: flnm
 118
 119   integer :: iproj, ierr
 120   real :: pl1, pl2, pl3, pl4, plon, plat, rota, phic
 121   real :: xl, xr, xb, xt, wl, wr, wb, wt, yb
 122   integer :: ml, ih, i, j
 123
 124   integer, parameter :: lwrk = 20000, liwk = 50000
 125   real, dimension(lwrk) :: rwrk
 126   integer, dimension(liwk) :: iwrk
 127
 128   integer :: ilev
 129   integer :: found_it
 130   character(len=8) :: hlev
 131
 132 ! declarations for windowing
 133   integer :: ioff, joff, i1, j1, ix, jx, funit
 134   real, allocatable,dimension(:,:)  :: scr2d
 135   logical :: is_used
 136   namelist /plotfmt/ ix, jx, ioff, joff
 137
 138 ! This version allows the plotting of subsets of a lat/lon grid (i.e. NCEP GFS).
 139 ! ix,jx are the dimensions of the subset. ioff,joff are the offsets from 1,1
 140
 141   ix = iz
 142   jx = jz
 143   ioff = 0
 144   joff = 0
 145
 146 ! Read parameters from Fortran namelist
 147   do funit=10,100
 148     inquire(unit=funit, opened=is_used)
 149     if (.not. is_used) exit
 150   end do
 151   open(funit,file='namelist.wps',status='old',form='formatted',err=1000)
 152   read(funit,plotfmt,iostat=found_it)
 153   close(funit)
 154   if(found_it .gt. 0 ) then
 155      print *,'error reading the plotfmt namelist record in namelist.wps'
 156      print *,'you may have: ix, jx, ioff, joff ONLY'
 157      stop 1234
 158   end if
 159
 160
 161 ! ioff = 250     ! e.g. east of the Philippines from 0.5 degree GFS
 162 ! joff = 140
 163 ! ix = 20
 164 ! jx = 20
 165
 166   if (ix+ioff .gt. iz .or. jx+joff .gt. jz) then
 167 !   print *,'map subset is too large. Setting to full domain'
 168     ix = iz
 169     jx = jz
 170     ioff = 0
 171     joff = 0
 172   endif
 173 ! compute upper left point for the map   (works for NCEP GFS and godas)
 174   pl1 = lat1 + (joff*dy)
 175   pl2 = lon1 + (ioff*dx)
 176
 177   allocate (scr2d(ix,jx))
 178
 179   do i = 1, ix
 180   do j = 1, jx
 181     i1 = i + ioff
 182     j1 = j + joff
 183     scr2d(i,j) = tcr2d(i1,j1)
 184   enddo
 185   enddo
 186
 187   select case (llflag)
 188   case (PROJ_LATLON)
 189      call fmtxyll(float(ix), float(jx), pl3, pl4, 'CE', pl1, pl2, &
 190           plon, truelat1, truelat2, dx, dy)
 191      plon = (pl2 + pl4) / 2.
 192      plat = 0.
 193      rota = 0.
 194      iproj=8
 195   case (PROJ_MERC)
 196      pl1 = lat1
 197      pl2 = lon1
 198      plon = 0.
 199      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ME', pl1, pl2, &
 200           plon, truelat1, truelat2, dx, dy)
 201      plat = 0.
 202      rota = 0
 203      iproj = 9
 204   case (PROJ_LC)
 205      pl1 = lat1
 206      pl2 = lon1
 207      plon = lov
 208      call fmtxyll(float(ix), float(jx), pl3, pl4, 'LC', pl1, pl2,&
 209           plon, truelat1, truelat2, dx, dy)
 210      plat = truelat1
 211      rota = truelat2
 212      iproj=3
 213 ! This never used to be a problem, but currently we seem to need
 214 ! truelat1 (in plat) differ from truelat2 (in rota) for the
 215 ! NCAR-Graphics map routines to work.  Maybe it's just a compiler
 216 ! thing.  So if the truelats are the same, we add an epsilon:
 217      if (abs(plat - rota) < 1.E-8) then
 218         plat = plat + 1.E-8
 219         rota = rota - 1.E-8
 220      endif
 221   case (PROJ_PS)
 222      print*, 'ix, jx = ', ix, jx
 223      print*, 'lat1, lon1 = ', lat1, lon1
 224      pl1 = lat1
 225      pl2 = lon1
 226      plon = lov
 227      plat = 90.
 228      print*, 'plon, plat = ', plon, plat
 229      phic = 90.
 230      rota = 0.
 231      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ST', pl1, pl2,&
 232           plon, truelat1, truelat2, dx, dy)
 233      iproj=1
 234      print*, pl1, pl2, pl3, pl4
 235   case default
 236      print*,'Unsupported map projection ',llflag,' in input'
 237      stop
 238   end select
 239
 240   call gsplci(2)   ! Use a different color for the map
 241   call supmap(iproj,plat,plon,rota,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 242   call gsplci(1)
 243 ! call supmap(iproj,plat+0.001,plon,rota-0.001,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 244   if (ierr.ne.0) then
 245      print*, 'supmap ierr = ', ierr
 246          stop
 247 !    stop
 248   endif
 249   call getset(xl,xr,xb,xt,wl,wr,wb,wt,ml)
 250
 251   write(hlev,'(I8)') ilev
 252
 253   call set(0., 1., 0., 1., 0., 1., 0., 1., 1)
 254   if ( xb .lt. .16 ) then
 255     yb = .16    ! xb depends on the projection, so fix yb and use it for labels
 256   else
 257     yb = xb
 258   endif
 259   call pchiqu(0.1, yb-0.05, hlev//'  '//field, .020, 0.0, -1.0)
 260   print*, field//'#'//units//'#'//trim(Desc)
 261 ! call pchiqu(0.1, xb-0.12, Desc, .012, 0.0, -1.0)
 262   hunit = '                                      '
 263   ih = 0
 264   do i = 1, len(units)
 265      if (units(i:i).eq.'{') then
 266         hunit(ih+1:ih+3) = '~S~'
 267         ih = ih + 3
 268         elseif (units(i:i).eq.'}') then
 269         hunit(ih+1:ih+3) = '~N~'
 270         ih = ih + 3
 271      else
 272         ih = ih + 1
 273         hunit(ih:ih) = units(i:i)
 274      endif
 275   enddo
 276   if ( ifv .le. 3 ) then
 277     tmp32 = 'MM5 intermediate format'
 278   else if ( ifv .eq. 4 ) then
 279     tmp32 = 'SI intermediate format'
 280   else if ( ifv .eq. 5 ) then
 281     tmp32 = 'WPS intermediate format'
 282   endif
 283   call pchiqu(0.1, yb-0.09, hunit, .015, 0.0, -1.0)
 284   call pchiqu(0.1, yb-0.12, Desc, .013, 0.0, -1.0)
 285   call pchiqu(0.6, yb-0.12, source, .013, 0.0, -1.0)
 286   call pchiqu(0.1, yb-0.15, tmp32, .013, 0.0, -1.0)
 287   call pchiqu(0.6, yb-0.15, flnm, .013, 0.0, -1.0)
 288
 289   call set(xl,xr,xb,xt,1.,float(ix),1.,float(jx),ml)
 290
 291   call CPSETI ('SET - Do-SET-Call Flag', 0)
 292   call CPSETR ('SPV - Special Value', -1.E30)
 293
 294   call cpseti('LLP', 3)
 295
 296   if (dy.lt.0.) then
 297      call array_flip(scr2d, ix, jx)
 298   endif
 299
 300   call cprect(scr2d,ix,ix,jx,rwrk,lwrk,iwrk,liwk)
 301   call cpcldr(scr2d,rwrk,iwrk)
 302   call cplbdr(scr2d,rwrk,iwrk)
 303
 304   deallocate (scr2d)
 305   call frame
 306   return
 307 1000 write(0,*) 'Error opening file namelist.wps, Stopping'
 308   stop 'namelist missing'
 309
 310 end subroutine plt2d
 311
 312 subroutine stopit
 313   call graceful_stop
 314 end
 315
 316 subroutine graceful_stop
 317   call gdawk(2)
 318   call gdawk(1)
 319   call gclwk(2)
 320   call gclwk(1)
 321   call gclks
 322   print*, 'Graceful Stop.'
 323   stop
 324 end subroutine graceful_stop
 325
 326 subroutine fmtxyll(x, y, xlat, xlon, project, glat1, glon1, gclon,&
 327      gtrue1, gtrue2, gdx, gdy)
 328   implicit none
 329
 330   real , intent(in) :: x, y, glat1, glon1, gtrue1, gtrue2, gdx, gdy, gclon
 331   character(len=2), intent(in) :: project
 332   real , intent(out) :: xlat, xlon
 333
 334   real :: gx1, gy1, gkappa
 335   real :: grrth = 6370.
 336
 337   real :: r, y1
 338   integer :: iscan, jscan
 339   real, parameter :: pi = 3.1415926534
 340   real, parameter :: degran = pi/180.
 341   real, parameter :: raddeg = 180./pi
 342   real :: gt
 343
 344   if (project.eq.'CE') then  ! Cylindrical Equidistant grid
 345
 346      xlat = glat1 + gdy*(y-1.)
 347      xlon = glon1 + gdx*(x-1.)
 348
 349   elseif (project == "ME") then
 350
 351      gt = grrth * cos(gtrue1*degran)
 352      xlon = glon1 + (gdx*(x-1.)/gt)*raddeg
 353      y1 = gt*alog((1.+sin(glat1*degran))/cos(glat1*degran))/gdy
 354      xlat = 90. - 2. * atan(exp(-gdy*(y+y1-1.)/gt))*raddeg
 355
 356   elseif (project.eq.'ST') then  ! Polar Stereographic grid
 357
 358      r = grrth/gdx * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
 359      gx1 = r * sind(glon1-gclon)
 360      gy1 = - r * cosd(glon1-gclon)
 361
 362      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 363      xlat = 90. - 2.*atan2d((r*gdx),(grrth*(1.+sind(gtrue1))))
 364      xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon
 365
 366   elseif (project.eq.'LC') then  ! Lambert-conformal grid
 367
 368      call glccone(gtrue1, gtrue2, 1, gkappa)
 369
 370      r = grrth/(gdx*gkappa)*sind(90.-gtrue1) * &
 371           (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
 372      gx1 =  r*sind(gkappa*(glon1-gclon))
 373      gy1 = -r*cosd(gkappa*(glon1-gclon))
 374
 375      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 376      xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
 377           ((r*gkappa*gdx)/(grrth*sind(90.-gtrue1)))**(1./gkappa))
 378      xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon
 379
 380   else
 381
 382      write(*,'("Unrecoginzed projection: ", A)') project
 383      write(*,'("Abort in FMTXYLL",/)')
 384      stop
 385
 386   endif
 387 contains
 388   real function sind(theta)
 389     real :: theta
 390     sind = sin(theta*degran)
 391   end function sind
 392   real function cosd(theta)
 393     real :: theta
 394     cosd = cos(theta*degran)
 395   end function cosd
 396   real function tand(theta)
 397     real :: theta
 398     tand = tan(theta*degran)
 399   end function tand
 400   real function atand(x)
 401     real :: x
 402     atand = atan(x)*raddeg
 403   end function atand
 404   real function atan2d(x,y)
 405     real :: x,y
 406     atan2d = atan2(x,y)*raddeg
 407   end function atan2d
 408
 409   subroutine glccone (fsplat,ssplat,sign,confac)
 410     implicit none
 411     real, intent(in) :: fsplat,ssplat
 412     integer, intent(in) :: sign
 413     real, intent(out) :: confac
 414     if (abs(fsplat-ssplat).lt.1.E-3) then
 415        confac = sind(fsplat)
 416     else
 417        confac = log10(cosd(fsplat))-log10(cosd(ssplat))
 418        confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
 419             log10(tand(45.-float(sign)*ssplat/2.)))
 420     endif
 421   end subroutine glccone
 422
 423 end subroutine fmtxyll
 424
 425 subroutine array_flip(array, ix, jx)
 426   implicit none
 427   integer :: ix, jx
 428   real , dimension(ix,jx) :: array
 429
 430   real, dimension(ix) :: hold
 431   integer :: i, j, jj
 432
 433   do j = 1, jx/2
 434      jj = jx+1-j
 435      hold = array(1:ix, j)
 436      array(1:ix,j) = array(1:ix,jj)
 437      array(1:ix,jj) = hold
 438   enddo
 439 end subroutine array_flip
 440