c*********************************************************************** subroutine NLTEdlvr09_CZALU(ig) c jul 2011 malv+fgg c*********************************************************************** implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! common variables and constants include 'nltedefs.h' include 'nlte_atm.h' include 'nlte_data.h' include 'nlte_results.h' include 'tcr_15um.h' include 'nlte_matrix.h' include 'nlte_rates.h' include 'nlte_curtis.h' c arguments integer ig !ADDED FOR TRACEBACK c local variables ! matrixes and vectors real*8 e110(nl), e210(nl), e310(nl), e410(nl) real*8 e121(nl), e112(nl) real*8 f1(nl,nl) real*8 cax1(nl,nl), cax2(nl,nl), cax3(nl,nl) real*8 v1(nl), v2(nl), v3(nl) real*8 alf11(nl,nl), alf12(nl,nl) real*8 alf21(nl,nl), alf31(nl,nl), alf41(nl,nl) real*8 a11(nl), a1112(nl,nl) real*8 a1121(nl,nl), a1131(nl,nl), a1141(nl,nl) real*8 a21(nl), a2131(nl,nl), a2141(nl,nl) real*8 a2111(nl,nl), a2112(nl,nl) real*8 a31(nl), a3121(nl,nl), a3141(nl,nl) real*8 a3111(nl,nl), a3112(nl,nl) real*8 a41(nl), a4121(nl,nl), a4131(nl,nl) real*8 a4111(nl,nl), a4112(nl,nl) real*8 a12(nl), a1211(nl,nl) real*8 a1221(nl,nl), a1231(nl,nl), a1241(nl,nl) real*8 aalf11(nl,nl),aalf21(nl,nl),aalf31(nl,nl),aalf41(nl,nl) real*8 aa11(nl), aa1121(nl,nl), aa1131(nl,nl), aa1141(nl,nl) real*8 aa21(nl), aa2111(nl,nl), aa2131(nl,nl), aa2141(nl,nl) real*8 aa31(nl), aa3111(nl,nl), aa3121(nl,nl), aa3141(nl,nl) real*8 aa41(nl), aa4111(nl,nl), aa4121(nl,nl), aa4131(nl,nl) real*8 aa12(nl) real*8 aa1211(nl,nl), aa1221(nl,nl), aa1231(nl,nl), aa1241(nl,nl) real*8 aa1112(nl,nl), aa2112(nl,nl), aa3112(nl,nl), aa4112(nl,nl) real*8 aaalf11(nl,nl),aaalf21(nl,nl),aaalf31(nl,nl),aaalf41(nl,nl) real*8 aaa11(nl),aaa1121(nl,nl),aaa1131(nl,nl),aaa1141(nl,nl) real*8 aaa21(nl),aaa2111(nl,nl),aaa2131(nl,nl),aaa2141(nl,nl) real*8 aaa31(nl),aaa3111(nl,nl),aaa3121(nl,nl),aaa3141(nl,nl) real*8 aaa41(nl),aaa4111(nl,nl),aaa4121(nl,nl),aaa4131(nl,nl) real*8 aaaalf11(nl,nl),aaaalf41(nl,nl) real*8 aaaa11(nl),aaaa1141(nl,nl) real*8 aaaa41(nl),aaaa4111(nl,nl) ! populations real*8 n10(nl), n11(nl) real*8 n20(nl), n21(nl) real*8 n30(nl), n31(nl) real*8 n40(nl), n41(nl) ! productions and loses real*8 d19a1,d19b1,d19c1 real*8 d19ap1,d19bp1,d19cp1 real*8 d19a2,d19b2,d19c2 real*8 d19ap2,d19bp2,d19cp2 real*8 d19a3,d19b3,d19c3 real*8 d19ap3,d19bp3,d19cp3 real*8 d19a4,d19b4,d19c4 real*8 d19ap4,d19bp4,d19cp4 real*8 l11, l12, l21, l31, l41 real*8 p11, p12, p21, p31, p41 real*8 p1112, p1211, p1221, p1231, p1241 real*8 p1121, p1131, p1141 real*8 p2111, p2112, p2131, p2141 real*8 p3111, p3112, p3121, p3141 real*8 p4111, p4112, p4121, p4131 real*8 ps11, ps21, ps31, ps41, ps12 real*8 pl11, pl12, pl21, pl31, pl41 c local constants and indexes integer ii ! decides if output of tv,hr integer icurt ! decides if read/comp c.matrix real*8 co2t real*8 ftest real*8 a11_einst(nl), a12_einst(nl) real*8 a21_einst(nl), a31_einst(nl), a41_einst(nl) real tsurf real*8 nu11, nu12, nu121, nu21, nu31, nu41 integer i, j, ik, isot , icurtishb integer i_by15sh, i_col020, i_col010636 c external functions and subroutines external planckdp real*8 planckdp ! subroutines called: ! mz4sub, dmzout, readc_mz4, mztf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! start program ii = 4 icurt = 1 call zero4v( aa11, aa21, aa31, aa41, nl) call zero4m( aa1121, aa1131, aa1141, aalf11, nl) call zero4m( aa2111, aa2131, aa2141, aalf21, nl) call zero4m( aa3111, aa3121, aa3141, aalf31, nl) call zero4m( aa4111, aa4121, aa4131, aalf41, nl) call zero4m( aa1112, aa2112, aa3112, aa4112, nl) call zero4m( aa1211, aa1221, aa1231, aa1241, nl) call zero3v( aaa41, aaa31, aaa11, nl ) call zero3m( aaa4111, aaa4131, aaalf41, nl) call zero3m( aaa3111, aaa3141, aaalf31, nl) call zero3m( aaa1131, aaa1141, aaalf11, nl) call zero2v( aaaa11, aaaa41, nl ) call zero2m( aaaa1141, aaaalf11, nl) call zero2m( aaaa4111, aaaalf41, nl) !write (*,*) ' --- c z a simple --- input_cza : ', input_cza call zero3v (vt11,vt12,vt13,nl) call zero3v (vt21,vt22,vt23,nl) call zero3v (vt31,vt32,vt33,nl) call zero3v (vt41,vt42,vt43,nl) call zero3v (hr110,hr121,hr132,nl) call zero3v (hr210,hr221,hr232,nl) call zero3v (hr310,hr321,hr332,nl) call zero3v (hr410,hr421,hr432,nl) call zero3v (sl110,sl121,sl132,nl) call zero3v (sl210,sl221,sl232,nl) call zero3v (sl310,sl321,sl332,nl) call zero3v (sl410,sl421,sl432,nl) call zero4v (el11,el21,el31,el41,nl) call zero4v (e110,e210,e310,e410,nl) call zero3v (el12,e121,e112,nl) call zero3m (cax1,cax2,cax3,nl) call zerom (f1,nl) call zero3v (v1,v2,v3,nl) call zero4m (alf11,alf21,alf31,alf41,nl) call zerom (alf12,nl) call zero2v (a11,a12,nl) call zero3v (a21,a31,a41,nl) call zero3m (a1121,a1131,a1141,nl) call zerom (a1112,nl) call zero3m (a1221,a1231,a1241,nl) call zerom (a1211,nl) call zero2m (a2111,a2112,nl) call zero2m (a2131,a2141,nl) call zero2m (a3111,a3112,nl) call zero2m (a3121,a3141,nl) call zero2m (a4111,a4112,nl) call zero2m (a4121,a4131,nl) call zero4v (n11,n21,n31,n41,nl) nu11 = nu(1,1) nu12 = nu(1,2) nu121 = nu12-nu11 nu21 = nu(2,1) nu31 = nu(3,1) nu41 = nu(4,1) ftest = 1.d0 i_by15sh = 1 i_col020 = 1 i_col010636 = 1 101 format(a1) 180 format(a80) c establishing molecular populations needed as input do i=1,nl n10(i) = dble( co2(i) * imr(1) ) n20(i) = dble( co2(i) * imr(2) ) n30(i) = dble( co2(i) * imr(3) ) n40(i) = dble( co2(i) * imr(4) ) if ( input_cza.ge.1 ) then n11(i) = n10(i) *2.d0 *exp( dble(-ee*nu(1,1))/v626t1(i) ) n21(i) = n20(i) *2.d0 *exp( dble(-ee*nu(2,1))/v628t1(i) ) n31(i) = n30(i) *2.d0* exp( dble(-ee*nu(3,1))/v636t1(i) ) n41(i) = n40(i) *2.d0* exp( dble(-ee*nu(4,1))/v627t1(i) ) end if enddo cc cc curtis matrix calculation cc if ( input_cza.ge.1 ) then if (itt_cza.eq.15 ) then call cm15um_hb_simple ( ig,icurt ) elseif (itt_cza.eq.13) then call mztvc_626fh(ig) endif endif do 4,i=nl,1,-1 !---------------------------------------------- co2t = dble ( co2(i) *(imr(1)+imr(3)+imr(2)+imr(4)) ) call getk ( t(i) ) ps11 = 0.d0 ps21 = 0.d0 ps31 = 0.d0 ps41 = 0.d0 ps12 = 0.d0 ! V-T productions and losses V-T isot = 1 d19b1 = dble(k19ba(isot)*co2t+k19bb(isot)*n2(i)) @ + dble(k19bc(isot)*co(i)) d19c1 = dble(k19ca(isot)*co2t+k19cb(isot)*n2(i)) @ + dble(k19cc(isot)*co(i)) d19bp1 = dble( k19bap(isot)*co2t + k19bbp(isot)*n2(i) ) @ + dble( k19bcp(isot)*co(i) ) d19cp1 = dble( k19cap(isot)*co2t + k19cbp(isot)*n2(i) ) @ + dble( k19ccp(isot)*co(i) ) isot = 2 d19c2 =dble(k19ca(isot)*co2t+k19cb(isot)*n2(i)) @ + dble(k19cc(isot)*co(i)) d19cp2 =dble( k19cap(isot)*co2t + k19cbp(isot)*n2(i) ) @ + dble( k19ccp(isot)*co(i) ) isot = 3 d19c3 =dble(k19ca(isot)*co2t+k19cb(isot)*n2(i)) @ + dble(k19cc(isot)*co(i)) d19cp3 =dble( k19cap(isot)*co2t + k19cbp(isot)*n2(i) ) @ + dble( k19ccp(isot)*co(i) ) isot = 4 d19c4 =dble(k19ca(isot)*co2t+k19cb(isot)*n2(i)) @ + dble(k19cc(isot)*co(i)) d19cp4 =dble( k19cap(isot)*co2t + k19cbp(isot)*n2(i) ) @ + dble(k19ccp(isot)*co(i) ) ! l11 = d19c1 + k20c(1)*dble(o3p(i)) p11 = ( d19cp1 + k20cp(1)*dble(o3p(i)) ) * n10(i) l21 = d19c2 + k20c(2)*dble(o3p(i)) p21 = ( d19cp2 + k20cp(2)*dble(o3p(i)) ) *n20(i) l31 = d19c3 + k20c(3)*dble(o3p(i)) p31 = ( d19cp3 + k20cp(3)*dble(o3p(i)) ) *n30(i) l41 = d19c4 + k20c(4)*dble(o3p(i)) p41 = ( d19cp4 + k20cp(4)*dble(o3p(i)) ) *n40(i) ! Addition of V-V l11 = l11 + k21cp(2)*n20(i) + k21cp(3)*n30(i) + k21cp(4)*n40(i) p1121 = k21c(2) * n10(i) p1131 = k21c(3) * n10(i) p1141 = k21c(4) * n10(i) ! l21 = l21 + k21c(2)*n10(i) + k23k21c*n30(i) + k24k21c*n40(i) p2111 = k21cp(2) * n20(i) p2131 = k23k21cp * n20(i) p2141 = k24k21cp * n20(i) ! l31 = l31 + k21c(3)*n10(i) + k23k21cp*n20(i) + k34k21c*n40(i) p3111 = k21cp(3)* n30(i) p3121 = k23k21c * n30(i) p3141 = k34k21cp* n30(i) ! l41 = l41 + k21c(4)*n10(i) + k24k21cp*n20(i) + k34k21cp*n30(i) p4111 = k21cp(4)* n40(i) p4121 = k24k21c * n40(i) p4131 = k34k21c * n40(i) if ( input_cza.ge.1 ) then l12 = d19b1 @ + k20b(1)*dble(o3p(i)) @ + k21b(1)*n10(i) @ + k33c*( n20(i) + n30(i) + n40(i) ) p12 = k21bp(1)*n11(i) * n11(i) p1211 = d19bp1 + k20bp(1)*dble(o3p(i)) p1221 = k33cp(2)*n11(i) p1231 = k33cp(3)*n11(i) p1241 = k33cp(4)*n11(i) l11 = l11 + d19bp1 @ + k20bp(1)*dble(o3p(i)) @ + 2.d0 * k21bp(1) * n11(i) @ + k33cp(2)*n21(i) + k33cp(3)*n31(i) + k33cp(4)*n41(i) p1112 = d19b1 @ + k20b(1)*dble(o3p(i)) @ + 2.d0*k21b(1)*n10(i) @ + k33c*( n20(i) + n30(i) + n40(i) ) l21 = l21 + k33cp(2)*n11(i) p2112 = k33c*n20(i) l31 = l31 + k33cp(3)*n11(i) p3112 = k33c*n30(i) l41 = l41 + k33cp(4)*n11(i) p4112 = k33c*n40(i) end if ! Changes in local losses for ITT=13,15 cases a21_einst(i) = 1.3452d00 * 1.8 / 4.0 * taustar21(i) a31_einst(i) = 1.1878d00 * 1.8 / 4.0 * taustar31(i) a41_einst(i) = 1.2455d00 * 1.8 / 4.0 * taustar41(i) l21 = l21 + a21_einst(i) l31 = l31 + a31_einst(i) l41 = l41 + a41_einst(i) if (input_cza.ge.1 .and. itt_cza.eq.13) then a12_einst(i) = 4.35d00 / 3.0d0 * 1.8 / 4.0 * taustar12(i) l12=l12+a12_einst(i) endif if (itt_cza.eq.24) then a11_einst(i) = a11_einst(i) * 1.8 / 4.0 * taustar11(i) l11 = l11 + a11_einst(i) endif ! vectors and matrices for the formulation a11(i) = dble(gamma*nu11**3.) * 1.d0/2.d0 * (p11+ps11) / @ (n10(i)*l11) a1121(i,i) = dble((nu11/nu21))**3.d0 * n20(i)/n10(i) * p1121/l11 a1131(i,i) = dble((nu11/nu31))**3.d0 * n30(i)/n10(i) * p1131/l11 a1141(i,i) = dble((nu11/nu41))**3.d0 * n40(i)/n10(i) * p1141/l11 e110(i) = 2.d0* dble(vlight*nu11**2.) * 1.d0/2.d0 / @ ( n10(i) * l11 ) a21(i) = dble( gamma*nu21**3.) * 1.d0/2.d0 * @ (p21+ps21)/(n20(i)*l21) a2111(i,i) = dble((nu21/nu11))**3.d0 * n10(i)/n20(i) * p2111/l21 a2131(i,i) = dble((nu21/nu31))**3.d0 * n30(i)/n20(i) * p2131/l21 a2141(i,i) = dble((nu21/nu41))**3.d0 * n40(i)/n20(i) * p2141/l21 e210(i) = 2.d0*dble(vlight*nu21**2.) * 1.d0/2.d0 / @ ( n20(i) * l21 ) a31(i) = dble(gamma*nu31**3.) * 1.d0/2.d0 * (p31+ps31) / @ (n30(i)*l31) a3111(i,i) = dble((nu31/nu11))**3.d0 * n10(i)/n30(i) * p3111/l31 a3121(i,i) = dble((nu31/nu21))**3.d0 * n20(i)/n30(i) * p3121/l31 a3141(i,i) = dble((nu31/nu41))**3.d0 * n40(i)/n30(i) * p3141/l31 e310(i) = 2.d0*dble(vlight*nu31**2.) * 1.d0/2.d0 / @ ( n30(i) * l31 ) a41(i) = dble(gamma*nu41**3.) * 1.d0/2.d0 * (p41+ps41) / @ (n40(i)*l41) a4111(i,i) = dble((nu41/nu11))**3.d0 * n10(i)/n40(i) * p4111/l41 a4121(i,i) = dble((nu41/nu21))**3.d0 * n20(i)/n40(i) * p4121/l41 a4131(i,i) = dble((nu41/nu31))**3.d0 * n30(i)/n40(i) * p4131/l41 e410(i) = 2.d0*dble(vlight*nu41**2.) * 1.d0/2.d0 / @ ( n40(i) * l41 ) if (input_cza.ge.1) then a1112(i,i) = dble((nu11/nu121))**3.d0 * n11(i)/n10(i) * @ p1112/l11 a2112(i,i) = dble((nu21/nu121))**3.d0 * n11(i)/n20(i) * @ p2112/l21 a3112(i,i) = dble((nu31/nu121))**3.d0 * n11(i)/n30(i) * @ p3112/l31 a4112(i,i) = dble((nu41/nu121))**3.d0 * n11(i)/n40(i) * @ p4112/l41 e112(i) = -2.d0*dble(vlight*nu11**3.)/nu121 /2.d0 / @ ( n10(i)*l11 ) a12(i) = dble( gamma*nu121**3.) *2.d0/4.d0* (p12+ps12)/ @ (n11(i)*l12) a1211(i,i) = dble((nu121/nu11))**3.d0 * n10(i)/n11(i) * @ p1211/l12 a1221(i,i) = dble((nu121/nu21))**3.d0 * n20(i)/n11(i) * @ p1221/l12 a1231(i,i) = dble((nu121/nu31))**3.d0 * n30(i)/n11(i) * @ p1231/l12 a1241(i,i) = dble((nu121/nu41))**3.d0 * n40(i)/n11(i) * @ p1241/l12 e121(i) = 2.d0*dble(vlight*nu121**2.) *2.d0/4.d0 / @ ( n11(i) * l12 ) end if 4 continue !------------------------------------------------------- ! Change C.M. do i=1,nl do j=1,nl c210(i,j) = 0.0d0 c310(i,j) = 0.0d0 c410(i,j) = 0.0d0 end do end do if ( itt_cza.eq.13 ) then do i=1,nl do j=1,nl c121(i,j) = 0.0d0 end do end do endif !Añadido para hacer diagonal C121 ! if ( itt_cza.eq.15 ) then ! do i=1,nl ! do j=1,nl ! if(abs(i-j).eq.1.or.abs(i-j).eq.2) c121(i,j) = 0.0d0 ! end do ! end do ! endif if ( itt_cza.eq.24 ) then do i=1,nl do j=1,nl c110(i,j) = 0.0d0 end do end do endif ! Lower Boundary tsurf = t(1) + tsurf_excess do i=1,nl sl110(i) = sl110(i) + vc110(i) * planckdp( tsurf, nu11 ) sl210(i) = sl210(i) + vc210(i) * planckdp( tsurf, nu21 ) sl310(i) = sl310(i) + vc310(i) * planckdp( tsurf, nu31 ) sl410(i) = sl410(i) + vc410(i) * planckdp( tsurf, nu41 ) end do if (input_cza.ge.1) then do i=1,nl sl121(i) = sl121(i) + vc121(i) * planckdp( tsurf, nu121 ) end do endif !!!!!!!!!!!! Solucion del sistema !! Paso 0 : Calculo de los alphas alf11, alf21, alf31, alf41, alf12 call unit ( cax2, nl ) call diago ( cax1, e110, nl ) call mulmm ( cax3, cax1,c110, nl ) ! cax3=matmul(cax1,c110) call resmm ( alf11, cax2,cax3, nl ) call diago ( cax1, e210, nl ) call mulmm ( cax3, cax1,c210, nl ) ! cax3=matmul(cax1,c210) call resmm ( alf21, cax2,cax3, nl ) call diago ( cax1, e310, nl ) call mulmm ( cax3, cax1,c310, nl ) ! cax3=matmul(cax1,c310) call resmm ( alf31, cax2,cax3, nl ) ! call diago ( cax1, e410, nl ) call mulmm ( cax3, cax1,c410, nl ) ! cax3=matmul(cax1,c410) call resmm ( alf41, cax2,cax3, nl ) ! ! if(ig.eq.2223.and.input_cza.eq.1) then ! open(168,file='output_curtis_c121diagminus2.dat') ! do i=1,nl ! do j=1,nl ! write(168,*)i,j,c110(i,j),c121(i,j) ! enddo ! enddo ! close(168) ! open(178,file='output_taustar.dat') ! do i=1,nl ! write(178,*)i,taustar21(i),taustar31(i),taustar41(i) ! enddo ! close(178) ! endif if (input_cza.ge.1) then call diago ( cax1, e121, nl ) call mulmm ( cax3, cax1,c121, nl ) ! cax3=matmul(cax1,c121) call resmm ( alf12, cax2,cax3, nl ) endif !! Paso 1 : Calculo de vectores y matrices con 1 barra (aa***) if (input_cza.eq.0) then ! Skip paso 1, pues el12 no se calcula ! el11 call sypvvv( aa11, a11,e110,sl110, nl ) call samem( aa1121, a1121, nl ) call samem( aa1131, a1131, nl ) call samem( aa1141, a1141, nl ) call samem( aalf11, alf11, nl ) ! el21 call sypvvv( aa21, a21,e210,sl210, nl ) call samem( aa2111, a2111, nl ) call samem( aa2131, a2131, nl ) call samem( aa2141, a2141, nl ) call samem( aalf21, alf21, nl ) ! el31 call sypvvv( aa31, a31,e310,sl310, nl ) call samem( aa3111, a3111, nl ) call samem( aa3121, a3121, nl ) call samem( aa3141, a3141, nl ) call samem( aalf31, alf31, nl ) ! el41 call sypvvv( aa41, a41,e410,sl410, nl ) call samem( aa4111, a4111, nl ) call samem( aa4121, a4121, nl ) call samem( aa4131, a4131, nl ) call samem( aalf41, alf41, nl ) else ! (input_cza.ge.1) , FH ! call sypvvv( v1, a12,e121,sl121, nl ) ! a12 + e121 * sl121 ! aa11 call sypvvv( v2, a11,e110,sl110, nl ) call trucommvv( aa11 , alf12,a1112,v2, v1, nl ) ! aalf11 call invdiag( cax1, a1112, nl ) call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a1112) ! cax2=matmul(alf12,cax1) call mulmm( cax3, cax2, alf11, nl ) ! cax3=matmul(cax2,alf11) call resmm( aalf11, cax3, a1211, nl ) ! aa1121 call trucodiag(aa1121, alf12,a1112,a1121, a1221, nl) ! aa1131 call trucodiag(aa1131, alf12,a1112,a1131, a1231, nl) ! aa1141 call trucodiag(aa1141, alf12,a1112,a1141, a1241, nl) ! aa21 call sypvvv( v2, a21,e210,sl210, nl ) call trucommvv( aa21 , alf12,a2112,v2, v1, nl ) ! aalf21 call invdiag( cax1, a2112, nl ) call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a2112) ! cax2=matmul(alf12,cax1) call mulmm( cax3, cax2, alf21, nl ) ! cax3=matmul(cax2,alf21) call resmm( aalf21, cax3, a1221, nl ) ! aa2111 call trucodiag(aa2111, alf12,a2112,a2111, a1211, nl) ! aa2131 call trucodiag(aa2131, alf12,a2112,a2131, a1231, nl) ! aa2141 call trucodiag(aa2141, alf12,a2112,a2141, a1241, nl) ! aa31 call sypvvv( v2, a31,e310,sl310, nl ) call trucommvv( aa31 , alf12,a3112,v2, v1, nl ) ! aalf31 call invdiag( cax1, a3112, nl ) call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a3112) ! cax2=matmul(alf12,cax1) call mulmm( cax3, cax2, alf31, nl ) ! cax3=matmul(cax2,alf31) call resmm( aalf31, cax3, a1231, nl ) ! aa3111 call trucodiag(aa3111, alf12,a3112,a3111, a1211, nl) ! aa3121 call trucodiag(aa3121, alf12,a3112,a3121, a1221, nl) ! aa3141 call trucodiag(aa3141, alf12,a3112,a3141, a1241, nl) ! aa41 call sypvvv( v2, a41,e410,sl410, nl ) call trucommvv( aa41 , alf12,a4112,v2, v1, nl ) ! aalf41 call invdiag( cax1, a4112, nl ) call mulmm( cax2, alf12, cax1, nl ) ! alf12 * (1/a4112) ! cax2=matmul(alf12,cax1) call mulmm( cax3, cax2, alf41, nl ) ! cax3=matmul(cax2,alf41) call resmm( aalf41, cax3, a1241, nl ) ! aa4111 call trucodiag(aa4111, alf12,a4112,a4111, a1211, nl) ! aa4121 call trucodiag(aa4121, alf12,a4112,a4121, a1221, nl) ! aa4131 call trucodiag(aa4131, alf12,a4112,a4131, a1231, nl) endif ! Final caso input_cza.ge.1 !! Paso 2 : Calculo de vectores y matrices con 2 barras (aaa***) ! aaalf41 call invdiag( cax1, aa4121, nl ) call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a4121) ! cax2=matmul(aalf21,cax1) call mulmm( cax3, cax2, aalf41, nl ) ! cax3=matmul(cax2,aalf41) call resmm( aaalf41, cax3, aa2141, nl ) ! aaa41 call trucommvv(aaa41, aalf21,aa4121,aa41, aa21, nl) ! aaa4111 call trucodiag(aaa4111, aalf21,aa4121,aa4111, aa2111, nl) ! aaa4131 call trucodiag(aaa4131, aalf21,aa4121,aa4131, aa2131, nl) ! aaalf31 call invdiag( cax1, aa3121, nl ) call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a3121) ! cax2=matmul(aalf21,cax1) call mulmm( cax3, cax2, aalf31, nl ) ! cax3=matmul(cax2,aalf31) call resmm( aaalf31, cax3, aa2131, nl ) ! aaa31 call trucommvv(aaa31, aalf21,aa3121,aa31, aa21, nl) ! aaa3111 call trucodiag(aaa3111, aalf21,aa3121,aa3111, aa2111, nl) ! aaa3141 call trucodiag(aaa3141, aalf21,aa3121,aa3141, aa2141, nl) ! aaalf11 call invdiag( cax1, aa1121, nl ) call mulmm( cax2, aalf21, cax1, nl ) ! alf21 * (1/a1121) ! cax2=matmul(aalf21,cax1) call mulmm( cax3, cax2, aalf11, nl ) ! cax3=matmul(cax2,aalf11) call resmm( aaalf11, cax3, aa2111, nl ) ! aaa11 call trucommvv(aaa11, aalf21,aa1121,aa11, aa21, nl) ! aaa1131 call trucodiag(aaa1131, aalf21,aa1121,aa1131, aa2131, nl) ! aaa1141 call trucodiag(aaa1141, aalf21,aa1121,aa1141, aa2141, nl) !! Paso 3 : Calculo de vectores y matrices con 3 barras (aaaa***) ! aaaalf41 call invdiag( cax1, aaa4131, nl ) call mulmm( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131) ! cax2=matmul(aaalf31,cax1) call mulmm( cax3, cax2, aaalf41, nl ) ! cax3=matmul(cax2,aaalf41) call resmm( aaaalf41, cax3, aaa3141, nl ) ! aaaa41 call trucommvv(aaaa41, aaalf31,aaa4131,aaa41, aaa31, nl) ! aaaa4111 call trucodiag(aaaa4111, aaalf31,aaa4131,aaa4111,aaa3111, nl) ! aaaalf11 call invdiag( cax1, aaa1131, nl ) call mulmm( cax2, aaalf31, cax1, nl ) ! aaalf31 * (1/aaa4131) ! cax2=matmul(aaalf31,cax1) call mulmm( cax3, cax2, aaalf11, nl ) ! cax3=matmul(cax2,aaalf11) call resmm( aaaalf11, cax3, aaa3111, nl ) ! aaaa11 call trucommvv(aaaa11, aaalf31,aaa1131,aaa11, aaa31, nl) ! aaaa1141 call trucodiag(aaaa1141, aaalf31,aaa1131,aaa1141,aaa3141, nl) !! Paso 4 : Calculo de vectores y matrices finales y calculo de J1 call trucommvv(v1, aaaalf41,aaaa1141,aaaa11, aaaa41, nl) ! call invdiag( cax1, aaaa1141, nl ) call mulmm( cax2, aaaalf41, cax1, nl ) ! aaaalf41 * (1/aaaa1141) ! cax2=matmul(aaaalf41,cax1) call mulmm( cax3, cax2, aaaalf11, nl ) ! cax3=matmul(cax2,aaaalf11) call resmm( cax1, cax3, aaaa4111, nl ) ! call LUdec ( el11, cax1, v1, nl, nl2 ) ! Solucion para el41 call sypvmv( v1, aaaa41, aaaa4111,el11, nl ) call LUdec ( el41, aaaalf41, v1, nl, nl2 ) ! Solucion para el31 call sypvmv( v2, aaa31, aaa3111,el11, nl ) call sypvmv( v1, v2, aaa3141,el41, nl ) call LUdec ( el31, aaalf31, v1, nl, nl2 ) ! Solucion para el21 call sypvmv( v3, aa21, aa2111,el11, nl ) call sypvmv( v2, v3, aa2131,el31, nl ) call sypvmv( v1, v2, aa2141,el41, nl ) call LUdec ( el21, aalf21, v1, nl, nl2 ) !!! el11(1) = planckdp( t(1), nu11 ) el21(1) = planckdp( t(1), nu21 ) el31(1) = planckdp( t(1), nu31 ) el41(1) = planckdp( t(1), nu41 ) el11(nl) = 2.d0 * el11(nl-1) - el11(nl2) el21(nl) = 2.d0 * el21(nl-1) - el21(nl2) el31(nl) = 2.d0 * el31(nl-1) - el31(nl2) el41(nl) = 2.d0 * el41(nl-1) - el41(nl2) call mulmv ( v1, c110,el11, nl ) call sumvv ( hr110, v1,sl110, nl ) ! Solucion para el12 if (input_cza.ge.1) then call sypvmv( v1, a12, a1211,el11, nl ) call sypvmv( v3, v1, a1221,el21, nl ) call sypvmv( v2, v3, a1231,el31, nl ) call sypvmv( v1, v2, a1241,el41, nl ) call LUdec ( el12, alf12, v1, nl, nl2 ) el12(1) = planckdp( t(1), nu121 ) el12(nl) = 2.d0 * el12(nl-1) - el12(nl2) if (itt_cza.eq.15) then call mulmv ( v1, c121,el12, nl ) call sumvv ( hr121, v1,sl121, nl ) endif end if if (input_cza.lt.1) then do i=1,nl pl11 = el11(i)/dble( gamma * nu11**3.0d0 * 1./2. / n10(i) ) pl21 = el21(i)/dble( gamma * nu21**3.0d0 * 1./2. / n20(i) ) pl31 = el31(i)/dble( gamma * nu31**3.0d0 * 1./2. / n30(i) ) pl41 = el41(i)/dble( gamma * nu41**3.0d0 * 1./2. / n40(i) ) vt11(i) = dble(-ee*nu11) / log( abs(pl11) / (2.0d0*n10(i)) ) vt21(i) = dble(-ee*nu21) / log( abs(pl21) / (2.0d0*n20(i)) ) vt31(i) = dble(-ee*nu31) / log( abs(pl31) / (2.0d0*n30(i)) ) vt41(i) = dble(-ee*nu41) / log( abs(pl41) / (2.0d0*n40(i)) ) hr210(i) = sl210(i) - hplanck*vlight*nu21 * a21_einst(i)*pl21 hr310(i) = sl310(i) - hplanck*vlight*nu31 * a31_einst(i)*pl31 hr410(i) = sl410(i) - hplanck*vlight*nu41 * a41_einst(i)*pl41 ! hr410(i) = 0. enddo call dinterconnection ( v626t1, vt11 ) call dinterconnection ( v628t1, vt21 ) call dinterconnection ( v636t1, vt31 ) call dinterconnection ( v627t1, vt41 ) else do i=1,nl pl21 = el21(i)/dble( gamma * nu21**3.0d0 * 1./2. / n20(i) ) pl31 = el31(i)/dble( gamma * nu31**3.0d0 * 1./2. / n30(i) ) pl41 = el41(i)/dble( gamma * nu41**3.0d0 * 1./2. / n40(i) ) hr210(i) = sl210(i) - hplanck*vlight*nu21 * a21_einst(i)*pl21 hr310(i) = sl310(i) - hplanck*vlight*nu31 * a31_einst(i)*pl31 hr410(i) = sl410(i) - hplanck*vlight*nu41 * a41_einst(i)*pl41 ! hr410(i) = 0. if (itt_cza.eq.13) then pl12 = el12(i)/dble( gamma * nu121**3.0d0 * 2./4. / n11(i) ) hr121(i) = - hplanck*vlight * nu121 * a12_einst(i) * pl12 hr121(i) = hr121(i) + sl121(i) endif enddo endif ! K/Dday do i=1,nl hr110(i)=hr110(i)*( hrkday_factor(i) / nt(i) ) hr210(i)=hr210(i)*( hrkday_factor(i) / nt(i) ) hr310(i)=hr310(i)*( hrkday_factor(i) / nt(i) ) hr410(i)=hr410(i)*( hrkday_factor(i) / nt(i) ) hr121(i)=hr121(i)*( hrkday_factor(i) / nt(i) ) end do c output !codigo = codeout !call dmzout_tv ( 1 ) !call dmzout_hr ( 1 ) c final subrutina return end