c*********************************************************************** subroutine NLTEdlvr09_TCOOL (ngridgcm,n_gcm, @ p_gcm, t_gcm, z_gcm, @ co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm, @ q15umco2_gcm ) c jul 2011 malv+fgg c*********************************************************************** implicit none include "dimensions.h" include "dimphys.h" include 'nlte_paramdef.h' include 'nlte_commons.h' include "chimiedata.h" include "conc.h" c Arguments integer n_gcm,ngridgcm real p_gcm(ngridgcm,n_gcm), t_gcm(ngridgcm,n_gcm) real co2vmr_gcm(ngridgcm,n_gcm), n2vmr_gcm(ngridgcm,n_gcm) real covmr_gcm(ngridgcm,n_gcm), o3pvmr_gcm(ngridgcm,n_gcm) real q15umco2_gcm(ngridgcm,n_gcm) real z_gcm(ngridgcm,n_gcm) c local variables and constants integer iz, i, j, k, l, ig,istyle real*8 q15umco2_nl(nl) real*8 zld(nl), zgcmd(n_gcm) real*8 auxdgcm(n_gcm) real p_ig(n_gcm),z_ig(n_gcm) real t_ig(n_gcm) real co2_ig(n_gcm),n2_ig(n_gcm),co_ig(n_gcm),o3p_ig(n_gcm) real mmean_ig(n_gcm),cpnew_ig(n_gcm) c********************************************************************** do ig=1,ngridgcm do l=1,n_gcm p_ig(l)=p_gcm(ig,l) t_ig(l)=t_gcm(ig,l) co2_ig(l)=co2vmr_gcm(ig,l) n2_ig(l)=n2vmr_gcm(ig,l) o3p_ig(l)=o3pvmr_gcm(ig,l) co_ig(l)=covmr_gcm(ig,l) z_ig(l)=z_gcm(ig,l)/1000. mmean_ig(l)=mmean(ig,l) cpnew_ig(l)=cpnew(ig,l) enddo call NLTEdlvr09_ZGRID (n_gcm, @ p_ig, t_ig, z_ig, @ co2_ig,n2_ig,co_ig, $ o3p_ig , mmean_ig, cpnew_ig) c And sets zero to all Curtis Matrixes and Escape Transmissions call leetvt call zero3m (c110,cup110,cdw110, nl) call zero2v (taugr110,vc110, nl) if (itt_cza.eq.24) then call mzescape ( ig,taustar11,tauinf110,tauii110, @ 1, 1,irw_mztf,imu ) istyle=2 call mzescape_normaliz ( taustar11, istyle ) else call mztud (ig, c110,cup110,cdw110, vc110,taugr110, @ 1, 1, irw_mztf, imu, 0,0,0 ) endif call mztvc (ig,vc210, 1, 2, irw_mztf, imu, 0,0,0 ) call mztvc (ig,vc310, 1, 3, irw_mztf, imu, 0,0,0 ) call mztvc (ig,vc410, 1, 4, irw_mztf, imu, 0,0,0 ) call mzescape_fb (ig) input_cza = 0 call NLTEdlvr09_CZALU(ig) if (itt_cza.ne.24) then call mzescape_fh (ig) input_cza = 1 call NLTEdlvr09_CZALU(ig) endif c total cooling rate c smoothing and c interpolation back to original Pgrid c do i = 1, nl q15umco2_nl(i) = hr110(i) + hr210(i) + hr310(i) + hr410(i) @ + hr121(i) enddo do i=1,nl zld(i) = - dble ( alog(pl(i)) ) enddo do i=1,n_gcm zgcmd(i) = - dble( alog(p_gcm(ig,i)) ) enddo call zerov( auxdgcm, n_gcm ) call interdp_limits @ (auxdgcm,zgcmd,n_gcm,jlowerboundary,jtopboundary, @ q15umco2_nl,zld,nl,1,nl,1) call suaviza ( auxdgcm, n_gcm, 1, zgcmd ) do i=1,n_gcm q15umco2_gcm(ig,i) = sngl( auxdgcm(i) ) enddo enddo c end subroutine return end c*********************************************************************** subroutine NLTEdlvr09_ZGRID (n_gcm, @ p_gcm, t_gcm, z_gcm, @ co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm ,mmean_gcm, @ cpnew_gcm) c jul 2011 malv+fgg First version c*********************************************************************** implicit none include "dimensions.h" include "dimphys.h" include 'nlte_paramdef.h' include 'nlte_commons.h' include 'chimiedata.h' include 'conc.h' c Arguments integer n_gcm real p_gcm(n_gcm), t_gcm(n_gcm) real co2vmr_gcm(n_gcm), n2vmr_gcm(n_gcm) real covmr_gcm(n_gcm), o3pvmr_gcm(n_gcm) real z_gcm(n_gcm) real mmean_gcm(n_gcm) real cpnew_gcm(n_gcm) c local variables integer i, j , iz ! real distancia, meanm, gz, Hkm real zmin, zmax, deltazz, deltazzy real nt_gcm(n_gcm) real mmean_nlte(n_gcm),cpnew_nlte(n_gcm) c functions external hrkday_convert real hrkday_convert c*********************************************************************** ! Define working grid for MZ1D model (NL, ZL, ZMIN) ! y otro mas fino para M.Curtis (NZ, ZX, ZXMIN = ZMIN ! Para ello hace falta una z de ref del GCM, que voy a suponer la inferior ! Primero, construimos escala z_gcm ! z_gcm (1) = zmin_gcm ! [km] !write (*,*) ' iz, p, g, H, z =', 1, p_gcm(1), z_gcm(1) ! do iz = 2, n_gcm ! do iz=1,n_gcm ! z_gcm(iz)=zlay(iz)/1.e3 ! meanm = ( co2vmr_gcm(iz)*44. + o3pvmr_gcm(iz)*16. ! @ + n2vmr_gcm(iz)*28. + covmr_gcm(iz)*28. ) ! meanm = meanm / n_avog ! distancia = ( radio + z_gcm(iz-1) )*1.e5 ! gz = gg * masa / ( distancia * distancia ) ! Hkm = 0.5*( t_gcm(iz)+t_gcm(iz-1) ) / ( meanm * gz ) ! Hkm = kboltzman * Hkm *1e-5 ! [km] ! z_gcm(iz) = z_gcm(iz-1) - Hkm * log( p_gcm(iz)/p_gcm(iz-1) ) !write (*,*) iz, p_gcm(iz), gz, Hkm, z_gcm(iz) ! enddo ! Segundo, definimos los límites del modelo, entre las 2 presiones clave ! Bottom boundary for NLTE model : Pbottom=2e-2mb=1.974e-5 atm jlowerboundary = 1 do while ( p_gcm(jlowerboundary) .gt. Pbottom_atm ) jlowerboundary = jlowerboundary + 1 enddo zmin = z_gcm(jlowerboundary) ! write (*,*) ' jlowerboundary, Pmin, zmin =', ! @ jlowerboundary, p_gcm(jlowerboundary), zmin ! Top boundary for NLTE model : Ptop=2e-7mb = 1.974e-5 atm jtopboundary = jlowerboundary do while ( p_gcm(jtopboundary) .gt. Ptop_atm ) jtopboundary = jtopboundary + 1 enddo zmax = z_gcm(jtopboundary) ! write (*,*) ' jtopboundary, Pmax, zmax =', ! @ jtopboundary, p_gcm(jtopboundary),zmax deltaz = (zmax-zmin) / (nl-1) do i=1,nl zl(i) = zmin + (i-1) * deltaz enddo ! write (*,*) ' ZL grid: dz,zmin,zmax ', deltaz, zl(1),zl(nl) ! Creamos el perfil interpolando call intersp ( pl,zl,nl, p_gcm,z_gcm,n_gcm, 2) ! [atm] call intersp ( t,zl,nl, t_gcm,z_gcm,n_gcm, 1) do i = 1, n_gcm nt_gcm(i) = 7.339e+21 * p_gcm(i) / t_gcm(i) ! [cm-3] enddo call intersp ( nt,zl,nl, nt_gcm,z_gcm,n_gcm, 2) call intersp (co2vmr,zl,nl, co2vmr_gcm,z_gcm,n_gcm, 1) call intersp ( n2vmr,zl,nl, n2vmr_gcm,z_gcm,n_gcm, 1) call intersp ( covmr,zl,nl, covmr_gcm,z_gcm,n_gcm, 1) call intersp (o3pvmr,zl,nl, o3pvmr_gcm,z_gcm,n_gcm, 1) call intersp (mmean_nlte,zl,nl,mmean_gcm,z_gcm,n_gcm,1) call intersp (cpnew_nlte,zl,nl,cpnew_gcm,z_gcm,n_gcm,1) do i = 1, nl co2(i) = nt(i) * co2vmr(i) n2(i) = nt(i) * n2vmr(i) co(i) = nt(i) * covmr(i) o3p(i) = nt(i) * o3pvmr(i) ! hrkday_factor(i) = hrkday_convert( t(i), ! @ co2vmr(i), o3pvmr(i), n2vmr(i), covmr(i) ) hrkday_factor(i) = hrkday_convert(mmean_nlte(i) & ,cpnew_nlte(i)) enddo c Fine grid for transmittance calculations deltazy = (zmax-zmin) / (nzy-1) do i=1,nzy zy(i) = zmin + (i-1) * deltazy enddo ! write (*,*) ' ZY grid: nzy,dzy,zmin,zmax ', ! @ nzy, deltazy, zy(1),zy(nzy) call intersp ( py,zy,nzy, p_gcm,z_gcm,n_gcm, 2) ! [atm] call intersp ( ty,zy,nzy, t_gcm,z_gcm,n_gcm, 1) call intersp ( nty,zy,nzy, nt_gcm,z_gcm,n_gcm, 2) call intersp ( co2y,zy,nzy, co2vmr_gcm,z_gcm,n_gcm, 1) do i=1,nzy co2y(i) = co2y(i) * nty(i) enddo c end return end c*********************************************************************** subroutine NLTEdlvr09_CZALU(ig) c jul 2011 malv+fgg c*********************************************************************** implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! common variables and constants include 'nlte_paramdef.h' include 'nlte_commons.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 c*********************************************************************** c hrkday_convert.f c c fortran function that returns the factor for conversion from c hr' [erg s-1 cm-3] to hr [ k day-1 ] c c mar 2010 fgg adapted to GCM c jan 99 malv add o2 as major component. c ago 98 malv also returns cp_avg,pm_avg c jul 98 malv first version. c*********************************************************************** function hrkday_convert @ ( mmean_nlte,cpmean_nlte ) implicit none include 'comcstfi.h' include 'param.h' c argumentos real mmean_nlte,cpmean_nlte real hrkday_convert ccccccccccccccccccccccccccccccccccccc hrkday_convert = daysec * n_avog / & ( cpmean_nlte * 1.e4 * mmean_nlte ) c end return end c*********************************************************************** subroutine sypvvv(a,b,c,d,n) c a(i)=b(i)+c(i)*d(i) c jul 2011 malv+fgg c*********************************************************************** real*8 a(n),b(n),c(n),d(n) integer n,i do 1,i=2,n-1 a(i)= b(i) + c(i) * d(i) 1 continue a(1) = 0.0d0 a(n) = 0.0d0 return end c*********************************************************************** subroutine sypvmv(v,u,c,w,n) c inputs: matriz diagonal c , vectores u,w c output: vector v c Operacion a realizar: v = u + c * w c jul 2011 malv+fgg c*********************************************************************** real*8 v(n),u(n),c(n,n),w(n) integer n,i do 1,i=2,n-1 v(i)= u(i) + c(i,i) * w(i) 1 continue v(1) = 0.0d0 v(n) = 0.0d0 return end c*********************************************************************** subroutine trucommvv(v,b,c,u,w,n) c inputs: matrices b,c , vectores u,w c output: vector v c Operacion a realizar: v = b * c^(-1) * u + w c La matriz c va a ser invertida c c es diagonal, b no c Aprovechamos esa condicion para invertir c, y acelerar el calculo c jul 2011 malv+fgg c*********************************************************************** real*8 v(n),b(n,n),c(n,n),u(n),w(n), sum integer n,i,j,k do 1,i=2,n-1 sum=0.0d0 do 2,j=2,n-1 sum=sum+ (b(i,j)) * (u(j)/c(j,j)) 2 continue v(i) = sum + w(i) 1 continue v(1) = 0.d0 v(n) = 0.d0 return end c*********************************************************************** subroutine trucodiag(a,b,c,d,e,n) c inputs: matrices b,c,d,e c output: matriz diagonal a c Operacion a realizar: a = b * c^(-1) * d + e c La matriz c va a ser invertida c Todas las matrices de entrada son diagonales excepto b c Aprovechamos esa condicion para invertir c, acelerar el calculo, y c ademas, para forzar que a sea diagonal c jul 2011 malv+fgg c*********************************************************************** real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n), sum integer n,i,j,k do 1,i=2,n-1 sum=0.0d0 do 2,j=2,n-1 sum=sum+ (b(i,j)) * (d(j,j)/c(j,j)) 2 continue a(i,i) = sum + e(i,i) 1 continue do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end c*********************************************************************** subroutine invdiag(a,b,n) c inverse of a diagonal matrix c jul 2011 malv c*********************************************************************** implicit none integer n,i,j,k real*8 a(n,n),b(n,n) do 1,i=2,n-1 do 2,j=2,n-1 if (i.eq.j) then a(i,j) = 1.d0/b(i,i) else a(i,j)=0.0d0 end if 2 continue 1 continue do k=1,n a(n,k) = 0.0d0 a(1,k) = 0.0d0 a(k,1) = 0.0d0 a(k,n) = 0.0d0 end do return end