C----------------------------------------------------------------------- C C orb.F C C Subroutines dealing with the solar zenith angle and the earth's C orbit. C C---------------------------Subroutines contained----------------------- C C orb_params --- Calculate the orbital parameters for a given C situation/year. C orb_decl ----- Calculate the solar declination angle and C Earth/Sun distance factor for a given C time of the year. C orb_print ---- Print out information on the orbital parameters C to use. C C---------------------------Code history-------------------------------- C C Original version: Erik Kluzek C Date: Oct/1997 C C----------------------------------------------------------------------- C C Version information: C C CVS: $Id: orb.F,v 1.7.4.1 1998/08/17 21:40:32 zender Exp $ C CVS: $Source: /fs/cgd/csm/models/CVS.REPOS/shared/csm_share/Attic/orb.F,v $ C CVS: $Name: ccm3_6_brnchT_crm2_1_2 $ C C----------------------------------------------------------------------- C======================================================================= subroutine orb_params( iyear_AD , eccen , obliq, mvelp, 1 $ obliqr , lambm0 , mvelpp, log_print ) C----------------------------------------------------------------------- C C Calculate earth's orbital parameters using Dave Threshers C formula which came from Berger, Andre. 1978 C "A Simple Algorithm to Compute Long-Term Variations C of Daily Insolation". Contribution 18, Institute of Astronomy and C Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, C Belgium. C C---------------------------Code history-------------------------------- C C Original Author: Erik Kluzek C Date: Oct/97 C C----------------------------------------------------------------------- implicit none C---------------------------- Input Arguments -------------------------- real eccen, ! Earth's orbital eccentricity $ obliq, ! Earth's obliquity in degree's $ mvelp ! Earth's moving vernal equinox longitude integer iyear_AD ! Year to calculate orbit for.. logical log_print! Flag to print-out status information or not. ! (This turns off ALL status printing including) ! (error messages.) C--------------------------- Output Arguments -------------------------- real $ obliqr, ! Earth's obliquity in radians $ lambm0, ! Mean longitude of perihelion at the c ! vernal equinox (radians) $ mvelpp ! Earth's moving vernal equinox longitude c ! of perihelion plus pi (radians) C------------------------------Parameters------------------------------- C C Parameters for calculating earth's orbital characteristics C integer poblen, ! number of elements in the series to calc obliquity $ pecclen, ! number of elements in the series to calc eccentricity $ pmvelen ! number of elements in the series to calc vernal equinox parameter( poblen = 47, pecclen = 19, pmvelen = 78 ) real $ psecdeg, ! arc seconds to degrees conversion $ degrad, ! degrees to radians conversion factor $ obamp(poblen), ! amplitudes for obliquity cosine series $ obrate(poblen), ! rates for obliquity cosine series $ obphas(poblen), ! phases for obliquity cosine series $ ecamp(pecclen), ! amplitudes for eccentricity/fvelp cosine/sine series $ ecrate(pecclen), ! rates for eccentricity/fvelp cosine/sine series $ ecphas(pecclen), ! phases for eccentricity/fvelp cosine/sine series $ mvamp(pmvelen), ! amplitudes for mvelp sine series $ mvrate(pmvelen), ! rates for mvelp sine series $ mvphas(pmvelen), ! phases for mvelp sine series $ yb4_1950AD ! number of years before 1950 AD C parameter(psecdeg = 1./3600.) C C Cosine series data for computation of obliquity: C amplitude (arc seconds), rate (arc seconds/year), phase (degrees). C data obamp /-2462.2214466, -857.3232075, -629.3231835, $ -414.2804924, -311.7632587, 308.9408604, $ -162.5533601, -116.1077911, 101.1189923, $ -67.6856209, 24.9079067, 22.5811241, $ -21.1648355, -15.6549876, 15.3936813, $ 14.6660938, -11.7273029, 10.2742696, $ 6.4914588, 5.8539148, -5.4872205, $ -5.4290191, 5.1609570, 5.0786314, $ -4.0735782, 3.7227167, 3.3971932, $ -2.8347004, -2.6550721, -2.5717867, $ -2.4712188, 2.4625410, 2.2464112, $ -2.0755511, -1.9713669, -1.8813061, $ -1.8468785, 1.8186742, 1.7601888, $ -1.5428851, 1.4738838, -1.4593669, $ 1.4192259, -1.1818980, 1.1756474, $ -1.1316126, 1.0896928/ C data obrate /31.609974, 32.620504, 24.172203, $ 31.983787, 44.828336, 30.973257, $ 43.668246, 32.246691, 30.599444, $ 42.681324, 43.836462, 47.439436, $ 63.219948, 64.230478, 1.010530, $ 7.437771, 55.782177, 0.373813, $ 13.218362, 62.583231, 63.593761, $ 76.438310, 45.815258, 8.448301, $ 56.792707, 49.747842, 12.058272, $ 75.278220, 65.241008, 64.604291, $ 1.647247, 7.811584, 12.207832, $ 63.856665, 56.155990, 77.448840, $ 6.801054, 62.209418, 20.656133, $ 48.344406, 55.145460, 69.000539, $ 11.071350, 74.291298, 11.047742, $ 0.636717, 12.844549/ C data obphas /251.9025, 280.8325, 128.3057, $ 292.7252, 15.3747, 263.7951, $ 308.4258, 240.0099, 222.9725, $ 268.7809, 316.7998, 319.6024, $ 143.8050, 172.7351, 28.9300, $ 123.5968, 20.2082, 40.8226, $ 123.4722, 155.6977, 184.6277, $ 267.2772, 55.0196, 152.5268, $ 49.1382, 204.6609, 56.5233, $ 200.3284, 201.6651, 213.5577, $ 17.0374, 164.4194, 94.5422, $ 131.9124, 61.0309, 296.2073, $ 135.4894, 114.8750, 247.0691, $ 256.6114, 32.1008, 143.6804, $ 16.8784, 160.6835, 27.5932, $ 348.1074, 82.6496/ C C Cosine/sine series data for computation of eccentricity and C fixed vernal equinox longitude of perihelion (fvelp): C amplitude, rate (arc seconds/year), phase (degrees). C data ecamp /0.01860798, 0.01627522, -0.01300660, $ 0.00988829, -0.00336700, 0.00333077, $ -0.00235400, 0.00140015, 0.00100700, $ 0.00085700, 0.00064990, 0.00059900, $ 0.00037800, -0.00033700, 0.00027600, $ 0.00018200, -0.00017400, -0.00012400, $ 0.00001250/ C data ecrate /4.2072050, 7.3460910, 17.8572630, $ 17.2205460, 16.8467330, 5.1990790, $ 18.2310760, 26.2167580, 6.3591690, $ 16.2100160, 3.0651810, 16.5838290, $ 18.4939800, 6.1909530, 18.8677930, $ 17.4255670, 6.1860010, 18.4174410, $ 0.6678630/ C data ecphas /28.620089, 193.788772, 308.307024, $ 320.199637, 279.376984, 87.195000, $ 349.129677, 128.443387, 154.143880, $ 291.269597, 114.860583, 332.092251, $ 296.414411, 145.769910, 337.237063, $ 152.092288, 126.839891, 210.667199, $ 72.108838/ C C Sine series data for computation of moving vernal equinox C longitude of perihelion: C amplitude (arc seconds), rate (arc seconds/year), phase (degrees). C data mvamp /7391.0225890, 2555.1526947, 2022.7629188, $ -1973.6517951, 1240.2321818, 953.8679112, $ -931.7537108, 872.3795383, 606.3544732, $ -496.0274038, 456.9608039, 346.9462320, $ -305.8412902, 249.6173246, -199.1027200, $ 191.0560889, -175.2936572, 165.9068833, $ 161.1285917, 139.7878093, -133.5228399, $ 117.0673811, 104.6907281, 95.3227476, $ 86.7824524, 86.0857729, 70.5893698, $ -69.9719343, -62.5817473, 61.5450059, $ -57.9364011, 57.1899832, -57.0236109, $ -54.2119253, 53.2834147, 52.1223575, $ -49.0059908, -48.3118757, -45.4191685, $ -42.2357920, -34.7971099, 34.4623613, $ -33.8356643, 33.6689362, -31.2521586, $ -30.8798701, 28.4640769, -27.1960802, $ 27.0860736, -26.3437456, 24.7253740, $ 24.6732126, 24.4272733, 24.0127327, $ 21.7150294, -21.5375347, 18.1148363, $ -16.9603104, -16.1765215, 15.5567653, $ 15.4846529, 15.2150632, 14.5047426, $ -14.3873316, 13.1351419, 12.8776311, $ 11.9867234, 11.9385578, 11.7030822, $ 11.6018181, -11.2617293, -10.4664199, $ 10.4333970, -10.2377466, 10.1934446, $ -10.1280191, 10.0289441, -10.0034259/ C data mvrate /31.609974, 32.620504, 24.172203, $ 0.636717, 31.983787, 3.138886, $ 30.973257, 44.828336, 0.991874, $ 0.373813, 43.668246, 32.246691, $ 30.599444, 2.147012, 10.511172, $ 42.681324, 13.650058, 0.986922, $ 9.874455, 13.013341, 0.262904, $ 0.004952, 1.142024, 63.219948, $ 0.205021, 2.151964, 64.230478, $ 43.836462, 47.439436, 1.384343, $ 7.437771, 18.829299, 9.500642, $ 0.431696, 1.160090, 55.782177, $ 12.639528, 1.155138, 0.168216, $ 1.647247, 10.884985, 5.610937, $ 12.658184, 1.010530, 1.983748, $ 14.023871, 0.560178, 1.273434, $ 12.021467, 62.583231, 63.593761, $ 76.438310, 4.280910, 13.218362, $ 17.818769, 8.359495, 56.792707, $ 8.448301, 1.978796, 8.863925, $ 0.186365, 8.996212, 6.771027, $ 45.815258, 12.002811, 75.278220, $ 65.241008, 18.870667, 22.009553, $ 64.604291, 11.498094, 0.578834, $ 9.237738, 49.747842, 2.147012, $ 1.196895, 2.133898, 0.173168/ C data mvphas /251.9025, 280.8325, 128.3057, $ 348.1074, 292.7252, 165.1686, $ 263.7951, 15.3747, 58.5749, $ 40.8226, 308.4258, 240.0099, $ 222.9725, 106.5937, 114.5182, $ 268.7809, 279.6869, 39.6448, $ 126.4108, 291.5795, 307.2848, $ 18.9300, 273.7596, 143.8050, $ 191.8927, 125.5237, 172.7351, $ 316.7998, 319.6024, 69.7526, $ 123.5968, 217.6432, 85.5882, $ 156.2147, 66.9489, 20.2082, $ 250.7568, 48.0188, 8.3739, $ 17.0374, 155.3409, 94.1709, $ 221.1120, 28.9300, 117.1498, $ 320.5095, 262.3602, 336.2148, $ 233.0046, 155.6977, 184.6277, $ 267.2772, 78.9281, 123.4722, $ 188.7132, 180.1364, 49.1382, $ 152.5268, 98.2198, 97.4808, $ 221.5376, 168.2438, 161.1199, $ 55.0196, 262.6495, 200.3284, $ 201.6651, 294.6547, 99.8233, $ 213.5577, 154.1631, 232.7153, $ 138.3034, 204.6609, 106.5938, $ 250.4676, 332.3345, 27.3039/ C C---------------------- Parameters ------------------------------------- #include <orb.h> C---------------------------Local variables----------------------------- C integer i ! Index for series summations real obsum, ! Obliquity series summation $ cossum, ! Cosine series summation for eccentricity/fvelp $ sinsum, ! Sine series summation for eccentricity/fvelp $ fvelp, ! Fixed vernal equinox longitude of perihelion $ mvsum, ! mvelp series summation $ beta, ! Intermediate argument for lambm0 $ years ! Years to time of interest (negative = past; C ! positive = future) real eccen2 ! eccentricity squared real eccen3 ! eccentricity cubed real pi ! pi C----------------------------------------------------------------------- C C radinp and algorithms below will need a degrees to radians conversion C factor. C if ( log_print ) then print *,'(orb_params) Calculate characteristics of the orbit: ' print *,'(orb_params) CVS revision: $Revision: 1.7.4.1 $' print *,'(orb_params) CVS Tag: $Name: ccm3_6_brnchT_crm2_1_2 $' end if pi = 4.*atan(1.) degrad = pi/180. C C Check for flag to use input orbit parameters C if ( iyear_AD .eq. ORB_NOT_YEAR_BASED ) then c c Check input obliq, eccen, and mvelp to ensure reasonable c if( obliq .eq. ORB_UNDEF_REAL )then if ( log_print ) then print *, '(orb_params) Have to specify orbital parameters!:' print *, 'Either set: iyear_AD, or obliq, eccen, and mvelp:' print *, 'iyear_AD is the year to simulate the orbit for ' > , '(ie. 1950): ' print *, 'obliq, eccen, mvelp specify the orbit directly: ' print *, 'The AMIP II settings (for a 1995 orbit) are: ' print *, 'obliq = 23.4441' print *, 'eccen = 0.016715' print *, 'mvelp = 102.7' end if stop 999 else if ( log_print ) then print *, '(orb_params) Use input orbital parameters: ' end if if( (obliq.lt.ORB_OBLIQ_MIN).or.(obliq.gt.ORB_OBLIQ_MAX) ) then if ( log_print ) then print *, '(orb_params): Input obliquity unreasonable: ' > , obliq end if stop 999 end if if( (eccen.lt.ORB_ECCEN_MIN).or.(eccen.gt.ORB_ECCEN_MAX) ) then if ( log_print ) then print *, '(orb_params): Input eccentricity unreasonable: ' > , eccen end if stop 999 end if if( (mvelp.lt.ORB_MVELP_MIN).or.(mvelp.gt.ORB_MVELP_MAX) ) then if ( log_print ) then print *, '(orb_params): Input mvelp unreasonable: ' > , mvelp end if stop 999 end if eccen2 = eccen*eccen eccen3 = eccen2*eccen else C C Otherwise calculate based on years before present C if ( log_print ) then print *, '(orb_params): Calculate orbit for year: ' > , iyear_AD end if yb4_1950AD = 1950.0 - float(iyear_AD) if ( abs(yb4_1950AD) .gt. 1000000.0 )then if ( log_print ) then print *, '(orb_params) orbit only valid for years +-1000000' print *, '(orb_params) Relative to 1950 AD' print *, '(orb_params) # of years before 1950: ', yb4_1950AD print *, '(orb_params) Year to simulate was : ', iyear_AD end if stop 999 end if C C C The following calculates the earth's obliquity, orbital eccentricity C (and various powers of it) and vernal equinox mean longitude of C perihelion for years in the past (future = negative of years past), C using constants (see parameter section) given in the program of: C C Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations C of Daily Insolation. Contribution 18, Institute of Astronomy and C Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium. C C and formulas given in the paper (where less precise constants are also C given): C C Berger, Andre. 1978. Long-Term Variations of Daily Insolation and C Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367 C C The algorithm is valid only to 1,000,000 years past or hence. C For a solution valid to 5-10 million years past see the above author. C Algorithm below is better for years closer to present than is the C 5-10 million year solution. C C Years to time of interest must be negative of years before present C (1950) in formulas that follow. C years = - yb4_1950AD C C In the summations below, cosine or sine arguments, which end up in C degrees, must be converted to radians via multiplication by degrad. C C Summation of cosine series for obliquity (epsilon in Berger 1978) in C degrees. Convert the amplitudes and rates, which are in arc seconds, into C degrees via multiplication by psecdeg (arc seconds to degrees conversion C factor). For obliq, first term is Berger 1978's epsilon star; second C term is series summation in degrees. C obsum = 0.0 do i = 1, poblen obsum = obsum + $ obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + $ obphas(i))*degrad) end do obliq = 23.320556 + obsum C C Summation of cosine and sine series for computation of eccentricity C (eccen; e in Berger 1978) and fixed vernal equinox longitude of perihelion C (fvelp; pi in Berger 1978), which is used for computation of moving vernal C equinox longitude of perihelion. Convert the rates, which are in arc C seconds, into degrees via multiplication by psecdeg. C cossum = 0.0 do i = 1, pecclen cossum = cossum + $ ecamp(i)*cos((ecrate(i)*psecdeg*years + $ ecphas(i))*degrad) end do C sinsum = 0.0 do i = 1, pecclen sinsum = sinsum + $ ecamp(i)*sin((ecrate(i)*psecdeg*years + $ ecphas(i))*degrad) end do C C Use summations to calculate eccentricity C eccen2 = cossum*cossum + sinsum*sinsum eccen = sqrt(eccen2) eccen3 = eccen2*eccen C C A series of cases for fvelp, which is in radians. C if (abs(cossum) .le. 1.0E-8) then if (sinsum .eq. 0.0) then fvelp = 0.0 else if (sinsum .lt. 0.0) then fvelp = 1.5*pi else if (sinsum .gt. 0.0) then fvelp = .5*pi endif else if (cossum .lt. 0.0) then fvelp = atan(sinsum/cossum) + pi else if (cossum .gt. 0.0) then if (sinsum .lt. 0.0) then fvelp = atan(sinsum/cossum) + 2.0*pi else fvelp = atan(sinsum/cossum) endif endif C C Summation of sine series for computation of moving vernal equinox longitude C of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp, C first term is fvelp in degrees; second term is Berger 1978's psi bar times C years and in degrees; third term is Berger 1978's zeta; fourth term is C series summation in degrees. Convert the amplitudes and rates, which are C in arc seconds, into degrees via multiplication by psecdeg. Series summation C plus second and third terms constitute Berger 1978's psi, which is the C general precession. C mvsum = 0.0 do i = 1, pmvelen mvsum = mvsum + $ mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + $ mvphas(i))*degrad) end do mvelp = fvelp/degrad + 50.439273*psecdeg*years + 3.392506 $ + mvsum C C Cases to make sure mvelp is between 0 and 360. C do while (mvelp .lt. 0.0) mvelp = mvelp + 360.0 end do do while (mvelp .ge. 360.0) mvelp = mvelp - 360.0 end do end if ! end of test on whether to calculate or use input orbital params C C Orbit needs the obliquity in radians C obliqr = obliq*degrad C C 180 degrees must be added to mvelp since observations are made from the C earth and the sun is considered (wrongly for the algorithm) to go around C the earth. For a more graphic explanation see Appendix B in: C C A. Berger, M. Loutre and C. Tricot. 1993. Insolation and Earth's Orbital C Periods. J. of Geophysical Research 98:10,341-10,362. C C Additionally, orbit will need this value in radians. So mvelp becomes C mvelpp (mvelp plus pi) C mvelpp = (mvelp + 180.)*degrad C C Set up an argument used several times in lambm0 calculation ahead. C beta = sqrt(1. - eccen2) C C The mean longitude at the vernal equinox (lambda m nought in Berger C 1978; in radians) is calculated from the following formula given in C Berger 1978. At the vernal equinox the true longitude (lambda in Berger C 1978) is 0. C lambm0 = 2.*((.5*eccen + .125*eccen3)*(1. + beta)*sin(mvelpp) $ - .25*eccen2*(.5 + beta)*sin(2.*mvelpp) $ + .125*eccen3*(1./3. + beta)*sin(3.*mvelpp)) C if ( log_print ) then write(6,*)'(orb_params) ' $ , '------ Computed Orbital Parameters ------' write(6,*)'(orb_params) Eccentricity = ', eccen write(6,*)'(orb_params) Obliquity (deg) = ', obliq write(6,*)'(orb_params) Obliquity (rad) = ',obliqr write(6,*)'(orb_params) Long of perh(deg) = ',mvelp write(6,*)'(orb_params) Long of perh(rad) = ',mvelpp write(6,*)'(orb_params) Long at v.e.(rad) = ',lambm0 write(6,*)'(orb_params) ' $ , '-----------------------------------------' end if C C return end C======================================================================= subroutine orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr , 2 $ delta ,eccf) C----------------------------------------------------------------------- C C Compute earth/orbit parameters using formula suggested by C Duane Thresher. C C---------------------------Code history-------------------------------- C C Original version: Erik Kluzek C Date: Oct/1997 C C----------------------------------------------------------------------- implicit none C------------------------------Arguments-------------------------------- C Input arguments C real calday, ! Calendar day, including fraction $ eccen, ! Eccentricity $ obliqr, ! Earth's obliquity in radians $ lambm0, ! Mean longitude of perihelion at the C ! vernal equinox (radians) $ mvelpp ! Earth's moving vernal equinox longitude C ! of perihelion plus pi (radians) C C Output arguments C real delta, ! Solar declination angle in radians $ eccf ! Earth-sun distance factor ( i.e. (1/r)**2 ) C C---------------------------Local variables----------------------------- C real $ pie, ! pi $ dayspy ! days per year data dayspy / 365. / real ve ! Calday of vernal equinox parameter( ve = 80.5 ) ! correct for Jan 1 = calday 1 parameter( pie = 3.14159265358979323846 ) C real lambm, ! Lambda m, earth's mean longitude of perihelion (radians) $ lmm, ! Intermediate argument involving lambm $ lamb, ! Lambda, the earth's longitude of perihelion $ invrho, ! Inverse normalized sun/earth distance $ sinl ! Sine of lmm C C Compute eccentricity factor and solar declination using C day value where a round day (such as 213.0) refers to 0z at C Greenwich longitude. C C Use formulas from Berger, Andre 1978: Long-Term Variations of Daily C Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci. C 35:2362-2367. C C To get the earth's true longitude (position in orbit; lambda in Berger 1978), C which is necessary to find the eccentricity factor and declination, C must first calculate the mean longitude (lambda m in Berger 1978) at C the present day. This is done by adding to lambm0 (the mean longitude C at the vernal equinox, set as March 21 at noon, when lambda = 0; in radians) C an increment (delta lambda m in Berger 1978) that is the number of C days past or before (a negative increment) the vernal equinox divided by C the days in a model year times the 2*pi radians in a complete orbit. C lambm = lambm0 + (calday - ve)*2.*pie/dayspy lmm = lambm - mvelpp C C The earth's true longitude, in radians, is then found from C the formula in Berger 1978: C sinl = sin(lmm) lamb = lambm + eccen*(2.*sinl $ + eccen*(1.25*sin(2.*lmm) $ + eccen*((13.0/12.0)*sin(3.*lmm) - 0.25*sinl))) C C Using the obliquity, eccentricity, moving vernal equinox longitude of C perihelion (plus), and earth's true longitude, the declination (delta) C and the normalized earth/sun distance (rho in Berger 1978; actually inverse C rho will be used), and thus the eccentricity factor (eccf), can be calculated C from formulas given in Berger 1978. C invrho = (1. + eccen*cos(lamb - mvelpp)) $ / (1. - eccen*eccen) C C Set solar declination and eccentricity factor C delta = asin(sin(obliqr)*sin(lamb)) eccf = invrho*invrho C return C end C======================================================================= subroutine orb_print( iyear_AD, eccen, obliq, mvelp ) C----------------------------------------------------------------------- C C Print out the information on the input orbital characteristics C C---------------------------Code history-------------------------------- C C Original version: Erik Kluzek C Date: Oct/1997 C C----------------------------------------------------------------------- implicit none C---------------------------Parameters---------------------------------- #include <orb.h> C---------------------------Input Arguments----------------------------- real eccen ! Earth's eccentricity factor (unitless) (typically 0 to 0.1) real obliq ! Earth's obliquity angle (degree's) (-90 to +90) (typically 22-26) real mvelp ! Earth's moving vernal equinox at perhelion (degree's) (0 to 360.0) integer iyear_AD ! Year (AD) to simulate above earth's orbital parameters for C----------------------------------------------------------------------- C if ( iyear_AD .eq. ORB_NOT_YEAR_BASED )then if ( obliq .eq. ORB_UNDEF_REAL )then write(6,*) 'Orbit parameters not set!' else write(6,*) 'Orbital parameters: ' write(6,*) 'Obliquity (degree): ', obliq write(6,*) 'Eccentricity (unitless): ', eccen write(6,*) 'Long. of moving Perhelion (deg): ', mvelp end if else if ( iyear_AD .gt. 0 )then write(6,*) 'Orbital parameters calculated for given year: ' $ , iyear_AD, ' AD' else write(6,*) 'Orbital parameters calculated for given year: ' $ , iyear_AD, ' BC' end if end if c return end C=======================================================================