SUBROUTINE calchim(ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim, . ctemp,cplay,cplev, . dqyc) c------------------------------------------------- c c Introduction d une routine chimique c c Auteur: S. Lebonnois, 01/2000 | 09/2003 c adaptation pour Titan 3D: 02/2009 c c------------------------------------------------- c use dimphy implicit none #include "dimensions.h" #include "clesphys.h" #include "paramet.h" #include "YOMCST.h" #include "titan_for.h" !!! doit etre en accord avec titan.h #include "aerprod.h" c Arguments c --------- INTEGER ny ! nb de composes (nqmax-nmicro) REAL qy_c(jjm+1,klev,NC) ! Especes chimiques apres adv.+diss. character*10 nomqy_c(NC+1) ! Noms des especes chimiques REAL declin_rad,ls_rad ! declinaison et long solaire en radians REAL dtchim ! pas de temps chimie REAL ctemp(jjm+1,klev) ! Temperature REAL cplay(jjm+1,klev) ! pression (Pa) REAL cplev(jjm+1,klev) ! pression intercouches (Pa) REAL dqyc(jjm+1,klev,NC) ! Tendances especes chimiques c Local variables : c ----------------- c variables envoyees dans la chimie: double precision integer i,j,l,ic REAL temp_c(klev),press_c(klev) ! T,p(mbar) a 1 lat donnee REAL declin_c ! declinaison en degres REAL cqy(klev,NC) ! frac mol qui seront modifiees REAL surfhaze(klev) REAL cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4) REAL rinter(klev),nb(klev) REAL a,b character str1*1,str2*2 integer ntotftop parameter (ntotftop=50) integer nftop(ntotftop),isaisonflux REAL fluxtop(NC),latit,tabletmp(ntotftop),factflux character*10 nomsftop(ntotftop+1) character*20 formt1,formt2,formt3 c variables locales initialisees au premier appel LOGICAL firstcal DATA firstcal/.true./ SAVE firstcal REAL mass(NC),duree REAL tablefluxtop(NC,jjm+1,5) REAL botCH4 DATA botCH4/0.05/ REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:) integer reactif(5,NR),nom_prod(NC),nom_perte(NC) integer prod(200,NC),perte(2,200,NC) SAVE mass,tablefluxtop,botCH4 SAVE reactif,nom_prod,nom_perte,prod,perte c----------------------------------------------------------------------- c*********************************************************************** c c Initialisations : c ---------------- duree = dtchim ! passage en real*4 pour appel a routines C IF (firstcal) THEN print*,'CHIMIE, premier appel' c ************************************ c Au premier appel, initialisation de toutes les variables c pour les routines de la chimie. c ************************************ allocate(krpd(15,ND+1,klev,jjm+1),krate(klev,NR)) c Verification dimension verticale: coherence titan_for.h et klev c -------------------------------- if (klev.ne.NLEV) then print*,'PROBLEME de coherence dans la dimension verticale:' . ,klev,NLEV STOP endif c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro c ---------------------------------- if (ny.ne.NC) then print*,'PROBLEME de coherence dans le nombre de composes:' . ,ny,NC STOP endif c calcul de temp_c, densites et press_c a l'equateur: c -------------------------------------------------- print*,'pression, densites et temp a l equateur (chimie):' print*,'level, press_c, nb, temp_c' DO l=1,klev c temp_c (K): temp_c(l) = ctemp(jjm/2+1,l) c press_c (mbar): press_c(l) = cplay(jjm/2+1,l)/100. c nb (cm-3): nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l)) print*,l,press_c(l),nb(l),temp_c(l) ENDDO c caracteristiques des composes: c ----------------------------- mass = 0.0 call comp(nomqy_c,mass) print*,' Mass' do ic=1,NC print*,nomqy_c(ic),mass(ic) enddo c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON... c flux dans la derniere couche:(issu du modele 1D, a 500 km) c ----------------------------- c !! en cm-2.s-1 !! c ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54 c c Lecture des tables de flux en fonction lat et saison c ---------------------------------------------------- WRITE(str2(1:2),'(i2.2)') ntotftop formt1 = '(A10,'//str2//'(3X,A10))' formt2 = '(F12.6,'//str2//'(2X,F13.6))' formt3 = '(F13.6,'//str2//'(2X,F13.6))' do j=1,jjp1 do ic=1,NC do l=1,5 tablefluxtop(ic,j,l) = 0. enddo enddo enddo do l=1,4 WRITE(str1(1:1),'(i1.1)') l c hx1 -> Ls=RPI/2 / hx4 -> Ls=0 open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1) c open(10,file="FLUXTOP/flux500.hs"//str1) c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ... read(10,formt1) nomsftop do j=1,ntotftop do i=1,10 !justification a gauche... if (nomsftop(j+1)(i:i).ne.' ') then nomsftop(j) = nomsftop(j+1)(i:) goto 100 endif enddo 100 continue c print*,nomsftop(j) nftop(j) = 0 do i=1,NC if (nomqy_c(i).eq.nomsftop(j)) then nftop(j) = i endif enddo c if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j)) enddo c lecture des flux. Les formats sont un peu alambiques... c a revoir eventuellement do j=1,jjp1/2+1 read(10,formt2) latit,tabletmp do i=1,ntotftop if (nftop(i).ne.0) then tablefluxtop(nftop(i),j,l+1) = tabletmp(i) endif enddo enddo do j=jjp1/2+2,jjp1 read(10,formt3) latit,tabletmp do i=1,ntotftop if (nftop(i).ne.0) then tablefluxtop(nftop(i),j,l+1) = tabletmp(i) endif enddo enddo close(10) enddo ! l do j=1,jjp1 do ic=1,NC tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5) enddo enddo c Stockage des composes utilises dans la prod d'aerosols c (aerprod=1) et dans H-> H2 (htoh2=1): utilaer c ! decalage de 1 car utilise dans le c ! c ------------------------------------------ c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10 c ------------------------------------------ do ic=1,NC c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!! c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c if (nomqy_c(ic).eq."CH4") then c do l=1,llm c do j=1,ip1jmp1 c if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015 c enddo c enddo c endif c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (nomqy_c(ic).eq."C4H2") then utilaer(10) = ic-1 endif if (nomqy_c(ic).eq."HCN") then utilaer(6) = ic-1 do j=1,jjp1 do i=1,5 tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6. enddo enddo endif if (nomqy_c(ic).eq."HC3N") then utilaer(7) = ic-1 endif if (nomqy_c(ic).eq."NCCN") then utilaer(14) = ic-1 endif if (nomqy_c(ic).eq."CH3CN") then utilaer(15) = ic-1 utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle endif if (nomqy_c(ic).eq."H") then utilaer(1) = ic-1 endif if (nomqy_c(ic).eq."H2") then utilaer(2) = ic-1 endif if (nomqy_c(ic).eq."C2H2") then utilaer(3) = ic-1 do j=1,jjp1 do i=1,5 c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3 tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10. enddo enddo endif if (nomqy_c(ic).eq."AC6H6") then utilaer(13) = ic-1 endif if (nomqy_c(ic).eq."C2H3CN") then utilaer(16) = ic-1 endif if (nomqy_c(ic).eq."C2") then utilaer(4) = ic-1 endif if (nomqy_c(ic).eq."C2H") then utilaer(5) = ic-1 endif if (nomqy_c(ic).eq."C3N") then utilaer(8) = ic-1 endif if (nomqy_c(ic).eq."H2CN") then utilaer(9) = ic-1 endif if (nomqy_c(ic).eq."C4H3") then utilaer(11) = ic-1 endif if (nomqy_c(ic).eq."AC6H5") then utilaer(12) = ic-1 endif if (nomqy_c(ic).eq."C2H6") then do j=1,jjp1 do i=1,5 c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640. tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200. enddo enddo endif if (nomqy_c(ic).eq."C3H8") then do j=1,jjp1 do i=1,5 c tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.) tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.) enddo enddo endif if (nomqy_c(ic).eq."C4H10") then do j=1,jjp1 do i=1,5 tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.) enddo enddo endif c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c if ((nomqy_c(ic).eq."HC3N").or. c $ (nomqy_c(ic).eq."C3N")) then c DO j=1,ip1jmp1 c do l=1,34 ! p>~ 1 mbar c qy_c(j,l,ic) = 1.e-30 c enddo c ENDDO c endif c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! enddo c taux de photodissociations: c -------------------------- call disso(krpd,jjp1) c reactions chimiques: c ------------------- call chimie(nomqy_c,nb,temp_c,krate,reactif, . nom_perte,nom_prod,perte,prod) c print*,'nom_prod, nom_perte:' c do ic=1,NC c print*,nom_prod(ic),nom_perte(ic) c enddo c print*,'premieres prod, perte(1:reaction,2:compagnon):' c do ic=1,NC c print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic) c enddo c l = klev-3 c print*,'krate a p=',press_c(l),' reactifs et produits:' c do ic=1,ND+1 c print*,ic,krpd(7,ic,l,4)*nb(l)," ", c . nomqy_c(reactif(1,ic)+1), c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) c enddo c do ic=ND+2,NR c print*,ic,krate(l,ic)," ", c . nomqy_c(reactif(1,ic)+1), c . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), c . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) c enddo ENDIF ! premier appel c*********************************************************************** c----------------------------------------------------------------------- c calcul declin_c (en degres) c --------------------------- declin_c = declin_rad*180./RPI c print*,'declinaison en degre=',declin_c c----------------------------------------------------------------------- c c calcul du facteur d'interpolation entre deux saisons pour flux c -------------------------------------------------------------- isaisonflux = int(ls_rad*2./RPI)+1 if (isaisonflux.eq.5) then isaisonflux = 1 endif factflux = (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/ . (sin(isaisonflux*RPI/2.)-sin((isaisonflux-1)*RPI/2.)) c*********************************************************************** c*********************************************************************** c c BOUCLE SUR LES LATITUDES c DO j=1,jjp1 c*********************************************************************** c*********************************************************************** c----------------------------------------------------------------------- c c Temperature, pression (mbar), densite (cm-3) c ------------------------------------------- DO l=1,klev c temp_c (K): temp_c(l) = ctemp(j,l) c press_c (mbar): press_c(l) = cplay(j,l)/100. c nb (cm-3): nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l)) ENDDO c----------------------------------------------------------------------- c c Distances radiales (intercouches, en km) c ---------------------------------------- rinter(1) = RA/1000. do l=2,klev c A REVOIR: g doit varier avec r ! rinter(l) = rinter(l-1) + . (cplev(j,l-1)-cplev(j,l))/cplay(j,l)*(RD*temp_c(l)/RG)/1000. c print*,rinter(l) enddo c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c FLUX AU TOP: interpolation en saison et c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C) c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do ic=1,NC fluxtop(ic) = tablefluxtop(ic,j,isaisonflux) *(1-factflux) . + tablefluxtop(ic,j,isaisonflux+1)* factflux fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5) enddo c TEST: sortie de l'un des flux avec saveLs et factflux c dans un fichier special pour tracer avec gnuplot c if (j.eq.2) then c open(unit=11,file="flux_80N.txt",status='old',position='append') c write(11,*) ls_rad*180./RPI, factflux, c . ((rinter(llm)-rinter(llm-1))*1e5), fluxtop c write(11,*) ls_rad*180./RPI, factflux, c . fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5), !C2H2 c . fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5) !HCN c close(11) c endif c if (j.eq.17) then c open(unit=11,file="flux_eq.txt",status='old',position='append') c write(11,*) ls_rad*180./RPI, factflux, c . ((rinter(llm)-rinter(llm-1))*1e5), fluxtop c write(11,*) ls_rad*180./RPI, factflux, c . fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5), !C2H2 c . fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5) !HCN c close(11) c endif c FIN TEST: sortie de l'un des flux avec ls et factflux c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (firstcal.and.(j.eq.1)) then print*,'Alt, densites et temp au pole (chimie):' print*,'level, z_bas, nb, temp_c' do l=1,klev print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l) enddo endif if (firstcal.and.(j.eq.jjm/2)) then c print*,'g,mugaz' c print*,g,mugaz print*,'Alt, densites et temp a l equateur (chimie):' print*,'level, z_bas, nb, temp_c' do l=1,klev print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l) enddo endif c----------------------------------------------------------------------- c c composition c ------------ do ic=1,NC do l=1,klev cqy(l,ic) = qy_c(j,l,ic) enddo enddo c----------------------------------------------------------------------- c c total haze area (um2/cm3) c ------------------------- if (htoh2.eq.1) then do l=1,klev ! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!! surfhaze(l) = psurfhaze(j,klev+1-l) c if (j.eq.25) c . print*,'psurfhaze en um2/cm3:',surfhaze(l) enddo endif c----------------------------------------------------------------------- c c Appel de chimietitan c -------------------- call gptitan(jjp1,rinter,temp_c,nb, $ nomqy_c,cqy,fluxtop, $ declin_c,duree,(j-1),mass, $ botCH4,krpd,krate,reactif, $ nom_prod,nom_perte,prod,perte, $ aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, $ htoh2,surfhaze) c if ( j.eq.jjm/2 ) c $ print*,cqy(1,1),cqy(klev,1),cqy(1,2),cqy(klev,2) c if ( j.eq.jjm/2 ) c $ print*,qy_c(j,1,1),qy_c(j,klev,1),qy_c(j,1,2),qy_c(j,klev,2) c stop c Tendances composition c --------------------- do ic=1,NC do l=1,klev dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim ! en /s enddo enddo c----------------------------------------------------------------------- c c production aer c -------------- if (aerprod.eq.1) then do ic=1,4 do l=1,klev prodaer(j,l,ic) = cprodaer(l,ic) maer(j,l,ic) = cmaer(l,ic) csn(j,l,ic) = ccsn(l,ic) csh(j,l,ic) = ccsh(l,ic) enddo enddo endif c*********************************************************************** c*********************************************************************** c c FIN: BOUCLE SUR LES LATITUDES c ENDDO c*********************************************************************** c*********************************************************************** firstcal = .false. RETURN END