source: trunk/MESOSCALE/LMD_MM_MARS/SRC/WPS/util/src/plotfmt.F90 @ 3568

program pltfmt

  use read_met_module

  implicit none
!
! Utility program to plot up the files created by pregrid.
! Uses NCAR graphics routines.  If you don't have NCAR Graphics, you're 
! out of luck.
!
   INTEGER :: istatus, version, nx, ny, iproj
   integer :: idum, ilev
   REAL :: xfcst, xlvl, startlat, startlon, starti, startj, &
           deltalat, deltalon, dx, dy, xlonc, truelat1, truelat2, earth_radius
   REAL, POINTER, DIMENSION(:,:) :: slab
   LOGICAL :: is_wind_grid_rel

   CHARACTER ( LEN =132 )            :: flnm

   CHARACTER ( LEN = 24 )            :: hdate
   CHARACTER ( LEN =  9 )            :: field
   CHARACTER ( LEN = 25 )            :: units
   CHARACTER ( LEN = 46 )            :: desc
   CHARACTER ( LEN = 32 )            :: map_source

!
!   Set up the graceful stop (Sun, SGI, DEC).
!
   integer, external :: graceful_stop
   integer, external :: signal
   integer :: iii

  iii = signal(2, graceful_stop, -1)

  call getarg(1,flnm)

  call gopks(6,idum)
  call gopwk(1,55,1)
  call gopwk(2,56,3)
  call gacwk(1)
  call gacwk(2)
  call pcseti('FN', 21)
  call pcsetc('FC', '~')

  call gscr(1,0, 1.000, 1.000, 1.000)
  call gscr(1,1, 0.000, 0.000, 0.000)
  call gscr(1,2, 0.900, 0.600, 0.600)

   CALL read_met_init(TRIM(flnm), .true., '0000-00-00_00', istatus)

   IF ( istatus == 0 ) THEN

      CALL  read_next_met_field(version, field, hdate, xfcst, xlvl, units, desc, &
                          iproj, startlat, startlon, starti, startj, deltalat, &
                          deltalon, dx, dy, xlonc, truelat1, truelat2, earth_radius, &
                          nx, ny, map_source, &
                          slab, is_wind_grid_rel, istatus)

      DO WHILE (istatus == 0)

         ilev = nint(xlvl)

         if (iproj == PROJ_LATLON) then
            call plt2d(slab, nx, ny, iproj, &
                       startlat, startlon, deltalon, deltalat, xlonc, truelat1, truelat2, &
                       field, ilev, units, version, desc, map_source)
         else
            call plt2d(slab, nx, ny, iproj, &
                       startlat, startlon, dx, dy, xlonc, truelat1, truelat2, &
                       field, ilev, units, version, desc, map_source)
         end if


         IF (ASSOCIATED(slab)) DEALLOCATE(slab)

         CALL  read_next_met_field(version, field, hdate, xfcst, xlvl, units, desc, &
                             iproj, startlat, startlon, starti, startj, deltalat, &
                             deltalon, dx, dy, xlonc, truelat1, truelat2, earth_radius, &
                             nx, ny, map_source, &
                             slab, is_wind_grid_rel, istatus)
      END DO

      CALL read_met_close()

   ELSE

      print *, 'File = ',TRIM(flnm)
      print *, 'Problem with input file, I can''t open it'
      STOP

   END IF

  call stopit

end program pltfmt

subroutine plt2d(scr2d, ix, jx, llflag, &
     lat1, lon1, dx, dy, lov, truelat1, truelat2, &
     field, ilev, units, ifv, Desc, source)

  use misc_definitions_module

  implicit none

  integer :: llflag
  integer :: ix, jx, ifv
  real, dimension(ix,jx) :: scr2d(ix,jx)
  real :: lat1, lon1, lov, truelat1, truelat2
  real :: dx, dy
  character(len=*) :: field
  character(len=*) ::  units
  character(len=*) :: Desc
  character(len=*) :: source
  character(len=30) :: hunit
  character(len=32) :: tmp32

  integer :: iproj, ierr
  real :: pl1, pl2, pl3, pl4, plon, plat, rota, phic
  real :: xl, xr, xb, xt, wl, wr, wb, wt
  integer :: ml, ih, i

  integer, parameter :: lwrk = 20000, liwk = 50000
  real, dimension(lwrk) :: rwrk
  integer, dimension(liwk) :: iwrk

  integer :: ilev
  character(len=8) :: hlev

  select case (llflag)
  case (PROJ_LATLON)
     pl1 = lat1
     pl2 = lon1
     call fmtxyll(float(ix), float(jx), pl3, pl4, 'CE', pl1, pl2, &
          plon, truelat1, truelat2, dx, dy)
     plon = (pl2 + pl4) / 2.
     plat = 0.
     rota = 0.
     iproj=8
  case (PROJ_MERC)
     pl1 = lat1
     pl2 = lon1
     plon = 0.
     call fmtxyll(float(ix), float(jx), pl3, pl4, 'ME', pl1, pl2, &
          plon, truelat1, truelat2, dx, dy)
     plat = 0.
     rota = 0
     iproj = 9
  case (PROJ_LC)
     pl1 = lat1
     pl2 = lon1
     plon = lov
     call fmtxyll(float(ix), float(jx), pl3, pl4, 'LC', pl1, pl2,&
          plon, truelat1, truelat2, dx, dy)
     plat = truelat1
     rota = truelat2
     iproj=3
! This never used to be a problem, but currently we seem to need
! truelat1 (in plat) differ from truelat2 (in rota) for the 
! NCAR-Graphics map routines to work.  Maybe it's just a compiler
! thing.  So if the truelats are the same, we add an epsilon:
     if (abs(plat - rota) < 1.E-8) then
        plat = plat + 1.E-8
        rota = rota - 1.E-8
     endif
  case (PROJ_PS)
     print*, 'ix, jx = ', ix, jx
     print*, 'lat1, lon1 = ', lat1, lon1
     pl1 = lat1
     pl2 = lon1
     plon = lov
     plat = 90.
     print*, 'plon, plat = ', plon, plat
     phic = 90.
     rota = 0.
     call fmtxyll(float(ix), float(jx), pl3, pl4, 'ST', pl1, pl2,&
          plon, truelat1, truelat2, dx, dy)
     iproj=1
     print*, pl1, pl2, pl3, pl4
  case default
     print*,'Unsupported map projection ',llflag,' in input'
     stop
  end select

  call supmap(iproj,plat,plon,rota,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
! call supmap(iproj,plat+0.001,plon,rota-0.001,pl1,pl2,pl3,pl4,2,30,4,0,ierr)
  if (ierr.ne.0) then
     print*, 'supmap ierr = ', ierr
         stop
!    stop
  endif
  call getset(xl,xr,xb,xt,wl,wr,wb,wt,ml)

  write(hlev,'(I8)') ilev

  call set(0., 1., 0., 1., 0., 1., 0., 1., 1)
  call pchiqu(0.1, xb-0.04, hlev//'  '//field, .020, 0.0, -1.0)
  print*, field//'#'//units//'#'//trim(Desc)
! call pchiqu(0.1, xb-0.12, Desc, .012, 0.0, -1.0)
  hunit = '                                      '
  ih = 0
  do i = 1, len(units)
     if (units(i:i).eq.'{') then
        hunit(ih+1:ih+3) = '~S~'
        ih = ih + 3
        elseif (units(i:i).eq.'}') then
        hunit(ih+1:ih+3) = '~N~'
        ih = ih + 3
     else
        ih = ih + 1
        hunit(ih:ih) = units(i:i)
     endif
  enddo
  if ( ifv .le. 3 ) then
    tmp32 = 'MM5 intermediate format'
  else if ( ifv .eq. 4 ) then
    tmp32 = 'SI intermediate format'
  else if ( ifv .eq. 5 ) then
    tmp32 = 'WPS intermediate format'
  endif
  call pchiqu(0.1, xb-0.08, hunit, .015, 0.0, -1.0)
  call pchiqu(0.1, xb-0.12, Desc, .013, 0.0, -1.0)
  call pchiqu(0.6, xb-0.12, source, .013, 0.0, -1.0)
  call pchiqu(0.1, xb-0.15, tmp32, .013, 0.0, -1.0)

  call set(xl,xr,xb,xt,1.,float(ix),1.,float(jx),ml)

  call CPSETI ('SET - Do-SET-Call Flag', 0)
  call CPSETR ('SPV - Special Value', -1.E30)

  if (dy.lt.0.) then
     call array_flip(scr2d, ix, jx)
  endif

  call cpseti('LLP', 3)
  call cprect(scr2d,ix,ix,jx,rwrk,lwrk,iwrk,liwk)
  call cpcldr(scr2d,rwrk,iwrk)
  call cplbdr(scr2d,rwrk,iwrk)

  if (dy.lt.0.) then
     call array_flip(scr2d, ix, jx)
  endif

  call frame

end subroutine plt2d

subroutine stopit
  call graceful_stop
end

subroutine graceful_stop
  call gdawk(2)
  call gdawk(1)
  call gclwk(2)
  call gclwk(1)
  call gclks
  print*, 'Graceful Stop.'
  stop
end subroutine graceful_stop

subroutine fmtxyll(x, y, xlat, xlon, project, glat1, glon1, gclon,&
     gtrue1, gtrue2, gdx, gdy)
  implicit none

  real , intent(in) :: x, y, glat1, glon1, gtrue1, gtrue2, gdx, gdy, gclon
  character(len=2), intent(in) :: project
  real , intent(out) :: xlat, xlon

  real :: gx1, gy1, gkappa
  real :: grrth = 6370.

  real :: r, y1
  integer :: iscan, jscan
  real, parameter :: pi = 3.1415926534
  real, parameter :: degran = pi/180.
  real, parameter :: raddeg = 180./pi
  real :: gt

  if (project.eq.'CE') then  ! Cylindrical Equidistant grid

     xlat = glat1 + gdy*(y-1.)
     xlon = glon1 + gdx*(x-1.)
     
  elseif (project == "ME") then

     gt = grrth * cos(gtrue1*degran)
     xlon = glon1 + (gdx*(x-1.)/gt)*raddeg
     y1 = gt*alog((1.+sin(glat1*degran))/cos(glat1*degran))/gdy
     xlat = 90. - 2. * atan(exp(-gdy*(y+y1-1.)/gt))*raddeg

  elseif (project.eq.'ST') then  ! Polar Stereographic grid

     r = grrth/gdx * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
     gx1 = r * sind(glon1-gclon)
     gy1 = - r * cosd(glon1-gclon)

     r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
     xlat = 90. - 2.*atan2d((r*gdx),(grrth*(1.+sind(gtrue1))))
     xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon

  elseif (project.eq.'LC') then  ! Lambert-conformal grid

     call glccone(gtrue1, gtrue2, 1, gkappa)

     r = grrth/(gdx*gkappa)*sind(90.-gtrue1) * &
          (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
     gx1 =  r*sind(gkappa*(glon1-gclon))
     gy1 = -r*cosd(gkappa*(glon1-gclon))

     r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
     xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
          ((r*gkappa*gdx)/(grrth*sind(90.-gtrue1)))**(1./gkappa))
     xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon

  else

     write(*,'("Unrecoginzed projection: ", A)') project
     write(*,'("Abort in FMTXYLL",/)')
     stop

  endif
contains
  real function sind(theta)
    real :: theta
    sind = sin(theta*degran)
  end function sind
  real function cosd(theta)
    real :: theta
    cosd = cos(theta*degran)
  end function cosd
  real function tand(theta)
    real :: theta
    tand = tan(theta*degran)
  end function tand
  real function atand(x)
    real :: x
    atand = atan(x)*raddeg
  end function atand
  real function atan2d(x,y)
    real :: x,y
    atan2d = atan2(x,y)*raddeg
  end function atan2d

  subroutine glccone (fsplat,ssplat,sign,confac)
    implicit none
    real, intent(in) :: fsplat,ssplat
    integer, intent(in) :: sign
    real, intent(out) :: confac
    if (abs(fsplat-ssplat).lt.1.E-3) then
       confac = sind(fsplat)
    else
       confac = log10(cosd(fsplat))-log10(cosd(ssplat))
       confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
            log10(tand(45.-float(sign)*ssplat/2.)))
    endif
  end subroutine glccone

end subroutine fmtxyll

subroutine array_flip(array, ix, jx)
  implicit none
  integer :: ix, jx
  real , dimension(ix,jx) :: array

  real, dimension(ix) :: hold
  integer :: i, j, jj

  do j = 1, jx/2
     jj = jx+1-j
     hold = array(1:ix, j)
     array(1:ix,j) = array(1:ix,jj)
     array(1:ix,jj) = hold
  enddo
end subroutine array_flip