SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad, & dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc,NLEV) !-------------------------------------------------------------------- ! ! Introduction d une routine chimique ! ! Auteurs: S. Lebonnois, 01/2000 | 09/2003 ! adaptation pour Titan 3D: 02/2009 ! adaptation pour // : 04/2013 ! extension chimie jusqu a 1300km : 10/2013 ! ! J. Vatant d'Ollone, 02/2017 ! + adaptation pour le nouveau titan issu du generic ! !--------------------------------------------------------------------- ! use dimphy use datafile_mod, only: datadir use comcstfi_mod, only: rad, pi, kbol use moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat use mod_grid_phy_lmdz, only: nbp_lat implicit none ! Parameters ( dans common_mod in old titan) ! + doivent etre en accord avec titan.h ! ------------------------------------------------- INTEGER,PARAMETER :: NC = 44 ! Nb de composes dans la chimie INTEGER,PARAMETER :: ND = 54 ! Nb de photodissociations INTEGER,PARAMETER :: NR = 377 ! Nb de reactions dans la chimie INTEGER,PARAMETER :: NLRT = 650 ! Pour l'UV 650 niveaux de 2 km ! Dummy parameters without microphysics ! + Just here to keep the whole stuff running without modif C sources ! --------------------------------------------------------------------- INTEGER :: utilaer(16) INTEGER :: aerprod = 0 INTEGER :: htoh2 = 0 ! Arguments ! --------- INTEGER nlon ! nb of horiz points INTEGER ny ! nb de composes (nqmax-nmicro) REAL*8 qy_c(nlon,klev,NC) ! Especes chimiques apres adv.+diss. character*10 nomqy_c(NC+1) ! Noms des especes chimiques REAL*8 declin_rad,ls_rad ! declinaison et long solaire en radians REAL*8 dtchim ! pas de temps chimie REAL*8 ctemp(nlon,klev) ! Temperature REAL*8 cplay(nlon,klev) ! pression (Pa) REAL*8 cplev(nlon,klev+1) ! pression intercouches (Pa) REAL*8 czlay(nlon,klev) ! altitude (m) REAL*8 czlev(nlon,klev+1) ! altitude intercouches (m) REAL*8 dqyc(nlon,klev,NC) ! Tendances especes chimiques INTEGER NLEV ! nbp_lev+70 - a changer si =/ 55 layers ? ! Local variables : ! ----------------- integer i,j,l,ic,jm1 ! variables envoyees dans la chimie: double precision REAL*8 temp_c(NLEV) REAL*8 press_c(NLEV) ! T,p(mbar) a 1 lat donnee REAL*8 cqy(NLEV,NC) ! frac mol qui seront modifiees REAL*8 cqy0(NLEV,NC) ! frac mol avant modif REAL*8 surfhaze(NLEV) REAL*8 cprodaer(NLEV,4),cmaer(NLEV,4) REAL*8 ccsn(NLEV,4),ccsh(NLEV,4) ! rmil: milieu de couche, grille pour K,D,p,ct,T,y ! rinter: intercouche (grille RA dans la chimie) REAL*8 rmil(NLEV),rinter(NLEV),nb(NLEV) REAL*8,save,allocatable :: kedd(:) REAL*8 a,b character str1*1,str2*2 REAL*8 latit character*20 formt1,formt2,formt3 ! variables locales initialisees au premier appel LOGICAL firstcall DATA firstcall/.true./ SAVE firstcall integer dec,declinint,ialt REAL*8 declin_c ! declinaison en degres real*8 factalt,factdec,krpddec,krpddecp1,krpddecm1 real*8 duree REAL*8,save :: mass(NC) REAL*8,save,allocatable :: md(:,:) REAL*8,save :: botCH4 DATA botCH4/0.05/ real*8,save :: r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver REAL*8,save,allocatable :: krpd(:,:,:,:),krate(:,:) integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC) integer,save :: prod(200,NC),perte(2,200,NC) character dummy*30,fich*7,ficdec*15,curdec*15,name*10 real*8 ficalt,oldq,newq,zalt logical okfic !----------------------------------------------------------------------- !*********************************************************************** ! ! Initialisations : ! ---------------- duree = dtchim ! passage en real*4 pour appel a routines C IF (firstcall) THEN print*,'CHIMIE, premier appel' ! ************************************ ! Au premier appel, initialisation de toutes les variables ! pour les routines de la chimie. ! ************************************ allocate(kedd(NLEV)) allocate(krpd(15,ND+1,NLRT,nbp_lat),krate(NLEV,NR),md(NLEV,NC)) ! Verification du nombre de composes: coherence common_mod et nqmax-nmicro ! ---------------------------------- if (ny.ne.NC) then print*,'PROBLEME de coherence dans le nombre de composes:',ny,NC STOP endif ! calcul de temp_c, densites et press_c en moyenne planetaire : ! -------------------------------------------------------------- print*,'pression, densites et temp (init chimie):' print*,'level, press_c, nb, temp_c' DO l=1,klev rinter(l) = (zlevmoy(l)+rad)/1000. rmil(l) = (zlaymoy(l)+rad)/1000. ! temp_c (K): temp_c(l) = tmoy(l) ! press_c (mbar): press_c(l) = playmoy(l)/100. ! nb (cm-3): nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l)) print*,l,rmil(l)-rad/1000.,press_c(l),nb(l),temp_c(l) ENDDO rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000. ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km. do l=klev+2,NLEV rinter(l) = rinter(klev+1) & + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1) rmil(l-1) = (rinter(l-1)+rinter(l))/2. enddo rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2. ! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire ! lecture de tcp.ver, une seule fois ! remplissage r1d,t1d,ct1d,p1d open(11,file=TRIM(datadir)//'/tcp.ver',status='old') read(11,*) dummy do i=1,131 read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i) !print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i) enddo close(11) ! extension pour klev+1 a NLEV avec tcp.ver ! la jonction klev/klev+1 est brutale... Tant pis ? ialt=1 do l=klev+1,NLEV zalt=rmil(l)-rad/1000. do i=ialt,130 if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then ialt=i endif enddo factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt)) press_c(l) = exp( log(p1d(ialt)) * (1-factalt) & + log(p1d(ialt+1)) * factalt ) nb(l) = exp( log(ct1d(ialt)) * (1-factalt) & + log(ct1d(ialt+1)) * factalt ) temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt ! print*,l,zalt,press_c(l),nb(l),temp_c(l) enddo ! caracteristiques des composes: ! ----------------------------- mass(:) = 0.0 call comp(nomqy_c,nb,temp_c,mass,md) print*,' Mass' do ic=1,NC print*,nomqy_c(ic),mass(ic) ! if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic) enddo ! taux de photodissociations: ! -------------------------- call disso(krpd,nbp_lat) ! reactions chimiques: ! ------------------- call chimie(nomqy_c,nb,temp_c,krate,reactif, & nom_perte,nom_prod,perte,prod) ! print*,'nom_prod, nom_perte:' ! do ic=1,NC ! print*,nom_prod(ic),nom_perte(ic) ! enddo ! print*,'premieres prod, perte(1:reaction,2:compagnon):' ! do ic=1,NC ! print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic) ! enddo ! l = klev-3 ! print*,'krate a p=',press_c(l),' reactifs et produits:' ! do ic=1,ND+1 ! print*,ic,krpd(7,ic,l,4)*nb(l)," ", ! . nomqy_c(reactif(1,ic)+1), ! . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), ! . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) ! enddo ! do ic=ND+2,NR ! print*,ic,krate(l,ic)," ", ! . nomqy_c(reactif(1,ic)+1), ! . nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1), ! . nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1) ! enddo ! init kedd ! --------- ! kedd stays constant with time and space ! => linked to pressure rather than altitude... kedd(:) = 5e2 ! valeur mise par defaut ! UNITE ? doit etre ok pour gptitan do l=1,NLEV zalt=rmil(l)-rad/1000. ! z en km if ((zalt.ge.250.).and.(zalt.le.400.)) then kedd(l) = 10.**(3.+(zalt-250.)/50.) ! 1E3 at 250 km ! 1E6 at 400 km elseif ((zalt.gt.400.).and.(zalt.le.600.)) then kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.) ! 2E7 at 600 km elseif ((zalt.gt.600.).and.(zalt.le.900.)) then kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.) ! 1E8 above 900 km elseif ( zalt.gt.900. ) then kedd(l) = 1.e8 endif enddo ENDIF ! premier appel !*********************************************************************** !----------------------------------------------------------------------- ! calcul declin_c (en degres) ! --------------------------- declin_c = declin_rad*180./pi ! print*,'declinaison en degre=',declin_c !*********************************************************************** !*********************************************************************** ! ! BOUCLE SUR LES LATITUDES ! ! * Permet de faire le calcul une seule fois par lat ! DO j=1,nlon if (j.eq.1) then jm1=1 else jm1=j-1 endif if((j.eq.1).or.(klat(j).ne.klat(jm1))) then !*********************************************************************** !*********************************************************************** !----------------------------------------------------------------------- ! ! Distance radiale (intercouches, en km) ! ---------------------------------------- do l=1,klev rinter(l) = (rad+czlev(j,l))/1000. rmil(l) = (rad+czlay(j,l))/1000. ! print*,rinter(l) enddo rinter(klev+1)=(rad+czlev(j,klev+1))/1000. ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km. do l=klev+2,NLEV rinter(l) = rinter(klev+1) & + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1) rmil(l-1) = (rinter(l-1)+rinter(l))/2. enddo rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2. !----------------------------------------------------------------------- ! ! Temperature, pression (mbar), densite (cm-3) ! ------------------------------------------- DO l=1,klev ! temp_c (K): temp_c(l) = ctemp(j,l) ! press_c (mbar): press_c(l) = cplay(j,l)/100. ! nb (cm-3): nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l)) ENDDO ! extension pour klev+1 a NLEV avec tcp.ver ialt=1 do l=klev+1,NLEV zalt=rmil(l)-rad/1000. do i=ialt,130 if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then ialt=i endif enddo factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt)) press_c(l) = exp( log(p1d(ialt)) * (1-factalt) & + log(p1d(ialt+1)) * factalt ) nb(l) = exp( log(ct1d(ialt)) * (1-factalt) & + log(ct1d(ialt+1)) * factalt ) temp_c(l) = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt enddo !----------------------------------------------------------------------- ! ! passage krpd => krate ! --------------------- ! initialisation krate pour dissociations if ((declin_c*10+267).lt.14.) then declinint = 0 dec = 0 else if ((declin_c*10+267).gt.520.) then declinint = 14 dec = 534 else declinint = 1 dec = 27 do while( (declin_c*10+267).ge.real(dec+20) ) dec = dec+40 declinint = declinint+1 enddo endif endif if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then factdec = ( declin_c - (dec-267)/10. ) / 4. else factdec = ( declin_c - (dec-267)/10. ) / 2.7 endif do l=1,NLEV ! INTERPOL EN ALT POUR k (krpd tous les deux km) ialt = int((rmil(l)-rad/1000.)/2.)+1 factalt = (rmil(l)-rad/1000.)/2.-(ialt-1) do i=1,ND+1 krpddecm1 = krpd(declinint ,i,ialt ,klat(j)) * (1-factalt) & + krpd(declinint ,i,ialt+1,klat(j)) * factalt krpddec = krpd(declinint+1,i,ialt ,klat(j)) * (1-factalt) & + krpd(declinint+1,i,ialt+1,klat(j)) * factalt krpddecp1 = krpd(declinint+2,i,ialt ,klat(j)) * (1-factalt) & + krpd(declinint+2,i,ialt+1,klat(j)) * factalt ! ND+1 correspond a la dissociation de N2 par les GCR if ( factdec.lt.0. ) then krate(l,i) = krpddecm1 * abs(factdec) & + krpddec * ( 1 + factdec) endif if ( factdec.gt.0. ) then krate(l,i) = krpddecp1 * factdec & + krpddec * ( 1 - factdec) endif if ( factdec.eq.0. ) krate(l,i) = krpddec enddo enddo !----------------------------------------------------------------------- ! ! composition ! ------------ do ic=1,NC do l=1,klev cqy(l,ic) = qy_c(j,l,ic) cqy0(l,ic) = cqy(l,ic) enddo enddo ! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV WRITE(str2,'(i2.2)') klat(j) fich = "comp_"//str2//".dat" inquire (file=fich,exist=okfic) if (okfic) then open(11,file=fich,status='old') ! premiere ligne=declin read(11,'(A15)') ficdec write(curdec,'(E15.5)') declin_c ! si la declin est la meme: ce fichier a deja ete reecrit ! on lit la colonne t-1 au lieu de la colonne t ! (cas d une bande de latitude partagee par 2 procs) do ic=1,NC read(11,'(A10)') name if (name.ne.nomqy_c(ic)) then print*,"probleme lecture ",fich print*,name,nomqy_c(ic) stop endif if (ficdec.eq.curdec) then do l=klev+1,NLEV read(11,*) ficalt,cqy(l,ic),newq enddo else do l=klev+1,NLEV read(11,*) ficalt,oldq,cqy(l,ic) enddo endif enddo close(11) else ! le fichier n'est pas la do ic=1,NC do l=klev+1,NLEV cqy(l,ic)=cqy(klev,ic) enddo enddo endif cqy0 = cqy !----------------------------------------------------------------------- ! ! Appel de chimietitan ! -------------------- call gptitan(rinter,temp_c,nb, & nomqy_c,cqy, & duree,(klat(j)-1),mass,md, & kedd,botCH4,krate,reactif, & nom_prod,nom_perte,prod,perte, & aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh, & htoh2,surfhaze) ! Tendances composition ! --------------------- do ic=1,NC do l=1,klev dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! en /s enddo enddo !----------------------------------------------------------------------- ! ! sauvegarde compo sur NLEV ! ------------------------- ! dans fichier compo_klat(j) (01 à 49) open(11,file=fich,status='replace') ! premiere ligne=declin write(11,'(E15.5)') declin_c do ic=1,NC write(11,'(A10)') nomqy_c(ic) do l=klev+1,NLEV write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-rad/1000.,cqy0(l,ic),cqy(l,ic) enddo enddo close(11) !*********************************************************************** !*********************************************************************** ! FIN: BOUCLE SUR LES LATITUDES else ! same latitude, we don't do calculations again, only adjust longitudinal variations dqyc(j,:,:) = dqyc(jm1,:,:)/qy_c(jm1,:,:)*qy_c(j,:,:) endif ENDDO !*********************************************************************** !*********************************************************************** firstcall = .false. end subroutine calchim