20020131: radius of earth bug in imar2gm.f -- second notice

David,

I updated the imar2gm.f recently to add the TANC projection as well as
fix an outstating bug in PS. The PS section I fixed does use the
radius of earth from the MciDAS file. I can make the other RADIUS to re
changes. There could still be some discrepancy in the lack of eccentricity
use however.

Steve Chiswell
Unidata User SUpport



>From: David Ovens <address@hidden>
>Organization: UCAR/Unidata
>Keywords: 200201312147.g0VLlkx05878

>GEMPAK Support,
>
>I was disappointed to see that my bug fixes from Dec. 2000 for
>gemlib/im/imar2gm.f have not fully made it into GEMPAK 5.6.E.1 or
>previous GEMPAK versions.  The crux of our differences now are:
>
> 1) in the actual use of the radius of the earth that
>    is in the McIDAS file (re) vs GEMPAK's value (RADIUS), and
> 2) the location of re in the ianav array for projections other than
>    'MERC' -- I have re = ianav(8) for everything but 'MERC' where I
>    have ianav(7), whereas you have re = ianav(7) for everything.
>
>Unless, McIDAS code has changed, I think your version is still in
>error.  Below you will find an excerpt from my earlier email and
>below that my newest version of imar2gm.f.
>
>Please let me know what you think, I'd like to get this incorporated
>into the standard release.  This is the only way that my McIDAS area
>files of topography jibe with the map database.
>
>David
>-- 
>
>David Ovens            e-mail: address@hidden
>(206) 685-8108          plan: Real-time MM5 forecasting for Pacific Northwest
>Research Meteorologist
>Dept of Atmospheric Sciences, Box 351640
>University of Washington 
>Seattle, WA  98195
>
>
>>From ovens Thu Dec 28 06:58:44 2000
>Subject: I found and fixed a bug in GEMPAK imar2gm.f
>To: address@hidden, justin (Justin Sharp), harry (Harry Edmon)
>Date: Thu, 28 Dec 2000 06:58:44 -0800 (PST)
>X-Mailer: ELM [version 2.5 PL2]
>Content-Length: 16228     
>Status: OR
>
>Hello,
>
>I've been attempting to modify Justin Sharp's excellent program that
>makes McIdas area files out of 30-second terrain data.  (A program that
>he has placed on CD with the terrain data set).  It has always worked
>fine with the RECTilinear projection of McIdas, but when I tried to
>use IMGREMAP as in the following code extract
>  $imgremap = "imgremap.k MYDATA/IMAGES.$area MYDATA/IMAGES.9998 PRO=MERC LATL
> ON=$center_lat $center_lon  RES=1 SIZE=ALL"; 
>I was having trouble with the GEMPAK map being offset with the
>terrain.  In McIdas, however, I could not discern any problem with the
>remapped file. 
>
>Well, I have gone to the McIdas source and to the GEMPAK source and I
>have been able to fix the problem in GEMPAK-5.6.A by correcting a
>discrepancy between the Earth radius used in GEMPAK and that used in
>the McIdas file.  I have modified $GEMPAK/source/gemlib/im/imar2gm.f
>as you will see below to replace RADIUS with re (obtained from the
>navigation header of the file).
>
>I have tested this with MERC, LAMB, and PS terrain/McIdas files and it
>all seems to line up fine using sfmap.
>
>David
>========================= imar2gm.f ====================
>       SUBROUTINE IM_AR2GM  ( imgfil, rnvblk, iret )
>C************************************************************************
>C* IM_AR2GM                                                            *
>C*                                                                     *
>C* This subroutine reads the header information from a MCIDAS AREA file
>       *
>C* and sets the navigation                                             *
>C*                                                                     *
>C* IM_AR2GM  ( IMGFIL, RNVBLK, IRET )                                  *
>C*                                                                     *
>C* Input parameters:                                                   *
>C*     IMGFIL          CHAR*           Image file                      *
>C*                                                                     *
>C* Output parameters:                                                  *
>C*     RNVBLK (LLNNAV) REAL            Navigation block                *
>C*     IRET            INTEGER         Return code                     *
>C*                                       0 = normal return             *
>C*                                      -2 = Error opening/reading file*
>C*                                      -3 = Invalid image file format *
>C*                                      -4 = Invalid image file header *
>C*                                      -5 = Invalid image navigation  *
>C**                                                                    *
>C* Log:                                                                
>       *
>C* J. Cowie/COMET       2/95                                           *
>C* C. Lin/EAI           6/95   Change input from file name to lunmf    *
>C*                             Change GG_* to IM_* for projection setup*
>C* J. Cowie/COMET      10/95   Set xy image scaling to 1.0             *
>C* J. Cowie/COMET       4/96   Add calibration type, check for 'RAW '  *
>C* T. Lee/GSC           7/96   Renamed IM_ARHD; Cleaned up & eliminated*
>C*                             map projection calls                    *
>C* D. Keiser/GSC        8/96   Added RETURN after call to GSATMG       *
>C* T. Lee/GSC           9/96   Used RADIUS and floating point          *
>C* D.W.Plummer/NCEP    10/96   Fix RADR processing bug, rem to RADIUS  *
>C* T. Lee/GSC          11/96   Declared timfil; Used meter everywhere  *
>C* J. Cowie/COMET       1/97   Changed IMGDEF common variables names   *
>C* S. Chiswell/UNIDATA  1/97   Add extra scaling chck for RECT images  *
>C* J. Cowie/COMET       2/97   Fixed RADR off-by-one error             *
>C* J. Cowie/COMET      12/97   Changed jyear to 1900 + current year    *
>C* J. Cowie/COMET      12/97   Added imradf, removed immode, imraw     *
>C* S. Danz/AWC          05/99   Fix to RECT off-by-one error            *
>C* S. Jacobs/NCEP       6/00   Used correct image name in GSATMG call  *
>C* S. Danz/AWC         10/00   Fix RECT setting center longitude       *
>C* J. Cowie/COMET      12/00   Added MCRADR projection type            *
>C* D. Ovens/UW          12/00   Fix MERC,LAMB,PS Earth Radius mismatches*
>C* D.W.Plummer/NCEP     9/01   Set imbswp=0 when no byte swapping      *
>C* S. Chiswell/Unidata 10/01   Fix PS dlat calculation                 *
>C************************************************************************
>       INCLUDE         'GEMPRM.PRM'
>       INCLUDE         'IMGDEF.CMN'
>       INCLUDE         'AREAFL.CMN'
>C*
>       CHARACTER*(*)   imgfil
>       REAL            rnvblk (*)
>C*
>       CHARACTER       navtyp*4, proj*3
>       INTEGER         iarray (704), idtarr (3)        
>C*
>       EQUIVALENCE     (iarray (1),  iadir (1))
>       EQUIVALENCE     (iarray (65), ianav (1))
>C------------------------------------------------------------------------
>       iret = 0
>       sign = 1.
>C
>C*     Open the file and read AREA DIR and NAV blocks
>C
>       CALL FL_DOPN  ( imgfil, 704, .false., lunmf, iret )
>       CALL FL_READ  ( lunmf, 1, 704, iarray, ier )
>       IF  ( ier .ne. 0 )  THEN
>               iret = -1
>               RETURN
>       END IF
>       CALL FL_CLOS ( lunmf, ier )
>C
>C*     Check to see if the header items needs byte-swapping.
>C
>       IF ( iadir (2) .ne. 4 ) THEN
>C      
>           imbswp = 1
>C
>C*         Swap bytes in AREA DIR except comment.
>C
>           ier = MV_SWP4 ( 24, iadir (  1 ), iadir (  1 ) )
>           ier = MV_SWP4 ( 19, iadir ( 33 ), iadir ( 33 ) )
>           ier = MV_SWP4 (  3, iadir ( 54 ), iadir ( 54 ) )
>           ier = MV_SWP4 (  7, iadir ( 58 ), iadir ( 58 ) )
>C
>C*         Swap bytes in NAV block except 1st word (nav label).
>C
>           ier = MV_SWP4 ( 640-1, ianav ( 2 ) , ianav ( 2 ) )
>C          
>       ELSE
>C
>           imbswp = 0
>C
>       END IF
>C
>C*     Determine calibration type      Calibrations:   'RAW'
>C*                                                     'BRIT'
>C*                                                     'TEMP'
>C
>       CALL ST_ITOC ( iadir (53), 1, cmcalb, ier )
>C
>C*     Parse AREA data header. Check for radar image.
>C
>       imsorc  = iadir ( 3 )
>       IF ( imsorc .eq. 7 ) then
>           imradf = 1
>           imtype = iadir ( 22 )
>         ELSE
>           imradf = 0
>           imtype = iadir ( 19 )
>       END IF
>C      
>       jyear = ( iadir ( 4 ) / 1000 ) + 1900
>       jday  = MOD ( iadir ( 4 ), 1000 )
>       CALL TI_JTOI ( jyear, jday, idtarr, ier )
>       imdate = idtarr ( 1 ) * 10000 + idtarr ( 2 ) * 100 + idtarr ( 3 )
>       imtime = iadir (  5 )
>       iullne = iadir (  6 )
>       iulele = iadir (  7 )
>       imnlin = iadir (  9 )
>       imnpix = iadir ( 10 )
>       imdpth = iadir ( 11 )
>       rmyres = iadir ( 12 )
>       rmxres = iadir ( 13 )
>       imnchl = iadir ( 14 )
>       imprsz = iadir ( 15 )
>       imdcsz = iadir ( 49 )
>       imclsz = iadir ( 50 )
>       imlvsz = iadir ( 51 )
>       imvald = iadir ( 36 )
>       icaloff = iadir ( 63 )
>       if(icaloff.gt.0) CALL IM_CALIB(imgfil,icaloff,iret)
>       rmxysc = 1.0
>C
>C*      Determine offset to start of data
>C
>       IF ( iadir ( 34 ) .gt. 0 ) THEN
>            imdoff = iadir ( 34 )
>        ELSE
>            iret = -4
>            RETURN
>        END IF
>        imldat = imnlin * ( imnpix * imdpth * imnchl + imprsz )
>C
>C*      Set image subset to full image for now.
>C
>        imleft = 1
>        imtop  = 1
>        imrght = imnpix
>        imbot  = imnlin
>C
>C*     Set navigation.
>C
>C*     Standard latitude.
>C
>       angdd  = ianav ( 4 ) / 10000
>       angmm  = ianav ( 4 ) / 100 - angdd * 100
>       angss  = ianav ( 4 ) - angdd * 10000 - angmm * 100
>       stdlat = angdd + angmm / 60. + angss / 3600.
>       phi0   = stdlat
>       phi0r  = phi0 * DTR
>       IF ( phi0r .lt. 0. ) sign = - 1.
>C
>C*     Normal longitude. If west positive, make west negative.
>C
>       angdd = ianav ( 6 ) / 10000
>       angmm = ianav ( 6 ) / 100 - angdd * 100
>       angss = ianav ( 6 ) - angdd * 10000 - angmm * 100
>       clon  = angdd + angmm / 60. + angss / 3600.
>       IF ( ianav ( 10 ) .ge. 0 ) clon = - clon
>C
>C*      Set pixel/grid spacing and earth radius.
>C
>       dx    = ianav ( 5 ) * iadir ( 13 )
>cdo    re    = ianav ( 7 )
>       re    = ianav ( 8 )
>       if (re .le. 6200000.) re = RADIUS
>C
>C*     Set pixel/grid dimension.
>C
>        rgx =   imnpix
>        rgy =   imnlin
>C
>C*      Location of pole point (rxp, ryp); (1,1) at lower-left corner.
>C
>        rxp = float ( ianav ( 3 ) - iadir ( 7 ) ) / iadir ( 13 ) + 1.
>       ryp = rgy - float ( ianav ( 2 ) - iadir ( 6 ) ) / iadir ( 12 )
>C
>       CALL ST_ITOC ( ianav, 1, navtyp, ier )
>C
>C*     ******************************************
>C*     *** Satellite projection (MCIDAS GOES) ***
>C*     ******************************************
>C
>       IF  ( navtyp .eq. 'GOES' .or. navtyp .eq. 'GVAR' .or.
>     +       navtyp .eq. 'RADR')  THEN
>C
>           CALL GSATMG ( imgfil, iadir, ianav, imleft, imtop, imrght,
>     +                   imbot, ier )
>           RETURN
>C
>C*     ***************************
>C*     *** Polar Stereographic ***
>C*     ***************************
>C
>         ELSE IF  ( navtyp .eq. 'PS' )  THEN
>C
>C*         Compute lat/lon of the lower-left corner.
>C
>
>C          squished projection, need (12) LINERES
>           dy = ianav ( 5 ) * rmyres
>
>
>           dxp = 1. - rxp
>           dyp = 1. - ryp
>           alpha = 1. + SIN ( ABS ( phi0r ) )
>           rm  = dy * SQRT ( dxp * dxp + dyp * dyp ) / alpha
>cdo        dlat1 = sign * ( HALFPI - 2. * ATAN ( rm / RADIUS ) ) * RTD
>           dlat1 = sign * ( HALFPI - 2. * ATAN ( rm / re ) ) * RTD
>           IF ( dyp .ne. 0. ) THEN
>               dyp   = - dyp * sign
>               thta  = ATAN2 ( dxp, dyp ) * RTD
>               dlon1 = clon + thta
>               CALL PRNLON ( 1, dlon1, ier )
>             ELSE
>               dlon1 = clon
>           END IF
>C
>C*         Compute lat/lon of the upper-right corner.
>C
>           dxp = rgx - rxp
>           dyp = rgy - ryp
>           rm  = dy * SQRT ( dxp * dxp + dyp * dyp ) / alpha
>cdo        dlat2 = sign * ( HALFPI - 2. * ATAN ( rm / RADIUS ) ) * RTD
>           dlat2 = sign * ( HALFPI - 2. * ATAN ( rm / re ) ) * RTD
>           IF ( dyp .ne. 0. ) THEN
>               dyp = - dyp * sign
>               thta = ATAN2 ( dxp, dyp ) * RTD
>               dlon2 = clon + thta
>               CALL PRNLON ( 1, dlon2, iret )
>             ELSE
>               dlon2 = clon
>           END IF
>C
>C*         Set map projection angles. Angle 2 is the central longitude.
>C
>           IF ( stdlat .ge. 0. ) THEN
>               dangl1 = 90.
>             ELSE
>               dangl1 = - 90.
>           END IF
>           dangl2 = clon
>           dangl3 = 0.
>           proj   = 'STR'
>C
>C*       ****************
>C*       *** Mercator ***    
>C*       ****************
>C
>         ELSE IF  ( navtyp .eq. 'MERC' )  THEN
>C
>C*      Reset earth radius and set eccentricity.
>C
>            re    = ianav ( 7 )
>            if (re .le. 6200000.) re = RADIUS
>            ecc   = ianav ( 8 ) / 1000000.
>C
>C*         Compute lat/lon of the lower-left corner.
>C
>           dxp   = 1. - rxp
>           dyp   = 1. - ryp
>           rm    = dx * dyp
>cdo        rcos  = RADIUS * COS ( phi0r )
>           rcos  = re * COS ( phi0r )
>           arg   = EXP ( rm / rcos )
>           dlat1 = ( 2. * ATAN ( arg ) - HALFPI ) * RTD
>           dlon1 = clon + ( dx * dxp / rcos ) * RTD
>           CALL PRNLON ( 1, dlon1, ier )
>C
>C*         Compute lat/lon of the upper-right corner point.
>C
>           dxp = rgx - rxp
>            dyp = rgy - ryp
>           rm  = dx * dyp
>           arg = EXP ( rm / rcos )
>           dlat2 = ( 2. * ATAN ( arg ) - HALFPI ) * RTD
>           dlon2 = clon + ( dx * dxp / rcos ) * RTD
>           CALL PRNLON ( 1, dlon2, ier )
>C
>C*         Set the map projection and angles.
>C
>           dangl1 = 0.
>           dangl2 = clon
>           dangl3 = 0.
>           proj   = 'MER'
>C
>C*       *************************
>C*       *** Lambert Conformal ***
>C*       *************************
>C
>         ELSE IF  ( navtyp .eq. 'LAMB' )  THEN
>C
>C*         Standard latitude #1
C
>           std1 = stdlat
>C
>C*         Standard latitude #2
>C
>           angdd = ianav ( 5 ) / 10000
>           angmm = ianav ( 5 ) / 100 - angdd * 100
>           angss = ianav ( 5 ) - angdd * 10000 - angmm * 100
>           std2  = angdd + angmm / 60. + angss / 3600.
>C
>C*         Normal longitude. If west positive, make west negative.
>C
>           angdd = ianav ( 7 ) / 10000
>           angmm = ianav ( 7 ) / 100 - angdd * 100
>           angss = ianav ( 7 ) - angdd * 10000 - angmm * 100
>           clon  = angdd + angmm / 60. + angss / 3600.
>           IF ( ianav ( 10 ) .ge. 0 ) clon = - clon
>C
>C*         Compute pixel/grid spacing and colatitudes.
>C
>           dx   = ianav ( 6 ) * iadir ( 13 )
>           psi1 = HALFPI - ABS ( std1 ) * DTR
>           psi2 = HALFPI - ABS ( std2 ) * DTR
>C
>C*          Compute cone constant.
>C
>           IF ( psi1 .eq. psi2 ) THEN
>               ccone = COS ( psi1 )
>             ELSE
>               tmp1 = ALOG ( SIN ( psi2 ) / SIN ( psi1 ) )
>               tmp2 = ALOG ( TAN ( psi2 / 2. ) / TAN ( psi1 / 2. ) )
>               ccone  = tmp1 / tmp2
>           END IF
>C 
>C*          Compute lat/lon of the lower-left corner. Sign = 1/-1 
>C*         denotes NH/SH.
>C
>           dxp = 1. - rxp
>           dyp = 1. - ryp
>           rm  = dx * SQRT ( dxp * dxp + dyp * dyp )
>cdo        tmp = ccone / ( RADIUS * SIN ( psi1 ) ) 
>           tmp = ccone / ( re * SIN ( psi1 ) ) 
>           arg = ( rm * tmp ) ** ( 1. / ccone ) * TAN ( psi1 / 2. )
>           dlat1 = sign * ( HALFPI - 2. * ATAN ( arg ) ) * RTD
>C
>           IF ( dyp .ne. 0. ) THEN
>               dyp   = - dyp
>               thta  = ATAN2 ( dxp, dyp ) * RTD / ccone
>               dlon1 = clon + thta
>               CALL PRNLON ( 1, dlon1, ier ) 
>              ELSE
>               dlon1 = clon
>           END IF
>C
>C*         Compute lat/lon of the upper-right corner.
>C
>           dxp = rgx - rxp
>            dyp = rgy - ryp
>           rm  = dx * SQRT ( dxp * dxp + dyp * dyp )
>           arg = ( rm * tmp ) ** ( 1. / ccone ) * TAN ( psi1 / 2. )
>           dlat2 = sign * ( HALFPI - 2. * ATAN ( arg ) ) * RTD
>           IF ( dyp .ne. 0. ) THEN
>               dyp   = - dyp
>               thta  = ATAN2 ( dxp, dyp ) * RTD / ccone
>               dlon2 = clon + thta
>               CALL PRNLON ( 1, dlon2, ier )
>             ELSE
>               dlon2 = clon
>           END IF
>C
>C*         Northern or Southern Hemisphere projection
>C
>           IF ( ( std1 + std2 ) .gt. 0. ) THEN
>               proj = 'LCC'
>             ELSE
>               proj = 'SCC'
>           END IF
>           dangl1 = std1
>           dangl2 = clon
>           dangl3 = std2
>C
>C
>C*       *******************
>C*       *** Rectilinear ***
>C*       *******************
>C
>         ELSE IF ( navtyp .eq. 'RECT' ) THEN
>C
>C*         Row/column and lat/lon of reference point.
>C
>           rrow = ianav ( 2 )
>           rcol = ianav ( 4 )
>           dlat = ianav ( 3 ) / 10000.
>           dlon = ianav ( 5 ) / 10000.
>C
>C*         Determine longitude convention.
>C
>           IF ( ianav ( 11 ) .ge. 0 ) sign = - 1.
>           dlon = sign * dlon
>C
>C*         Lat/lon resolution (degrees/imageline).
>C*         Check for scaling factor. If it is not present use the
>C*         default scaling.
>C
>           IF  ( ianav ( 14 ) .eq. 0 )  THEN
>               reslat = ianav ( 6 ) / 10000.
>             ELSE
>               reslat = ianav ( 6 ) / ( 10. ** ianav (14) )
>           END IF
>           IF  ( ianav ( 15 ) .eq. 0 )  THEN
>               reslon = ianav ( 7 ) / 10000.
>             ELSE
>               reslon = ianav ( 7 ) / ( 10. ** ianav (15) )
>           END IF
>C
>C*         Determine upper-left lat/lon
>C*         (McIDAS image full coordinates 1,1).
>C
>           uldlat = dlat + reslat * ( rrow - iullne )
>           uldlon = dlon + reslon * sign * ( rcol - iulele )
>C
>C*         Lat/lon for lower-left corner.
>C
>           dlat1 = uldlat - ( reslat * (imnlin * rmyres - 1))
>           dlon1 = uldlon
>C
>C*         Lat/lon for upper-right corner.
>C
>           dlat2 = uldlat
>           dlon2 = uldlon - sign * ( reslon * (imnpix * rmxres - 1) )
>C
>C*         Set angle1 to the lon/lat aspect ratio.
>C*         Set angle2 to the central longitude.
>C
>           dangl1 = reslon / reslat
>           dangl2 = uldlon -
>     +              sign * ( reslon * ( (imnpix * rmxres) / 2 - 1 ) )
>           dangl3 = 0.
>           proj   = 'MCD' 
>C
>C*       Not valid
>C
>         ELSE
>             iret = -5
>             RETURN
>       END IF
>C
>C*      Set navigation blocks.
>C
>       CALL GR_MNAV  ( proj, imnpix, imnlin, dlat1, dlon1, dlat2,
>     +                  dlon2, dangl1, dangl2, dangl3, .true.,
>     +                  rnvblk, ier )
>C
>C*     Set map projection.
>C
>       CALL GSMPRJ ( proj, dangl1, dangl2, dangl3,
>     +         dlat1, dlon1, dlat2, dlon2, ier )
>        IF ( ier .ne. 0 ) iret = -5
>C
>       RETURN
>       END
>