WPS/util/src/plotfmt.F90

   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)
  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)
  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, fg_data%map_source)
  70          else
  71             call plt2d(fg_data%slab, fg_data%nx, fg_data%ny, fg_data%iproj, &
  72                        fg_data%startlat, fg_data%startlon, fg_data%dx, fg_data%dy, fg_data%xlonc, &
  73                        fg_data%truelat1, fg_data%truelat2, fg_data%field, ilev, fg_data%units, &
  74                        fg_data%version, fg_data%desc, fg_data%map_source)
  75          end if
  76
  77
  78          IF (ASSOCIATED(fg_data%slab)) DEALLOCATE(fg_data%slab)
  79
  80          CALL  read_next_met_field(fg_data, istatus)
  81       END DO
  82
  83       CALL read_met_close()
  84
  85    ELSE
  86
  87       print *, 'File = ',TRIM(flnm)
  88       print *, 'Problem with input file, I can''t open it'
  89       STOP
  90
  91    END IF
  92
  93   call stopit
  94
  95 end program plotfmt
  96
  97 subroutine plt2d(scr2d, ix, jx, llflag, &
  98      lat1, lon1, dx, dy, lov, truelat1, truelat2, &
  99      field, ilev, units, ifv, Desc, source)
 100
 101   use misc_definitions_module
 102
 103   implicit none
 104
 105   integer :: llflag
 106   integer :: ix, jx, ifv
 107   real, dimension(ix,jx) :: scr2d(ix,jx)
 108   real :: lat1, lon1, lov, truelat1, truelat2
 109   real :: dx, dy
 110   character(len=*) :: field
 111   character(len=*) ::  units
 112   character(len=*) :: Desc
 113   character(len=*) :: source
 114   character(len=30) :: hunit
 115   character(len=32) :: tmp32
 116
 117   integer :: iproj, ierr
 118   real :: pl1, pl2, pl3, pl4, plon, plat, rota, phic
 119   real :: xl, xr, xb, xt, wl, wr, wb, wt, yb
 120   integer :: ml, ih, i
 121
 122   integer, parameter :: lwrk = 20000, liwk = 50000
 123   real, dimension(lwrk) :: rwrk
 124   integer, dimension(liwk) :: iwrk
 125
 126   integer :: ilev
 127   character(len=8) :: hlev
 128
 129   select case (llflag)
 130   case (PROJ_LATLON)
 131      pl1 = lat1
 132      pl2 = lon1
 133      call fmtxyll(float(ix), float(jx), pl3, pl4, 'CE', pl1, pl2, &
 134           plon, truelat1, truelat2, dx, dy)
 135      plon = (pl2 + pl4) / 2.
 136      plat = 0.
 137      rota = 0.
 138      iproj=8
 139   case (PROJ_MERC)
 140      pl1 = lat1
 141      pl2 = lon1
 142      plon = 0.
 143      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ME', pl1, pl2, &
 144           plon, truelat1, truelat2, dx, dy)
 145      plat = 0.
 146      rota = 0
 147      iproj = 9
 148   case (PROJ_LC)
 149      pl1 = lat1
 150      pl2 = lon1
 151      plon = lov
 152      call fmtxyll(float(ix), float(jx), pl3, pl4, 'LC', pl1, pl2,&
 153           plon, truelat1, truelat2, dx, dy)
 154      plat = truelat1
 155      rota = truelat2
 156      iproj=3
 157 ! This never used to be a problem, but currently we seem to need
 158 ! truelat1 (in plat) differ from truelat2 (in rota) for the
 159 ! NCAR-Graphics map routines to work.  Maybe it's just a compiler
 160 ! thing.  So if the truelats are the same, we add an epsilon:
 161      if (abs(plat - rota) < 1.E-8) then
 162         plat = plat + 1.E-8
 163         rota = rota - 1.E-8
 164      endif
 165   case (PROJ_PS)
 166      print*, 'ix, jx = ', ix, jx
 167      print*, 'lat1, lon1 = ', lat1, lon1
 168      pl1 = lat1
 169      pl2 = lon1
 170      plon = lov
 171      plat = 90.
 172      print*, 'plon, plat = ', plon, plat
 173      phic = 90.
 174      rota = 0.
 175      call fmtxyll(float(ix), float(jx), pl3, pl4, 'ST', pl1, pl2,&
 176           plon, truelat1, truelat2, dx, dy)
 177      iproj=1
 178      print*, pl1, pl2, pl3, pl4
 179   case default
 180      print*,'Unsupported map projection ',llflag,' in input'
 181      stop
 182   end select
 183
 184   call supmap(iproj,plat,plon,rota,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 185 ! call supmap(iproj,plat+0.001,plon,rota-0.001,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
 186   if (ierr.ne.0) then
 187      print*, 'supmap ierr = ', ierr
 188          stop
 189 !    stop
 190   endif
 191   call getset(xl,xr,xb,xt,wl,wr,wb,wt,ml)
 192
 193   write(hlev,'(I8)') ilev
 194
 195   call set(0., 1., 0., 1., 0., 1., 0., 1., 1)
 196   if ( xb .lt. .16 ) then
 197     yb = .16    ! xb depends on the projection, so fix yb and use it for labels
 198   else
 199     yb = xb
 200   endif
 201   call pchiqu(0.1, yb-0.05, hlev//'  '//field, .020, 0.0, -1.0)
 202   print*, field//'#'//units//'#'//trim(Desc)
 203 ! call pchiqu(0.1, xb-0.12, Desc, .012, 0.0, -1.0)
 204   hunit = '                                      '
 205   ih = 0
 206   do i = 1, len(units)
 207      if (units(i:i).eq.'{') then
 208         hunit(ih+1:ih+3) = '~S~'
 209         ih = ih + 3
 210         elseif (units(i:i).eq.'}') then
 211         hunit(ih+1:ih+3) = '~N~'
 212         ih = ih + 3
 213      else
 214         ih = ih + 1
 215         hunit(ih:ih) = units(i:i)
 216      endif
 217   enddo
 218   if ( ifv .le. 3 ) then
 219     tmp32 = 'MM5 intermediate format'
 220   else if ( ifv .eq. 4 ) then
 221     tmp32 = 'SI intermediate format'
 222   else if ( ifv .eq. 5 ) then
 223     tmp32 = 'WPS intermediate format'
 224   endif
 225   call pchiqu(0.1, yb-0.09, hunit, .015, 0.0, -1.0)
 226   call pchiqu(0.1, yb-0.12, Desc, .013, 0.0, -1.0)
 227   call pchiqu(0.6, yb-0.12, source, .013, 0.0, -1.0)
 228   call pchiqu(0.1, yb-0.15, tmp32, .013, 0.0, -1.0)
 229
 230   call set(xl,xr,xb,xt,1.,float(ix),1.,float(jx),ml)
 231
 232   call CPSETI ('SET - Do-SET-Call Flag', 0)
 233   call CPSETR ('SPV - Special Value', -1.E30)
 234
 235   if (dy.lt.0.) then
 236      call array_flip(scr2d, ix, jx)
 237   endif
 238
 239   call cpseti('LLP', 3)
 240   call cprect(scr2d,ix,ix,jx,rwrk,lwrk,iwrk,liwk)
 241   call cpcldr(scr2d,rwrk,iwrk)
 242   call cplbdr(scr2d,rwrk,iwrk)
 243
 244   if (dy.lt.0.) then
 245      call array_flip(scr2d, ix, jx)
 246   endif
 247
 248   call frame
 249
 250 end subroutine plt2d
 251
 252 subroutine stopit
 253   call graceful_stop
 254 end
 255
 256 subroutine graceful_stop
 257   call gdawk(2)
 258   call gdawk(1)
 259   call gclwk(2)
 260   call gclwk(1)
 261   call gclks
 262   print*, 'Graceful Stop.'
 263   stop
 264 end subroutine graceful_stop
 265
 266 subroutine fmtxyll(x, y, xlat, xlon, project, glat1, glon1, gclon,&
 267      gtrue1, gtrue2, gdx, gdy)
 268   implicit none
 269
 270   real , intent(in) :: x, y, glat1, glon1, gtrue1, gtrue2, gdx, gdy, gclon
 271   character(len=2), intent(in) :: project
 272   real , intent(out) :: xlat, xlon
 273
 274   real :: gx1, gy1, gkappa
 275   real :: grrth = 6370.
 276
 277   real :: r, y1
 278   integer :: iscan, jscan
 279   real, parameter :: pi = 3.1415926534
 280   real, parameter :: degran = pi/180.
 281   real, parameter :: raddeg = 180./pi
 282   real :: gt
 283
 284   if (project.eq.'CE') then  ! Cylindrical Equidistant grid
 285
 286      xlat = glat1 + gdy*(y-1.)
 287      xlon = glon1 + gdx*(x-1.)
 288
 289   elseif (project == "ME") then
 290
 291      gt = grrth * cos(gtrue1*degran)
 292      xlon = glon1 + (gdx*(x-1.)/gt)*raddeg
 293      y1 = gt*alog((1.+sin(glat1*degran))/cos(glat1*degran))/gdy
 294      xlat = 90. - 2. * atan(exp(-gdy*(y+y1-1.)/gt))*raddeg
 295
 296   elseif (project.eq.'ST') then  ! Polar Stereographic grid
 297
 298      r = grrth/gdx * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
 299      gx1 = r * sind(glon1-gclon)
 300      gy1 = - r * cosd(glon1-gclon)
 301
 302      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 303      xlat = 90. - 2.*atan2d((r*gdx),(grrth*(1.+sind(gtrue1))))
 304      xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon
 305
 306   elseif (project.eq.'LC') then  ! Lambert-conformal grid
 307
 308      call glccone(gtrue1, gtrue2, 1, gkappa)
 309
 310      r = grrth/(gdx*gkappa)*sind(90.-gtrue1) * &
 311           (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
 312      gx1 =  r*sind(gkappa*(glon1-gclon))
 313      gy1 = -r*cosd(gkappa*(glon1-gclon))
 314
 315      r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
 316      xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
 317           ((r*gkappa*gdx)/(grrth*sind(90.-gtrue1)))**(1./gkappa))
 318      xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon
 319
 320   else
 321
 322      write(*,'("Unrecoginzed projection: ", A)') project
 323      write(*,'("Abort in FMTXYLL",/)')
 324      stop
 325
 326   endif
 327 contains
 328   real function sind(theta)
 329     real :: theta
 330     sind = sin(theta*degran)
 331   end function sind
 332   real function cosd(theta)
 333     real :: theta
 334     cosd = cos(theta*degran)
 335   end function cosd
 336   real function tand(theta)
 337     real :: theta
 338     tand = tan(theta*degran)
 339   end function tand
 340   real function atand(x)
 341     real :: x
 342     atand = atan(x)*raddeg
 343   end function atand
 344   real function atan2d(x,y)
 345     real :: x,y
 346     atan2d = atan2(x,y)*raddeg
 347   end function atan2d
 348
 349   subroutine glccone (fsplat,ssplat,sign,confac)
 350     implicit none
 351     real, intent(in) :: fsplat,ssplat
 352     integer, intent(in) :: sign
 353     real, intent(out) :: confac
 354     if (abs(fsplat-ssplat).lt.1.E-3) then
 355        confac = sind(fsplat)
 356     else
 357        confac = log10(cosd(fsplat))-log10(cosd(ssplat))
 358        confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
 359             log10(tand(45.-float(sign)*ssplat/2.)))
 360     endif
 361   end subroutine glccone
 362
 363 end subroutine fmtxyll
 364
 365 subroutine array_flip(array, ix, jx)
 366   implicit none
 367   integer :: ix, jx
 368   real , dimension(ix,jx) :: array
 369
 370   real, dimension(ix) :: hold
 371   integer :: i, j, jj
 372
 373   do j = 1, jx/2
 374      jj = jx+1-j
 375      hold = array(1:ix, j)
 376      array(1:ix,j) = array(1:ix,jj)
 377      array(1:ix,jj) = hold
 378   enddo
 379 end subroutine array_flip
 380