c*********************************************************************** c mzescape.f c*********************************************************************** c c program for calculating atmospheric escape functions, from c a calculation of transmittances and derivatives of these ones subroutine mzescape(ig,taustar,tauinf,tauii, ib,isot, iirw,iimu) c jul 2011 malv+fgg adapted to LMD-MGCM c nov 99 malv adapt mztf to compute taustar (pg.23b-ma c nov 98 malv allow for overlaping in the lorentz line c jan 98 malv version for mz1d. based on curtis/mztf.for c 17-jul-96 mlp&crs change the calculation of mr. c evitar: divide por cero. anhadiendo: ff c oct-92 malv correct s(t) dependence for all histogr bands c june-92 malv proper lower levels for laser bands c may-92 malv new temperature dependence for laser bands c @ 991 malv boxing for the averaged absorber amount and t c ? malv extension up to 200 km altitude in mars c 13-nov-86 mlp include the temperature weighted to match c the eqw in the strong doppler limit. c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig ! ADDED FOR TRACEBACK real*8 taustar(nl) ! o real*8 tauinf(nl) ! o real*8 tauii(nl) ! o integer ib ! i integer isot ! i integer iirw ! i integer iimu ! i c local variables and constants integer i, in, ir, im, k,j integer nmu parameter (nmu = 8) ! real*8 tauinf(nl) real*8 con(nzy), coninf real*8 c1, c2, ccc real*8 t1, t2 real*8 p1, p2 real*8 mr1, mr2 real*8 st1, st2 real*8 c1box(70), c2box(70) real*8 ff ! to avoid too small numbers real*8 tvtbs(nzy) real*8 st, beta, ts, eqwmu real*8 mu(nmu), amu(nmu) real*8 zld(nl), zyd(nzy) real*8 correc real deltanux ! width of vib-rot band (cm-1) character isotcode*2 real*8 maximum real*8 csL, psL, Desp, wsL ! for Strong Lorentz limit c formats 111 format(a1) 112 format(a2) 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** c some needed values ! rl=sqrt(log(2.d0)) ! pi2 = 3.14159265358989d0 beta = 1.8d0 ! imrco = 0.9865 c esto es para que las subroutines de mztfsub calculen we c de la forma apropiada para mztf, no para fot icls=icls_mztf c codigos para filenames ! if (isot .eq. 1) isotcode = '26' ! if (isot .eq. 2) isotcode = '28' ! if (isot .eq. 3) isotcode = '36' ! if (isot .eq. 4) isotcode = '27' ! if (isot .eq. 5) isotcode = '62' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then !! encode(2,101,ibcode1)ib ! write ( ibcode1, 101) ib ! else !! encode(2,202,ibcode2)ib ! write (ibcode2, 202) ib ! endif ! write (*,'( 30h calculating curtis matrix : ,2x, ! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot c integration in angle !!!!!!!!!!!!!!!!!!!! c------- diffusivity approx. if (iimu.eq.1) then ! write (*,*) ' diffusivity approx. beta = ',beta mu(1) = 1.0d0 amu(1)= 1.0d0 c-------data for 8 points integration elseif (iimu.eq.4) then write (*,*)' 4 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.339981043584856)/2.0d0 mu(2)=(1.0d0-0.339981043584856)/2.0d0 mu(3)=(1.0d0+0.861136311594053)/2.0d0 mu(4)=(1.0d0-0.861136311594053)/2.0d0 amu(1)=0.652145154862546 amu(2)=amu(1) amu(3)=0.347854845137454 amu(4)=amu(3) beta=1.0d0 c-------data for 8 points integration elseif(iimu.eq.8) then write (*,*)' 8 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.183434642495650)/2.0d0 mu(2)=(1.0d0-0.183434642495650)/2.0d0 mu(3)=(1.0d0+0.525532409916329)/2.0d0 mu(4)=(1.0d0-0.525532409916329)/2.0d0 mu(5)=(1.0d0+0.796666477413627)/2.0d0 mu(6)=(1.0d0-0.796666477413627)/2.0d0 mu(7)=(1.0d0+0.960289856497536)/2.0d0 mu(8)=(1.0d0-0.960289856497536)/2.0d0 amu(1)=0.362683783378362 amu(2)=amu(1) amu(3)=0.313706645877887 amu(4)=amu(3) amu(5)=0.222381034453374 amu(6)=amu(5) amu(7)=0.101228536290376 amu(8)=amu(7) beta=1.0d0 end if c!!!!!!!!!!!!!!!!!!!!!!! ccc ccc determine abundances included in the absorber amount ccc c first, set up the grid ready for interpolation. do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble(zl(i)) enddo c vibr. temp of the bending mode : if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state) c por similitud a la que se hace en cza.for do i=1,nzy if (isot.eq.5) then con(i) = dble( coy(i) * imrco ) else con(i) = dble( co2y(i) * imr(isot) ) correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) ) con(i) = con(i) * ( 1.d0 - correc ) endif c----------------------------------------------------------------------- c mlp & cristina. 17 july 1996 c change the calculation of mr. it is used for calculating partial press c alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin c collisions with other co2 isotopes (including the major one, 626) c as if they were with n2. assuming mr as co2/nt, we consider collisions c of type 628-626 as of 626-626 instead of as 626-n2. c mrx(i)=con(i)/ntx(i) ! old malv ! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs c jan 98: c esta modif de mlp implica anular el correc (deberia revisar esto) mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98 c----------------------------------------------------------------------- end do ! como beta y 1.d5 son comunes a todas las weighted absorber amounts, ! los simplificamos: ! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) ) coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) ) ! write (*,*) ' coninf =', coninf ccc ccc temp dependence of the band strength and ccc nlte correction factor for the absorber amount ccc call mztf_correccion ( coninf, con, ib, isot, 0 ) ccc ccc reads histogrammed spectral data (strength for lte and vmr=1) ccc !hfile1 = dirspec//'hi'//dn ! Ya no distinguimos entre d/n !! hfile1 = dirspec//'hid' ! (see why in his.for) ! hfile1='hid' !! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' ! else ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' ! endif !write (*,*) ' /MZESCAPE/ hisfile: ', hisfile ! the argument to rhist is to make this compatible with mztf_comp.f, ! which is a useful modification of mztf.f (to change strengths of bands ! call rhist (1.0) if(ib.eq.1) then if(isot.eq.1) then !Case 1 mm=mm_c1 nbox=nbox_c1 tmin=tmin_c1 tmax=tmax_c1 do i=1,nbox_max no(i)=no_c1(i) dist(i)=dist_c1(i) do j=1,nhist sk1(j,i)=sk1_c1(j,i) xls1(j,i)=xls1_c1(j,i) xln1(j,i)=xln1_c1(j,i) xld1(j,i)=xld1_c1(j,i) enddo enddo do j=1,nhist thist(j)=thist_c1(j) enddo else if(isot.eq.2) then !Case 2 mm=mm_c2 nbox=nbox_c2 tmin=tmin_c2 tmax=tmax_c2 do i=1,nbox_max no(i)=no_c2(i) dist(i)=dist_c2(i) do j=1,nhist sk1(j,i)=sk1_c2(j,i) xls1(j,i)=xls1_c2(j,i) xln1(j,i)=xln1_c2(j,i) xld1(j,i)=xld1_c2(j,i) enddo enddo do j=1,nhist thist(j)=thist_c2(j) enddo else if(isot.eq.3) then !Case 3 mm=mm_c3 nbox=nbox_c3 tmin=tmin_c3 tmax=tmax_c3 do i=1,nbox_max no(i)=no_c3(i) dist(i)=dist_c3(i) do j=1,nhist sk1(j,i)=sk1_c3(j,i) xls1(j,i)=xls1_c3(j,i) xln1(j,i)=xln1_c3(j,i) xld1(j,i)=xld1_c3(j,i) enddo enddo do j=1,nhist thist(j)=thist_c3(j) enddo else if(isot.eq.4) then !Case 4 mm=mm_c4 nbox=nbox_c4 tmin=tmin_c4 tmax=tmax_c4 do i=1,nbox_max no(i)=no_c4(i) dist(i)=dist_c4(i) do j=1,nhist sk1(j,i)=sk1_c4(j,i) xls1(j,i)=xls1_c4(j,i) xln1(j,i)=xln1_c4(j,i) xld1(j,i)=xld1_c4(j,i) enddo enddo do j=1,nhist thist(j)=thist_c4(j) enddo else write(*,*)'isot must be 2,3 or 4 for ib=1!!' write(*,*)'stop at mzescape/312' stop endif else if (ib.eq.2) then if(isot.eq.1) then !Case 5 mm=mm_c5 nbox=nbox_c5 tmin=tmin_c5 tmax=tmax_c5 do i=1,nbox_max no(i)=no_c5(i) dist(i)=dist_c5(i) do j=1,nhist sk1(j,i)=sk1_c5(j,i) xls1(j,i)=xls1_c5(j,i) xln1(j,i)=xln1_c5(j,i) xld1(j,i)=xld1_c5(j,i) enddo enddo do j=1,nhist thist(j)=thist_c5(j) enddo else write(*,*)'isot must be 1 for ib=2!!' write(*,*)'stop at mzescape/336' stop endif else if (ib.eq.3) then if(isot.eq.1) then !Case 6 mm=mm_c6 nbox=nbox_c6 tmin=tmin_c6 tmax=tmax_c6 do i=1,nbox_max no(i)=no_c6(i) dist(i)=dist_c6(i) do j=1,nhist sk1(j,i)=sk1_c6(j,i) xls1(j,i)=xls1_c6(j,i) xln1(j,i)=xln1_c6(j,i) xld1(j,i)=xld1_c6(j,i) enddo enddo do j=1,nhist thist(j)=thist_c6(j) enddo else write(*,*)'isot must be 1 for ib=3!!' write(*,*)'stop at mzescape/360' stop endif else if (ib.eq.4) then if(isot.eq.1) then !Case 7 mm=mm_c7 nbox=nbox_c7 tmin=tmin_c7 tmax=tmax_c7 do i=1,nbox_max no(i)=no_c7(i) dist(i)=dist_c7(i) do j=1,nhist sk1(j,i)=sk1_c7(j,i) xls1(j,i)=xls1_c7(j,i) xln1(j,i)=xln1_c7(j,i) xld1(j,i)=xld1_c7(j,i) enddo enddo do j=1,nhist thist(j)=thist_c7(j) enddo else write(*,*)'isot must be 1 for ib=4!!' write(*,*)'stop at mzescape/384' stop endif else write(*,*)'ib must be 1,2,3 or 4!!' write(*,*)'stop at mzescape/389' endif if (isot.ne.5) deltanux = deltanu(isot,ib) if (isot.eq.5) deltanux = deltanuco c****** c****** calculation of tauinf(nl) c****** call initial ff=1.0e10 do i=nl,1,-1 if(i.eq.nl)then call intz (zl(i),c2,p2,mr2,t2, con) do kr=1,nbox ta(kr)=t2 end do ! write (*,*) ' i, t2 =', i, t2 call interstrength (st2,t2,ka,ta) aa = p2 * coninf * mr2 * (st2 * ff) bb = p2 * coninf * st2 cc = coninf * st2 dd = t2 * coninf * st2 do kr=1,nbox ccbox(kr) = coninf * ka(kr) ddbox(kr) = t2 * ccbox(kr) ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) else call intz (zl(i),c1,p1,mr1,t1, con) do kr=1,nbox ta(kr)=t1 end do ! write (*,*) ' i, t1 =', i, t1 call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 ccc = ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + @ ( c1box(kr) + c2box(kr) )/2.d0 ddbox(kr) = ddbox(kr) + @ ( t1*c1box(kr)+t2*c2box(kr) )/2.d0 end do mr2 = mr1 c2=c1 do kr=1,nbox c2box(kr) = c1box(kr) end do t2=t1 p2=p1 end if pt = bb / cc pp = aa / (cc*ff) ! ta=dd/cc ! tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do ! write (*,*) ' i, ts =', i, ts call interstrength(st,ts,ka,ta) ! call intershape(alsa,alna,alda,tdop) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mzescape/480',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw (ig,iirw, 0, csL,psL, Desp, wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do ! tauinf(i) = exp( - eqwmu / dble(deltanux) ) tauinf(i) = 1.d0 - eqwmu / dble(deltanux) if (tauinf(i).lt.0.d0) tauinf(i) = 0.0d0 if (i.eq.nl) then taustar(i) = 0.0d0 else taustar(i) = dble(deltanux) * (tauinf(i+1)-tauinf(i)) ! ~ / ( beta * cc * 1.d5 ) ~ / ( beta * ccc * 1.d5 ) endif end do ! i continue c****** c****** calculation of tau(in,ir) for n<=r c****** do 1 in=1,nl-1 call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox c1box(kr) = c1 * ka(kr) * dble(deltaz) end do c1 = c1 * st1 * dble(deltaz) call intz (zl(in+1), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox c2box(kr) = c2 * ka(kr) * dble(deltaz) end do c2 = c2 * st2 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + (c1box(kr)+c2box(kr))/2.d0 ddbox(kr) = ddbox(kr) + (t1*c1box(kr)+t2*c2box(kr))/2.d0 end do mr1=mr2 t1=t2 c1=c2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength(st,ts,ka,ta) call intershape(alsa,alna,alda,ta) eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mzescape/566',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw (ig,iirw, 0, csL,psL, Desp, wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tauii(in) = exp( - eqwmu / dble(deltanux) ) !write (*,*) 'i,tauii=',in,tauii(in) 1 continue tauii(nl) = 1.0d0 c end return end c*********************************************************************** c mzescape_normaliz.f c*********************************************************************** c c program for correcting some strange values and for normalizing c the atmospheric escape functions computed by mzescape_15um.f c possibilities according to istyle (see mzescape_15um.f). c subroutine mzescape_normaliz ( taustar, istyle ) c dic 99 malv first version c jul 2011 malv+fgg Adapted to LMD-MGCM c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real*8 taustar(nl) ! o integer istyle ! i c local variables and constants integer i, imaximum real*8 maximum c*********************************************************************** ! ! correcting strange values at top, eliminating local maxima, etc... ! taustar(nl) = taustar(nl-1) if ( istyle .eq. 1 ) then imaximum = nl maximum = taustar(nl) do i=1,nl-1 if (taustar(i).gt.maximum) taustar(i) = taustar(nl) enddo elseif ( istyle .eq. 2 ) then imaximum = nl maximum = taustar(nl) do i=nl-1,1,-1 if (taustar(i).gt.maximum) then maximum = taustar(i) imaximum = i endif enddo do i=imaximum,nl if (taustar(i).lt.maximum) taustar(i) = maximum enddo endif ! ! normalizing ! do i=1,nl taustar(i) = taustar(i) / maximum enddo c end return end c*********************************************************************** c mzescape_fb.f c*********************************************************************** subroutine mzescape_fb(ig) c computes the escape functions of the most important 15um bands c this calls mzescape ( taustar,tauinf,tauii, ib,isot, iirw,iimu c nov 99 malv based on cm15um_fb.f c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c local variables integer i, ib, ik, istyle integer ig !ADDED FOR TRACEBACK real*8 tau_factor real*8 aux(nl), aux2(nl), aux3(nl) c*********************************************************************** call mzescape (ig,taustar21,tauinf210,tauii210,1,2,irw_mztf,imu) call mzescape (ig,taustar31,tauinf310,tauii310,1,3,irw_mztf,imu) call mzescape (ig,taustar41,tauinf410,tauii410,1,4,irw_mztf,imu) istyle = 2 call mzescape_normaliz ( taustar21, istyle ) call mzescape_normaliz ( taustar31, istyle ) call mzescape_normaliz ( taustar41, istyle ) c end return end c*********************************************************************** c mzescape_fh.f c*********************************************************************** subroutine mzescape_fh(ig) c jul 2011 malv+fgg c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c local variables integer i, ib, ik, istyle integer ig ! ADDED FOR TRACEBACK real*8 tau_factor real*8 aux(nl), aux2(nl), aux3(nl) c*********************************************************************** call zero4v( aux, taustar12,tauinf121,tauii121, nl) do ik=1,3 ib=ik+1 call mzescape ( ig,aux,aux2,aux3, ib, 1,irw_mztf,imu ) tau_factor = 1.d0 if (ik.eq.1) tau_factor = dble(667.75/618.03) if (ik.eq.3) tau_factor = dble(667.75/720.806) do i=1,nl taustar12(i) = taustar12(i) + aux(i) * tau_factor tauinf121(i) = tauinf121(i) + aux2(i) * tau_factor tauii121(i) = tauii121(i) + aux3(i) * tau_factor enddo enddo istyle = 2 call mzescape_normaliz ( taustar12, istyle ) c end return end c*********************************************************************** c mztud.f c*********************************************************************** subroutine mztud ( ig,cf,cfup,cfdw,vc,taugr, ib,isot, @ iirw,iimu,itauout,icfout,itableout ) c program for calculating atmospheric transmittances c to be used in the calculation of curtis matrix coefficients c i*out = 1 output of data c i*out = 0 no output c itableout = 30 output de toda la C.M. y el VC y las poblaciones de los c estados 626(020), esta opcion nueva se añade porque c itableout=1 saca o bien solamente de 5 en 5 capas c o bien los elementos de C.M. desde una cierta capa c (consultese elimin_mz1d.f que es quien lo hace); lo c de las poblaciones (020) lo hace mztf_correcion.f c jul 2011 malv+fgg Adapted to LMD-MGCM c jan 07 malv Add new vertical fine grid zy, similar to zx c sep-oct 01 malv update for fluxes for hb and fb, adapt to Linux c nov 98 mavl allow for overlaping in the lorentz line c jan 98 malv version for mz1d. based on curtis/mztf.for c 17-jul-96 mlp&crs change the calculation of mr. c evitar: divide por cero. anhadiendo: ff c oct-92 malv correct s(t) dependence for all histogr bands c june-92 malv proper lower levels for laser bands c may-92 malv new temperature dependence for laser bands c @ 991 malv boxing for the averaged absorber amount and t c ? malv extension up to 200 km altitude in mars c 13-nov-86 mlp include the temperature weighted to match c the eqw in the strong doppler limit. c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig !ADDED FOR TRACEBACK real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o real*8 vc(nl), taugr(nl) ! o integer ib ! i integer isot ! i integer iirw ! i integer iimu ! i integer itauout ! i integer icfout ! i integer itableout ! i c local variables and constants integer i, in, ir, im, k,j integer nmu parameter (nmu = 8) real*8 tau(nl,nl) real*8 tauinf(nl) real*8 con(nzy), coninf real*8 c1, c2 real*8 t1, t2 real*8 p1, p2 real*8 mr1, mr2 real*8 st1, st2 real*8 c1box(70), c2box(70) real*8 ff ! to avoid too small numbers real*8 tvtbs(nzy) real*8 st, beta, ts, eqwmu real*8 mu(nmu), amu(nmu) real*8 zld(nl), zyd(nzy) real*8 correc real deltanux ! width of vib-rot band (cm-1) character isotcode*2 integer idummy real*8 Desp,wsL c formats 111 format(a1) 112 format(a2) 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** c some needed values ! rl=sqrt(log(2.d0)) ! pi2 = 3.14159265358989d0 beta = 1.8d0 ! beta = 1.0d0 idummy = 0 Desp = 0.0d0 wsL = 0.0d0 ! write (*,*) ' MZTUD/ iirw = ', iirw c esto es para que las subroutines de mztfsub calculen we c de la forma apropiada para mztf, no para fot icls=icls_mztf c codigos para filenames ! if (isot .eq. 1) isotcode = '26' ! if (isot .eq. 2) isotcode = '28' ! if (isot .eq. 3) isotcode = '36' ! if (isot .eq. 4) isotcode = '27' ! if (isot .eq. 5) isotcode = '62' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (ibcode1,101) ib ! else ! write (ibcode2,202) ib ! endif ! write (*,'( 30h calculating curtis matrix : ,2x, ! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot c integration in angle !!!!!!!!!!!!!!!!!!!! c------- diffusivity approx. if (iimu.eq.1) then ! write (*,*) ' diffusivity approx. beta = ',beta mu(1) = 1.0d0 amu(1)= 1.0d0 c-------data for 8 points integration elseif (iimu.eq.4) then write (*,*)' 4 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.339981043584856)/2.0d0 mu(2)=(1.0d0-0.339981043584856)/2.0d0 mu(3)=(1.0d0+0.861136311594053)/2.0d0 mu(4)=(1.0d0-0.861136311594053)/2.0d0 amu(1)=0.652145154862546 amu(2)=amu(1) amu(3)=0.347854845137454 amu(4)=amu(3) beta=1.0d0 c-------data for 8 points integration elseif(iimu.eq.8) then write (*,*)' 8 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.183434642495650)/2.0d0 mu(2)=(1.0d0-0.183434642495650)/2.0d0 mu(3)=(1.0d0+0.525532409916329)/2.0d0 mu(4)=(1.0d0-0.525532409916329)/2.0d0 mu(5)=(1.0d0+0.796666477413627)/2.0d0 mu(6)=(1.0d0-0.796666477413627)/2.0d0 mu(7)=(1.0d0+0.960289856497536)/2.0d0 mu(8)=(1.0d0-0.960289856497536)/2.0d0 amu(1)=0.362683783378362 amu(2)=amu(1) amu(3)=0.313706645877887 amu(4)=amu(3) amu(5)=0.222381034453374 amu(6)=amu(5) amu(7)=0.101228536290376 amu(8)=amu(7) beta=1.0d0 end if c!!!!!!!!!!!!!!!!!!!!!!! ccc ccc determine abundances included in the absorber amount ccc c first, set up the grid ready for interpolation. do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble(zl(i)) enddo c vibr. temp of the bending mode : if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) !if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 ) c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state) c por similitud a la que se hace en cza.for ; esto solo se hace para CO2 !write (*,*) 'imr(isot) = ', isot, imr(isot) do i=1,nzy if (isot.eq.5) then con(i) = dble( coy(i) * imrco ) else con(i) = dble( co2y(i) * imr(isot) ) correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) ) con(i) = con(i) * ( 1.d0 - correc ) ! write (*,*) ' iz, correc, co2y(i), con(i) =', ! @ i,correc,co2y(i),con(i) endif !----------------------------------------------------------------- ! mlp & cristina. 17 july 1996 change the calculation of mr. ! it is used for calculating partial press ! alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) ! for an isotope, if mr is obtained by ! co2*imr(iso)/nt ! we are considerin collisions with other co2 isotopes ! (including the major one, 626) as if they were with n2. ! assuming mr as co2/nt, we consider collisions ! of type 628-626 as of 626-626 instead of as 626-n2. ! mrx(i)=con(i)/ntx(i) ! old malv ! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs ! jan 98: ! esta modif de mlp implica anular el correc (deberia revisar esto) mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98 !----------------------------------------------------------------- end do ! como beta y 1.d5 son comunes a todas las weighted absorber amounts, ! los simplificamos: ! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) ) !write (*,*) ' con(nz), con(nz-1) =', con(nz), con(nz-1) coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) ) !write (*,*) ' coninf =', coninf ccc ccc temp dependence of the band strength and ccc nlte correction factor for the absorber amount ccc call mztf_correccion ( coninf, con, ib, isot, itableout ) ccc ccc reads histogrammed spectral data (strength for lte and vmr=1) ccc !hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto ! hfile1 = dirspec//'hid' !(see why in his.for) ! hfile1='hid' !! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' ! else ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' ! endif if(ib.eq.1) then if(isot.eq.1) then !Case 1 mm=mm_c1 nbox=nbox_c1 tmin=tmin_c1 tmax=tmax_c1 do i=1,nbox_max no(i)=no_c1(i) dist(i)=dist_c1(i) do j=1,nhist sk1(j,i)=sk1_c1(j,i) xls1(j,i)=xls1_c1(j,i) xln1(j,i)=xln1_c1(j,i) xld1(j,i)=xld1_c1(j,i) enddo enddo do j=1,nhist thist(j)=thist_c1(j) enddo else if(isot.eq.2) then !Case 2 mm=mm_c2 nbox=nbox_c2 tmin=tmin_c2 tmax=tmax_c2 do i=1,nbox_max no(i)=no_c2(i) dist(i)=dist_c2(i) do j=1,nhist sk1(j,i)=sk1_c2(j,i) xls1(j,i)=xls1_c2(j,i) xln1(j,i)=xln1_c2(j,i) xld1(j,i)=xld1_c2(j,i) enddo enddo do j=1,nhist thist(j)=thist_c2(j) enddo else if(isot.eq.3) then !Case 3 mm=mm_c3 nbox=nbox_c3 tmin=tmin_c3 tmax=tmax_c3 do i=1,nbox_max no(i)=no_c3(i) dist(i)=dist_c3(i) do j=1,nhist sk1(j,i)=sk1_c3(j,i) xls1(j,i)=xls1_c3(j,i) xln1(j,i)=xln1_c3(j,i) xld1(j,i)=xld1_c3(j,i) enddo enddo do j=1,nhist thist(j)=thist_c3(j) enddo else if(isot.eq.4) then !Case 4 mm=mm_c4 nbox=nbox_c4 tmin=tmin_c4 tmax=tmax_c4 do i=1,nbox_max no(i)=no_c4(i) dist(i)=dist_c4(i) do j=1,nhist sk1(j,i)=sk1_c4(j,i) xls1(j,i)=xls1_c4(j,i) xln1(j,i)=xln1_c4(j,i) xld1(j,i)=xld1_c4(j,i) enddo enddo do j=1,nhist thist(j)=thist_c4(j) enddo else write(*,*)'isot must be 2,3 or 4 for ib=1!!' write(*,*)'stop at mztud/324' stop endif else if (ib.eq.2) then if(isot.eq.1) then !Case 5 mm=mm_c5 nbox=nbox_c5 tmin=tmin_c5 tmax=tmax_c5 do i=1,nbox_max no(i)=no_c5(i) dist(i)=dist_c5(i) do j=1,nhist sk1(j,i)=sk1_c5(j,i) xls1(j,i)=xls1_c5(j,i) xln1(j,i)=xln1_c5(j,i) xld1(j,i)=xld1_c5(j,i) enddo enddo do j=1,nhist thist(j)=thist_c5(j) enddo else write(*,*)'isot must be 1 for ib=2!!' write(*,*)'stop at mztud/348' stop endif else if (ib.eq.3) then if(isot.eq.1) then !Case 6 mm=mm_c6 nbox=nbox_c6 tmin=tmin_c6 tmax=tmax_c6 do i=1,nbox_max no(i)=no_c6(i) dist(i)=dist_c6(i) do j=1,nhist sk1(j,i)=sk1_c6(j,i) xls1(j,i)=xls1_c6(j,i) xln1(j,i)=xln1_c6(j,i) xld1(j,i)=xld1_c6(j,i) enddo enddo do j=1,nhist thist(j)=thist_c6(j) enddo else write(*,*)'isot must be 1 for ib=3!!' write(*,*)'stop at mztud/372' stop endif else if (ib.eq.4) then if(isot.eq.1) then !Case 7 mm=mm_c7 nbox=nbox_c7 tmin=tmin_c7 tmax=tmax_c7 do i=1,nbox_max no(i)=no_c7(i) dist(i)=dist_c7(i) do j=1,nhist sk1(j,i)=sk1_c7(j,i) xls1(j,i)=xls1_c7(j,i) xln1(j,i)=xln1_c7(j,i) xld1(j,i)=xld1_c7(j,i) enddo enddo do j=1,nhist thist(j)=thist_c7(j) enddo else write(*,*)'isot must be 1 for ib=4!!' write(*,*)'stop at mztud/396' stop endif else write(*,*)'ib must be 1,2,3 or 4!!' write(*,*)'stop at mztud/401' endif ! write (*,*) 'hisfile: ', hisfile ! the argument to rhist is to make this compatible with mztf_comp.f, ! which is a useful modification of mztf.f (to change strengths of bands ! call rhist (1.0) if (isot.ne.5) deltanux = deltanu(isot,ib) if (isot.eq.5) deltanux = deltanuco c****** c****** calculation of tauinf(nl) c****** call initial ff=1.0e10 do i=nl,1,-1 if(i.eq.nl)then call intz (zl(i),c2,p2,mr2,t2, con) do kr=1,nbox ta(kr)=t2 end do ! write (*,*) ' i, t2 =', i, t2 call interstrength (st2,t2,ka,ta) aa = p2 * coninf * mr2 * (st2 * ff) bb = p2 * coninf * st2 cc = coninf * st2 dd = t2 * coninf * st2 do kr=1,nbox ccbox(kr) = coninf * ka(kr) ddbox(kr) = t2 * ccbox(kr) ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) else call intz (zl(i),c1,p1,mr1,t1, con) do kr=1,nbox ta(kr)=t1 end do ! write (*,*) ' i, t1 =', i, t1 call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + @ ( c1box(kr) + c2box(kr) )/2.d0 ddbox(kr) = ddbox(kr) + @ ( t1*c1box(kr)+t2*c2box(kr) )/2.d0 end do mr2 = mr1 c2=c1 do kr=1,nbox c2box(kr) = c1box(kr) end do t2=t1 p2=p1 end if pt = bb / cc pp = aa / (cc*ff) ! ta=dd/cc ! tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do ! write (*,*) ' i, ts =', i, ts call interstrength(st,ts,ka,ta) ! call intershape(alsa,alna,alda,tdop) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path uapp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mztud/504',ua(kr),ccbox(kr), $ ka(kr),beta,mu(im),kr,im,i,nl call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tauinf(i) = exp( - eqwmu / dble(deltanux) ) end do ! if ( isot.eq.1 .and. ib.eq.2 ) then ! write (*,*) ' tauinf(nl) = ', tauinf(nl) ! write (*,*) ' tauinf(1) = ', tauinf(1) ! endif c****** c****** calculation of tau(in,ir) for n<=r c****** do 1 in=1,nl-1 call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) do 2 ir=in,nl-1 if (ir.eq.in) then tau(in,ir) = 1.d0 goto 2 end if call intz (zl(ir), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5 c2 = c2 * st2 * dble(deltaz) c aa = aa + ( p1*mr1*c1 + p2*mr2*c2 ) / 2.d0 aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0 ddbox(kr) = ddbox(kr) + $ ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0 end do mr1=mr2 t1=t2 c1=c2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) * ta=dd/cc * tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength(st,ts,ka,ta) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tau(in,ir) = exp( - eqwmu / dble(deltanux) ) 2 continue 1 continue ! if ( isot.eq.1 .and. ib.eq.2 ) then ! write (*,*) ' tau(1,*) , *=1,20 ' ! write (*,*) ( sngl(tau(1,k)), k=1,20 ) ! endif c********** c********** calculation of tau(in,ir) for n>r c********** in=nl call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) do 4 ir=in-1,1,-1 call intz (zl(ir), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0 ddbox(kr) = ddbox(kr) + $ ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0 end do mr1=mr2 c1=c2 t1=t2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) ts = dd / cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength (st,ts,ka,ta) call intershape (alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mztud/691',ua(kr),ccbox(kr), $ ka(kr),beta,mu(im),kr,im,i,nl call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tau(in,ir) = exp( - eqwmu / dble(deltanux) ) 4 continue c c due to the simmetry of the transmittances c do in=nl-1,2,-1 do ir=in-1,1,-1 tau(in,ir) = tau(ir,in) end do end do ccc ccc writing out transmittances ccc if (itauout.eq.1) then ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat', ! @ access='sequential', form='unformatted' ) ! else ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat', ! @ access='sequential', form='unformatted' ) ! endif ! write(1) dummy ! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)' ! do in=1,nl ! write (1) tauinf(in), ( tau(in,ir), ir=1,nl ) ! end do ! close(unit=1) elseif (itauout.eq.2) then ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat') ! else ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat') ! endif ! !write(1,*) dummy ! !write(1,*) 'tij for curtis matrix calculations ' ! !write(1,*)' cira mars model atmosphere ' ! !write(1,*)' beta= ',beta,'deltanu= ',deltanux ! write(1,*) nl ! write(1,*) ! @ ' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)' ! do in=1,nl ! write (1,*) tauinf(in) ! write (1,*) (tau(in,ir), ir=1,nl) ! end do ! close(unit=1) ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (*,'(1x, 31htransmitances written out in: ,a22)') ! @ 'taul'//isotcode//dn//ibcode1 ! else ! write (*,'(1x, 31htransmitances written out in: ,a22)') ! @ 'taul'//isotcode//dn//ibcode2 ! endif end if c cleaning of transmittances ! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy, ! @ isotcode,dn,ibcode2) c construction of the curtis matrix call mzcud ( tauinf,tau, cf,cfup,cfdw, vc,taugr, @ ib,isot,icfout,itableout ) c end return end c*********************************************************************** c mzcud.f c*********************************************************************** subroutine mzcud( tauinf,tau, c,cup,cdw,vc,taugr, @ ib,isot,icfout,itableout ) c old times mlp first version of mzcf c a.k.murphy method to avoid extrapolation in the curtis matrix c feb-89 malv AKM method to avoid extrapolation in C.M. c 25-sept-96 cristina dejar las matrices en doble precision c jan 98 malv version para mz1d c oct 01 malv update version for fluxes for hb and fb c jul 2011 malv+fgg Adapted to LMD-MGCM c*********************************************************************** implicit none include 'comcstfi.h' include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o real*8 vc(nl), taugr(nl) ! o real*8 tau(nl,nl) ! i real*8 tauinf(nl) ! i integer ib ! i integer isot ! i integer icfout, itableout ! i c external external bandid character*2 bandid c local variables integer i, in, ir, iw, itblout real*8 cfup(nl,nl), cfdw(nl,nl) real*8 a(nl,nl), cf(nl,nl) character isotcode*2, bcode*2 c formats 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** if (isot.eq.1) isotcode = '26' if (isot.eq.2) isotcode = '28' if (isot.eq.3) isotcode = '36' if (isot.eq.4) isotcode = '27' if (isot.eq.5) isotcode = 'co' bcode = bandid( ib ) ! write (*,*) ' ' do in=1,nl do ir=1,nl cf(in,ir) = 0.0d0 cfup(in,ir) = 0.0d0 cfdw(in,ir) = 0.0d0 c(in,ir) = 0.0d0 cup(in,ir) = 0.0d0 cdw(in,ir) = 0.0d0 a(in,ir) = 0.0d0 end do vc(in) = 0.0d0 taugr(in) = 0.0d0 end do c the next lines are a reduced and equivalent way of calculating c the c(in,ir) elements for n=2,nl1 and r=1,nl c do in=2,nl1 c do ir=1,nl c if(ir.eq.1)then c c(in,ir)=tau(in-1,1)-tau(in+1,1) c elseif(ir.eq.nl)then c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1) c else c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir) c end if c c(in,ir)=c(in,ir)*pi*deltanu(ib)/(2.*deltaz*1.0e5) c end do c end do c go to 1000 c calculation of the matrix cfup(nl,nl) cfup(1,1) = 1.d0 - tau(1,1) do in=2,nl do ir=1,in if (ir.eq.1) then cfup(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.in) then cfup(in,ir) = 1.d0 - tau(in,ir-1) else cfup(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do ! contribution to upwards fluxes from bb at bottom : do in=1,nl taugr(in) = tau(in,1) enddo c calculation of the matrix cfdw(nl,nl) cfdw(nl,nl) = 1.d0 - tauinf(nl) do in=1,nl-1 do ir=in,nl if (ir.eq.in) then cfdw(in,ir) = 1.d0 - tau(in,ir) elseif (ir.eq.nl) then cfdw(in,ir) = tau(in,ir-1) - tauinf(in) else cfdw(in,ir) = tau(in,ir-1) - tau(in,ir) end if end do end do c calculation of the matrix cf(nl,nl) do in=1,nl do ir=1,nl if (ir.eq.1) then ! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below) ! cf(in,ir) = tau(in,ir) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below) cf(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.nl) then cf(in,ir) = tauinf(in) - tau(in,ir-1) else cf(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do c definition of the a(nl,nl) matrix do in=2,nl-1 do ir=1,nl if (ir.eq.in+1) a(in,ir) = -1.d0 if (ir.eq.in-1) a(in,ir) = +1.d0 a(in,ir) = a(in,ir) / ( 2.d0*deltaz*1.d5 ) end do end do ! this is not needed anymore in the akm scheme ! a(1,1) = +3.d0 ! a(1,2) = -4.d0 ! a(1,3) = +1.d0 ! a(nl,nl) = -3.d0 ! a(nl,nl1) = +4.d0 ! a(nl,nl2) = -1.d0 c calculation of the final curtis matrix ("reduced" by murphy's method) if (isot.ne.5) then do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib) cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib) cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib) end do taugr(in) = taugr(in) * pi*deltanu(isot,ib) end do else do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanuco enddo taugr(in) = taugr(in) * pi*deltanuco enddo endif do in=2,nl-1 do ir=1,nl do i=1,nl ! only c contains the matrix a. matrixes cup,cdw dont because ! these two will be used for flux calculations, not ! only for flux divergencies c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir) ! from this matrix we will extract (see below) the ! nl2 x nl2 "core" for the "reduced" final curtis matrix. end do cup(in,ir) = cfup(in,ir) cdw(in,ir) = cfdw(in,ir) end do ! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above) !vc(in) = c(in,1) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above) if (isot.ne.5) then vc(in) = pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) else vc(in) = pi*deltanuco/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) endif end do 5 continue ! write (*,*) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl) ! call elimin_dibuja(c,nl,itableout) c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa): c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw) iw = nan if (isot.eq.4) iw = 5 ! eliminates values < 1.d-19 if (itableout.eq.30) then itblout = 0 else itblout = itableout endif call elimin_mz1d (c,vc,0,iw,itblout,nw) ! upper boundary condition ! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==> do in=2,nl-1 c(in,nl-2) = c(in,nl-2) - c(in,nl) c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl) cup(in,nl-2) = cup(in,nl-2) - cup(in,nl) cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl) cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl) cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl) end do ! j(nl) = j(nl1) ==> ! do in=2,nl1 ! c(in,nl1) = c(in,nl1) + c(in,nl) ! end do ! 1000 continue if (icfout.eq.1) then ! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then ! codmatrx = codmatrx_fb ! else ! codmatrx = codmatrx_hot ! end if ! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! ibcode2 = '0'//ibcode1 ! else ! write ( ibcode2, 202) ib ! endif ! open ( 1, access='sequential', form='unformatted', file= ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 2, access='sequential', form='unformatted', file= ! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 3, access='sequential', form='unformatted', file= ! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 4, access='sequential', form='unformatted', file= ! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat') ! write(1) dummy ! write(1) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' ! do in=2,nl-1 ! write(1) vc(in), (c(in,ir) , ir=2,nl-1 ) !! write (*,*) in, vc(in) ! end do ! write(2) dummy ! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(2) ( cup(in,ir) , ir=1,nl ) ! end do ! write(3) dummy ! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(3) (cdw(in,ir) , ir=1,nl ) ! end do ! write(4) dummy ! write(4)' format: (taugr(n), n=1,nl)' ! do in=1,nl ! write(4) (taugr(in), ir=1,nl ) ! end do ! !write (*,*) ' Last value in file: ', taugr(nl) ! write (*,'(1x,30hcurtis matrix written out in: ,a,a,a,a)' ) ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat', ! @ dircurtis//'cflgr'//isotcode//dn//ibcode2//codmatrx//'.dat' ! close (1) ! close (2) ! close (3) ! close (4) else ! write (*,*) ' no curtis matrix output file ', char(10) end if if (itableout.eq.30) then ! Force output of C.M. in ascii file ! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then ! codmatrx = codmatrx_fb ! else ! codmatrx = codmatrx_hot ! end if ! if (ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! ibcode2 = '0'//ibcode1 ! else ! write ( ibcode2, 202) ib ! endif ! open (10, file= ! & dircurtis//'table'//isotcode//dn//ibcode2//codmatrx//'.dat') ! write(10,*) nl, ' = number of layers ' ! write(10,*) ' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' ! do in=2,nl-1 ! write(10,*) vc(in), (c(in,ir) , ir=2,nl-1 ) ! enddo ! close (10) endif c end return end c*********************************************************************** c mztvc c*********************************************************************** subroutine mztvc ( ig,vc, ib,isot, @ iirw,iimu,itauout,icfout,itableout ) c jul 2011 malv+fgg c*********************************************************************** implicit none include 'comcstfi.h' include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig ! ADDED FOR TRACEBACK real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o real*8 vc(nl), taugr(nl) ! o integer ib ! i integer isot ! i integer iirw ! i integer iimu ! i integer itauout ! i integer icfout ! i integer itableout ! i c local variables and constants integer i, in, ir, im, k ,j integer nmu parameter (nmu = 8) real*8 tau(nl,nl) real*8 tauinf(nl) real*8 con(nzy), coninf real*8 c1, c2 real*8 t1, t2 real*8 p1, p2 real*8 mr1, mr2 real*8 st1, st2 real*8 c1box(70), c2box(70) real*8 ff ! to avoid too small numbers real*8 tvtbs(nzy) real*8 st, beta, ts, eqwmu real*8 mu(nmu), amu(nmu) real*8 zld(nl), zyd(nzy) real*8 correc real deltanux ! width of vib-rot band (cm-1) character isotcode*2 integer idummy real*8 Desp,wsL c formats 111 format(a1) 112 format(a2) 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** c some needed values ! rl=sqrt(log(2.d0)) ! pi2 = 3.14159265358989d0 beta = 1.8d0 ! beta = 1.0d0 idummy = 0 Desp = 0.0d0 wsL = 0.0d0 !write (*,*) ' MZTUD/ iirw = ', iirw c esto es para que las subroutines de mztfsub calculen we c de la forma apropiada para mztf, no para fot icls=icls_mztf c codigos para filenames ! if (isot .eq. 1) isotcode = '26' ! if (isot .eq. 2) isotcode = '28' ! if (isot .eq. 3) isotcode = '36' ! if (isot .eq. 4) isotcode = '27' ! if (isot .eq. 5) isotcode = '62' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (ibcode1,101) ib ! else ! write (ibcode2,202) ib ! endif ! write (*,'( 30h calculating curtis matrix : ,2x, ! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot c integration in angle !!!!!!!!!!!!!!!!!!!! c------- diffusivity approx. if (iimu.eq.1) then ! write (*,*) ' diffusivity approx. beta = ',beta mu(1) = 1.0d0 amu(1)= 1.0d0 c-------data for 8 points integration elseif (iimu.eq.4) then write (*,*)' 4 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.339981043584856)/2.0d0 mu(2)=(1.0d0-0.339981043584856)/2.0d0 mu(3)=(1.0d0+0.861136311594053)/2.0d0 mu(4)=(1.0d0-0.861136311594053)/2.0d0 amu(1)=0.652145154862546 amu(2)=amu(1) amu(3)=0.347854845137454 amu(4)=amu(3) beta=1.0d0 c-------data for 8 points integration elseif(iimu.eq.8) then write (*,*)' 8 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.183434642495650)/2.0d0 mu(2)=(1.0d0-0.183434642495650)/2.0d0 mu(3)=(1.0d0+0.525532409916329)/2.0d0 mu(4)=(1.0d0-0.525532409916329)/2.0d0 mu(5)=(1.0d0+0.796666477413627)/2.0d0 mu(6)=(1.0d0-0.796666477413627)/2.0d0 mu(7)=(1.0d0+0.960289856497536)/2.0d0 mu(8)=(1.0d0-0.960289856497536)/2.0d0 amu(1)=0.362683783378362 amu(2)=amu(1) amu(3)=0.313706645877887 amu(4)=amu(3) amu(5)=0.222381034453374 amu(6)=amu(5) amu(7)=0.101228536290376 amu(8)=amu(7) beta=1.0d0 end if c!!!!!!!!!!!!!!!!!!!!!!! ccc ccc determine abundances included in the absorber amount ccc c first, set up the grid ready for interpolation. do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble(zl(i)) enddo c vibr. temp of the bending mode : if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) !if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 ) c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state) c por similitud a la que se hace en cza.for ; esto solo se hace para CO2 !write (*,*) 'imr(isot) = ', isot, imr(isot) do i=1,nzy if (isot.eq.5) then con(i) = dble( coy(i) * imrco ) else con(i) = dble( co2y(i) * imr(isot) ) correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) ) con(i) = con(i) * ( 1.d0 - correc ) ! write (*,*) ' iz, correc, co2y(i), con(i) =', ! @ i,correc,co2y(i),con(i) endif !----------------------------------------------------------------- ! mlp & cristina. 17 july 1996 change the calculation of mr. ! it is used for calculating partial press ! alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) ! for an isotope, if mr is obtained by ! co2*imr(iso)/nt ! we are considerin collisions with other co2 isotopes ! (including the major one, 626) as if they were with n2. ! assuming mr as co2/nt, we consider collisions ! of type 628-626 as of 626-626 instead of as 626-n2. ! mrx(i)=con(i)/ntx(i) ! old malv ! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs ! jan 98: ! esta modif de mlp implica anular el correc (deberia revisar esto) mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98 !----------------------------------------------------------------- end do ! como beta y 1.d5 son comunes a todas las weighted absorber amounts, ! los simplificamos: ! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) ) !write (*,*) ' con(nz), con(nz-1) =', con(nz), con(nz-1) coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) ) !write (*,*) ' coninf =', coninf ccc ccc temp dependence of the band strength and ccc nlte correction factor for the absorber amount ccc call mztf_correccion ( coninf, con, ib, isot, itableout ) ccc ccc reads histogrammed spectral data (strength for lte and vmr=1) ccc !hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto !! hfile1 = dirspec//'hid' !(see why in his.for) ! hfile1='hid' !! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' ! else ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' ! endif ! write (*,*) 'hisfile: ', hisfile ! the argument to rhist is to make this compatible with mztf_comp.f, ! which is a useful modification of mztf.f (to change strengths of bands ! call rhist (1.0) if(ib.eq.1) then if(isot.eq.1) then !Case 1 mm=mm_c1 nbox=nbox_c1 tmin=tmin_c1 tmax=tmax_c1 do i=1,nbox_max no(i)=no_c1(i) dist(i)=dist_c1(i) do j=1,nhist sk1(j,i)=sk1_c1(j,i) xls1(j,i)=xls1_c1(j,i) xln1(j,i)=xln1_c1(j,i) xld1(j,i)=xld1_c1(j,i) enddo enddo do j=1,nhist thist(j)=thist_c1(j) enddo else if(isot.eq.2) then !Case 2 mm=mm_c2 nbox=nbox_c2 tmin=tmin_c2 tmax=tmax_c2 do i=1,nbox_max no(i)=no_c2(i) dist(i)=dist_c2(i) do j=1,nhist sk1(j,i)=sk1_c2(j,i) xls1(j,i)=xls1_c2(j,i) xln1(j,i)=xln1_c2(j,i) xld1(j,i)=xld1_c2(j,i) enddo enddo do j=1,nhist thist(j)=thist_c2(j) enddo else if(isot.eq.3) then !Case 3 mm=mm_c3 nbox=nbox_c3 tmin=tmin_c3 tmax=tmax_c3 do i=1,nbox_max no(i)=no_c3(i) dist(i)=dist_c3(i) do j=1,nhist sk1(j,i)=sk1_c3(j,i) xls1(j,i)=xls1_c3(j,i) xln1(j,i)=xln1_c3(j,i) xld1(j,i)=xld1_c3(j,i) enddo enddo do j=1,nhist thist(j)=thist_c3(j) enddo else if(isot.eq.4) then !Case 4 mm=mm_c4 nbox=nbox_c4 tmin=tmin_c4 tmax=tmax_c4 do i=1,nbox_max no(i)=no_c4(i) dist(i)=dist_c4(i) do j=1,nhist sk1(j,i)=sk1_c4(j,i) xls1(j,i)=xls1_c4(j,i) xln1(j,i)=xln1_c4(j,i) xld1(j,i)=xld1_c4(j,i) enddo enddo do j=1,nhist thist(j)=thist_c4(j) enddo else write(*,*)'isot must be 2,3 or 4 for ib=1!!' write(*,*)'stop at mztvc/310' stop endif else if (ib.eq.2) then if(isot.eq.1) then !Case 5 mm=mm_c5 nbox=nbox_c5 tmin=tmin_c5 tmax=tmax_c5 do i=1,nbox_max no(i)=no_c5(i) dist(i)=dist_c5(i) do j=1,nhist sk1(j,i)=sk1_c5(j,i) xls1(j,i)=xls1_c5(j,i) xln1(j,i)=xln1_c5(j,i) xld1(j,i)=xld1_c5(j,i) enddo enddo do j=1,nhist thist(j)=thist_c5(j) enddo else write(*,*)'isot must be 1 for ib=2!!' write(*,*)'stop at mztvc/334' stop endif else if (ib.eq.3) then if(isot.eq.1) then !Case 6 mm=mm_c6 nbox=nbox_c6 tmin=tmin_c6 tmax=tmax_c6 do i=1,nbox_max no(i)=no_c6(i) dist(i)=dist_c6(i) do j=1,nhist sk1(j,i)=sk1_c6(j,i) xls1(j,i)=xls1_c6(j,i) xln1(j,i)=xln1_c6(j,i) xld1(j,i)=xld1_c6(j,i) enddo enddo do j=1,nhist thist(j)=thist_c6(j) enddo else write(*,*)'isot must be 1 for ib=3!!' write(*,*)'stop at mztvc/358' stop endif else if (ib.eq.4) then if(isot.eq.1) then !Case 7 mm=mm_c7 nbox=nbox_c7 tmin=tmin_c7 tmax=tmax_c7 do i=1,nbox_max no(i)=no_c7(i) dist(i)=dist_c7(i) do j=1,nhist sk1(j,i)=sk1_c7(j,i) xls1(j,i)=xls1_c7(j,i) xln1(j,i)=xln1_c7(j,i) xld1(j,i)=xld1_c7(j,i) enddo enddo do j=1,nhist thist(j)=thist_c7(j) enddo else write(*,*)'isot must be 1 for ib=4!!' write(*,*)'stop at mztvc/382' stop endif else write(*,*)'ib must be 1,2,3 or 4!!' write(*,*)'stop at mztvc/387' endif c****** c****** calculation of tau(1,ir) for 1<=r c****** call initial ff=1.0e10 in=1 tau(in,1) = 1.d0 call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox c1box(kr) = c1 * ka(kr) * dble(deltaz) end do c1 = c1 * st1 * dble(deltaz) do 2 ir=2,nl call intz (zl(ir), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox c2box(kr) = c2 * ka(kr) * dble(deltaz) end do c2 = c2 * st2 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0 ddbox(kr) = ddbox(kr) + $ ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0 end do mr1=mr2 t1=t2 c1=c2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength(st,ts,ka,ta) call intershape(alsa,alna,alda,ta) eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tau(in,ir) = exp( - eqwmu / dble(deltanu(isot,ib)) ) 2 continue c c due to the simmetry of the transmittances c do in=nl,2,-1 tau(in,1) = tau(1,in) end do vc(1) = 0.0d0 vc(nl) = 0.0d0 do in=2,nl-1 ! poner aqui nl-1 luego vc(in) = pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) end do c end return end c*********************************************************************** c mztvc_626fh.F c*********************************************************************** subroutine mztvc_626fh(ig) c jul 2011 malv+fgg c*********************************************************************** implicit none !!!!!!!!!!!!!!!!!!!!!!! ! common variables & constants include 'nlte_paramdef.h' include 'nlte_commons.h' !!!!!!!!!!!!!!!!!!!!!!! ! arguments integer ig ! ADDED FOR TRACEBACK !!!!!!!!!!!!!!!!!!!!!!! ! local variables real*4 cdummy(nl,nl), csngl(nl,nl) real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl) real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor integer itauout,icfout,itableout, interpol,ismooth, isngldble integer i,j,ik,ist,isot,ib,itt !character bandcode*2 character isotcode*2 !character codmatrx_hot*5 !!!!!!!!!!!!!!!!!!!!!!! ! external functions external bandid character*2 bandid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! subroutines called: ! mz4sub, dmzout, readc_mz4, readcupdw, mztf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! formatos 132 format(i2) ************************************************************************ ************************************************************************ isngldble = 1 ! =1 --> dble precission fileroot = 'cfl' ist = 1 isot = 26 write (isotcode,132) isot call zerov( vc121, nl ) do 11, ik=1,3 ib=ik+1 call mztvc (ig,v1, ib, 1, irw_mztf, imu, 0,0,0 ) do i=1,nl if(ik.eq.1)then vc_factor = dble(667.75/618.03) elseif(ik.eq.2)then vc_factor = 1.d0 elseif(ik.eq.3)then vc_factor = dble(667.75/720.806) end if vc121(i) = vc121(i) + v1(i) * vc_factor end do 11 continue return end c*********************************************************************** c mztf_correccion c*********************************************************************** subroutine mztf_correccion (coninf, con, ib, isot, icurt_pop) c including the dependence of the absort. coeff. on temp., vibr. temp., c function, etc.., when neccessary. imr is already corrected in his.for c we follow pg.39b-43a (l5): c tvt1 is the vibr temp of the upper level c tvt is the vibr temp of the transition itself c tvtbs is the vibr temp of the bending mode (used in qv) c for fundamental bands, they are not used at the moment. c for the 15 fh and sh bands, only tvt0 is used at the moment. c for the laser band, all of them are used following pg. 41a -l5- : c we need s(z) and we can read s(tk) from the histogram (also called c what we have to calculate now is the factor s(z)/s(tk) or following c l5 notebook notation, s_nlte/s_lte. c s_nlte/s_lte = xfactor = xlower * xqv * xes c icurt_pop = 30 -> Output of populations of the 0200,0220,1000 states c = otro -> no output of these populations c oct 92 malv c jan 98 malv version for mz1d c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ib, isot integer icurt_pop ! output of Fermi states population real*8 con(nzy), coninf ! local variables integer i real*8 tvt0(nzy),tvt1(nzy),tvtbs(nzy), zld(nl),zyd(nzy) real xalfa, xbeta, xtv1000, xtv0200, xtv0220, xfactor real xqv, xnu_trans, xtv_trans, xes, xlower c*********************************************************************** xfactor = 1.0 do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble( zl(i) ) end do ! tvtbs is the bending mode of the molecule. used in xqv. if (isot.eq.1) call interdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp (tvtbs,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp (tvtbs,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp (tvtbs,zyd,nzy, v627t1,zld,nl, 1 ) if (isot.eq.5) call interdp (tvtbs,zyd,nzy, vcot1,zld,nl, 1 ) ! tvt0 is the lower level of the transition. used in xlower. if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4 .or. ib.eq.15) then if (isot.eq.1) call interdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 ) if (isot.eq.2) call interdp (tvt0,zyd,nzy, v628t1,zld,nl, 1 ) if (isot.eq.3) call interdp (tvt0,zyd,nzy, v636t1,zld,nl, 1 ) if (isot.eq.4) call interdp (tvt0,zyd,nzy, v627t1,zld,nl, 1 ) elseif (ib.eq.6 .or. ib.eq.8 .or. ib.eq.10 @ .or. ib.eq.13 .or. ib.eq.14 @ .or. ib.eq.17 .or. ib.eq.19 .or. ib.eq.20) then if (isot.eq.1) call interdp ( tvt0,zyd,nzy, v626t2,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvt0,zyd,nzy, v628t2,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvt0,zyd,nzy, v636t2,zld,nl, 1 ) if (isot.eq.4) then call interdp ( tvt0,zyd,nzy, v627t2,zld,nl, 1 ) endif else do i=1,nzy tvt0(i) = dble( ty(i) ) end do end if c tvt is the vt of the transition. used in xes. c since xes=1.0 except for the laser bands, tvt is only needed for them c but it is actually calculated from the tv of the upper and lower level c of the transition. hence, only tvt1 remains to be read for the laser b c tvt1 is the upper level of the transition. if (ib.eq.13 .or. ib.eq.14) then if (isot.eq.1) call interdp ( tvt1,zyd,nzy, v626t4,zld,nl, 1 ) if (isot.eq.2) call interdp ( tvt1,zyd,nzy, v628t4,zld,nl, 1 ) if (isot.eq.3) call interdp ( tvt1,zyd,nzy, v636t4,zld,nl, 1 ) if (isot.eq.4) call interdp ( tvt1,zyd,nzy, v627t4,zld,nl, 1 ) end if c here we weight the absorber amount by a factor which compensate the l c value of the strength read from hitran. we use that factor in order t c correct the product s*m when we later multiply those two variables. ! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then ! open (30, file='020populations.dat') ! write (30,*) ' z tv(020) tv0200 tv0220 tv1000 ' ! endif do i=1,nzy if (isot.eq.1) then !!! vt of the 3 levels in (020) (see pag. 36a-sn1 for this) xalfa = 1.d0/2.d0* exp( dble(-ee*(nu12_1000-nu(1,2))/ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu12_0200-nu(1,2))/ty(i)) ) xtv0200 = dble( - ee * nu12_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) xtv0220 = dble( - ee * nu(1,2) ) / @ ( log( 1.d0/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) xtv1000 = dble( - ee * nu12_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - @ dble(ee*nu(1,2))/tvt0(i) ) !!! correccion 8-Nov-04 (see pag.9b-Marte4-) xtv0200 = dble( - ee * nu12_0200 / @ (log(4.*xbeta/(1.d0+xalfa+xbeta))-ee*nu(1,2)/tvt0(i))) xtv0220 = dble( - ee * nu(1,2) / @ ( log(2./(1.d0+xalfa+xbeta)) - ee*nu(1,2)/tvt0(i) ) ) xtv1000 = dble( - ee * nu12_1000 / @ ( log(4.*xalfa/(1.d0+xalfa+xbeta))-ee*nu(1,2)/tvt0(i))) ! if ( icurt_pop.eq.30 ) then ! write (30,'( 1x,f7.2, 3x,f8.3, 2x,3(1x,f8.3) )') ! @ zx(i), tvt0(i), xtv0200, xtv0220, xtv1000 ! endif !!! xlower and xes for the band if (ib.eq.19) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xes = 1.0d0 elseif (ib.eq.17) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xes = 1.0d0 elseif (ib.eq.20) then xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0220 ) ) xes = 1.0d0 elseif (ib.eq.14) then xlower = exp( dble(ee*nu12_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(1,4)-nu12_1000 ) xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)- @ nu12_1000/xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu12_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble(nu(1,4)-nu12_0200) xtv_trans = xnu_trans / dble(nu(1,4)/tvt1(i)- @ nu12_0200/xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*667.3801/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*667.3801/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.2) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu22_1000-nu(2,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu22_0200-nu(2,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu22_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu22_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(2,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu22_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble(nu(2,4)-nu22_1000) xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu22_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(2,4)-nu22_0200 ) xtv_trans = xnu_trans / dble(nu(2,4)/tvt1(i)-nu22_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*662.3734/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*662.3734/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.3) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu32_1000-nu(3,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu32_0200-nu(3,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu32_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu32_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(3,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu32_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(3,4)-nu32_1000 ) xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu32_0200) * @ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(3,4)-nu32_0200 ) xtv_trans = xnu_trans / dble(nu(3,4)/tvt1(i)-nu32_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*648.4784/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*648.4784/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.4) then xalfa = 1.d0/2.d0* exp( dble(-ee*(nu42_1000-nu(4,2))/ @ ty(i)) ) xbeta = 1.d0/2.d0* exp( dble(-ee*(nu42_0200-nu(4,2))/ @ ty(i)) ) xtv0200 = dble( - ee * nu42_0200 ) / @ ( log( xbeta/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/ @ tvt0(i) ) xtv1000 = dble( - ee * nu42_1000 ) / @ ( log( xalfa/(1.d0+xalfa+xbeta) ) - dble(ee*nu(4,2))/ @ tvt0(i) ) if (ib.eq.14) then xlower = exp( dble(ee*nu42_1000) * @ ( 1.d0/dble(ty(i))-1.d0/xtv1000 ) ) xnu_trans = dble( nu(4,4)-nu42_1000 ) xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_1000/ @ xtv1000) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) elseif (ib.eq.13) then xlower = exp( dble(ee*nu42_0200) * $ ( 1.d0/dble(ty(i))-1.d0/xtv0200 ) ) xnu_trans = dble( nu(4,4)-nu42_0200 ) xtv_trans = xnu_trans / dble(nu(4,4)/tvt1(i)-nu42_0200/ @ xtv0200) xes = (1.d0-exp( dble(-ee*xnu_trans/xtv_trans) )) / @ (1.d0-exp( dble(-ee*xnu_trans/ty(i)) )) else xlower = exp( dble(ee*elow(isot,ib)) * @ ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) ) xes = 1.0d0 end if xqv = (1.d0-exp( dble(-ee*664.7289/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*664.7289/ty(i)) )) xfactor = xlower * xqv**2.d0 * xes elseif (isot.eq.5 .and. ib.eq.1) then xlower = 1.d0 xes = 1.0d0 xqv = (1.d0-exp( dble(-ee*nuco_10/tvtbs(i)) )) / @ (1.d0-exp( dble(-ee*nuco_10/ty(i)) )) xfactor = xlower * xqv * xes end if con(i) = con(i) * xfactor if (i.eq.nzy) coninf = coninf * xfactor end do ! if ( isot.eq.1 .and. icurt_pop.eq.30 ) then ! close (30) ! endif return end c*********************************************************************** c mztf.f c*********************************************************************** c c program for calculating atmospheric transmittances c to be used in the calculation of curtis matrix coefficients subroutine mztf ( ig,cf,cfup,cfdw,vc,taugr, ib,isot, @ iirw,iimu,itauout,icfout,itableout ) c i*out = 1 output of data c i*out = 0 no output c c jul 2011 malv+fgg adapted to LMD-MGCM c nov 98 mavl allow for overlaping in the lorentz line c jan 98 malv version for mz1d. based on curtis/mztf.for c 17-jul-96 mlp&crs change the calculation of mr. c evitar: divide por cero. anhadiendo: ff c oct-92 malv correct s(t) dependence for all histogr bands c june-92 malv proper lower levels for laser bands c may-92 malv new temperature dependence for laser bands c @ 991 malv boxing for the averaged absorber amount and t c ? malv extension up to 200 km altitude in mars c 13-nov-86 mlp include the temperature weighted to match c the eqw in the strong doppler limit. c*********************************************************************** implicit none include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments integer ig !ADDED FOR TRACEBACK real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o. real*8 vc(nl), taugr(nl) ! o integer ib ! i integer isot ! i integer iirw ! i integer iimu ! i integer itauout ! i integer icfout ! i integer itableout ! i c local variables and constants integer i, in, ir, im, k ,j integer nmu parameter (nmu = 8) real*8 tau(nl,nl) real*8 tauinf(nl) real*8 con(nzy), coninf real*8 c1, c2 real*8 t1, t2 real*8 p1, p2 real*8 mr1, mr2 real*8 st1, st2 real*8 c1box(70), c2box(70) real*8 ff ! to avoid too small numbers real*8 tvtbs(nzy) real*8 st, beta, ts, eqwmu real*8 mu(nmu), amu(nmu) real*8 zld(nl), zyd(nzy) real*8 correc real deltanux ! width of vib-rot band (cm-1) ! character isotcode*2 integer idummy real*8 Desp,wsL c formats ! 111 format(a1) ! 112 format(a2) 101 format(i1) 202 format(i2) ! 180 format(a80) ! 181 format(a80) c*********************************************************************** c some needed values ! rl=sqrt(log(2.d0)) ! pi2 = 3.14159265358989d0 beta = 1.8d0 idummy = 0 Desp = 0.d0 wsL = 0.d0 c esto es para que las subroutines de mztfsub calculen we c de la forma apropiada para mztf, no para fot icls=icls_mztf c codigos para filenames ! if (isot .eq. 1) isotcode = '26' ! if (isot .eq. 2) isotcode = '28' ! if (isot .eq. 3) isotcode = '36' ! if (isot .eq. 4) isotcode = '27' ! if (isot .eq. 5) isotcode = '62' ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (ibcode1,101) ib ! else ! write (ibcode2,202) ib ! endif ! write (*,'( 30h calculating curtis matrix : ,2x, ! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot c integration in angle !!!!!!!!!!!!!!!!!!!! c------- diffusivity approx. if (iimu.eq.1) then ! write (*,*) ' diffusivity approx. beta = ',beta mu(1) = 1.0d0 amu(1)= 1.0d0 c-------data for 8 points integration elseif (iimu.eq.4) then write (*,*)' 4 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.339981043584856)/2.0d0 mu(2)=(1.0d0-0.339981043584856)/2.0d0 mu(3)=(1.0d0+0.861136311594053)/2.0d0 mu(4)=(1.0d0-0.861136311594053)/2.0d0 amu(1)=0.652145154862546 amu(2)=amu(1) amu(3)=0.347854845137454 amu(4)=amu(3) beta=1.0d0 c-------data for 8 points integration elseif(iimu.eq.8) then write (*,*)' 8 points for the gauss-legendre angle quadrature.' mu(1)=(1.0d0+0.183434642495650)/2.0d0 mu(2)=(1.0d0-0.183434642495650)/2.0d0 mu(3)=(1.0d0+0.525532409916329)/2.0d0 mu(4)=(1.0d0-0.525532409916329)/2.0d0 mu(5)=(1.0d0+0.796666477413627)/2.0d0 mu(6)=(1.0d0-0.796666477413627)/2.0d0 mu(7)=(1.0d0+0.960289856497536)/2.0d0 mu(8)=(1.0d0-0.960289856497536)/2.0d0 amu(1)=0.362683783378362 amu(2)=amu(1) amu(3)=0.313706645877887 amu(4)=amu(3) amu(5)=0.222381034453374 amu(6)=amu(5) amu(7)=0.101228536290376 amu(8)=amu(7) beta=1.0d0 end if c!!!!!!!!!!!!!!!!!!!!!!! ccc ccc determine abundances included in the absorber amount ccc c first, set up the grid ready for interpolation. do i=1,nzy zyd(i) = dble(zy(i)) enddo do i=1,nl zld(i) = dble(zl(i)) enddo c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state) c por similitud a la que se hace en cza.for do i=1,nzy if (isot.eq.5) then con(i) = dble( coy(i) * imrco ) else con(i) = dble( co2y(i) * imr(isot) ) c vibr. temp of the bending mode : if(isot.eq.1) call interdp(tvtbs,zyd,nzy,v626t1,zld,nl,1) if(isot.eq.2) call interdp(tvtbs,zyd,nzy,v628t1,zld,nl,1) if(isot.eq.3) call interdp(tvtbs,zyd,nzy,v636t1,zld,nl,1) if(isot.eq.4) call interdp(tvtbs,zyd,nzy,v627t1,zld,nl,1) correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) ) con(i) = con(i) * ( 1.d0 - correc ) endif c----------------------------------------------------------------------- c mlp & cristina. 17 july 1996 c change the calculation of mr. it is used for calculating partial press c alpha = alpha(self,co2)*pp +alpha(n2)*(pt-pp) c for an isotope, if mr is obtained by co2*imr(iso)/nt we are considerin c collisions with other co2 isotopes (including the major one, 626) c as if they were with n2. assuming mr as co2/nt, we consider collisions c of type 628-626 as of 626-626 instead of as 626-n2. c mrx(i)=con(i)/ntx(i) ! old malv ! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs c jan 98: c esta modif de mlp implica anular el correc (deberia revisar esto) mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98 c----------------------------------------------------------------------- end do ! como beta y 1.d5 son comunes a todas las weighted absorber amounts, ! los simplificamos: ! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) ) coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) ) ! write (*,*) ' coninf =', coninf ccc ccc temp dependence of the band strength and ccc nlte correction factor for the absorber amount ccc call mztf_correccion ( coninf, con, ib, isot, itableout ) ccc ccc reads histogrammed spectral data (strength for lte and vmr=1) ccc !hfile1 = dirspec//'hi'//dn ! ya no distinguimos entre d/n !! hfile1 = dirspec//'hid' ! (see why in his.for) ! hfile='hid' !! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his' ! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his' ! ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat' ! else ! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat' ! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat' ! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat' ! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat' ! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat' ! endif ! write (*,*) 'hisfile: ', hisfile ! the argument to rhist is to make this compatible with mztf_comp.f, ! which is a useful modification of mztf.f (to change strengths of bands ! call rhist (1.0) if(ib.eq.1) then if(isot.eq.1) then !Case 1 mm=mm_c1 nbox=nbox_c1 tmin=tmin_c1 tmax=tmax_c1 do i=1,nbox_max no(i)=no_c1(i) dist(i)=dist_c1(i) do j=1,nhist sk1(j,i)=sk1_c1(j,i) xls1(j,i)=xls1_c1(j,i) xln1(j,i)=xln1_c1(j,i) xld1(j,i)=xld1_c1(j,i) enddo enddo do j=1,nhist thist(j)=thist_c1(j) enddo else if(isot.eq.2) then !Case 2 mm=mm_c2 nbox=nbox_c2 tmin=tmin_c2 tmax=tmax_c2 do i=1,nbox_max no(i)=no_c2(i) dist(i)=dist_c2(i) do j=1,nhist sk1(j,i)=sk1_c2(j,i) xls1(j,i)=xls1_c2(j,i) xln1(j,i)=xln1_c2(j,i) xld1(j,i)=xld1_c2(j,i) enddo enddo do j=1,nhist thist(j)=thist_c2(j) enddo else if(isot.eq.3) then !Case 3 mm=mm_c3 nbox=nbox_c3 tmin=tmin_c3 tmax=tmax_c3 do i=1,nbox_max no(i)=no_c3(i) dist(i)=dist_c3(i) do j=1,nhist sk1(j,i)=sk1_c3(j,i) xls1(j,i)=xls1_c3(j,i) xln1(j,i)=xln1_c3(j,i) xld1(j,i)=xld1_c3(j,i) enddo enddo do j=1,nhist thist(j)=thist_c3(j) enddo else if(isot.eq.4) then !Case 4 mm=mm_c4 nbox=nbox_c4 tmin=tmin_c4 tmax=tmax_c4 do i=1,nbox_max no(i)=no_c4(i) dist(i)=dist_c4(i) do j=1,nhist sk1(j,i)=sk1_c4(j,i) xls1(j,i)=xls1_c4(j,i) xln1(j,i)=xln1_c4(j,i) xld1(j,i)=xld1_c4(j,i) enddo enddo do j=1,nhist thist(j)=thist_c4(j) enddo else write(*,*)'isot must be 2,3 or 4 for ib=1!!' write(*,*)'stop at mztf_overlap/317' stop endif else if (ib.eq.2) then if(isot.eq.1) then !Case 5 mm=mm_c5 nbox=nbox_c5 tmin=tmin_c5 tmax=tmax_c5 do i=1,nbox_max no(i)=no_c5(i) dist(i)=dist_c5(i) do j=1,nhist sk1(j,i)=sk1_c5(j,i) xls1(j,i)=xls1_c5(j,i) xln1(j,i)=xln1_c5(j,i) xld1(j,i)=xld1_c5(j,i) enddo enddo do j=1,nhist thist(j)=thist_c5(j) enddo else write(*,*)'isot must be 1 for ib=2!!' write(*,*)'stop at mztf_overlap/341' stop endif else if (ib.eq.3) then if(isot.eq.1) then !Case 6 mm=mm_c6 nbox=nbox_c6 tmin=tmin_c6 tmax=tmax_c6 do i=1,nbox_max no(i)=no_c6(i) dist(i)=dist_c6(i) do j=1,nhist sk1(j,i)=sk1_c6(j,i) xls1(j,i)=xls1_c6(j,i) xln1(j,i)=xln1_c6(j,i) xld1(j,i)=xld1_c6(j,i) enddo enddo do j=1,nhist thist(j)=thist_c6(j) enddo else write(*,*)'isot must be 1 for ib=3!!' write(*,*)'stop at mztf_overlap/365' stop endif else if (ib.eq.4) then if(isot.eq.1) then !Case 7 mm=mm_c7 nbox=nbox_c7 tmin=tmin_c7 tmax=tmax_c7 do i=1,nbox_max no(i)=no_c7(i) dist(i)=dist_c7(i) do j=1,nhist sk1(j,i)=sk1_c7(j,i) xls1(j,i)=xls1_c7(j,i) xln1(j,i)=xln1_c7(j,i) xld1(j,i)=xld1_c7(j,i) enddo enddo do j=1,nhist thist(j)=thist_c7(j) enddo else write(*,*)'isot must be 1 for ib=4!!' write(*,*)'stop at mztf_overlap/389' stop endif else write(*,*)'ib must be 1,2,3 or 4!!' write(*,*)'stop at mztf_overlap/394' endif if (isot.ne.5) deltanux = deltanu(isot,ib) if (isot.eq.5) deltanux = deltanuco c****** c****** calculation of tauinf(nl) c****** call initial ff=1.0e10 do i=nl,1,-1 if(i.eq.nl)then call intz (zl(i),c2,p2,mr2,t2, con) do kr=1,nbox ta(kr)=t2 end do ! write (*,*) ' i, t2 =', i, t2 call interstrength (st2,t2,ka,ta) aa = p2 * coninf * mr2 * (st2 * ff) bb = p2 * coninf * st2 cc = coninf * st2 dd = t2 * coninf * st2 do kr=1,nbox ccbox(kr) = coninf * ka(kr) ddbox(kr) = t2 * ccbox(kr) ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) else call intz (zl(i),c1,p1,mr1,t1, con) do kr=1,nbox ta(kr)=t1 end do ! write (*,*) ' i, t1 =', i, t1 call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) )/2.d0 ddbox(kr) = ddbox(kr) + (t1*c1box(kr)+t2*c2box(kr))/2.d0 end do mr2 = mr1 c2=c1 do kr=1,nbox c2box(kr) = c1box(kr) end do t2=t1 p2=p1 end if pt = bb / cc pp = aa / (cc*ff) ! ta=dd/cc ! tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do ! write (*,*) ' i, ts =', i, ts call interstrength(st,ts,ka,ta) ! call intershape(alsa,alna,alda,tdop) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mztf_overlap/483',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tauinf(i) = exp( - eqwmu / dble(deltanux) ) end do ! i continue ! if ( isot.eq.1 .and. ib.eq.2 ) then ! write (*,*) ' tauinf(nl) = ', tauinf(nl) ! write (*,*) ' tauinf(1) = ', tauinf(1) ! endif c****** c****** calculation of tau(in,ir) for n<=r c****** do 1 in=1,nl-1 call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) do 2 ir=in,nl-1 if (ir.eq.in) then tau(in,ir) = 1.d0 goto 2 end if call intz (zl(ir), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.e5 c2 = c2 * st2 * dble(deltaz) c aa = aa + ( p1*mr1*c1 + p2*mr2*c2 ) / 2.d0 aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0 ddbox(kr) = ddbox(kr) + $ ( t1*c1box(kr) + t2*c2box(kr) ) /2.d0 end do mr1=mr2 t1=t2 c1=c2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) * ta=dd/cc * tdop = ta ts = dd/cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength(st,ts,ka,ta) call intershape(alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mztf_overlap/581',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tau(in,ir) = exp( - eqwmu / dble(deltanux) ) 2 continue 1 continue ! if ( isot.eq.1 .and. ib.eq.2 ) then ! write (*,*) ' tau(1,*) , *=1,20 ' ! write (*,*) ( sngl(tau(1,k)), k=1,20 ) ! endif c********** c********** calculation of tau(in,ir) for n>r c********** in=nl call initial call intz (zl(in), c1,p1,mr1,t1, con) do kr=1,nbox ta(kr) = t1 end do call interstrength (st1,t1,ka,ta) do kr=1,nbox ! c1box(kr) = c1 * ka(kr) * beta * dble(deltaz) * 1.d5 c1box(kr) = c1 * ka(kr) * dble(deltaz) end do ! c1 = c1 * st1 * beta * dble(deltaz) * 1.d5 c1 = c1 * st1 * dble(deltaz) do 4 ir=in-1,1,-1 call intz (zl(ir), c2,p2,mr2,t2, con) do kr=1,nbox ta(kr) = t2 end do call interstrength (st2,t2,ka,ta) do kr=1,nbox ! c2box(kr) = c2 * ka(kr) * beta * dble(deltaz) * 1.d5 c2box(kr) = c2 * ka(kr) * dble(deltaz) end do ! c2 = c2 * st2 * beta * dble(deltaz) * 1.d5 c2 = c2 * st2 * dble(deltaz) aa = aa + ( p1*mr1*(c1*ff) + p2*mr2*(c2*ff)) / 2.d0 bb = bb + ( p1*c1 + p2*c2 ) / 2.d0 cc = cc + ( c1 + c2 ) / 2.d0 dd = dd + ( t1*c1 + t2*c2 ) / 2.d0 do kr=1,nbox ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0 ddbox(kr) = ddbox(kr) + ( t1*c1box(kr) + t2*c2box(kr) )/2.d0 end do mr1=mr2 c1=c2 t1=t2 p1=p2 do kr=1,nbox c1box(kr) = c2box(kr) end do pt = bb / cc pp = aa / (cc * ff) ts = dd / cc do kr=1,nbox ta(kr) = ddbox(kr) / ccbox(kr) end do call interstrength (st,ts,ka,ta) call intershape (alsa,alna,alda,ta) * ua = cc/st c next loop calculates the eqw for an especified path ua,pp,pt,ta eqwmu = 0.0d0 do im = 1,iimu eqw=0.0d0 do kr=1,nbox ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im) if(ua(kr).lt.0.)write(*,*)'mztf_overlap/674',ua(kr), $ ccbox(kr),ka(kr),beta,mu(im),kr,im,i,nl call findw(ig,iirw, idummy,c1,p1,Desp,wsL) if ( i_supersat .eq. 0 ) then eqw=eqw+no(kr)*w else eqw = w + (no(kr)-1) * ( asat_box*dist(kr) ) endif end do eqwmu = eqwmu + eqw * mu(im)*amu(im) end do tau(in,ir) = exp( - eqwmu / dble(deltanux) ) 4 continue c c due to the simmetry of the transmittances c do in=nl-1,2,-1 do ir=in-1,1,-1 tau(in,ir) = tau(ir,in) end do end do ccc ccc writing out transmittances ccc if (itauout.eq.1) then ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat', ! @ access='sequential', form='unformatted' ) ! else ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat', ! @ access='sequential', form='unformatted' ) ! endif ! write(1) dummy ! write(1)' format: (tauinf(n),(tau(n,r),r=1,nl),n=1,nl)' ! do in=1,nl ! write (1) tauinf(in), ( tau(in,ir), ir=1,nl ) ! end do ! close(unit=1) elseif (itauout.eq.2) then ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode1//'.dat') ! else ! open( 1, file= ! @ dircurtis//'taul'//isotcode//dn//ibcode2//'.dat') ! endif ! !write(1,*) dummy ! !write(1,*) 'tij for curtis matrix calculations ' ! !write(1,*)' cira mars model atmosphere ' ! write(1,*)' beta= ',beta,'deltanu= ',deltanux ! write(1,*)' number of elements (in,ir)= ',nl,nl ! write(1,*)' format: (tauinf(in),(tau(in,ir),ir=1,nl),in=1,nl)' ! do in=1,nl ! write (1,*) tauinf(in) ! do ir=1,nl ! write(1,*) tau(in,ir) ! end do ! end do ! close(unit=1) ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (*,'(1x, 31htransmitances written out in: ,a22)') ! @ 'taul'//isotcode//dn//ibcode1 ! else ! write (*,'(1x, 31htransmitances written out in: ,a22)') ! @ 'taul'//isotcode//dn//ibcode2 ! endif end if c cleaning of transmittances ! call elimin_tau(tau,tauinf,nl,nan,itableout,nw,dummy, ! @ isotcode,dn,ibcode2) c construction of the curtis matrix call mzcf ( tauinf,tau, cf,cfup,cfdw, vc,taugr, @ ib,isot,icfout,itableout ) c end return end c*********************************************************************** c mzcf c*********************************************************************** subroutine mzcf( tauinf,tau, c,cup,cdw,vc,taugr, @ ib,isot,icfout,itableout ) c a.k.murphy method to avoid extrapolation in the curtis matrix c feb-89 m. angel granada c 25-sept-96 cristina dejar las matrices en doble precision c jan 98 malv version para mz1d c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none include 'comcstfi.h' include 'nlte_paramdef.h' include 'nlte_commons.h' c arguments real*8 c(nl,nl), cup(nl,nl), cdw(nl,nl) ! o real*8 vc(nl), taugr(nl) ! o real*8 tau(nl,nl) ! i real*8 tauinf(nl) ! i integer ib ! i integer isot ! i integer icfout, itableout ! i c external external bandid character*2 bandid c local variables integer i, in, ir, iw real*8 cfup(nl,nl), cfdw(nl,nl) real*8 a(nl,nl), cf(nl,nl) character isotcode*2, bcode*2 c formats 101 format(i1) 202 format(i2) 180 format(a80) 181 format(a80) c*********************************************************************** if (isot.eq.1) isotcode = '26' if (isot.eq.2) isotcode = '28' if (isot.eq.3) isotcode = '36' if (isot.eq.4) isotcode = '27' if (isot.eq.5) isotcode = 'co' bcode = bandid( ib ) ! write (*,*) ' ' do in=1,nl do ir=1,nl cf(in,ir) = 0.0d0 cfup(in,ir) = 0.0d0 cfdw(in,ir) = 0.0d0 c(in,ir) = 0.0d0 cup(in,ir) = 0.0d0 cdw(in,ir) = 0.0d0 a(in,ir) = 0.0d0 end do vc(in) = 0.0d0 taugr(in) = 0.0d0 end do c the next lines are a reduced and equivalent way of calculating c the c(in,ir) elements for n=2,nl1 and r=1,nl c do in=2,nl1 c do ir=1,nl c if(ir.eq.1)then c c(in,ir)=tau(in-1,1)-tau(in+1,1) c elseif(ir.eq.nl)then c c(in,ir)=tau(in+1,nl1)-tauinf(in+1)-tau(in-1,nl1)+tauinf(in-1) c else c c(in,ir)=tau(in+1,ir-1)-tau(in+1,ir)-tau(in-1,ir-1)+tau(in-1,ir) c end if c c(in,ir)=c(in,ir)*pi*deltanu(ib)/(2.*deltaz*1.0e5) c end do c end do c go to 1000 c calculation of the matrix cfup(nl,nl) cfup(1,1) = 1.d0 - tau(1,1) do in=2,nl do ir=1,in if (ir.eq.1) then cfup(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.in) then cfup(in,ir) = 1.d0 - tau(in,ir-1) else cfup(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do ! contribution to upwards fluxes from bb at bottom : do in=1,nl taugr(in) = tau(in,1) enddo c calculation of the matrix cfdw(nl,nl) cfdw(nl,nl) = 1.d0 - tauinf(nl) do in=1,nl-1 do ir=in,nl if (ir.eq.in) then cfdw(in,ir) = 1.d0 - tau(in,ir) elseif (ir.eq.nl) then cfdw(in,ir) = tau(in,ir-1) - tauinf(in) else cfdw(in,ir) = tau(in,ir-1) - tau(in,ir) end if end do end do c calculation of the matrix cf(nl,nl) do in=1,nl do ir=1,nl if (ir.eq.1) then ! version con l_bb(tg) = l_bb(t(1))=j(1) (see also vc below) ! cf(in,ir) = tau(in,ir) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see also vc below) cf(in,ir) = tau(in,ir) - tau(in,1) elseif (ir.eq.nl) then cf(in,ir) = tauinf(in) - tau(in,ir-1) else cf(in,ir) = tau(in,ir) - tau(in,ir-1) end if end do end do c definition of the a(nl,nl) matrix do in=2,nl-1 do ir=1,nl if (ir.eq.in+1) a(in,ir) = -1.d0 if (ir.eq.in-1) a(in,ir) = +1.d0 a(in,ir) = a(in,ir) / ( 2.d0*deltaz*1.d5 ) end do end do ! this is not needed anymore in the akm scheme ! a(1,1) = +3.d0 ! a(1,2) = -4.d0 ! a(1,3) = +1.d0 ! a(nl,nl) = -3.d0 ! a(nl,nl1) = +4.d0 ! a(nl,nl2) = -1.d0 c calculation of the final curtis matrix ("reduced" by murphy's method) if (isot.ne.5) then do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanu(isot,ib) cfup(in,ir) = cfup(in,ir) * pi*deltanu(isot,ib) cfdw(in,ir) = cfdw(in,ir) * pi*deltanu(isot,ib) end do taugr(in) = taugr(in) * pi*deltanu(isot,ib) end do else do in=1,nl do ir=1,nl cf(in,ir) = cf(in,ir) * pi*deltanuco enddo taugr(in) = taugr(in) * pi*deltanuco enddo endif do in=2,nl-1 do ir=1,nl do i=1,nl ! only c contains the matrix a. matrixes cup,cdw dont because ! these two will be used for flux calculations, not ! only for flux divergencies c(in,ir) = c(in,ir) + a(in,i) * cf(i,ir) ! from this matrix we will extract (see below) the ! nl2 x nl2 "core" for the "reduced" final curtis matrix. end do cup(in,ir) = cfup(in,ir) cdw(in,ir) = cfdw(in,ir) end do ! version con l_bb(tg) = l_bb(t(1))=j(1) (see cf above) !vc(in) = c(in,1) ! version con l_bb(tg) =/= l_bb(t(1))=j(1) (see cf above) vc(in) = pi*deltanu(isot,ib)/( 2.d0*deltaz*1.d5 ) * @ ( tau(in-1,1) - tau(in+1,1) ) end do 5 continue ! write (*,*) 'mztf/1/ c(2,*) =', (c(2,i), i=1,nl) ! call elimin_dibuja(c,nl,itableout) c ventana del smoothing de c es nw=3 y de vc es 5 (puesto en lisa): c subroutine elimin_mz4(c,vc,ilayer,nl,nan,iw, nw) iw = nan if (isot.eq.4) iw = 5 call elimin_mz1d (c,vc,0,iw,itableout,nw) ! upper boundary condition ! j'(nl) = j'(nl1) ==> j(nl) = 2j(nl1) - j(nl2) ==> do in=2,nl-1 c(in,nl-2) = c(in,nl-2) - c(in,nl) c(in,nl-1) = c(in,nl-1) + 2.d0*c(in,nl) cup(in,nl-2) = cup(in,nl-2) - cup(in,nl) cup(in,nl-1) = cup(in,nl-1) + 2.d0*cup(in,nl) cdw(in,nl-2) = cdw(in,nl-2) - cdw(in,nl) cdw(in,nl-1) = cdw(in,nl-1) + 2.d0*cdw(in,nl) end do ! j(nl) = j(nl1) ==> ! do in=2,nl1 ! c(in,nl1) = c(in,nl1) + c(in,nl) ! end do ! 1000 continue if (icfout.eq.1) then ! if (ib.eq.1 .or. ib.eq.12 .or. ib.eq.16 .or. ib.eq.18) then ! codmatrx = codmatrx_fb ! else ! codmatrx = codmatrx_hot ! end if ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! open ( 1, access='sequential', form='unformatted', file= ! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat') ! open ( 2, access='sequential', form='unformatted', file= ! @ dircurtis//'cflup'//isotcode//dn//ibcode1//codmatrx//'.dat') ! open ( 3, access='sequential', form='unformatted', file= ! @ dircurtis//'cfldw'//isotcode//dn//ibcode1//codmatrx//'.dat') ! else ! open ( 1, access='sequential', form='unformatted', file= ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 2, access='sequential', form='unformatted', file= ! @ dircurtis//'cflup'//isotcode//dn//ibcode2//codmatrx//'.dat') ! open ( 3, access='sequential', form='unformatted', file= ! @ dircurtis//'cfldw'//isotcode//dn//ibcode2//codmatrx//'.dat') ! endif ! write(1) dummy ! write(1)' format: (vc(n),(ch(n,r),r=2,nl-1),n=2,nl-1)' ! do in=2,nl-1 ! write(1) vc(in), (c(in,ir) , ir=2,nl-1 ) ! es mas importante la precision que ocupar mucho espacio asi que ! escribiremos las matrices en doble precision y por tanto en ! [lib]readc_mz4.for no hay que reconvertirlas a doble precision. ! ch is stored in single prec. to save storage space. ! end do ! write(2) dummy ! write(2)' format: (cfup(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(2) ( cup(in,ir) , ir=1,nl ) ! end do ! write(3) dummy ! write(3)' format: (cfdw(n,r),r=1,nl), n=1,nl)' ! do in=1,nl ! write(3) (cdw(in,ir) , ir=1,nl ) ! end do ! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5 ! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then ! write (*,'(1x,30hcurtis matrix written out in: ,a50)' ) ! @ dircurtis//'cfl'//isotcode//dn//ibcode1//codmatrx//'.dat' ! else ! write (*,'(1x,30hcurtis matrix written out in: ,a50)' ) ! @ dircurtis//'cfl'//isotcode//dn//ibcode2//codmatrx//'.dat' ! endif else ! write (*,*) ' no curtis matrix output file ', char(10) end if c end return end c*********************************************************************** c cm15um_hb_simple c*********************************************************************** subroutine cm15um_hb_simple (ig,icurt) c computing the curtix matrixes for the 15 um hot bands c (las de las bandas fudnamentales las calcula cm15um_fb) c jan 98 malv version de mod3/cm_15um.f para mz1d c jul 2011 malv+fgg adapted to LMD-MGCM c*********************************************************************** implicit none !!!!!!!!!!!!!!!!!!!!!!! ! common variables & constants include 'nlte_paramdef.h' include 'nlte_commons.h' !!!!!!!!!!!!!!!!!!!!!!! ! arguments integer ig ! ADDED FOR TRACEBACK integer icurt ! icurt=0,1,2 ! new calculations? (see caa.f heads) !!!!!!!!!!!!!!!!!!!!!!! ! local variables real*4 cdummy(nl,nl), csngl(nl,nl) real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl) real*8 v1(nl), v2(nl), v3(nl), cm_factor, vc_factor integer itauout,icfout,itableout, interpol,ismooth, isngldble integer i,j,ik,ist,isot,ib,itt !character bandcode*2 character isotcode*2 character codmatrx_hot*5 !!!!!!!!!!!!!!!!!!!!!!! ! external functions external bandid character*2 bandid !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! subroutines called: ! mz4sub, dmzout, readc_mz4, readcupdw, mztf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! formatos 132 format(i2) ************************************************************************ ************************************************************************ call zerom (c121,nl) call zerov (vc121,nl) call zerom (cup121,nl) call zerom (cdw121,nl) call zerov (taugr121,nl) itauout = 0 ! =1 --> with output of tau icfout = 0 ! =1 --> with output of cf itableout = 0 ! =1 --> with output of table of taus isngldble = 1 ! =1 --> dble precission codmatrx_hot=' ' if (icurt.eq.2) then icfout=1 elseif (icurt.eq.0) then write (*,'(a,a$)') @ ' hot bands. code for old matrixes (5 chars): ' read (*,'(a)') codmatrx_hot endif fileroot = 'cfl' ! ====================== curtis matrixes for fh bands ================== ! una piedra en el camino ... ! write (*,*) ' cm15um_hb/1 ' ccc if ( input_cza.ge.1 ) then ccc if (icurt.eq.2) then write (*,'(a,a$)') @ ' new calculation of curt. mat. for fh bands.', @ ' code for new matrixes : ' read (*,'(a)') codmatrx_hot elseif (icurt.eq.0) then write (*,'(a,a$)') @ ' reading in curt. mat. for fh bands.', @ ' code for old matrixes : ' read (*,'(a)') codmatrx_hot else ! write (*,'(a)') ! @ ' new calculation of curt. mat. for fh bands.' end if ! fh bands for the 626 isotope ================================- ist = 1 isot = 26 ! encode (2,132,isotcode) isot write (isotcode,132) isot do 11, ik=1,3 ib=ik+1 if (icurt.gt.0) then call zero3m (cax1,cax2,cax3,nl) ! una piedra en el camino ... !write (*,*) ' cm15um_hb/11 ' !write (*,*) ' ib, ist, irw, imu =', ib, ist, irw_mztf, imu call mztf ( ig,cax1,cax2,cax3,v1,v2, ib,ist,irw_mztf,imu, @ itauout,icfout,itableout) ! else ! bandcode = bandid(ib) ! filend=isotcode//dn//bandcode//codmatrx_hot !! write (*,*) char(9), fileroot//filend ! call zero3m (cax1,cax2,cax3,nl) ! call readcud_mz1d ( cax1,cax2,cax3, v1, v2, ! @ fileroot,filend, csngl, nl,nan,0,isngldble) end if c calculating the total c121(n,r) matrix for the first hot band do i=1,nl if ( ib .eq. 4 ) then ! write (*,*) ' ' ! write (*,*) i, ' ib,ist, altura :', ib,ist, zl(i) endif ! if ( v1(i) .le. 1.d-99 ) v1(i) = 0.0d0 ! if ( v2(i) .le. 1.d-99 ) v2(i) = 0.0d0 if(ik.eq.1)then cm_factor = (dble(618.03/667.75))**2.d0* @ exp( dble(ee*(667.75-618.03)/t(i)) ) vc_factor = dble(667.75/618.03) elseif(ik.eq.2)then cm_factor = 1.d0 vc_factor = 1.d0 elseif(ik.eq.3)then cm_factor = ( dble(720.806/667.75) )**2.d0* @ exp( dble(ee*(667.75-720.806)/t(i)) ) vc_factor = dble(667.75/720.806) end if do j=1,nl ! if (cax1(i,j) .le. 1.d-99 ) cax1(i,j) = 0.0d0 ! if (cax2(i,j) .le. 1.d-99 ) cax2(i,j) = 0.0d0 ! if (cax3(i,j) .le. 1.d-99 ) cax3(i,j) = 0.0d0 c121(i,j) = c121(i,j) + cax1(i,j) * cm_factor cup121(i,j) = cup121(i,j) + cax2(i,j) * cm_factor cdw121(i,j) = cdw121(i,j) + cax3(i,j) * cm_factor end do ! write (*,*) ' i =', i ! write (*,*) ' vc_factor =', vc_factor ! write (*,*) ' v1 =', v1(i) ! write (*,*) ' v2 =', v2(i) ! write (*,*) vc121(i), taugr121(i) ! write (*,*) v1(i) * vc_factor ! write (*,*) vc121(i) + v1(i) * vc_factor vc121(i) = vc121(i) + v1(i) * vc_factor ! write (*,*) v2(i) * vc_factor ! write (*,*) taugr121(i) + v2(i) * vc_factor taugr121(i) = taugr121(i) + v2(i) * vc_factor end do 11 continue ccc end if ccc return end