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