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