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, lexist
 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   inquire(file='namelist.wps', exist = lexist)
 152   if (lexist) then
 153     open(funit,file='namelist.wps',status='old',form='formatted',err=1000)
 154     read(funit,plotfmt,iostat=found_it)
 155     close(funit)
 156     if(found_it .gt. 0 ) then
 157        print *,'error reading the plotfmt namelist record in namelist.wps'
 158        print *,'only ix, jx, ioff, joff are recognized.'
 159        stop 'namelist problem'
 160     end if
 161   end if
 162
 163
 164 ! ioff = 250     ! e.g. east of the Philippines from 0.5 degree GFS
 165 ! joff = 140
 166 ! ix = 20
 167 ! jx = 20
 168
 169   if (ix+ioff .gt. iz .or. jx+joff .gt. jz) then
 170 !   print *,'map subset is too large. Setting to full domain'
 171     ix = iz
 172     jx = jz
 173     ioff = 0
 174     joff = 0
 175   endif
 176 ! compute upper left point for the map   (works for NCEP GFS and godas)
 177   pl1 = lat1 + (joff*dy)
 178   pl2 = lon1 + (ioff*dx)
 179
 180   allocate (scr2d(ix,jx))
 181
 182   do i = 1, ix
 183   do j = 1, jx
 184     i1 = i + ioff
 185     j1 = j + joff
 186     scr2d(i,j) = tcr2d(i1,j1)
 187   enddo
 188   enddo
 189
 190   select case (llflag)
 191   case (PROJ_LATLON)
 192      call fmtxyll(float(ix), float(jx), pl3, pl4, 'CE', pl1, pl2, &
 193           plon, truelat1, truelat2, dx, dy)
 194      plon = (pl2 + pl4) / 2.
 195      plat = 0.
 196      rota = 0.
 197      iproj=8
 198   case (PROJ_MERC)
 199      pl1 = lat1
 200      pl2 = lon1
 201      plon = 0.
 202      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ME', pl1, pl2, &
 203           plon, truelat1, truelat2, dx, dy)
 204      plat = 0.
 205      rota = 0
 206      iproj = 9
 207   case (PROJ_LC)
 208      pl1 = lat1
 209      pl2 = lon1
 210      plon = lov
 211      call fmtxyll(float(ix), float(jx), pl3, pl4, 'LC', pl1, pl2,&
 212           plon, truelat1, truelat2, dx, dy)
 213      plat = truelat1
 214      rota = truelat2
 215      iproj=3
 216 ! This never used to be a problem, but currently we seem to need
 217 ! truelat1 (in plat) differ from truelat2 (in rota) for the
 218 ! NCAR-Graphics map routines to work.  Maybe it's just a compiler
 219 ! thing.  So if the truelats are the same, we add an epsilon:
 220      if (abs(plat - rota) < 1.E-8) then
 221         plat = plat + 1.E-8
 222         rota = rota - 1.E-8
 223      endif
 224   case (PROJ_PS)
 225      print*, 'ix, jx = ', ix, jx
 226      print*, 'lat1, lon1 = ', lat1, lon1
 227      pl1 = lat1
 228      pl2 = lon1
 229      plon = lov
 230      if (truelat1 .lt. 0. ) then
 231        plat = -90.
 232      else
 233        plat = 90.
 234      endif
 235      print*, 'plon, plat = ', plon, plat
 236      rota = 0.
 237      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ST', pl1, pl2,&
 238           plon, truelat1, truelat2, dx, dy)
 239      iproj=1
 240      print*, pl1, pl2, pl3, pl4
 241   case default
 242      print*,'Unsupported map projection ',llflag,' in input'
 243      stop
 244   end select
 245
 246   call gsplci(2)   ! Use a different color for the map
 247 ! jlts = 2 means corner points are provided in p1,p2,pl3,pl4
 248   call supmap(iproj,plat,plon,rota,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 249   call gsplci(1)
 250 ! call supmap(iproj,plat+0.001,plon,rota-0.001,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 251   if (ierr.ne.0) then
 252      print*, 'supmap ierr = ', ierr
 253          stop
 254 !    stop
 255   endif
 256   call getset(xl,xr,xb,xt,wl,wr,wb,wt,ml)
 257
 258   write(hlev,'(I8)') ilev
 259
 260   call set(0., 1., 0., 1., 0., 1., 0., 1., 1)
 261   if ( xb .lt. .16 ) then
 262     yb = .16    ! xb depends on the projection, so fix yb and use it for labels
 263   else
 264     yb = xb
 265   endif
 266   call pchiqu(0.1, yb-0.05, hlev//'  '//field, .020, 0.0, -1.0)
 267   print*, field//'#'//units//'#'//trim(Desc)
 268 ! call pchiqu(0.1, xb-0.12, Desc, .012, 0.0, -1.0)
 269   hunit = '                                      '
 270   ih = 0
 271   do i = 1, len(units)
 272      if (units(i:i).eq.'{') then
 273         hunit(ih+1:ih+3) = '~S~'
 274         ih = ih + 3
 275         elseif (units(i:i).eq.'}') then
 276         hunit(ih+1:ih+3) = '~N~'
 277         ih = ih + 3
 278      else
 279         ih = ih + 1
 280         hunit(ih:ih) = units(i:i)
 281      endif
 282   enddo
 283   if ( ifv .le. 3 ) then
 284     tmp32 = 'MM5 intermediate format'
 285   else if ( ifv .eq. 4 ) then
 286     tmp32 = 'SI intermediate format'
 287   else if ( ifv .eq. 5 ) then
 288     tmp32 = 'WPS intermediate format'
 289   endif
 290   call pchiqu(0.1, yb-0.09, hunit, .015, 0.0, -1.0)
 291   call pchiqu(0.1, yb-0.12, Desc, .013, 0.0, -1.0)
 292   call pchiqu(0.6, yb-0.12, source, .013, 0.0, -1.0)
 293   call pchiqu(0.1, yb-0.15, tmp32, .013, 0.0, -1.0)
 294   call pchiqu(0.6, yb-0.15, flnm, .013, 0.0, -1.0)
 295
 296   call set(xl,xr,xb,xt,1.,float(ix),1.,float(jx),ml)
 297
 298   call CPSETI ('SET - Do-SET-Call Flag', 0)
 299   call CPSETR ('SPV - Special Value', -1.E30)
 300
 301   call cpseti('LLP', 3)
 302
 303   if (dy.lt.0.) then
 304      call array_flip(scr2d, ix, jx)
 305   endif
 306
 307   call cprect(scr2d,ix,ix,jx,rwrk,lwrk,iwrk,liwk)
 308   call cpcldr(scr2d,rwrk,iwrk)
 309   call cplbdr(scr2d,rwrk,iwrk)
 310
 311   deallocate (scr2d)
 312   call frame
 313   return
 314 1000 write(0,*) 'Error opening file namelist.wps, Stopping'
 315   stop 'namelist missing'
 316
 317 end subroutine plt2d
 318
 319 subroutine stopit
 320   call graceful_stop
 321 end
 322
 323 subroutine graceful_stop
 324   call gdawk(2)
 325   call gdawk(1)
 326   call gclwk(2)
 327   call gclwk(1)
 328   call gclks
 329   print*, 'Graceful Stop.'
 330   stop
 331 end subroutine graceful_stop
 332
 333 subroutine fmtxyll(x, y, xlat, xlon, project, glat1, glon1, gclon,&
 334      gtrue1, gtrue2, gdx, gdy)
 335   implicit none
 336
 337   real , intent(in) :: x, y, glat1, glon1, gtrue1, gtrue2, gdx, gdy, gclon
 338   character(len=2), intent(in) :: project
 339   real , intent(out) :: xlat, xlon
 340
 341   real :: gx1, gy1, gkappa
 342   real :: grrth = 6370., phist, phemi
 343
 344   real :: r, y1
 345   integer :: iscan, jscan
 346   real, parameter :: pi = 3.1415926534
 347   real, parameter :: degran = pi/180.
 348   real, parameter :: raddeg = 180./pi
 349   real :: gt
 350
 351   if (project.eq.'CE') then  ! Cylindrical Equidistant grid
 352
 353      xlat = glat1 + gdy*(y-1.)
 354      xlon = glon1 + gdx*(x-1.)
 355
 356   elseif (project == "ME") then
 357
 358      gt = grrth * cos(gtrue1*degran)
 359      xlon = glon1 + (gdx*(x-1.)/gt)*raddeg
 360      y1 = gt*alog((1.+sin(glat1*degran))/cos(glat1*degran))/gdy
 361      xlat = 90. - 2. * atan(exp(-gdy*(y+y1-1.)/gt))*raddeg
 362
 363   elseif (project.eq.'ST') then  ! Polar Stereographic grid
 364
 365      if (gtrue1 .lt. 0.)then
 366        phemi = -1.
 367      else
 368        phemi = 1.
 369      endif
 370      r = grrth/gdx * tand((phemi*90. - glat1)/2.) * (1.+ phemi * sind(gtrue1))
 371      gx1 = r * sind(glon1-gclon)
 372      gy1 = - r * cosd(glon1-gclon)
 373
 374      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 375      xlat = phemi*90. - 2.*atan2d((r*gdx),(grrth*(1.+ phemi*sind(gtrue1))))
 376      xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon
 377
 378   elseif (project.eq.'LC') then  ! Lambert-conformal grid
 379
 380      call glccone(gtrue1, gtrue2, 1, gkappa)
 381
 382      r = grrth/(gdx*gkappa)*sind(90.-gtrue1) * &
 383           (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
 384      gx1 =  r*sind(gkappa*(glon1-gclon))
 385      gy1 = -r*cosd(gkappa*(glon1-gclon))
 386
 387      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 388      xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
 389           ((r*gkappa*gdx)/(grrth*sind(90.-gtrue1)))**(1./gkappa))
 390      xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon
 391
 392   else
 393
 394      write(*,'("Unrecoginzed projection: ", A)') project
 395      write(*,'("Abort in FMTXYLL",/)')
 396      stop
 397
 398   endif
 399 contains
 400   real function sind(theta)
 401     real :: theta
 402     sind = sin(theta*degran)
 403   end function sind
 404   real function cosd(theta)
 405     real :: theta
 406     cosd = cos(theta*degran)
 407   end function cosd
 408   real function tand(theta)
 409     real :: theta
 410     tand = tan(theta*degran)
 411   end function tand
 412   real function atand(x)
 413     real :: x
 414     atand = atan(x)*raddeg
 415   end function atand
 416   real function atan2d(x,y)
 417     real :: x,y
 418     atan2d = atan2(x,y)*raddeg
 419   end function atan2d
 420
 421   subroutine glccone (fsplat,ssplat,sign,confac)
 422     implicit none
 423     real, intent(in) :: fsplat,ssplat
 424     integer, intent(in) :: sign
 425     real, intent(out) :: confac
 426     if (abs(fsplat-ssplat).lt.1.E-3) then
 427        confac = sind(fsplat)
 428     else
 429        confac = log10(cosd(fsplat))-log10(cosd(ssplat))
 430        confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
 431             log10(tand(45.-float(sign)*ssplat/2.)))
 432     endif
 433   end subroutine glccone
 434
 435 end subroutine fmtxyll
 436
 437 subroutine array_flip(array, ix, jx)
 438   implicit none
 439   integer :: ix, jx
 440   real , dimension(ix,jx) :: array
 441
 442   real, dimension(ix) :: hold
 443   integer :: i, j, jj
 444
 445   do j = 1, jx/2
 446      jj = jx+1-j
 447      hold = array(1:ix, j)
 448      array(1:ix,j) = array(1:ix,jj)
 449      array(1:ix,jj) = hold
 450   enddo
 451 end subroutine array_flip
 452