c*********************************************************************** c File with all subroutines required by mztf c Subroutines previously included in mztfsub_overlap.F c c jan 98 malv basado en mztfsub_solar c jul 2011 malv+fgg adapted to LMD-MGCM c c contiene: c initial c intershape c interstrength c intz c rhist c we c simrul c fi c f c findw c voigtf c*********************************************************************** c **************************************************************** subroutine initial c ma & crs !evita troubles 16-july-96 c **************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c local variables integer i c *************** eqw = 0.0d00 aa = 0.0d00 bb = 0.0d00 cc = 0.0d00 dd = 0.0d00 do i=1,nbox ua(i) = 0.0d0 ccbox(i) = 0.0d0 ddbox(i) = 0.0d0 end do return end c ********************************************************************** subroutine intershape(alsx,alnx,adx,xtemp) c interpolates the line shape parameters at a temperature xtemp from c input histogram data. c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max), & xtemp(nbox_max) c local variables integer i, k c *********** ! write (*,*) 'intershape xtemp =', xtemp do 1, k=1,nbox if (xtemp(k).gt.tmax) then write (*,*) ' WARNING ! Tpath,tmax= ',xtemp(k),tmax xtemp(k) = tmax endif if (xtemp(k).lt.tmin) then write (*,*) ' WARNING ! Tpath,tmin= ',xtemp(k),tmin xtemp(k) = tmin endif i = 1 do while (i.le.mm) i = i + 1 if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles alsx(k)=xls1(i,k) !16-july-1996 alnx(k)=xln1(i,k) adx(k)=xld1(i,k) goto 1 elseif ( thist(i) .le. xtemp(k) ) then alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) + @ xls1(i-1,k)*(xtemp(k)-thist(i)) )) / $ (thist(i-1)-thist(i)) alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) + @ xln1(i-1,k)*(xtemp(k)-thist(i)) )) / $ (thist(i-1)-thist(i)) adx(k) = (( xld1(i,k)*(thist(i-1)-xtemp(k)) + @ xld1(i-1,k)*(xtemp(k)-thist(i)) )) / $ (thist(i-1)-thist(i)) goto 1 end if end do write (*,*) @ ' error in xtemp(k). it should be between tmin and tmax' 1 continue return end c ********************************************************************** subroutine interstrength (stx, ts, sx, xtemp) c interpolates the line strength at a temperature xtemp from c input histogram data. c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real*8 stx ! output, total band strength real*8 ts ! input, temp for stx real*8 sx(nbox_max) ! output, strength for each box real*8 xtemp(nbox_max) ! input, temp for sx c local variables integer i, k c *********** do 1, k=1,nbox ! if(xtemp(k).lt.ts)then ! write(*,*)'***********************' ! write(*,*)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k ! write(*,*)'***********************' ! endif if (xtemp(k).gt.tmax) xtemp(k) = tmax if (xtemp(k).lt.tmin) xtemp(k) = tmin i = 1 do while (i.le.mm-1) i = i + 1 ! write(*,*)'mztfsub_overlap/136',i,xtemp(k),thist(i) if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then sx(k) = sk1(i,k) ! write(*,*)'mztfsub_overlap/139',sx(k),k,i goto 1 elseif ( thist(i) .le. xtemp(k) ) then sx(k) = ( sk1(i,k)*(thist(i-1)-xtemp(k)) + sk1(i-1,k)* @ (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i)) ! write(*,*)'mztfsub_overlap/144',sx(k),k,i goto 1 end if end do write (*,*) ' error in xtemp(kr) =', xtemp(k), @ '. it should be between ' write (*,*) ' tmin =',tmin, ' and tmax =',tmax stop 1 continue stx = 0.d0 if (ts.gt.tmax) ts = dble( tmax ) if (ts.lt.tmin) ts = dble( tmin ) i = 1 do while (i.le.mm-1) i = i + 1 ! write(*,*)'mztfsub_overlap/160',i,ts,thist(i) if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then do k=1,nbox stx = stx + no(k) * sk1(i,k) ! write(*,*)'mztfsub_overlap/164',stx end do return elseif ( thist(i) .le. ts ) then do k=1,nbox stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) + @ sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i)) ! write(*,*)'mztfsub_overlap/171',stx end do ! stop return end if end do return end c ********************************************************************** subroutine intz(h,aco2,ap,amr,at, con) c return interp. concentration, pressure,mixing ratio and temperature c for a input height h c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real h ! i real*8 con(nzy) ! i real*8 aco2, ap, at, amr ! o c local variables integer k c ************ if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then write (*,*) ' zp= ',h,' zy(1)= ',zy(1) stop'from intz: error in interpolation, z < minimum height' elseif (h.gt.zy(nzy)) then write (*,*) ' zp= ',h,' zy(nzy)= ',zy(nzy) stop'from intz: error in interpolation, z > maximum height' end if if (h.eq.zy(nzy)) then ap = dble( py(nzy) ) aco2= con(nzy) at = dble( ty(nzy) ) amr = dble( mr(nzy) ) return end if do k=1,nzy-1 if( abs( h-zy(k) ).le.( 1.e-5 ) ) then ap = dble( py(k) ) aco2= con(k) at = dble( ty(k) ) amr = dble( mr(k) ) return elseif(h.gt.zy(k).and.h.lt.zy(k+1))then ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * @ (h-zy(k)) / (zy(k+1)-zy(k)) ) ) aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) * @ (h-zy(k)) / (zy(k+1)-zy(k)) ) at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/ @ (zy(k+1)-zy(k)) ) amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/ @ (zy(k+1)-zy(k)) ) return end if end do return end c ********************************************************************** real*8 function iaa_we(ig,me,pe,plaux,idummy,nt_local,p_local, $ Desp,wsL) c icls=5 -->para mztf c icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2) c calculates an approximate equivalent width for an error estimate. c c ioverlap = 0 ....... no correction for overlaping c 1 ....... "lisat" first correction (see overlap_box. c 2 ....... " " " plus "supersaturation" c idummy=0 do nothing c 1 write out some values for diagnostics c 2 correct the Strong Lorentz behaviour for SZA>90 c 3 casos 1 & 2 c malv nov-98 add overlaping's corrections c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig ! ADDED FOR TRACEBACK real*8 me ! I. path's absorber amount real*8 pe ! I. path's presion total real*8 plaux ! I. path's partial pressure of CO2 real*8 nt_local ! I. needed for strong limit of Lorentz profil real*8 p_local ! I. " " " integer idummy ! I. indica varias opciones real*8 wsL ! O. need this for strong Lorentz correction real*8 Desp ! I. need this for strong Lorentz correction c local variables integer i real*8 y,x,wl,wd real*8 cn(0:7),dn(0:7) real*8 pi, xx real*8 f_sat_box real*8 dv_sat_box, dv_corte_box real*8 area_core_box, area_wing_box real*8 wlgood , parentesis , xlor real*8 wsl_grad c data blocks data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2, @ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4, @ 3.42209457833d-5,-1.28380804108d-6/ data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1, @ 8.21896973657d-1,-2.5222672453,6.1007027481, @ -8.51001627836,4.6535116765/ c *********** c equivalent width of atmospheric line. pi = acos(-1.d0) if ( idummy.gt.9 ) @ write (*,*) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), plaux y=ka(kr)*me ! x=y/(2.0*pi*(alsa(kr)*pl+alna(kr)*(pe-pl))) x=y/(2.0d0*pi* alsa(kr)*plaux) !+alna(kr)*(pe-pl))) ! Strong limit of Lorentz profile: WL = 2 SQRT( S * m * alsa*pl ) ! Para anular esto, comentar las siguientes 5 lineas ! if ( x .gt. 1.e6 ) then ! wl = 2.0*sqrt( y * alsa(kr)*pl ) ! else ! wl=y/sqrt(1.0d0+pi*x/2.0d0) ! endif wl=y/sqrt(1.0d0+pi*x/2.0d0) if (wl .le. 0.d0) then write(*,*)'mztfsub_overlap/496',ig,y,ka(kr),me,kr stop'WE/Lorentz EQW zero or negative!/498' !,ig endif if ( idummy.gt.9 ) @ write (*,*) ' y, x =', y, x xlor = x if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then ! en caso que estemos en el regimen ! Strong Lorentz y la presion local ! vaya disminuyendo, corregimos la EQW ! con un gradiente analitico (notebook) wsL = 2.0*sqrt( y * alsa(kr)*plaux ) wsl_grad = - 2.d0 * ka(kr)*alsa(kr) * nt_local*p_local / wsL wlgood = w_strongLor_prev(kr) + wsl_grad * Desp if (idummy.eq.12) @ write (*,*) ' W(wrong), W_SL, W_SL prev, W_SL corrected=', @ wl, wsL, w_strongLor_prev(kr), wlgood wl = wlgood endif ! wsL = wl pero esto no lo hacemos todavia, porque necesitamos ! el valor que ahora mismo tiene wsL para corregir la ! expresion R&W below ! write (*,*) 'WE arguments me,pe,pl =', me,pe,pl ! write (*,*) 'WE/ wl,ka(kr),alsa(kr) =', ! @ wl, ka(kr),alsa(kr) !>>>>>>> 500 format (a,i3,3(2x,1pe15.8)) 600 format (a,2(2x,1pe16.9)) 700 format (a,3(1x,1pe16.9)) ! if (kr.eq.8 .or. kr.eq.13) then ! write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl ! write (*,700) ' /aln,als,d_x=', alna(kr),alsa(kr), ! @ 2.0*pi*( alsa(kr)*pl + alna(kr)*(pe-pl) ) ! write (*,600) ' /alsa*p_CO2, alna*p_n2 :', ! @ alsa(kr)*pl, alna(kr)*(pe-pl) ! write (*,600) ' a*p, S =', ! @ alsa(kr)*pl + alna(kr)*(pe-pl) , ka(kr) ! write (*,600) ' /S*m, x =', y, x ! write (*,600) ' /aprox, WL =', ! @ 2.*sqrt( y*( alsa(kr)*pl+alna(kr)*(pe-pl) ) ), WL ! endif ! ! corrections to lorentz eqw due to overlaping and super-saturation ! i_supersat = 0 if ( icls.eq.5 .and. ioverlap.gt.0 ) then ! for the moment, only consider overlaping for mztf.f, not fot.f ! definition of saturation in the lisat model ! asat_box = 0.99d0 f_sat_box = 2.d0 * x xx = f_sat_box / log( 1./(1-asat_box) ) if ( xx .lt. 1.0d0 ) then dv_sat_box = 0.0d0 asat_box = 1.0d0 - exp( - f_sat_box ) else dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 ) ! approximation: only use of alsa in mars and venus endif ! area of saturated line ! area_core_box = 2.d0 * dv_sat_box * asat_box area_wing_box = 0.5d0 * ( wl - area_core_box ) dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box ! super-saturation or simple overlaping? ! ! i_supersat = 0 xx = dv_sat_box - ( 0.5d0 * dist(kr) ) if ( xx .ge. 0.0 ! definition of supersaturation @ .and. dv_sat_box .gt. 0.0 ! definition of saturation @ .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line @ ! and not too far apart @ then i_supersat = 1 else ! no super-saturation, then use "lisat + first correction", i.e., ! correct for line products ! wl = wl endif end if ! end of overlaping loop if (icls.eq.2) then iaa_we = wl return endif cc doppler limit: if ( idummy.gt.9 ) @ write (*,*) ' S*m, alf_dop =', y, alda(kr)*sqrt(pi) x = y / (alda(kr)*sqrt(pi)) if ( x.lt.1.e-10 ) then ! to avoid underflow wd = y else wd=alda(kr)*sqrt(4.0*pi*x**2*(1.0+log(1.0+x))/(4.0+pi*x**2)) endif if ( idummy.gt.9 ) @ write (*,*) ' wd =', wd cc doppler weak limit c wd = ka(kr) * me cc good doppler if(icls.eq.5) then !para mztf !write (*,*) 'para mztf, icls=',icls if (x.lt.5.) then wd = 0.d0 do i=0,7 wd = wd + cn(i) * x**i end do wd = alda(kr) * x * sqrt(pi) * wd elseif (x.gt.5.) then wd = 0.d0 do i=0,7 wd = wd + dn(i) / (log(x))**i end do wd = alda(kr) * sqrt(log(x)) * wd else stop ' x should not be less than zero' end if end if if ( i_supersat .eq. 0 ) then parentesis = wl**2+wd**2-(wd*wl/y)**2 ! changed +(wd*wl/y)**2 to -...14-3-84 if ( parentesis .lt. 0.0 ) then if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then parentesis = wl**2+wd**2-(wd*wsL/y)**2 ! este cambio puede ser necesario cuando se hace ! correccion Strong Lor, para evitar valores ! negativos del parentesis en sqrt( ) else stop ' WE/ Error en las EQW wl,wl,y ' endif endif iaa_we = sqrt( parentesis ) ! write (*,*) ' from iaa_we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd) ! write (*,*) ' from iaa_we: we', iaa_we else iaa_we = wl ! if there is supersaturation we can ignore wd completely; ! mztf.f will compute the eqw of the whole box afterwards endif if (icls.eq.3) iaa_we = wd if ( idummy.gt.9 ) @ write (*,*) ' wl,wd,w =', wl,wd,iaa_we wsL = wl return end c ********************************************************************** real*8 function simrul(a,b,fsim,c,acc) c adaptively integrates fsim from a to b, within the criterion acc. c ********************************************************************** implicit none real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim real*8 s1(70),s2(70),s3(70) real*8 c1, c2 integer*4 m,n,j res=0. c=0. m=0 n=0 j=30 g0=fsim(a) g2=fsim((a+b)/2.) g4=fsim(b) a0=(b-a)*(g0+4.0*g2+g4)/2.0 1 d=2.0**n h=(b-a)/(4.0*d) x=a+(4.0*m+1.0)*h g1=fsim(x) g3=fsim(x+2.0*h) a1=h*(g0+4.0*g1+g2) a2=h*(g2+4.0*g3+g4) if ( abs(a1+a2-a0).gt.(acc/d)) goto 2 res=res+(16.0*(a1+a2)-a0)/45.0 m=m+1 c=a+m*(b-a)/d 6 if (m.eq.(2*(m/2))) goto 4 if ((m.ne.1).or.(n.ne.0)) goto 5 8 simrul=res return 2 m=2*m n=n+1 if (n.gt.j) goto 3 a0=a1 s1(n)=a2 s2(n)=g3 s3(n)=g4 g4=g2 g2=g1 goto 1 3 c1=c-(b-a)/d c2=c+(b-a)/d write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2) 7 format(2x,'17hsimrule fails at ',/,3e15.6,/,3e15.6) goto 8 5 a0=s1(n) g0=g4 g2=s2(n) g4=s3(n) goto 1 4 m=m/2 n=n-1 goto 6 end c ********************************************************************** subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL) c this routine sets up accuracy criteria and calls simrule between limit c that depend on the number of atmospheric and cell paths. it gives eqw. c Add correction for EQW in Strong Lorentz regime and SZA>90 c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig ! ADDED FOR TRACEBACK integer iirw integer idummy ! I. indica varias opciones real*8 c1 ! I. needed for strong limit of Lorentz profil real*8 p1 ! I. " " " real*8 wsL ! O. need this for strong Lorentz correction real*8 Desp ! I. need this for strong Lorentz correction c local variables real*8 ept,eps,xa real*8 acc, c real*8 iaa_we real*8 iaa_f, iaa_fi, simrul external iaa_f,iaa_fi c ********** *********** ********* if(icls.eq.5) then !para mztf ! if(ig.eq.1682)write(*,*)'mztfsub_overlap/768',ua(kr),iirw if (iirw.eq.2) then !iirw=icf=2 ==> we use the w&r formula w = iaa_we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL ) return end if ept=iaa_we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL) else !para fot if (iirw.eq.2) then ! icf=2 ==> we use the w&r formula w = iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) return end if ept=iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL) end if c the next block is a modification to avoid nul we. c this situation appears for weak lines and low path temperature, but c there is not any loss of accuracy. first july 1986 if (ept.eq.0.) then ! for weak lines sometimes we=0 ept=1.0e-18 write (*,*) 'ept =',ept write (*,*) 'from we: we=0.0' return end if acc = 4.d0 acc = 10.d0**(-acc) eps = acc * ept !accuracy 10-4 atmospheric eqw. xa=0.5*ept/iaa_f(0.d0) !width of doppler shifted atmospheric line. w = 2.0*( simrul(0.0d0,xa,iaa_f,c,eps) . + simrul(0.1d0,1.0/xa,iaa_fi,c,eps) ) !no shift. return end c ********************************************************************** double precision function iaa_fi(y) c returns the value of f(1/y) c ********************************************************************** implicit none real*8 iaa_f, y iaa_fi=iaa_f(1.0/y)/y**2 return end c ********************************************************************** double precision function iaa_f(nuaux) c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof c ********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' double precision tra,xa,ya,za,yy,nuaux double precision voigtf tra=1.0d0 yy=1.0d0/alda(kr) xa=nuaux*yy ya= ( alsa(kr)*pp + alna(kr)*(pt-pp) ) * yy za=ka(kr)*yy if(icls.eq.5) then !para mztf ! write (*,*) 'icls=',icls tra=za*ua(kr)*voigtf(sngl(xa),sngl(ya)) else tra=za*sl_ua*voigtf(sngl(xa),sngl(ya)) end if if (tra.gt.50.0) then tra=1.0 !2.0e-22 overflow cut-off. else if (tra.gt.1.0e-4) then tra=1.0-exp(-tra) end if iaa_f=tra return end c ********************************************************************** double precision function voigtf(x1,y) c computes voigt function for any value of x1 and any +ve value of y. c where possible uses modified lorentz and modified doppler approximatio c otherwise uses a rearranged rybicki routine. c c(n) = exp(-(n/h)**2)/(pi*sqrt(pi)), with h = 2.5 . c accurate to better than 1 in 10000. c ********************************************************************** implicit none real x1, y real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c integer*4 n dimension c(10) complex xp,xpp,z data c(1)/0.15303405/ data c(2)/0.94694928e-1/ data c(3)/0.42549174e-1/ data c(4)/0.13882935e-1/ data c(5)/0.32892528e-2/ data c(6)/0.56589906e-3/ data c(7)/0.70697890e-4/ data c(8)/0.64135678e-5/ data c(9)/0.42249221e-6/ data c(10)/0.20209868e-7/ x=abs(x1) if (x.gt.7.2) goto 1 if ((y+x*0.3).gt.5.4) goto 1 if (y.gt.0.01) goto 3 if (x.lt.2.1) goto 2 goto 3 c here uses modified lorentz approx. 1 xx=x*x xxyy=xx+y*y b=xx/xxyy voigtf=y*(1.+(2.*b-0.5+(0.75-(9.-12.*b)*b)/xxyy)/ * xxyy)/(xxyy*3.141592654) return c here uses modified doppler approx. 2 xx=x*x voigtf=0.56418958*exp(-xx)*(1.-y*(1.-0.5*y)*(1.1289-xx*(1.1623+ * xx*(0.080812+xx*(0.13854-xx*(0.033605-0.0073972*xx)))))) return c here uses a rearranged rybicki routine. 3 xh=2.5*x xhxh=xh*xh yh=2.5*y yhyh=yh*yh f1=xhxh+yhyh f2=f1-0.5*yhyh if (y.lt.0.1) goto 20 p=-y*7.8539816 !7.8539816=2.5*pi q=x*7.8539816 xpp=cmplx(p,q) z=cexp(xpp) d1=xh*aimag(z) d2=-d1 d3=yh*(1.-real(z)) d4=-d3+2.*yh voig=0.17958712*(d1+d3)/f1 goto 30 20 p=x*7.8539816 q=y*7.8539816 xp=cmplx(p,q) z=ccos(xp) d1=xh*aimag(z) d2=-d1 d3=yh*(1.-real(z)) d4=-d3+2.*yh voig=0.56418958*exp(y*y-x*x)*cos(2.*x*y)+0.17958712*(d1+d3)/f1 30 xn=0. do 55 n=1,10,2 xn=xn+1. xnxn=xn*xn g1=xh-xn g2=g1*(xh+xn) g3=0.5*g2*g2 voig=voig+c(n)*(d2*(g2+yhyh)+d4*(f1+xnxn))/ & (g3+yhyh*(f2+xnxn)) xn=xn+1. xnxn=xn*xn g1=xh-xn g2=g1*(xh+xn) g3=0.5*g2*g2 voig=voig+c(n+1)*(d1*(g2+yhyh)+d3*(f1+xnxn))/ @ (g3+yhyh*(f2+xnxn)) 55 continue voigtf=voig return end c ********************************************************************** c elimin_mz1d.F (includes smooth_cf) c ************************************************************************ subroutine elimin_mz1d (c,vc, ilayer,nanaux,itblout, nwaux) c Eliminate anomalous negative numbers in c(nl,nl) according to "nanaux": c nanaux = 0 -> no eliminate c @ -> eliminate all numbers with absol.value eliminate all anomalous negative numbers in c(n,r). c 3 -> eliminate all anomalous negative numbers far from the main c diagonal. c 8 -> eliminate all non-zero numbers outside the main diagonal, c and the contibution of lower boundary. c 9 -> eliminate all non-zero numbers outside the main diagonal. c 4 -> hace un smoothing cuando la distancia de separacion entre c el valor maximo y el minimo de cf > 50 capas. c 5 -> elimina valores menores que 1.0d-19 c 6 -> incluye los dos casos 4 y 5 c 7 -> llama a lisa: smooth con width=nw & elimina mejorado c 78-> incluye los dos casos 7 y 8 c 79-> incluye los dos casos 7 y 9 c itblout (itableout in calling program) is the option for writing c out or not the purged c(n,r) matrix: c itblout = 0 -> no write c = 1 -> write out in curtis***.out according to ilayer c ilayer is the index for the layer selected to write out the matrix: c ilayer = 0 => matrix elements written out cover all the altitudes c with 5 layers steps c > 0 => " " " " are c(ilayer,*) c NOTA: c EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ) c UTILIZANDO itableout=30 EN MZTUD c jul 2011 malv+fgg adapted to LMD-MGCM c Sep-04 FGG+MALV Correct include and call parameters c cristina 25-sept-1996 y 27-ene-1997 c JAN 98 MALV Version for mz1d c ************************************************************************ implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' integer nanaux,j,i,itblout,kk,k,ir,in integer ilayer,jmin, jmax,np,nwaux,ntimes,ntimes2 !* real*8 c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini real*8 c(nl,nl), vc(nl), amax, cmax, cmin, mini real*8 aux(nl), auxs(nl) character layercode*3 ntimes=0 ntimes2=0 ! type *,'from elimin_mz4: nan, nw',nan,nw if (nanaux .eq. 0) goto 200 if(nanaux.eq.1)then do i=1,nl amax=1.0d-36 do j=1,nl if(abs(c(i,j)).gt.amax)amax=abs(c(i,j)) end do do j=1,nl if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0 end do enddo elseif(nanaux.eq.2)then do i=1,nl do j=1,nl if( ( j.le.(i-2) .or. j.gt.(i+2) ).and. @ ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0 end do enddo elseif(nanaux.eq.3)then do i=1,nl do j=1,nl if (abs(i-j).ge.10) c(i,j)=0.0d0 end do enddo elseif(nanaux.eq.8)then do i=1,nl do j=1,i-1 c(i,j)=0.0d0 enddo do j=i+1,nl c(i,j)=0.0d0 enddo vc(i)= 0.d0 enddo elseif(nanaux.eq.9)then do i=1,nl do j=1,i-1 c(i,j)=0.0d0 enddo do j=i+1,nl c(i,j)=0.0d0 enddo enddo ! elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then ! call lisa(c, vc, nl, nw) end if if(nanaux.eq.78)then do i=1,nl do j=1,i-1 c(i,j)=0.0d0 enddo do j=i+1,nl c(i,j)=0.0d0 enddo vc(i)= 0.d0 enddo endif if(nanaux.eq.79)then do i=1,nl do j=1,i-1 c(i,j)=0.0d0 enddo do j=i+1,nl c(i,j)=0.0d0 enddo enddo endif if(nanaux.eq.5.or.nanaux.eq.6)then do i=1,nl mini = 1.0d-19 do j=1,nl if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then ntimes2=ntimes2+1 end if if ( abs(c(i,j)).le.mini) c(i,j)=0.d0 end do enddo end if if(nanaux.eq.4.or.nanaux.eq.6)then do i=1,nl do j=1,nl aux(j)=c(i,j) auxs(j)=c(i,j) end do !call maxdp_2(aux,nl,cmax,jmax) cmax=maxval(aux) jmax=maxloc(aux,dim=1) if(abs(jmax-i).ge.50) then call smooth_cf(aux,auxs,i,nl,3) !!!call smooth_cf(aux,auxs,i,nl,5) ntimes=ntimes+1 end if do j=1,nl c(i,j)=auxs(j) end do end do end if ! type *, 'elimin_mz4: c(n,r) procesed for elimination. ' ! type *, ' ' ! if(nan.eq.4.or.nan.eq.6) type *, ' call smoothing:',ntimes ! if(nan.eq.5.or.nan.eq.6) type *, ' call elimina: ',ntimes2 ! if(nan.eq.7) type *, ' from elimin: lisa w=',nw ! type *, ' ' 200 continue c writting out of c(n,r) in ascii file ! if(itblout.eq.1) then ! if (ilayer.eq.0) then ! open (unit=2, status='new', ! @ file=dircurtis//'curtis_gnu.out', recl=1024) ! write(2,'(a)') ! @ ' curtis matrix: table with 1.e+7 * acf(n,r) ' ! write(2,114) 'n,r', ( i, i=nl,1,-5 ) ! do in=nl,1,-5 ! write(2,*) ! write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 ) ! end do ! close(2) ! write (*,*) ' ' ! write (*,*) ' curtis.out has been created. ' ! write (*,*) ' ' ! else ! write (layercode,132) ilayer ! open (2, status='new', ! @ file=dircurtis//'curtis'//layercode//'.out') ! write(2,'(a)') ! @ ' curtis matrix: table with 1.e+7 * acf(n,r) ' ! write(2,116) ' layer x c(',layercode, ! @ ',x) c(x,', layercode,')' ! do in=nl,1,-1 ! if (c(ilayer,ilayer).ne.0.d0) then ! write(2,117) in, c(ilayer,in), c(in,ilayer), ! @ c(ilayer,in)/c(ilayer,ilayer), ! @ c(in,ilayer)/c(ilayer,ilayer) ! else ! write(2,118) in, c(ilayer,in), c(in,ilayer) ! end if ! end do ! close(2) ! write (*,*) ' ' ! write (*,*) dircurtis//'curtis'//layercode//'.out', ! @ ' has been created.' ! write (*,*) ' ' ! end if ! elseif(itblout.eq.0)then ! continue ! else ! write (*,*) ' error from elimin: ', ! @ ' itblout should be 1 or 0; itblout= ',itblout ! stop ! end if return 112 format(10x,10(i3,9x)) 113 format(1x,i3,2x,9(1pe9.2,2x)) 114 format(1x,a3, 11(8x,i3)) 115 format( 1x,i3, 2x, 11(1pe10.3)) 116 format( 1x,a17,a2,a18,a2,a1 ) 117 format( 3x,i3, 4(8x,1pe10.3) ) 118 format( 3x,i3, 2(8x,1pe10.3) ) 120 format( 1x,i3, 1x,i3, 2x, 11(1pe10.3)) 132 format(i3) ! cambio: los formatos 114, 115 , 117 y 118 ! cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3 ! y ahora en vez de 11 capas de 5 en 5, hay 28 ! end c************************************************************************** subroutine smooth_cf( c, cs, i, nl, w ) c hace un smoothing de c(i,*), de la contribucion de todas las capas c menos de la capa en cuestion, la i. c opcion w (width): el tamanho de la ventana del smoothing. c output values: cs c************************************************************************** implicit none integer j,np,i,nl,w real*8 c(nl), cs(nl) if(w.eq.0) then do j=1,nl cs(j)=c(j) end do elseif(w.eq.3) then ! write (*,*) 'smoothing w=3' do j=1,i-4 if(j.eq.1) then cs(j)=c(j) else cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1)) end if end do do j=i+4,nl-1 if(j.eq.nl) then cs(j)=c(j) else cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1)) end if end do elseif(w.eq.5) then ! type *,'smoothing w=5' do j=3,i-4 if(j.eq.1) then cs(j)=c(j) else cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2)) end if end do do j=i+4,nl-2 if(j.eq.nl) then cs(j)=c(j) else cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2)) end if end do end if return end c***************************************************************************** c suaviza c***************************************************************************** c subroutine suaviza ( x, n, ismooth, y ) c c x - input and return values c y - auxiliary vector c ismooth = 0 --> no smoothing is performed c ismooth = 1 --> weak smoothing (5 points, centred weighted) c ismooth = 2 --> normal smoothing (3 points, evenly weighted) c ismooth = 3 --> strong smoothing (5 points, evenly weighted) c malv august 1991 c***************************************************************************** implicit none integer n, imax, imin, i, ismooth real*8 x(n), y(n) c***************************************************************************** imin=1 imax=n if (ismooth.eq.0) then return elseif (ismooth.eq.1) then ! 5 points, with central weighting do i=imin,imax if(i.eq.imin)then y(i)=x(imin) elseif(i.eq.imax)then y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + & x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0 else y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0 end if end do elseif (ismooth.eq.2) then ! 3 points, evenly spaced do i=imin,imax if(i.eq.imin)then y(i)=x(imin) elseif(i.eq.imax)then y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0 else y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 end if end do elseif (ismooth.eq.3) then ! 5 points, evenly spaced do i=imin,imax if(i.eq.imin)then y(i) = x(imin) elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 elseif(i.eq.imax)then y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0 else y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0 end if end do else write (*,*) ' Error in suaviza.f Wrong ismooth value.' stop endif c rehago el cambio, para devolver x(i) do i=imin,imax x(i)=y(i) end do return end c***************************************************************************** c LUdec.F (includes lubksb_dp and ludcmp_dp subroutines) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Solution of linear equation without inverting matrix c using LU decomposition: c AA * xx = bb AA, bb: known c xx: to be found c AA and bb are not modified in this subroutine c c MALV , Sep 2007 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine LUdec(xx,aa,bb,m,n) implicit none ! Arguments integer,intent(in) :: m, n real*8,intent(in) :: aa(m,m), bb(m) real*8,intent(out) :: xx(m) ! Local variables real*8 a(n,n), b(n), x(n), d integer i, j, indx(n) ! Subrutinas utilizadas ! ludcmp_dp, lubksb_dp !!!!!!!!!!!!!!! Comienza el programa !!!!!!!!!!!!!! do i=1,n b(i) = bb(i+1) do j=1,n a(i,j) = aa(i+1,j+1) enddo enddo ! Descomposicion de auxm1 call ludcmp_dp ( a, n, n, indx, d) ! Sustituciones foward y backwards para hallar la solucion do i=1,n x(i) = b(i) enddo call lubksb_dp( a, n, n, indx, x ) do i=1,n xx(i+1) = x(i) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine ludcmp_dp(a,n,np,indx,d) c jul 2011 malv+fgg implicit none integer,intent(in) :: n, np real*8,intent(inout) :: a(np,np) real*8,intent(out) :: d integer,intent(out) :: indx(n) integer i, j, k, imax real*8,parameter :: tiny=1.0d-20 real*8 vv(n), aamax, sum, dum d=1.0d0 do 12 i=1,n aamax=0.0d0 do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.0) then write(*,*) 'ludcmp_dp: singular matrix!' stop endif vv(i)=1.0d0/aamax 12 continue do 19 j=1,n if (j.gt.1) then do 14 i=1,j-1 sum=a(i,j) if (i.gt.1)then do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum endif 14 continue endif aamax=0.0d0 do 16 i=j,n sum=a(i,j) if (j.gt.1)then do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum endif dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(j.ne.n)then if(a(j,j).eq.0.0)a(j,j)=tiny dum=1.0d0/a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue if(a(n,n).eq.0.0)a(n,n)=tiny return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine lubksb_dp(a,n,np,indx,b) c jul 2011 malv+fgg implicit none integer,intent(in) :: n,np real*8,intent(in) :: a(np,np) integer,intent(in) :: indx(n) real*8,intent(out) :: b(n) real*8 sum integer ii, ll, i, j ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0)then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.0) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) if(i.lt.n)then do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue endif b(i)=sum/a(i,i) 14 continue return end c***************************************************************************** c intersp c *********************************************************************** subroutine intersp(yy,zz,m,y,z,n,opt) c interpolation soubroutine. input values: y(n) at z(n). c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic c jul 2011 malv+fgg c *********************************************************************** implicit none integer n,m,i,j,opt real zz(m),yy(m),z(n),y(n) real zmin,zzmin,zmax,zzmax ! write(*,*) ' interpolating' ! call minsp(z,n,zmin) zmin=minval(z) ! call minsp(zz,m,zzmin) zzmin=minval(zz) ! call maxsp(z,n,zmax) zmax=maxval(z) ! call maxsp(zz,m,zzmax) zzmax=maxval(zz) if(zzmin.lt.zmin)then write(*,*) 'from interp: new variable out of limits' write(*,*) zzmin,'must be .ge. ',zmin stop ! elseif(zzmax.gt.zmax)then ! write(*,*)'from interp: new variable out of limits' ! write(*,*)zzmax, 'must be .le. ',zmax ! stop end if do 1,i=1,m do 2,j=1,n-1 if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 2 continue c in this case (zz(m).ge.z(n)) and j leaves the loop with j=n-1+1=n if(opt.eq.1)then yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1)) elseif(opt.eq.2)then if(y(n).eq.0.0.or.y(n-1).eq.0.0)then yy(i)=0.0 else yy(i)=exp(log(y(n-1))+log(y(n)/y(n-1))* @ (zz(i)-z(n-1))/(z(n)-z(n-1))) end if else write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt end if goto 1 3 continue if(opt.eq.1)then yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) elseif(opt.eq.2)then if(y(j+1).eq.0.0.or.y(j).eq.0.0)then yy(i)=0.0 else yy(i)=exp(log(y(j))+log(y(j+1)/y(j))* @ (zz(i)-z(j))/(z(j+1)-z(j))) end if else write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt end if 1 continue return end c***************************************************************************** c interdp c *********************************************************************** subroutine interdp(yy,zz,m,y,z,n,opt) c interpolation soubroutine. input values: y(n) at z(n). c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic c jul 2011: malv+fgg Adapted to LMD-MGCM c *********************************************************************** implicit none integer n,m,i,j,opt real*8 zz(m),yy(m),z(n),y(n), zmin,zzmin,zmax,zzmax ! write (*,*) ' d interpolating ' ! call mindp (z,n,zmin) zmin=minval(z) ! call mindp (zz,m,zzmin) zzmin=minval(zz) ! call maxdp (z,n,zmax) zmax=maxval(z) ! call maxdp (zz,m,zzmax) zzmax=maxval(zz) if(zzmin.lt.zmin)then write (*,*) 'from d interp: new variable out of limits' write (*,*) zzmin,'must be .ge. ',zmin stop ! elseif(zzmax.gt.zmax)then ! write (*,*) 'from interp: new variable out of limits' ! write (*,*) zzmax, 'must be .le. ',zmax ! stop end if do 1,i=1,m do 2,j=1,n-1 if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 2 continue c in this case (zz(m).eq.z(n)) and j leaves the loop with j=n-1+1=n if(opt.eq.1)then yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1)) elseif(opt.eq.2)then if(y(n).eq.0.0d0.or.y(n-1).eq.0.0d0)then yy(i)=0.0d0 else yy(i)=dexp(dlog(y(n-1))+dlog(y(n)/y(n-1))* @ (zz(i)-z(n-1))/(z(n)-z(n-1))) end if else write (*,*) @ ' from d interp error: opt must be 1 or 2, opt= ',opt end if goto 1 3 continue if(opt.eq.1)then yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) ! write (*,*) ' ' ! write (*,*) ' z(j),z(j+1) =', z(j),z(j+1) ! write (*,*) ' t(j),t(j+1) =', y(j),y(j+1) ! write (*,*) ' zz, tt = ', zz(i), yy(i) elseif(opt.eq.2)then if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then yy(i)=0.0d0 else yy(i)=dexp(dlog(y(j))+dlog(y(j+1)/y(j))* @ (zz(i)-z(j))/(z(j+1)-z(j))) end if else write (*,*) ' from interp error: opt must be 1 or 2, opt= ', @ opt end if 1 continue return end c***************************************************************************** c interdp_limits.F c *********************************************************************** subroutine interdp_limits ( yy,zz,m, i1,i2, y,z,n, j1,j2, opt) c Interpolation soubroutine. c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n c Input values: y(n) , z(n) (solo se usan los valores entre j1,j2) c zz(m) (solo se necesita entre i1,i2) c Output values: yy(m) (solo se calculan entre i1,i2) c Options: opt=1 -> lineal ,, opt=2 -> logarithmic c Difference with interdp: c here interpolation proceeds between indexes i1,i2 only c if i1=1 & i2=m, both subroutines are exactly the same c thus previous calls to interdp or interdp2 could be easily replaced c JAN 98 MALV Version for mz1d c jul 2011 malv+fgg Adapted to LMD-MGCM c *********************************************************************** implicit none ! Arguments integer n,m ! I. Dimensions integer i1, i2, j1, j2, opt ! I real*8 zz(m),yy(m) ! O real*8 z(n),y(n) ! I ! Local variables integer i,j real*8 zmin,zzmin,zmax,zzmax c ******************************* ! type *, ' d interpolating ' ! call mindp_limits (z,n,zmin, j1,j2) zmin=minval(z(j1:j2)) ! call mindp_limits (zz,m,zzmin, i1,i2) zzmin=minval(zz(i1:i2)) ! call maxdp_limits (z,n,zmax, j1,j2) zmax=maxval(z(j1:j2)) ! call maxdp_limits (zz,m,zzmax, i1,i2) zzmax=maxval(zz(i1:i2)) if(zzmin.lt.zmin)then write (*,*) 'from d interp: new variable out of limits' write (*,*) zzmin,'must be .ge. ',zmin stop ! elseif(zzmax.gt.zmax)then ! type *,'from interp: new variable out of limits' ! type *,zzmax, 'must be .le. ',zmax ! stop end if do 1,i=i1,i2 do 2,j=j1,j2-1 if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3 2 continue c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2 if(opt.eq.1)then yy(i)=y(j2-1)+(y(j2)-y(j2-1))* $ (zz(i)-z(j2-1))/(z(j2)-z(j2-1)) elseif(opt.eq.2)then if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then yy(i)=0.0d0 else yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))* @ (zz(i)-z(j2-1))/(z(j2)-z(j2-1))) end if else write (*,*) ' d interp : opt must be 1 or 2, opt= ',opt end if goto 1 3 continue if(opt.eq.1)then yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j)) ! type *, ' ' ! type *, ' z(j),z(j+1) =', z(j),z(j+1) ! type *, ' t(j),t(j+1) =', y(j),y(j+1) ! type *, ' zz, tt = ', zz(i), yy(i) elseif(opt.eq.2)then if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then yy(i)=0.0d0 else yy(i)=exp(log(y(j))+log(y(j+1)/y(j))* @ (zz(i)-z(j))/(z(j+1)-z(j))) end if else write (*,*) ' interp : opt must be 1 or 2, opt= ',opt end if 1 continue return end c***************************************************************************** c Subroutines previously included in tcrco2_subrut.F c*********************************************************************** c tcrco2_subrut.f c c jan 98 malv version for mz1d. copied from solar10/mz4sub.f c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** ************************************************************************ subroutine dinterconnection ( v, vt ) * input: vib. temp. from che*.for programs, vt(nl) * output: test vibrational temp. for other che*.for, v(nl) ! iconex_smooth=1 ==> with smoothing ! iconex_smooth=0 ==> without smoothing ! iconex_tk=40 ==> with forced lte up to 40 km ! iconex_tk=20 ==> with forced lte up to 20 km ************************************************************************ implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c argumentos real*8 vt(nl), v(nl) c local variables integer i c ************* do i=1,nl v(i) = vt(i) end do ! lo siguiente se utilizaba en solar10, pero es mejor introducirlo en ! la driver. por ahora no lo uso todavia. ! call fluctua(v,iconex_fluctua) ! call smooth_nl(v,iconex_smooth,nl) ! call forzar_tk(v,iconex_tk) return end c*********************************************************************** function planckdp(tp,xnu) c returns the black body function at wavenumber xnu and temperature t. c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' ! common/datis/ pi, vlight, ee, hplanck, gamma, ab, ! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg ! real*8 pi, vlight, ee, hplanck, gamma, ab, ! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg real*8 planckdp real*8 xnu real tp planckdp = gamma*xnu**3.0 / exp( ee*xnu/dble(tp) ) !erg cm-2.sr-1/cm-1. return end c **************************************************************** function bandid (ib) c returns the 2 character code of the band c **************************************************************** implicit none integer ib character*2 bandid 132 format(i2) ! encode (2,132,bandid) ib write ( bandid, 132) ib if ( ib .eq. 1 ) bandid = '01' if ( ib .eq. 2 ) bandid = '02' if ( ib .eq. 3 ) bandid = '03' if ( ib .eq. 4 ) bandid = '04' if ( ib .eq. 5 ) bandid = '05' if ( ib .eq. 6 ) bandid = '06' if ( ib .eq. 7 ) bandid = '07' if ( ib .eq. 8 ) bandid = '08' if ( ib .eq. 9 ) bandid = '09' if ( ib .eq. 0 ) bandid = '00' c end return end c***************************************************************************** c Subroutines previously included in mat_oper.F c***************************************************************************** c set of subroutines for the cz*.for programs: ! subroutine unit(a,n) ! subroutine diago(a,v,n) diagonal matrix with v ! subroutine invdiag(a,b,n) inverse of diagonal matrix ! subroutine sypvvv(a,b,c,d,n) suma y prod de 3 vectores, muy comun ! subroutine sypvmv(v,w,b,u,n) suma y prod de 3 vectores, muy comun ! subroutine mulmvv(w,b,u,v,n) prod matriz vector vector ! subroutine muymvv(w,b,u,v,n) prod matriz (inv.vector) vector ! subroutine samem (a,m,n) ! subroutine mulmv(a,b,c,n) ! subroutine mulmm(a,b,c,n) ! subroutine resmm(a,b,c,n) ! subroutine mulvv(a,b,c,n) ! subroutine sumvv(a,b,c,n) ! subroutine zerom(a,n) ! subroutine zero4m(a,b,c,d,n) ! subroutine zero3m(a,b,c,n) ! subroutine zero2m(a,b,n) ! subroutine zerov(a,n) ! subroutine zero4v(a,b,c,d,n) ! subroutine zero3v(a,b,c,n) ! subroutine zero2v(a,b,n) ! ! ! May-05 Sustituimos todos los zerojt de cristina por las subrutinas ! genericas zerov*** ! c *********************************************************************** subroutine unit(a,n) c store the unit value in the diagonal of a c *********************************************************************** real*8 a(n,n) integer n,i,j,k do 1,i=2,n-1 do 2,j=2,n-1 if(i.eq.j) then a(i,j) = 1.d0 else a(i,j)=0.0d0 end if 2 continue 1 continue do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c *********************************************************************** subroutine diago(a,v,n) c store the vector v in the diagonal elements of the square matrix a c *********************************************************************** implicit none integer n,i,j,k real*8 a(n,n),v(n) do 1,i=2,n-1 do 2,j=2,n-1 if(i.eq.j) then a(i,j) = v(i) else a(i,j)=0.0d0 end if 2 continue 1 continue do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c *********************************************************************** subroutine samem (a,m,n) c store the matrix m in the matrix a c *********************************************************************** real*8 a(n,n),m(n,n) integer n,i,j,k do 1,i=2,n-1 do 2,j=2,n-1 a(i,j) = m(i,j) 2 continue 1 continue do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c *********************************************************************** subroutine mulmv(a,b,c,n) c do a(i)=b(i,j)*c(j). a, b, and c must be distint c *********************************************************************** implicit none integer n,i,j real*8 a(n),b(n,n),c(n),sum do 1,i=2,n-1 sum=0.0d0 do 2,j=2,n-1 sum=sum+ (b(i,j)) * (c(j)) 2 continue a(i)=sum 1 continue a(1) = 0.0d0 a(n) = 0.0d0 return end c *********************************************************************** subroutine mulmm(a,b,c,n) c *********************************************************************** real*8 a(n,n), b(n,n), c(n,n) integer n,i,j,k ! do i=2,n-1 ! do j=2,n-1 ! a(i,j)= 0.d00 ! do k=2,n-1 ! a(i,j) = a(i,j) + b(i,k) * c(k,j) ! end do ! end do ! end do do j=2,n-1 do i=2,n-1 a(i,j)=0.d0 enddo do k=2,n-1 do i=2,n-1 a(i,j)=a(i,j)+b(i,k)*c(k,j) enddo enddo enddo do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c *********************************************************************** subroutine resmm(a,b,c,n) c *********************************************************************** real*8 a(n,n), b(n,n), c(n,n) integer n,i,j,k do i=2,n-1 do j=2,n-1 a(i,j)= b(i,j) - c(i,j) end do end do do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c *********************************************************************** subroutine sumvv(a,b,c,n) c a(i)=b(i)+c(i) c *********************************************************************** implicit none integer n,i real*8 a(n),b(n),c(n) do 1,i=2,n-1 a(i)= (b(i)) + (c(i)) 1 continue a(1) = 0.0d0 a(n) = 0.0d0 return end c *********************************************************************** subroutine zerom(a,n) c a(i,j)= 0.0 c *********************************************************************** implicit none integer n,i,j real*8 a(n,n) do 1,i=1,n do 2,j=1,n a(i,j) = 0.0d0 2 continue 1 continue return end c *********************************************************************** subroutine zero4m(a,b,c,d,n) c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 c *********************************************************************** real*8 a(n,n), b(n,n), c(n,n), d(n,n) integer n,i,j do 1,i=1,n do 2,j=1,n a(i,j) = 0.0d0 b(i,j) = 0.0d0 c(i,j) = 0.0d0 d(i,j) = 0.0d0 2 continue 1 continue return end c *********************************************************************** subroutine zero3m(a,b,c,n) c a(i,j) = b(i,j) = c(i,j) = 0.0 c ********************************************************************** real*8 a(n,n), b(n,n), c(n,n) integer n,i,j do 1,i=1,n do 2,j=1,n a(i,j) = 0.0d0 b(i,j) = 0.0d0 c(i,j) = 0.0d0 2 continue 1 continue return end c *********************************************************************** subroutine zero2m(a,b,n) c a(i,j) = b(i,j) = 0.0 c *********************************************************************** real*8 a(n,n), b(n,n) integer n,i,j do 1,i=1,n do 2,j=1,n a(i,j) = 0.0d0 b(i,j) = 0.0d0 2 continue 1 continue return end c *********************************************************************** subroutine zerov(a,n) c a(i)= 0.0 c *********************************************************************** real*8 a(n) integer n,i do 1,i=1,n a(i) = 0.0d0 1 continue return end c *********************************************************************** subroutine zero4v(a,b,c,d,n) c a(i) = b(i) = c(i) = d(i,j) = 0.0 c *********************************************************************** real*8 a(n), b(n), c(n), d(n) integer n,i do 1,i=1,n a(i) = 0.0d0 b(i) = 0.0d0 c(i) = 0.0d0 d(i) = 0.0d0 1 continue return end c *********************************************************************** subroutine zero3v(a,b,c,n) c a(i) = b(i) = c(i) = 0.0 c *********************************************************************** real*8 a(n), b(n), c(n) integer n,i do 1,i=1,n a(i) = 0.0d0 b(i) = 0.0d0 c(i) = 0.0d0 1 continue return end c *********************************************************************** subroutine zero2v(a,b,n) c a(i) = b(i) = 0.0 c *********************************************************************** real*8 a(n), b(n) integer n,i do 1,i=1,n a(i) = 0.0d0 b(i) = 0.0d0 1 continue return end