00001
00002
00003
00004
00005
00006 MODULE shr_orb_mod
00007
00008 use shr_kind_mod
00009 use shr_sys_mod
00010 use shr_const_mod
00011 use shr_log_mod, only: s_loglev => shr_log_Level
00012 use shr_log_mod, only: s_logunit => shr_log_Unit
00013
00014 IMPLICIT none
00015
00016
00017
00018
00019 public :: shr_orb_cosz
00020 public :: shr_orb_params
00021 public :: shr_orb_decl
00022 public :: shr_orb_print
00023
00024 real (SHR_KIND_R8),public,parameter :: SHR_ORB_UNDEF_REAL = 1.e36_SHR_KIND_R8
00025 integer(SHR_KIND_IN),public,parameter :: SHR_ORB_UNDEF_INT = 2000000000
00026
00027
00028
00029
00030 private
00031
00032 real (SHR_KIND_R8),parameter :: pi = SHR_CONST_PI
00033 real (SHR_KIND_R8),parameter :: SHR_ORB_ECCEN_MIN = 0.0_SHR_KIND_R8
00034 real (SHR_KIND_R8),parameter :: SHR_ORB_ECCEN_MAX = 0.1_SHR_KIND_R8
00035 real (SHR_KIND_R8),parameter :: SHR_ORB_OBLIQ_MIN = -90.0_SHR_KIND_R8
00036 real (SHR_KIND_R8),parameter :: SHR_ORB_OBLIQ_MAX = +90.0_SHR_KIND_R8
00037 real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MIN = 0.0_SHR_KIND_R8
00038 real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MAX = 360.0_SHR_KIND_R8
00039
00040
00041 CONTAINS
00042
00043
00044 real(SHR_KIND_R8) FUNCTION shr_orb_cosz(jday,lat,lon,declin)
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 real (SHR_KIND_R8),intent(in) :: jday
00060 real (SHR_KIND_R8),intent(in) :: lat
00061 real (SHR_KIND_R8),intent(in) :: lon
00062 real (SHR_KIND_R8),intent(in) :: declin
00063
00064
00065
00066 shr_orb_cosz = sin(lat)*sin(declin) - &
00067 & cos(lat)*cos(declin)*cos(jday*2.0_SHR_KIND_R8*pi + lon)
00068
00069 END FUNCTION shr_orb_cosz
00070
00071
00072
00073 SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , &
00074 & obliqr , lambm0 , mvelpp, log_print )
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 real (SHR_KIND_R8),intent(inout) :: eccen
00092 real (SHR_KIND_R8),intent(inout) :: obliq
00093 real (SHR_KIND_R8),intent(inout) :: mvelp
00094 integer(SHR_KIND_IN),intent(in) :: iyear_AD
00095 logical ,intent(in) :: log_print
00096 real (SHR_KIND_R8),intent(out) :: obliqr
00097 real (SHR_KIND_R8),intent(out) :: lambm0
00098
00099 real (SHR_KIND_R8),intent(out) :: mvelpp
00100
00101
00102
00103 integer(SHR_KIND_IN),parameter :: poblen =47
00104 integer(SHR_KIND_IN),parameter :: pecclen=19
00105 integer(SHR_KIND_IN),parameter :: pmvelen=78
00106 real (SHR_KIND_R8),parameter :: psecdeg = 1.0_SHR_KIND_R8/3600.0_SHR_KIND_R8
00107
00108 real (SHR_KIND_R8) :: degrad = pi/180._SHR_KIND_R8
00109 real (SHR_KIND_R8) :: yb4_1950AD
00110
00111
00112
00113
00114 real (SHR_KIND_R8), parameter :: obamp(poblen) =
00115 (/ -2462.2214466_SHR_KIND_R8, -857.3232075_SHR_KIND_R8, -629.3231835_SHR_KIND_R8, &
00116 -414.2804924_SHR_KIND_R8, -311.7632587_SHR_KIND_R8, 308.9408604_SHR_KIND_R8, &
00117 -162.5533601_SHR_KIND_R8, -116.1077911_SHR_KIND_R8, 101.1189923_SHR_KIND_R8, &
00118 -67.6856209_SHR_KIND_R8, 24.9079067_SHR_KIND_R8, 22.5811241_SHR_KIND_R8, &
00119 -21.1648355_SHR_KIND_R8, -15.6549876_SHR_KIND_R8, 15.3936813_SHR_KIND_R8, &
00120 14.6660938_SHR_KIND_R8, -11.7273029_SHR_KIND_R8, 10.2742696_SHR_KIND_R8, &
00121 6.4914588_SHR_KIND_R8, 5.8539148_SHR_KIND_R8, -5.4872205_SHR_KIND_R8, &
00122 -5.4290191_SHR_KIND_R8, 5.1609570_SHR_KIND_R8, 5.0786314_SHR_KIND_R8, &
00123 -4.0735782_SHR_KIND_R8, 3.7227167_SHR_KIND_R8, 3.3971932_SHR_KIND_R8, &
00124 -2.8347004_SHR_KIND_R8, -2.6550721_SHR_KIND_R8, -2.5717867_SHR_KIND_R8, &
00125 -2.4712188_SHR_KIND_R8, 2.4625410_SHR_KIND_R8, 2.2464112_SHR_KIND_R8, &
00126 -2.0755511_SHR_KIND_R8, -1.9713669_SHR_KIND_R8, -1.8813061_SHR_KIND_R8, &
00127 -1.8468785_SHR_KIND_R8, 1.8186742_SHR_KIND_R8, 1.7601888_SHR_KIND_R8, &
00128 -1.5428851_SHR_KIND_R8, 1.4738838_SHR_KIND_R8, -1.4593669_SHR_KIND_R8, &
00129 1.4192259_SHR_KIND_R8, -1.1818980_SHR_KIND_R8, 1.1756474_SHR_KIND_R8, &
00130 -1.1316126_SHR_KIND_R8, 1.0896928_SHR_KIND_R8/)
00131
00132 real (SHR_KIND_R8), parameter :: obrate(poblen) =
00133 (/ 31.609974_SHR_KIND_R8, 32.620504_SHR_KIND_R8, 24.172203_SHR_KIND_R8, &
00134 31.983787_SHR_KIND_R8, 44.828336_SHR_KIND_R8, 30.973257_SHR_KIND_R8, &
00135 43.668246_SHR_KIND_R8, 32.246691_SHR_KIND_R8, 30.599444_SHR_KIND_R8, &
00136 42.681324_SHR_KIND_R8, 43.836462_SHR_KIND_R8, 47.439436_SHR_KIND_R8, &
00137 63.219948_SHR_KIND_R8, 64.230478_SHR_KIND_R8, 1.010530_SHR_KIND_R8, &
00138 7.437771_SHR_KIND_R8, 55.782177_SHR_KIND_R8, 0.373813_SHR_KIND_R8, &
00139 13.218362_SHR_KIND_R8, 62.583231_SHR_KIND_R8, 63.593761_SHR_KIND_R8, &
00140 76.438310_SHR_KIND_R8, 45.815258_SHR_KIND_R8, 8.448301_SHR_KIND_R8, &
00141 56.792707_SHR_KIND_R8, 49.747842_SHR_KIND_R8, 12.058272_SHR_KIND_R8, &
00142 75.278220_SHR_KIND_R8, 65.241008_SHR_KIND_R8, 64.604291_SHR_KIND_R8, &
00143 1.647247_SHR_KIND_R8, 7.811584_SHR_KIND_R8, 12.207832_SHR_KIND_R8, &
00144 63.856665_SHR_KIND_R8, 56.155990_SHR_KIND_R8, 77.448840_SHR_KIND_R8, &
00145 6.801054_SHR_KIND_R8, 62.209418_SHR_KIND_R8, 20.656133_SHR_KIND_R8, &
00146 48.344406_SHR_KIND_R8, 55.145460_SHR_KIND_R8, 69.000539_SHR_KIND_R8, &
00147 11.071350_SHR_KIND_R8, 74.291298_SHR_KIND_R8, 11.047742_SHR_KIND_R8, &
00148 0.636717_SHR_KIND_R8, 12.844549_SHR_KIND_R8/)
00149
00150 real (SHR_KIND_R8), parameter :: obphas(poblen) =
00151 (/ 251.9025_SHR_KIND_R8, 280.8325_SHR_KIND_R8, 128.3057_SHR_KIND_R8, &
00152 292.7252_SHR_KIND_R8, 15.3747_SHR_KIND_R8, 263.7951_SHR_KIND_R8, &
00153 308.4258_SHR_KIND_R8, 240.0099_SHR_KIND_R8, 222.9725_SHR_KIND_R8, &
00154 268.7809_SHR_KIND_R8, 316.7998_SHR_KIND_R8, 319.6024_SHR_KIND_R8, &
00155 143.8050_SHR_KIND_R8, 172.7351_SHR_KIND_R8, 28.9300_SHR_KIND_R8, &
00156 123.5968_SHR_KIND_R8, 20.2082_SHR_KIND_R8, 40.8226_SHR_KIND_R8, &
00157 123.4722_SHR_KIND_R8, 155.6977_SHR_KIND_R8, 184.6277_SHR_KIND_R8, &
00158 267.2772_SHR_KIND_R8, 55.0196_SHR_KIND_R8, 152.5268_SHR_KIND_R8, &
00159 49.1382_SHR_KIND_R8, 204.6609_SHR_KIND_R8, 56.5233_SHR_KIND_R8, &
00160 200.3284_SHR_KIND_R8, 201.6651_SHR_KIND_R8, 213.5577_SHR_KIND_R8, &
00161 17.0374_SHR_KIND_R8, 164.4194_SHR_KIND_R8, 94.5422_SHR_KIND_R8, &
00162 131.9124_SHR_KIND_R8, 61.0309_SHR_KIND_R8, 296.2073_SHR_KIND_R8, &
00163 135.4894_SHR_KIND_R8, 114.8750_SHR_KIND_R8, 247.0691_SHR_KIND_R8, &
00164 256.6114_SHR_KIND_R8, 32.1008_SHR_KIND_R8, 143.6804_SHR_KIND_R8, &
00165 16.8784_SHR_KIND_R8, 160.6835_SHR_KIND_R8, 27.5932_SHR_KIND_R8, &
00166 348.1074_SHR_KIND_R8, 82.6496_SHR_KIND_R8/)
00167
00168
00169
00170
00171
00172 real (SHR_KIND_R8), parameter :: ecamp (pecclen) =
00173 (/ 0.01860798_SHR_KIND_R8, 0.01627522_SHR_KIND_R8, -0.01300660_SHR_KIND_R8, &
00174 0.00988829_SHR_KIND_R8, -0.00336700_SHR_KIND_R8, 0.00333077_SHR_KIND_R8, &
00175 -0.00235400_SHR_KIND_R8, 0.00140015_SHR_KIND_R8, 0.00100700_SHR_KIND_R8, &
00176 0.00085700_SHR_KIND_R8, 0.00064990_SHR_KIND_R8, 0.00059900_SHR_KIND_R8, &
00177 0.00037800_SHR_KIND_R8, -0.00033700_SHR_KIND_R8, 0.00027600_SHR_KIND_R8, &
00178 0.00018200_SHR_KIND_R8, -0.00017400_SHR_KIND_R8, -0.00012400_SHR_KIND_R8, &
00179 0.00001250_SHR_KIND_R8/)
00180
00181 real (SHR_KIND_R8), parameter :: ecrate(pecclen) =
00182 (/ 4.2072050_SHR_KIND_R8, 7.3460910_SHR_KIND_R8, 17.8572630_SHR_KIND_R8, &
00183 17.2205460_SHR_KIND_R8, 16.8467330_SHR_KIND_R8, 5.1990790_SHR_KIND_R8, &
00184 18.2310760_SHR_KIND_R8, 26.2167580_SHR_KIND_R8, 6.3591690_SHR_KIND_R8, &
00185 16.2100160_SHR_KIND_R8, 3.0651810_SHR_KIND_R8, 16.5838290_SHR_KIND_R8, &
00186 18.4939800_SHR_KIND_R8, 6.1909530_SHR_KIND_R8, 18.8677930_SHR_KIND_R8, &
00187 17.4255670_SHR_KIND_R8, 6.1860010_SHR_KIND_R8, 18.4174410_SHR_KIND_R8, &
00188 0.6678630_SHR_KIND_R8/)
00189
00190 real (SHR_KIND_R8), parameter :: ecphas(pecclen) =
00191 (/ 28.620089_SHR_KIND_R8, 193.788772_SHR_KIND_R8, 308.307024_SHR_KIND_R8, &
00192 320.199637_SHR_KIND_R8, 279.376984_SHR_KIND_R8, 87.195000_SHR_KIND_R8, &
00193 349.129677_SHR_KIND_R8, 128.443387_SHR_KIND_R8, 154.143880_SHR_KIND_R8, &
00194 291.269597_SHR_KIND_R8, 114.860583_SHR_KIND_R8, 332.092251_SHR_KIND_R8, &
00195 296.414411_SHR_KIND_R8, 145.769910_SHR_KIND_R8, 337.237063_SHR_KIND_R8, &
00196 152.092288_SHR_KIND_R8, 126.839891_SHR_KIND_R8, 210.667199_SHR_KIND_R8, &
00197 72.108838_SHR_KIND_R8/)
00198
00199
00200
00201
00202 real (SHR_KIND_R8), parameter :: mvamp (pmvelen) =
00203 (/ 7391.0225890_SHR_KIND_R8, 2555.1526947_SHR_KIND_R8, 2022.7629188_SHR_KIND_R8, &
00204 -1973.6517951_SHR_KIND_R8, 1240.2321818_SHR_KIND_R8, 953.8679112_SHR_KIND_R8, &
00205 -931.7537108_SHR_KIND_R8, 872.3795383_SHR_KIND_R8, 606.3544732_SHR_KIND_R8, &
00206 -496.0274038_SHR_KIND_R8, 456.9608039_SHR_KIND_R8, 346.9462320_SHR_KIND_R8, &
00207 -305.8412902_SHR_KIND_R8, 249.6173246_SHR_KIND_R8, -199.1027200_SHR_KIND_R8, &
00208 191.0560889_SHR_KIND_R8, -175.2936572_SHR_KIND_R8, 165.9068833_SHR_KIND_R8, &
00209 161.1285917_SHR_KIND_R8, 139.7878093_SHR_KIND_R8, -133.5228399_SHR_KIND_R8, &
00210 117.0673811_SHR_KIND_R8, 104.6907281_SHR_KIND_R8, 95.3227476_SHR_KIND_R8, &
00211 86.7824524_SHR_KIND_R8, 86.0857729_SHR_KIND_R8, 70.5893698_SHR_KIND_R8, &
00212 -69.9719343_SHR_KIND_R8, -62.5817473_SHR_KIND_R8, 61.5450059_SHR_KIND_R8, &
00213 -57.9364011_SHR_KIND_R8, 57.1899832_SHR_KIND_R8, -57.0236109_SHR_KIND_R8, &
00214 -54.2119253_SHR_KIND_R8, 53.2834147_SHR_KIND_R8, 52.1223575_SHR_KIND_R8, &
00215 -49.0059908_SHR_KIND_R8, -48.3118757_SHR_KIND_R8, -45.4191685_SHR_KIND_R8, &
00216 -42.2357920_SHR_KIND_R8, -34.7971099_SHR_KIND_R8, 34.4623613_SHR_KIND_R8, &
00217 -33.8356643_SHR_KIND_R8, 33.6689362_SHR_KIND_R8, -31.2521586_SHR_KIND_R8, &
00218 -30.8798701_SHR_KIND_R8, 28.4640769_SHR_KIND_R8, -27.1960802_SHR_KIND_R8, &
00219 27.0860736_SHR_KIND_R8, -26.3437456_SHR_KIND_R8, 24.7253740_SHR_KIND_R8, &
00220 24.6732126_SHR_KIND_R8, 24.4272733_SHR_KIND_R8, 24.0127327_SHR_KIND_R8, &
00221 21.7150294_SHR_KIND_R8, -21.5375347_SHR_KIND_R8, 18.1148363_SHR_KIND_R8, &
00222 -16.9603104_SHR_KIND_R8, -16.1765215_SHR_KIND_R8, 15.5567653_SHR_KIND_R8, &
00223 15.4846529_SHR_KIND_R8, 15.2150632_SHR_KIND_R8, 14.5047426_SHR_KIND_R8, &
00224 -14.3873316_SHR_KIND_R8, 13.1351419_SHR_KIND_R8, 12.8776311_SHR_KIND_R8, &
00225 11.9867234_SHR_KIND_R8, 11.9385578_SHR_KIND_R8, 11.7030822_SHR_KIND_R8, &
00226 11.6018181_SHR_KIND_R8, -11.2617293_SHR_KIND_R8, -10.4664199_SHR_KIND_R8, &
00227 10.4333970_SHR_KIND_R8, -10.2377466_SHR_KIND_R8, 10.1934446_SHR_KIND_R8, &
00228 -10.1280191_SHR_KIND_R8, 10.0289441_SHR_KIND_R8, -10.0034259_SHR_KIND_R8/)
00229
00230 real (SHR_KIND_R8), parameter :: mvrate(pmvelen) =
00231 (/ 31.609974_SHR_KIND_R8, 32.620504_SHR_KIND_R8, 24.172203_SHR_KIND_R8, &
00232 0.636717_SHR_KIND_R8, 31.983787_SHR_KIND_R8, 3.138886_SHR_KIND_R8, &
00233 30.973257_SHR_KIND_R8, 44.828336_SHR_KIND_R8, 0.991874_SHR_KIND_R8, &
00234 0.373813_SHR_KIND_R8, 43.668246_SHR_KIND_R8, 32.246691_SHR_KIND_R8, &
00235 30.599444_SHR_KIND_R8, 2.147012_SHR_KIND_R8, 10.511172_SHR_KIND_R8, &
00236 42.681324_SHR_KIND_R8, 13.650058_SHR_KIND_R8, 0.986922_SHR_KIND_R8, &
00237 9.874455_SHR_KIND_R8, 13.013341_SHR_KIND_R8, 0.262904_SHR_KIND_R8, &
00238 0.004952_SHR_KIND_R8, 1.142024_SHR_KIND_R8, 63.219948_SHR_KIND_R8, &
00239 0.205021_SHR_KIND_R8, 2.151964_SHR_KIND_R8, 64.230478_SHR_KIND_R8, &
00240 43.836462_SHR_KIND_R8, 47.439436_SHR_KIND_R8, 1.384343_SHR_KIND_R8, &
00241 7.437771_SHR_KIND_R8, 18.829299_SHR_KIND_R8, 9.500642_SHR_KIND_R8, &
00242 0.431696_SHR_KIND_R8, 1.160090_SHR_KIND_R8, 55.782177_SHR_KIND_R8, &
00243 12.639528_SHR_KIND_R8, 1.155138_SHR_KIND_R8, 0.168216_SHR_KIND_R8, &
00244 1.647247_SHR_KIND_R8, 10.884985_SHR_KIND_R8, 5.610937_SHR_KIND_R8, &
00245 12.658184_SHR_KIND_R8, 1.010530_SHR_KIND_R8, 1.983748_SHR_KIND_R8, &
00246 14.023871_SHR_KIND_R8, 0.560178_SHR_KIND_R8, 1.273434_SHR_KIND_R8, &
00247 12.021467_SHR_KIND_R8, 62.583231_SHR_KIND_R8, 63.593761_SHR_KIND_R8, &
00248 76.438310_SHR_KIND_R8, 4.280910_SHR_KIND_R8, 13.218362_SHR_KIND_R8, &
00249 17.818769_SHR_KIND_R8, 8.359495_SHR_KIND_R8, 56.792707_SHR_KIND_R8, &
00250 8.448301_SHR_KIND_R8, 1.978796_SHR_KIND_R8, 8.863925_SHR_KIND_R8, &
00251 0.186365_SHR_KIND_R8, 8.996212_SHR_KIND_R8, 6.771027_SHR_KIND_R8, &
00252 45.815258_SHR_KIND_R8, 12.002811_SHR_KIND_R8, 75.278220_SHR_KIND_R8, &
00253 65.241008_SHR_KIND_R8, 18.870667_SHR_KIND_R8, 22.009553_SHR_KIND_R8, &
00254 64.604291_SHR_KIND_R8, 11.498094_SHR_KIND_R8, 0.578834_SHR_KIND_R8, &
00255 9.237738_SHR_KIND_R8, 49.747842_SHR_KIND_R8, 2.147012_SHR_KIND_R8, &
00256 1.196895_SHR_KIND_R8, 2.133898_SHR_KIND_R8, 0.173168_SHR_KIND_R8/)
00257
00258 real (SHR_KIND_R8), parameter :: mvphas(pmvelen) =
00259 (/ 251.9025_SHR_KIND_R8, 280.8325_SHR_KIND_R8, 128.3057_SHR_KIND_R8, &
00260 348.1074_SHR_KIND_R8, 292.7252_SHR_KIND_R8, 165.1686_SHR_KIND_R8, &
00261 263.7951_SHR_KIND_R8, 15.3747_SHR_KIND_R8, 58.5749_SHR_KIND_R8, &
00262 40.8226_SHR_KIND_R8, 308.4258_SHR_KIND_R8, 240.0099_SHR_KIND_R8, &
00263 222.9725_SHR_KIND_R8, 106.5937_SHR_KIND_R8, 114.5182_SHR_KIND_R8, &
00264 268.7809_SHR_KIND_R8, 279.6869_SHR_KIND_R8, 39.6448_SHR_KIND_R8, &
00265 126.4108_SHR_KIND_R8, 291.5795_SHR_KIND_R8, 307.2848_SHR_KIND_R8, &
00266 18.9300_SHR_KIND_R8, 273.7596_SHR_KIND_R8, 143.8050_SHR_KIND_R8, &
00267 191.8927_SHR_KIND_R8, 125.5237_SHR_KIND_R8, 172.7351_SHR_KIND_R8, &
00268 316.7998_SHR_KIND_R8, 319.6024_SHR_KIND_R8, 69.7526_SHR_KIND_R8, &
00269 123.5968_SHR_KIND_R8, 217.6432_SHR_KIND_R8, 85.5882_SHR_KIND_R8, &
00270 156.2147_SHR_KIND_R8, 66.9489_SHR_KIND_R8, 20.2082_SHR_KIND_R8, &
00271 250.7568_SHR_KIND_R8, 48.0188_SHR_KIND_R8, 8.3739_SHR_KIND_R8, &
00272 17.0374_SHR_KIND_R8, 155.3409_SHR_KIND_R8, 94.1709_SHR_KIND_R8, &
00273 221.1120_SHR_KIND_R8, 28.9300_SHR_KIND_R8, 117.1498_SHR_KIND_R8, &
00274 320.5095_SHR_KIND_R8, 262.3602_SHR_KIND_R8, 336.2148_SHR_KIND_R8, &
00275 233.0046_SHR_KIND_R8, 155.6977_SHR_KIND_R8, 184.6277_SHR_KIND_R8, &
00276 267.2772_SHR_KIND_R8, 78.9281_SHR_KIND_R8, 123.4722_SHR_KIND_R8, &
00277 188.7132_SHR_KIND_R8, 180.1364_SHR_KIND_R8, 49.1382_SHR_KIND_R8, &
00278 152.5268_SHR_KIND_R8, 98.2198_SHR_KIND_R8, 97.4808_SHR_KIND_R8, &
00279 221.5376_SHR_KIND_R8, 168.2438_SHR_KIND_R8, 161.1199_SHR_KIND_R8, &
00280 55.0196_SHR_KIND_R8, 262.6495_SHR_KIND_R8, 200.3284_SHR_KIND_R8, &
00281 201.6651_SHR_KIND_R8, 294.6547_SHR_KIND_R8, 99.8233_SHR_KIND_R8, &
00282 213.5577_SHR_KIND_R8, 154.1631_SHR_KIND_R8, 232.7153_SHR_KIND_R8, &
00283 138.3034_SHR_KIND_R8, 204.6609_SHR_KIND_R8, 106.5938_SHR_KIND_R8, &
00284 250.4676_SHR_KIND_R8, 332.3345_SHR_KIND_R8, 27.3039_SHR_KIND_R8/)
00285
00286
00287 integer(SHR_KIND_IN) :: i
00288 real (SHR_KIND_R8) :: obsum
00289 real (SHR_KIND_R8) :: cossum
00290 real (SHR_KIND_R8) :: sinsum
00291 real (SHR_KIND_R8) :: fvelp
00292 real (SHR_KIND_R8) :: mvsum
00293 real (SHR_KIND_R8) :: beta
00294 real (SHR_KIND_R8) :: years
00295 real (SHR_KIND_R8) :: eccen2
00296 real (SHR_KIND_R8) :: eccen3
00297
00298
00299 character(*),parameter :: svnID = "SVN " //
00300 "$Id: shr_orb_mod.F90 6752 2007-10-04 21:02:15Z jwolfe $"
00301 character(*),parameter :: svnURL = "SVN <unknown URL>"
00302
00303
00304 character(len=*),parameter :: F00 = "('(shr_orb_params) ',4a)"
00305 character(len=*),parameter :: F01 = "('(shr_orb_params) ',a,i9)"
00306 character(len=*),parameter :: F02 = "('(shr_orb_params) ',a,f6.3)"
00307 character(len=*),parameter :: F03 = "('(shr_orb_params) ',a,es14.6)"
00308
00309
00310
00311
00312 if ( log_print .and. s_loglev > 0 ) then
00313 write(s_logunit,F00) 'Calculate characteristics of the orbit:'
00314 write(s_logunit,F00) svnID
00315 write(s_logunit,F00) svnURL
00316 end if
00317
00318
00319
00320 IF ( iyear_AD == SHR_ORB_UNDEF_INT ) THEN
00321
00322
00323
00324 if( obliq == SHR_ORB_UNDEF_REAL )then
00325 write(s_logunit,F00) 'Have to specify orbital parameters:'
00326 write(s_logunit,F00) 'Either set: iyear_AD, OR [obliq, eccen, and mvelp]:'
00327 write(s_logunit,F00) 'iyear_AD is the year to simulate orbit for (ie. 1950): '
00328 write(s_logunit,F00) 'obliq, eccen, mvelp specify the orbit directly:'
00329 write(s_logunit,F00) 'The AMIP II settings (for a 1995 orbit) are: '
00330 write(s_logunit,F00) ' obliq = 23.4441'
00331 write(s_logunit,F00) ' eccen = 0.016715'
00332 write(s_logunit,F00) ' mvelp = 102.7'
00333 call shr_sys_abort()
00334 else if ( log_print ) then
00335 write(s_logunit,F00) 'Use input orbital parameters: '
00336 end if
00337 if( (obliq < SHR_ORB_OBLIQ_MIN).or.(obliq > SHR_ORB_OBLIQ_MAX) ) then
00338 write(s_logunit,F03) 'Input obliquity unreasonable: ', obliq
00339 call shr_sys_abort()
00340 end if
00341 if( (eccen < SHR_ORB_ECCEN_MIN).or.(eccen > SHR_ORB_ECCEN_MAX) ) then
00342 write(s_logunit,F03) 'Input eccentricity unreasonable: ', eccen
00343 call shr_sys_abort()
00344 end if
00345 if( (mvelp < SHR_ORB_MVELP_MIN).or.(mvelp > SHR_ORB_MVELP_MAX) ) then
00346 write(s_logunit,F03) 'Input mvelp unreasonable: ' , mvelp
00347 call shr_sys_abort()
00348 end if
00349 eccen2 = eccen*eccen
00350 eccen3 = eccen2*eccen
00351
00352 ELSE
00353
00354 if ( log_print .and. s_loglev > 0) then
00355 write(s_logunit,F01) 'Calculate orbit for year: ' , iyear_AD
00356 end if
00357 yb4_1950AD = 1950.0_SHR_KIND_R8 - real(iyear_AD,SHR_KIND_R8)
00358 if ( abs(yb4_1950AD) .gt. 1000000.0_SHR_KIND_R8 )then
00359 write(s_logunit,F00) 'orbit only valid for years+-1000000'
00360 write(s_logunit,F00) 'Relative to 1950 AD'
00361 write(s_logunit,F03) '# of years before 1950: ',yb4_1950AD
00362 write(s_logunit,F01) 'Year to simulate was : ',iyear_AD
00363 call shr_sys_abort()
00364 end if
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 years = - yb4_1950AD
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 obsum = 0.0_SHR_KIND_R8
00401 do i = 1, poblen
00402 obsum = obsum + obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + &
00403 & obphas(i))*degrad)
00404 end do
00405 obliq = 23.320556_SHR_KIND_R8 + obsum
00406
00407
00408
00409
00410
00411
00412
00413 cossum = 0.0_SHR_KIND_R8
00414 do i = 1, pecclen
00415 cossum = cossum+ecamp(i)*cos((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
00416 end do
00417
00418 sinsum = 0.0_SHR_KIND_R8
00419 do i = 1, pecclen
00420 sinsum = sinsum+ecamp(i)*sin((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
00421 end do
00422
00423
00424
00425 eccen2 = cossum*cossum + sinsum*sinsum
00426 eccen = sqrt(eccen2)
00427 eccen3 = eccen2*eccen
00428
00429
00430
00431 if (abs(cossum) .le. 1.0E-8_SHR_KIND_R8) then
00432 if (sinsum .eq. 0.0_SHR_KIND_R8) then
00433 fvelp = 0.0_SHR_KIND_R8
00434 else if (sinsum .lt. 0.0_SHR_KIND_R8) then
00435 fvelp = 1.5_SHR_KIND_R8*pi
00436 else if (sinsum .gt. 0.0_SHR_KIND_R8) then
00437 fvelp = .5_SHR_KIND_R8*pi
00438 endif
00439 else if (cossum .lt. 0.0_SHR_KIND_R8) then
00440 fvelp = atan(sinsum/cossum) + pi
00441 else if (cossum .gt. 0.0_SHR_KIND_R8) then
00442 if (sinsum .lt. 0.0_SHR_KIND_R8) then
00443 fvelp = atan(sinsum/cossum) + 2.0_SHR_KIND_R8*pi
00444 else
00445 fvelp = atan(sinsum/cossum)
00446 endif
00447 endif
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 mvsum = 0.0_SHR_KIND_R8
00459 do i = 1, pmvelen
00460 mvsum = mvsum + mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + &
00461 & mvphas(i))*degrad)
00462 end do
00463 mvelp = fvelp/degrad + 50.439273_SHR_KIND_R8*psecdeg*years + 3.392506_SHR_KIND_R8 + mvsum
00464
00465
00466
00467 do while (mvelp .lt. 0.0_SHR_KIND_R8)
00468 mvelp = mvelp + 360.0_SHR_KIND_R8
00469 end do
00470 do while (mvelp .ge. 360.0_SHR_KIND_R8)
00471 mvelp = mvelp - 360.0_SHR_KIND_R8
00472 end do
00473
00474 END IF
00475
00476
00477
00478 obliqr = obliq*degrad
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 mvelpp = (mvelp + 180._SHR_KIND_R8)*degrad
00491
00492
00493
00494 beta = sqrt(1._SHR_KIND_R8 - eccen2)
00495
00496
00497
00498
00499
00500
00501 lambm0 = 2._SHR_KIND_R8*((.5_SHR_KIND_R8*eccen + .125_SHR_KIND_R8*eccen3)*(1._SHR_KIND_R8 + beta)*sin(mvelpp) &
00502 & - .250_SHR_KIND_R8*eccen2*(.5_SHR_KIND_R8 + beta)*sin(2._SHR_KIND_R8*mvelpp) &
00503 & + .125_SHR_KIND_R8*eccen3*(1._SHR_KIND_R8/3._SHR_KIND_R8 + beta)*sin(3._SHR_KIND_R8*mvelpp))
00504
00505 if ( log_print ) then
00506 write(s_logunit,F03) '------ Computed Orbital Parameters ------'
00507 write(s_logunit,F03) 'Eccentricity = ',eccen
00508 write(s_logunit,F03) 'Obliquity (deg) = ',obliq
00509 write(s_logunit,F03) 'Obliquity (rad) = ',obliqr
00510 write(s_logunit,F03) 'Long of perh(deg) = ',mvelp
00511 write(s_logunit,F03) 'Long of perh(rad) = ',mvelpp
00512 write(s_logunit,F03) 'Long at v.e.(rad) = ',lambm0
00513 write(s_logunit,F03) '-----------------------------------------'
00514 end if
00515
00516 END SUBROUTINE shr_orb_params
00517
00518
00519
00520 SUBROUTINE shr_orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr ,delta ,eccf)
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 real (SHR_KIND_R8),intent(in) :: calday
00536 real (SHR_KIND_R8),intent(in) :: eccen
00537 real (SHR_KIND_R8),intent(in) :: obliqr
00538 real (SHR_KIND_R8),intent(in) :: lambm0
00539
00540 real (SHR_KIND_R8),intent(in) :: mvelpp
00541
00542 real (SHR_KIND_R8),intent(out) :: delta
00543 real (SHR_KIND_R8),intent(out) :: eccf
00544
00545
00546 real (SHR_KIND_R8),parameter :: dayspy = 365.0_SHR_KIND_R8
00547 real (SHR_KIND_R8),parameter :: ve = 80.5_SHR_KIND_R8
00548
00549
00550 real (SHR_KIND_R8) :: lambm
00551 real (SHR_KIND_R8) :: lmm
00552 real (SHR_KIND_R8) :: lamb
00553 real (SHR_KIND_R8) :: invrho
00554 real (SHR_KIND_R8) :: sinl
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 lambm = lambm0 + (calday - ve)*2._SHR_KIND_R8*pi/dayspy
00574 lmm = lambm - mvelpp
00575
00576
00577
00578
00579 sinl = sin(lmm)
00580 lamb = lambm + eccen*(2._SHR_KIND_R8*sinl + eccen*(1.25_SHR_KIND_R8*sin(2._SHR_KIND_R8*lmm) &
00581 & + eccen*((13.0_SHR_KIND_R8/12.0_SHR_KIND_R8)*sin(3._SHR_KIND_R8*lmm) - 0.25_SHR_KIND_R8*sinl)))
00582
00583
00584
00585
00586
00587
00588
00589 invrho = (1._SHR_KIND_R8 + eccen*cos(lamb - mvelpp)) / (1._SHR_KIND_R8 - eccen*eccen)
00590
00591
00592
00593 delta = asin(sin(obliqr)*sin(lamb))
00594 eccf = invrho*invrho
00595
00596 return
00597
00598 END SUBROUTINE shr_orb_decl
00599
00600
00601
00602 SUBROUTINE shr_orb_print( iyear_AD, eccen, obliq, mvelp )
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 integer(SHR_KIND_IN),intent(in) :: iyear_AD
00617 real (SHR_KIND_R8),intent(in) :: eccen
00618
00619 real (SHR_KIND_R8),intent(in) :: obliq
00620
00621 real (SHR_KIND_R8),intent(in) :: mvelp
00622
00623
00624 character(len=*),parameter :: F00 = "('(shr_orb_print) ',4a)"
00625 character(len=*),parameter :: F01 = "('(shr_orb_print) ',a,i9.4)"
00626 character(len=*),parameter :: F02 = "('(shr_orb_print) ',a,f6.3)"
00627 character(len=*),parameter :: F03 = "('(shr_orb_print) ',a,es14.6)"
00628
00629
00630 if (s_loglev > 0) then
00631 if ( iyear_AD .ne. SHR_ORB_UNDEF_INT ) then
00632 if ( iyear_AD > 0 ) then
00633 write(s_logunit,F01) 'Orbital parameters calculated for year: AD ',iyear_AD
00634 else
00635 write(s_logunit,F01) 'Orbital parameters calculated for year: BC ',iyear_AD
00636 end if
00637 else if ( obliq /= SHR_ORB_UNDEF_REAL ) then
00638 write(s_logunit,F03) 'Orbital parameters: '
00639 write(s_logunit,F03) 'Obliquity (degree): ', obliq
00640 write(s_logunit,F03) 'Eccentricity (unitless): ', eccen
00641 write(s_logunit,F03) 'Long. of moving Perhelion (deg): ', mvelp
00642 else
00643 write(s_logunit,F03) 'Orbit parameters not set!'
00644 end if
00645 endif
00646
00647 END SUBROUTINE shr_orb_print
00648
00649
00650 END MODULE shr_orb_mod