| 1 | c********************************************************************** |
|---|
| 2 | |
|---|
| 3 | c Includes the following old 1-D model files/subroutines |
|---|
| 4 | |
|---|
| 5 | c -MZTCRSUB_dlvr11.f |
|---|
| 6 | c *dinterconnection |
|---|
| 7 | c *planckd |
|---|
| 8 | c *leetvt |
|---|
| 9 | c -MZTFSUB_dlvr11_02.f |
|---|
| 10 | c *initial |
|---|
| 11 | c *intershphunt |
|---|
| 12 | c *interstrhunt |
|---|
| 13 | c *intzhunt |
|---|
| 14 | c *intzhunt_cts |
|---|
| 15 | c *rhist |
|---|
| 16 | c *we_clean |
|---|
| 17 | c *mztf_correccion |
|---|
| 18 | c *mzescape_normaliz |
|---|
| 19 | c *mzescape_normaliz_02 |
|---|
| 20 | c -interdpESCTVCISO_dlvr11.f |
|---|
| 21 | c -hunt_cts.f |
|---|
| 22 | c -huntdp.f |
|---|
| 23 | c -hunt.f |
|---|
| 24 | c -interdp_limits.f |
|---|
| 25 | c -interhunt2veces.f |
|---|
| 26 | c -interhunt5veces.f |
|---|
| 27 | c -interhuntdp3veces.f |
|---|
| 28 | c -interhuntdp4veces.f |
|---|
| 29 | c -interhuntdp.f |
|---|
| 30 | c -interhunt.f |
|---|
| 31 | c -interhuntlimits2veces.f |
|---|
| 32 | c -interhuntlimits5veces.f |
|---|
| 33 | c -interhuntlimits.f |
|---|
| 34 | c -lubksb_dp.f |
|---|
| 35 | c -ludcmp_dp.f |
|---|
| 36 | c -LUdec.f |
|---|
| 37 | c -mat_oper.f |
|---|
| 38 | c *unit |
|---|
| 39 | c *diago |
|---|
| 40 | c *invdiag |
|---|
| 41 | c *samem |
|---|
| 42 | c *mulmv |
|---|
| 43 | c *trucodiag |
|---|
| 44 | c *trucommvv |
|---|
| 45 | c *sypvmv |
|---|
| 46 | c *mulmm |
|---|
| 47 | c *resmm |
|---|
| 48 | c *sumvv |
|---|
| 49 | c *sypvvv |
|---|
| 50 | c *zerom |
|---|
| 51 | c *zero4m |
|---|
| 52 | c *zero3m |
|---|
| 53 | c *zero2m |
|---|
| 54 | c *zerov |
|---|
| 55 | c *zero4v |
|---|
| 56 | c *zero3v |
|---|
| 57 | c *zero2v |
|---|
| 58 | c -suaviza.f |
|---|
| 59 | |
|---|
| 60 | c********************************************************************** |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | c *** Old MZTCRSUB_dlvr11.f *** |
|---|
| 64 | |
|---|
| 65 | !************************************************************************ |
|---|
| 66 | |
|---|
| 67 | ! subroutine dinterconnection ( v, vt ) |
|---|
| 68 | |
|---|
| 69 | |
|---|
| 70 | ************************************************************************ |
|---|
| 71 | |
|---|
| 72 | ! implicit none |
|---|
| 73 | ! include 'nlte_paramdef.h' |
|---|
| 74 | |
|---|
| 75 | c argumentos |
|---|
| 76 | ! real*8 vt(nl), v(nl) |
|---|
| 77 | |
|---|
| 78 | c local variables |
|---|
| 79 | ! integer i |
|---|
| 80 | |
|---|
| 81 | c ************* |
|---|
| 82 | ! |
|---|
| 83 | ! do i=1,nl |
|---|
| 84 | ! v(i) = vt(i) |
|---|
| 85 | ! end do |
|---|
| 86 | |
|---|
| 87 | ! return |
|---|
| 88 | ! end |
|---|
| 89 | |
|---|
| 90 | c*********************************************************************** |
|---|
| 91 | function planckdp(tp,xnu) |
|---|
| 92 | c*********************************************************************** |
|---|
| 93 | |
|---|
| 94 | implicit none |
|---|
| 95 | |
|---|
| 96 | include 'nlte_paramdef.h' |
|---|
| 97 | |
|---|
| 98 | real*8 planckdp |
|---|
| 99 | real*8 xnu |
|---|
| 100 | real tp |
|---|
| 101 | |
|---|
| 102 | planckdp = gamma*xnu**3.0d0 / exp( ee*xnu/dble(tp) ) |
|---|
| 103 | !erg cm-2.sr-1/cm-1. |
|---|
| 104 | |
|---|
| 105 | c end |
|---|
| 106 | return |
|---|
| 107 | end |
|---|
| 108 | |
|---|
| 109 | c*********************************************************************** |
|---|
| 110 | subroutine leetvt |
|---|
| 111 | |
|---|
| 112 | c*********************************************************************** |
|---|
| 113 | |
|---|
| 114 | implicit none |
|---|
| 115 | |
|---|
| 116 | include 'nlte_paramdef.h' |
|---|
| 117 | include 'nlte_commons.h' |
|---|
| 118 | |
|---|
| 119 | c local variables |
|---|
| 120 | integer i |
|---|
| 121 | real*8 zld(nl), zyd(nzy) |
|---|
| 122 | real*8 xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy) |
|---|
| 123 | |
|---|
| 124 | c*********************************************************************** |
|---|
| 125 | |
|---|
| 126 | do i=1,nzy |
|---|
| 127 | zyd(i) = dble(zy(i)) |
|---|
| 128 | xvt11(i)= dble( ty(i) ) |
|---|
| 129 | xvt21(i)= dble( ty(i) ) |
|---|
| 130 | xvt31(i)= dble( ty(i) ) |
|---|
| 131 | xvt41(i)= dble( ty(i) ) |
|---|
| 132 | end do |
|---|
| 133 | |
|---|
| 134 | do i=1,nl |
|---|
| 135 | zld(i) = dble( zl(i) ) |
|---|
| 136 | enddo |
|---|
| 137 | call interhuntdp4veces ( v626t1,v628t1,v636t1,v627t1, zld,nl, |
|---|
| 138 | $ xvt11, xvt21, xvt31, xvt41, zyd,nzy, 1 ) |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | c end |
|---|
| 142 | return |
|---|
| 143 | end |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | c *** MZTFSUB_dlvr11_02.f *** |
|---|
| 147 | |
|---|
| 148 | |
|---|
| 149 | c **************************************************************** |
|---|
| 150 | subroutine initial |
|---|
| 151 | |
|---|
| 152 | c **************************************************************** |
|---|
| 153 | |
|---|
| 154 | implicit none |
|---|
| 155 | |
|---|
| 156 | include 'nlte_paramdef.h' |
|---|
| 157 | include 'nlte_commons.h' |
|---|
| 158 | |
|---|
| 159 | c local variables |
|---|
| 160 | integer i |
|---|
| 161 | |
|---|
| 162 | c *************** |
|---|
| 163 | |
|---|
| 164 | eqw = 0.0d00 |
|---|
| 165 | aa = 0.0d00 |
|---|
| 166 | cc = 0.0d00 |
|---|
| 167 | dd = 0.0d00 |
|---|
| 168 | |
|---|
| 169 | do i=1,nbox |
|---|
| 170 | ccbox(i) = 0.0d0 |
|---|
| 171 | ddbox(i) = 0.0d0 |
|---|
| 172 | end do |
|---|
| 173 | |
|---|
| 174 | return |
|---|
| 175 | end |
|---|
| 176 | |
|---|
| 177 | c ********************************************************************** |
|---|
| 178 | |
|---|
| 179 | subroutine intershphunt (i, alsx,adx,xtemp) |
|---|
| 180 | |
|---|
| 181 | c ********************************************************************** |
|---|
| 182 | |
|---|
| 183 | implicit none |
|---|
| 184 | |
|---|
| 185 | include 'nlte_paramdef.h' |
|---|
| 186 | include 'nlte_commons.h' |
|---|
| 187 | |
|---|
| 188 | c arguments |
|---|
| 189 | real*8 alsx(nbox_max),adx(nbox_max) ! Output |
|---|
| 190 | real*8 xtemp(nbox_max) ! Input |
|---|
| 191 | integer i ! I , O |
|---|
| 192 | |
|---|
| 193 | c local variables |
|---|
| 194 | integer k |
|---|
| 195 | real*8 factor |
|---|
| 196 | real*8 temperatura ! para evitar valores ligeramnt out of limits |
|---|
| 197 | |
|---|
| 198 | c *********** |
|---|
| 199 | |
|---|
| 200 | do 1, k=1,nbox_max |
|---|
| 201 | temperatura = xtemp(k) |
|---|
| 202 | if (abs(xtemp(k)-thist(1)).le.0.01d0) then |
|---|
| 203 | temperatura=thist(1) |
|---|
| 204 | elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then |
|---|
| 205 | temperatura=thist(nhist) |
|---|
| 206 | endif |
|---|
| 207 | call huntdp ( thist,nhist, temperatura, i ) |
|---|
| 208 | if ( i.eq.0 .or. i.eq.nhist ) then |
|---|
| 209 | write (*,*) ' HUNT/ Limits input grid:', |
|---|
| 210 | @ thist(1),thist(nhist) |
|---|
| 211 | write (*,*) ' HUNT/ location in new grid:', xtemp(k) |
|---|
| 212 | stop ' INTERSHP/ Interpolation error. T out of Histogram.' |
|---|
| 213 | endif |
|---|
| 214 | factor = 1.d0 / (thist(i+1)-thist(i)) |
|---|
| 215 | alsx(k) = (( xls1(i,k)*(thist(i+1)-xtemp(k)) + |
|---|
| 216 | @ xls1(i+1,k)*(xtemp(k)-thist(i)) )) * factor |
|---|
| 217 | adx(k) = (( xld1(i,k)*(thist(i+1)-xtemp(k)) + |
|---|
| 218 | @ xld1(i+1,k)*(xtemp(k)-thist(i)) )) * factor |
|---|
| 219 | 1 continue |
|---|
| 220 | |
|---|
| 221 | return |
|---|
| 222 | end |
|---|
| 223 | |
|---|
| 224 | c ********************************************************************** |
|---|
| 225 | |
|---|
| 226 | subroutine interstrhunt (i, stx, ts, sx, xtemp ) |
|---|
| 227 | |
|---|
| 228 | c ********************************************************************** |
|---|
| 229 | |
|---|
| 230 | implicit none |
|---|
| 231 | |
|---|
| 232 | include 'nlte_paramdef.h' |
|---|
| 233 | include 'nlte_commons.h' |
|---|
| 234 | |
|---|
| 235 | c arguments |
|---|
| 236 | real*8 stx ! output, total band strength |
|---|
| 237 | real*8 ts ! input, temp for stx |
|---|
| 238 | real*8 sx(nbox_max) ! output, strength for each box |
|---|
| 239 | real*8 xtemp(nbox_max) ! input, temp for sx |
|---|
| 240 | integer i |
|---|
| 241 | |
|---|
| 242 | c local variables |
|---|
| 243 | integer k |
|---|
| 244 | real*8 factor |
|---|
| 245 | real*8 temperatura |
|---|
| 246 | |
|---|
| 247 | c *********** |
|---|
| 248 | |
|---|
| 249 | do 1, k=1,nbox |
|---|
| 250 | temperatura = xtemp(k) |
|---|
| 251 | if (abs(xtemp(k)-thist(1)).le.0.01d0) then |
|---|
| 252 | temperatura=thist(1) |
|---|
| 253 | elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then |
|---|
| 254 | temperatura=thist(nhist) |
|---|
| 255 | endif |
|---|
| 256 | call huntdp ( thist,nhist, temperatura, i ) |
|---|
| 257 | if ( i.eq.0 .or. i.eq.nhist ) then |
|---|
| 258 | write(*,*)'HUNT/ Limits input grid:',thist(1),thist(nhist) |
|---|
| 259 | write(*,*)'HUNT/ location in new grid:',xtemp(k) |
|---|
| 260 | stop'INTERSTR/1/ Interpolation error. T out of Histogram.' |
|---|
| 261 | endif |
|---|
| 262 | factor = 1.d0 / (thist(i+1)-thist(i)) |
|---|
| 263 | sx(k) = ( sk1(i,k) * (thist(i+1)-xtemp(k)) |
|---|
| 264 | @ + sk1(i+1,k) * (xtemp(k)-thist(i)) ) * factor |
|---|
| 265 | 1 continue |
|---|
| 266 | |
|---|
| 267 | |
|---|
| 268 | temperatura = ts |
|---|
| 269 | if (abs(ts-thist(1)).le.0.01d0) then |
|---|
| 270 | temperatura=thist(1) |
|---|
| 271 | elseif (abs(ts-thist(nhist)).le.0.01d0) then |
|---|
| 272 | temperatura=thist(nhist) |
|---|
| 273 | endif |
|---|
| 274 | call huntdp ( thist,nhist, temperatura, i ) |
|---|
| 275 | if ( i.eq.0 .or. i.eq.nhist ) then |
|---|
| 276 | write (*,*) ' HUNT/ Limits input grid:', |
|---|
| 277 | @ thist(1),thist(nhist) |
|---|
| 278 | write (*,*) ' HUNT/ location in new grid:', ts |
|---|
| 279 | stop ' INTERSTR/2/ Interpolat error. T out of Histogram.' |
|---|
| 280 | endif |
|---|
| 281 | factor = 1.d0 / (thist(i+1)-thist(i)) |
|---|
| 282 | stx = 0.d0 |
|---|
| 283 | do k=1,nbox |
|---|
| 284 | stx = stx + no(k) * ( sk1(i,k)*(thist(i+1)-ts) + |
|---|
| 285 | @ sk1(i+1,k)*(ts-thist(i)) ) * factor |
|---|
| 286 | end do |
|---|
| 287 | |
|---|
| 288 | |
|---|
| 289 | return |
|---|
| 290 | end |
|---|
| 291 | |
|---|
| 292 | c ********************************************************************** |
|---|
| 293 | |
|---|
| 294 | subroutine intzhunt (k, h, aco2,ap,amr,at, con) |
|---|
| 295 | |
|---|
| 296 | c k lleva la posicion de la ultima llamada a intz , necesario para |
|---|
| 297 | c que esto represente una aceleracion real. |
|---|
| 298 | c ********************************************************************** |
|---|
| 299 | |
|---|
| 300 | implicit none |
|---|
| 301 | include 'nlte_paramdef.h' |
|---|
| 302 | include 'nlte_commons.h' |
|---|
| 303 | |
|---|
| 304 | c arguments |
|---|
| 305 | real h ! i |
|---|
| 306 | real*8 con(nzy) ! i |
|---|
| 307 | real*8 aco2, ap, at, amr ! o |
|---|
| 308 | integer k ! i |
|---|
| 309 | |
|---|
| 310 | c local variables |
|---|
| 311 | real factor |
|---|
| 312 | |
|---|
| 313 | c ************ |
|---|
| 314 | |
|---|
| 315 | call hunt ( zy,nzy, h, k ) |
|---|
| 316 | factor = (h-zy(k)) / (zy(k+1)-zy(k)) |
|---|
| 317 | ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * factor ) ) |
|---|
| 318 | aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor) |
|---|
| 319 | aco2 = exp( aco2 ) |
|---|
| 320 | at = dble( ty(k) + (ty(k+1)-ty(k)) * factor ) |
|---|
| 321 | amr = dble( mr(k) + (mr(k+1)-mr(k)) * factor ) |
|---|
| 322 | |
|---|
| 323 | |
|---|
| 324 | return |
|---|
| 325 | end |
|---|
| 326 | |
|---|
| 327 | c ********************************************************************** |
|---|
| 328 | |
|---|
| 329 | subroutine intzhunt_cts (k, h, nzy_cts_real, |
|---|
| 330 | @ aco2,ap,amr,at, con) |
|---|
| 331 | |
|---|
| 332 | c k lleva la posicion de la ultima llamada a intz , necesario para |
|---|
| 333 | c que esto represente una aceleracion real. |
|---|
| 334 | c ********************************************************************** |
|---|
| 335 | |
|---|
| 336 | implicit none |
|---|
| 337 | include 'nlte_paramdef.h' |
|---|
| 338 | include 'nlte_commons.h' |
|---|
| 339 | |
|---|
| 340 | c arguments |
|---|
| 341 | real h ! i |
|---|
| 342 | real*8 con(nzy_cts) ! i |
|---|
| 343 | real*8 aco2, ap, at, amr ! o |
|---|
| 344 | integer k ! i |
|---|
| 345 | integer nzy_cts_real ! i |
|---|
| 346 | |
|---|
| 347 | c local variables |
|---|
| 348 | real factor |
|---|
| 349 | |
|---|
| 350 | c ************ |
|---|
| 351 | |
|---|
| 352 | call hunt_cts ( zy_cts,nzy_cts, nzy_cts_real, h, k ) |
|---|
| 353 | factor = (h-zy_cts(k)) / (zy_cts(k+1)-zy_cts(k)) |
|---|
| 354 | ap = dble( exp( log(py_cts(k)) + |
|---|
| 355 | @ log(py_cts(k+1)/py_cts(k)) * factor ) ) |
|---|
| 356 | aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor) |
|---|
| 357 | aco2 = exp( aco2 ) |
|---|
| 358 | at = dble( ty_cts(k) + (ty_cts(k+1)-ty_cts(k)) * factor ) |
|---|
| 359 | amr = dble( mr_cts(k) + (mr_cts(k+1)-mr_cts(k)) * factor ) |
|---|
| 360 | |
|---|
| 361 | |
|---|
| 362 | return |
|---|
| 363 | end |
|---|
| 364 | |
|---|
| 365 | |
|---|
| 366 | c ********************************************************************** |
|---|
| 367 | |
|---|
| 368 | real*8 function we_clean ( y,pl, xalsa, xalda ) |
|---|
| 369 | |
|---|
| 370 | c ********************************************************************** |
|---|
| 371 | |
|---|
| 372 | implicit none |
|---|
| 373 | |
|---|
| 374 | include 'nlte_paramdef.h' |
|---|
| 375 | |
|---|
| 376 | c arguments |
|---|
| 377 | real*8 y ! I. path's absorber amount * strength |
|---|
| 378 | real*8 pl ! I. path's partial pressure of CO2 |
|---|
| 379 | real*8 xalsa ! I. Self lorentz linewidth for 1 isot & 1 box |
|---|
| 380 | real*8 xalda ! I. Doppler linewidth " " |
|---|
| 381 | |
|---|
| 382 | c local variables |
|---|
| 383 | integer i |
|---|
| 384 | real*8 x,wl,wd,wvoigt |
|---|
| 385 | real*8 cn(0:7),dn(0:7) |
|---|
| 386 | real*8 factor, denom |
|---|
| 387 | real*8 pi, pi2, sqrtpi |
|---|
| 388 | |
|---|
| 389 | c data blocks |
|---|
| 390 | data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2, |
|---|
| 391 | @ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4, |
|---|
| 392 | @ 3.42209457833d-5,-1.28380804108d-6/ |
|---|
| 393 | data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1, |
|---|
| 394 | @ 8.21896973657d-1,-2.5222672453,6.1007027481, |
|---|
| 395 | @ -8.51001627836,4.6535116765/ |
|---|
| 396 | |
|---|
| 397 | c *********** |
|---|
| 398 | |
|---|
| 399 | pi = 3.141592 |
|---|
| 400 | pi2= 6.28318531 |
|---|
| 401 | sqrtpi = 1.77245385 |
|---|
| 402 | |
|---|
| 403 | x=y / ( pi2 * xalsa*pl ) |
|---|
| 404 | |
|---|
| 405 | |
|---|
| 406 | c Lorentz |
|---|
| 407 | wl=y/sqrt(1.0d0+pi*x/2.0d0) |
|---|
| 408 | |
|---|
| 409 | c Doppler |
|---|
| 410 | x = y / (xalda*sqrtpi) |
|---|
| 411 | if (x .lt. 5.0d0) then |
|---|
| 412 | wd = cn(0) |
|---|
| 413 | factor = 1.d0 |
|---|
| 414 | do i=1,7 |
|---|
| 415 | factor = factor * x |
|---|
| 416 | wd = wd + cn(i) * factor |
|---|
| 417 | end do |
|---|
| 418 | wd = xalda * x * sqrtpi * wd |
|---|
| 419 | else |
|---|
| 420 | wd = dn(0) |
|---|
| 421 | factor = 1.d0 / log(x) |
|---|
| 422 | denom = 1.d0 |
|---|
| 423 | do i=1,7 |
|---|
| 424 | denom = denom * factor |
|---|
| 425 | wd = wd + dn(i) * denom |
|---|
| 426 | end do |
|---|
| 427 | wd = xalda * sqrt(log(x)) * wd |
|---|
| 428 | end if |
|---|
| 429 | |
|---|
| 430 | c Voigt |
|---|
| 431 | wvoigt = wl*wl + wd*wd - (wd*wl/y)*(wd*wl/y) |
|---|
| 432 | |
|---|
| 433 | if ( wvoigt .lt. 0.0d0 ) then |
|---|
| 434 | write (*,*) ' Subroutine WE/ Error in Voift EQS calculation' |
|---|
| 435 | write (*,*) ' WL, WD, X, Y = ', wl, wd, x, y |
|---|
| 436 | stop ' ERROR : Imaginary EQW. Revise spectral data. ' |
|---|
| 437 | endif |
|---|
| 438 | |
|---|
| 439 | we_clean = sqrt( wvoigt ) |
|---|
| 440 | |
|---|
| 441 | |
|---|
| 442 | return |
|---|
| 443 | end |
|---|
| 444 | |
|---|
| 445 | |
|---|
| 446 | c *********************************************************************** |
|---|
| 447 | |
|---|
| 448 | subroutine mztf_correccion (coninf, con, ib ) |
|---|
| 449 | |
|---|
| 450 | c *********************************************************************** |
|---|
| 451 | |
|---|
| 452 | implicit none |
|---|
| 453 | |
|---|
| 454 | include 'nlte_paramdef.h' |
|---|
| 455 | include 'nlte_commons.h' |
|---|
| 456 | |
|---|
| 457 | c arguments |
|---|
| 458 | integer ib |
|---|
| 459 | real*8 con(nzy), coninf |
|---|
| 460 | |
|---|
| 461 | ! local variables |
|---|
| 462 | integer i, isot |
|---|
| 463 | real*8 tvt0(nzy), tvtbs(nzy), zld(nl),zyd(nzy) |
|---|
| 464 | real*8 xqv, xes, xlower, xfactor |
|---|
| 465 | |
|---|
| 466 | c ********* |
|---|
| 467 | |
|---|
| 468 | isot = 1 |
|---|
| 469 | nu11 = dble( nu(1,1) ) |
|---|
| 470 | |
|---|
| 471 | do i=1,nzy |
|---|
| 472 | zyd(i) = dble(zy(i)) |
|---|
| 473 | enddo |
|---|
| 474 | do i=1,nl |
|---|
| 475 | zld(i) = dble( zl(i) ) |
|---|
| 476 | end do |
|---|
| 477 | |
|---|
| 478 | ! tvtbs |
|---|
| 479 | call interhuntdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) |
|---|
| 480 | |
|---|
| 481 | ! tvt0 |
|---|
| 482 | if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4) then |
|---|
| 483 | call interhuntdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 ) |
|---|
| 484 | else |
|---|
| 485 | do i=1,nzy |
|---|
| 486 | tvt0(i) = dble( ty(i) ) |
|---|
| 487 | end do |
|---|
| 488 | end if |
|---|
| 489 | |
|---|
| 490 | c factor |
|---|
| 491 | do i=1,nzy |
|---|
| 492 | |
|---|
| 493 | xlower = exp( ee*dble(elow(isot,ib)) * |
|---|
| 494 | @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) |
|---|
| 495 | xes = 1.0d0 |
|---|
| 496 | xqv = ( 1.d0-exp( -ee*nu11/tvtbs(i) ) ) / |
|---|
| 497 | @ (1.d0-exp( -ee*nu11/dble(ty(i)) )) |
|---|
| 498 | xfactor = xlower * xqv**2.d0 * xes |
|---|
| 499 | |
|---|
| 500 | con(i) = con(i) * xfactor |
|---|
| 501 | if (i.eq.nzy) coninf = coninf * xfactor |
|---|
| 502 | |
|---|
| 503 | end do |
|---|
| 504 | |
|---|
| 505 | |
|---|
| 506 | return |
|---|
| 507 | end |
|---|
| 508 | |
|---|
| 509 | |
|---|
| 510 | c *********************************************************************** |
|---|
| 511 | |
|---|
| 512 | subroutine mzescape_normaliz ( taustar, istyle ) |
|---|
| 513 | |
|---|
| 514 | c *********************************************************************** |
|---|
| 515 | |
|---|
| 516 | implicit none |
|---|
| 517 | include 'nlte_paramdef.h' |
|---|
| 518 | |
|---|
| 519 | c arguments |
|---|
| 520 | real*8 taustar(nl) ! o |
|---|
| 521 | integer istyle ! i |
|---|
| 522 | |
|---|
| 523 | c local variables and constants |
|---|
| 524 | integer i, imaximum |
|---|
| 525 | real*8 maximum |
|---|
| 526 | |
|---|
| 527 | c *************** |
|---|
| 528 | |
|---|
| 529 | taustar(nl) = taustar(nl-1) |
|---|
| 530 | |
|---|
| 531 | if ( istyle .eq. 1 ) then |
|---|
| 532 | imaximum = nl |
|---|
| 533 | maximum = taustar(nl) |
|---|
| 534 | do i=1,nl-1 |
|---|
| 535 | if (taustar(i).gt.maximum) taustar(i) = taustar(nl) |
|---|
| 536 | enddo |
|---|
| 537 | elseif ( istyle .eq. 2 ) then |
|---|
| 538 | imaximum = nl |
|---|
| 539 | maximum = taustar(nl) |
|---|
| 540 | do i=nl-1,1,-1 |
|---|
| 541 | if (taustar(i).gt.maximum) then |
|---|
| 542 | maximum = taustar(i) |
|---|
| 543 | imaximum = i |
|---|
| 544 | endif |
|---|
| 545 | enddo |
|---|
| 546 | do i=imaximum,nl |
|---|
| 547 | if (taustar(i).lt.maximum) taustar(i) = maximum |
|---|
| 548 | enddo |
|---|
| 549 | endif |
|---|
| 550 | |
|---|
| 551 | do i=1,nl |
|---|
| 552 | taustar(i) = taustar(i) / maximum |
|---|
| 553 | enddo |
|---|
| 554 | |
|---|
| 555 | |
|---|
| 556 | c end |
|---|
| 557 | return |
|---|
| 558 | end |
|---|
| 559 | |
|---|
| 560 | c *********************************************************************** |
|---|
| 561 | |
|---|
| 562 | subroutine mzescape_normaliz_02 ( taustar, nn, istyle ) |
|---|
| 563 | |
|---|
| 564 | c *********************************************************************** |
|---|
| 565 | |
|---|
| 566 | implicit none |
|---|
| 567 | |
|---|
| 568 | c arguments |
|---|
| 569 | real*8 taustar(nn) ! i,o |
|---|
| 570 | integer istyle ! i |
|---|
| 571 | integer nn ! i |
|---|
| 572 | |
|---|
| 573 | c local variables and constants |
|---|
| 574 | integer i, imaximum |
|---|
| 575 | real*8 maximum |
|---|
| 576 | |
|---|
| 577 | c *************** |
|---|
| 578 | |
|---|
| 579 | taustar(nn) = taustar(nn-1) |
|---|
| 580 | |
|---|
| 581 | if ( istyle .eq. 1 ) then |
|---|
| 582 | imaximum = nn |
|---|
| 583 | maximum = taustar(nn) |
|---|
| 584 | do i=1,nn-1 |
|---|
| 585 | if (taustar(i).gt.maximum) taustar(i) = taustar(nn) |
|---|
| 586 | enddo |
|---|
| 587 | elseif ( istyle .eq. 2 ) then |
|---|
| 588 | imaximum = nn |
|---|
| 589 | maximum = taustar(nn) |
|---|
| 590 | do i=nn-1,1,-1 |
|---|
| 591 | if (taustar(i).gt.maximum) then |
|---|
| 592 | maximum = taustar(i) |
|---|
| 593 | imaximum = i |
|---|
| 594 | endif |
|---|
| 595 | enddo |
|---|
| 596 | do i=imaximum,nn |
|---|
| 597 | if (taustar(i).lt.maximum) taustar(i) = maximum |
|---|
| 598 | enddo |
|---|
| 599 | endif |
|---|
| 600 | |
|---|
| 601 | do i=1,nn |
|---|
| 602 | taustar(i) = taustar(i) / maximum |
|---|
| 603 | enddo |
|---|
| 604 | |
|---|
| 605 | |
|---|
| 606 | c end |
|---|
| 607 | return |
|---|
| 608 | end |
|---|
| 609 | |
|---|
| 610 | |
|---|
| 611 | c *** interdp_ESCTVCISO_dlvr11.f *** |
|---|
| 612 | |
|---|
| 613 | c*********************************************************************** |
|---|
| 614 | |
|---|
| 615 | subroutine interdp_ESCTVCISO |
|---|
| 616 | |
|---|
| 617 | c*********************************************************************** |
|---|
| 618 | |
|---|
| 619 | implicit none |
|---|
| 620 | |
|---|
| 621 | include 'nlte_paramdef.h' |
|---|
| 622 | include 'nlte_commons.h' |
|---|
| 623 | |
|---|
| 624 | c local variables |
|---|
| 625 | integer i |
|---|
| 626 | real*8 lnpnb(nl) |
|---|
| 627 | |
|---|
| 628 | |
|---|
| 629 | c*********************************************************************** |
|---|
| 630 | |
|---|
| 631 | c Use pressure in the NLTE grid but in log and in nb |
|---|
| 632 | do i=1,nl |
|---|
| 633 | lnpnb(i) = log( dble( pl(i) * 1013.25 * 1.e6) ) |
|---|
| 634 | enddo |
|---|
| 635 | |
|---|
| 636 | c Interpolations |
|---|
| 637 | |
|---|
| 638 | call interhuntdp3veces |
|---|
| 639 | @ ( taustar21,taustar31,taustar41, lnpnb, nl, |
|---|
| 640 | @ tstar21tab,tstar31tab,tstar41tab, lnpnbtab, nztabul, |
|---|
| 641 | @ 1 ) |
|---|
| 642 | |
|---|
| 643 | call interhuntdp3veces ( vc210,vc310,vc410, lnpnb, nl, |
|---|
| 644 | @ vc210tab,vc310tab,vc410tab, lnpnbtab, nztabul, 2 ) |
|---|
| 645 | |
|---|
| 646 | c end |
|---|
| 647 | return |
|---|
| 648 | end |
|---|
| 649 | |
|---|
| 650 | |
|---|
| 651 | c *** hunt_cts.f *** |
|---|
| 652 | |
|---|
| 653 | cccc |
|---|
| 654 | SUBROUTINE hunt_cts(xx,n,n_cts,x,jlo) |
|---|
| 655 | c |
|---|
| 656 | c La dif con hunt es el uso de un indice superior (n_cts) mas bajito que (n) |
|---|
| 657 | c |
|---|
| 658 | c Arguments |
|---|
| 659 | INTEGER jlo ! O |
|---|
| 660 | INTEGER n ! I |
|---|
| 661 | INTEGER n_cts ! I |
|---|
| 662 | REAL xx(n) ! I |
|---|
| 663 | REAL x ! I |
|---|
| 664 | |
|---|
| 665 | c Local variables |
|---|
| 666 | INTEGER inc,jhi,jm |
|---|
| 667 | LOGICAL ascnd |
|---|
| 668 | c |
|---|
| 669 | cccc |
|---|
| 670 | c |
|---|
| 671 | ascnd=xx(n_cts).ge.xx(1) |
|---|
| 672 | if(jlo.le.0.or.jlo.gt.n_cts)then |
|---|
| 673 | jlo=0 |
|---|
| 674 | jhi=n_cts+1 |
|---|
| 675 | goto 3 |
|---|
| 676 | endif |
|---|
| 677 | inc=1 |
|---|
| 678 | if(x.ge.xx(jlo).eqv.ascnd)then |
|---|
| 679 | 1 jhi=jlo+inc |
|---|
| 680 | ! write (*,*) jlo |
|---|
| 681 | if(jhi.gt.n_cts)then |
|---|
| 682 | jhi=n_cts+1 |
|---|
| 683 | ! write (*,*) jhi-1 |
|---|
| 684 | else if(x.ge.xx(jhi).eqv.ascnd)then |
|---|
| 685 | jlo=jhi |
|---|
| 686 | inc=inc+inc |
|---|
| 687 | ! write (*,*) jlo |
|---|
| 688 | goto 1 |
|---|
| 689 | endif |
|---|
| 690 | else |
|---|
| 691 | jhi=jlo |
|---|
| 692 | 2 jlo=jhi-inc |
|---|
| 693 | ! write (*,*) jlo |
|---|
| 694 | if(jlo.lt.1)then |
|---|
| 695 | jlo=0 |
|---|
| 696 | else if(x.lt.xx(jlo).eqv.ascnd)then |
|---|
| 697 | jhi=jlo |
|---|
| 698 | inc=inc+inc |
|---|
| 699 | goto 2 |
|---|
| 700 | endif |
|---|
| 701 | endif |
|---|
| 702 | 3 if(jhi-jlo.eq.1)then |
|---|
| 703 | if(x.eq.xx(n_cts))jlo=n_cts-1 |
|---|
| 704 | if(x.eq.xx(1))jlo=1 |
|---|
| 705 | ! write (*,*) jlo |
|---|
| 706 | return |
|---|
| 707 | endif |
|---|
| 708 | jm=(jhi+jlo)/2 |
|---|
| 709 | if(x.ge.xx(jm).eqv.ascnd)then |
|---|
| 710 | jlo=jm |
|---|
| 711 | else |
|---|
| 712 | jhi=jm |
|---|
| 713 | endif |
|---|
| 714 | ! write (*,*) jhi-1 |
|---|
| 715 | goto 3 |
|---|
| 716 | c |
|---|
| 717 | END |
|---|
| 718 | |
|---|
| 719 | |
|---|
| 720 | c *** huntdp.f *** |
|---|
| 721 | |
|---|
| 722 | cccc |
|---|
| 723 | SUBROUTINE huntdp(xx,n,x,jlo) |
|---|
| 724 | c |
|---|
| 725 | c Arguments |
|---|
| 726 | INTEGER jlo ! O |
|---|
| 727 | INTEGER n ! I |
|---|
| 728 | REAL*8 xx(n) ! I |
|---|
| 729 | REAL*8 x ! I |
|---|
| 730 | |
|---|
| 731 | c Local variables |
|---|
| 732 | INTEGER inc,jhi,jm |
|---|
| 733 | LOGICAL ascnd |
|---|
| 734 | c |
|---|
| 735 | cccc |
|---|
| 736 | c |
|---|
| 737 | ascnd=xx(n).ge.xx(1) |
|---|
| 738 | if(jlo.le.0.or.jlo.gt.n)then |
|---|
| 739 | jlo=0 |
|---|
| 740 | jhi=n+1 |
|---|
| 741 | goto 3 |
|---|
| 742 | endif |
|---|
| 743 | inc=1 |
|---|
| 744 | if(x.ge.xx(jlo).eqv.ascnd)then |
|---|
| 745 | 1 jhi=jlo+inc |
|---|
| 746 | if(jhi.gt.n)then |
|---|
| 747 | jhi=n+1 |
|---|
| 748 | else if(x.ge.xx(jhi).eqv.ascnd)then |
|---|
| 749 | jlo=jhi |
|---|
| 750 | inc=inc+inc |
|---|
| 751 | goto 1 |
|---|
| 752 | endif |
|---|
| 753 | else |
|---|
| 754 | jhi=jlo |
|---|
| 755 | 2 jlo=jhi-inc |
|---|
| 756 | if(jlo.lt.1)then |
|---|
| 757 | jlo=0 |
|---|
| 758 | else if(x.lt.xx(jlo).eqv.ascnd)then |
|---|
| 759 | jhi=jlo |
|---|
| 760 | inc=inc+inc |
|---|
| 761 | goto 2 |
|---|
| 762 | endif |
|---|
| 763 | endif |
|---|
| 764 | 3 if(jhi-jlo.eq.1)then |
|---|
| 765 | if(x.eq.xx(n))jlo=n-1 |
|---|
| 766 | if(x.eq.xx(1))jlo=1 |
|---|
| 767 | return |
|---|
| 768 | endif |
|---|
| 769 | jm=(jhi+jlo)/2 |
|---|
| 770 | if(x.ge.xx(jm).eqv.ascnd)then |
|---|
| 771 | jlo=jm |
|---|
| 772 | else |
|---|
| 773 | jhi=jm |
|---|
| 774 | endif |
|---|
| 775 | goto 3 |
|---|
| 776 | c |
|---|
| 777 | END |
|---|
| 778 | |
|---|
| 779 | |
|---|
| 780 | c *** hunt.f *** |
|---|
| 781 | |
|---|
| 782 | cccc |
|---|
| 783 | SUBROUTINE hunt(xx,n,x,jlo) |
|---|
| 784 | c |
|---|
| 785 | c Arguments |
|---|
| 786 | INTEGER jlo ! O |
|---|
| 787 | INTEGER n ! I |
|---|
| 788 | REAL xx(n) ! I |
|---|
| 789 | REAL x ! I |
|---|
| 790 | |
|---|
| 791 | c Local variables |
|---|
| 792 | INTEGER inc,jhi,jm |
|---|
| 793 | LOGICAL ascnd |
|---|
| 794 | c |
|---|
| 795 | cccc |
|---|
| 796 | c |
|---|
| 797 | ascnd=xx(n).ge.xx(1) |
|---|
| 798 | if(jlo.le.0.or.jlo.gt.n)then |
|---|
| 799 | jlo=0 |
|---|
| 800 | jhi=n+1 |
|---|
| 801 | goto 3 |
|---|
| 802 | endif |
|---|
| 803 | inc=1 |
|---|
| 804 | if(x.ge.xx(jlo).eqv.ascnd)then |
|---|
| 805 | 1 jhi=jlo+inc |
|---|
| 806 | ! write (*,*) jlo |
|---|
| 807 | if(jhi.gt.n)then |
|---|
| 808 | jhi=n+1 |
|---|
| 809 | ! write (*,*) jhi-1 |
|---|
| 810 | else if(x.ge.xx(jhi).eqv.ascnd)then |
|---|
| 811 | jlo=jhi |
|---|
| 812 | inc=inc+inc |
|---|
| 813 | ! write (*,*) jlo |
|---|
| 814 | goto 1 |
|---|
| 815 | endif |
|---|
| 816 | else |
|---|
| 817 | jhi=jlo |
|---|
| 818 | 2 jlo=jhi-inc |
|---|
| 819 | ! write (*,*) jlo |
|---|
| 820 | if(jlo.lt.1)then |
|---|
| 821 | jlo=0 |
|---|
| 822 | else if(x.lt.xx(jlo).eqv.ascnd)then |
|---|
| 823 | jhi=jlo |
|---|
| 824 | inc=inc+inc |
|---|
| 825 | goto 2 |
|---|
| 826 | endif |
|---|
| 827 | endif |
|---|
| 828 | 3 if(jhi-jlo.eq.1)then |
|---|
| 829 | if(x.eq.xx(n))jlo=n-1 |
|---|
| 830 | if(x.eq.xx(1))jlo=1 |
|---|
| 831 | ! write (*,*) jlo |
|---|
| 832 | return |
|---|
| 833 | endif |
|---|
| 834 | jm=(jhi+jlo)/2 |
|---|
| 835 | if(x.ge.xx(jm).eqv.ascnd)then |
|---|
| 836 | jlo=jm |
|---|
| 837 | else |
|---|
| 838 | jhi=jm |
|---|
| 839 | endif |
|---|
| 840 | ! write (*,*) jhi-1 |
|---|
| 841 | goto 3 |
|---|
| 842 | c |
|---|
| 843 | END |
|---|
| 844 | |
|---|
| 845 | |
|---|
| 846 | c *** interdp_limits.f *** |
|---|
| 847 | |
|---|
| 848 | c *********************************************************************** |
|---|
| 849 | |
|---|
| 850 | subroutine interdp_limits ( yy, zz, m, i1,i2, |
|---|
| 851 | @ y, z, n, j1,j2, opt) |
|---|
| 852 | |
|---|
| 853 | c Interpolation soubroutine. |
|---|
| 854 | c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m |
|---|
| 855 | c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n |
|---|
| 856 | c Input values: y(n) , z(n) (solo se usarann los valores entre j1,j2) |
|---|
| 857 | c zz(m) (solo se necesita entre i1,i2) |
|---|
| 858 | c Output values: yy(m) (solo se calculan entre i1,i2) |
|---|
| 859 | c Options: opt=1 -> lineal ,, opt=2 -> logarithmic |
|---|
| 860 | c Difference with interdp: |
|---|
| 861 | c here interpolation proceeds between indexes i1,i2 only |
|---|
| 862 | c if i1=1 & i2=m, both subroutines are exactly the same |
|---|
| 863 | c thus previous calls to interdp or interdp2 could be easily replaced |
|---|
| 864 | |
|---|
| 865 | c JAN 98 MALV Version for mz1d |
|---|
| 866 | c *********************************************************************** |
|---|
| 867 | |
|---|
| 868 | implicit none |
|---|
| 869 | |
|---|
| 870 | ! Arguments |
|---|
| 871 | integer n,m ! I. Dimensions |
|---|
| 872 | integer i1, i2, j1, j2, opt ! I |
|---|
| 873 | real*8 zz(m) ! I |
|---|
| 874 | real*8 yy(m) ! O |
|---|
| 875 | real*8 z(n),y(n) ! I |
|---|
| 876 | |
|---|
| 877 | ! Local variables |
|---|
| 878 | integer i,j |
|---|
| 879 | real*8 zmin,zzmin,zmax,zzmax |
|---|
| 880 | |
|---|
| 881 | c ******************************* |
|---|
| 882 | |
|---|
| 883 | ! write (*,*) ' d interpolating ' |
|---|
| 884 | ! call mindp_limits (z,n,zmin, j1,j2) |
|---|
| 885 | ! call mindp_limits (zz,m,zzmin, i1,i2) |
|---|
| 886 | ! call maxdp_limits (z,n,zmax, j1,j2) |
|---|
| 887 | ! call maxdp_limits (zz,m,zzmax, i1,i2) |
|---|
| 888 | zmin=minval(z(j1:j2)) |
|---|
| 889 | zzmin=minval(zz(i1:i2)) |
|---|
| 890 | zmax=maxval(z(j1:j2)) |
|---|
| 891 | zzmax=maxval(zz(i1:i2)) |
|---|
| 892 | |
|---|
| 893 | if(zzmin.lt.zmin)then |
|---|
| 894 | write (*,*) 'from d interp: new variable out of limits' |
|---|
| 895 | write (*,*) zzmin,'must be .ge. ',zmin |
|---|
| 896 | stop |
|---|
| 897 | ! elseif(zzmax.gt.zmax)then |
|---|
| 898 | ! type *,'from interp: new variable out of limits' |
|---|
| 899 | ! type *,zzmax, 'must be .le. ',zmax |
|---|
| 900 | ! stop |
|---|
| 901 | end if |
|---|
| 902 | |
|---|
| 903 | do 1,i=i1,i2 |
|---|
| 904 | |
|---|
| 905 | do 2,j=j1,j2-1 |
|---|
| 906 | if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 |
|---|
| 907 | 2 continue |
|---|
| 908 | c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2 |
|---|
| 909 | if(opt.eq.1)then |
|---|
| 910 | yy(i)=y(j2-1)+(y(j2)-y(j2-1))*(zz(i)-z(j2-1))/ |
|---|
| 911 | $ (z(j2)-z(j2-1)) |
|---|
| 912 | elseif(opt.eq.2)then |
|---|
| 913 | if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then |
|---|
| 914 | yy(i)=0.0d0 |
|---|
| 915 | else |
|---|
| 916 | yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))* |
|---|
| 917 | @ (zz(i)-z(j2-1))/(z(j2)-z(j2-1))) |
|---|
| 918 | end if |
|---|
| 919 | else |
|---|
| 920 | write (*,*) ' d interp : opt must be 1 or 2, opt= ',opt |
|---|
| 921 | end if |
|---|
| 922 | goto 1 |
|---|
| 923 | 3 continue |
|---|
| 924 | if(opt.eq.1)then |
|---|
| 925 | yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 926 | ! type *, ' ' |
|---|
| 927 | ! type *, ' z(j),z(j+1) =', z(j),z(j+1) |
|---|
| 928 | ! type *, ' t(j),t(j+1) =', y(j),y(j+1) |
|---|
| 929 | ! type *, ' zz, tt = ', zz(i), yy(i) |
|---|
| 930 | elseif(opt.eq.2)then |
|---|
| 931 | if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then |
|---|
| 932 | yy(i)=0.0d0 |
|---|
| 933 | else |
|---|
| 934 | yy(i)=exp(log(y(j))+log(y(j+1)/y(j))* |
|---|
| 935 | @ (zz(i)-z(j))/(z(j+1)-z(j))) |
|---|
| 936 | end if |
|---|
| 937 | else |
|---|
| 938 | write (*,*) ' interp : opt must be 1 or 2, opt= ',opt |
|---|
| 939 | end if |
|---|
| 940 | 1 continue |
|---|
| 941 | return |
|---|
| 942 | end |
|---|
| 943 | |
|---|
| 944 | |
|---|
| 945 | |
|---|
| 946 | c *** interhunt2veces.f *** |
|---|
| 947 | |
|---|
| 948 | c *********************************************************************** |
|---|
| 949 | |
|---|
| 950 | subroutine interhunt2veces ( y1,y2, zz,m, |
|---|
| 951 | @ x1,x2, z,n, opt) |
|---|
| 952 | |
|---|
| 953 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 954 | c input values: y(n) at z(n) |
|---|
| 955 | c output values: yy(m) at zz(m) |
|---|
| 956 | c options: 1 -> lineal |
|---|
| 957 | c 2 -> logarithmic |
|---|
| 958 | c *********************************************************************** |
|---|
| 959 | |
|---|
| 960 | implicit none |
|---|
| 961 | |
|---|
| 962 | ! Arguments |
|---|
| 963 | integer n,m,opt ! I |
|---|
| 964 | real zz(m),z(n) ! I |
|---|
| 965 | real y1(m),y2(m) ! O |
|---|
| 966 | real x1(n),x2(n) ! I |
|---|
| 967 | |
|---|
| 968 | |
|---|
| 969 | ! Local variables |
|---|
| 970 | integer i, j |
|---|
| 971 | real factor |
|---|
| 972 | real zaux |
|---|
| 973 | |
|---|
| 974 | !!!! |
|---|
| 975 | |
|---|
| 976 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 977 | |
|---|
| 978 | do 1,i=1,m ! |
|---|
| 979 | |
|---|
| 980 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 981 | zaux = zz(i) |
|---|
| 982 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 983 | zaux=z(1) |
|---|
| 984 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 985 | zaux=z(n) |
|---|
| 986 | endif |
|---|
| 987 | call hunt ( z,n, zaux, j ) |
|---|
| 988 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 989 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 990 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 991 | stop ' interhunt2/ Interpolat error. zz out of limits.' |
|---|
| 992 | endif |
|---|
| 993 | |
|---|
| 994 | ! Perform interpolation |
|---|
| 995 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 996 | if (opt.eq.1) then |
|---|
| 997 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 998 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 999 | else |
|---|
| 1000 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1001 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1002 | end if |
|---|
| 1003 | |
|---|
| 1004 | 1 continue |
|---|
| 1005 | |
|---|
| 1006 | return |
|---|
| 1007 | end |
|---|
| 1008 | |
|---|
| 1009 | |
|---|
| 1010 | c *** interhunt5veces.f *** |
|---|
| 1011 | |
|---|
| 1012 | c *********************************************************************** |
|---|
| 1013 | |
|---|
| 1014 | subroutine interhunt5veces ( y1,y2,y3,y4,y5, zz,m, |
|---|
| 1015 | @ x1,x2,x3,x4,x5, z,n, opt) |
|---|
| 1016 | |
|---|
| 1017 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1018 | c input values: y(n) at z(n) |
|---|
| 1019 | c output values: yy(m) at zz(m) |
|---|
| 1020 | c options: 1 -> lineal |
|---|
| 1021 | c 2 -> logarithmic |
|---|
| 1022 | c *********************************************************************** |
|---|
| 1023 | |
|---|
| 1024 | implicit none |
|---|
| 1025 | ! Arguments |
|---|
| 1026 | integer n,m,opt ! I |
|---|
| 1027 | real zz(m),z(n) ! I |
|---|
| 1028 | real y1(m),y2(m),y3(m),y4(m),y5(m) ! O |
|---|
| 1029 | real x1(n),x2(n),x3(n),x4(n),x5(n) ! I |
|---|
| 1030 | |
|---|
| 1031 | |
|---|
| 1032 | ! Local variables |
|---|
| 1033 | integer i, j |
|---|
| 1034 | real factor |
|---|
| 1035 | real zaux |
|---|
| 1036 | |
|---|
| 1037 | !!!! |
|---|
| 1038 | |
|---|
| 1039 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1040 | |
|---|
| 1041 | do 1,i=1,m ! |
|---|
| 1042 | |
|---|
| 1043 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1044 | zaux = zz(i) |
|---|
| 1045 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 1046 | zaux=z(1) |
|---|
| 1047 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 1048 | zaux=z(n) |
|---|
| 1049 | endif |
|---|
| 1050 | call hunt ( z,n, zaux, j ) |
|---|
| 1051 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1052 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1053 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1054 | stop ' interhunt5/ Interpolat error. zz out of limits.' |
|---|
| 1055 | endif |
|---|
| 1056 | |
|---|
| 1057 | ! Perform interpolation |
|---|
| 1058 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1059 | if (opt.eq.1) then |
|---|
| 1060 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1061 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1062 | y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor |
|---|
| 1063 | y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor |
|---|
| 1064 | y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor |
|---|
| 1065 | else |
|---|
| 1066 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1067 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1068 | y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor ) |
|---|
| 1069 | y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor ) |
|---|
| 1070 | y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor ) |
|---|
| 1071 | end if |
|---|
| 1072 | |
|---|
| 1073 | 1 continue |
|---|
| 1074 | |
|---|
| 1075 | return |
|---|
| 1076 | end |
|---|
| 1077 | |
|---|
| 1078 | |
|---|
| 1079 | |
|---|
| 1080 | c *** interhuntdp3veces.f *** |
|---|
| 1081 | |
|---|
| 1082 | c *********************************************************************** |
|---|
| 1083 | |
|---|
| 1084 | subroutine interhuntdp3veces ( y1,y2,y3, zz,m, |
|---|
| 1085 | @ x1,x2,x3, z,n, opt) |
|---|
| 1086 | |
|---|
| 1087 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1088 | c input values: x(n) at z(n) |
|---|
| 1089 | c output values: y(m) at zz(m) |
|---|
| 1090 | c options: opt = 1 -> lineal |
|---|
| 1091 | c opt=/=1 -> logarithmic |
|---|
| 1092 | c *********************************************************************** |
|---|
| 1093 | ! Arguments |
|---|
| 1094 | integer n,m,opt ! I |
|---|
| 1095 | real*8 zz(m),z(n) ! I |
|---|
| 1096 | real*8 y1(m),y2(m),y3(m) ! O |
|---|
| 1097 | real*8 x1(n),x2(n),x3(n) ! I |
|---|
| 1098 | |
|---|
| 1099 | |
|---|
| 1100 | ! Local variables |
|---|
| 1101 | integer i, j |
|---|
| 1102 | real*8 factor |
|---|
| 1103 | real*8 zaux |
|---|
| 1104 | |
|---|
| 1105 | !!!! |
|---|
| 1106 | |
|---|
| 1107 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1108 | |
|---|
| 1109 | do 1,i=1,m ! |
|---|
| 1110 | |
|---|
| 1111 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1112 | zaux = zz(i) |
|---|
| 1113 | if (abs(zaux-z(1)).le.0.01d0) then |
|---|
| 1114 | zaux=z(1) |
|---|
| 1115 | elseif (abs(zaux-z(n)).le.0.01d0) then |
|---|
| 1116 | zaux=z(n) |
|---|
| 1117 | endif |
|---|
| 1118 | call huntdp ( z,n, zaux, j ) |
|---|
| 1119 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1120 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1121 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1122 | stop ' INTERHUNTDP3/ Interpolat error. zz out of limits.' |
|---|
| 1123 | endif |
|---|
| 1124 | |
|---|
| 1125 | ! Perform interpolation |
|---|
| 1126 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1127 | if (opt.eq.1) then |
|---|
| 1128 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1129 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1130 | y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor |
|---|
| 1131 | else |
|---|
| 1132 | y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor ) |
|---|
| 1133 | y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor ) |
|---|
| 1134 | y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor ) |
|---|
| 1135 | end if |
|---|
| 1136 | |
|---|
| 1137 | 1 continue |
|---|
| 1138 | |
|---|
| 1139 | return |
|---|
| 1140 | end |
|---|
| 1141 | |
|---|
| 1142 | |
|---|
| 1143 | c *** interhuntdp4veces.f *** |
|---|
| 1144 | |
|---|
| 1145 | c *********************************************************************** |
|---|
| 1146 | |
|---|
| 1147 | subroutine interhuntdp4veces ( y1,y2,y3,y4, zz,m, |
|---|
| 1148 | @ x1,x2,x3,x4, z,n, opt) |
|---|
| 1149 | |
|---|
| 1150 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1151 | c input values: x1(n),x2(n),x3(n),x4(n) at z(n) |
|---|
| 1152 | c output values: y1(m),y2(m),y3(m),y4(m) at zz(m) |
|---|
| 1153 | c options: 1 -> lineal |
|---|
| 1154 | c 2 -> logarithmic |
|---|
| 1155 | c *********************************************************************** |
|---|
| 1156 | |
|---|
| 1157 | implicit none |
|---|
| 1158 | |
|---|
| 1159 | ! Arguments |
|---|
| 1160 | integer n,m,opt ! I |
|---|
| 1161 | real*8 zz(m),z(n) ! I |
|---|
| 1162 | real*8 y1(m),y2(m),y3(m),y4(m) ! O |
|---|
| 1163 | real*8 x1(n),x2(n),x3(n),x4(n) ! I |
|---|
| 1164 | |
|---|
| 1165 | |
|---|
| 1166 | ! Local variables |
|---|
| 1167 | integer i, j |
|---|
| 1168 | real*8 factor |
|---|
| 1169 | real*8 zaux |
|---|
| 1170 | |
|---|
| 1171 | !!!! |
|---|
| 1172 | |
|---|
| 1173 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1174 | |
|---|
| 1175 | do 1,i=1,m ! |
|---|
| 1176 | |
|---|
| 1177 | ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1178 | zaux = zz(i) |
|---|
| 1179 | if (abs(zaux-z(1)).le.0.01d0) then |
|---|
| 1180 | zaux=z(1) |
|---|
| 1181 | elseif (abs(zaux-z(n)).le.0.01d0) then |
|---|
| 1182 | zaux=z(n) |
|---|
| 1183 | endif |
|---|
| 1184 | call huntdp ( z,n, zaux, j ) |
|---|
| 1185 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1186 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1187 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1188 | stop ' INTERHUNTDP4/ Interpolat error. zz out of limits.' |
|---|
| 1189 | endif |
|---|
| 1190 | |
|---|
| 1191 | ! Perform interpolation |
|---|
| 1192 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1193 | if (opt.eq.1) then |
|---|
| 1194 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1195 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1196 | y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor |
|---|
| 1197 | y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor |
|---|
| 1198 | else |
|---|
| 1199 | y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor ) |
|---|
| 1200 | y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor ) |
|---|
| 1201 | y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor ) |
|---|
| 1202 | y4(i) = dexp( dlog(x4(j)) + dlog(x4(j+1)/x4(j)) * factor ) |
|---|
| 1203 | end if |
|---|
| 1204 | |
|---|
| 1205 | 1 continue |
|---|
| 1206 | |
|---|
| 1207 | return |
|---|
| 1208 | end |
|---|
| 1209 | |
|---|
| 1210 | |
|---|
| 1211 | c *** interhuntdp.f *** |
|---|
| 1212 | |
|---|
| 1213 | c *********************************************************************** |
|---|
| 1214 | |
|---|
| 1215 | subroutine interhuntdp ( y1, zz,m, |
|---|
| 1216 | @ x1, z,n, opt) |
|---|
| 1217 | |
|---|
| 1218 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1219 | c input values: x1(n) at z(n) |
|---|
| 1220 | c output values: y1(m) at zz(m) |
|---|
| 1221 | c options: 1 -> lineal |
|---|
| 1222 | c 2 -> logarithmic |
|---|
| 1223 | c *********************************************************************** |
|---|
| 1224 | |
|---|
| 1225 | implicit none |
|---|
| 1226 | |
|---|
| 1227 | ! Arguments |
|---|
| 1228 | integer n,m,opt ! I |
|---|
| 1229 | real*8 zz(m),z(n) ! I |
|---|
| 1230 | real*8 y1(m) ! O |
|---|
| 1231 | real*8 x1(n) ! I |
|---|
| 1232 | |
|---|
| 1233 | |
|---|
| 1234 | ! Local variables |
|---|
| 1235 | integer i, j |
|---|
| 1236 | real*8 factor |
|---|
| 1237 | real*8 zaux |
|---|
| 1238 | |
|---|
| 1239 | !!!! |
|---|
| 1240 | |
|---|
| 1241 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1242 | |
|---|
| 1243 | do 1,i=1,m ! |
|---|
| 1244 | |
|---|
| 1245 | ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1246 | zaux = zz(i) |
|---|
| 1247 | if (abs(zaux-z(1)).le.0.01d0) then |
|---|
| 1248 | zaux=z(1) |
|---|
| 1249 | elseif (abs(zaux-z(n)).le.0.01d0) then |
|---|
| 1250 | zaux=z(n) |
|---|
| 1251 | endif |
|---|
| 1252 | call huntdp ( z,n, zaux, j ) |
|---|
| 1253 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1254 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1255 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1256 | stop ' INTERHUNT/ Interpolat error. zz out of limits.' |
|---|
| 1257 | endif |
|---|
| 1258 | |
|---|
| 1259 | ! Perform interpolation |
|---|
| 1260 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1261 | if (opt.eq.1) then |
|---|
| 1262 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1263 | else |
|---|
| 1264 | y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor ) |
|---|
| 1265 | end if |
|---|
| 1266 | |
|---|
| 1267 | 1 continue |
|---|
| 1268 | |
|---|
| 1269 | return |
|---|
| 1270 | end |
|---|
| 1271 | |
|---|
| 1272 | |
|---|
| 1273 | c *** interhunt.f *** |
|---|
| 1274 | |
|---|
| 1275 | c*********************************************************************** |
|---|
| 1276 | |
|---|
| 1277 | subroutine interhunt ( y1, zz,m, |
|---|
| 1278 | @ x1, z,n, opt) |
|---|
| 1279 | |
|---|
| 1280 | c interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1281 | c input values: x1(n) at z(n) |
|---|
| 1282 | c output values: y1(m) at zz(m) |
|---|
| 1283 | c options: 1 -> lineal |
|---|
| 1284 | c 2 -> logarithmic |
|---|
| 1285 | c*********************************************************************** |
|---|
| 1286 | |
|---|
| 1287 | implicit none |
|---|
| 1288 | |
|---|
| 1289 | ! Arguments |
|---|
| 1290 | integer n,m,opt ! I |
|---|
| 1291 | real zz(m),z(n) ! I |
|---|
| 1292 | real y1(m) ! O |
|---|
| 1293 | real x1(n) ! I |
|---|
| 1294 | |
|---|
| 1295 | |
|---|
| 1296 | ! Local variables |
|---|
| 1297 | integer i, j |
|---|
| 1298 | real factor |
|---|
| 1299 | real zaux |
|---|
| 1300 | |
|---|
| 1301 | !!!! |
|---|
| 1302 | |
|---|
| 1303 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1304 | |
|---|
| 1305 | do 1,i=1,m ! |
|---|
| 1306 | |
|---|
| 1307 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1308 | zaux = zz(i) |
|---|
| 1309 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 1310 | zaux=z(1) |
|---|
| 1311 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 1312 | zaux=z(n) |
|---|
| 1313 | endif |
|---|
| 1314 | call hunt ( z,n, zaux, j ) |
|---|
| 1315 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1316 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1317 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1318 | stop ' interhunt/ Interpolat error. z out of limits.' |
|---|
| 1319 | endif |
|---|
| 1320 | |
|---|
| 1321 | ! Perform interpolation |
|---|
| 1322 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1323 | if (opt.eq.1) then |
|---|
| 1324 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1325 | else |
|---|
| 1326 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1327 | end if |
|---|
| 1328 | |
|---|
| 1329 | |
|---|
| 1330 | 1 continue |
|---|
| 1331 | |
|---|
| 1332 | return |
|---|
| 1333 | end |
|---|
| 1334 | |
|---|
| 1335 | |
|---|
| 1336 | c *** interhuntlimits2veces.f *** |
|---|
| 1337 | |
|---|
| 1338 | c*********************************************************************** |
|---|
| 1339 | |
|---|
| 1340 | subroutine interhuntlimits2veces |
|---|
| 1341 | @ ( y1,y2, zz,m, limite1,limite2, |
|---|
| 1342 | @ x1,x2, z,n, opt) |
|---|
| 1343 | |
|---|
| 1344 | c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1345 | c Input values: x1,x2(n) at z(n) |
|---|
| 1346 | c Output values: |
|---|
| 1347 | c y1,y2(m) at zz(m) pero solo entre los indices de zz |
|---|
| 1348 | c siguientes: [limite1,limite2] |
|---|
| 1349 | c Options: 1 -> linear in z and linear in x |
|---|
| 1350 | c 2 -> linear in z and logarithmic in x |
|---|
| 1351 | c 3 -> logarithmic in z and linear in x |
|---|
| 1352 | c 4 -> logarithmic in z and logaritmic in x |
|---|
| 1353 | c |
|---|
| 1354 | c NOTAS: Esta subrutina extiende y generaliza la usual |
|---|
| 1355 | c "interhunt5veces" en 2 direcciones: |
|---|
| 1356 | c - la condicion en los limites es que zz(limite1:limite2) |
|---|
| 1357 | c esté dentro de los limites de z (pero quizas no todo zz) |
|---|
| 1358 | c - se han añadido 3 opciones mas al caso de interpolacion |
|---|
| 1359 | c logaritmica, ahora se hace en log de z, de x o de ambos. |
|---|
| 1360 | c Notese que esta subrutina engloba a la interhunt5veces |
|---|
| 1361 | c ( esta es reproducible haciendo limite1=1 y limite2=m |
|---|
| 1362 | c y usando una de las 2 primeras opciones opt=1,2 ) |
|---|
| 1363 | c*********************************************************************** |
|---|
| 1364 | |
|---|
| 1365 | implicit none |
|---|
| 1366 | |
|---|
| 1367 | ! Arguments |
|---|
| 1368 | integer n,m,opt, limite1,limite2 ! I |
|---|
| 1369 | real zz(m),z(n) ! I |
|---|
| 1370 | real y1(m),y2(m) ! O |
|---|
| 1371 | real x1(n),x2(n) ! I |
|---|
| 1372 | |
|---|
| 1373 | |
|---|
| 1374 | ! Local variables |
|---|
| 1375 | integer i, j |
|---|
| 1376 | real factor |
|---|
| 1377 | real zaux |
|---|
| 1378 | |
|---|
| 1379 | !!!! |
|---|
| 1380 | |
|---|
| 1381 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1382 | |
|---|
| 1383 | do 1,i=limite1,limite2 |
|---|
| 1384 | |
|---|
| 1385 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1386 | zaux = zz(i) |
|---|
| 1387 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 1388 | zaux=z(1) |
|---|
| 1389 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 1390 | zaux=z(n) |
|---|
| 1391 | endif |
|---|
| 1392 | call hunt ( z,n, zaux, j ) |
|---|
| 1393 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1394 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1395 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1396 | stop ' interhuntlimits/ Interpolat error. z out of limits.' |
|---|
| 1397 | endif |
|---|
| 1398 | |
|---|
| 1399 | ! Perform interpolation |
|---|
| 1400 | if (opt.eq.1) then |
|---|
| 1401 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1402 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1403 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1404 | |
|---|
| 1405 | elseif (opt.eq.2) then |
|---|
| 1406 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1407 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1408 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1409 | elseif (opt.eq.3) then |
|---|
| 1410 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1411 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1412 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1413 | elseif (opt.eq.4) then |
|---|
| 1414 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1415 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1416 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1417 | end if |
|---|
| 1418 | |
|---|
| 1419 | |
|---|
| 1420 | 1 continue |
|---|
| 1421 | |
|---|
| 1422 | return |
|---|
| 1423 | end |
|---|
| 1424 | |
|---|
| 1425 | |
|---|
| 1426 | c *** interhuntlimits5veces.f *** |
|---|
| 1427 | |
|---|
| 1428 | c*********************************************************************** |
|---|
| 1429 | |
|---|
| 1430 | subroutine interhuntlimits5veces |
|---|
| 1431 | @ ( y1,y2,y3,y4,y5, zz,m, limite1,limite2, |
|---|
| 1432 | @ x1,x2,x3,x4,x5, z,n, opt) |
|---|
| 1433 | |
|---|
| 1434 | c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1435 | c Input values: x1,x2,..,x5(n) at z(n) |
|---|
| 1436 | c Output values: |
|---|
| 1437 | c y1,y2,...,y5(m) at zz(m) pero solo entre los indices de zz |
|---|
| 1438 | c siguientes: [limite1,limite2] |
|---|
| 1439 | c Options: 1 -> linear in z and linear in x |
|---|
| 1440 | c 2 -> linear in z and logarithmic in x |
|---|
| 1441 | c 3 -> logarithmic in z and linear in x |
|---|
| 1442 | c 4 -> logarithmic in z and logaritmic in x |
|---|
| 1443 | c |
|---|
| 1444 | c NOTAS: Esta subrutina extiende y generaliza la usual |
|---|
| 1445 | c "interhunt5veces" en 2 direcciones: |
|---|
| 1446 | c - la condicion en los limites es que zz(limite1:limite2) |
|---|
| 1447 | c esté dentro de los limites de z (pero quizas no todo zz) |
|---|
| 1448 | c - se han añadido 3 opciones mas al caso de interpolacion |
|---|
| 1449 | c logaritmica, ahora se hace en log de z, de x o de ambos. |
|---|
| 1450 | c Notese que esta subrutina engloba a la interhunt5veces |
|---|
| 1451 | c ( esta es reproducible haciendo limite1=1 y limite2=m |
|---|
| 1452 | c y usando una de las 2 primeras opciones opt=1,2 ) |
|---|
| 1453 | c*********************************************************************** |
|---|
| 1454 | |
|---|
| 1455 | implicit none |
|---|
| 1456 | |
|---|
| 1457 | ! Arguments |
|---|
| 1458 | integer n,m,opt, limite1,limite2 ! I |
|---|
| 1459 | real zz(m),z(n) ! I |
|---|
| 1460 | real y1(m),y2(m),y3(m),y4(m),y5(m) ! O |
|---|
| 1461 | real x1(n),x2(n),x3(n),x4(n),x5(n) ! I |
|---|
| 1462 | |
|---|
| 1463 | |
|---|
| 1464 | ! Local variables |
|---|
| 1465 | integer i, j |
|---|
| 1466 | real factor |
|---|
| 1467 | real zaux |
|---|
| 1468 | |
|---|
| 1469 | !!!! |
|---|
| 1470 | |
|---|
| 1471 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1472 | |
|---|
| 1473 | do 1,i=limite1,limite2 |
|---|
| 1474 | |
|---|
| 1475 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1476 | zaux = zz(i) |
|---|
| 1477 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 1478 | zaux=z(1) |
|---|
| 1479 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 1480 | zaux=z(n) |
|---|
| 1481 | endif |
|---|
| 1482 | call hunt ( z,n, zaux, j ) |
|---|
| 1483 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1484 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1485 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1486 | stop ' interhuntlimits/ Interpolat error. z out of limits.' |
|---|
| 1487 | endif |
|---|
| 1488 | |
|---|
| 1489 | ! Perform interpolation |
|---|
| 1490 | if (opt.eq.1) then |
|---|
| 1491 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1492 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1493 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1494 | y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor |
|---|
| 1495 | y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor |
|---|
| 1496 | y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor |
|---|
| 1497 | elseif (opt.eq.2) then |
|---|
| 1498 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1499 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1500 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1501 | y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor ) |
|---|
| 1502 | y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor ) |
|---|
| 1503 | y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor ) |
|---|
| 1504 | elseif (opt.eq.3) then |
|---|
| 1505 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1506 | y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor |
|---|
| 1507 | y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor |
|---|
| 1508 | y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor |
|---|
| 1509 | y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor |
|---|
| 1510 | y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor |
|---|
| 1511 | elseif (opt.eq.4) then |
|---|
| 1512 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1513 | y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor ) |
|---|
| 1514 | y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor ) |
|---|
| 1515 | y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor ) |
|---|
| 1516 | y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor ) |
|---|
| 1517 | y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor ) |
|---|
| 1518 | end if |
|---|
| 1519 | |
|---|
| 1520 | |
|---|
| 1521 | 1 continue |
|---|
| 1522 | |
|---|
| 1523 | return |
|---|
| 1524 | end |
|---|
| 1525 | |
|---|
| 1526 | |
|---|
| 1527 | |
|---|
| 1528 | c *** interhuntlimits.f *** |
|---|
| 1529 | |
|---|
| 1530 | c*********************************************************************** |
|---|
| 1531 | |
|---|
| 1532 | subroutine interhuntlimits ( y, zz,m, limite1,limite2, |
|---|
| 1533 | @ x, z,n, opt) |
|---|
| 1534 | |
|---|
| 1535 | c Interpolation soubroutine basada en Numerical Recipes HUNT.FOR |
|---|
| 1536 | c Input values: x(n) at z(n) |
|---|
| 1537 | c Output values: y(m) at zz(m) pero solo entre los indices de zz |
|---|
| 1538 | c siguientes: [limite1,limite2] |
|---|
| 1539 | c Options: 1 -> linear in z and linear in x |
|---|
| 1540 | c 2 -> linear in z and logarithmic in x |
|---|
| 1541 | c 3 -> logarithmic in z and linear in x |
|---|
| 1542 | c 4 -> logarithmic in z and logaritmic in x |
|---|
| 1543 | c |
|---|
| 1544 | c NOTAS: Esta subrutina extiende y generaliza la usual "interhunt" |
|---|
| 1545 | c en 2 direcciones: |
|---|
| 1546 | c - la condicion en los limites es que zz(limite1:limite2) |
|---|
| 1547 | c esté dentro de los limites de z (pero quizas no todo zz) |
|---|
| 1548 | c - se han añadido 3 opciones mas al caso de interpolacion |
|---|
| 1549 | c logaritmica, ahora se hace en log de z, de x o de ambos. |
|---|
| 1550 | c Notese que esta subrutina engloba a la usual interhunt |
|---|
| 1551 | c ( esta es reproducible haciendo limite1=1 y limite2=m |
|---|
| 1552 | c y usando una de las 2 primeras opciones opt=1,2 ) |
|---|
| 1553 | c*********************************************************************** |
|---|
| 1554 | |
|---|
| 1555 | implicit none |
|---|
| 1556 | |
|---|
| 1557 | ! Arguments |
|---|
| 1558 | integer n,m,opt, limite1,limite2 ! I |
|---|
| 1559 | real zz(m),z(n) ! I |
|---|
| 1560 | real y(m) ! O |
|---|
| 1561 | real x(n) ! I |
|---|
| 1562 | |
|---|
| 1563 | |
|---|
| 1564 | ! Local variables |
|---|
| 1565 | integer i, j |
|---|
| 1566 | real factor |
|---|
| 1567 | real zaux |
|---|
| 1568 | |
|---|
| 1569 | !!!! |
|---|
| 1570 | |
|---|
| 1571 | j = 1 ! initial first guess (=n/2 is anothr pssblty) |
|---|
| 1572 | |
|---|
| 1573 | do 1,i=limite1,limite2 |
|---|
| 1574 | |
|---|
| 1575 | ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)] |
|---|
| 1576 | zaux = zz(i) |
|---|
| 1577 | if (abs(zaux-z(1)).le.0.01) then |
|---|
| 1578 | zaux=z(1) |
|---|
| 1579 | elseif (abs(zaux-z(n)).le.0.01) then |
|---|
| 1580 | zaux=z(n) |
|---|
| 1581 | endif |
|---|
| 1582 | call hunt ( z,n, zaux, j ) |
|---|
| 1583 | if ( j.eq.0 .or. j.eq.n ) then |
|---|
| 1584 | write (*,*) ' HUNT/ Limits input grid:', z(1),z(n) |
|---|
| 1585 | write (*,*) ' HUNT/ location in new grid:', zz(i) |
|---|
| 1586 | stop ' interhuntlimits/ Interpolat error. z out of limits.' |
|---|
| 1587 | endif |
|---|
| 1588 | |
|---|
| 1589 | ! Perform interpolation |
|---|
| 1590 | if (opt.eq.1) then |
|---|
| 1591 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1592 | y(i) = x(j) + (x(j+1)-x(j)) * factor |
|---|
| 1593 | elseif (opt.eq.2) then |
|---|
| 1594 | factor = (zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1595 | y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor ) |
|---|
| 1596 | elseif (opt.eq.3) then |
|---|
| 1597 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1598 | y(i) = x(j) + (x(j+1)-x(j)) * factor |
|---|
| 1599 | elseif (opt.eq.4) then |
|---|
| 1600 | factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j))) |
|---|
| 1601 | y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor ) |
|---|
| 1602 | end if |
|---|
| 1603 | |
|---|
| 1604 | |
|---|
| 1605 | 1 continue |
|---|
| 1606 | |
|---|
| 1607 | return |
|---|
| 1608 | end |
|---|
| 1609 | |
|---|
| 1610 | |
|---|
| 1611 | c *** lubksb_dp.f *** |
|---|
| 1612 | |
|---|
| 1613 | subroutine lubksb_dp(a,n,np,indx,b) |
|---|
| 1614 | |
|---|
| 1615 | implicit none |
|---|
| 1616 | |
|---|
| 1617 | integer,intent(in) :: n,np |
|---|
| 1618 | real*8,intent(in) :: a(np,np) |
|---|
| 1619 | integer,intent(in) :: indx(n) |
|---|
| 1620 | real*8,intent(out) :: b(n) |
|---|
| 1621 | |
|---|
| 1622 | real*8 sum |
|---|
| 1623 | integer ii, ll, i, j |
|---|
| 1624 | |
|---|
| 1625 | ii=0 |
|---|
| 1626 | do 12 i=1,n |
|---|
| 1627 | ll=indx(i) |
|---|
| 1628 | sum=b(ll) |
|---|
| 1629 | b(ll)=b(i) |
|---|
| 1630 | if (ii.ne.0)then |
|---|
| 1631 | do 11 j=ii,i-1 |
|---|
| 1632 | sum=sum-a(i,j)*b(j) |
|---|
| 1633 | 11 continue |
|---|
| 1634 | else if (sum.ne.0.0) then |
|---|
| 1635 | ii=i |
|---|
| 1636 | endif |
|---|
| 1637 | b(i)=sum |
|---|
| 1638 | 12 continue |
|---|
| 1639 | do 14 i=n,1,-1 |
|---|
| 1640 | sum=b(i) |
|---|
| 1641 | if(i.lt.n)then |
|---|
| 1642 | do 13 j=i+1,n |
|---|
| 1643 | sum=sum-a(i,j)*b(j) |
|---|
| 1644 | 13 continue |
|---|
| 1645 | endif |
|---|
| 1646 | b(i)=sum/a(i,i) |
|---|
| 1647 | 14 continue |
|---|
| 1648 | return |
|---|
| 1649 | end |
|---|
| 1650 | |
|---|
| 1651 | |
|---|
| 1652 | c *** ludcmp_dp.f *** |
|---|
| 1653 | |
|---|
| 1654 | subroutine ludcmp_dp(a,n,np,indx,d) |
|---|
| 1655 | |
|---|
| 1656 | implicit none |
|---|
| 1657 | |
|---|
| 1658 | integer,intent(in) :: n, np |
|---|
| 1659 | real*8,intent(inout) :: a(np,np) |
|---|
| 1660 | real*8,intent(out) :: d |
|---|
| 1661 | integer,intent(out) :: indx(n) |
|---|
| 1662 | |
|---|
| 1663 | integer nmax, i, j, k, imax |
|---|
| 1664 | real*8 tiny |
|---|
| 1665 | parameter (nmax=100,tiny=1.0d-20) |
|---|
| 1666 | real*8 vv(nmax), aamax, sum, dum |
|---|
| 1667 | |
|---|
| 1668 | |
|---|
| 1669 | d=1.0d0 |
|---|
| 1670 | do 12 i=1,n |
|---|
| 1671 | aamax=0.0d0 |
|---|
| 1672 | do 11 j=1,n |
|---|
| 1673 | if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) |
|---|
| 1674 | 11 continue |
|---|
| 1675 | if (aamax.eq.0.0) then |
|---|
| 1676 | write(*,*) 'ludcmp_dp: singular matrix!' |
|---|
| 1677 | stop |
|---|
| 1678 | endif |
|---|
| 1679 | vv(i)=1.0d0/aamax |
|---|
| 1680 | 12 continue |
|---|
| 1681 | do 19 j=1,n |
|---|
| 1682 | if (j.gt.1) then |
|---|
| 1683 | do 14 i=1,j-1 |
|---|
| 1684 | sum=a(i,j) |
|---|
| 1685 | if (i.gt.1)then |
|---|
| 1686 | do 13 k=1,i-1 |
|---|
| 1687 | sum=sum-a(i,k)*a(k,j) |
|---|
| 1688 | 13 continue |
|---|
| 1689 | a(i,j)=sum |
|---|
| 1690 | endif |
|---|
| 1691 | 14 continue |
|---|
| 1692 | endif |
|---|
| 1693 | aamax=0.0d0 |
|---|
| 1694 | do 16 i=j,n |
|---|
| 1695 | sum=a(i,j) |
|---|
| 1696 | if (j.gt.1)then |
|---|
| 1697 | do 15 k=1,j-1 |
|---|
| 1698 | sum=sum-a(i,k)*a(k,j) |
|---|
| 1699 | 15 continue |
|---|
| 1700 | a(i,j)=sum |
|---|
| 1701 | endif |
|---|
| 1702 | dum=vv(i)*abs(sum) |
|---|
| 1703 | if (dum.ge.aamax) then |
|---|
| 1704 | imax=i |
|---|
| 1705 | aamax=dum |
|---|
| 1706 | endif |
|---|
| 1707 | 16 continue |
|---|
| 1708 | if (j.ne.imax)then |
|---|
| 1709 | do 17 k=1,n |
|---|
| 1710 | dum=a(imax,k) |
|---|
| 1711 | a(imax,k)=a(j,k) |
|---|
| 1712 | a(j,k)=dum |
|---|
| 1713 | 17 continue |
|---|
| 1714 | d=-d |
|---|
| 1715 | vv(imax)=vv(j) |
|---|
| 1716 | endif |
|---|
| 1717 | indx(j)=imax |
|---|
| 1718 | if(j.ne.n)then |
|---|
| 1719 | if(a(j,j).eq.0.0)a(j,j)=tiny |
|---|
| 1720 | dum=1.0d0/a(j,j) |
|---|
| 1721 | do 18 i=j+1,n |
|---|
| 1722 | a(i,j)=a(i,j)*dum |
|---|
| 1723 | 18 continue |
|---|
| 1724 | endif |
|---|
| 1725 | 19 continue |
|---|
| 1726 | if(a(n,n).eq.0.0)a(n,n)=tiny |
|---|
| 1727 | return |
|---|
| 1728 | end |
|---|
| 1729 | |
|---|
| 1730 | |
|---|
| 1731 | c *** LUdec.f *** |
|---|
| 1732 | |
|---|
| 1733 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1734 | c |
|---|
| 1735 | c Solution of linear equation without inverting matrix |
|---|
| 1736 | c using LU decomposition: |
|---|
| 1737 | c AA * xx = bb AA, bb: known |
|---|
| 1738 | c xx: to be found |
|---|
| 1739 | c AA and bb are not modified in this subroutine |
|---|
| 1740 | c |
|---|
| 1741 | c MALV , Sep 2007 |
|---|
| 1742 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1743 | |
|---|
| 1744 | subroutine LUdec(xx,aa,bb,m,n) |
|---|
| 1745 | |
|---|
| 1746 | implicit none |
|---|
| 1747 | |
|---|
| 1748 | ! Arguments |
|---|
| 1749 | integer,intent(in) :: m, n |
|---|
| 1750 | real*8,intent(in) :: aa(m,m), bb(m) |
|---|
| 1751 | real*8,intent(out) :: xx(m) |
|---|
| 1752 | |
|---|
| 1753 | |
|---|
| 1754 | ! Local variables |
|---|
| 1755 | real*8 a(n,n), b(n), x(n), d |
|---|
| 1756 | integer i, j, indx(n) |
|---|
| 1757 | |
|---|
| 1758 | |
|---|
| 1759 | ! Subrutinas utilizadas |
|---|
| 1760 | ! ludcmp_dp, lubksb_dp |
|---|
| 1761 | |
|---|
| 1762 | !!!!!!!!!!!!!!!Comienza el programa !!!!!!!!!!!!!! |
|---|
| 1763 | |
|---|
| 1764 | do i=1,n |
|---|
| 1765 | b(i) = bb(i+1) |
|---|
| 1766 | do j=1,n |
|---|
| 1767 | a(i,j) = aa(i+1,j+1) |
|---|
| 1768 | enddo |
|---|
| 1769 | enddo |
|---|
| 1770 | |
|---|
| 1771 | ! Descomposicion de auxm1 |
|---|
| 1772 | call ludcmp_dp ( a, n, n, indx, d) |
|---|
| 1773 | |
|---|
| 1774 | ! Sustituciones foward y backwards para hallar la solucion |
|---|
| 1775 | do i=1,n |
|---|
| 1776 | x(i) = b(i) |
|---|
| 1777 | enddo |
|---|
| 1778 | call lubksb_dp( a, n, n, indx, x ) |
|---|
| 1779 | |
|---|
| 1780 | do i=1,n |
|---|
| 1781 | xx(i+1) = x(i) |
|---|
| 1782 | enddo |
|---|
| 1783 | |
|---|
| 1784 | return |
|---|
| 1785 | end |
|---|
| 1786 | |
|---|
| 1787 | |
|---|
| 1788 | c *** mat_oper.f *** |
|---|
| 1789 | |
|---|
| 1790 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1791 | |
|---|
| 1792 | c *********************************************************************** |
|---|
| 1793 | subroutine unit(a,n) |
|---|
| 1794 | c store the unit value in the diagonal of a |
|---|
| 1795 | c *********************************************************************** |
|---|
| 1796 | implicit none |
|---|
| 1797 | real*8 a(n,n) |
|---|
| 1798 | integer n,i,j,k |
|---|
| 1799 | do 1,i=2,n-1 |
|---|
| 1800 | do 2,j=2,n-1 |
|---|
| 1801 | if(i.eq.j) then |
|---|
| 1802 | a(i,j) = 1.d0 |
|---|
| 1803 | else |
|---|
| 1804 | a(i,j)=0.0d0 |
|---|
| 1805 | end if |
|---|
| 1806 | 2 continue |
|---|
| 1807 | 1 continue |
|---|
| 1808 | do k=1,n |
|---|
| 1809 | a(n,k) = 0.0d0 |
|---|
| 1810 | a(1,k) = 0.0d0 |
|---|
| 1811 | a(k,1) = 0.0d0 |
|---|
| 1812 | a(k,n) = 0.0d0 |
|---|
| 1813 | end do |
|---|
| 1814 | return |
|---|
| 1815 | end |
|---|
| 1816 | |
|---|
| 1817 | c *********************************************************************** |
|---|
| 1818 | subroutine diago(a,v,n) |
|---|
| 1819 | c store the vector v in the diagonal elements of the square matrix a |
|---|
| 1820 | c *********************************************************************** |
|---|
| 1821 | implicit none |
|---|
| 1822 | |
|---|
| 1823 | integer n,i,j,k |
|---|
| 1824 | real*8 a(n,n),v(n) |
|---|
| 1825 | |
|---|
| 1826 | do 1,i=2,n-1 |
|---|
| 1827 | do 2,j=2,n-1 |
|---|
| 1828 | if(i.eq.j) then |
|---|
| 1829 | a(i,j) = v(i) |
|---|
| 1830 | else |
|---|
| 1831 | a(i,j)=0.0d0 |
|---|
| 1832 | end if |
|---|
| 1833 | 2 continue |
|---|
| 1834 | 1 continue |
|---|
| 1835 | do k=1,n |
|---|
| 1836 | a(n,k) = 0.0d0 |
|---|
| 1837 | a(1,k) = 0.0d0 |
|---|
| 1838 | a(k,1) = 0.0d0 |
|---|
| 1839 | a(k,n) = 0.0d0 |
|---|
| 1840 | end do |
|---|
| 1841 | return |
|---|
| 1842 | end |
|---|
| 1843 | |
|---|
| 1844 | c *********************************************************************** |
|---|
| 1845 | subroutine invdiag(a,b,n) |
|---|
| 1846 | c inverse of a diagonal matrix |
|---|
| 1847 | c *********************************************************************** |
|---|
| 1848 | implicit none |
|---|
| 1849 | |
|---|
| 1850 | integer n,i,j,k |
|---|
| 1851 | real*8 a(n,n),b(n,n) |
|---|
| 1852 | |
|---|
| 1853 | do 1,i=2,n-1 |
|---|
| 1854 | do 2,j=2,n-1 |
|---|
| 1855 | if (i.eq.j) then |
|---|
| 1856 | a(i,j) = 1.d0/b(i,i) |
|---|
| 1857 | else |
|---|
| 1858 | a(i,j)=0.0d0 |
|---|
| 1859 | end if |
|---|
| 1860 | 2 continue |
|---|
| 1861 | 1 continue |
|---|
| 1862 | do k=1,n |
|---|
| 1863 | a(n,k) = 0.0d0 |
|---|
| 1864 | a(1,k) = 0.0d0 |
|---|
| 1865 | a(k,1) = 0.0d0 |
|---|
| 1866 | a(k,n) = 0.0d0 |
|---|
| 1867 | end do |
|---|
| 1868 | return |
|---|
| 1869 | end |
|---|
| 1870 | |
|---|
| 1871 | |
|---|
| 1872 | c *********************************************************************** |
|---|
| 1873 | subroutine samem (a,m,n) |
|---|
| 1874 | c store the matrix m in the matrix a |
|---|
| 1875 | c *********************************************************************** |
|---|
| 1876 | implicit none |
|---|
| 1877 | real*8 a(n,n),m(n,n) |
|---|
| 1878 | integer n,i,j,k |
|---|
| 1879 | do 1,i=2,n-1 |
|---|
| 1880 | do 2,j=2,n-1 |
|---|
| 1881 | a(i,j) = m(i,j) |
|---|
| 1882 | 2 continue |
|---|
| 1883 | 1 continue |
|---|
| 1884 | do k=1,n |
|---|
| 1885 | a(n,k) = 0.0d0 |
|---|
| 1886 | a(1,k) = 0.0d0 |
|---|
| 1887 | a(k,1) = 0.0d0 |
|---|
| 1888 | a(k,n) = 0.0d0 |
|---|
| 1889 | end do |
|---|
| 1890 | return |
|---|
| 1891 | end |
|---|
| 1892 | |
|---|
| 1893 | |
|---|
| 1894 | c *********************************************************************** |
|---|
| 1895 | subroutine mulmv(a,b,c,n) |
|---|
| 1896 | c do a(i)=b(i,j)*c(j). a, b, and c must be distint |
|---|
| 1897 | c *********************************************************************** |
|---|
| 1898 | implicit none |
|---|
| 1899 | |
|---|
| 1900 | integer n,i,j |
|---|
| 1901 | real*8 a(n),b(n,n),c(n),sum |
|---|
| 1902 | |
|---|
| 1903 | do 1,i=2,n-1 |
|---|
| 1904 | sum=0.0d0 |
|---|
| 1905 | do 2,j=2,n-1 |
|---|
| 1906 | sum = sum + b(i,j) * c(j) |
|---|
| 1907 | 2 continue |
|---|
| 1908 | a(i)=sum |
|---|
| 1909 | 1 continue |
|---|
| 1910 | a(1) = 0.0d0 |
|---|
| 1911 | a(n) = 0.0d0 |
|---|
| 1912 | return |
|---|
| 1913 | end |
|---|
| 1914 | |
|---|
| 1915 | |
|---|
| 1916 | c *********************************************************************** |
|---|
| 1917 | subroutine trucodiag(a,b,c,d,e,n) |
|---|
| 1918 | c inputs: matrices b,c,d,e |
|---|
| 1919 | c output: matriz diagonal a |
|---|
| 1920 | c Operacion a realizar: a = b * c^(-1) * d + e |
|---|
| 1921 | c La matriz c va a ser invertida |
|---|
| 1922 | c Todas las matrices de entrada son diagonales excepto b |
|---|
| 1923 | c Aprovechamos esa condicion para invertir c, acelerar el calculo, y |
|---|
| 1924 | c ademas, para forzar que a sea diagonal |
|---|
| 1925 | c *********************************************************************** |
|---|
| 1926 | implicit none |
|---|
| 1927 | real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n), sum |
|---|
| 1928 | integer n,i,j,k |
|---|
| 1929 | do 1,i=2,n-1 |
|---|
| 1930 | sum=0.0d0 |
|---|
| 1931 | do 2,j=2,n-1 |
|---|
| 1932 | sum=sum+ b(i,j) * d(j,j)/c(j,j) |
|---|
| 1933 | 2 continue |
|---|
| 1934 | a(i,i) = sum + e(i,i) |
|---|
| 1935 | 1 continue |
|---|
| 1936 | do k=1,n |
|---|
| 1937 | a(n,k) = 0.0d0 |
|---|
| 1938 | a(1,k) = 0.0d0 |
|---|
| 1939 | a(k,1) = 0.0d0 |
|---|
| 1940 | a(k,n) = 0.0d0 |
|---|
| 1941 | end do |
|---|
| 1942 | return |
|---|
| 1943 | end |
|---|
| 1944 | |
|---|
| 1945 | |
|---|
| 1946 | c *********************************************************************** |
|---|
| 1947 | subroutine trucommvv(v,b,c,u,w,n) |
|---|
| 1948 | c inputs: matrices b,c , vectores u,w |
|---|
| 1949 | c output: vector v |
|---|
| 1950 | c Operacion a realizar: v = b * c^(-1) * u + w |
|---|
| 1951 | c La matriz c va a ser invertida |
|---|
| 1952 | c c es diagonal, b no |
|---|
| 1953 | c Aprovechamos esa condicion para invertir c, y acelerar el calculo |
|---|
| 1954 | c *********************************************************************** |
|---|
| 1955 | implicit none |
|---|
| 1956 | real*8 v(n),b(n,n),c(n,n),u(n),w(n), sum |
|---|
| 1957 | integer n,i,j |
|---|
| 1958 | do 1,i=2,n-1 |
|---|
| 1959 | sum=0.0d0 |
|---|
| 1960 | do 2,j=2,n-1 |
|---|
| 1961 | sum=sum+ b(i,j) * u(j)/c(j,j) |
|---|
| 1962 | 2 continue |
|---|
| 1963 | v(i) = sum + w(i) |
|---|
| 1964 | 1 continue |
|---|
| 1965 | v(1) = 0.d0 |
|---|
| 1966 | v(n) = 0.d0 |
|---|
| 1967 | return |
|---|
| 1968 | end |
|---|
| 1969 | |
|---|
| 1970 | |
|---|
| 1971 | c *********************************************************************** |
|---|
| 1972 | subroutine sypvmv(v,u,c,w,n) |
|---|
| 1973 | c inputs: matriz diagonal c , vectores u,w |
|---|
| 1974 | c output: vector v |
|---|
| 1975 | c Operacion a realizar: v = u + c * w |
|---|
| 1976 | c *********************************************************************** |
|---|
| 1977 | implicit none |
|---|
| 1978 | real*8 v(n),u(n),c(n,n),w(n) |
|---|
| 1979 | integer n,i |
|---|
| 1980 | do 1,i=2,n-1 |
|---|
| 1981 | v(i)= u(i) + c(i,i) * w(i) |
|---|
| 1982 | 1 continue |
|---|
| 1983 | v(1) = 0.0d0 |
|---|
| 1984 | v(n) = 0.0d0 |
|---|
| 1985 | return |
|---|
| 1986 | end |
|---|
| 1987 | |
|---|
| 1988 | |
|---|
| 1989 | c *********************************************************************** |
|---|
| 1990 | subroutine sumvv(a,b,c,n) |
|---|
| 1991 | c a(i)=b(i)+c(i) |
|---|
| 1992 | c *********************************************************************** |
|---|
| 1993 | implicit none |
|---|
| 1994 | |
|---|
| 1995 | integer n,i |
|---|
| 1996 | real*8 a(n),b(n),c(n) |
|---|
| 1997 | |
|---|
| 1998 | do 1,i=2,n-1 |
|---|
| 1999 | a(i)= b(i) + c(i) |
|---|
| 2000 | 1 continue |
|---|
| 2001 | a(1) = 0.0d0 |
|---|
| 2002 | a(n) = 0.0d0 |
|---|
| 2003 | return |
|---|
| 2004 | end |
|---|
| 2005 | |
|---|
| 2006 | |
|---|
| 2007 | c *********************************************************************** |
|---|
| 2008 | subroutine sypvvv(a,b,c,d,n) |
|---|
| 2009 | c a(i)=b(i)+c(i)*d(i) |
|---|
| 2010 | c *********************************************************************** |
|---|
| 2011 | implicit none |
|---|
| 2012 | real*8 a(n),b(n),c(n),d(n) |
|---|
| 2013 | integer n,i |
|---|
| 2014 | do 1,i=2,n-1 |
|---|
| 2015 | a(i)= b(i) + c(i) * d(i) |
|---|
| 2016 | 1 continue |
|---|
| 2017 | a(1) = 0.0d0 |
|---|
| 2018 | a(n) = 0.0d0 |
|---|
| 2019 | return |
|---|
| 2020 | end |
|---|
| 2021 | |
|---|
| 2022 | |
|---|
| 2023 | c *********************************************************************** |
|---|
| 2024 | ! subroutine zerom(a,n) |
|---|
| 2025 | c a(i,j)= 0.0 |
|---|
| 2026 | c *********************************************************************** |
|---|
| 2027 | ! implicit none |
|---|
| 2028 | ! integer n,i,j |
|---|
| 2029 | ! real*8 a(n,n) |
|---|
| 2030 | |
|---|
| 2031 | ! do 1,i=1,n |
|---|
| 2032 | ! do 2,j=1,n |
|---|
| 2033 | ! a(i,j) = 0.0d0 |
|---|
| 2034 | ! 2 continue |
|---|
| 2035 | ! 1 continue |
|---|
| 2036 | ! return |
|---|
| 2037 | ! end |
|---|
| 2038 | |
|---|
| 2039 | |
|---|
| 2040 | c *********************************************************************** |
|---|
| 2041 | subroutine zero4m(a,b,c,d,n) |
|---|
| 2042 | c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 |
|---|
| 2043 | c *********************************************************************** |
|---|
| 2044 | implicit none |
|---|
| 2045 | real*8 a(n,n), b(n,n), c(n,n), d(n,n) |
|---|
| 2046 | integer n |
|---|
| 2047 | a(1:n,1:n)=0.d0 |
|---|
| 2048 | b(1:n,1:n)=0.d0 |
|---|
| 2049 | c(1:n,1:n)=0.d0 |
|---|
| 2050 | d(1:n,1:n)=0.d0 |
|---|
| 2051 | ! do 1,i=1,n |
|---|
| 2052 | ! do 2,j=1,n |
|---|
| 2053 | ! a(i,j) = 0.0d0 |
|---|
| 2054 | ! b(i,j) = 0.0d0 |
|---|
| 2055 | ! c(i,j) = 0.0d0 |
|---|
| 2056 | ! d(i,j) = 0.0d0 |
|---|
| 2057 | ! 2 continue |
|---|
| 2058 | ! 1 continue |
|---|
| 2059 | return |
|---|
| 2060 | end |
|---|
| 2061 | |
|---|
| 2062 | |
|---|
| 2063 | c *********************************************************************** |
|---|
| 2064 | subroutine zero3m(a,b,c,n) |
|---|
| 2065 | c a(i,j) = b(i,j) = c(i,j) = 0.0 |
|---|
| 2066 | c ********************************************************************** |
|---|
| 2067 | implicit none |
|---|
| 2068 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 2069 | integer n |
|---|
| 2070 | a(1:n,1:n)=0.d0 |
|---|
| 2071 | b(1:n,1:n)=0.d0 |
|---|
| 2072 | c(1:n,1:n)=0.d0 |
|---|
| 2073 | ! do 1,i=1,n |
|---|
| 2074 | ! do 2,j=1,n |
|---|
| 2075 | ! a(i,j) = 0.0d0 |
|---|
| 2076 | ! b(i,j) = 0.0d0 |
|---|
| 2077 | ! c(i,j) = 0.0d0 |
|---|
| 2078 | ! 2 continue |
|---|
| 2079 | ! 1 continue |
|---|
| 2080 | return |
|---|
| 2081 | end |
|---|
| 2082 | |
|---|
| 2083 | |
|---|
| 2084 | c *********************************************************************** |
|---|
| 2085 | subroutine zero2m(a,b,n) |
|---|
| 2086 | c a(i,j) = b(i,j) = 0.0 |
|---|
| 2087 | c *********************************************************************** |
|---|
| 2088 | implicit none |
|---|
| 2089 | real*8 a(n,n), b(n,n) |
|---|
| 2090 | integer n |
|---|
| 2091 | a(1:n,1:n)=0.d0 |
|---|
| 2092 | b(1:n,1:n)=0.d0 |
|---|
| 2093 | ! do 1,i=1,n |
|---|
| 2094 | ! do 2,j=1,n |
|---|
| 2095 | ! a(i,j) = 0.0d0 |
|---|
| 2096 | ! b(i,j) = 0.0d0 |
|---|
| 2097 | ! 2 continue |
|---|
| 2098 | ! 1 continue |
|---|
| 2099 | return |
|---|
| 2100 | end |
|---|
| 2101 | |
|---|
| 2102 | |
|---|
| 2103 | c *********************************************************************** |
|---|
| 2104 | ! subroutine zerov(a,n) |
|---|
| 2105 | c a(i)= 0.0 |
|---|
| 2106 | c *********************************************************************** |
|---|
| 2107 | ! implicit none |
|---|
| 2108 | ! real*8 a(n) |
|---|
| 2109 | ! integer n,i |
|---|
| 2110 | ! do 1,i=1,n |
|---|
| 2111 | ! a(i) = 0.0d0 |
|---|
| 2112 | ! 1 continue |
|---|
| 2113 | ! return |
|---|
| 2114 | ! end |
|---|
| 2115 | |
|---|
| 2116 | |
|---|
| 2117 | c *********************************************************************** |
|---|
| 2118 | subroutine zero4v(a,b,c,d,n) |
|---|
| 2119 | c a(i) = b(i) = c(i) = d(i,j) = 0.0 |
|---|
| 2120 | c *********************************************************************** |
|---|
| 2121 | implicit none |
|---|
| 2122 | real*8 a(n), b(n), c(n), d(n) |
|---|
| 2123 | integer n |
|---|
| 2124 | a(1:n)=0.d0 |
|---|
| 2125 | b(1:n)=0.d0 |
|---|
| 2126 | c(1:n)=0.d0 |
|---|
| 2127 | d(1:n)=0.d0 |
|---|
| 2128 | ! do 1,i=1,n |
|---|
| 2129 | ! a(i) = 0.0d0 |
|---|
| 2130 | ! b(i) = 0.0d0 |
|---|
| 2131 | ! c(i) = 0.0d0 |
|---|
| 2132 | ! d(i) = 0.0d0 |
|---|
| 2133 | ! 1 continue |
|---|
| 2134 | return |
|---|
| 2135 | end |
|---|
| 2136 | |
|---|
| 2137 | |
|---|
| 2138 | c *********************************************************************** |
|---|
| 2139 | subroutine zero3v(a,b,c,n) |
|---|
| 2140 | c a(i) = b(i) = c(i) = 0.0 |
|---|
| 2141 | c *********************************************************************** |
|---|
| 2142 | implicit none |
|---|
| 2143 | real*8 a(n), b(n), c(n) |
|---|
| 2144 | integer n |
|---|
| 2145 | a(1:n)=0.d0 |
|---|
| 2146 | b(1:n)=0.d0 |
|---|
| 2147 | c(1:n)=0.d0 |
|---|
| 2148 | ! do 1,i=1,n |
|---|
| 2149 | ! a(i) = 0.0d0 |
|---|
| 2150 | ! b(i) = 0.0d0 |
|---|
| 2151 | ! c(i) = 0.0d0 |
|---|
| 2152 | ! 1 continue |
|---|
| 2153 | return |
|---|
| 2154 | end |
|---|
| 2155 | |
|---|
| 2156 | |
|---|
| 2157 | c *********************************************************************** |
|---|
| 2158 | subroutine zero2v(a,b,n) |
|---|
| 2159 | c a(i) = b(i) = 0.0 |
|---|
| 2160 | c *********************************************************************** |
|---|
| 2161 | implicit none |
|---|
| 2162 | real*8 a(n), b(n) |
|---|
| 2163 | integer n |
|---|
| 2164 | a(1:n)=0.d0 |
|---|
| 2165 | b(1:n)=0.d0 |
|---|
| 2166 | ! do 1,i=1,n |
|---|
| 2167 | ! a(i) = 0.0d0 |
|---|
| 2168 | ! b(i) = 0.0d0 |
|---|
| 2169 | ! 1 continue |
|---|
| 2170 | return |
|---|
| 2171 | end |
|---|
| 2172 | |
|---|
| 2173 | c *********************************************************************** |
|---|
| 2174 | |
|---|
| 2175 | |
|---|
| 2176 | c**************************************************************************** |
|---|
| 2177 | |
|---|
| 2178 | c *** suaviza.f *** |
|---|
| 2179 | |
|---|
| 2180 | c***************************************************************************** |
|---|
| 2181 | c |
|---|
| 2182 | subroutine suaviza ( x, n, ismooth, y ) |
|---|
| 2183 | c |
|---|
| 2184 | c x - input and return values |
|---|
| 2185 | c y - auxiliary vector |
|---|
| 2186 | c ismooth = 0 --> no smoothing is performed |
|---|
| 2187 | c ismooth = 1 --> weak smoothing (5 points, centred weighted) |
|---|
| 2188 | c ismooth = 2 --> normal smoothing (3 points, evenly weighted) |
|---|
| 2189 | c ismooth = 3 --> strong smoothing (5 points, evenly weighted) |
|---|
| 2190 | |
|---|
| 2191 | |
|---|
| 2192 | c august 1991 |
|---|
| 2193 | c***************************************************************************** |
|---|
| 2194 | |
|---|
| 2195 | implicit none |
|---|
| 2196 | |
|---|
| 2197 | integer n, imax, imin, i, ismooth |
|---|
| 2198 | real*8 x(n), y(n) |
|---|
| 2199 | c***************************************************************************** |
|---|
| 2200 | |
|---|
| 2201 | imin=1 |
|---|
| 2202 | imax=n |
|---|
| 2203 | |
|---|
| 2204 | if (ismooth.eq.0) then |
|---|
| 2205 | |
|---|
| 2206 | return |
|---|
| 2207 | |
|---|
| 2208 | elseif (ismooth.eq.1) then ! 5 points, with central weighting |
|---|
| 2209 | |
|---|
| 2210 | do i=imin,imax |
|---|
| 2211 | if(i.eq.imin)then |
|---|
| 2212 | y(i)=x(imin) |
|---|
| 2213 | elseif(i.eq.imax)then |
|---|
| 2214 | y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 |
|---|
| 2215 | elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then |
|---|
| 2216 | y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + |
|---|
| 2217 | @ x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0 |
|---|
| 2218 | else |
|---|
| 2219 | y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0 |
|---|
| 2220 | end if |
|---|
| 2221 | end do |
|---|
| 2222 | |
|---|
| 2223 | elseif (ismooth.eq.2) then ! 3 points, evenly spaced |
|---|
| 2224 | |
|---|
| 2225 | do i=imin,imax |
|---|
| 2226 | if(i.eq.imin)then |
|---|
| 2227 | y(i)=x(imin) |
|---|
| 2228 | elseif(i.eq.imax)then |
|---|
| 2229 | y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 |
|---|
| 2230 | else |
|---|
| 2231 | y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 |
|---|
| 2232 | end if |
|---|
| 2233 | end do |
|---|
| 2234 | |
|---|
| 2235 | elseif (ismooth.eq.3) then ! 5 points, evenly spaced |
|---|
| 2236 | |
|---|
| 2237 | do i=imin,imax |
|---|
| 2238 | if(i.eq.imin)then |
|---|
| 2239 | y(i) = x(imin) |
|---|
| 2240 | elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then |
|---|
| 2241 | y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 |
|---|
| 2242 | elseif(i.eq.imax)then |
|---|
| 2243 | y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0 |
|---|
| 2244 | else |
|---|
| 2245 | y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0 |
|---|
| 2246 | end if |
|---|
| 2247 | end do |
|---|
| 2248 | |
|---|
| 2249 | else |
|---|
| 2250 | |
|---|
| 2251 | write (*,*) ' Error in suaviza.f Wrong ismooth value.' |
|---|
| 2252 | stop |
|---|
| 2253 | |
|---|
| 2254 | endif |
|---|
| 2255 | |
|---|
| 2256 | c rehago el cambio, para devolver x(i) |
|---|
| 2257 | do i=imin,imax |
|---|
| 2258 | x(i)=y(i) |
|---|
| 2259 | end do |
|---|
| 2260 | |
|---|
| 2261 | return |
|---|
| 2262 | end |
|---|
| 2263 | |
|---|
| 2264 | |
|---|
| 2265 | c *********************************************************************** |
|---|
| 2266 | subroutine mulmmf90(a,b,c,n) |
|---|
| 2267 | c *********************************************************************** |
|---|
| 2268 | implicit none |
|---|
| 2269 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 2270 | integer n |
|---|
| 2271 | |
|---|
| 2272 | a=matmul(b,c) |
|---|
| 2273 | a(1,:)=0.d0 |
|---|
| 2274 | a(:,1)=0.d0 |
|---|
| 2275 | a(n,:)=0.d0 |
|---|
| 2276 | a(:,n)=0.d0 |
|---|
| 2277 | |
|---|
| 2278 | return |
|---|
| 2279 | end |
|---|
| 2280 | |
|---|
| 2281 | |
|---|
| 2282 | c *********************************************************************** |
|---|
| 2283 | subroutine resmmf90(a,b,c,n) |
|---|
| 2284 | c *********************************************************************** |
|---|
| 2285 | implicit none |
|---|
| 2286 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 2287 | integer n |
|---|
| 2288 | |
|---|
| 2289 | a=b-c |
|---|
| 2290 | a(1,:)=0.d0 |
|---|
| 2291 | a(:,1)=0.d0 |
|---|
| 2292 | a(n,:)=0.d0 |
|---|
| 2293 | a(:,n)=0.d0 |
|---|
| 2294 | |
|---|
| 2295 | return |
|---|
| 2296 | end |
|---|
| 2297 | |
|---|
| 2298 | |
|---|
| 2299 | c******************************************************************* |
|---|
| 2300 | |
|---|
| 2301 | subroutine gethist_03 (ihist) |
|---|
| 2302 | |
|---|
| 2303 | c******************************************************************* |
|---|
| 2304 | |
|---|
| 2305 | implicit none |
|---|
| 2306 | |
|---|
| 2307 | include 'nlte_paramdef.h' |
|---|
| 2308 | include 'nlte_commons.h' |
|---|
| 2309 | |
|---|
| 2310 | |
|---|
| 2311 | c arguments |
|---|
| 2312 | integer ihist |
|---|
| 2313 | |
|---|
| 2314 | c local variables |
|---|
| 2315 | integer j, r, mm |
|---|
| 2316 | real*8 xx |
|---|
| 2317 | |
|---|
| 2318 | c *************** |
|---|
| 2319 | |
|---|
| 2320 | nbox = nbox_stored(ihist) |
|---|
| 2321 | do j=1,mm_stored(ihist) |
|---|
| 2322 | thist(j) = thist_stored(ihist,j) |
|---|
| 2323 | do r=1,nbox_stored(ihist) |
|---|
| 2324 | no(r) = no_stored(ihist,r) |
|---|
| 2325 | sk1(j,r) = sk1_stored(ihist,j,r) |
|---|
| 2326 | xls1(j,r) = xls1_stored(ihist,j,r) |
|---|
| 2327 | xld1(j,r) = xld1_stored(ihist,j,r) |
|---|
| 2328 | enddo |
|---|
| 2329 | enddo |
|---|
| 2330 | |
|---|
| 2331 | |
|---|
| 2332 | return |
|---|
| 2333 | end |
|---|
| 2334 | |
|---|
| 2335 | |
|---|
| 2336 | c ******************************************************************* |
|---|
| 2337 | |
|---|
| 2338 | subroutine rhist_03 (ihist) |
|---|
| 2339 | |
|---|
| 2340 | c ******************************************************************* |
|---|
| 2341 | |
|---|
| 2342 | implicit none |
|---|
| 2343 | |
|---|
| 2344 | include 'nlte_paramdef.h' |
|---|
| 2345 | include 'nlte_commons.h' |
|---|
| 2346 | |
|---|
| 2347 | |
|---|
| 2348 | c arguments |
|---|
| 2349 | integer ihist |
|---|
| 2350 | |
|---|
| 2351 | c local variables |
|---|
| 2352 | integer j, r, mm |
|---|
| 2353 | real*8 xx |
|---|
| 2354 | |
|---|
| 2355 | c *************** |
|---|
| 2356 | |
|---|
| 2357 | open(unit=3,file=hisfile,status='old') |
|---|
| 2358 | |
|---|
| 2359 | read(3,*) |
|---|
| 2360 | read(3,*) |
|---|
| 2361 | read(3,*) mm_stored(ihist) |
|---|
| 2362 | read(3,*) |
|---|
| 2363 | read(3,*) nbox_stored(ihist) |
|---|
| 2364 | read(3,*) |
|---|
| 2365 | |
|---|
| 2366 | if ( nbox_stored(ihist) .gt. nbox_max ) then |
|---|
| 2367 | write (*,*) ' nbox too large in input file ', hisfile |
|---|
| 2368 | stop ' Check maximum number nbox_max in mz1d.par ' |
|---|
| 2369 | endif |
|---|
| 2370 | |
|---|
| 2371 | do j=1,mm_stored(ihist) |
|---|
| 2372 | read(3,*) thist_stored(ihist,j) |
|---|
| 2373 | do r=1,nbox_stored(ihist) |
|---|
| 2374 | read(3,*) no_stored(ihist,r), |
|---|
| 2375 | & sk1_stored(ihist,j,r), |
|---|
| 2376 | & xls1_stored(ihist,j,r), |
|---|
| 2377 | & xx, |
|---|
| 2378 | & xld1_stored(ihist,j,r) |
|---|
| 2379 | enddo |
|---|
| 2380 | |
|---|
| 2381 | enddo |
|---|
| 2382 | |
|---|
| 2383 | close(unit=3) |
|---|
| 2384 | |
|---|
| 2385 | |
|---|
| 2386 | return |
|---|
| 2387 | end |
|---|
