[gradsusr] Regarding uv2trw function

Mike Fiorino michael.fiorino at noaa.gov
Mon Oct 26 11:42:32 EDT 2020


oops...forgot to attach the fortran...

hi Raghu,

sorry for the confusion about this extension...(I cc'd gradsusr)

attached is my version of the test script that does what you're looking 
for (mf-test.gs)

I've also attached the fortran code (ftn_uv2trw.F) and a couple of plots

a couple of notes...

1) the center in your script is way off from the wind speed minimum.  
the script shows how to get that center interactively by clicking on the 
wind speed plot.  this will make a very big difference in both the 
tangential and radial winds.  using the wind speed minimum, the diff is 
< 5 kts, but with your center > 80 kts!

2) the doc is correct...you used lon,lat but should be lat,lon  see below

3) I'm running this version on CentOs 7:

Config: v2.2.1.oga.1 little-endian readline grib2 netcdf hdf4-sds hdf5 
opendap-grids,stn athena geotiff shapefile
Grid Analysis and Display System (GrADS) Version 2.2.1.oga.1

you should update...

thanks for the enquiry...  had some GrADS fun this morning ;)

CC4N & best /R Mike


ga-> d uv2trw

Error from uv2trw: Too few args

           Usage:  uv2trw(u,v,latc,lonc,[iout]))

returns the tangential velocity (iout=1; the default) relative to a 
fixed lonc/latc
if no lonc/latc is provided, then the center of the grid is assumed the 
fixed point

Examples:

    ga-> d uv2trw(ua,va,latc,lonc,3)     (return the xy vcomp tangential 
wind
    ga-> d uv2trw(ua,va,latc,lonc,2)     (return the xy ucomp tangential 
wind
    ga-> d uv2trw(ua,va,latc,lonc,1)     (return the total tangential wind
    ga-> d uv2trw(ua,va,latc,lonc,-1)    (return the total radial wind
    ga-> d uv2trw(ua,va,latc,lonc,-2)    (return the xy ucomp radial wind
    ga-> d uv2trw(ua,va,latc,lonc,-3)    (return the xy vcomp radial wind




On 10/26/2020 01:29, Raghavendra Ashrit wrote:
> Dear Michael
>
> I am writing to you since I could not find any discussion or reply on
> gradsuser mail group.
>
> I am trying to use a function uv2trw developed and implemented by you in
> opengrads for a recent super cyclone 'Amphan' case over India. I am using
> opengrads, with all env correctly specified.
> ------------------
> Grid Analysis and Display System (GrADS) Version 2.0.2.oga.1
> Copyright (c) 1988-2011 by Brian Doty and the
> Institute for Global Environment and Society (IGES)
> -------------------------
> 'q udf' lists out all the stuff including the following
> uv2trw Find radial/tangential velocity f_uv2trw@^libmf.gex
>
> When I try to use it the function to obtain tangential wind for a level
> (pl see test.gs and the uvt.nc shared on google drive)
> https://drive.google.com/drive/folders/1AD1U-h-vWHBr-YAcOW-ko0DHax3hvYXP?usp=sharing
> I just get v wind at that level. Same with radial wind. Let me know if I
> am missing something. If possible put me in touch with who can assist.
>
> Raghu
>
>
>
>



-------------- next part --------------
function main(args)

rc=gsfallow(on)
rc=const()

# -- open wind file
#
'sdfopen uvt.nc'

latb=19
late=23
lonb=85
lone=91

'set lat 'latb' 'late
'set lon 'lonb' 'lone
'set z 20'

# -- tracker center?
#
latc=20.28
lonc=87.66

getwmin=0

# -- wind min center
#
wlonc=87.84
wlatc=20.88

vskip=3

# -- plot the wind [kts]

'u=u*'_ms2kt
'v=v*'_ms2kt
'w=mag(u,v)'

'set gxout contour'
'set cint 10'
'd w'

'q w2xy 'lonc' 'latc
xtc=subwrd(result,3)
ytc=subwrd(result,6)

print 'CCCCC 'xtc' 'ytc
'draw wxsym 41 'xtc' 'ytc' 0.50 2'

# -- plot wind min interactive when getwmin=1
#
'q w2xy 'wlonc' 'wlatc
xtcw=subwrd(result,3)
ytcw=subwrd(result,6)
'draw wxsym 41 'xtcw' 'ytcw' 0.35 1'

# -- get wind min interactively when getwmin=1
#
if(getwmin)

  'q xy2w 'xw' 'yw

  wlonc=subwrd(result,3)
  wlatc=subwrd(result,6)

  print 'wind min lon: 'wlonc' lat: 'wlatc
  'q pos'
  'quit'

endif

# -- wind barbs

'set gxout barb'
'set digsiz 0.05'
'd u;skip(v,'vskip')'

# -- title and output
#
'draw title total wind [kts]'
'gxprint totwind.png x1024 y768'

'q pos'
'c'

# -- now get the rad/tan wind


# -- tangential wind speed
#
'vtw=uv2trw(u,v,'wlatc','wlonc',1)'
'vt=uv2trw(u,v,'latc','lonc',1)'

'set gxout contour'
'd vtw'

# -- u/v components of the tangential wind
#
'uvt=uv2trw(u,v,'wlatc','wlonc',2)'
'vvt=uv2trw(u,v,'wlatc','wlonc',3)'

'set gxout barb'
'd uvt;skip(vvt,3)'
'draw title tangential wind [kts]'
'gxprint tanwind.png x1024 y768'

'q pos'
'c'


# -- radial wind
#
'vrw=uv2trw(u,v,'wlatc','wlonc',-1)'
'set gxout contour'
'd vrw'

# -- u/v comp of radial wind
#
'uvr=uv2trw(u,v,'wlatc','wlonc',-2)'
'vvr=uv2trw(u,v,'wlatc','wlonc',-3)'

'set gxout vector'
'd uvr;skip(vvr,2)'

'draw title radial wind [kts]'
'gxprint radwind.png x1024 y768'

'q pos'
'quit'


'set gxout contour'
'set ccolor rainbow'
'd vt'

'set gxout vector'
'set arrscl 0.5 30'
'set ccolor 2'
'd uvt;vvt'


'q pos'
'quit'

'set ccolor 3'
'set arrscl 0.5 30'
'd uvr;vvr'

'set ccolor 1'
'set arrscl 0.5 30'
'd u;v'

return


function getloc(result)

return xmin

-------------- next part --------------
A non-text attachment was scrubbed...
Name: totwind.png
Type: image/png
Size: 188078 bytes
Desc: not available
URL: <http://gradsusr.org/pipermail/gradsusr/attachments/20201026/016319f7/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: tanwind.png
Type: image/png
Size: 195206 bytes
Desc: not available
URL: <http://gradsusr.org/pipermail/gradsusr/attachments/20201026/016319f7/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: radwind.png
Type: image/png
Size: 164642 bytes
Desc: not available
URL: <http://gradsusr.org/pipermail/gradsusr/attachments/20201026/016319f7/attachment-0005.png>
-------------- next part --------------
      subroutine ftnuv2trw(u, v, lat, lon, dum, undef, 
     $     ni, nj, centlat, centlon, 
     $     iout, rc)

      USE gex
      USE const

      implicit integer(i-n)
      implicit real(kind=iGaKind) (a-h,o-z)

      real(kind=iGaKind) ::  u(ni,nj),v(ni,nj)
      real(kind=iGaKind) ::  dum(ni,nj)
      real(kind=iGaKind) ::  lon(ni),lat(nj)
      
      integer rc

C         initialize dum to input field
C

      do i=1,ni
        xlon=lon(i)

        do j=1,nj

          xlat=lat(j)
          udat=u(i,j)
          vdat=v(i,j)

          call getvrvt (centlon,centlat,xlon,xlat,
     $         udat,vdat,
     $         vr,vt,uvr,vvr,uvt,vvt,
     $         igvtret)

          if(iout == 3) dum(i,j)=vvt
          if(iout == 2) dum(i,j)=uvt
          if(iout == 1) dum(i,j)=vt

          if(iout == -1) dum(i,j)=vr
          if(iout == -2) dum(i,j)=uvr
          if(iout == -3) dum(i,j)=vvr

        enddo
      enddo
      
      do i=1,ni
        do j=1,nj
          u(i,j)=dum(i,j)
        end do
      end do

      rc=0
      return
      end



c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
      subroutine getvrvt (centlon,centlat,xlon,xlat,
     &     udat,vdat,
     &     vr,vt,uvr,vvr,uvt,vvt,
     &     igvtret)
c
c     ABSTRACT: This subroutine takes as input a u-wind and v-wind value
c     at an input (xlon,xlat) location and returns the tangential and
c     radial wind components relative to the input center lat/lon 
c     position (centlon,centlat).  The only trick to this whole 
c     subroutine is figuring out the angle from the center point to the
c     data point, and we do this by creating a triangle with the leg 
c     from the center point to the data point being the hypotenuse.
c
c     NOTE: All longitudes must be in positive degrees east (0-360) !!!
c
c     INPUT:
c     centlon  Longitude of center point
c     centlat  Latitude  of center point
c     xlon     Longitude of pt at which vr & vt will be computed
c     xlat     Latitude  of pt at which vr & vt will be computed
c     udat     u-value of wind at the point (xlon,xlat) 
c     vdat     v-value of wind at the point (xlon,xlat) 
c
c     OUTPUT:
c     vr      Radial wind component at (xlon,xlat) wrt (centlon,centlat)
c     vt      Tang   wind component at (xlon,xlat) wrt (centlon,centlat)
c     igvtret Return code from this subroutine

      USE const
      USE gex

      implicit integer(i-n)
      implicit real(kind=iGaKind) (a-h,o-z)

      logical verb

      verb=.false.
CCC      verb=.true.

      call calcdist1(centlon,centlat,xlon,xlat,hyp_dist)

      xlatdiff = abs(centlat - xlat)
      xlondiff = abs(centlon - xlon)

      
      if (xlondiff == 0 .and. xlatdiff > 0) then

        if (centlat > xlat) psi = 180   ! pt directly south of ctr
        if (centlat < xlat) psi = 0     ! pt directly north of ctr

      else if (xlondiff > 0 .and. xlatdiff == 0) then

        if (centlon > xlon) psi = 270   ! pt directly west of ctr
        if (centlon < xlon) psi = 90    ! pt directly east of ctr
C         
C         exactly on the center
C
      else if (hyp_dist ==  0) then

        vr=0.0
        vt=0.0
        uvr=0.0
        vvr=0.0
        uvt=0.0
        vvt=0.0
        return

      else

        ! This next part figures out the angle from the center point
        ! (centlon,centlat) to the data point (xlon,xlat).  It does 
        ! this by setting up a triangle and then using inverse trig
        ! functions to get the angle.  Since this is a kludgy way to
        ! do it that doesn't account for the curvature of the earth,
        ! we'll do it 2 ways, using asin and then acos, then take the
        ! average of those 2 for the angle.  hyp_dist, calculated just
        ! above, is the distance from the center pt to the data pt.

        opp_dist  = xlatdiff/360. * earthcircum
        sin_value = opp_dist / hyp_dist

        if (sin_value > 1.0) then

          sin_value=1.0

          if(verb) then
            print *,' '
            print *,'           !!! In getvrvt, sin_value > 1, setting to 1.'
            print *,'           !!! opp_dist= ',opp_dist,' hyp_dist= ',hyp_dist
            print *,'           !!! sin_value = ',sin_value
            print *,'           !!! centlon= ',centlon,' centlat= ',centlat
            print *,'           !!! xlon=    ',xlon,' xlat=    ',xlat
            print *,' '
          endif

        endif

        sin_angle = asin(sin_value) * rad2deg

        call calcdist1(centlon,centlat,xlon,centlat,adj_dist)
        cos_value = adj_dist / hyp_dist

        if (cos_value > 1.0) then

          cos_value = 1.0

          if(verb) then
            print *,' '
            print *,'           !!! In getvrvt, cos_value > 1, setting to 1.'
            print *,'           !!! adj_dist= ',adj_dist,' hyp_dist= ',hyp_dist
            print *,'           !!! cos_value = ',cos_value
            print *,'           !!! centlon= ',centlon,' centlat= ',centlat
            print *,'           !!! xlon=    ',xlon,' xlat=    ',xlat
            print *,' '
          endif

        endif

        cos_angle = acos(cos_value) * rad2deg

        theta = 0.5 * (sin_angle + cos_angle)

        ! The previous lines of code just calculated an angle between
        ! 0 and 90.  This next if structure adjusts that angle to 
        ! instead be between 0 and 360.

        if      (centlat <= xlat .and. centlon <= xlon) then
          psi = 90 - theta
        else if (centlat >  xlat .and. centlon <= xlon) then
          psi = 90 + theta
        else if (centlat >= xlat .and. centlon >= xlon) then
          psi = 270 - theta
        else if (centlat <  xlat .and. centlon >= xlon) then
          psi = 270 + theta
        endif

      endif

      psi=psi*deg2rad
      uvrcomp = udat * sin(psi)
      vvrcomp = vdat * cos(psi)
      vr      = uvrcomp + vvrcomp

      uvtcomp = (-udat) * cos(psi)
      vvtcomp = vdat    * sin(psi)
      vt      = uvtcomp + vvtcomp

      vvr = vr * cos (psi)
      uvr = vr * sin (psi)
    
      vvt = -vt * cos (psi + pi*0.5)
      uvt = -vt * sin (psi + pi*0.5)

      return 
      end


      subroutine calcdist1(rlonb,rlatb,rlonc,rlatc,xdist)

c
c     ABSTRACT: This subroutine computes the distance between two 
c               lat/lon points by using spherical coordinates to 
c               calculate the great circle distance between the points.
c                       Figure out the angle (a) between pt.B and pt.C,
c             N. Pole   then figure out how much of a % of a great 
c               x       circle distance that angle represents.
c              / \
c            b/   \     cos(a) = (cos b)(cos c) + (sin b)(sin c)(cos A)
c            /     \
c        pt./<--A-->\c     NOTE: The latitude arguments passed to the
c        B /         \           subr are the actual lat vals, but in
c                     \          the calculation we use 90-lat.
c               a      \
c                       \pt.  NOTE: You may get strange results if you:
c                         C    (1) use positive values for SH lats AND
c                              you try computing distances across the 
c                              equator, or (2) use lon values of 0 to
c                              -180 for WH lons AND you try computing
c                              distances across the 180E meridian.
c    
c     NOTE: In the diagram above, (a) is the angle between pt. B and
c     pt. C (with pt. x as the vertex), and (A) is the difference in
c     longitude (in degrees, absolute value) between pt. B and pt. C.
c
c     !!! NOTE !!! -- THE PARAMETER earthcircum IS DEFINED (AS OF THE 
c     ORIGINAL WRITING OF THIS SYSTEM) IN KM, NOT M, SO BE AWARE THAT
c     THE DISTANCE RETURNED FROM THIS SUBROUTINE IS ALSO IN KM.
c         
      USE gex
      USE const
      

      implicit integer(i-n)
      implicit real(kind=iGaKind) (a-h,o-z)

      if (rlatb < 0.0 .or. rlatc < 0.0) then
        pole = -90.
      else
        pole = 90.
      endif

      distlatb = (pole - rlatb) * deg2rad
      distlatc = (pole - rlatc) * deg2rad
      difflon  = abs( (rlonb - rlonc)*deg2rad )

      cosanga = ( cos(distlatb) * cos(distlatc) + 
     &            sin(distlatb) * sin(distlatc) * cos(difflon))
 
c     This next check of cosanga is needed since I have had ACOS crash
c     when calculating the distance between 2 identical points (should
c     = 0), but the input for ACOS was just slightly over 1
c     (e.g., 1.00000000007), due to (I'm guessing) rounding errors.

      if (cosanga > 1.0) then
        cosanga = 1.0
      endif

      degrees    = acos(cosanga) / deg2rad
      circ_fract = degrees / 360.
      xdist      = circ_fract * earthcircum
c
c     NOTE: whether this subroutine returns the value of the distance
c           in km or m depends on the scale of the parameter earthcircum. 
c           At the original writing of this subroutine (7/97), earthcircum
c           was given in km.
c
      return
      end


More information about the gradsusr mailing list