program magmap c======================================================================= c UNIT 0: MAIN PROGRAM c======================================================================= c Driver program for making maps of magnetic fields or evaluating c field elements at specified points based upon spherical harmonic c models. Runs under the command language paradigm. c The program is built from 7 units: c UNIT 0: Main program c UNIT M: Mapping and plotting c UNIT I: Input commands c UNIT S: Spherical harmonic input and normalization c UNIT K: Decoding routines for Command Line Interface c UNIT Y: Generate and sum Ylm c UNIT P: Single point evaluation c parameter (inmx=200) parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /ndict/ iecho,nin,memory,istate(inmx) common /prolix/ iprint c ce 1000 call inpsh call getdat call point call mapper goto 1000 end c_______________________________________________________________________ subroutine projg(iproj, geo, u) c$$$$ calls nothing c Applies the forwd map projection from geographic coordinates geo c to the point u in R**2 on the map c dimension u(2), geo(2), x(3) common /orient/ center(2),O(3,3) data deg/0.01745329/,degby2/0.008726646/,pi/3.1415926/ ce if (geo(1) .gt. 90.0) then u(1)=999.0 return endif c cendg=center(2) c goto (1000, 2000, 3000, 4000, 5000), iproj c c Rectilinear projection: u(1)=longitude, u(2)=latitude c 0 <= u(1) <= 360; -90 <= u(2) <= +90 c 1000 u(2)=geo(1) u(1)=mod(geo(2) - cendg + 540.0, 360.0) - 180.0 return c c Lambert equal area: before rotations (0,0) in u plane is (0,0,1) in c R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 4 is mapped from the (punctured) sphere c Map onto unit sphere 2000 x(3)=sin(deg*geo(1)) rho =cos(deg*geo(1)) x(1)=rho*cos(deg*geo(2)) x(2)=rho*sin(deg*geo(2)) c Apply coordinate rotation O**T x if (O(3,3) .ge. 2.0) then x1=x(1) x2=x(2) x3=x(3) else x1 =O(1,1)*x(1) + O(2,1)*x(2) + O(3,1)*x(3) x2 =O(1,2)*x(1) + O(2,2)*x(2) + O(3,2)*x(3) x3 =O(1,3)*x(1) + O(2,3)*x(2) + O(3,3)*x(3) endif c c Map into Lambert plane: rho = 2*sin(theta/2) rho=sqrt(2.0*(1.0 -x3)/(x1**2 + x2**2)) u(1)= rho*x2 u(2)=-rho*x1 return c c Hammer-Aitoff equal-area projection: unit sphere is mapped on c the region u(1)**2/4 + u(2)**2 < 1 3000 clat=cos(deg*geo(1)) half=degby2*(mod(geo(2) - cendg + 540.0, 360.0) - 180.0) R=sqrt(1.0/(abs(1.0 + clat*cos(half)) + 1.0e-15)) u(1)=2.0*R*clat*sin(half) u(2)=R*sin(deg*geo(1)) return c c Orthogonal projection: before rotations (0,0) in u plane is (0,0,1) c in R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 1 is mapped into the (punctured) sphere. c Map onto unit sphere 4000 x(3)=sin(deg*geo(1)) rho =cos(deg*geo(1)) x(1)=rho*cos(deg*geo(2)) x(2)=rho*sin(deg*geo(2)) c Apply coordinate rotation O**T x if (O(3,3) .ge. 2.0) then x1=x(1) x2=x(2) x3=x(3) else x1 =O(1,1)*x(1) + O(2,1)*x(2) + O(3,1)*x(3) x2 =O(1,2)*x(1) + O(2,2)*x(2) + O(3,2)*x(3) x3 =O(1,3)*x(1) + O(2,3)*x(2) + O(3,3)*x(3) endif c Project upper hemisphere onto equatorial plane: rho = sin(theta) if (x3 .ge. 0.0) then u(1)= x2 u(2)=-x1 else u(1)=2.0 u(2)=2.0 endif return c c Mercator projection: u(1)=longitude (radians), u(2)=M c -pi <= u(1) <= pi; -pi/2 <= u(2) <= +pi/2 c 5000 u(2)=log(abs(tan(degby2*(geo(1)+90.0)))) u(1)=(mod(geo(2) - cendg + 540.0, 360.0) - 180.0)*deg if (abs(u(2)) .gt. pi) u(1)=999.0 return c end c_______________________________________________________________________ blockdata coast c Loads common /world/ with a coastline data set. Each real word c contains a lat-long pair recovered by: c lat=c(j)/10; long=1000*mod(abs(c(j),1.0) c Breaks signaled by lat=99 c common /world/ c(3217),nworld data nworld/3217/ c data (c(j),j=1,130)/94.28,90.2823,84.283,90.2837,95.2842, $100.284,110.285,115.287,119.288,110.288,121.29,115.291, $110.292,104.292,100.298,94.2989,86.2992,79.301,70.3016, $63.3026,59.3032,56.3064,51.3074,46.3081,42.3086,38.3092, $30.3093,20.3096,16.3102,8.31,0.3093,-9.3083,0.3091, $-5.3118,-10.314,-14.315,-20.3156,-27.316,-23.317,-29.318, $-33.3209,-40.3219,-46.3227,-51.3235,-59.325,-70.3252, $-80.3251,-89.3249,-99.3242,-107.323,-110.323,-114.323, $-120.322,-127.322,-132.321,-140.321,-150.321,-160.321, $-170.321,-177.321,-183.32,-191.321,-200.32,-210.319, $-220.319,-230.318,-234.315,-237.314,-245.313,-250.312, $-254.312,-263.311,-270.312,-280.311,-286.311,-292.31, $-300.31,-309.309,-317.309,-321.308,-328.308,-334.307, $-340.307,-347.306,-342.302,-349.302,-360.303,-370.303, $-376.303,-380.303,-384.302,-389.3,-394.298,-402.298, $-406.298,-411.297,-407.295,-412.295,-420.295,-427.296, $-436.295,-443.295,-449.295,-454.293,-460.292,-466.293, $-470.293,-476.294,-481.294,-484.293,-490.293,-496.292, $-501.292,990.37,-521.3,-516.299,-512.299,-516.302,-520.301, $-521.3,990.37,-501.292,-504.291,-510.291,-516.291,-522.292, $-530.292,-535.292,-540.293,-545.294/ data (c(j),j=131,260)/-550.294,-554.291,-551.29,-544.289, $-537.288,-530.287,-520.287,-512.286,-505.286,-497.286, $-489.286,-480.286,-470.286,-466.286,-460.286,-453.287, $-446.287,-440.287,-431.287,-421.288,-425.287,-433.286, $-427.286,-420.286,-413.286,-401.286,-391.287,-382.286, $-372.286,-360.287,-350.288,-340.288,-330.288,-320.288, $-310.288,-302.288,-296.289,-289.288,-283.289,-273.289, $-261.289,-250.289,-240.289,-230.289,-220.29,-210.29, $-200.29,-190.29,-180.289,-173.289,-165.288,-160.287, $-156.286,-152.285,-143.284,-135.284,-126.283,-117.283, $-109.282,-100.282,-90.2813,-80.2809,-70.2803,-64.2799, $-59.2788,-50.2788,-41.2787,-35.2795,-29.2801,-22.2791, $-10.2791,0.2798,10.2805,16.2812,27.2817,37.283,46.2828, $56.2827,66.2828,73.2822,83.2818,89.2813,84.28,990.37, $310.033,320.035,330.035,340.036,347.036,356.036,367.036, $361.034,360.033,990.37,350.033,345.033,350.034,355.035, $350.033,990.37,360.033,366.032,361.031,365.029,374.028, $380.027,388.027,394.026,402.027,410.03,420.034,405.038, $405.039,420.041,435.039,450.038,454.037,472.038,464.034, $452.036,450.034,457.033,461.033,466.031,464.03,435.029, $410.029,407.027,401.024,405.023/ data (c(j),j=261,390)/396.023,389.024,385.024,379.025, $372.024,373.023,990.37,450.047,464.049,479.05,479.052, $457.053,450.051,420.052,413.052,417.054,383.053,364.054, $361.052,376.049,398.05,435.047,450.047,990.37,356.024, $350.025,354.025,356.024,990.37,373.023,364.023,369.022, $366.022,370.022,379.021,383.021,390.021,399.02,406.02, $417.02,423.019,429.018,434.017,438.016,446.015,450.015, $456.014,450.013,443.012,439.013,435.014,429.014,423.015, $419.015,416.016,412.017,408.018,402.019,406.017,399.017, $395.017,390.017,383.017,379.016,387.016,990.37,381.016, $375.013,369.015,366.015,370.015,376.015,381.016,990.37, $387.016,397.016,405.015,410.014,414.013,420.012,427.011, $430.011,990.37,429.01,425.009,418.009,414.009,422.01, $429.01,990.37,412.009,407.01,400.01,391.01,400.009, $406.008,409.008,412.009,990.37,400.003,393.003,398.003, $400.003,990.37,430.011,440.01,444.009,438.008,431.007, $434.005,430.003,421.003,417.003,413.002,401,396.36,389, $383.36,377.359,369.358,364.355,361.355,990.37,104.107, $96.1061,90.1053,86.1049,94.1046,101.105,107.104/ data (c(j),j=391,520)/110.103,121.103,128.102,135.101, $120.1,110.099,100.099,92.0993,83.1001,74.1005,69.1009, $63.1024,53.103,43.1034,33.1033,23.1039,14.1041,20.1027, $25.1019,35.1012,43.1006,59.1004,70.0999,78.0993,82.0984, $91.0983,100.099,111.099,123.099,133.099,144.098,153.098, $163.098,170.097,163.096,990.37,135.093,125.093,117.093, $126.092,135.093,990.37,163.096,158.095,171.094,183.094, $190.094,196.093,202.093,211.092,220.092,230.091,220.09, $213.087,208.087,200.086,190.084,183.084,176.083,170.082, $164.082,159.081,153.08,144.08,134.08,126.08,117.08,106.08, $96.0789,90.078,82.0774,990.37,100.08,97.0805,89.081, $77.0816,69.0817,63.081,73.0796,87.0797,100.08,990.37, $82.0774,90.0764,102.076,113.076,122.075,133.075,144.074, $152.074,160.073,170.073,180.073,190.073,200.072,210.073, $217.072,224.072,214.072,210.071,220.069,225.069,230.07, $238.068,241.067,249.067,255.067,251.062,255.059,263.057, $270.057,264.055,271.053,277.052,286.051,294.051,301.05, $307.049,300.049,294.048,286.049,279.049,270.05,263.05, $256.05,253.051,261.051,256.052,250.052,254.052/ data (c(j),j=521,650)/241.052,250.055,257.056,263.057, $252.057,243.057,238.058,228.059,224.06,214.059,204.059, $190.058,180.057,170.055,163.053,156.053,152.052,149.051, $140.049,133.047,129.045,134.043,143.043,153.043,164.043, $173.042,180.042,192.041,201.04,208.039,217.039,225.039, $235.039,241.038,251.037,260.037,268.036,276.036,282.035, $288.033,300.033,990.37,300.033,290.033,280.034,270.034, $260.035,250.035,240.036,229.036,222.037,210.037,200.037, $189.038,183.038,174.039,163.039,156.04,150.041,146.041, $139.042,130.043,121.044,116.043,107.044,104.045,108.046, $111.048,116.05,120.051,109.051,100.051,90.0509,80.0502, $70.0497,60.0494,50.0486,40.0478,31.047,23.046,17.045, $10.044,0.0432,-10.0422,-20.0413,-26.0406,-37.04,-47.0395, $-55.0392,-62.039,-70.0398,-80.0397,-90.0398,-100.04, $-107.041,-118.041,-129.041,-140.041,-149.041,-156.041, $-163.04,990.37,-120.049,-128.05,-140.05,-150.05,-160.05, $-156.05,-165.05,-173.05,-184.049,-194.049,-206.049, $-215.048,-225.048,-235.048,-244.047,-252.047,-256.045, $-253.045,-248.044,-239.044,-229.043,-219.043,-212.044, $-206.044,-198.045,-187.044,-174.044,-163.045/ data (c(j),j=651,780)/ -156.047,-150.048,-136.048,-132.049, $-120.049,990.37,-163.04,-171.039,-179.037,-187.036, $-194.036,-200.035,-206.035,-214.035,-223.036,-230.036, $-240.036,-247.035,-254.033,-260.033,-270.033,-280.033, $-290.032,-300.031,-310.031,-316.03,-324.029,-332.028, $-337.027,-342.025,-349.02,-344.019,-338.018,-328.018, $-324.018,-315.018,-307.017,-297.017,-288.017,-281.016, $-272.015,-260.015,-250.015,-241.014,-230.014,-223.015, $-215.014,-206.013,-197.013,-188.012,-178.012,-167.012, $-154.012,-144.012,-135.013,-127.013,-124.014,-114.014, $-104.014,-98.0132,-90.0129,-82.0133,-70.0129,-61.0123, $-50.012,-40.0112,-30.0103,-20.0095,-8.0089,-1.0094,10.0098, $20.0101,31.0103,40.0097,46.009,50.0056,60.0051,64.004, $59.001,56,52.359,48.358,52.356,45.3529,51.3507,59.35, $66.349,71.348,80.3472,90.3469,100.346,110.345,120.345, $124.343,134.343,144.343,150.343,157.344,167.344,176.344, $188.344,196.344,203.344,211.343,221.343,230.344,240.344, $249.345,260.346,266.346,273.347,280.347,288.349,291.35, $296.35,302.351,308.35,316.35,321.351,329.351,333.351, $341.353,350.354,355.354,359.354,355.355,360.359,366.001, $369.003,373.01/ data (c(j),j=781,910)/ 366.011,362.011,356.011,350.011, $346.011,342.01,336.011,332.011,327.014,323.015,315.016, $310.018,306.019,303.019,310.02,320.02,327.021,323.023, $316.025,311.028,310.033,990.37,687.049,691.049,695.049, $692.05,688.05,687.049,990.37,800.048,804.047,809.051, $806.052,802.049,800.048,990.37,806.057,811.055,815.057, $811.059,807.058,806.057,990.37,803.062,808.06,803.062, $990.37,810.066,807.064,812.063,810.066,990.37,704.057, $708.054,712.054,721.052,726.053,731.053,737.054,742.055, $751.056,754.058,758.059,762.062,765.066,768.067,762.068, $758.066,755.063,751.062,745.059,740.058,734.057,730.056, $723.057,714.056,709.057,704.057,990.37,680.05,683.051, $689.054,683.054,680.053,675.053,678.054,682.055,685.057, $688.059,683.059,687.061,693.06,698.06,703.059,698.061, $693.064,690.066,686.067,681.068,687.069,690.068,693.068, $700.067,705.067,710.067,713.068,720.069,729.07,734.07, $730.07,726.073,717.072,709.073,699.073,689.073,684.074, $676.073,670.072,667.071,660.072,667.073,672.074,680.075, $685.074,688.075,700.074,706.074,713.073,719.073,725.075/ data (c(j),j=911,1040)/ 730.074,723.076,713.075,710.078, $713.077,719.076,723.077,720.081,717.081,723.082,730.081, $736.08,739.087,743.086,749.087,755.089,758.093,762.099, $764.1,990.37,788.098,794.094,798.094,802.092,805.092, $808.093,812.096,808.097,805.097,802.097,798.1,792.1, $788.099,788.098,990.37,779.099,784.1,789.101,793.102, $790.104,784.105,781.101,779.099,990.37,764.1,771.101, $774.102,771.104,767.107,763.113,759.114,753.114,749.113, $746.112,741.11,736.108,729.106,736.11,740.11,736.112, $741.113,745.112,741.113,736.113,732.119,728.123,732.123, $736.123,732.128,728.13,724.13,721.129,716.129,711.13, $707.131,715.132,719.133,715.133,719.14,724.14,990.37, $746.139,751.137,757.137,760.138,756.141,750.14,746.139, $990.37,750.144,753.142,758.143,767.142,759.143,756.145, $752.145,750.144,990.37,747.15,750.147,754.147,751.151, $747.15,990.37,732.144,736.141,739.142,734.144,732.144, $990.37,724.14,728.141,725.145,722.149,716.15,710.152, $707.159,703.16,698.159,694.161,699.168,694.169,691.168, $687.17,990.37,716.24,712.237,715.236,719.235,725.235/ data (c(j),j=1041,1170)/ 729.235,735.235,738.236,744.235, $741.243,737.244,733.243,729.242,726.241,719.239,716.24, $990.37,693.244,700.243,703.245,710.241,713.243,722.241, $728.242,732.244,726.245,730.247,727.249,723.248,727.25, $720.251,716.252,724.252,731.252,736.253,732.255,728.255, $724.255,719.255,715.256,710.255,706.256,701.259,696.259, $692.258,686.257,690.254,693.253,688.253,685.25,691.246, $693.244,990.37,688.243,683.246,679.245,676.247,679.25, $676.251,669.253,678.252,682.253,686.252,682.255,677.258, $683.261,687.261,693.262,698.262,695.263,690.264,686.265, $680.264,672.264,680.265,686.266,693.266,698.264,702.263, $706.264,713.264,718.265,724.265,733.264,740.265,732.269, $726.268,719.266,716.266,707.267,703.268,696.268,689.269, $682.27,990.37,717.264,712.262,718.26,722.259,725.258, $729.258,725.259,732.26,736.259,739.261,736.263,731.262, $726.263,717.264,990.37,746.265,750.263,756.265,752.267, $746.266,746.265,990.37,750.262,756.26,760.256,766.255, $758.258,764.258,757.262,750.262,990.37,750.254,746.249, $749.248,757.243,762.244,755.249,758.25,764.25,758.252/ data (c(j),j=1171,1300)/ 755.255,750.254,990.37,760.241, $757.24,761.237,764.238,770.24,773.241,770.244,765.244, $762.242,767.242,760.241,990.37,778.247,774.25,778.247, $990.37,782.25,785.249,782.25,990.37,785.255,789.256, $793.254,788.257,785.26,781.261,784.257,785.255,990.37, $780.262,788.262,785.265,782.265,778.264,780.262,990.37, $781.271,784.267,789.266,794.268,800.264,810.265,813.265, $808.269,805.27,797.273,793.275,788.273,781.271,990.37, $765.27,769.271,773.273,778.272,774.274,781.275,787.273, $794.275,798.273,803.274,807.281,814.27,818.269,822.273, $826.277,829.28,826.291,829.295,825.297,821.298,817.296, $813.295,809.293,805.291,800.289,796.288,792.284,780.284, $774.282,770.281,765.282,762.28,765.27,990.37,766.269, $762.271,756.271,750.28,745.28,754.268,763.267,767.263, $771.264,765.267,766.269,990.37,721.27,733.271,737.274, $732.275,735.276,728.284,725.284,720.286,716.286,708.29, $703.292,697.293,692.292,684.292,680.293,673.296,670.297, $665.299,660.298,654.297,647.296,653.295,659.294,664.293, $656.292,650.293,647.294,640.295,624.295,629.294,634.292/ data (c(j),j=1301,1430)/ 630.292,626.293,618.294,622.291, $627.29,636.288,640.287,645.286,641.283,648.282,653.283, $661.286,664.286,671.288,680.287,684.286,688.285,693.284, $702.282,698.281,702.273,707.271,710.272,721.27,990.37, $683.284,675.285,671.284,677.283,683.284,990.37,658.274, $654.275,650.277,647.277,639.279,636.28,640.277,635.276, $631.275,636.274,641.274,645.274,653.274,658.274,990.37, $-782.183,-778.2,-772.202,-768.209,-765.21,-762.211, $-756.216,-750.221,-745.224,-742.229,-738.236,-731.237, $-734.239,-738.247,-744.252,-738.253,-746.255,-750.259, $-743.259,-740.259,-735.259,-732.26,-728.258,-724.258, $-721.258,-717.26,-720.264,-726.265,-729.266,-726.268, $-731.271,-737.278,-731.28,-726.282,-729.285,-726.288, $-720.285,-717.284,-711.283,-714.285,-710.285,-706.286, $-700.285,-694.284,-691.287,-686.289,-689.29,-694.29, $-690.292,-683.293,-677.292,-670.291,-667.292,-670.293, $-666.294,-662.294,-654.296,-649.296,-646.298,-640.299, $-636.301,-630.302,-635.303,-639.301,-644.301,-648.301, $-654.301,-658.3,-663.3,-668.299,-674.299,-679.299,-684.299, $-689.299,-694.299,-699.299,-703.299,-712.3,-716.3,-721.3, $-725.3,-730.3,-734.3,-741.299/ data (c(j),j=1431,1560)/ -746.299,-751.299,-754.3,-761.303, $-764.305,-769.308,-773.311,-777.312,-780.315,-776.318, $-780.32,-777.324,-772.326,-766.33,-762.332,-757.333, $-751.335,-744.336,-740.338,-737.339,-730.34,-727.341, $-723.344,-720.345,-714.348,-710.349,-704.351,-700.358, $-690.36,-696,-699.001,-695.013,-691.015,-696.016,-699.02, $-696.029,-691.033,-686.034,-689.035,-695.037,-699.038, $-695.04,-689.04,-682.041,-677.044,-673.046,-670.049, $-663.05,-660.051,-657.053,-660.055,-665.057,-670.057, $-673.059,-677.066,-688.07,-684.072,-689.074,-693.077, $-689.078,-685.079,-681.079,-675.081,-671.082,-663.081, $-659.082,-664.083,-661.087,-665.089,-668.089,-665.092, $-661.095,-653.096,-650.096,-654.098,-658.099,-653.101, $-659.105,-664.108,-669.109,-664.11,-660.111,-657.112, $-660.115,-664.115,-668.117,-659.121,-663.123,-660.123, $-666.127,-670.128,-656.131,-659.131,-663.138,-667.14, $-671.143,-676.146,-679.147,-683.148,-690.156,-693.159, $-700.161,-704.161,-709.164,-706.167,-710.168,-717.17, $-723.171,-729.17,-734.169,-740.168,-743.167,-748.167, $-756.164,-766.164,-775.164,-783.165,-780.167,-776.167, $-779.178,990.37,682.27,691.271,687.272,676.272,671.273, $676.274,681.273,686.274,691.275/ data (c(j),j=1561,1690)/ 697.275,691.278,684.279,677.278, $667.278,662.276,665.275,662.274,654.273,650.273,641.272, $635.27,630.269,620.267,610.266,600.265,590.265,586.266, $577.267,569.268,565.272,559.273,553.275,540.278,529.278, $519.279,510.28,516.281,522.282,530.281,538.281,545.28, $553.283,558.283,566.284,574.283,580.282,584.281,591.282, $597.283,603.282,607.282,614.282,623.282,616.288,610.289, $600.29,591.291,587.291,581.292,584.293,589.294,598.295, $603.295,596.296,588.297,577.298,568.299,561.298,556.299, $551.3,547.302,539.302,535.303,524.304,518.304,513.303, $513.302,990.37,509.303,514.304,505.304,496.303,499.304, $495.304,484.306,476.306,480.307,472.307,465.307,476.306, $468.305,474.305,479.301,483.302,490.302,501.303,509.303, $990.37,498.296,494.298,490.299,498.296,990.37,513.302, $507.301,500.3,496.293,491.292,485.291,480.29,474.29, $480.291,485.292,489.293,479.295,469.295,460.296,456.297, $460.299,469.3,462.3,458.3,454.3,448.298,443.296,436.295, $443.294,449.295,456.296,450.294,444.292,437.291,428.289, $420.289,416.29,408.287,398.286,390.285,394.285/ data (c(j),j=1691,1820)/ 384.285,377.284,370.284,379.284, $383.284,374.284,370.284,366.284,358.284,351.284,346.283, $339.282,334.281,325.28,320.279,310.279,300.279,290.279, $284.28,274.28,266.28,258.28,250.28,255.279,266.278, $278.277,285.278,292.277,299.276,295.275,300.274,290.271, $297.268,294.268,290.265,286.264,282.263,270.263,260.263, $250.262,240.262,230.262,220.262,215.263,210.263,200.263, $190.264,182.266,187.268,196.27,209.27,212.271,206.273, $197.273,181.272,186.272,173.272,165.272,158.271,154.276, $150.277,140.277,130.277,120.277,111.276,100.277,92.278, $88.279,94.28,990.37,219.275,226.276,230.278,225.28, $217.282,212.284,206.285,200.286,205.283,212.282,219.28, $224.279,220.277,216.276,219.275,990.37,196.287,192.291, $183.292,174.288,180.288,184.286,192.287,196.287,990.37, $84.28,75.2803,78.279,82.2786,88.2764,92.2766,97.276, $100.275,110.274,116.274,122.273,127.273,132.272,139.27, $146.268,152.267,157.267,161.266,158.264,161.262,167.261, $173.259,176.259,181.258,187.256,193.256,197.255,205.254, $214.255,220.254,228.254,234.254,240.253,247.252,254.252/ data (c(j),j=1821,1950)/ 261.251,265.251,269.25,274.25, $281.249,286.248,293.248,300.247,310.247,314.247,317.246, $314.245,307.245,300.246,294.246,289.247,283.247,277.248, $271.248,266.248,260.249,250.249,245.249,236.25,231.25, $236.25,242.249,247.248,256.248,261.248,266.247,270.247, $280.245,287.246,293.245,296.245,305.244,313.244,321.243, $328.243,334.242,342.241,352.239,360.238,369.238,380.237, $385.237,390.236,398.236,406.236,414.236,422.236,428.235, $438.236,448.236,458.236,467.236,475.236,483.235,490.235, $494.234,499.233,507.232,513.232,519.232,527.231,534.23, $542.23,549.229,990.37,189.204,200.204,193.205,187.204, $189.204,990.37,540.228,531.227,525.228,520.229,528.228, $535.228,540.228,990.37,549.229,554.227,562.227,568.226, $576.225,582.224,586.222,592.221,597.22,602.216,607.214, $599.211,596.21,600.208,609.209,613.209,609.208,603.208, $597.207,590.206,589.207,990.37,584.207,580.207,575.205, $567.206,573.207,579.207,583.208,584.207,990.37,589.207, $582.206,578.205,574.204,570.203,566.202,559.201,556.2, $549.197,545.196,551.197,554.198,558.198,563.2,567.201/ data (c(j),j=1951,2080)/ 575.202,580.202,589.203,586.202, $589.202,584.201,590.2,585.198,593.198,602.198,596.197, $598.196,990.37,600.194,600.194,990.37,598.196,603.195, $609.195,614.194,618.194,624.195,629.195,633.196,638.199, $645.199,648.199,643.197,651.193,654.193,659.193,663.194, $660.198,663.199,667.197,671.197,675.196,678.195,683.193, $687.194,690.196,695.197,700.197,703.198,708.201,713.203, $708.204,712.205,708.205,705.208,702.213,698.218,694.221, $691.222,696.226,699.229,703.23,697.231,701.232,706.232, $701.233,697.233,701.235,696.236,692.241,688.243,990.37, $687.17,690.171,700.17,697.174,694.178,990.37,708.179, $713.179,708.179,990.37,694.178,689.18,990.37,686.181, $680.183,676.184,666.185,670.186,665.189,660.19,656.189, $649.188,646.187,641.187,644.186,648.184,653.184,662.181, $656.181,651.18,990.37,646.179,649.177,645.176,642.178, $631.179,626.18,622.176,617.174,610.172,605.171,599.17, $604.169,598.166,604.166,598.165,589.163,584.162,577.162, $573.163,566.163,559.163,563.163,554.162,547.162,540.16, $531.16,522.158,515.158,509.157,520.156,530.156,540.156/ data (c(j),j=2081,2210)/ 550.155,560.156,569.157,575.157, $579.157,584.159,589.16,595.16,599.161,604.162,608.164, $613.164,622.164,626.164,614.162,609.161,606.16,613.16, $619.16,616.159,619.159,614.157,607.157,603.155,595.154, $590.154,595.151,591.146,594.146,590.142,584.141,577.14, $571.139,563.138,557.137,551.136,547.135,538.137,542.138, $536.138,542.139,538.14,533.141,522.141,516.141,510.141, $500.141,490.14,990.37,490.142,500.142,511.142,521.142, $533.142,541.142,534.143,523.143,512.144,501.144,488.145, $493.144,484.143,476.143,470.143,463.144,467.143,460.142, $467.142,475.142,482.142,490.142,990.37,490.14,480.14, $470.139,460.138,452.137,443.136,435.135,428.134,433.132, $427.131,420.13,410.13,405.129,400.128,990.37,452.142, $448.143,443.143,436.145,429.145,421.143,426.142,418.141, $422.14,428.14,433.14,437.141,444.142,452.142,990.37, $415.141,407.142,400.142,390.142,384.142,381.141,373.141, $366.141,358.141,350.14,357.14,353.139,347.139,352.139, $347.138,343.137,339.136,334.136,340.135,347.135,339.132, $335.132,328.132,320.132,311.131,320.13,326.131,331.13/ data (c(j),j=2211,2340)/ 337.13,343.131,349.132,356.133, $360.136,369.137,374.137,367.137,371.138,378.139,385.14, $394.14,404.14,411.14,415.141,990.37,340.133,344.134, $338.135,333.134,327.133,334.132,340.133,990.37,400.128, $391.128,384.129,377.129,370.13,360.13,352.129,344.127, $353.126,359.127,368.126,378.126,388.125,394.125,398.125, $393.122,388.122,392.122,401.122,406.122,410.122,400.12, $392.119,386.118,380.118,374.119,370.12,378.121,374.122, $369.122,365.121,360.12,356.12,350.119,342.12,330.121, $319.122,309.122,304.121,300.122,290.122,281.121,270.12, $259.12,252.119,243.118,236.117,230.116,226.114,219.113, $214.111,205.111,210.11,990.37,200.11,209.111,193.111, $186.11,183.11,187.109,194.109,200.109,200.11,990.37, $210.11,217.109,210.107,200.106,190.106,180.106,173.107, $164.108,152.109,140.109,130.109,120.109,111.109,108.108, $104.107,990.37,220.121,229.12,237.12,243.12,251.121, $241.122,230.121,220.121,990.37,186.121,172.123,162.122, $154.121,141.122,129.124,134.123,139.123,131.123,139.122, $135.121,126.122,123.121,135.12,146.12,153.12,162.12/ data (c(j),j=2341,2470)/ 170.12,180.12,186.121,990.37, $125.124,113.126,104.125,100.125,96.124,92.123,98.1224, $106.123,119.122,115.123,110.123,119.124,125.124,990.37, $98.1255,89.1264,79.1267,70.1266,64.1263,71.126,62.126, $56.1255,60.1247,64.1241,70.1241,73.1244,79.1236,75.1235, $78.1227,69.122,80.1221,86.123,82.1242,86.1247,90.1249, $98.1255,990.37,-13.1307,-8.131,-4.1321,-8.134,-20.134, $-30.1349,-26.1363,-20.1371,-14.1378,-18.1387,-22.1398, $-26.1414,-30.1425,-34.1434,-38.1445,-43.1453,-48.146, $-54.1458,-60.1476,-64.148,-69.1472,-76.1477,-85.1484, $-90.149,-95.1496,-100.15,-106.15,-102.15,-94.147,-81.1461, $-77.1451,-83.1431,-90.1434,-94.1426,-90.141,-81.14, $-74.1383,-62.1385,-55.1382,-49.1373,-46.1356,-41.1347, $-34.1336,-40.1332,-34.1327,-28.1318,-24.133,-20.1337, $-16.1316,-13.1307,990.37,-55.1483,-50.1511,-41.1515, $-44.1523,-56.1519,-61.151,-55.1483,990.37,-28.1508, $-33.1522,-41.1532,-49.1529,-39.1524,-35.152,-28.1508, $990.37,-55.1548,-60.1554,-68.1558,-63.155,-55.1548,990.37, $-100.16,-94.1596,-100.161,-100.16,990.37,-203.164,-207.165, $-216.166,-223.167,-214.165,-203.164,990.37,-123.137, $-131.136,-143.136,-149.136/ data (c(j),j=2471,2600)/ -154.136,-160.137,-166.138,-170.139, $-174.14,-170.141,-160.141,-150.142,-140.142,-130.142, $-120.142,-110.142,-107.143,-117.143,-126.143,-133.144, $-143.144,-150.145,-160.145,-170.146,-180.146,-190.146, $-194.147,-200.148,-204.149,-211.149,-220.149,-224.151, $-233.151,-240.151,-247.152,-256.153,-264.153,-271.153, $-281.154,-290.154,-300.153,-310.153,-320.152,-328.152, $-333.151,-341.151,-350.151,-360.15,-369.15,-376.15, $-379.149,-386.147,-389.146,-384.145,-380.145,-384.144, $-388.144,-384.143,-380.14,-375.14,-370.14,-360.14,-354.139, $-348.139,-341.138,-351.138,-340.137,-334.138,-326.138, $-334.137,-340.136,-350.136,-337.135,-330.134,-323.134, $-315.131,-320.128,-326.125,-330.124,-337.124,-344.119, $-350.118,-341.115,-335.115,-331.116,-321.116,-314.116, $-303.115,-291.115,-280.114,-270.114,-260.113,-255.113, $-261.114,-250.114,-240.113,-234.114,-225.113,-219.114, $-223.114,-216.115,-210.116,-206.117,-202.118,-199.119, $-190.121,-180.122,-171.122,-163.123,-171.123,-160.123, $-153.125,-145.125,-139.126,-148.128,-138.13,-133.13, $-127.13,-124.131,-116.132,-112.132,-119.133,-125.136, $-120.136,-123.137,990.37,-408.145,-412.147,-419.148, $-428.148,-436.147,-429.146,-421.145,-415.145/ data (c(j),j=2601,2730)/ -408.145,990.37,-179.177,-172.178, $-179.179,-179.177,990.37,-345.173,-351.174,-359.175, $-369.175,-365.175,-376.176,-379.177,-374.178,-386.178, $-390.178,-397.177,-406.176,-411.176,-416.175,-407.175, $-401.175,-393.174,-380.175,-374.175,-367.174,-360.174, $-354.173,-345.173,990.37,-407.173,-412.173,-418.174, $-425.174,-432.173,-439.173,-450.171,-459.171,-465.169, $-461.167,-450.167,-440.168,-431.17,-426.171,-419.171, $-410.172,-407.173,990.37,0.128,10.1277,20.1282,10.1283, $6.1287,2.129,-9.1283,0.128,990.37,-32.128,-38.1308, $-32.128,990.37,-38.1266,-31.126,-36.1272,-38.1266,990.37, $-69.1341,-60.134,-55.1345,-62.1347,-69.1341,990.37, $-103.124,-93.124,-87.1252,-84.1271,-90.1261,-94.1254, $-101.125,-103.124,990.37,-97.119,-94.12,-101.121,-97.119, $990.37,-87.12,-83.1206,-80.123,-87.1228,-90.1217,-87.12, $990.37,-90.1159,-81.1161,-86.1191,-90.118,-90.1159,990.37, $-83.1156,-86.1143,-83.1131,-80.1106,-77.1096,-70.1063, $-67.1054,-60.106,-68.1085,-62.1107,-66.1117,-70.1124, $-77.113,-83.1156,990.37,0.1201,-8.1203,-13.1208,-8.1217, $-15.1223,-19.1218,-24.122,-31.1225,-37.1222,-41.123, $-46.1233,-53.1231,-46.1223,-49.122,-40.1216/ data (c(j),j=2731,2860)/ -35.121,-26.1211,-30.1204,-40.1204, $-47.1205,-55.1205,-46.1198,-35.1196,-28.1189,-20.1194, $-11.1195,-7.1199,0.12,10.1203,16.1251,5.1247,0.1201, $990.37,85.1173,94.1181,103.119,112.12,102.119,92.1185, $85.1173,990.37,0.1091,10.109,20.1095,16.1111,23.1113, $30.1118,40.1136,46.114,50.115,60.116,70.1169,65.1176, $60.1182,53.1192,50.1183,44.1186,35.1175,23.118,14.1184, $9.119,0.1176,-10.117,-16.1163,-21.1167,-31.116,-38.1154, $-34.1145,-30.113,-34.1117,-27.1116,-16.1099,-12.11,-4.1094, $0.1091,990.37,0.1035,-6.1032,-10.1042,-19.1047,-15.1054, $-30.1065,-40.1059,-50.1058,-58.1057,-49.1034,-40.1024, $-30.1015,-20.1009,-10.1004,0.0997,10.099,18.0989,23.0977, $31.0974,39.0967,44.0958,56.0951,50.098,40.0988,35.0995, $29.1002,21.1009,16.1019,11.1026,5.1036,0.1035,990.37, $361.355,364.354,370.354,380.351,385.351,394.351,402.351, $410.351,423.351,431.351,434.352,443.359,450.359,460.359, $469.358,474.358,485.355,496.358,493.36,990.37,513.35, $521.35,524.35,531.351,534.35,539.35,545.352,551.352, $546.355,542.355,532.354,526.354,521.354,517.352,513.35, $990.37,500.354/ data (c(j),j=2861,2990)/ 506.355,509.356,514.357,518.355, $528.356,532.356,540.357,544.356,549.357,546.355,550.355, $554.355,558.355,553.355,559.354,564.355,569.354,574.353, $571.354,577.354,585.355,579.356,576.356,571.358,566.358, $562.357,558.358,547.359,543.36,536,531,523.002,519.002, $516.001,512.002,508.001,504.358,500.355,500.354,990.37, $493.36,497,501.002,507.002,513.004,517.005,526.005, $521.005,526.006,533.006,538.009,547.009,554.009,565.008, $570.009,574.01,571.011,563.011,558.01,554.011,550.011, $557.011,549.012,545.012,551.01,544.01,539.011,546.017, $542.019,549.02,553.022,560.021,567.021,575.022,568.024, $572.025,582.024,586.023,591.024,594.025,599.029,604.029, $600.024,606.022,617.022,627.021,633.023,639.024,645.025, $651.025,656.025,651.022,646.021,637.021,629.019,623.018, $617.018,610.017,605.018,599.019,593.019,588.018,575.017, $990.37,579.019,573.019,569.019,576.019,579.019,990.37, $575.017,565.016,560.016,553.014,562.013,570.012,580.012, $591.011,595.011,589.01,583.009,579.008,587.006,596.006, $606.006,615.005,621.006,625.007,629.008,636.01/ data (c(j),j=2991,3120)/ 644.011,649.012,658.013,666.014, $671.015,678.016,683.016,689.016,686.017,689.018,694.018, $698.02,702.023,706.024,709.027,706.03,990.37,765.017, $770.015,774.014,780.014,783.016,787.012,791.011,796.011, $792.016,799.016,794.019,790.02,786.022,782.022,778.024, $774.024,778.022,784.02,780.019,774.018,765.017,990.37, $794.021,798.022,802.018,798.027,794.026,791.024,794.023, $794.021,990.37,706.03,700.03,697.031,693.033,689.036, $697.037,692.039,680.04,676.041,670.041,664.041,661.04, $664.035,657.035,648.035,644.035,639.037,644.037,648.037, $645.04,654.04,658.041,665.042,661.044,669.045,680.044, $686.043,678.047,674.045,670.046,667.046,676.048,680.05, $990.37,599.316,606.314,610.312,615.311,622.31,629.309, $636.309,641.308,651.308,656.307,659.306,664.307,670.306, $675.306,681.307,684.308,689.309,693.309,697.309,700.308, $697.308,693.307,701.305,707.305,710.308,713.308,717.307, $720.305,723.305,727.305,730.305,736.304,741.304,744.303, $748.303,754.301,758.301,762.298,766.291,769.289,773.291, $777.289,780.287,785.287,788.29,793.294,797.295,800.295/ data (c(j),j=3121,3217)/ 803.292,807.294,811.296,814.298, $819.299,816.301,819.301,823.305,820.308,825.308,822.312, $818.314,823.313,827.312,830.312,833.32,837.331,833.333, $830.336,827.339,823.338,820.328,817.326,820.329,817.333, $820.336,816.338,812.336,816.34,811.347,807.346,801.343, $796.342,791.341,781.341,775.341,772.342,768.342,762.338, $757.341,753.342,750.343,745.341,742.341,738.34,734.34, $729.338,724.338,721.338,724.336,720.337,716.338,704.338, $712.336,705.335,702.334,697.337,690.335,687.334,682.331, $672.327,667.326,662.325,655.323,649.32,642.32,636.32, $631.319,621.318,614.318,608.317,601.317,599.316,990.37, $633.341,637.339,646.338,650.338,653.338,658.336,663.337, $658.339,655.339,652.339,656.339,660.34,657.341,660.341, $665.344,662.345,657.346,651.346,646.346,641.345,637.343, $633.341,990.37/ end c_____________________________________________________________________ c======================================================================= c UNIT M: MAPPING AND PLOTTING c======================================================================= subroutine atlas(kind) c$$$$ calls getchr getarr getone rotate mult c Determines map projection (kind) and prepares various constants c for it in common /orient/. c Finds the array size and proportions and the appropriate step c size and limits for the projection; returned to caller in /scan/. c kind=1: Radian projection c kind=2: Lambert equal-area c kind=3: Hammer-Aitoff equal area (A) c kind=4: Orthographic projection c kind=5: Mercator projection c character*1 type dimension unit(3),axel(3),O1(3,3),O2(3,3),ulim(5) common /scan/ m,n,umin,vmin,du,factor,hash common /orient/ center(2),O(3,3) common /prolix/ iprint save type, em,iadd c data type /'R'/, em/40/, deg/0.01745329252/, axel(3)/0/ data unit /0,0,1/, enlarg/1/, pi/3.14159265/, iadd/0/ ce call getchr('projection', type, nf) kind=index('rlaomRLAOM', type) if (kind .eq. 0) then write(*,'(3a)')'>>> Unknown map projection: ',type, $ '. Magmap substitutes proj=r' type='r' kind=1 endif if (kind .gt. 5) kind=kind-5 c call getarr('center', center, 2, ncen) if (ncen .eq. 1) center(2)=center(1) c call getone('dimensions', em, ignor) c call getone('magnify', enlarg, ignor) factor=1.0/max(1.0, enlarg) if (enlarg .gt. 1.0 .and.iprint.ge.0) $write(*,'(/a,f8.2)') $' Map enlarged by the factor ',enlarg hash=kind+1/(factor+100/(center(1)+100/(center(2)+.123))) c c Alternative specification: directly by a window call getarr('window', ulim, 5, jwin) if (jwin .eq. 5 .and. ulim(1).lt.ulim(2) .and. $ ulim(3).lt.ulim(4) .and. ulim(5).gt.0.0) then umin=ulim(1) vmin=ulim(3) du=ulim(5) m=1 + nint((ulim(4)-vmin)/du) n=1 + nint((ulim(2)-umin)/du) iadd=5 hash=-kind+1/(umin+1/(vmin+1/(du+.567))) elseif (jwin .eq. 0) then iadd=0 endif igo=kind + iadd c c Prepare the projections goto (1000, 2000, 3000, 4000, 5000, $ 1010, 2010, 3010, 4010, 5010), igo c c kind=1: radian projection 1000 m=nint(em) n=2*m - 1 umin=-180.0*factor vmin=-90.0*factor du=-umin/(m-1) c 1010 if (ncen .eq. 1) center(2)=center(1) if (iprint.ge.0) $write(*,'(/a/a,f8.1)')' Rectilinear projection', $' Central meridian ',center(2) return c c kind=2: Lambert equal-area 2000 m=nint(em) n=m umin=-1.4142*factor vmin=umin du=-2.0*umin/(m-1) c 2010 if (center(1) .ge. 90.0) then O(3,3)=2.0 else c Build the rotation matrix used in projx; arrange meridian through c the map center to be vertical on the page; if center is one of the c poles, the Greenwich meridian is vertical axel(1)= sin(deg*center(2)) axel(2)=-cos(deg*center(2)) call rotate(center(2), unit, O1) call rotate(center(1)-90.0, axel, O2) call mult(3,3,3, O2,O1,O) endif c if (iprint.ge.0) $write(*,'(/a/a,2f8.1)')' Lambert equal-area projection', $' Geographic coordinates of map center ',center c return c c kind=3: Hammer-Aitoff equal area 3000 m=nint(em) n=2*m - 1 umin=-2.0*factor vmin=-factor du=-umin/(m-1) O(1,1)= cos(deg*center(2)) O(2,2)=-sin(deg*center(2)) 3010 if (iprint.ge.0) $write(*,'(/a/a,f8.1)')' Hammer-Aitoff equal-area projection', $' Central meridian ',center(2) return c c kind=4: Orthographic projection 4000 m=nint(em) n=m umin=-factor vmin=umin du=-2.0*umin/(m-1) c 4010 if (center(1) .ge. 90.0) then O(3,3)=2.0 else c Build the rotation matrix used in projx; arrange meridian through c the map center to be vertical on the page; if center is one of the c poles, the Greenwich meridian is vertical axel(1)= sin(deg*center(2)) axel(2)=-cos(deg*center(2)) call rotate(center(2), unit, O1) call rotate(center(1)-90.0, axel, O2) call mult(3,3,3, O2,O1,O) endif c if (iprint.ge.0) $write(*,'(/a/a,2f8.1)')' Orthographic projection', $' Geographic coordinates of map center ',center return c c kind=5: Mercator projection 5000 m=nint(em) n=1.5*m - 1 du=2.0*pi*factor/(n-1) umin=-pi*factor vmin=-0.5*du*(m-1) c 5010 if (iprint.ge.0) $write(*,'(/a/a,f8.1)')' Mercator projection', $' Central meridian ',center(2) return c end c_______________________________________________________________________ blockdata middle common /orient/ center(2),O(3,3) data center/90,0/, O(3,3)/2/ end c_______________________________________________________________________ subroutine projx(iproj, u, x) c$$$$ calls nothing c Applies the inverse map projection from a point u in R**2 on the map c plane to the unit vector in R**3 in the unit sphere dimension u(2), x(3) common /orient/ center(2),O(3,3) data deg/0.01745329/ c goto (1000, 2000, 3000, 4000, 5000), iproj c c Rectilinear projection: u(1)=longitude, u(2)=latitude c 0 <= u(1) <= 360; -90 <= u(2) <= +90 c 1000 cosla=cos(deg*u(2)) x(1)=cosla*cos(deg*(u(1)+center(2))) x(2)=cosla*sin(deg*(u(1)+center(2))) x(3)=sin(deg*u(2)) return c c Lambert equal area: before rotations (0,0) in u plane is (0,0,1) in c R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 4 is mapped into the (punctured) sphere 2000 sq=u(1)**2 + u(2)**2 R=sqrt(1.0 - 0.25*sq) x(1)=-R*u(2) x(2)= R*u(1) x(3)=1.0 - 0.5*sq if (O(3,3) .ge. 2.0) return c Apply coordinate rotation x1 =O(1,1)*x(1) + O(1,2)*x(2) + O(1,3)*x(3) x2 =O(2,1)*x(1) + O(2,2)*x(2) + O(2,3)*x(3) x(3)=O(3,1)*x(1) + O(3,2)*x(2) + O(3,3)*x(3) x(2)=x2 x(1)=x1 return c c Hammer-Aitoff equal-area projection: u(1)**2/4 + u(2)**2 < 1 is c mapped onto the unit sphere 3000 q=sqrt(2.0 - u(2)**2 - 0.25*u(1)**2) x(3)=q*u(2) cosla=sqrt(1.0 - x(3)**2) sinh=0.5*u(1)*q/(cosla + 1.0e-10) x2=2.0*sinh*sqrt(1.0 - min(1.0, sinh**2))*cosla x1=(1.0 - 2.0*sinh**2)*cosla c Apply coordinate rotation - only about z axis x(1)= x1*O(1,1) + x2*O(2,2) x(2)=-x1*O(2,2) + x2*O(1,1) return c c Orthogonal projection: before rotations (0,0) in u plane is (0,0,1) c in R**3, and negative u(2) is the Greenwich meridian. The disk c u(1)**2 + u(2)**2 < 1 is mapped into the (punctured) sphere 4000 R=sqrt(u(1)**2 + u(2)**2) if (R .gt. 1.0) then c (u, v) lies off the sphere; place image on limb x(1)=-u(2)/R x(2)= u(1)/R x(3)= 0.0 R=1.0 else x(1)=-u(2) x(2)= u(1) x(3)= sqrt(1.0 - R**2) endif if (O(3,3) .ge. 2.0) return c Apply coordinate rotation x1 =O(1,1)*x(1) + O(1,2)*x(2) + O(1,3)*x(3) x2 =O(2,1)*x(1) + O(2,2)*x(2) + O(2,3)*x(3) x(3)=O(3,1)*x(1) + O(3,2)*x(2) + O(3,3)*x(3) x(2)=x2 x(1)=x1 return c c Mercator projection: u(1)=longitude (radians), u(2)=M c -pi <= u(1) <= pi; -pi <= u(2) <= +pi c 5000 cosla=1.0/cosh(u(2)) x(1)=cosla*cos(u(1)+deg*center(2)) x(2)=cosla*sin(u(1)+deg*center(2)) x(3)=tanh(u(2)) return c end c_______________________________________________________________________ subroutine rotate(theta, en, a) c$$$$ calls nothing c Finds the 3 x 3 rotation matrix a that transforms vectors rotating c them about the axis en by an angle theta (degrees) anticlockwise c looking from origin along the axis. c c Rotation about the axis n, where n is unit by the angle A. c written intelligibly the matrix is really c c R = I cos A + n n (1-cos A) + sin A n x c c where n x means "n cross". dimension en(3), a(3,3) data rad/0.01745329252/ c d=sqrt(en(1)**2 + en(2)**2 + en(3)**2) if (d .eq. 0) write(*, *)'>>> Zero length rotation axis' c cost=cos (rad*theta) sint=sin (rad*theta)/d sinq=2.0*(sin(rad*theta/2.0)/d)**2 c a(1,1)=cost + sinq*en(1)*en(1) a(1,2)= sinq*en(1)*en(2) - sint*en(3) a(1,3)= sinq*en(1)*en(3) + sint*en(2) a(2,1)= sinq*en(2)*en(1) + sint*en(3) a(2,2)=cost + sinq*en(2)*en(2) a(2,3)= sinq*en(2)*en(3) - sint*en(1) a(3,1)= sinq*en(3)*en(1) - sint*en(2) a(3,2)= sinq*en(3)*en(2) + sint*en(1) a(3,3)=cost + sinq*en(3)*en(3) c return end c______________________________________________________________________ subroutine mult(m,n,k, a,b,c) c$$$$ calls nothing c Elementary matrix multiplication: c = a b where array sizes are c those appearing in the dimension statement below dimension a(m,n), b(n,k), c(m,k) c do 1300 i=1, m do 1200 j=1, k c(i,j)=0.0 do 1100 l=1, n c(i,j)=c(i,j) + a(i,l)*b(l,j) 1100 continue 1200 continue 1300 continue c return end c______________________________________________________________________ subroutine mapper c$$$$ calls atlas getchr getarr getone shsum projx plot c Calculates a field model on a regular array, which is in a specified c map plane. writes the required field element to a diskfile parameter (pi=3.14159265358979, deg=180.0/pi) parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) character*80 name, oldnam, code*1, type*1, nrm*1 complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /scan/ m,n,umin,vmin,du,factor,hash common /orient/ center(2),O(3,3) common /prolix/ iprint dimension u(2),r(3),fld(500) save code, shrnk, oldnam, kount, radius c data code/'R'/, bdum/1.0e8/, shrnk/1.0/, kount/0/ data type/'N'/, name,oldnam/'fort.9',' '/, radius/1.0/ c c Check to see if a map is needed ce call getchr('projection', type, nf) if (type .eq. 'N' .or. type.eq.'n') then write(*,*)' No map required in this run' return endif kount=kount + 1 c call getchr('normalize', nrm, ignor) call getone('igrf', yr, jgrf) if (nrm.eq.'S' .or. jgrf.gt. 0) then shrnk=1.0e-3 write(*,'(/1x,a)') $ 'Gauss coeffs assmed to be in nT; then mapped fields are in muT' endif c c Get the file name for array call getchr('output', name, nnam) if (name .eq. oldnam) then write(*,'(2a/a)') $'>>> Warning - Previous array of field values will be overwritten' $ ,'>>> A new filename should be set with an output command' elseif (nnam .le. 0) then if (kount .gt. 1) write(name,'(a5,i2)') 'fort.',88+kount write(*,*)'>>> No output file assigned; default to: ',name endif close (unit=9) open(unit=9, file=name) c c Get output array size and projection information; found in /scan/ call atlas(iproj) if (iprint.ge.0) $write(*,'(/a,2i4)')' Dimension of output array',m,n c c Get the field element desired; map onto an index call getchr('field', code, nf) le=index('RXYZFHVID',code) if (le .eq. 0) then write(*,'(3a)')'>>> Unknown field element: ',code, $ '. Magmap assumes field=R' code='R' le=1 endif if (iprint.ge.0) $write(*,'(/(a,a))')' Map of field element: ',code, $ ' written to output file: ',name(1:35) if (iprint.ge.0) $write(*,'(a,2f8.3,a,2f8.3,a,f8.3)')' Map plane = (', $umin,umin +(n-1)*du,')x(',vmin,vmin +(m-1)*du,'); dx=dy=',du c bmin=1e8 bmax=-bmin call getone('radius', radius, ignor) c c Scan through u plane: vmin,umin, etc come through /scan/ from atlas do 1500 i=1, m u(2)=vmin + (m-i)*du do 1400 j=1, n u(1)=umin + (j-1)*du call projx(iproj, u, r) c Insert dummy value for illegal coordinates if (r(1)**2+r(2)**2+r(3)**2 .eq. 0.0) then fld(j)=bdum goto 1400 endif c Sum SH series r(1)=radius*r(1) r(2)=radius*r(2) r(3)=radius*r(3) c call shsum(r, lmax, blm, b1,b2,b3,V) if (lext.gt. 0) then call shsum(r, -lext, alm, exb1,exb2,exb3,exV) b1=b1 + exb1 b2=b2 + exb2 b3=b3 + exb3 V =V + exV endif c c Select the kind of field to be mapped goto (1010,1020,1030,1040,1050,1060,1070,1080,1090),le c Radial field desired: 'R' 1010 fld(j)=b1 goto 1400 c North field desired: 'X' 1020 fld(j)=-b2 goto 1400 c East field desired: 'Y' 1030 fld(j)=b3 goto 1400 c Downward field desired: 'Z' 1040 fld(j)=-b1 goto 1400 c Total field desired: 'T' 1050 fld(j)=sqrt(b1**2 + b2**2 + b3**2) goto 1400 c Horizontal field desired: 'H' 1060 fld(j)=sqrt(b2**2 + b3**2) goto 1400 c Scalar potential V/a: 'V' 1070 fld(j)=V goto 1400 c Inclination desired: 'I' 1080 fld(j)=atan(-b1/sqrt(b2**2+b3**2))*deg/shrnk goto 1400 c Declination desired: 'D' 1090 fld(j)=atan2(b3, -b2)*deg/shrnk c 1400 continue write(9,'(1p,g12.4)') (shrnk *fld(j),j=1,n) do 1450 j=1, n if (fld(j) .ne. bdum) then bmin=min(fld(j), bmin) bmax=max(fld(j), bmax) endif 1450 continue 1500 continue close(unit=9) c c Inform user about sizes if (iprint.ge.0) $write(*,'(/3(a,1p,g12.3))') $' Field values between',bmin*shrnk,' and',bmax*shrnk,' muT' c c Display results if requested call plot(m,n, iproj, name) c return end c_______________________________________________________________________ subroutine launch(lk, script) c$$$$ calls system c CAUTION: the call to system is not a part of the ANSI standard and c may have to be deleted if taken out of the local environment. c The calls to launch can be safely deleted if necessary. c Runs the plotfile in the file script written by magmap. character*64 script c if (script .eq. ' ') return close(unit=10) c Contour and color are plot utilities written by R.L.Parker c available at http://igppweb.ucsd.edu/~parker if (lk.eq.0) call system('sleep 2; contour <'//script) if (lk.eq.1) call system('sleep 2; color <'//script) c Any PostScript screen display can be substituted for gv call system('gv mypost &') return end c_______________________________________________________________________ subroutine getdat c$$$$ calls getchr getone getarr serial bang c Input routine. Accepts SH coeffs from a diskfile of the scalar c magnetic potential V. The default is fully normalized complex c values, integrating to unity over the sphere. The coeffs are stored c this way. These are the coeffs loaded into common /coeff/. c c Each input line should contain l, m, real, imag. c Or equivalently l, m, cosine coeff, sine coeff. c Negative m coeffs must not be provided: they are assumed to be c those appropriate to real omega. c The coeffs may be Schmidt normalized referred to the surface, c or one of 3 fully normalized conventions c c c parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) character*80 name, style,line, norm*2 dimension sum(4),w(4) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /prolix/ iprint common /earth/ rad(4) save name, norm, radius, year c data name/' '/, norm/' '/, radius/1/, pi/3.14159265/ data factor/1/, year/0/, minyr/1899/ c c c Check first if internally stored IGRF model is requested call getone('igrf', year, inunit) if (inunit.ge.1) then if (1900.0 .le.year .and. year.le. 2010.0) then inunit=nint(year) norm='S' if (iprint.ge.0) $ write(*,'(a,i5/2a)')' Field evaluated for IGRF ',inunit, $ ' Using IAGA IGRF-11: ', $ '' goto 1010 elseif (year .ne. 0.0) then write(*,'(a,f7.1)')'>>> Inappropriate year for IGRF model:', $ year stop endif endif c c Ascertain the name of the input file call getchr('file', name, nfile) if (nfile .le. 0 .and. name(1:1) .eq. ' ') then write(*,*) $ '>>> Program cannot continue without filename for SH coeffs' stop elseif (name(1:1) .eq. '*') then if (iprint.ge.0) write(*,'(/a)') $ ' SH coefficients input from the terminal' inunit=5 else if (iprint.ge.0) write(*,'(/a,a)') $ ' SH coeffs from file: ',name(1:51) inunit=7 open (unit=7, status='OLD', err=2100, file=name) endif if (nfile .le. 0 .and. name(1:1) .ne. ' ') then if (iprint.ge.0) write(*,'(a)') $ ' Program continues with SH coefficients previously input' return endif c call getchr('normalize', norm, ignor) if (norm .eq. ' ') then write(*,*)'>>> Program cannot continue without SH normalization' stop endif if (index('S F1F2G1U ',norm) .eq. 0) then write(*,*)'>>> Normalization must be one of: S F1 F2 G1 U' stop endif if (iprint.ge.0) $write(*,'(a,a)')' Normalization of coeffs: ',norm c 1010 call getone('radius', radius, nrad) if (radius .eq. 0) radius=0.54701 rad(4)=radius c c Set the earth dimensions: call getarr('earth', rad, 3, nrad) if (iprint.ge.0) write(*,'(a,f9.3)') $' Equatorial radius (km) ',rad(1), $' Polar radius (km) ',rad(2), $' Ref radius aref (km) ',rad(3), $' Evaluation r/aref ',rad(4) c c Zero the arrays do 1100 n=1, mxb blm(n) = 0.0 if (n.le.mxa) alm(n)=0.0 1100 continue c c In contrast to other parameters, lmax is reset to maximum value c each time through the SH input process elmax=sqrt(2*mxb + 0.25) - 0.5 nexmx=sqrt(2*mxa + 0.25) - 0.5 call getone('lmax', elmax, nfound) lmax=elmax call getone('scale', factor, ignor) if (factor .ne. 1.0 .and.iprint.ge.0) $write(*,'(a,1p,g11.4)')' Coeffs scaled by',factor c c Check normalization f=sqrt(2.0*3.1415926535) rt2=sqrt(2.0) ff=factor c Norm is F2 if (norm .eq. 'F2') ff=factor*f*rt2 c c Check for serial data; read in file here call getchr('serial', style, nser) if (nser .ge. 0 .or. inunit .ge. minyr) then call serial(inunit, 0, ll, mm, glm, hlm) if (inunit .ge. minyr) style=' ' inunit=0 endif c call getone('list', dum, nlist) c lmx=0 kmx=0 bmx=0 if (nlist.ge.0 .and. inunit.ne.5 .and.iprint.ge.0) $write(*,'(a)') ' ', ' Table of coefficients as read' if (inunit.eq.5 .and.iprint.ge.0) $write(*,'(/(a))') ' Enter SH Coefficients NOW:', $' Enter l, m, glm, hlm each on a new line; l < -499 terminates' c c c Input the SH file one line at a time. do 1200 k=1, mxb c c Read the line of values keying on inuit for source if (inunit.eq. 7) then 1175 read (7,'(a80)', end=1250) line if (line(1:1) .eq. '%') then if (nlist.ge. 0 .and. iprint.ge. 0) write(*,'(a)') line goto 1175 endif read (line, *, err=2000) ldeg,mm, glm, hlm elseif (inunit.eq. 5) then read (*,*, err=2000) ldeg,mm, glm, hlm elseif (inunit.eq. 0) then call serial(0, k, ll, mm, glm, hlm) ldeg=ll endif ll=abs(ldeg) c if (nlist.ge. 0 .and. iprint.ge. 0) $ write(*,'(2i5,2f12.2)') ldeg,mm,glm,hlm if (ldeg.le. -500) goto 1250 if (ll.gt. lmax .or. -ldeg.gt.nexmx) goto 1200 lmx=max(lmx, ldeg) kmx=max(kmx, -ldeg) c Convert Schmidt semi-normalized to fully normalized coefficients if (norm .eq. 'S') then vl=f*(-1.0)**mm /sqrt(2.0*ll+1.0) glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =rt2*glm c Convert unnormalized to fully normalized coefficients elseif (norm .eq. 'U') then vl=f*(-1.0)**mm*sqrt(bang(ll+mm)/(bang(ll-mm)*(4.0*ll+2.0))) glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =2.0*glm c Norm is G1 elseif (norm .eq. 'G1') then vl=f*(-1.0)**mm glm = vl*glm hlm =-vl*hlm if (mm .eq. 0) glm =rt2*glm endif c Scale and store in temp array j=1 + mm + (ll*(ll+1))/2 if (ldeg.ge. 0) blm(j)=ff*cmplx(glm, hlm) if (ldeg.lt. 0) alm(j)=ff*cmplx(glm, hlm) bmx=max(abs(blm(j)), bmx) 1200 continue c c 1250 if (inunit.eq.7) close (unit=7) if (bmx .eq. 0) bmx=1 lmax=lmx lext=kmx if (iprint.ge.0) write(*,'(a,i5)') $' Maximum degree of internal coefficients:',lmax, $' Max degree of external coeffs:',lext if (iprint.ge.0) $write(*,'(a,f9.5)')' Normalized radius r/aref = ',radius c Check if dipole term(s) are to be dropped for this read only trm=3.0 call getone('nondipole', trm, nondip) if (nondip .ge. 0) then if (iprint.ge.0) then if (trm.eq. 3)write(*,*)'Dipole terms deleted if present' if (trm.ne. 3)write(*,*)'Axial dipole term deleted' endif blm(2)=0.0 if (trm.eq. 3.0) blm(3)=0.0 endif c c List the fully normalized coefficients if requested call getone('list', dum, nlist) if (nlist .ge. 0) then if (iprint.ge.0) write(*,'(/a,f9.5,a,f9.3,a)') $ ' Fully normalized SH coefficients referred to radius',radius, $ ' = ',radius*rad(3),' km' bscale=10.0**(-3*int((log10(bmx)-4.0)/3.0)) if (bscale .ne. 1.0 .and.iprint.ge.0) $ write(*,'(a, 1p,g10.2,a)') $ '(Values in the table must be multiplied by ',1.0/bscale,')' do 1350 l=0, lmax j=1 + (l*(l+1))/2 if (iprint.ge.0)write(*,110) l,(blm(k)*bscale,k=j, j+l) 110 format(i4,6f12.2:/(f16.2,5f12.2:)) 1350 continue if (lext.gt. 0) then write(*,'(a)')' External field coefficients' do 1360 l=1, lext j=1 + (l*(l+1))/2 if (iprint.ge.0) write(*,110) l,(alm(k)*bscale,k=j, j+l) 1360 continue endif endif c c Print weighted sums call getone('sums', dum, nesum) if (nesum .lt. 0) return j=0 do 1600 ns=1, 4 sum(ns)=0 1600 continue do 1800 l=1, lmax el=l w(1)=el+1 w(2)=(el+1)**2 w(3)=(el+1)*(2*el+1)**2*(2*el+3)/el w(4)=el*(el+1)**3 do 1700 m=0, l j=j + 1 bmagsq=2.0*abs(blm(j))**2 if (m .eq. 0) bmagsq=0.5*bmagsq do 1650 ns=1, 4 sum(ns)=sum(ns) + w(ns)*bmagsq 1650 continue 1700 continue 1800 continue c if (iprint.ge.0) write(*,'(/a, 1p, g12.4, a)') $' Field energy:',sum(1), '*a**3/(2*mu0)', $' RMS radial field:',sqrt(sum(2)/(4*pi)), ' ', $' Lower bound on ohmic heating:',sum(3), '*a/(sigma*mu0**2)', $' RMS grad Br:', sqrt(sum(4)/(4*pi)),'/a' return c c Error return 2000 write(*,*)'>>> Unable to read input line from disk' write(*,*)'>>> Last successful l, m:',ll,mm stop 2100 write(*,*)'>>> Cannot open the file ',name stop end c_______________________________________________________________________ function bang(k) c$$$$ calls nothing c Factorial of integer; no checks for negative argument or overflow bang=1.0 do 1100 j=1, k bang=bang * j 1100 continue return end c_______________________________________________________________________ blockdata tellus c Set default values in km for the equatorial, polar and reference c radii of the earth common /earth/ rad(4) data rad/ 6378.17, 6356.91, 6371.20, 1.0 / c end c_______________________________________________________________________ subroutine serial(inunit,k, l,m, br,bi) c$$$$ calls getchr getone c Reads serial data or returns values upon request. c If inunit = 0 returns kth pair (br, bi) in SH coefficient list and c If 0 < inunit < minyr reads from disk file to eof c If minyr < inunit uses inunit as year for IGRF field model c then corresponding degree and order are l, m. c c If nondipole field is specified, routine assumes series starts at l=2 c parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) parameter (lmmx=2*mxb) character*20 style, line*80 common /igrf/ ibeg(24),ghlm(20*130+3*208) dimension b1(lmmx) save nondip, nozer, kmx, b1 data minyr/1899/ c c c Either fill b1 array or return the l,m pair if (inunit .gt. 0) then c c nondip=2 if dipoles are dropped, 0 otherwise dum=3.0 call getone('nondipole', dum, nondip) if (dum .ne. 3.0) write(*,*) $ '>>> All l=1 terms deleted in serial input mode' nondip=2 + max(-2, nondip) c c Either read data from a file or from igrf common if (inunit .lt. minyr) then c Reads coefficients serially, absent l, m labels, saved in b1 c until needed i1=1 do 1500 l=1, lmmx 1100 read (inunit,'(a)',end=1600) line if (line(1:1) .eq. '%') goto 1100 read (line, *, end=1200, err=5000) (b1(i),i=i1, lmmx) 1200 i1=i 1500 continue 1600 kmx=i1-1 c c Get igrf coefficients from common statement else k1=1 + (inunit - 1900)/5 ibgk=ibeg(k1) + 2*nondip kend= ibeg(k1+1)-1 do 2100 ik=ibgk, kend b1(ik - ibgk +1)=ghlm(ik) 2100 continue kmx=ibeg(k1+1) - ibgk endif c c Coefficients entered into b1 array if (inunit.eq.7) close (unit=inunit) call getchr('serial', style, ignor) nozer=min(1, index(style, 'out')) else c c c Retrieve SH pairs from b1. First find l, m; then work out index in b1 kk=k + nondip l=sqrt(0.251+2*kk) - 0.5 m=kk - (l*(l+1))/2 n=2*(kk - nondip) - 1 - nozer*(l - nondip/2) if (m .eq. 0) n=n+nozer if (n .gt. kmx) l=-999 br=b1(n) bi=b1(n+1)*min(1,m) endif return c c Error return 5000 write(*,*)'>>> Error in serial data on line ', l write(*,*) line stop end c_______________________________________________________________________ subroutine inpsh c$$$$ calls clear icheck minlen common /prolix/ iprint c Input routine for commands to read SH coeffs and plot maps c c Reads from the standard input until eof or code "continue" or "quit". c Saves lines in the input store /store/ for later retrieval by c getarr, getone or getchr. c c Prints a glossary upon request c parameter (inmx=200) character *80 input(inmx),line,cmd*4 common /store/ input common /ndict/ iecho,nin,memory,istate(inmx) c write(*,'(a)') ' ', $' =================',' ', $'Enter commands for spherical harmonic coefficients (? for help)' c do 1500 l=nin+1, inmx read(*,'(80a)', end=3000) line if (line(1:5) .eq. 'conti' .or. line(1:4) .eq. 'exec') goto 2000 c List a glossary of codes if (line(1:1) .eq. '?') write(*, '(a/a/(2x,a))') $'Enter commands from the following list:', $' (M) means mandatory information',' ', $'?: Remind me of the command list again', $'execute: Quit reading commands and begin calculations', $'quit: End calculations', $'file filename: Input the filename of coeffs', $'serial [without zeros] SH coeffs in natural order without l,m', $'normalization X: X=S Schmidt, X=F1 fully normalized, integral=1', $' X=F2 integral=4*pi, X=G1 Gravity normal, U=unnormalized (M)', $'igrf year: Use internal IGRF coeffs in place of file input', $'lmax: Accept coefficients up to this degree', $'nondipole: Delete dipole terms from all subsequent input', $'list: List SH file as read and in F1 normalization', $'radius r: Calculate field & plot map on sphere radius = r*aref', $'earth a,b,aref: Equatorial, polar, ref radius (km)', $' ---------', $'output filename: Diskfile to receive map array', $'silent: Suppress terminal output henceforth (except error mssg)', $'dimensions m, n: Array dimensions' if (line(1:1) .eq. '?') write(*, '(2x,a)') $'center lat, long: Geographic coordinates of map center', $'projection p: p=l Lambert equal area, p=r radian, p=a Aitoff, ', $' p=o orthographic, m=Mercator ', $'noplot: Save map array but do not plot it', $'magnification f: Draw map with scale f times bigger than normal', $'coast [file]: Lat-long pairs for coast; blank =internal dataset', $'symbols file type ht [color=X]: File of lat-longs for symbols', $'field E: E=one of (X, Y, Z, H, F, R, V/a, I, D)', $'scale s: Multiply SH coeffs by s before plotting', $'site lat, long, alt: Evaluate field elements at site with given', $' geographic coords and altitude in meters', $'sums: Print certain weighted sums of coeffs', $'contour [file]: Prepare a contour script in file and run it', $'color [file]: Prepare a color contour script in file and run it', $'clear cmd: Remove all instances of cmd from stack', $' ' if (line(1:1) .ne. ' ') then if (line(1:4) .eq. 'quit') goto 3000 if (line(1:4) .eq. 'echo') iecho=-iecho if (line(1:4) .eq. 'sile') iprint=-1 c Alternative file designators suppress each other if (line(1:4) .eq. 'file') call clear('igrf') if (line(1:4) .eq. 'igrf') call clear('file') c Alternative plotting formats suppress each other if (line(1:4) .eq. 'cont') call clear('colo') if (line(1:4) .eq. 'colo') call clear('cont') c nin=nin + 1 call icheck(line) input(nin)=line istate(nin)=0 endif 1500 continue write(*,'(a,i5,a/a)')'>>> Command stack full: only',inmx, $' allowed','>>> Recompile MAGMAP with a larger inmx' stop c 2000 if (iprint.ge.0) then write(*, '(5x,a)') ' ','===================' do 2200 i=1, nin k=minlen(input(i)) write (*,'(5x,a)') input(i)(1:k) 2200 continue write(*,'(5x,a)') ' ','===================',' ' else write(*,'(/a/a)')' Magmap screen output suppressed', $ '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~',' ' endif call getchr('clear', cmd, klr) c Remove cleared commands if (klr .gt. 0) call clear(cmd) if (klr .gt. 0) call clear('clear') c return c 3000 write(*,*) ' ',' Magmap run complete' call plot(0,0,0,' ') stop end c______________________________________________________________________ function minlen(str) c$$$$ calls nothing c Returns length of the string str with trailing blanks removed. c If str is blank returns 1 character*(*) str c minlen=1 n=len(str) do 1100 k=n, 1, -1 minlen=k if (str(k:k) .ne. ' ' .or. k.eq. 1) return 1100 continue end c_________________________________________________________________ subroutine icheck(line) c$$$$ calls nothing c Checks the command string line (1:4) against a list; appends a c warning if command is not present in the list character*80 line, com*4 ce com=line(1:4) if (0 .eq. $+index('cent clea coas colo cont dime eart fiel file igrf', com) $+index('list lmax magn nond norm outp plot proj radi scal', com) $+index('seri site sums symb wind nopl', com) $) line=line(1:20)//' <<<<<<< Unrecognized command' return end c______________________________________________________________________ blockdata iounit parameter (inmx=200) common /prolix/ iprint common /ndict/ iecho,nin,memory,istate(inmx) data iprint /1/, iecho/-1/, nin/0/ end c______________________________________________________________________ blockdata igref c IAGA IGRF-11: 1900-2010. All degree 10, until 2000, then deg 13 c See common /igrf/ ibeg(24), $y1900(130), y1905(130), y1910(130), y1915(130), y1920(130), $y1925(130), y1930(130), y1935(130), y1940(130), y1945(130), $y1950(130), y1955(130), y1960(130), y1965(130), y1970(130), $y1975(130), y1980(130), y1985(130), y1990(130), y1995(130), $y2000(208), y2005(208), y2010(208) data ibeg/1,131,261,391,521,651,781,911,1041,1171,1301, $1431,1561,1691,1821,1951,2081,2211,2341, $2471,2601,2809,3017, 3225/ c 1995 2000 2005 2010 c data y1900/ $-31543,0,-2298,5922,-677,0,2905,-1061,924,1121, $1022,0,-1469,-330,1256,3,572,523,876,0,628,195, $660,-69,-361,-210,134,-75,-184,0,328,-210,264,53, $5,-33,-86,-124,-16,3,63,0,61,-9,-11,83,-217,2, $-58,-35,59,36,-90,-69,70,0,-55,-45,0,-13,34, $-10,-41,-1,-21,28,18,-12,6,-22,11,0,8,8,-4, $-14,-9,7,1,-13,2,5,-9,16,5,-5,8,-18,8,0,10, $-20,1,14,-11,5,12,-3,1,-2,-2,8,2,10,-1,-2,-1, $2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2, $4,2,0,0,-6/ data y1905/ $-31464,0,-2298,5909,-728,0,2928,-1086,1041,1065, $1037,0,-1494,-357,1239,34,635,480,880,0,643,203, $653,-77,-380,-201,146,-65,-192,0,328,-193,259,56, $-1,-32,-93,-125,-26,11,62,0,60,-7,-11,86,-221, $4,-57,-32,57,32,-92,-67,70,0,-54,-46,0,-14,33, $-11,-41,0,-20,28,18,-12,6,-22,11,0,8,8,-4,-15, $-9,7,1,-13,2,5,-8,16,5,-5,8,-18,8,0,10,-20, $1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2,-1,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2,4, $2,0,0,-6/ data y1910/ $-31354,0,-2297,5898,-769,0,2948,-1128,1176,1000, $1058,0,-1524,-389,1223,62,705,425,884,0,660,211, $644,-90,-400,-189,160,-55,-201,0,327,-172,253,57, $-9,-33,-102,-126,-38,21,62,0,58,-5,-11,89,-224, $5,-54,-29,54,28,-95,-65,71,0,-54,-47,1,-14,32, $-12,-40,1,-19,28,18,-13,6,-22,11,0,8,8,-4,-15, $-9,6,1,-13,2,5,-8,16,5,-5,8,-18,8,0,10,-20, $1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2,-1,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,2,4, $2,0,0,-6/ data y1915/ $-31212,0,-2306,5875,-802,0,2956,-1191,1309,917, $1084,0,-1559,-421,1212,84,778,360,887,0,678,218, $631,-109,-416,-173,178,-51,-211,0,327,-148,245, $58,-16,-34,-111,-126,-51,32,61,0,57,-2,-10,93, $-228,8,-51,-26,49,23,-98,-62,72,0,-54,-48,2, $-14,31,-12,-38,2,-18,28,19,-15,6,-22,11,0,8,8, $-4,-15,-9,6,2,-13,3,5,-8,16,6,-5,8,-18,8,0, $10,-20,1,14,-11,5,12,-3,1,-2,-2,8,2,10,0,-2, $-1,2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2, $1,4,2,0,0,-6/ data y1920/ $-31060,0,-2317,5845,-839,0,2959,-1259,1407,823, $1111,0,-1600,-445,1205,103,839,293,889,0,695,220, $616,-134,-424,-153,199,-57,-221,0,326,-122,236, $58,-23,-38,-119,-125,-62,43,61,0,55,0,-10,96, $-233,11,-46,-22,44,18,-101,-57,73,0,-54,-49,2, $-14,29,-13,-37,4,-16,28,19,-16,6,-22,11,0,7,8, $-3,-15,-9,6,2,-14,4,5,-7,17,6,-5,8,-19,8,0, $10,-20,1,14,-11,5,12,-3,1,-2,-2,9,2,10,0,-2, $-1,2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2, $1,4,3,0,0,-6/ data y1925/ $-30926,0,-2318,5817,-893,0,2969,-1334,1471,728, $1140,0,-1645,-462,1202,119,881,229,891,0,711,216, $601,-163,-426,-130,217,-70,-230,0,326,-96,226,58, $-28,-44,-125,-122,-69,51,61,0,54,3,-9,99,-238, $14,-40,-18,39,13,-103,-52,73,0,-54,-50,3,-14, $27,-14,-35,5,-14,29,19,-17,6,-21,11,0,7,8,-3, $-15,-9,6,2,-14,4,5,-7,17,7,-5,8,-19,8,0,10, $-20,1,14,-11,5,12,-3,1,-2,-2,9,2,10,0,-2,-1, $2,-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,1, $4,3,0,0,-6/ data y1930/ $-30805,0,-2316,5808,-951,0,2980,-1424,1517,644, $1172,0,-1692,-480,1205,133,907,166,896,0,727,205, $584,-195,-422,-109,234,-90,-237,0,327,-72,218,60, $-32,-53,-131,-118,-74,58,60,0,53,4,-9,102,-242, $19,-32,-16,32,8,-104,-46,74,0,-54,-51,4,-15,25, $-14,-34,6,-12,29,18,-18,6,-20,11,0,7,8,-3,-15, $-9,5,2,-14,5,5,-6,18,8,-5,8,-19,8,0,10,-20, $1,14,-12,5,12,-3,1,-2,-2,9,3,10,0,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-2,1,4, $3,0,0,-6/ data y1935/ $-30715,0,-2306,5812,-1018,0,2984,-1520,1550,586, $1206,0,-1740,-494,1215,146,918,101,903,0,744,188, $565,-226,-415,-90,249,-114,-241,0,329,-51,211,64, $-33,-64,-136,-115,-76,64,59,0,53,4,-8,104,-246, $25,-25,-15,25,4,-106,-40,74,0,-53,-52,4,-17,23, $-14,-33,7,-11,29,18,-19,6,-19,11,0,7,8,-3,-15, $-9,5,1,-15,6,5,-6,18,8,-5,7,-19,8,0,10,-20, $1,15,-12,5,11,-3,1,-3,-2,9,3,11,0,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-1,2,4, $3,0,0,-6/ data y1940/ $-30654,0,-2292,5821,-1106,0,2981,-1614,1566,528, $1240,0,-1790,-499,1232,163,916,43,914,0,762,169, $550,-252,-405,-72,265,-141,-241,0,334,-33,208,71, $-33,-75,-141,-113,-76,69,57,0,54,4,-7,105,-249, $33,-18,-15,18,0,-107,-33,74,0,-53,-52,4,-18,20, $-14,-31,7,-9,29,17,-20,5,-19,11,0,7,8,-3,-14, $-10,5,1,-15,6,5,-5,19,9,-5,7,-19,8,0,10,-21, $1,15,-12,5,11,-3,1,-3,-2,9,3,11,1,-2,-2,2, $-3,0,-4,2,2,1,-5,2,-2,6,6,-4,4,0,0,-1,2,4, $3,0,0,-6/ data y1945/ $-30594,0,-2285,5810,-1244,0,2990,-1702,1578,477, $1282,0,-1834,-499,1255,186,913,-11,944,0,776,144, $544,-276,-421,-55,304,-178,-253,0,346,-12,194,95, $-20,-67,-142,-119,-82,82,59,0,57,6,6,100,-246, $16,-25,-9,21,-16,-104,-39,70,0,-40,-45,0,-18,0, $2,-29,6,-10,28,15,-17,29,-22,13,0,7,12,-8,-21, $-5,-12,9,-7,7,2,-10,18,7,3,2,-11,5,0,-21,-27, $1,17,-11,29,3,-9,16,4,-3,9,-4,6,-3,1,-4,8, $-3,0,11,5,1,1,2,-20,-5,-1,-1,-6,8,6,-1,-4, $-3,-2,5,0,-2,-2/ data y1950/ $-30554,0,-2250,5815,-1341,0,2998,-1810,1576,381, $1297,0,-1889,-476,1274,206,896,-46,954,0,792,136, $528,-278,-408,-37,303,-210,-240,0,349,3,211,103, $-20,-87,-147,-122,-76,80,54,0,57,-1,4,99,-247, $33,-16,-12,12,-12,-105,-30,65,0,-55,-35,2,-17, $1,0,-40,10,-7,36,5,-18,19,-16,22,0,15,5,-4, $-22,-1,0,11,-21,15,-8,-13,17,5,-4,-1,-17,3,0, $-7,-24,-1,19,-25,12,10,2,5,2,-5,8,-2,8,3,-11, $8,-7,-8,0,4,13,-1,-2,13,-10,-4,2,4,-3,12,6, $3,-3,2,6,10,11,3,8/ data y1955/ $-30500,0,-2215,5820,-1440,0,3003,-1898,1581,291, $1302,0,-1944,-462,1288,216,882,-83,958,0,796,133, $510,-274,-397,-23,290,-230,-229,0,360,15,230,110, $-23,-98,-152,-121,-69,78,47,0,57,-9,3,96,-247, $48,-8,-16,7,-12,-107,-24,65,0,-56,-50,2,-24,10, $-4,-32,8,-11,28,9,-20,18,-18,11,0,9,10,-6,-15, $-14,5,6,-23,10,3,-7,23,6,-4,9,-13,4,0,9,-11, $-4,12,-5,7,2,6,4,-2,1,10,2,7,2,-6,5,5,-3,0, $-5,-4,-1,0,2,-8,-3,-2,7,-4,4,1,-2,-3,6,7,-2, $-1,0,-3/ data y1960/ $-30421,0,-2169,5791,-1555,0,3002,-1967,1590,206, $1302,0,-1992,-414,1289,224,878,-130,957,0,800, $135,504,-278,-394,3,269,-255,-222,0,362,16,242, $125,-26,-117,-156,-114,-63,81,46,0,58,-10,1,99, $-237,60,-1,-20,-2,-11,-113,-17,67,0,-56,-55,5, $-28,15,-6,-32,7,-7,23,17,-18,8,-17,15,0,6,11, $-4,-14,-11,7,2,-18,10,4,-5,23,10,1,8,-20,4,0, $6,-18,0,12,-9,2,1,0,4,-3,-1,9,-2,8,3,0,-1, $5,1,0,-3,4,4,1,0,0,-1,2,4,-5,6,1,1,-1,-1, $6,2,0,0,-7/ data y1965/ $-30334,0,-2119,5776,-1662,0,2997,-2016,1594,114, $1297,0,-2038,-404,1292,240,856,-165,957,0,804, $148,479,-269,-390,13,252,-269,-219,0,358,19,254, $128,-31,-126,-157,-97,-62,81,45,0,61,-11,8,100, $-228,68,4,-32,1,-8,-111,-7,75,0,-57,-61,4,-27, $13,-2,-26,6,-6,26,13,-23,1,-12,13,0,5,7,-4, $-12,-14,9,0,-16,8,4,-1,24,11,-3,4,-17,8,0,10, $-22,2,15,-13,7,10,-4,-1,-5,-1,10,5,10,1,-4, $-2,1,-2,0,-3,2,2,1,-5,2,-2,6,4,-4,4,0,0,-2, $2,3,2,0,0,-6/ data y1970/ $-30220,0,-2068,5737,-1781,0,3000,-2047,1611,25, $1287,0,-2091,-366,1278,251,838,-196,952,0,800, $167,461,-266,-395,26,234,-279,-216,0,359,26,262, $139,-42,-139,-160,-91,-56,83,43,0,64,-12,15,100, $-212,72,2,-37,3,-6,-112,1,72,0,-57,-70,1,-27, $14,-4,-22,8,-2,23,13,-23,-2,-11,14,0,6,7,-2, $-15,-13,6,-3,-17,5,6,0,21,11,-6,3,-16,8,0,10, $-21,2,16,-12,6,10,-4,-1,-5,0,10,3,11,1,-2,-1, $1,-3,0,-3,1,2,1,-5,3,-1,4,6,-4,4,0,1,-1,0, $3,3,1,-1,-4/ data y1975/ $-30100,0,-2013,5675,-1902,0,3010,-2067,1632,-68, $1276,0,-2144,-333,1260,262,830,-223,946,0,791, $191,438,-265,-405,39,216,-288,-218,0,356,31,264, $148,-59,-152,-159,-83,-49,88,45,0,66,-13,28,99, $-198,75,1,-41,6,-4,-111,11,71,0,-56,-77,1,-26, $16,-5,-14,10,0,22,12,-23,-5,-12,14,0,6,6,-1, $-16,-12,4,-8,-19,4,6,0,18,10,-10,1,-17,7,0, $10,-21,2,16,-12,7,10,-4,-1,-5,-1,10,4,11,1, $-3,-2,1,-3,0,-3,1,2,1,-5,3,-2,4,5,-4,4,-1, $1,-1,0,3,3,1,-1,-5/ data y1980/ $-29992,0,-1956,5604,-1997,0,3027,-2129,1663,-200, $1281,0,-2180,-336,1251,271,833,-252,938,0,782, $212,398,-257,-419,53,199,-297,-218,0,357,46,261, $150,-74,-151,-162,-78,-48,92,48,0,66,-15,42,93, $-192,71,4,-43,14,-2,-108,17,72,0,-59,-82,2,-27, $21,-5,-12,16,1,18,11,-23,-2,-10,18,0,6,7,0, $-18,-11,4,-7,-22,4,9,3,16,6,-13,-1,-15,5,0, $10,-21,1,16,-12,9,9,-5,-3,-6,-1,9,7,10,2,-6, $-5,2,-4,0,-4,1,2,0,-5,3,-2,6,5,-4,3,0,1,-1, $2,4,3,0,0,-6/ data y1985/-29873,0,-1905,5500,-2072,0,3044,-2197,1687,-306, $1296,0,-2208,-310,1247,284,829,-297,936,0,780,232,361,-249, $-424,69,170,-297,-214,0,355,47,253,150,-93,-154,-164,-75,-46, $95,53,0,65,-16,51,88,-185,69,4,-48,16,-1,-102,21,74,0,-62,-83, $3,-27,24,-2,-6,20,4,17,10,-23,0,-7,21,0,6,8,0,-19,-11,5,-9, $-23,4,11,4,14,4,-15,-4,-11,5,0,10,-21,1,15,-12,9,9,-6,-3,-6, $-1,9,7,9,1,-7,-5,2,-4,0,-4,1,3,0,-5,3,-2,6,5,-4,3,0,1,-1,2,4, $3,0,0,-6/ data y1990/-29775,0,-1848,5406,-2131,0,3059,-2279,1686,-373, $1314,0,-2239,-284,1248,293,802,-352,939,0,780,247,325,-240, $-423,84,141,-299,-214,0,353,46,245,154,-109,-153,-165,-69,-36, $97,61,0,65,-16,59,82,-178,69,3,-52,18,1,-96,24,77,0,-64,-80,2, $-26,26,0,-1,21,5,17,9,-23,0,-4,23,0,5,10,-1,-19,-10,6,-12,-22, $3,12,4,12,2,-16,-6,-10,4,0,9,-20,1,15,-12,11,9,-7,-4,-7,-2,9, $7,8,1,-7,-6,2,-3,0,-4,2,2,1,-5,3,-2,6,4,-4,3,0,1,-2,3,3,3,-1, $0,-6/ data y1995/-29692,0,-1784,5306,-2200,0,3070,-2366,1681,-413, $1335,0,-2267,-262,1249,302,759,-427,940,0,780,262,290,-236,-418, $97,122,-306,-214,0,352,46,235,165,-118,-143,-166,-55,-17,107, $68,0,67,-17,68,72,-170,67,-1,-58,19,1,-93,36,77,0,-72,-69,1, $-25,28,4,5,24,4,17,8,-24,-2,-6,25,0,6,11,-6,-21,-9,8,-14,-23, $9,15,6,11,-5,-16,-7,-4,4,0,9,-20,3,15,-10,12,8,-6,-8,-8,-1,8, $10,5,-2,-8,-8,3,-3,0,-6,1,2,0,-4,4,-1,5,4,-5,2,-1,2,-2,5,1,1, $-2,0,-7/ data y2000/-29619.4,0,-1728.2,5186.1,-2267.7,0,3068.4,-2481.6, $1670.9,-458.0,1339.6,0,-2288.0,-227.6,1252.1,293.4,714.5,-491.1, $932.3,0,786.8,272.6,250.0,-231.9,-403.0,119.8,111.3,-303.8,-218.8, $0,351.4,43.8,222.3,171.9,-130.4,-133.1,-168.6,-39.3,-12.9,106.3, $72.3,0,68.2,-17.4,74.2,63.7,-160.9,65.1,-5.9,-61.2,16.9,0.7,-90.4, $43.8,79.0,0,-74.0,-64.6,0.0,-24.2,33.3,6.2,9.1,24.0,6.9,14.8,7.3, $-25.4,-1.2,-5.8,24.4,0,6.6,11.9,-9.2,-21.5,-7.9,8.5,-16.6,-21.5, $9.1,15.5,7.0,8.9,-7.9,-14.9,-7.0,-2.1,5.0,0,9.4,-19.7,3.0,13.4, $-8.4,12.5,6.3,-6.2,-8.9,-8.4,-1.5,8.4,9.3,3.8,-4.3,-8.2,-8.2, $4.8,-2.6,0,-6.0,1.7,1.7,0.0,-3.1,4.0,-0.5,4.9,3.7,-5.9,1.0,-1.2, $2.0,-2.9,4.2,0.2,0.3,-2.2,-1.1,-7.4,2.7,0,-1.7,0.1,-1.9,1.3,1.5, $-0.9,-0.1,-2.6,0.1,0.9,-0.7,-0.7,0.7,-2.8,1.7,-0.9,0.1,-1.2,1.2, $-1.9,4.0,-0.9,-2.2,0,-0.3,-0.4,0.2,0.3,0.9,2.5,-0.2,-2.6,0.9, $0.7,-0.5,0.3,0.3,0.0,-0.3,0.0,-0.4,0.3,-0.1,-0.9,-0.2,-0.4,-0.4, $0.8,-0.2,0,-0.9,-0.9,0.3,0.2,0.1,1.8,-0.4,-0.4,1.3,-1.0,-0.4, $-0.1,0.7,0.7,-0.4,0.3,0.3,0.6,-0.1,0.3,0.4,-0.2,0.0,-0.5,0.1,-0.9/ data y2005/-29554.63,0,-1669.05,5077.99,-2337.24,0,3047.69, $-2594.50,1657.76,-515.43,1336.30,0,-2305.83,-198.86,1246.39, $269.72,672.51,-524.72,920.55,0,797.96,282.07,210.65,-225.23, $-379.86,145.15,100.00,-305.36,-227.00,0,354.41,42.72,208.95, $180.25,-136.54,-123.45,-168.05,-19.57,-13.55,103.85,73.60, $0,69.56,-20.33,76.74,54.75,-151.34,63.63,-14.58,-63.53,14.58, $0.24,-86.36,50.94,79.88,0,-74.46,-61.14,-1.65,-22.57,38.73,6.82, $12.30,25.35,9.37,10.93,5.42,-26.32,1.94,-4.64,24.80,0,7.62,11.20, $-11.73,-20.88,-6.88,9.83,-18.11,-19.71,10.17,16.22,9.36,7.61, $-11.25,-12.76,-4.87,-0.06,5.58,0,9.76,-20.11,3.58,12.69,-6.94, $12.67,5.01,-6.72,-10.76,-8.16,-1.25,8.10,8.76,2.92,-6.66,-7.73, $-9.22,6.01,-2.17,0,-6.12,2.19,1.42,0.10,-2.35,4.46,-0.15,4.76, $3.06,-6.58,0.29,-1.01,2.06,-3.47,3.77,-0.86,-0.21,-2.31,-2.09, $-7.93,2.95,0,-1.60,0.26,-1.88,1.44,1.44,-0.77,-0.31,-2.27,0.29, $0.90,-0.79,-0.58,0.53,-2.69,1.80,-1.08,0.16,-1.58,0.96,-1.90, $3.99,-1.39,-2.15,0,-0.29,-0.55,0.21,0.23,0.89,2.38,-0.38,-2.63, $0.96,0.61,-0.30,0.40,0.46,0.01,-0.35,0.02,-0.36,0.28,0.08,-0.87, $-0.49,-0.34,-0.08,0.88,-0.16,0,-0.88,-0.76,0.30,0.33,0.28,1.72, $-0.43,-0.54,1.18,-1.07,-0.37,-0.04,0.75,0.63,-0.26,0.21,0.35, $0.53,-0.05,0.38,0.41,-0.22,-0.10,-0.57,-0.18,-0.82/ data y2010/-29496.5,0,-1585.9,4945.1,-2396.6,0,3026.0,-2707.7, $1668.6,-575.4,1339.7,0,-2326.3,-160.5,1231.7,251.7,634.2,-536.8, $912.6,0,809.0,286.4,166.6,-211.2,-357.1,164.4,89.7,-309.2,-231.1, $0,357.2,44.7,200.3,188.9,-141.2,-118.1,-163.1,0.1,-7.7,100.9, $72.8,0,68.6,-20.8,76.0,44.2,-141.4,61.5,-22.9,-66.3,13.1,3.1, $-77.9,54.9,80.4,0,-75.0,-57.8,-4.7,-21.2,45.3,6.6,14.0,24.9, $10.4,7.0,1.6,-27.7,4.9,-3.4,24.3,0,8.2,10.9,-14.5,-20.0,-5.7, $11.9,-19.3,-17.4,11.6,16.7,10.9,7.1,-14.1,-10.8,-3.7,1.7,5.4, $0,9.4,-20.5,3.4,11.6,-5.3,12.8,3.1,-7.2,-12.4,-7.4,-0.8,8.0, $8.4,2.2,-8.4,-6.1,-10.1,7.0,-2.0,0,-6.3,2.8,0.9,-0.1,-1.1,4.7, $-0.2,4.4,2.5,-7.2,-0.3,-1.0,2.2,-4.0,3.1,-2.0,-1.0,-2.0,-2.8, $-8.3,3.0,0,-1.5,0.1,-2.1,1.7,1.6,-0.6,-0.5,-1.8,0.5,0.9,-0.8, $-0.4,0.4,-2.5,1.8,-1.3,0.2,-2.1,0.8,-1.9,3.8,-1.8,-2.1,0,-0.2, $-0.8,0.3,0.3,1.0,2.2,-0.7,-2.5,0.9,0.5,-0.1,0.6,0.5,0.0,-0.4, $0.1,-0.4,0.3,0.2,-0.9,-0.8,-0.2,0.0,0.8,-0.2,0,-0.9,-0.8,0.3, $0.3,0.4,1.7,-0.4,-0.6,1.1,-1.2,-0.3,-0.1,0.8,0.5,-0.2,0.1,0.4, $0.5,0.0,0.4,0.4,-0.2,-0.3,-0.5,-0.3,-0.8/ end c______________________________________________________________ c======================================================================= c UNIT K: DECODING ROUTINES FOR COMMAND KIT c======================================================================= subroutine getarr(code, values, nwant, nfound) c$$$$ calls nothing c c Extracts an array of numbers from the input store. That store is c a large array in common /store/ which has been filled earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c values the real output array of values found. c nwant the maximum number of numbers expected . c nfound the number of numbers actually found in the input. c If the line contains fewer than nwant values, this is the c value returned in nfound. If an error is discovered c nfound=-n, where n is the number of numbers properly c decoded. If there are no numbers after the codeword c nfound=0. Finally, if the code is absent from the store c nfound=-99 and the array is left undisturbed. c parameter (inmx=200) character *80 input(inmx),line common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) dimension values(*) character local*40, code*4, char*80 c c Read the store in normal order ce do 1010 lin=nin, 1, -1 line=input(lin) c Check code and if line is fresh or old if (code .eq. line(1:4)) then if (iecho .ge. 1) write(*,'(2a)')'==> ',line c Increment read count; copy tail into char, discarding stuff after % kread(lin)=kread(lin) + 1 n1=index(line, ' ')+1 n2=index(line, '%')-1 if (n2 .le. 0) n2=80 char=line(n1:n2) kn=n2 - n1 + 1 goto 1020 endif 1010 continue c Code word not found nfound=-99 return c 1020 continue k1=1 c Up to nwant numbers are sought do 1800 nval=1, nwant do 1100 k=k1, kn if (char(k:k) .ne. ' ') goto 1200 1100 continue nfound=nval-1 return 1200 do 1300 l=k, kn if (char(l:l).eq. ',' .or. char(l:l).eq.' ') goto 1500 1300 continue 1500 if (index(' Ee+-.0123456789', char(k:k)) .eq. 0) goto 2000 if (k-l+1 .gt. 40) goto 2010 k2=l+1 local=char(k:l-1) read (local, '(f40.0)', err=2000) values(nval) 1800 k1=k2 nval=1 - nwant 1900 nfound=1 - nval return c 2000 write(*,*) $'>>> Unreadable numbers in this input line:',line goto 1900 2010 write(*,*) $'>>> Program cannot read a number of more than 40 characters:', $char(k:l-1) goto 1900 end c______________________________________________________________ subroutine getone(code, value, nfound) c$$$$ calls getarr c c Extracts a single number from the input store. That store is c a large array in common /store/ which has been filled earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned for numbers. c value the real output variable containing the desired number. c nfound is 1 if a number is successfully read in; it is zero c the number is absent or unreadable. nfound = -99 if the c code is absent from the input store. c character*4 code dimension v(1) ce call getarr(code, v, 1, nfound) if (nfound .eq. 1) value=v(1) return end c______________________________________________________________ subroutine getchr(code, char, nfound) c$$$$ calls nothing c Extracts a single character variable from the input store. That c store is a large array in common /store/ which has been filled c earlier c c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned. c char the character output variable containing the desired string. c nfound is 1 if a string is successfully read in; it is zero c the line is blank after the code. nfound = -99 if the c code is absent from the input store. c parameter (inmx=200) character *80 input(inmx),line, char*(*), code*4 common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) c Inspect the store in normal order ce do 1010 lin=nin, 1, -1 line=input(lin) c Check code and if line is fresh or old if (code .eq. line(1:4)) then if (iecho .ge. 1) write(*,'(2a)')'==> ',line c Increment read count kread(lin)=kread(lin) + 1 goto 1020 endif 1010 continue c Code word not found nfound=-99 return c 1020 continue c n1=index(line, ' ')+1 n2=index(line, '%')-1 if (n2 .le. 0) n2=80 c Save in char everything after 1st blank to just before % do 1200 k=n1, n2 if (line(k:k) .ne. ' ') then char=line(k:n2) nfound=1 return endif 1200 continue c c Blank field after code word nfound=0 return end c______________________________________________________________ subroutine clear(code) c$$$$ calls nothing c c Removes all command entries identified by code. c code A 4-bye identifying code. Only lines in the input store c beginning with this code are scanned. c if byte 5 is '1' clear only one command not all c parameter (inmx=200) character *80 input(inmx),line, code*5 common /store/ input common /ndict/ iecho,nin,memory,kread(inmx) c c Read the store in normal order ce do 1010 lin= nin, 1, -1 line=input(lin) c Check code and clear the line if present through stack if (code(1:4) .eq. line(1:4)) then input(lin)=' cleared' if (code(5:5).eq.'1') return endif 1010 continue return end c______________________________________________________________ c======================================================================= c UNIT Y: GENERATE AND SUM Ylm c======================================================================= subroutine shsum(ra, lmx,blm, br,bth,bla,V) c$$$$ calls allm c Evaluates a vector-valued function -grad V, and V, where V is a c harmonic function with internal (lmx>0) or external (lmx<0) c sources and has complex spherical harmonic coefficients blm, c fully normalized. Thus c when lmx>0: V = sum{l,m} [b_lm*(1/r)**(l+1) * Ylm] c when lmx<0: V = sum{l,m} [b_lm*r**l * Ylm] c The coefficients are stored serially c c ra 3-vector defining obs point c lmx spherical harmonic degree of series. c if lmx < 0 fields are of external origin, and l(max)=-lmx! c blm complex coefficients in a fully normalized spherical harmonic c series for V. Stored serially. c values are given only for 0 .le. m .le. lmax. c br returned radial component of -grad V. c bth returned value of theta (colatitude) component of -grad V c bla returned value of lambda (longitude) component of -grad V c V returned value of scalar potential V c parameter (mxdeg=96) c dimension r(3),ra(3) complex dx,slm,ex(mxdeg), blm(*) double precision plm,dlm,cosco dimension plm(mxdeg+1),dlm(mxdeg+1) c data near/0/ c ce lmax=abs(lmx) if (lmax .ge. mxdeg) then write(*,'(a,i5/a)')'>>> lmax too large for shsum:', lmax, $ '>>> Reduce lmax or recompile shsum with bigger mxdeg' stop endif c Form unit vector rr=sqrt(ra(1)**2 + ra(2)**2 + ra(3)**2) r(1)=ra(1)/rr r(2)=ra(2)/rr r(3)=ra(3)/rr c sinth=sqrt(r(1)**2 + r(2)**2) if (sinth .eq. 0.0) then r(1)=2.0e-6 sinth=r(1) if (near .lt. 1) write(*,'(a)') $ '>>> One or more observer sites exactly at pole:', $ '>>> Magmap moves such sites 7 meters' near=near + 1 elseif (sinth .lt. 2.0e-6) then r(1)=r(1)*2.0e-6/sinth r(2)=r(2)*2.0e-6/sinth sinth=2.0e-6 endif c cosco=r(3) if (abs(r(3)).gt. 0.99) $ cosco=sqrt(1.0d00 - sinth**2)*sign(1.0, r(3)) c c Calculate and save exp(i*m*phi) ex(1)=1.0 dx=cmplx(r(1), r(2))/sinth do 1300 m=1, lmax ex(m+1)=dx*ex(m) 1300 continue c c Sum the series for -grad V and V br =0.0 bth=0.0 bla=0.0 V=0.0 j=0 do 1700 l=0, lmax call allm(l, cosco, plm,dlm) c Chooses external/internal part by sign of lmx n=l + 2 if (lmx.lt. 0) n=1 - l rpow= rr**n do 1600 m=0, l j=j + 1 slm=ex(m+1)*blm(j) / rpow em=m V =V + rr* plm(m+1)*real(slm) br =br + (n - 1.0) * plm(m+1)*real(slm) bth=bth - dlm(m+1)*real(slm) bla=bla - plm(m+1)*real(cmplx(0.0, em)*slm)/sinth 1600 continue 1700 continue c return end c_______________________________________________________________________ subroutine allm(l, amu, plm,dlm) c$$$$ calls nothing c For degree l finds fully-normalized associated Legendre functions c of order m= 0, 1, 2, ... l for the argument amu=cos(theta) and c corresponding theta derivatives (theta=colatitude). c Results of order m held in plm(m+1), dlm(m+1) c c plm(m+1)=P{l,m}(x)=t(m)*(-1)**m*sqrt((2*l+1)*(l-m)!/(l+m)!)*Plm(x) c where c Plm(x)=(1/2**l*l!)*(1-x**2)**(m/2)*(d/dx)**(l+m)(x**2-1)**l c and t(0)=1/sqrt(4*pi), t(m)=1/sqrt(2*pi), m>0. c c Note that 0 =< l =< 100, unless parameter lmx is changed c implicit double precision (a-h, o-z) parameter (lmx=100, nmx=2*lmx, twopi=2*3.14159265358979324d0) dimension plm(l+1),dlm(l+1),rootn(0:nmx),tlm(lmx+1) save init,rootn data init/0/ c c Fill array with square-roots of the integers ce if (init .eq. 0) then do 1000 n=0, nmx en=n rootn(n)=sqrt(en) 1000 continue init=1 endif c c Check degree constraints if (l .gt. lmx) then write(*,*)'>>> Degree too great for ALLM: l >',lmx stop elseif (l .le. 0) then plm(1)=1.0/sqrt(2.0*twopi) return endif c c Produce normalization constant c=(l + 0.5d0)/twopi do 1100 lx=1, l c=c - c/(2*lx) 1100 continue s=sqrt(1.0 - amu**2) twocot =-2.0*amu/s tlm(l+1)= 2.0*sqrt(c)*(-s)**l tlm(l) = tlm(l+1)*twocot*l/sqrt(2d0*l) c Recur downwards to m=0 do 1500 m=l-2, 0, -1 tlm(m+1)=((m+1)*twocot*tlm(m+2) - $ rootn(l+m+2)*rootn(l-m-1)*tlm(m+3))/(rootn(l+m+1)*rootn(l-m)) 1500 continue c c Calculate the derivative wrt theta cot=amu/s do 1600 m=0, l-1 plm(m+1)=tlm(m+1) dlm(m+1)=rootn(l+m+1)*rootn(l-m)*tlm(m+2)+m*cot*tlm(m+1) 1600 continue dlm(l+1)=l*cot*tlm(l+1) plm(l+1)=tlm(l+1) plm(1)=tlm(1)/2.0 dlm(1)=dlm(1)/2.0 return end c_______________________________________________________________________ c======================================================================= c UNIT P: SINGLE POINT EVALUATION c======================================================================= subroutine point c$$$$ calls getarr shsum clear c Calculates a field model at a specified geographical point parameter (mxdeg=96, mxb=(mxdeg+1)*(mxdeg+2)/2, mxa=mxb/10) complex blm,alm common /coeff/ lmax,lext,blm(mxb),alm(mxa) common /prolix/ iprint common /earth/ rad(4) dimension site(3),r(3) character*64 loci c data radian/0.01745329/ c ce a=1000.0*rad(1) b=1000.0*rad(2) if (rad(4) .le. 0.7) then a=1000.0*rad(3) b=a endif if (a .eq. b) then if (iprint.ge.0) write(*,'(a)')' ',' NOTE:', $ ' Ellipticity is always omitted for core evaluations', $ ' Hence geocentric coordinate system is in use',' ' if (iprint.ge.0) $ write(*,'(a,f9.1)')' Evaluation radius (km)',0.001*a*rad(4) endif c call getchr('site', loci, nsite) if (nsite .le. 0) return open(unit=12, file=loci, status='OLD',err=1100) if (iprint.ge.0) $write(*,'(a)')' ',' Evaluation sites read from:',loci nsite=99 goto 1200 1100 continue 1200 do 1800 k=1, 550000 kk=k if (nsite .ne. 99) then call getarr('site', site, 3, nsite) c Clear just one site command call clear('site1') height=0 if (nsite .eq. 3) height=site(3) else read (12,*, end=1810, err=2000) site(1),site(2),site(3) height=site(3) endif if (nsite .le. 0) return if (k .eq. 1 .and.iprint.ge.0) write(*,'(a)') $ ' Site values automatically sent to file fort.9. ', $ ' Each line in the file consists of:', $ ' lat, long, height, X, Y, Z, F, D, I',' ' alat=site(1) along=site(2) c Coordinate transformations from Langel, pp264-267 in Jacobs c "GEOMAGNETISM", 1987 c c Convert geodetic coordinates to geocentric ones sinla=sin(radian*alat) cosla=cos(radian*alat) r0=sqrt((a*cosla)**2 + (b*sinla)**2) x1=(a**2/r0 + height)*cosla y1=(b**2/r0 + height)*sinla r1=sqrt(x1**2 + y1**2) rabs=r1/a cost=y1/r1 sint=x1/r1 c c Find field elements in geocentric frame c Recall rad(3) is the reference radius scale=rad(1)*rabs / rad(3) r(3)=scale*cost r(1)=scale*sint*cos(radian*along) r(2)=scale*sint*sin(radian*along) call shsum(r, lmax, blm, Br,Bt,Bl,V) if (lext.gt. 0) then call shsum(r, -lext, alm, exBr,exBt,exBl,exV) Br=Br + exBr Bt=Bt + exBt Bl=Bl + exBl V =V + exV endif c c Convert field elements to geodetic frame sins=sinla*sint - cos(radian*alat)*cost coss=sqrt(1.0 - sins**2) F=sqrt(Br**2 + Bt**2 + Bl**2) X=-Bt*coss - Br*sins Y=Bl Z=Bt*sins - Br*coss decl=atan2(Y, X+1e-18)/radian ainc=atan2(Z, sqrt(X**2 + Y**2))/radian c c Print it out if (F .lt. 100000.0) then if (iprint.ge.0) $ write(*,'(/(a18,f10.2,a6,f10.2,a6,f10.2,a))') $' Site: lat ',alat,' long',along,' alt',height,' m', $' Elements (nT): X', X,' Y', Y,' Z',Z,' ', $' F', F,' Decl', decl,' Inc',ainc,' ' else if (iprint.ge.0) $ write(*,'(/ a18,f12.2,a6,f12.2,a6,f12.2,a/ $ (+1p,a18,g12.4,a6,g12.4,a6,g12.4,a))') $' Site: lat ',alat,' long',along,' alt',height,' m', $' Elements (nT): X', X,' Y', Y,' Z',Z,' ', $' F', F,' Decl', decl,' Inc',ainc,' ' endif write(8,'(f6.2,f8.2,1p,e11.4,0p,4f9.1,2f8.2)') $ alat,along,height,X,Y,Z,F,decl,ainc 1800 continue write(*,*)'>>> To many sites: N >', kk 1810 close(unit=12) return c c Error out 2000 write(*,*)'>>> Error reading site file at line:',k close (unit=12) return end c_______________________________________________________________________ subroutine plot(m1,n1, iproj, name) c$$$$ calls launch getchr projg c writes a contour or color command file and executes it c Also writes the coast outline file in mapped coordinates c m1,n1 array size; iproj projection type; name filename of array save kgrph character*(*) name character*64 shore, oldsho, script, tscrip, costid*8 character*64 syline, syfile*10 dimension geo(2), u(2) common /scan/ m,n,umin,vmin,du,factor,hash common /world/ c(3217),nworld common /prolix/ iprint parameter (rd=.0174533, deg=rd) parameter (latsp=15,lonsp=30, lato=90-latsp,nlat=1+2*lato, $ nmerid=nlat*360/lonsp, nparal=361*(180/latsp-1) ) data bdum/1.0e8/, u1,u2/2*1e9/ data script/' '/,nshor/0/, oldhsh/0/,idpts/0/ data oldsho/' '/ c c Don't plot the map or create script call getchr('noplot', script, noplot) if (noplot.ge. 0) return c c Run plotfile in script - final act if (m1 .eq. 0) then call launch(kgrph, script) return endif c c Decide if a plot is called for call getchr('color', tscrip, nscr1) call getchr('plot', tscrip, nscr2) c Color map requested (kgrph=1) nscr=max(nscr1, min(0,nscr2)) lk=kgrph if (nscr .ge. 0) kgrph =1 if (nscr .lt. 0) then call getchr('contour', tscrip, nscr) c Contour map requested (kgrph=0) if (nscr .ge. 0) kgrph=0 endif c No plot required if (nscr .lt. 0) return c c if (nscr .eq. 0) tscrip='tmp.plt' c New plotfile requested if (tscrip .ne. script) then c Run preceeding plotfile then open a new plotfile call launch(lk, script) script=tscrip open(unit=10, file=script) endif c c Create a file of mapped coastline or grid to be superimposed c on contour map call getchr('coast', shore, kind) c if (kind .ge. 0) then if (kind .gt. 0) then if (shore .ne. 'grid') $ open(unit=8, file=shore, status='old',err=2200) if (shore .eq. 'grid') kind=2 endif c nshor=nshor+1 oldsho=shore write(costid,'(a,i2)') 'coast',nshor+9 open (unit=11, file=costid) c crit=0.2 if (iproj .eq. 1) crit=0.2/deg npts=0 do 1100 j=1, 100000 c Read external file if (kind .eq. 1) then read (8,*, err=2000, end=1110) geo c Unpack internal coastline elseif (kind.eq. 0) then if (j .gt. nworld) goto 1110 geo(1)=0.1*c(j) geo(2)=1000.0*mod(abs(c(j)), 1.0) c Generate parallels for grid elseif (kind.eq. 2) then geo(1)=-lato + latsp*((j-1)/361) geo(2)=-180.0 + mod(j-1,361) if (j.eq. nparal) kind=3 c Generate meridians for grid elseif (kind.eq. 3) then geo(1)=-lato + mod((j-1-nparal), nlat) geo(2)=-180 + lonsp*((j-1-nparal)/nlat) if (j.eq. nparal+nmerid) goto 1110 endif c Transform lat-long data to x-y plane call projg(iproj, geo, u) dist=abs(u1-u(1)) + abs(u2-u(2)) c If jump is too large insert a pen-up command if (dist .gt. crit .and. $ geo(1).lt.90.0 .and. u1.ne. 99.0) then npts=npts + 1 write(11,'(a)') '99 370' endif npts=npts + 1 write(11,'(2f9.4)') u u1=u(1) u2=u(2) 1100 continue c 1110 if (kind .eq. 1) close (unit=8) close (unit=11) if (iprint.ge.0) write(*,'(a,i6,a)') $ ' A coastline has been added to the map of',npts,' points' endif 1200 continue c c Create files of mapped sample points to be mapped and superimposed do 1400 nfile=1, 100 call getchr('symbols', syline, mysy) call clear( 'symbols') c Delete symbol file if blank field or when projection changes if (mysy.eq.0 .or. (oldhsh*mysy.lt.0 .and. hash.ne.oldhsh)) then write(*,'(/a)')' Symbol data removed from this plot' write(*,*) $'To retain symbol data in new projection repeat symbol command(s)' write(10,'(a)')'symbols' goto 1410 elseif (mysy .lt. 0) then goto 1410 elseif (mysy .gt. 0) then c Create new symbol file mapped to new coordinates isy=index(syline, ' ') close (unit=8) open (unit=8, file=syline(1:isy), status='OLD',err=1390) idpts=idpts+1 write(syfile,'(a,i2)') 'sym',idpts+9 close (unit=11) open (unit=11, file=syfile) write(*,'(2a)')' Transforming coords in file: ',syline(1:isy) c c Read from symbol file, map point and write it out to the new file do 1300 j=1, nworld read (8,*, err=2000, end=1310) geo call projg(iproj, geo, u) write(11,'(2f9.4)') u 1300 continue c 1310 close (unit=8) close (unit=11) write(10,'(3a)') 'symbols ',syfile,syline(isy : 64) endif goto 1400 c 1390 write(*,*)'>>> Unable to open file: ',syline(1:isy) 1400 continue c c Save id of current map projection 1410 oldhsh=hash c c Create a file to draw the map and execute it: kgrph=0 means c a contour map, =1 means color goto (1500, 1600), kgrph+1 c c Write command file for contour c Draw a bounding circle or ellipse on the map 1500 if (iproj .gt. 1 .and. iproj.ne.5) then write(10,'(a)') 'axes -1 1 -1 1','file *','read 181' write(10,'(10f8.4)') $ (cos(rd*j)/factor,sin(rd*j)/factor,j=0,360,2) endif write(10,'(a,a/a,3i5)') 'file ',name,'read ',m,n,1 write(10,'(a,2f9.1)') 'null ',1-bdum/1000, bdum/1000-1 write(10,'(a,4f9.4,a)')'axes ',umin,-umin,vmin,-vmin if (kind .ge. 0) write(10,'(a/a/a,i6)') $'Color 9','file '//costid, 'read ',npts write(10,'(a)') 'axes 0 0 0 0', 'Color 0','letter 0.06' write(10,'(a)') 'width 5','fancy 0','plot ' return c c This is the color map code c Draw a bounding circle or ellipse on the map 1600 if (iproj .gt. 1 .and. iproj.ne.5) $write(10,'(a)') 'border mask' write(10,'(a,a/a,3i5)') 'file ',name,'read ',m,n,1 write(10,'(a,4f9.4,a)')'axes ',umin,-umin,vmin,-vmin,' hidden' if (kind .ge. 0) write(10,'(a,a)') 'lines ',costid write(10,'(a)') 'width 5','outline 0','palette 2', $'plot ' return c c Error return 2000 write(*,'(2a/a,i6)')'>>> Unreadable data in file',shore, $ ' on line ',j return 2200 write(*,*)'>>> Unable to open file:',shore goto 1200 c end c_______________________________________________________________________