| 1 | c*********************************************************************** |
|---|
| 2 | c File with all subroutines required by mztf |
|---|
| 3 | c Subroutines previously included in mztfsub_overlap.F |
|---|
| 4 | c |
|---|
| 5 | c jan 98 malv basado en mztfsub_solar |
|---|
| 6 | c jul 2011 malv+fgg adapted to LMD-MGCM |
|---|
| 7 | c |
|---|
| 8 | c contiene: |
|---|
| 9 | c initial |
|---|
| 10 | c intershape |
|---|
| 11 | c interstrength |
|---|
| 12 | c intz |
|---|
| 13 | c rhist |
|---|
| 14 | c we |
|---|
| 15 | c simrul |
|---|
| 16 | c fi |
|---|
| 17 | c f |
|---|
| 18 | c findw |
|---|
| 19 | c voigtf |
|---|
| 20 | c*********************************************************************** |
|---|
| 21 | |
|---|
| 22 | c **************************************************************** |
|---|
| 23 | subroutine initial |
|---|
| 24 | |
|---|
| 25 | c ma & crs !evita troubles 16-july-96 |
|---|
| 26 | c **************************************************************** |
|---|
| 27 | |
|---|
| 28 | implicit none |
|---|
| 29 | |
|---|
| 30 | include 'nlte_paramdef.h' |
|---|
| 31 | include 'nlte_commons.h' |
|---|
| 32 | |
|---|
| 33 | c local variables |
|---|
| 34 | integer i |
|---|
| 35 | |
|---|
| 36 | c *************** |
|---|
| 37 | |
|---|
| 38 | eqw = 0.0d00 |
|---|
| 39 | aa = 0.0d00 |
|---|
| 40 | bb = 0.0d00 |
|---|
| 41 | cc = 0.0d00 |
|---|
| 42 | dd = 0.0d00 |
|---|
| 43 | |
|---|
| 44 | do i=1,nbox |
|---|
| 45 | ua(i) = 0.0d0 |
|---|
| 46 | ccbox(i) = 0.0d0 |
|---|
| 47 | ddbox(i) = 0.0d0 |
|---|
| 48 | end do |
|---|
| 49 | |
|---|
| 50 | return |
|---|
| 51 | end |
|---|
| 52 | |
|---|
| 53 | c ********************************************************************** |
|---|
| 54 | subroutine intershape(alsx,alnx,adx,xtemp) |
|---|
| 55 | c interpolates the line shape parameters at a temperature xtemp from |
|---|
| 56 | c input histogram data. |
|---|
| 57 | c ********************************************************************** |
|---|
| 58 | |
|---|
| 59 | implicit none |
|---|
| 60 | |
|---|
| 61 | include 'nlte_paramdef.h' |
|---|
| 62 | include 'nlte_commons.h' |
|---|
| 63 | |
|---|
| 64 | c arguments |
|---|
| 65 | real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max), |
|---|
| 66 | & xtemp(nbox_max) |
|---|
| 67 | |
|---|
| 68 | c local variables |
|---|
| 69 | integer i, k |
|---|
| 70 | |
|---|
| 71 | c *********** |
|---|
| 72 | |
|---|
| 73 | ! write (*,*) 'intershape xtemp =', xtemp |
|---|
| 74 | |
|---|
| 75 | do 1, k=1,nbox |
|---|
| 76 | if (xtemp(k).gt.tmax) then |
|---|
| 77 | write (*,*) ' WARNING ! Tpath,tmax= ',xtemp(k),tmax |
|---|
| 78 | xtemp(k) = tmax |
|---|
| 79 | endif |
|---|
| 80 | if (xtemp(k).lt.tmin) then |
|---|
| 81 | write (*,*) ' WARNING ! Tpath,tmin= ',xtemp(k),tmin |
|---|
| 82 | xtemp(k) = tmin |
|---|
| 83 | endif |
|---|
| 84 | |
|---|
| 85 | i = 1 |
|---|
| 86 | do while (i.le.mm) |
|---|
| 87 | i = i + 1 |
|---|
| 88 | |
|---|
| 89 | if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles |
|---|
| 90 | alsx(k)=xls1(i,k) !16-july-1996 |
|---|
| 91 | alnx(k)=xln1(i,k) |
|---|
| 92 | adx(k)=xld1(i,k) |
|---|
| 93 | goto 1 |
|---|
| 94 | elseif ( thist(i) .le. xtemp(k) ) then |
|---|
| 95 | alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) + |
|---|
| 96 | @ xls1(i-1,k)*(xtemp(k)-thist(i)) )) / |
|---|
| 97 | $ (thist(i-1)-thist(i)) |
|---|
| 98 | alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) + |
|---|
| 99 | @ xln1(i-1,k)*(xtemp(k)-thist(i)) )) / |
|---|
| 100 | $ (thist(i-1)-thist(i)) |
|---|
| 101 | adx(k) = (( xld1(i,k)*(thist(i-1)-xtemp(k)) + |
|---|
| 102 | @ xld1(i-1,k)*(xtemp(k)-thist(i)) )) / |
|---|
| 103 | $ (thist(i-1)-thist(i)) |
|---|
| 104 | goto 1 |
|---|
| 105 | end if |
|---|
| 106 | end do |
|---|
| 107 | write (*,*) |
|---|
| 108 | @ ' error in xtemp(k). it should be between tmin and tmax' |
|---|
| 109 | 1 continue |
|---|
| 110 | |
|---|
| 111 | return |
|---|
| 112 | end |
|---|
| 113 | c ********************************************************************** |
|---|
| 114 | subroutine interstrength (stx, ts, sx, xtemp) |
|---|
| 115 | c interpolates the line strength at a temperature xtemp from |
|---|
| 116 | c input histogram data. |
|---|
| 117 | c ********************************************************************** |
|---|
| 118 | |
|---|
| 119 | implicit none |
|---|
| 120 | |
|---|
| 121 | include 'nlte_paramdef.h' |
|---|
| 122 | include 'nlte_commons.h' |
|---|
| 123 | |
|---|
| 124 | c arguments |
|---|
| 125 | real*8 stx ! output, total band strength |
|---|
| 126 | real*8 ts ! input, temp for stx |
|---|
| 127 | real*8 sx(nbox_max) ! output, strength for each box |
|---|
| 128 | real*8 xtemp(nbox_max) ! input, temp for sx |
|---|
| 129 | |
|---|
| 130 | c local variables |
|---|
| 131 | integer i, k |
|---|
| 132 | |
|---|
| 133 | c *********** |
|---|
| 134 | |
|---|
| 135 | do 1, k=1,nbox |
|---|
| 136 | ! if(xtemp(k).lt.ts)then |
|---|
| 137 | ! write(*,*)'***********************' |
|---|
| 138 | ! write(*,*)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k |
|---|
| 139 | ! write(*,*)'***********************' |
|---|
| 140 | ! endif |
|---|
| 141 | if (xtemp(k).gt.tmax) xtemp(k) = tmax |
|---|
| 142 | if (xtemp(k).lt.tmin) xtemp(k) = tmin |
|---|
| 143 | i = 1 |
|---|
| 144 | do while (i.le.mm-1) |
|---|
| 145 | i = i + 1 |
|---|
| 146 | ! write(*,*)'mztfsub_overlap/136',i,xtemp(k),thist(i) |
|---|
| 147 | if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then |
|---|
| 148 | sx(k) = sk1(i,k) |
|---|
| 149 | ! write(*,*)'mztfsub_overlap/139',sx(k),k,i |
|---|
| 150 | goto 1 |
|---|
| 151 | elseif ( thist(i) .le. xtemp(k) ) then |
|---|
| 152 | sx(k) = ( sk1(i,k)*(thist(i-1)-xtemp(k)) + sk1(i-1,k)* |
|---|
| 153 | @ (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i)) |
|---|
| 154 | ! write(*,*)'mztfsub_overlap/144',sx(k),k,i |
|---|
| 155 | goto 1 |
|---|
| 156 | end if |
|---|
| 157 | end do |
|---|
| 158 | write (*,*) ' error in xtemp(kr) =', xtemp(k), |
|---|
| 159 | @ '. it should be between ' |
|---|
| 160 | write (*,*) ' tmin =',tmin, ' and tmax =',tmax |
|---|
| 161 | stop |
|---|
| 162 | 1 continue |
|---|
| 163 | |
|---|
| 164 | stx = 0.d0 |
|---|
| 165 | if (ts.gt.tmax) ts = dble( tmax ) |
|---|
| 166 | if (ts.lt.tmin) ts = dble( tmin ) |
|---|
| 167 | i = 1 |
|---|
| 168 | do while (i.le.mm-1) |
|---|
| 169 | i = i + 1 |
|---|
| 170 | ! write(*,*)'mztfsub_overlap/160',i,ts,thist(i) |
|---|
| 171 | if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then |
|---|
| 172 | do k=1,nbox |
|---|
| 173 | stx = stx + no(k) * sk1(i,k) |
|---|
| 174 | ! write(*,*)'mztfsub_overlap/164',stx |
|---|
| 175 | end do |
|---|
| 176 | return |
|---|
| 177 | elseif ( thist(i) .le. ts ) then |
|---|
| 178 | do k=1,nbox |
|---|
| 179 | stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) + |
|---|
| 180 | @ sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i)) |
|---|
| 181 | ! write(*,*)'mztfsub_overlap/171',stx |
|---|
| 182 | end do |
|---|
| 183 | ! stop |
|---|
| 184 | return |
|---|
| 185 | end if |
|---|
| 186 | end do |
|---|
| 187 | |
|---|
| 188 | return |
|---|
| 189 | end |
|---|
| 190 | |
|---|
| 191 | |
|---|
| 192 | c ********************************************************************** |
|---|
| 193 | subroutine intz(h,aco2,ap,amr,at, con) |
|---|
| 194 | c return interp. concentration, pressure,mixing ratio and temperature |
|---|
| 195 | c for a input height h |
|---|
| 196 | c ********************************************************************** |
|---|
| 197 | |
|---|
| 198 | implicit none |
|---|
| 199 | include 'nlte_paramdef.h' |
|---|
| 200 | include 'nlte_commons.h' |
|---|
| 201 | |
|---|
| 202 | c arguments |
|---|
| 203 | real h ! i |
|---|
| 204 | real*8 con(nzy) ! i |
|---|
| 205 | real*8 aco2, ap, at, amr ! o |
|---|
| 206 | |
|---|
| 207 | c local variables |
|---|
| 208 | integer k |
|---|
| 209 | |
|---|
| 210 | c ************ |
|---|
| 211 | |
|---|
| 212 | if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then |
|---|
| 213 | write (*,*) ' zp= ',h,' zy(1)= ',zy(1) |
|---|
| 214 | stop'from intz: error in interpolation, z < minimum height' |
|---|
| 215 | elseif (h.gt.zy(nzy)) then |
|---|
| 216 | write (*,*) ' zp= ',h,' zy(nzy)= ',zy(nzy) |
|---|
| 217 | stop'from intz: error in interpolation, z > maximum height' |
|---|
| 218 | end if |
|---|
| 219 | |
|---|
| 220 | if (h.eq.zy(nzy)) then |
|---|
| 221 | ap = dble( py(nzy) ) |
|---|
| 222 | aco2= con(nzy) |
|---|
| 223 | at = dble( ty(nzy) ) |
|---|
| 224 | amr = dble( mr(nzy) ) |
|---|
| 225 | return |
|---|
| 226 | end if |
|---|
| 227 | |
|---|
| 228 | do k=1,nzy-1 |
|---|
| 229 | if( abs( h-zy(k) ).le.( 1.e-5 ) ) then |
|---|
| 230 | ap = dble( py(k) ) |
|---|
| 231 | aco2= con(k) |
|---|
| 232 | at = dble( ty(k) ) |
|---|
| 233 | amr = dble( mr(k) ) |
|---|
| 234 | return |
|---|
| 235 | elseif(h.gt.zy(k).and.h.lt.zy(k+1))then |
|---|
| 236 | ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * |
|---|
| 237 | @ (h-zy(k)) / (zy(k+1)-zy(k)) ) ) |
|---|
| 238 | aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) * |
|---|
| 239 | @ (h-zy(k)) / (zy(k+1)-zy(k)) ) |
|---|
| 240 | at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/ |
|---|
| 241 | @ (zy(k+1)-zy(k)) ) |
|---|
| 242 | amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/ |
|---|
| 243 | @ (zy(k+1)-zy(k)) ) |
|---|
| 244 | return |
|---|
| 245 | end if |
|---|
| 246 | end do |
|---|
| 247 | |
|---|
| 248 | return |
|---|
| 249 | end |
|---|
| 250 | |
|---|
| 251 | |
|---|
| 252 | |
|---|
| 253 | c ********************************************************************** |
|---|
| 254 | real*8 function iaa_we(ig,me,pe,plaux,idummy,nt_local,p_local, |
|---|
| 255 | $ Desp,wsL) |
|---|
| 256 | c icls=5 -->para mztf |
|---|
| 257 | c icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2) |
|---|
| 258 | c calculates an approximate equivalent width for an error estimate. |
|---|
| 259 | c |
|---|
| 260 | c ioverlap = 0 ....... no correction for overlaping |
|---|
| 261 | c 1 ....... "lisat" first correction (see overlap_box. |
|---|
| 262 | c 2 ....... " " " plus "supersaturation" |
|---|
| 263 | |
|---|
| 264 | c idummy=0 do nothing |
|---|
| 265 | c 1 write out some values for diagnostics |
|---|
| 266 | c 2 correct the Strong Lorentz behaviour for SZA>90 |
|---|
| 267 | c 3 casos 1 & 2 |
|---|
| 268 | |
|---|
| 269 | c malv nov-98 add overlaping's corrections |
|---|
| 270 | c ********************************************************************** |
|---|
| 271 | |
|---|
| 272 | implicit none |
|---|
| 273 | |
|---|
| 274 | include 'nlte_paramdef.h' |
|---|
| 275 | include 'nlte_commons.h' |
|---|
| 276 | |
|---|
| 277 | c arguments |
|---|
| 278 | integer ig ! ADDED FOR TRACEBACK |
|---|
| 279 | real*8 me ! I. path's absorber amount |
|---|
| 280 | real*8 pe ! I. path's presion total |
|---|
| 281 | real*8 plaux ! I. path's partial pressure of CO2 |
|---|
| 282 | real*8 nt_local ! I. needed for strong limit of Lorentz profil |
|---|
| 283 | real*8 p_local ! I. " " " |
|---|
| 284 | integer idummy ! I. indica varias opciones |
|---|
| 285 | real*8 wsL ! O. need this for strong Lorentz correction |
|---|
| 286 | real*8 Desp ! I. need this for strong Lorentz correction |
|---|
| 287 | |
|---|
| 288 | c local variables |
|---|
| 289 | integer i |
|---|
| 290 | real*8 y,x,wl,wd |
|---|
| 291 | real*8 cn(0:7),dn(0:7) |
|---|
| 292 | real*8 pi, xx |
|---|
| 293 | real*8 f_sat_box |
|---|
| 294 | real*8 dv_sat_box, dv_corte_box |
|---|
| 295 | real*8 area_core_box, area_wing_box |
|---|
| 296 | real*8 wlgood , parentesis , xlor |
|---|
| 297 | real*8 wsl_grad |
|---|
| 298 | |
|---|
| 299 | |
|---|
| 300 | c data blocks |
|---|
| 301 | data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2, |
|---|
| 302 | @ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4, |
|---|
| 303 | @ 3.42209457833d-5,-1.28380804108d-6/ |
|---|
| 304 | data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1, |
|---|
| 305 | @ 8.21896973657d-1,-2.5222672453,6.1007027481, |
|---|
| 306 | @ -8.51001627836,4.6535116765/ |
|---|
| 307 | |
|---|
| 308 | c *********** |
|---|
| 309 | |
|---|
| 310 | c equivalent width of atmospheric line. |
|---|
| 311 | |
|---|
| 312 | pi = acos(-1.d0) |
|---|
| 313 | |
|---|
| 314 | if ( idummy.gt.9 ) |
|---|
| 315 | @ write (*,*) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), plaux |
|---|
| 316 | |
|---|
| 317 | y=ka(kr)*me |
|---|
| 318 | ! x=y/(2.0*pi*(alsa(kr)*pl+alna(kr)*(pe-pl))) |
|---|
| 319 | x=y/(2.0d0*pi* alsa(kr)*plaux) !+alna(kr)*(pe-pl))) |
|---|
| 320 | |
|---|
| 321 | ! Strong limit of Lorentz profile: WL = 2 SQRT( S * m * alsa*pl ) |
|---|
| 322 | ! Para anular esto, comentar las siguientes 5 lineas |
|---|
| 323 | ! if ( x .gt. 1.e6 ) then |
|---|
| 324 | ! wl = 2.0*sqrt( y * alsa(kr)*pl ) |
|---|
| 325 | ! else |
|---|
| 326 | ! wl=y/sqrt(1.0d0+pi*x/2.0d0) |
|---|
| 327 | ! endif |
|---|
| 328 | |
|---|
| 329 | wl=y/sqrt(1.0d0+pi*x/2.0d0) |
|---|
| 330 | |
|---|
| 331 | if (wl .le. 0.d0) then |
|---|
| 332 | write(*,*)'mztfsub_overlap/496',ig,y,ka(kr),me,kr |
|---|
| 333 | stop'WE/Lorentz EQW zero or negative!/498' !,ig |
|---|
| 334 | endif |
|---|
| 335 | |
|---|
| 336 | if ( idummy.gt.9 ) |
|---|
| 337 | @ write (*,*) ' y, x =', y, x |
|---|
| 338 | |
|---|
| 339 | xlor = x |
|---|
| 340 | if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then |
|---|
| 341 | ! en caso que estemos en el regimen |
|---|
| 342 | ! Strong Lorentz y la presion local |
|---|
| 343 | ! vaya disminuyendo, corregimos la EQW |
|---|
| 344 | ! con un gradiente analitico (notebook) |
|---|
| 345 | wsL = 2.0*sqrt( y * alsa(kr)*plaux ) |
|---|
| 346 | wsl_grad = - 2.d0 * ka(kr)*alsa(kr) * nt_local*p_local / wsL |
|---|
| 347 | wlgood = w_strongLor_prev(kr) + wsl_grad * Desp |
|---|
| 348 | if (idummy.eq.12) |
|---|
| 349 | @ write (*,*) ' W(wrong), W_SL, W_SL prev, W_SL corrected=', |
|---|
| 350 | @ wl, wsL, w_strongLor_prev(kr), wlgood |
|---|
| 351 | wl = wlgood |
|---|
| 352 | endif |
|---|
| 353 | ! wsL = wl pero esto no lo hacemos todavia, porque necesitamos |
|---|
| 354 | ! el valor que ahora mismo tiene wsL para corregir la |
|---|
| 355 | ! expresion R&W below |
|---|
| 356 | |
|---|
| 357 | ! write (*,*) 'WE arguments me,pe,pl =', me,pe,pl |
|---|
| 358 | ! write (*,*) 'WE/ wl,ka(kr),alsa(kr) =', |
|---|
| 359 | ! @ wl, ka(kr),alsa(kr) |
|---|
| 360 | |
|---|
| 361 | |
|---|
| 362 | !>>>>>>> |
|---|
| 363 | 500 format (a,i3,3(2x,1pe15.8)) |
|---|
| 364 | 600 format (a,2(2x,1pe16.9)) |
|---|
| 365 | 700 format (a,3(1x,1pe16.9)) |
|---|
| 366 | ! if (kr.eq.8 .or. kr.eq.13) then |
|---|
| 367 | ! write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl |
|---|
| 368 | ! write (*,700) ' /aln,als,d_x=', alna(kr),alsa(kr), |
|---|
| 369 | ! @ 2.0*pi*( alsa(kr)*pl + alna(kr)*(pe-pl) ) |
|---|
| 370 | ! write (*,600) ' /alsa*p_CO2, alna*p_n2 :', |
|---|
| 371 | ! @ alsa(kr)*pl, alna(kr)*(pe-pl) |
|---|
| 372 | ! write (*,600) ' a*p, S =', |
|---|
| 373 | ! @ alsa(kr)*pl + alna(kr)*(pe-pl) , ka(kr) |
|---|
| 374 | ! write (*,600) ' /S*m, x =', y, x |
|---|
| 375 | ! write (*,600) ' /aprox, WL =', |
|---|
| 376 | ! @ 2.*sqrt( y*( alsa(kr)*pl+alna(kr)*(pe-pl) ) ), WL |
|---|
| 377 | ! endif |
|---|
| 378 | ! |
|---|
| 379 | ! corrections to lorentz eqw due to overlaping and super-saturation |
|---|
| 380 | ! |
|---|
| 381 | |
|---|
| 382 | i_supersat = 0 |
|---|
| 383 | |
|---|
| 384 | if ( icls.eq.5 .and. ioverlap.gt.0 ) then |
|---|
| 385 | ! for the moment, only consider overlaping for mztf.f, not fot.f |
|---|
| 386 | |
|---|
| 387 | ! definition of saturation in the lisat model |
|---|
| 388 | ! |
|---|
| 389 | asat_box = 0.99d0 |
|---|
| 390 | f_sat_box = 2.d0 * x |
|---|
| 391 | xx = f_sat_box / log( 1./(1-asat_box) ) |
|---|
| 392 | if ( xx .lt. 1.0d0 ) then |
|---|
| 393 | dv_sat_box = 0.0d0 |
|---|
| 394 | asat_box = 1.0d0 - exp( - f_sat_box ) |
|---|
| 395 | else |
|---|
| 396 | dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 ) |
|---|
| 397 | ! approximation: only use of alsa in mars and venus |
|---|
| 398 | endif |
|---|
| 399 | |
|---|
| 400 | ! area of saturated line |
|---|
| 401 | ! |
|---|
| 402 | area_core_box = 2.d0 * dv_sat_box * asat_box |
|---|
| 403 | area_wing_box = 0.5d0 * ( wl - area_core_box ) |
|---|
| 404 | dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box |
|---|
| 405 | |
|---|
| 406 | ! super-saturation or simple overlaping? |
|---|
| 407 | ! |
|---|
| 408 | ! i_supersat = 0 |
|---|
| 409 | xx = dv_sat_box - ( 0.5d0 * dist(kr) ) |
|---|
| 410 | if ( xx .ge. 0.0 ! definition of supersaturation |
|---|
| 411 | @ .and. dv_sat_box .gt. 0.0 ! definition of saturation |
|---|
| 412 | @ .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line |
|---|
| 413 | @ ! and not too far apart |
|---|
| 414 | @ then |
|---|
| 415 | |
|---|
| 416 | i_supersat = 1 |
|---|
| 417 | |
|---|
| 418 | else |
|---|
| 419 | ! no super-saturation, then use "lisat + first correction", i.e., |
|---|
| 420 | ! correct for line products |
|---|
| 421 | ! |
|---|
| 422 | |
|---|
| 423 | wl = wl |
|---|
| 424 | |
|---|
| 425 | endif |
|---|
| 426 | |
|---|
| 427 | end if ! end of overlaping loop |
|---|
| 428 | |
|---|
| 429 | if (icls.eq.2) then |
|---|
| 430 | iaa_we = wl |
|---|
| 431 | return |
|---|
| 432 | endif |
|---|
| 433 | |
|---|
| 434 | cc doppler limit: |
|---|
| 435 | if ( idummy.gt.9 ) |
|---|
| 436 | @ write (*,*) ' S*m, alf_dop =', y, alda(kr)*sqrt(pi) |
|---|
| 437 | |
|---|
| 438 | x = y / (alda(kr)*sqrt(pi)) |
|---|
| 439 | if ( x.lt.1.e-10 ) then ! to avoid underflow |
|---|
| 440 | wd = y |
|---|
| 441 | else |
|---|
| 442 | wd=alda(kr)*sqrt(4.0*pi*x**2*(1.0+log(1.0+x))/(4.0+pi*x**2)) |
|---|
| 443 | endif |
|---|
| 444 | if ( idummy.gt.9 ) |
|---|
| 445 | @ write (*,*) ' wd =', wd |
|---|
| 446 | |
|---|
| 447 | cc doppler weak limit |
|---|
| 448 | c wd = ka(kr) * me |
|---|
| 449 | |
|---|
| 450 | cc good doppler |
|---|
| 451 | if(icls.eq.5) then !para mztf |
|---|
| 452 | !write (*,*) 'para mztf, icls=',icls |
|---|
| 453 | if (x.lt.5.) then |
|---|
| 454 | wd = 0.d0 |
|---|
| 455 | do i=0,7 |
|---|
| 456 | wd = wd + cn(i) * x**i |
|---|
| 457 | end do |
|---|
| 458 | wd = alda(kr) * x * sqrt(pi) * wd |
|---|
| 459 | elseif (x.gt.5.) then |
|---|
| 460 | wd = 0.d0 |
|---|
| 461 | do i=0,7 |
|---|
| 462 | wd = wd + dn(i) / (log(x))**i |
|---|
| 463 | end do |
|---|
| 464 | wd = alda(kr) * sqrt(log(x)) * wd |
|---|
| 465 | else |
|---|
| 466 | stop ' x should not be less than zero' |
|---|
| 467 | end if |
|---|
| 468 | end if |
|---|
| 469 | |
|---|
| 470 | |
|---|
| 471 | if ( i_supersat .eq. 0 ) then |
|---|
| 472 | |
|---|
| 473 | parentesis = wl**2+wd**2-(wd*wl/y)**2 |
|---|
| 474 | ! changed +(wd*wl/y)**2 to -...14-3-84 |
|---|
| 475 | |
|---|
| 476 | if ( parentesis .lt. 0.0 ) then |
|---|
| 477 | if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then |
|---|
| 478 | parentesis = wl**2+wd**2-(wd*wsL/y)**2 |
|---|
| 479 | ! este cambio puede ser necesario cuando se hace |
|---|
| 480 | ! correccion Strong Lor, para evitar valores |
|---|
| 481 | ! negativos del parentesis en sqrt( ) |
|---|
| 482 | else |
|---|
| 483 | stop ' WE/ Error en las EQW wl,wl,y ' |
|---|
| 484 | endif |
|---|
| 485 | endif |
|---|
| 486 | |
|---|
| 487 | iaa_we = sqrt( parentesis ) |
|---|
| 488 | ! write (*,*) ' from iaa_we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd) |
|---|
| 489 | ! write (*,*) ' from iaa_we: we', iaa_we |
|---|
| 490 | |
|---|
| 491 | else |
|---|
| 492 | |
|---|
| 493 | iaa_we = wl |
|---|
| 494 | ! if there is supersaturation we can ignore wd completely; |
|---|
| 495 | ! mztf.f will compute the eqw of the whole box afterwards |
|---|
| 496 | |
|---|
| 497 | endif |
|---|
| 498 | |
|---|
| 499 | if (icls.eq.3) iaa_we = wd |
|---|
| 500 | |
|---|
| 501 | if ( idummy.gt.9 ) |
|---|
| 502 | @ write (*,*) ' wl,wd,w =', wl,wd,iaa_we |
|---|
| 503 | |
|---|
| 504 | wsL = wl |
|---|
| 505 | |
|---|
| 506 | return |
|---|
| 507 | end |
|---|
| 508 | |
|---|
| 509 | |
|---|
| 510 | c ********************************************************************** |
|---|
| 511 | real*8 function simrul(a,b,fsim,c,acc) |
|---|
| 512 | c adaptively integrates fsim from a to b, within the criterion acc. |
|---|
| 513 | c ********************************************************************** |
|---|
| 514 | |
|---|
| 515 | implicit none |
|---|
| 516 | |
|---|
| 517 | real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim |
|---|
| 518 | real*8 s1(70),s2(70),s3(70) |
|---|
| 519 | real*8 c1, c2 |
|---|
| 520 | integer*4 m,n,j |
|---|
| 521 | |
|---|
| 522 | res=0. |
|---|
| 523 | c=0. |
|---|
| 524 | m=0 |
|---|
| 525 | n=0 |
|---|
| 526 | j=30 |
|---|
| 527 | g0=fsim(a) |
|---|
| 528 | g2=fsim((a+b)/2.) |
|---|
| 529 | g4=fsim(b) |
|---|
| 530 | a0=(b-a)*(g0+4.0*g2+g4)/2.0 |
|---|
| 531 | 1 d=2.0**n |
|---|
| 532 | h=(b-a)/(4.0*d) |
|---|
| 533 | x=a+(4.0*m+1.0)*h |
|---|
| 534 | g1=fsim(x) |
|---|
| 535 | g3=fsim(x+2.0*h) |
|---|
| 536 | a1=h*(g0+4.0*g1+g2) |
|---|
| 537 | a2=h*(g2+4.0*g3+g4) |
|---|
| 538 | if ( abs(a1+a2-a0).gt.(acc/d)) goto 2 |
|---|
| 539 | res=res+(16.0*(a1+a2)-a0)/45.0 |
|---|
| 540 | m=m+1 |
|---|
| 541 | c=a+m*(b-a)/d |
|---|
| 542 | 6 if (m.eq.(2*(m/2))) goto 4 |
|---|
| 543 | if ((m.ne.1).or.(n.ne.0)) goto 5 |
|---|
| 544 | 8 simrul=res |
|---|
| 545 | return |
|---|
| 546 | 2 m=2*m |
|---|
| 547 | n=n+1 |
|---|
| 548 | if (n.gt.j) goto 3 |
|---|
| 549 | a0=a1 |
|---|
| 550 | s1(n)=a2 |
|---|
| 551 | s2(n)=g3 |
|---|
| 552 | s3(n)=g4 |
|---|
| 553 | g4=g2 |
|---|
| 554 | g2=g1 |
|---|
| 555 | goto 1 |
|---|
| 556 | 3 c1=c-(b-a)/d |
|---|
| 557 | c2=c+(b-a)/d |
|---|
| 558 | write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) |
|---|
| 559 | write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) |
|---|
| 560 | 7 format(2x,'17hsimrule fails at ',/,3e15.6,/,3e15.6) |
|---|
| 561 | goto 8 |
|---|
| 562 | 5 a0=s1(n) |
|---|
| 563 | g0=g4 |
|---|
| 564 | g2=s2(n) |
|---|
| 565 | g4=s3(n) |
|---|
| 566 | goto 1 |
|---|
| 567 | 4 m=m/2 |
|---|
| 568 | n=n-1 |
|---|
| 569 | goto 6 |
|---|
| 570 | end |
|---|
| 571 | |
|---|
| 572 | c ********************************************************************** |
|---|
| 573 | subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL) |
|---|
| 574 | c this routine sets up accuracy criteria and calls simrule between limit |
|---|
| 575 | c that depend on the number of atmospheric and cell paths. it gives eqw. |
|---|
| 576 | |
|---|
| 577 | c Add correction for EQW in Strong Lorentz regime and SZA>90 |
|---|
| 578 | c ********************************************************************** |
|---|
| 579 | |
|---|
| 580 | implicit none |
|---|
| 581 | include 'nlte_paramdef.h' |
|---|
| 582 | include 'nlte_commons.h' |
|---|
| 583 | |
|---|
| 584 | c arguments |
|---|
| 585 | integer ig ! ADDED FOR TRACEBACK |
|---|
| 586 | integer iirw |
|---|
| 587 | integer idummy ! I. indica varias opciones |
|---|
| 588 | real*8 c1 ! I. needed for strong limit of Lorentz profil |
|---|
| 589 | real*8 p1 ! I. " " " |
|---|
| 590 | real*8 wsL ! O. need this for strong Lorentz correction |
|---|
| 591 | real*8 Desp ! I. need this for strong Lorentz correction |
|---|
| 592 | |
|---|
| 593 | c local variables |
|---|
| 594 | real*8 ept,eps,xa |
|---|
| 595 | real*8 acc, c |
|---|
| 596 | real*8 iaa_we |
|---|
| 597 | real*8 iaa_f, iaa_fi, simrul |
|---|
| 598 | |
|---|
| 599 | external iaa_f,iaa_fi |
|---|
| 600 | |
|---|
| 601 | c ********** *********** ********* |
|---|
| 602 | |
|---|
| 603 | if(icls.eq.5) then !para mztf |
|---|
| 604 | ! if(ig.eq.1682)write(*,*)'mztfsub_overlap/768',ua(kr),iirw |
|---|
| 605 | if (iirw.eq.2) then !iirw=icf=2 ==> we use the w&r formula |
|---|
| 606 | w = iaa_we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL ) |
|---|
| 607 | return |
|---|
| 608 | end if |
|---|
| 609 | ept=iaa_we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL) |
|---|
| 610 | else !para fot |
|---|
| 611 | if (iirw.eq.2) then ! icf=2 ==> we use the w&r formula |
|---|
| 612 | w = iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) |
|---|
| 613 | return |
|---|
| 614 | end if |
|---|
| 615 | ept=iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) |
|---|
| 616 | end if |
|---|
| 617 | |
|---|
| 618 | c the next block is a modification to avoid nul we. |
|---|
| 619 | c this situation appears for weak lines and low path temperature, but |
|---|
| 620 | c there is not any loss of accuracy. first july 1986 |
|---|
| 621 | if (ept.eq.0.) then ! for weak lines sometimes we=0 |
|---|
| 622 | ept=1.0e-18 |
|---|
| 623 | write (*,*) 'ept =',ept |
|---|
| 624 | write (*,*) 'from we: we=0.0' |
|---|
| 625 | return |
|---|
| 626 | end if |
|---|
| 627 | |
|---|
| 628 | acc = 4.d0 |
|---|
| 629 | acc = 10.d0**(-acc) |
|---|
| 630 | |
|---|
| 631 | eps = acc * ept !accuracy 10-4 atmospheric eqw. |
|---|
| 632 | xa=0.5*ept/iaa_f(0.d0) !width of doppler shifted atmospheric line. |
|---|
| 633 | w = 2.0*( simrul(0.0d0,xa,iaa_f,c,eps) |
|---|
| 634 | . + simrul(0.1d0,1.0/xa,iaa_fi,c,eps) ) |
|---|
| 635 | !no shift. |
|---|
| 636 | |
|---|
| 637 | return |
|---|
| 638 | end |
|---|
| 639 | |
|---|
| 640 | |
|---|
| 641 | c ********************************************************************** |
|---|
| 642 | double precision function iaa_fi(y) |
|---|
| 643 | c returns the value of f(1/y) |
|---|
| 644 | c ********************************************************************** |
|---|
| 645 | |
|---|
| 646 | implicit none |
|---|
| 647 | real*8 iaa_f, y |
|---|
| 648 | |
|---|
| 649 | iaa_fi=iaa_f(1.0/y)/y**2 |
|---|
| 650 | return |
|---|
| 651 | end |
|---|
| 652 | |
|---|
| 653 | |
|---|
| 654 | c ********************************************************************** |
|---|
| 655 | double precision function iaa_f(nuaux) |
|---|
| 656 | c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof |
|---|
| 657 | c ********************************************************************** |
|---|
| 658 | |
|---|
| 659 | implicit none |
|---|
| 660 | include 'nlte_paramdef.h' |
|---|
| 661 | include 'nlte_commons.h' |
|---|
| 662 | |
|---|
| 663 | double precision tra,xa,ya,za,yy,nuaux |
|---|
| 664 | double precision voigtf |
|---|
| 665 | tra=1.0d0 |
|---|
| 666 | |
|---|
| 667 | yy=1.0d0/alda(kr) |
|---|
| 668 | xa=nuaux*yy |
|---|
| 669 | ya= ( alsa(kr)*pp + alna(kr)*(pt-pp) ) * yy |
|---|
| 670 | za=ka(kr)*yy |
|---|
| 671 | |
|---|
| 672 | if(icls.eq.5) then !para mztf |
|---|
| 673 | ! write (*,*) 'icls=',icls |
|---|
| 674 | tra=za*ua(kr)*voigtf(sngl(xa),sngl(ya)) |
|---|
| 675 | else |
|---|
| 676 | tra=za*sl_ua*voigtf(sngl(xa),sngl(ya)) |
|---|
| 677 | end if |
|---|
| 678 | |
|---|
| 679 | if (tra.gt.50.0) then |
|---|
| 680 | tra=1.0 !2.0e-22 overflow cut-off. |
|---|
| 681 | else if (tra.gt.1.0e-4) then |
|---|
| 682 | tra=1.0-exp(-tra) |
|---|
| 683 | end if |
|---|
| 684 | |
|---|
| 685 | iaa_f=tra |
|---|
| 686 | return |
|---|
| 687 | end |
|---|
| 688 | |
|---|
| 689 | c ********************************************************************** |
|---|
| 690 | double precision function voigtf(x1,y) |
|---|
| 691 | c computes voigt function for any value of x1 and any +ve value of y. |
|---|
| 692 | c where possible uses modified lorentz and modified doppler approximatio |
|---|
| 693 | c otherwise uses a rearranged rybicki routine. |
|---|
| 694 | c c(n) = exp(-(n/h)**2)/(pi*sqrt(pi)), with h = 2.5 . |
|---|
| 695 | c accurate to better than 1 in 10000. |
|---|
| 696 | c ********************************************************************** |
|---|
| 697 | |
|---|
| 698 | implicit none |
|---|
| 699 | |
|---|
| 700 | real x1, y |
|---|
| 701 | real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig |
|---|
| 702 | |
|---|
| 703 | real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c |
|---|
| 704 | integer*4 n |
|---|
| 705 | |
|---|
| 706 | dimension c(10) |
|---|
| 707 | complex xp,xpp,z |
|---|
| 708 | |
|---|
| 709 | data c(1)/0.15303405/ |
|---|
| 710 | data c(2)/0.94694928e-1/ |
|---|
| 711 | data c(3)/0.42549174e-1/ |
|---|
| 712 | data c(4)/0.13882935e-1/ |
|---|
| 713 | data c(5)/0.32892528e-2/ |
|---|
| 714 | data c(6)/0.56589906e-3/ |
|---|
| 715 | data c(7)/0.70697890e-4/ |
|---|
| 716 | data c(8)/0.64135678e-5/ |
|---|
| 717 | data c(9)/0.42249221e-6/ |
|---|
| 718 | data c(10)/0.20209868e-7/ |
|---|
| 719 | |
|---|
| 720 | x=abs(x1) |
|---|
| 721 | if (x.gt.7.2) goto 1 |
|---|
| 722 | if ((y+x*0.3).gt.5.4) goto 1 |
|---|
| 723 | if (y.gt.0.01) goto 3 |
|---|
| 724 | if (x.lt.2.1) goto 2 |
|---|
| 725 | goto 3 |
|---|
| 726 | c here uses modified lorentz approx. |
|---|
| 727 | 1 xx=x*x |
|---|
| 728 | xxyy=xx+y*y |
|---|
| 729 | b=xx/xxyy |
|---|
| 730 | voigtf=y*(1.+(2.*b-0.5+(0.75-(9.-12.*b)*b)/xxyy)/ |
|---|
| 731 | * xxyy)/(xxyy*3.141592654) |
|---|
| 732 | return |
|---|
| 733 | c here uses modified doppler approx. |
|---|
| 734 | 2 xx=x*x |
|---|
| 735 | voigtf=0.56418958*exp(-xx)*(1.-y*(1.-0.5*y)*(1.1289-xx*(1.1623+ |
|---|
| 736 | * xx*(0.080812+xx*(0.13854-xx*(0.033605-0.0073972*xx)))))) |
|---|
| 737 | return |
|---|
| 738 | c here uses a rearranged rybicki routine. |
|---|
| 739 | 3 xh=2.5*x |
|---|
| 740 | xhxh=xh*xh |
|---|
| 741 | yh=2.5*y |
|---|
| 742 | yhyh=yh*yh |
|---|
| 743 | f1=xhxh+yhyh |
|---|
| 744 | f2=f1-0.5*yhyh |
|---|
| 745 | if (y.lt.0.1) goto 20 |
|---|
| 746 | p=-y*7.8539816 !7.8539816=2.5*pi |
|---|
| 747 | q=x*7.8539816 |
|---|
| 748 | xpp=cmplx(p,q) |
|---|
| 749 | z=cexp(xpp) |
|---|
| 750 | d1=xh*aimag(z) |
|---|
| 751 | d2=-d1 |
|---|
| 752 | d3=yh*(1.-real(z)) |
|---|
| 753 | d4=-d3+2.*yh |
|---|
| 754 | voig=0.17958712*(d1+d3)/f1 |
|---|
| 755 | goto 30 |
|---|
| 756 | 20 p=x*7.8539816 |
|---|
| 757 | q=y*7.8539816 |
|---|
| 758 | xp=cmplx(p,q) |
|---|
| 759 | z=ccos(xp) |
|---|
| 760 | d1=xh*aimag(z) |
|---|
| 761 | d2=-d1 |
|---|
| 762 | d3=yh*(1.-real(z)) |
|---|
| 763 | d4=-d3+2.*yh |
|---|
| 764 | voig=0.56418958*exp(y*y-x*x)*cos(2.*x*y)+0.17958712*(d1+d3)/f1 |
|---|
| 765 | 30 xn=0. |
|---|
| 766 | do 55 n=1,10,2 |
|---|
| 767 | xn=xn+1. |
|---|
| 768 | xnxn=xn*xn |
|---|
| 769 | g1=xh-xn |
|---|
| 770 | g2=g1*(xh+xn) |
|---|
| 771 | g3=0.5*g2*g2 |
|---|
| 772 | voig=voig+c(n)*(d2*(g2+yhyh)+d4*(f1+xnxn))/ |
|---|
| 773 | & (g3+yhyh*(f2+xnxn)) |
|---|
| 774 | xn=xn+1. |
|---|
| 775 | xnxn=xn*xn |
|---|
| 776 | g1=xh-xn |
|---|
| 777 | g2=g1*(xh+xn) |
|---|
| 778 | g3=0.5*g2*g2 |
|---|
| 779 | voig=voig+c(n+1)*(d1*(g2+yhyh)+d3*(f1+xnxn))/ |
|---|
| 780 | @ (g3+yhyh*(f2+xnxn)) |
|---|
| 781 | 55 continue |
|---|
| 782 | voigtf=voig |
|---|
| 783 | return |
|---|
| 784 | end |
|---|
| 785 | |
|---|
| 786 | |
|---|
| 787 | |
|---|
| 788 | c ********************************************************************** |
|---|
| 789 | c elimin_mz1d.F (includes smooth_cf) |
|---|
| 790 | c ************************************************************************ |
|---|
| 791 | subroutine elimin_mz1d (c,vc, ilayer,nanaux,itblout, nwaux) |
|---|
| 792 | |
|---|
| 793 | c Eliminate anomalous negative numbers in c(nl,nl) according to "nanaux": |
|---|
| 794 | |
|---|
| 795 | c nanaux = 0 -> no eliminate |
|---|
| 796 | c @ -> eliminate all numbers with absol.value<abs(max(c(n,r)))/300. |
|---|
| 797 | c 2 -> eliminate all anomalous negative numbers in c(n,r). |
|---|
| 798 | c 3 -> eliminate all anomalous negative numbers far from the main |
|---|
| 799 | c diagonal. |
|---|
| 800 | c 8 -> eliminate all non-zero numbers outside the main diagonal, |
|---|
| 801 | c and the contibution of lower boundary. |
|---|
| 802 | c 9 -> eliminate all non-zero numbers outside the main diagonal. |
|---|
| 803 | c 4 -> hace un smoothing cuando la distancia de separacion entre |
|---|
| 804 | c el valor maximo y el minimo de cf > 50 capas. |
|---|
| 805 | c 5 -> elimina valores menores que 1.0d-19 |
|---|
| 806 | c 6 -> incluye los dos casos 4 y 5 |
|---|
| 807 | c 7 -> llama a lisa: smooth con width=nw & elimina mejorado |
|---|
| 808 | c 78-> incluye los dos casos 7 y 8 |
|---|
| 809 | c 79-> incluye los dos casos 7 y 9 |
|---|
| 810 | |
|---|
| 811 | c itblout (itableout in calling program) is the option for writing |
|---|
| 812 | c out or not the purged c(n,r) matrix: |
|---|
| 813 | c itblout = 0 -> no write |
|---|
| 814 | c = 1 -> write out in curtis***.out according to ilayer |
|---|
| 815 | |
|---|
| 816 | c ilayer is the index for the layer selected to write out the matrix: |
|---|
| 817 | c ilayer = 0 => matrix elements written out cover all the altitudes |
|---|
| 818 | c with 5 layers steps |
|---|
| 819 | c > 0 => " " " " are c(ilayer,*) |
|---|
| 820 | c NOTA: |
|---|
| 821 | c EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ) |
|---|
| 822 | c UTILIZANDO itableout=30 EN MZTUD |
|---|
| 823 | |
|---|
| 824 | c jul 2011 malv+fgg adapted to LMD-MGCM |
|---|
| 825 | c Sep-04 FGG+MALV Correct include and call parameters |
|---|
| 826 | c cristina 25-sept-1996 y 27-ene-1997 |
|---|
| 827 | c JAN 98 MALV Version for mz1d |
|---|
| 828 | c ************************************************************************ |
|---|
| 829 | |
|---|
| 830 | implicit none |
|---|
| 831 | |
|---|
| 832 | include 'nlte_paramdef.h' |
|---|
| 833 | include 'nlte_commons.h' |
|---|
| 834 | |
|---|
| 835 | integer nanaux,j,i,itblout,kk,k,ir,in |
|---|
| 836 | integer ilayer,jmin, jmax,np,nwaux,ntimes,ntimes2 |
|---|
| 837 | !* real*8 c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini |
|---|
| 838 | real*8 c(nl,nl), vc(nl), amax, cmax, cmin, mini |
|---|
| 839 | real*8 aux(nl), auxs(nl) |
|---|
| 840 | character layercode*3 |
|---|
| 841 | |
|---|
| 842 | ntimes=0 |
|---|
| 843 | ntimes2=0 |
|---|
| 844 | ! type *,'from elimin_mz4: nan, nw',nan,nw |
|---|
| 845 | |
|---|
| 846 | if (nanaux .eq. 0) goto 200 |
|---|
| 847 | |
|---|
| 848 | if(nanaux.eq.1)then |
|---|
| 849 | do i=1,nl |
|---|
| 850 | amax=1.0d-36 |
|---|
| 851 | do j=1,nl |
|---|
| 852 | if(abs(c(i,j)).gt.amax)amax=abs(c(i,j)) |
|---|
| 853 | end do |
|---|
| 854 | do j=1,nl |
|---|
| 855 | if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0 |
|---|
| 856 | end do |
|---|
| 857 | enddo |
|---|
| 858 | elseif(nanaux.eq.2)then |
|---|
| 859 | do i=1,nl |
|---|
| 860 | do j=1,nl |
|---|
| 861 | if( ( j.le.(i-2) .or. j.gt.(i+2) ).and. |
|---|
| 862 | @ ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0 |
|---|
| 863 | end do |
|---|
| 864 | enddo |
|---|
| 865 | elseif(nanaux.eq.3)then |
|---|
| 866 | do i=1,nl |
|---|
| 867 | do j=1,nl |
|---|
| 868 | if (abs(i-j).ge.10) c(i,j)=0.0d0 |
|---|
| 869 | end do |
|---|
| 870 | enddo |
|---|
| 871 | elseif(nanaux.eq.8)then |
|---|
| 872 | do i=1,nl |
|---|
| 873 | do j=1,i-1 |
|---|
| 874 | c(i,j)=0.0d0 |
|---|
| 875 | enddo |
|---|
| 876 | do j=i+1,nl |
|---|
| 877 | c(i,j)=0.0d0 |
|---|
| 878 | enddo |
|---|
| 879 | vc(i)= 0.d0 |
|---|
| 880 | enddo |
|---|
| 881 | elseif(nanaux.eq.9)then |
|---|
| 882 | do i=1,nl |
|---|
| 883 | do j=1,i-1 |
|---|
| 884 | c(i,j)=0.0d0 |
|---|
| 885 | enddo |
|---|
| 886 | do j=i+1,nl |
|---|
| 887 | c(i,j)=0.0d0 |
|---|
| 888 | enddo |
|---|
| 889 | enddo |
|---|
| 890 | ! elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then |
|---|
| 891 | ! call lisa(c, vc, nl, nw) |
|---|
| 892 | end if |
|---|
| 893 | if(nanaux.eq.78)then |
|---|
| 894 | do i=1,nl |
|---|
| 895 | do j=1,i-1 |
|---|
| 896 | c(i,j)=0.0d0 |
|---|
| 897 | enddo |
|---|
| 898 | do j=i+1,nl |
|---|
| 899 | c(i,j)=0.0d0 |
|---|
| 900 | enddo |
|---|
| 901 | vc(i)= 0.d0 |
|---|
| 902 | enddo |
|---|
| 903 | endif |
|---|
| 904 | if(nanaux.eq.79)then |
|---|
| 905 | do i=1,nl |
|---|
| 906 | do j=1,i-1 |
|---|
| 907 | c(i,j)=0.0d0 |
|---|
| 908 | enddo |
|---|
| 909 | do j=i+1,nl |
|---|
| 910 | c(i,j)=0.0d0 |
|---|
| 911 | enddo |
|---|
| 912 | enddo |
|---|
| 913 | endif |
|---|
| 914 | |
|---|
| 915 | if(nanaux.eq.5.or.nanaux.eq.6)then |
|---|
| 916 | do i=1,nl |
|---|
| 917 | mini = 1.0d-19 |
|---|
| 918 | do j=1,nl |
|---|
| 919 | if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then |
|---|
| 920 | ntimes2=ntimes2+1 |
|---|
| 921 | end if |
|---|
| 922 | if ( abs(c(i,j)).le.mini) c(i,j)=0.d0 |
|---|
| 923 | end do |
|---|
| 924 | enddo |
|---|
| 925 | end if |
|---|
| 926 | |
|---|
| 927 | if(nanaux.eq.4.or.nanaux.eq.6)then |
|---|
| 928 | do i=1,nl |
|---|
| 929 | do j=1,nl |
|---|
| 930 | aux(j)=c(i,j) |
|---|
| 931 | auxs(j)=c(i,j) |
|---|
| 932 | end do |
|---|
| 933 | !call maxdp_2(aux,nl,cmax,jmax) |
|---|
| 934 | cmax=maxval(aux) |
|---|
| 935 | jmax=maxloc(aux,dim=1) |
|---|
| 936 | if(abs(jmax-i).ge.50) then |
|---|
| 937 | call smooth_cf(aux,auxs,i,nl,3) |
|---|
| 938 | !!!call smooth_cf(aux,auxs,i,nl,5) |
|---|
| 939 | ntimes=ntimes+1 |
|---|
| 940 | end if |
|---|
| 941 | do j=1,nl |
|---|
| 942 | c(i,j)=auxs(j) |
|---|
| 943 | end do |
|---|
| 944 | end do |
|---|
| 945 | end if |
|---|
| 946 | |
|---|
| 947 | ! type *, 'elimin_mz4: c(n,r) procesed for elimination. ' |
|---|
| 948 | ! type *, ' ' |
|---|
| 949 | ! if(nan.eq.4.or.nan.eq.6) type *, ' call smoothing:',ntimes |
|---|
| 950 | ! if(nan.eq.5.or.nan.eq.6) type *, ' call elimina: ',ntimes2 |
|---|
| 951 | ! if(nan.eq.7) type *, ' from elimin: lisa w=',nw |
|---|
| 952 | ! type *, ' ' |
|---|
| 953 | |
|---|
| 954 | |
|---|
| 955 | 200 continue |
|---|
| 956 | |
|---|
| 957 | c writting out of c(n,r) in ascii file |
|---|
| 958 | |
|---|
| 959 | ! if(itblout.eq.1) then |
|---|
| 960 | |
|---|
| 961 | ! if (ilayer.eq.0) then |
|---|
| 962 | |
|---|
| 963 | ! open (unit=2, status='new', |
|---|
| 964 | ! @ file=dircurtis//'curtis_gnu.out', recl=1024) |
|---|
| 965 | ! write(2,'(a)') |
|---|
| 966 | ! @ ' curtis matrix: table with 1.e+7 * acf(n,r) ' |
|---|
| 967 | ! write(2,114) 'n,r', ( i, i=nl,1,-5 ) |
|---|
| 968 | ! do in=nl,1,-5 |
|---|
| 969 | ! write(2,*) |
|---|
| 970 | ! write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 ) |
|---|
| 971 | ! end do |
|---|
| 972 | ! close(2) |
|---|
| 973 | |
|---|
| 974 | |
|---|
| 975 | ! write (*,*) ' ' |
|---|
| 976 | ! write (*,*) ' curtis.out has been created. ' |
|---|
| 977 | ! write (*,*) ' ' |
|---|
| 978 | |
|---|
| 979 | ! else |
|---|
| 980 | |
|---|
| 981 | ! write (layercode,132) ilayer |
|---|
| 982 | ! open (2, status='new', |
|---|
| 983 | ! @ file=dircurtis//'curtis'//layercode//'.out') |
|---|
| 984 | ! write(2,'(a)') |
|---|
| 985 | ! @ ' curtis matrix: table with 1.e+7 * acf(n,r) ' |
|---|
| 986 | ! write(2,116) ' layer x c(',layercode, |
|---|
| 987 | ! @ ',x) c(x,', layercode,')' |
|---|
| 988 | ! do in=nl,1,-1 |
|---|
| 989 | ! if (c(ilayer,ilayer).ne.0.d0) then |
|---|
| 990 | ! write(2,117) in, c(ilayer,in), c(in,ilayer), |
|---|
| 991 | ! @ c(ilayer,in)/c(ilayer,ilayer), |
|---|
| 992 | ! @ c(in,ilayer)/c(ilayer,ilayer) |
|---|
| 993 | ! else |
|---|
| 994 | ! write(2,118) in, c(ilayer,in), c(in,ilayer) |
|---|
| 995 | ! end if |
|---|
| 996 | ! end do |
|---|
| 997 | ! close(2) |
|---|
| 998 | ! write (*,*) ' ' |
|---|
| 999 | ! write (*,*) dircurtis//'curtis'//layercode//'.out', |
|---|
| 1000 | ! @ ' has been created.' |
|---|
| 1001 | ! write (*,*) ' ' |
|---|
| 1002 | |
|---|
| 1003 | ! end if |
|---|
| 1004 | |
|---|
| 1005 | ! elseif(itblout.eq.0)then |
|---|
| 1006 | |
|---|
| 1007 | ! continue |
|---|
| 1008 | |
|---|
| 1009 | ! else |
|---|
| 1010 | |
|---|
| 1011 | ! write (*,*) ' error from elimin: ', |
|---|
| 1012 | ! @ ' itblout should be 1 or 0; itblout= ',itblout |
|---|
| 1013 | ! stop |
|---|
| 1014 | |
|---|
| 1015 | ! end if |
|---|
| 1016 | |
|---|
| 1017 | return |
|---|
| 1018 | |
|---|
| 1019 | 112 format(10x,10(i3,9x)) |
|---|
| 1020 | 113 format(1x,i3,2x,9(1pe9.2,2x)) |
|---|
| 1021 | |
|---|
| 1022 | 114 format(1x,a3, 11(8x,i3)) |
|---|
| 1023 | 115 format( 1x,i3, 2x, 11(1pe10.3)) |
|---|
| 1024 | 116 format( 1x,a17,a2,a18,a2,a1 ) |
|---|
| 1025 | 117 format( 3x,i3, 4(8x,1pe10.3) ) |
|---|
| 1026 | 118 format( 3x,i3, 2(8x,1pe10.3) ) |
|---|
| 1027 | 120 format( 1x,i3, 1x,i3, 2x, 11(1pe10.3)) |
|---|
| 1028 | |
|---|
| 1029 | 132 format(i3) |
|---|
| 1030 | |
|---|
| 1031 | ! cambio: los formatos 114, 115 , 117 y 118 |
|---|
| 1032 | ! cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3 |
|---|
| 1033 | ! y ahora en vez de 11 capas de 5 en 5, hay 28 |
|---|
| 1034 | ! |
|---|
| 1035 | end |
|---|
| 1036 | c************************************************************************** |
|---|
| 1037 | subroutine smooth_cf( c, cs, i, nl, w ) |
|---|
| 1038 | c hace un smoothing de c(i,*), de la contribucion de todas las capas |
|---|
| 1039 | c menos de la capa en cuestion, la i. |
|---|
| 1040 | c opcion w (width): el tamanho de la ventana del smoothing. |
|---|
| 1041 | c output values: cs |
|---|
| 1042 | c************************************************************************** |
|---|
| 1043 | |
|---|
| 1044 | implicit none |
|---|
| 1045 | |
|---|
| 1046 | integer j,np,i,nl,w |
|---|
| 1047 | real*8 c(nl), cs(nl) |
|---|
| 1048 | |
|---|
| 1049 | if(w.eq.0) then |
|---|
| 1050 | do j=1,nl |
|---|
| 1051 | cs(j)=c(j) |
|---|
| 1052 | end do |
|---|
| 1053 | |
|---|
| 1054 | elseif(w.eq.3) then |
|---|
| 1055 | |
|---|
| 1056 | ! write (*,*) 'smoothing w=3' |
|---|
| 1057 | do j=1,i-4 |
|---|
| 1058 | if(j.eq.1) then |
|---|
| 1059 | cs(j)=c(j) |
|---|
| 1060 | else |
|---|
| 1061 | cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1)) |
|---|
| 1062 | end if |
|---|
| 1063 | end do |
|---|
| 1064 | do j=i+4,nl-1 |
|---|
| 1065 | if(j.eq.nl) then |
|---|
| 1066 | cs(j)=c(j) |
|---|
| 1067 | else |
|---|
| 1068 | cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1)) |
|---|
| 1069 | end if |
|---|
| 1070 | end do |
|---|
| 1071 | elseif(w.eq.5) then |
|---|
| 1072 | |
|---|
| 1073 | ! type *,'smoothing w=5' |
|---|
| 1074 | do j=3,i-4 |
|---|
| 1075 | if(j.eq.1) then |
|---|
| 1076 | cs(j)=c(j) |
|---|
| 1077 | else |
|---|
| 1078 | cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2)) |
|---|
| 1079 | end if |
|---|
| 1080 | end do |
|---|
| 1081 | do j=i+4,nl-2 |
|---|
| 1082 | if(j.eq.nl) then |
|---|
| 1083 | cs(j)=c(j) |
|---|
| 1084 | else |
|---|
| 1085 | cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2)) |
|---|
| 1086 | end if |
|---|
| 1087 | end do |
|---|
| 1088 | end if |
|---|
| 1089 | return |
|---|
| 1090 | end |
|---|
| 1091 | |
|---|
| 1092 | |
|---|
| 1093 | |
|---|
| 1094 | c***************************************************************************** |
|---|
| 1095 | c suaviza |
|---|
| 1096 | c***************************************************************************** |
|---|
| 1097 | c |
|---|
| 1098 | subroutine suaviza ( x, n, ismooth, y ) |
|---|
| 1099 | c |
|---|
| 1100 | c x - input and return values |
|---|
| 1101 | c y - auxiliary vector |
|---|
| 1102 | c ismooth = 0 --> no smoothing is performed |
|---|
| 1103 | c ismooth = 1 --> weak smoothing (5 points, centred weighted) |
|---|
| 1104 | c ismooth = 2 --> normal smoothing (3 points, evenly weighted) |
|---|
| 1105 | c ismooth = 3 --> strong smoothing (5 points, evenly weighted) |
|---|
| 1106 | |
|---|
| 1107 | |
|---|
| 1108 | c malv august 1991 |
|---|
| 1109 | c***************************************************************************** |
|---|
| 1110 | |
|---|
| 1111 | implicit none |
|---|
| 1112 | |
|---|
| 1113 | integer n, imax, imin, i, ismooth |
|---|
| 1114 | real*8 x(n), y(n) |
|---|
| 1115 | c***************************************************************************** |
|---|
| 1116 | |
|---|
| 1117 | imin=1 |
|---|
| 1118 | imax=n |
|---|
| 1119 | |
|---|
| 1120 | if (ismooth.eq.0) then |
|---|
| 1121 | |
|---|
| 1122 | return |
|---|
| 1123 | |
|---|
| 1124 | elseif (ismooth.eq.1) then ! 5 points, with central weighting |
|---|
| 1125 | |
|---|
| 1126 | do i=imin,imax |
|---|
| 1127 | if(i.eq.imin)then |
|---|
| 1128 | y(i)=x(imin) |
|---|
| 1129 | elseif(i.eq.imax)then |
|---|
| 1130 | y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 |
|---|
| 1131 | elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then |
|---|
| 1132 | y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + |
|---|
| 1133 | & x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0 |
|---|
| 1134 | else |
|---|
| 1135 | y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0 |
|---|
| 1136 | end if |
|---|
| 1137 | end do |
|---|
| 1138 | |
|---|
| 1139 | elseif (ismooth.eq.2) then ! 3 points, evenly spaced |
|---|
| 1140 | |
|---|
| 1141 | do i=imin,imax |
|---|
| 1142 | if(i.eq.imin)then |
|---|
| 1143 | y(i)=x(imin) |
|---|
| 1144 | elseif(i.eq.imax)then |
|---|
| 1145 | y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 |
|---|
| 1146 | else |
|---|
| 1147 | y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 |
|---|
| 1148 | end if |
|---|
| 1149 | end do |
|---|
| 1150 | |
|---|
| 1151 | elseif (ismooth.eq.3) then ! 5 points, evenly spaced |
|---|
| 1152 | |
|---|
| 1153 | do i=imin,imax |
|---|
| 1154 | if(i.eq.imin)then |
|---|
| 1155 | y(i) = x(imin) |
|---|
| 1156 | elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then |
|---|
| 1157 | y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 |
|---|
| 1158 | elseif(i.eq.imax)then |
|---|
| 1159 | y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0 |
|---|
| 1160 | else |
|---|
| 1161 | y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0 |
|---|
| 1162 | end if |
|---|
| 1163 | end do |
|---|
| 1164 | |
|---|
| 1165 | else |
|---|
| 1166 | |
|---|
| 1167 | write (*,*) ' Error in suaviza.f Wrong ismooth value.' |
|---|
| 1168 | stop |
|---|
| 1169 | |
|---|
| 1170 | endif |
|---|
| 1171 | |
|---|
| 1172 | c rehago el cambio, para devolver x(i) |
|---|
| 1173 | do i=imin,imax |
|---|
| 1174 | x(i)=y(i) |
|---|
| 1175 | end do |
|---|
| 1176 | |
|---|
| 1177 | return |
|---|
| 1178 | end |
|---|
| 1179 | |
|---|
| 1180 | |
|---|
| 1181 | |
|---|
| 1182 | |
|---|
| 1183 | c***************************************************************************** |
|---|
| 1184 | c LUdec.F (includes lubksb_dp and ludcmp_dp subroutines) |
|---|
| 1185 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1186 | c |
|---|
| 1187 | c Solution of linear equation without inverting matrix |
|---|
| 1188 | c using LU decomposition: |
|---|
| 1189 | c AA * xx = bb AA, bb: known |
|---|
| 1190 | c xx: to be found |
|---|
| 1191 | c AA and bb are not modified in this subroutine |
|---|
| 1192 | c |
|---|
| 1193 | c MALV , Sep 2007 |
|---|
| 1194 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1195 | |
|---|
| 1196 | subroutine LUdec(xx,aa,bb,m,n) |
|---|
| 1197 | |
|---|
| 1198 | implicit none |
|---|
| 1199 | |
|---|
| 1200 | ! Arguments |
|---|
| 1201 | integer,intent(in) :: m, n |
|---|
| 1202 | real*8,intent(in) :: aa(m,m), bb(m) |
|---|
| 1203 | real*8,intent(out) :: xx(m) |
|---|
| 1204 | |
|---|
| 1205 | |
|---|
| 1206 | ! Local variables |
|---|
| 1207 | real*8 a(n,n), b(n), x(n), d |
|---|
| 1208 | integer i, j, indx(n) |
|---|
| 1209 | |
|---|
| 1210 | |
|---|
| 1211 | ! Subrutinas utilizadas |
|---|
| 1212 | ! ludcmp_dp, lubksb_dp |
|---|
| 1213 | |
|---|
| 1214 | !!!!!!!!!!!!!!! Comienza el programa !!!!!!!!!!!!!! |
|---|
| 1215 | |
|---|
| 1216 | do i=1,n |
|---|
| 1217 | b(i) = bb(i+1) |
|---|
| 1218 | do j=1,n |
|---|
| 1219 | a(i,j) = aa(i+1,j+1) |
|---|
| 1220 | enddo |
|---|
| 1221 | enddo |
|---|
| 1222 | |
|---|
| 1223 | ! Descomposicion de auxm1 |
|---|
| 1224 | call ludcmp_dp ( a, n, n, indx, d) |
|---|
| 1225 | |
|---|
| 1226 | ! Sustituciones foward y backwards para hallar la solucion |
|---|
| 1227 | do i=1,n |
|---|
| 1228 | x(i) = b(i) |
|---|
| 1229 | enddo |
|---|
| 1230 | call lubksb_dp( a, n, n, indx, x ) |
|---|
| 1231 | |
|---|
| 1232 | do i=1,n |
|---|
| 1233 | xx(i+1) = x(i) |
|---|
| 1234 | enddo |
|---|
| 1235 | |
|---|
| 1236 | return |
|---|
| 1237 | end |
|---|
| 1238 | |
|---|
| 1239 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1240 | |
|---|
| 1241 | subroutine ludcmp_dp(a,n,np,indx,d) |
|---|
| 1242 | |
|---|
| 1243 | c jul 2011 malv+fgg |
|---|
| 1244 | |
|---|
| 1245 | implicit none |
|---|
| 1246 | |
|---|
| 1247 | integer,intent(in) :: n, np |
|---|
| 1248 | real*8,intent(inout) :: a(np,np) |
|---|
| 1249 | real*8,intent(out) :: d |
|---|
| 1250 | integer,intent(out) :: indx(n) |
|---|
| 1251 | |
|---|
| 1252 | integer i, j, k, imax |
|---|
| 1253 | real*8,parameter :: tiny=1.0d-20 |
|---|
| 1254 | real*8 vv(n), aamax, sum, dum |
|---|
| 1255 | |
|---|
| 1256 | |
|---|
| 1257 | d=1.0d0 |
|---|
| 1258 | do 12 i=1,n |
|---|
| 1259 | aamax=0.0d0 |
|---|
| 1260 | do 11 j=1,n |
|---|
| 1261 | if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) |
|---|
| 1262 | 11 continue |
|---|
| 1263 | if (aamax.eq.0.0) then |
|---|
| 1264 | write(*,*) 'ludcmp_dp: singular matrix!' |
|---|
| 1265 | stop |
|---|
| 1266 | endif |
|---|
| 1267 | vv(i)=1.0d0/aamax |
|---|
| 1268 | 12 continue |
|---|
| 1269 | do 19 j=1,n |
|---|
| 1270 | if (j.gt.1) then |
|---|
| 1271 | do 14 i=1,j-1 |
|---|
| 1272 | sum=a(i,j) |
|---|
| 1273 | if (i.gt.1)then |
|---|
| 1274 | do 13 k=1,i-1 |
|---|
| 1275 | sum=sum-a(i,k)*a(k,j) |
|---|
| 1276 | 13 continue |
|---|
| 1277 | a(i,j)=sum |
|---|
| 1278 | endif |
|---|
| 1279 | 14 continue |
|---|
| 1280 | endif |
|---|
| 1281 | aamax=0.0d0 |
|---|
| 1282 | do 16 i=j,n |
|---|
| 1283 | sum=a(i,j) |
|---|
| 1284 | if (j.gt.1)then |
|---|
| 1285 | do 15 k=1,j-1 |
|---|
| 1286 | sum=sum-a(i,k)*a(k,j) |
|---|
| 1287 | 15 continue |
|---|
| 1288 | a(i,j)=sum |
|---|
| 1289 | endif |
|---|
| 1290 | dum=vv(i)*abs(sum) |
|---|
| 1291 | if (dum.ge.aamax) then |
|---|
| 1292 | imax=i |
|---|
| 1293 | aamax=dum |
|---|
| 1294 | endif |
|---|
| 1295 | 16 continue |
|---|
| 1296 | if (j.ne.imax)then |
|---|
| 1297 | do 17 k=1,n |
|---|
| 1298 | dum=a(imax,k) |
|---|
| 1299 | a(imax,k)=a(j,k) |
|---|
| 1300 | a(j,k)=dum |
|---|
| 1301 | 17 continue |
|---|
| 1302 | d=-d |
|---|
| 1303 | vv(imax)=vv(j) |
|---|
| 1304 | endif |
|---|
| 1305 | indx(j)=imax |
|---|
| 1306 | if(j.ne.n)then |
|---|
| 1307 | if(a(j,j).eq.0.0)a(j,j)=tiny |
|---|
| 1308 | dum=1.0d0/a(j,j) |
|---|
| 1309 | do 18 i=j+1,n |
|---|
| 1310 | a(i,j)=a(i,j)*dum |
|---|
| 1311 | 18 continue |
|---|
| 1312 | endif |
|---|
| 1313 | 19 continue |
|---|
| 1314 | if(a(n,n).eq.0.0)a(n,n)=tiny |
|---|
| 1315 | return |
|---|
| 1316 | end |
|---|
| 1317 | |
|---|
| 1318 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1319 | |
|---|
| 1320 | subroutine lubksb_dp(a,n,np,indx,b) |
|---|
| 1321 | |
|---|
| 1322 | c jul 2011 malv+fgg |
|---|
| 1323 | |
|---|
| 1324 | implicit none |
|---|
| 1325 | |
|---|
| 1326 | integer,intent(in) :: n,np |
|---|
| 1327 | real*8,intent(in) :: a(np,np) |
|---|
| 1328 | integer,intent(in) :: indx(n) |
|---|
| 1329 | real*8,intent(out) :: b(n) |
|---|
| 1330 | |
|---|
| 1331 | real*8 sum |
|---|
| 1332 | integer ii, ll, i, j |
|---|
| 1333 | |
|---|
| 1334 | ii=0 |
|---|
| 1335 | do 12 i=1,n |
|---|
| 1336 | ll=indx(i) |
|---|
| 1337 | sum=b(ll) |
|---|
| 1338 | b(ll)=b(i) |
|---|
| 1339 | if (ii.ne.0)then |
|---|
| 1340 | do 11 j=ii,i-1 |
|---|
| 1341 | sum=sum-a(i,j)*b(j) |
|---|
| 1342 | 11 continue |
|---|
| 1343 | else if (sum.ne.0.0) then |
|---|
| 1344 | ii=i |
|---|
| 1345 | endif |
|---|
| 1346 | b(i)=sum |
|---|
| 1347 | 12 continue |
|---|
| 1348 | do 14 i=n,1,-1 |
|---|
| 1349 | sum=b(i) |
|---|
| 1350 | if(i.lt.n)then |
|---|
| 1351 | do 13 j=i+1,n |
|---|
| 1352 | sum=sum-a(i,j)*b(j) |
|---|
| 1353 | 13 continue |
|---|
| 1354 | endif |
|---|
| 1355 | b(i)=sum/a(i,i) |
|---|
| 1356 | 14 continue |
|---|
| 1357 | return |
|---|
| 1358 | end |
|---|
| 1359 | |
|---|
| 1360 | |
|---|
| 1361 | |
|---|
| 1362 | |
|---|
| 1363 | c***************************************************************************** |
|---|
| 1364 | c intersp |
|---|
| 1365 | c *********************************************************************** |
|---|
| 1366 | subroutine intersp(yy,zz,m,y,z,n,opt) |
|---|
| 1367 | c interpolation soubroutine. input values: y(n) at z(n). |
|---|
| 1368 | c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic |
|---|
| 1369 | |
|---|
| 1370 | c jul 2011 malv+fgg |
|---|
| 1371 | c *********************************************************************** |
|---|
| 1372 | |
|---|
| 1373 | implicit none |
|---|
| 1374 | |
|---|
| 1375 | integer n,m,i,j,opt |
|---|
| 1376 | real zz(m),yy(m),z(n),y(n) |
|---|
| 1377 | real zmin,zzmin,zmax,zzmax |
|---|
| 1378 | |
|---|
| 1379 | ! write(*,*) ' interpolating' |
|---|
| 1380 | ! call minsp(z,n,zmin) |
|---|
| 1381 | zmin=minval(z) |
|---|
| 1382 | ! call minsp(zz,m,zzmin) |
|---|
| 1383 | zzmin=minval(zz) |
|---|
| 1384 | ! call maxsp(z,n,zmax) |
|---|
| 1385 | zmax=maxval(z) |
|---|
| 1386 | ! call maxsp(zz,m,zzmax) |
|---|
| 1387 | zzmax=maxval(zz) |
|---|
| 1388 | |
|---|
| 1389 | if(zzmin.lt.zmin)then |
|---|
| 1390 | write(*,*) 'from interp: new variable out of limits' |
|---|
| 1391 | write(*,*) zzmin,'must be .ge. ',zmin |
|---|
| 1392 | stop |
|---|
| 1393 | ! elseif(zzmax.gt.zmax)then |
|---|
| 1394 | ! write(*,*)'from interp: new variable out of limits' |
|---|
| 1395 | ! write(*,*)zzmax, 'must be .le. ',zmax |
|---|
| 1396 | ! stop |
|---|
| 1397 | end if |
|---|
| 1398 | |
|---|
| 1399 | do 1,i=1,m |
|---|
| 1400 | |
|---|
| 1401 | do 2,j=1,n-1 |
|---|
| 1402 | if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 |
|---|
| 1403 | 2 continue |
|---|
| 1404 | c in this case (zz(m).ge.z(n)) and j leaves the loop with j=n-1+1=n |
|---|
| 1405 | if(opt.eq.1)then |
|---|
| 1406 | yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1)) |
|---|
| 1407 | elseif(opt.eq.2)then |
|---|
| 1408 | if(y(n).eq.0.0.or.y(n-1).eq.0.0)then |
|---|
| 1409 | yy(i)=0.0 |
|---|
| 1410 | else |
|---|
| 1411 | yy(i)=exp(log(y(n-1))+log(y(n)/y(n-1))* |
|---|
| 1412 | @ (zz(i)-z(n-1))/(z(n)-z(n-1))) |
|---|
| 1413 | end if |
|---|
| 1414 | else |
|---|
| 1415 | write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt |
|---|
| 1416 | end if |
|---|
| 1417 | goto 1 |
|---|
| 1418 | 3 continue |
|---|
| 1419 | if(opt.eq.1)then |
|---|
| 1420 | yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1421 | elseif(opt.eq.2)then |
|---|
| 1422 | if(y(j+1).eq.0.0.or.y(j).eq.0.0)then |
|---|
| 1423 | yy(i)=0.0 |
|---|
| 1424 | else |
|---|
| 1425 | yy(i)=exp(log(y(j))+log(y(j+1)/y(j))* |
|---|
| 1426 | @ (zz(i)-z(j))/(z(j+1)-z(j))) |
|---|
| 1427 | end if |
|---|
| 1428 | else |
|---|
| 1429 | write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt |
|---|
| 1430 | end if |
|---|
| 1431 | 1 continue |
|---|
| 1432 | |
|---|
| 1433 | return |
|---|
| 1434 | end |
|---|
| 1435 | |
|---|
| 1436 | |
|---|
| 1437 | |
|---|
| 1438 | c***************************************************************************** |
|---|
| 1439 | c interdp |
|---|
| 1440 | c *********************************************************************** |
|---|
| 1441 | subroutine interdp(yy,zz,m,y,z,n,opt) |
|---|
| 1442 | c interpolation soubroutine. input values: y(n) at z(n). |
|---|
| 1443 | c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic |
|---|
| 1444 | c jul 2011: malv+fgg Adapted to LMD-MGCM |
|---|
| 1445 | c *********************************************************************** |
|---|
| 1446 | implicit none |
|---|
| 1447 | integer n,m,i,j,opt |
|---|
| 1448 | real*8 zz(m),yy(m),z(n),y(n), zmin,zzmin,zmax,zzmax |
|---|
| 1449 | |
|---|
| 1450 | ! write (*,*) ' d interpolating ' |
|---|
| 1451 | ! call mindp (z,n,zmin) |
|---|
| 1452 | zmin=minval(z) |
|---|
| 1453 | ! call mindp (zz,m,zzmin) |
|---|
| 1454 | zzmin=minval(zz) |
|---|
| 1455 | ! call maxdp (z,n,zmax) |
|---|
| 1456 | zmax=maxval(z) |
|---|
| 1457 | ! call maxdp (zz,m,zzmax) |
|---|
| 1458 | zzmax=maxval(zz) |
|---|
| 1459 | |
|---|
| 1460 | if(zzmin.lt.zmin)then |
|---|
| 1461 | write (*,*) 'from d interp: new variable out of limits' |
|---|
| 1462 | write (*,*) zzmin,'must be .ge. ',zmin |
|---|
| 1463 | stop |
|---|
| 1464 | ! elseif(zzmax.gt.zmax)then |
|---|
| 1465 | ! write (*,*) 'from interp: new variable out of limits' |
|---|
| 1466 | ! write (*,*) zzmax, 'must be .le. ',zmax |
|---|
| 1467 | ! stop |
|---|
| 1468 | end if |
|---|
| 1469 | |
|---|
| 1470 | do 1,i=1,m |
|---|
| 1471 | |
|---|
| 1472 | do 2,j=1,n-1 |
|---|
| 1473 | if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 |
|---|
| 1474 | 2 continue |
|---|
| 1475 | c in this case (zz(m).eq.z(n)) and j leaves the loop with j=n-1+1=n |
|---|
| 1476 | if(opt.eq.1)then |
|---|
| 1477 | yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1)) |
|---|
| 1478 | elseif(opt.eq.2)then |
|---|
| 1479 | if(y(n).eq.0.0d0.or.y(n-1).eq.0.0d0)then |
|---|
| 1480 | yy(i)=0.0d0 |
|---|
| 1481 | else |
|---|
| 1482 | yy(i)=dexp(dlog(y(n-1))+dlog(y(n)/y(n-1))* |
|---|
| 1483 | @ (zz(i)-z(n-1))/(z(n)-z(n-1))) |
|---|
| 1484 | end if |
|---|
| 1485 | else |
|---|
| 1486 | write (*,*) |
|---|
| 1487 | @ ' from d interp error: opt must be 1 or 2, opt= ',opt |
|---|
| 1488 | end if |
|---|
| 1489 | goto 1 |
|---|
| 1490 | 3 continue |
|---|
| 1491 | if(opt.eq.1)then |
|---|
| 1492 | yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1493 | ! write (*,*) ' ' |
|---|
| 1494 | ! write (*,*) ' z(j),z(j+1) =', z(j),z(j+1) |
|---|
| 1495 | ! write (*,*) ' t(j),t(j+1) =', y(j),y(j+1) |
|---|
| 1496 | ! write (*,*) ' zz, tt = ', zz(i), yy(i) |
|---|
| 1497 | elseif(opt.eq.2)then |
|---|
| 1498 | if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then |
|---|
| 1499 | yy(i)=0.0d0 |
|---|
| 1500 | else |
|---|
| 1501 | yy(i)=dexp(dlog(y(j))+dlog(y(j+1)/y(j))* |
|---|
| 1502 | @ (zz(i)-z(j))/(z(j+1)-z(j))) |
|---|
| 1503 | end if |
|---|
| 1504 | else |
|---|
| 1505 | write (*,*) ' from interp error: opt must be 1 or 2, opt= ', |
|---|
| 1506 | @ opt |
|---|
| 1507 | end if |
|---|
| 1508 | 1 continue |
|---|
| 1509 | return |
|---|
| 1510 | end |
|---|
| 1511 | |
|---|
| 1512 | |
|---|
| 1513 | c***************************************************************************** |
|---|
| 1514 | c interdp_limits.F |
|---|
| 1515 | c *********************************************************************** |
|---|
| 1516 | |
|---|
| 1517 | subroutine interdp_limits ( yy,zz,m, i1,i2, y,z,n, j1,j2, opt) |
|---|
| 1518 | |
|---|
| 1519 | c Interpolation soubroutine. |
|---|
| 1520 | c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m |
|---|
| 1521 | c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n |
|---|
| 1522 | c Input values: y(n) , z(n) (solo se usan los valores entre j1,j2) |
|---|
| 1523 | c zz(m) (solo se necesita entre i1,i2) |
|---|
| 1524 | c Output values: yy(m) (solo se calculan entre i1,i2) |
|---|
| 1525 | c Options: opt=1 -> lineal ,, opt=2 -> logarithmic |
|---|
| 1526 | c Difference with interdp: |
|---|
| 1527 | c here interpolation proceeds between indexes i1,i2 only |
|---|
| 1528 | c if i1=1 & i2=m, both subroutines are exactly the same |
|---|
| 1529 | c thus previous calls to interdp or interdp2 could be easily replaced |
|---|
| 1530 | |
|---|
| 1531 | c JAN 98 MALV Version for mz1d |
|---|
| 1532 | c jul 2011 malv+fgg Adapted to LMD-MGCM |
|---|
| 1533 | c *********************************************************************** |
|---|
| 1534 | |
|---|
| 1535 | implicit none |
|---|
| 1536 | |
|---|
| 1537 | ! Arguments |
|---|
| 1538 | integer n,m ! I. Dimensions |
|---|
| 1539 | integer i1, i2, j1, j2, opt ! I |
|---|
| 1540 | real*8 zz(m),yy(m) ! O |
|---|
| 1541 | real*8 z(n),y(n) ! I |
|---|
| 1542 | |
|---|
| 1543 | ! Local variables |
|---|
| 1544 | integer i,j |
|---|
| 1545 | real*8 zmin,zzmin,zmax,zzmax |
|---|
| 1546 | |
|---|
| 1547 | c ******************************* |
|---|
| 1548 | |
|---|
| 1549 | ! type *, ' d interpolating ' |
|---|
| 1550 | ! call mindp_limits (z,n,zmin, j1,j2) |
|---|
| 1551 | zmin=minval(z(j1:j2)) |
|---|
| 1552 | ! call mindp_limits (zz,m,zzmin, i1,i2) |
|---|
| 1553 | zzmin=minval(zz(i1:i2)) |
|---|
| 1554 | ! call maxdp_limits (z,n,zmax, j1,j2) |
|---|
| 1555 | zmax=maxval(z(j1:j2)) |
|---|
| 1556 | ! call maxdp_limits (zz,m,zzmax, i1,i2) |
|---|
| 1557 | zzmax=maxval(zz(i1:i2)) |
|---|
| 1558 | |
|---|
| 1559 | if(zzmin.lt.zmin)then |
|---|
| 1560 | write (*,*) 'from d interp: new variable out of limits' |
|---|
| 1561 | write (*,*) zzmin,'must be .ge. ',zmin |
|---|
| 1562 | stop |
|---|
| 1563 | ! elseif(zzmax.gt.zmax)then |
|---|
| 1564 | ! type *,'from interp: new variable out of limits' |
|---|
| 1565 | ! type *,zzmax, 'must be .le. ',zmax |
|---|
| 1566 | ! stop |
|---|
| 1567 | end if |
|---|
| 1568 | |
|---|
| 1569 | do 1,i=i1,i2 |
|---|
| 1570 | |
|---|
| 1571 | do 2,j=j1,j2-1 |
|---|
| 1572 | if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 |
|---|
| 1573 | 2 continue |
|---|
| 1574 | c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2 |
|---|
| 1575 | if(opt.eq.1)then |
|---|
| 1576 | yy(i)=y(j2-1)+(y(j2)-y(j2-1))* |
|---|
| 1577 | $ (zz(i)-z(j2-1))/(z(j2)-z(j2-1)) |
|---|
| 1578 | elseif(opt.eq.2)then |
|---|
| 1579 | if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then |
|---|
| 1580 | yy(i)=0.0d0 |
|---|
| 1581 | else |
|---|
| 1582 | yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))* |
|---|
| 1583 | @ (zz(i)-z(j2-1))/(z(j2)-z(j2-1))) |
|---|
| 1584 | end if |
|---|
| 1585 | else |
|---|
| 1586 | write (*,*) ' d interp : opt must be 1 or 2, opt= ',opt |
|---|
| 1587 | end if |
|---|
| 1588 | goto 1 |
|---|
| 1589 | 3 continue |
|---|
| 1590 | if(opt.eq.1)then |
|---|
| 1591 | yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) |
|---|
| 1592 | ! type *, ' ' |
|---|
| 1593 | ! type *, ' z(j),z(j+1) =', z(j),z(j+1) |
|---|
| 1594 | ! type *, ' t(j),t(j+1) =', y(j),y(j+1) |
|---|
| 1595 | ! type *, ' zz, tt = ', zz(i), yy(i) |
|---|
| 1596 | elseif(opt.eq.2)then |
|---|
| 1597 | if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then |
|---|
| 1598 | yy(i)=0.0d0 |
|---|
| 1599 | else |
|---|
| 1600 | yy(i)=exp(log(y(j))+log(y(j+1)/y(j))* |
|---|
| 1601 | @ (zz(i)-z(j))/(z(j+1)-z(j))) |
|---|
| 1602 | end if |
|---|
| 1603 | else |
|---|
| 1604 | write (*,*) ' interp : opt must be 1 or 2, opt= ',opt |
|---|
| 1605 | end if |
|---|
| 1606 | 1 continue |
|---|
| 1607 | return |
|---|
| 1608 | end |
|---|
| 1609 | |
|---|
| 1610 | |
|---|
| 1611 | |
|---|
| 1612 | |
|---|
| 1613 | c***************************************************************************** |
|---|
| 1614 | c Subroutines previously included in tcrco2_subrut.F |
|---|
| 1615 | c*********************************************************************** |
|---|
| 1616 | c tcrco2_subrut.f |
|---|
| 1617 | c |
|---|
| 1618 | c jan 98 malv version for mz1d. copied from solar10/mz4sub.f |
|---|
| 1619 | c jul 2011 malv+fgg adapted to LMD-MGCM |
|---|
| 1620 | c*********************************************************************** |
|---|
| 1621 | |
|---|
| 1622 | ************************************************************************ |
|---|
| 1623 | |
|---|
| 1624 | subroutine dinterconnection ( v, vt ) |
|---|
| 1625 | |
|---|
| 1626 | * input: vib. temp. from che*.for programs, vt(nl) |
|---|
| 1627 | * output: test vibrational temp. for other che*.for, v(nl) |
|---|
| 1628 | ! iconex_smooth=1 ==> with smoothing |
|---|
| 1629 | ! iconex_smooth=0 ==> without smoothing |
|---|
| 1630 | ! iconex_tk=40 ==> with forced lte up to 40 km |
|---|
| 1631 | ! iconex_tk=20 ==> with forced lte up to 20 km |
|---|
| 1632 | ************************************************************************ |
|---|
| 1633 | |
|---|
| 1634 | implicit none |
|---|
| 1635 | include 'nlte_paramdef.h' |
|---|
| 1636 | include 'nlte_commons.h' |
|---|
| 1637 | |
|---|
| 1638 | c argumentos |
|---|
| 1639 | real*8 vt(nl), v(nl) |
|---|
| 1640 | |
|---|
| 1641 | c local variables |
|---|
| 1642 | integer i |
|---|
| 1643 | |
|---|
| 1644 | c ************* |
|---|
| 1645 | |
|---|
| 1646 | do i=1,nl |
|---|
| 1647 | v(i) = vt(i) |
|---|
| 1648 | end do |
|---|
| 1649 | |
|---|
| 1650 | ! lo siguiente se utilizaba en solar10, pero es mejor introducirlo en |
|---|
| 1651 | ! la driver. por ahora no lo uso todavia. |
|---|
| 1652 | ! call fluctua(v,iconex_fluctua) |
|---|
| 1653 | ! call smooth_nl(v,iconex_smooth,nl) |
|---|
| 1654 | ! call forzar_tk(v,iconex_tk) |
|---|
| 1655 | |
|---|
| 1656 | return |
|---|
| 1657 | end |
|---|
| 1658 | |
|---|
| 1659 | c*********************************************************************** |
|---|
| 1660 | function planckdp(tp,xnu) |
|---|
| 1661 | c returns the black body function at wavenumber xnu and temperature t. |
|---|
| 1662 | c*********************************************************************** |
|---|
| 1663 | |
|---|
| 1664 | implicit none |
|---|
| 1665 | |
|---|
| 1666 | include 'nlte_paramdef.h' |
|---|
| 1667 | include 'nlte_commons.h' |
|---|
| 1668 | |
|---|
| 1669 | ! common/datis/ pi, vlight, ee, hplanck, gamma, ab, |
|---|
| 1670 | ! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg |
|---|
| 1671 | ! real*8 pi, vlight, ee, hplanck, gamma, ab, |
|---|
| 1672 | ! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg |
|---|
| 1673 | |
|---|
| 1674 | real*8 planckdp |
|---|
| 1675 | real*8 xnu |
|---|
| 1676 | real tp |
|---|
| 1677 | |
|---|
| 1678 | planckdp = gamma*xnu**3.0 / exp( ee*xnu/dble(tp) ) |
|---|
| 1679 | !erg cm-2.sr-1/cm-1. |
|---|
| 1680 | |
|---|
| 1681 | return |
|---|
| 1682 | end |
|---|
| 1683 | |
|---|
| 1684 | c **************************************************************** |
|---|
| 1685 | function bandid (ib) |
|---|
| 1686 | c returns the 2 character code of the band |
|---|
| 1687 | c **************************************************************** |
|---|
| 1688 | implicit none |
|---|
| 1689 | |
|---|
| 1690 | integer ib |
|---|
| 1691 | character*2 bandid |
|---|
| 1692 | |
|---|
| 1693 | 132 format(i2) |
|---|
| 1694 | ! encode (2,132,bandid) ib |
|---|
| 1695 | write ( bandid, 132) ib |
|---|
| 1696 | |
|---|
| 1697 | if ( ib .eq. 1 ) bandid = '01' |
|---|
| 1698 | if ( ib .eq. 2 ) bandid = '02' |
|---|
| 1699 | if ( ib .eq. 3 ) bandid = '03' |
|---|
| 1700 | if ( ib .eq. 4 ) bandid = '04' |
|---|
| 1701 | if ( ib .eq. 5 ) bandid = '05' |
|---|
| 1702 | if ( ib .eq. 6 ) bandid = '06' |
|---|
| 1703 | if ( ib .eq. 7 ) bandid = '07' |
|---|
| 1704 | if ( ib .eq. 8 ) bandid = '08' |
|---|
| 1705 | if ( ib .eq. 9 ) bandid = '09' |
|---|
| 1706 | if ( ib .eq. 0 ) bandid = '00' |
|---|
| 1707 | |
|---|
| 1708 | c end |
|---|
| 1709 | return |
|---|
| 1710 | end |
|---|
| 1711 | |
|---|
| 1712 | |
|---|
| 1713 | |
|---|
| 1714 | c***************************************************************************** |
|---|
| 1715 | c Subroutines previously included in mat_oper.F |
|---|
| 1716 | c***************************************************************************** |
|---|
| 1717 | c set of subroutines for the cz*.for programs: |
|---|
| 1718 | ! subroutine unit(a,n) |
|---|
| 1719 | ! subroutine diago(a,v,n) diagonal matrix with v |
|---|
| 1720 | ! subroutine invdiag(a,b,n) inverse of diagonal matrix |
|---|
| 1721 | ! subroutine sypvvv(a,b,c,d,n) suma y prod de 3 vectores, muy comun |
|---|
| 1722 | ! subroutine sypvmv(v,w,b,u,n) suma y prod de 3 vectores, muy comun |
|---|
| 1723 | ! subroutine mulmvv(w,b,u,v,n) prod matriz vector vector |
|---|
| 1724 | ! subroutine muymvv(w,b,u,v,n) prod matriz (inv.vector) vector |
|---|
| 1725 | ! subroutine samem (a,m,n) |
|---|
| 1726 | ! subroutine mulmv(a,b,c,n) |
|---|
| 1727 | ! subroutine mulmm(a,b,c,n) |
|---|
| 1728 | ! subroutine resmm(a,b,c,n) |
|---|
| 1729 | ! subroutine mulvv(a,b,c,n) |
|---|
| 1730 | ! subroutine sumvv(a,b,c,n) |
|---|
| 1731 | ! subroutine zerom(a,n) |
|---|
| 1732 | ! subroutine zero4m(a,b,c,d,n) |
|---|
| 1733 | ! subroutine zero3m(a,b,c,n) |
|---|
| 1734 | ! subroutine zero2m(a,b,n) |
|---|
| 1735 | ! subroutine zerov(a,n) |
|---|
| 1736 | ! subroutine zero4v(a,b,c,d,n) |
|---|
| 1737 | ! subroutine zero3v(a,b,c,n) |
|---|
| 1738 | ! subroutine zero2v(a,b,n) |
|---|
| 1739 | |
|---|
| 1740 | ! |
|---|
| 1741 | ! |
|---|
| 1742 | ! May-05 Sustituimos todos los zerojt de cristina por las subrutinas |
|---|
| 1743 | ! genericas zerov*** |
|---|
| 1744 | ! |
|---|
| 1745 | c *********************************************************************** |
|---|
| 1746 | subroutine unit(a,n) |
|---|
| 1747 | c store the unit value in the diagonal of a |
|---|
| 1748 | c *********************************************************************** |
|---|
| 1749 | real*8 a(n,n) |
|---|
| 1750 | integer n,i,j,k |
|---|
| 1751 | do 1,i=2,n-1 |
|---|
| 1752 | do 2,j=2,n-1 |
|---|
| 1753 | if(i.eq.j) then |
|---|
| 1754 | a(i,j) = 1.d0 |
|---|
| 1755 | else |
|---|
| 1756 | a(i,j)=0.0d0 |
|---|
| 1757 | end if |
|---|
| 1758 | 2 continue |
|---|
| 1759 | 1 continue |
|---|
| 1760 | do k=1,n |
|---|
| 1761 | a(n,k) = 0.0d0 |
|---|
| 1762 | a(1,k) = 0.0d0 |
|---|
| 1763 | a(k,1) = 0.0d0 |
|---|
| 1764 | a(k,n) = 0.0d0 |
|---|
| 1765 | end do |
|---|
| 1766 | return |
|---|
| 1767 | end |
|---|
| 1768 | |
|---|
| 1769 | c *********************************************************************** |
|---|
| 1770 | subroutine diago(a,v,n) |
|---|
| 1771 | c store the vector v in the diagonal elements of the square matrix a |
|---|
| 1772 | c *********************************************************************** |
|---|
| 1773 | implicit none |
|---|
| 1774 | |
|---|
| 1775 | integer n,i,j,k |
|---|
| 1776 | real*8 a(n,n),v(n) |
|---|
| 1777 | |
|---|
| 1778 | do 1,i=2,n-1 |
|---|
| 1779 | do 2,j=2,n-1 |
|---|
| 1780 | if(i.eq.j) then |
|---|
| 1781 | a(i,j) = v(i) |
|---|
| 1782 | else |
|---|
| 1783 | a(i,j)=0.0d0 |
|---|
| 1784 | end if |
|---|
| 1785 | 2 continue |
|---|
| 1786 | 1 continue |
|---|
| 1787 | do k=1,n |
|---|
| 1788 | a(n,k) = 0.0d0 |
|---|
| 1789 | a(1,k) = 0.0d0 |
|---|
| 1790 | a(k,1) = 0.0d0 |
|---|
| 1791 | a(k,n) = 0.0d0 |
|---|
| 1792 | end do |
|---|
| 1793 | return |
|---|
| 1794 | end |
|---|
| 1795 | |
|---|
| 1796 | c *********************************************************************** |
|---|
| 1797 | subroutine samem (a,m,n) |
|---|
| 1798 | c store the matrix m in the matrix a |
|---|
| 1799 | c *********************************************************************** |
|---|
| 1800 | real*8 a(n,n),m(n,n) |
|---|
| 1801 | integer n,i,j,k |
|---|
| 1802 | do 1,i=2,n-1 |
|---|
| 1803 | do 2,j=2,n-1 |
|---|
| 1804 | a(i,j) = m(i,j) |
|---|
| 1805 | 2 continue |
|---|
| 1806 | 1 continue |
|---|
| 1807 | do k=1,n |
|---|
| 1808 | a(n,k) = 0.0d0 |
|---|
| 1809 | a(1,k) = 0.0d0 |
|---|
| 1810 | a(k,1) = 0.0d0 |
|---|
| 1811 | a(k,n) = 0.0d0 |
|---|
| 1812 | end do |
|---|
| 1813 | return |
|---|
| 1814 | end |
|---|
| 1815 | c *********************************************************************** |
|---|
| 1816 | subroutine mulmv(a,b,c,n) |
|---|
| 1817 | c do a(i)=b(i,j)*c(j). a, b, and c must be distint |
|---|
| 1818 | c *********************************************************************** |
|---|
| 1819 | implicit none |
|---|
| 1820 | |
|---|
| 1821 | integer n,i,j |
|---|
| 1822 | real*8 a(n),b(n,n),c(n),sum |
|---|
| 1823 | |
|---|
| 1824 | do 1,i=2,n-1 |
|---|
| 1825 | sum=0.0d0 |
|---|
| 1826 | do 2,j=2,n-1 |
|---|
| 1827 | sum=sum+ (b(i,j)) * (c(j)) |
|---|
| 1828 | 2 continue |
|---|
| 1829 | a(i)=sum |
|---|
| 1830 | 1 continue |
|---|
| 1831 | a(1) = 0.0d0 |
|---|
| 1832 | a(n) = 0.0d0 |
|---|
| 1833 | return |
|---|
| 1834 | end |
|---|
| 1835 | |
|---|
| 1836 | c *********************************************************************** |
|---|
| 1837 | subroutine mulmm(a,b,c,n) |
|---|
| 1838 | c *********************************************************************** |
|---|
| 1839 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 1840 | integer n,i,j,k |
|---|
| 1841 | |
|---|
| 1842 | ! do i=2,n-1 |
|---|
| 1843 | ! do j=2,n-1 |
|---|
| 1844 | ! a(i,j)= 0.d00 |
|---|
| 1845 | ! do k=2,n-1 |
|---|
| 1846 | ! a(i,j) = a(i,j) + b(i,k) * c(k,j) |
|---|
| 1847 | ! end do |
|---|
| 1848 | ! end do |
|---|
| 1849 | ! end do |
|---|
| 1850 | do j=2,n-1 |
|---|
| 1851 | do i=2,n-1 |
|---|
| 1852 | a(i,j)=0.d0 |
|---|
| 1853 | enddo |
|---|
| 1854 | do k=2,n-1 |
|---|
| 1855 | do i=2,n-1 |
|---|
| 1856 | a(i,j)=a(i,j)+b(i,k)*c(k,j) |
|---|
| 1857 | enddo |
|---|
| 1858 | enddo |
|---|
| 1859 | enddo |
|---|
| 1860 | do k=1,n |
|---|
| 1861 | a(n,k) = 0.0d0 |
|---|
| 1862 | a(1,k) = 0.0d0 |
|---|
| 1863 | a(k,1) = 0.0d0 |
|---|
| 1864 | a(k,n) = 0.0d0 |
|---|
| 1865 | end do |
|---|
| 1866 | |
|---|
| 1867 | return |
|---|
| 1868 | end |
|---|
| 1869 | |
|---|
| 1870 | c *********************************************************************** |
|---|
| 1871 | subroutine resmm(a,b,c,n) |
|---|
| 1872 | c *********************************************************************** |
|---|
| 1873 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 1874 | integer n,i,j,k |
|---|
| 1875 | |
|---|
| 1876 | do i=2,n-1 |
|---|
| 1877 | do j=2,n-1 |
|---|
| 1878 | a(i,j)= b(i,j) - c(i,j) |
|---|
| 1879 | end do |
|---|
| 1880 | end do |
|---|
| 1881 | do k=1,n |
|---|
| 1882 | a(n,k) = 0.0d0 |
|---|
| 1883 | a(1,k) = 0.0d0 |
|---|
| 1884 | a(k,1) = 0.0d0 |
|---|
| 1885 | a(k,n) = 0.0d0 |
|---|
| 1886 | end do |
|---|
| 1887 | |
|---|
| 1888 | return |
|---|
| 1889 | end |
|---|
| 1890 | |
|---|
| 1891 | c *********************************************************************** |
|---|
| 1892 | subroutine sumvv(a,b,c,n) |
|---|
| 1893 | c a(i)=b(i)+c(i) |
|---|
| 1894 | c *********************************************************************** |
|---|
| 1895 | implicit none |
|---|
| 1896 | |
|---|
| 1897 | integer n,i |
|---|
| 1898 | real*8 a(n),b(n),c(n) |
|---|
| 1899 | |
|---|
| 1900 | do 1,i=2,n-1 |
|---|
| 1901 | a(i)= (b(i)) + (c(i)) |
|---|
| 1902 | 1 continue |
|---|
| 1903 | a(1) = 0.0d0 |
|---|
| 1904 | a(n) = 0.0d0 |
|---|
| 1905 | return |
|---|
| 1906 | end |
|---|
| 1907 | |
|---|
| 1908 | c *********************************************************************** |
|---|
| 1909 | subroutine zerom(a,n) |
|---|
| 1910 | c a(i,j)= 0.0 |
|---|
| 1911 | c *********************************************************************** |
|---|
| 1912 | |
|---|
| 1913 | implicit none |
|---|
| 1914 | |
|---|
| 1915 | integer n,i,j |
|---|
| 1916 | real*8 a(n,n) |
|---|
| 1917 | |
|---|
| 1918 | do 1,i=1,n |
|---|
| 1919 | do 2,j=1,n |
|---|
| 1920 | a(i,j) = 0.0d0 |
|---|
| 1921 | 2 continue |
|---|
| 1922 | 1 continue |
|---|
| 1923 | return |
|---|
| 1924 | end |
|---|
| 1925 | |
|---|
| 1926 | c *********************************************************************** |
|---|
| 1927 | subroutine zero4m(a,b,c,d,n) |
|---|
| 1928 | c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 |
|---|
| 1929 | c *********************************************************************** |
|---|
| 1930 | real*8 a(n,n), b(n,n), c(n,n), d(n,n) |
|---|
| 1931 | integer n,i,j |
|---|
| 1932 | do 1,i=1,n |
|---|
| 1933 | do 2,j=1,n |
|---|
| 1934 | a(i,j) = 0.0d0 |
|---|
| 1935 | b(i,j) = 0.0d0 |
|---|
| 1936 | c(i,j) = 0.0d0 |
|---|
| 1937 | d(i,j) = 0.0d0 |
|---|
| 1938 | 2 continue |
|---|
| 1939 | 1 continue |
|---|
| 1940 | return |
|---|
| 1941 | end |
|---|
| 1942 | |
|---|
| 1943 | c *********************************************************************** |
|---|
| 1944 | subroutine zero3m(a,b,c,n) |
|---|
| 1945 | c a(i,j) = b(i,j) = c(i,j) = 0.0 |
|---|
| 1946 | c ********************************************************************** |
|---|
| 1947 | real*8 a(n,n), b(n,n), c(n,n) |
|---|
| 1948 | integer n,i,j |
|---|
| 1949 | do 1,i=1,n |
|---|
| 1950 | do 2,j=1,n |
|---|
| 1951 | a(i,j) = 0.0d0 |
|---|
| 1952 | b(i,j) = 0.0d0 |
|---|
| 1953 | c(i,j) = 0.0d0 |
|---|
| 1954 | 2 continue |
|---|
| 1955 | 1 continue |
|---|
| 1956 | return |
|---|
| 1957 | end |
|---|
| 1958 | |
|---|
| 1959 | c *********************************************************************** |
|---|
| 1960 | subroutine zero2m(a,b,n) |
|---|
| 1961 | c a(i,j) = b(i,j) = 0.0 |
|---|
| 1962 | c *********************************************************************** |
|---|
| 1963 | real*8 a(n,n), b(n,n) |
|---|
| 1964 | integer n,i,j |
|---|
| 1965 | do 1,i=1,n |
|---|
| 1966 | do 2,j=1,n |
|---|
| 1967 | a(i,j) = 0.0d0 |
|---|
| 1968 | b(i,j) = 0.0d0 |
|---|
| 1969 | 2 continue |
|---|
| 1970 | 1 continue |
|---|
| 1971 | return |
|---|
| 1972 | end |
|---|
| 1973 | c *********************************************************************** |
|---|
| 1974 | subroutine zerov(a,n) |
|---|
| 1975 | c a(i)= 0.0 |
|---|
| 1976 | c *********************************************************************** |
|---|
| 1977 | real*8 a(n) |
|---|
| 1978 | integer n,i |
|---|
| 1979 | do 1,i=1,n |
|---|
| 1980 | a(i) = 0.0d0 |
|---|
| 1981 | 1 continue |
|---|
| 1982 | return |
|---|
| 1983 | end |
|---|
| 1984 | c *********************************************************************** |
|---|
| 1985 | subroutine zero4v(a,b,c,d,n) |
|---|
| 1986 | c a(i) = b(i) = c(i) = d(i,j) = 0.0 |
|---|
| 1987 | c *********************************************************************** |
|---|
| 1988 | real*8 a(n), b(n), c(n), d(n) |
|---|
| 1989 | integer n,i |
|---|
| 1990 | do 1,i=1,n |
|---|
| 1991 | a(i) = 0.0d0 |
|---|
| 1992 | b(i) = 0.0d0 |
|---|
| 1993 | c(i) = 0.0d0 |
|---|
| 1994 | d(i) = 0.0d0 |
|---|
| 1995 | 1 continue |
|---|
| 1996 | return |
|---|
| 1997 | end |
|---|
| 1998 | c *********************************************************************** |
|---|
| 1999 | subroutine zero3v(a,b,c,n) |
|---|
| 2000 | c a(i) = b(i) = c(i) = 0.0 |
|---|
| 2001 | c *********************************************************************** |
|---|
| 2002 | real*8 a(n), b(n), c(n) |
|---|
| 2003 | integer n,i |
|---|
| 2004 | do 1,i=1,n |
|---|
| 2005 | a(i) = 0.0d0 |
|---|
| 2006 | b(i) = 0.0d0 |
|---|
| 2007 | c(i) = 0.0d0 |
|---|
| 2008 | 1 continue |
|---|
| 2009 | return |
|---|
| 2010 | end |
|---|
| 2011 | c *********************************************************************** |
|---|
| 2012 | subroutine zero2v(a,b,n) |
|---|
| 2013 | c a(i) = b(i) = 0.0 |
|---|
| 2014 | c *********************************************************************** |
|---|
| 2015 | real*8 a(n), b(n) |
|---|
| 2016 | integer n,i |
|---|
| 2017 | do 1,i=1,n |
|---|
| 2018 | a(i) = 0.0d0 |
|---|
| 2019 | b(i) = 0.0d0 |
|---|
| 2020 | 1 continue |
|---|
| 2021 | return |
|---|
| 2022 | end |
|---|
| 2023 | |
|---|
| 2024 | |
|---|
