Changeset 2433 for trunk/LMDZ.MARS/libf/aeronomars
- Timestamp:
- Nov 14, 2020, 4:51:23 PM (4 years ago)
- Location:
- trunk/LMDZ.MARS/libf/aeronomars
- Files:
-
- 4 edited
-
calchim_mod.F90 (modified) (1 diff)
-
photochemistry.F90 (modified) (11 diffs)
-
photolysis_mod.F90 (modified) (4 diffs)
-
photolysis_online.F (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
r2321 r2433 718 718 719 719 call photochemistry(nlayer,nq,nbq,ionchem,nb_reaction_3_max, & 720 nb_reaction_4_max, nb_phot_max, nphotion,&721 jonline, ig,lswitch,zycol,szacol,ptimestep,&720 nb_reaction_4_max,nphot,nb_phot_max,nphotion, & 721 jonline,ig,lswitch,zycol,szacol,ptimestep, & 722 722 zpress,zlocal,ztemp,ztelec,zdens,zmmean, & 723 723 dist_sol,zday,surfdust1d,surfice1d, & -
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
r2321 r2433 3 3 ! Photochemical routine 4 4 ! 5 ! Author : Franck Lefevre6 ! ------ 5 ! Authors: Franck Lefevre, Francisco Gonzalez-Galindo 6 ! ------- 7 7 ! 8 ! Version: 27/04/20178 ! Version: 14/11/2020 9 9 ! 10 10 ! ASIS scheme : for details on the method see … … 14 14 15 15 subroutine photochemistry(nlayer, nq, nesp, ionchem, nb_reaction_3_max, & 16 nb_reaction_4_max, n b_phot_max, nphotion,&16 nb_reaction_4_max, nphot, nb_phot_max, nphotion, & 17 17 jonline, ig, lswitch, zycol, sza, ptimestep, press, & 18 18 alt, temp, temp_elect, dens, zmmean, & … … 38 38 integer, intent(in) :: nb_reaction_4_max 39 39 ! number of bimolecular reactions 40 integer, intent(in) :: nphot 41 ! number of photodissociations 40 42 integer, intent(in) :: nb_phot_max 41 43 ! number of reactions treated numerically as photodissociations … … 121 123 integer,parameter :: i_elec = 32 122 124 125 !Tracer indexes for photionization coeffs 126 127 integer,parameter :: induv_co2 = 1 128 integer,parameter :: induv_o2 = 2 129 integer,parameter :: induv_o = 3 130 integer,parameter :: induv_h2o = 4 131 integer,parameter :: induv_h2 = 5 132 integer,parameter :: induv_h2o2= 6 133 integer,parameter :: induv_o3 = 7 134 integer,parameter :: induv_n2 = 8 135 integer,parameter :: induv_n = 9 136 integer,parameter :: induv_no = 10 137 integer,parameter :: induv_co = 11 138 integer,parameter :: induv_h = 12 139 integer,parameter :: induv_no2 = 13 140 123 141 integer :: ilay 124 142 … … 193 211 tau, sza, dist_sol, v_phot) 194 212 213 !Calculation of photoionization rates, if needed 195 214 if (jionos .and. ionchem) then 196 215 call jthermcalc_e107(ig,nlayer,2,c,nesp,temp,alt,sza,zday) 197 216 do ilay=1,lswitch-1 198 217 call phdisrate(ig,nlayer,2,sza,ilay) 199 enddo 200 v_phot(:,14)=jion(1,:,1) 201 v_phot(:,15)=jion(1,:,2) 202 v_phot(:,16)=jion(1,:,2) 203 v_phot(:,17)=jion(1,:,3) 204 v_phot(:,18)=jion(1,:,3) 205 v_phot(:,19)=jion(1,:,4) 206 v_phot(:,20)=jion(1,:,4) 207 v_phot(:,21)=jion(2,:,1) 208 v_phot(:,22)=jion(3,:,1) 209 v_phot(:,23)=jion(10,:,1) 210 v_phot(:,24)=jion(11,:,1) 211 v_phot(:,25)=jion(11,:,2) 212 v_phot(:,26)=jion(11,:,2) 213 v_phot(:,27)=jion(8,:,1) 214 v_phot(:,28)=jion(8,:,2) 215 v_phot(:,29)=jion(8,:,2) 216 v_phot(:,30)=jion(9,:,1) 217 v_phot(:,31)=jion(12,:,1) 218 endif 218 end do 219 !CO2 photoionization 220 v_phot(:,nphot+ 1) = jion(induv_co2,:,1) 221 v_phot(:,nphot+ 2) = jion(induv_co2,:,2) 222 v_phot(:,nphot+ 3) = jion(induv_co2,:,2) 223 v_phot(:,nphot+ 4) = jion(induv_co2,:,3) 224 v_phot(:,nphot+ 5) = jion(induv_co2,:,3) 225 v_phot(:,nphot+ 6) = jion(induv_co2,:,4) 226 v_phot(:,nphot+ 7) = jion(induv_co2,:,4) 227 !O2 photoionization 228 v_phot(:,nphot+ 8) = jion(induv_o2,:,1) 229 !O photoionization 230 v_phot(:,nphot+ 9) = jion(induv_o,:,1) 231 !NO photoionization 232 v_phot(:,nphot+10) = jion(induv_no,:,1) 233 !CO photoionization 234 v_phot(:,nphot+11) = jion(induv_co,:,1) 235 v_phot(:,nphot+12) = jion(induv_co,:,2) 236 v_phot(:,nphot+13) = jion(induv_co,:,2) 237 !N2 photoionization 238 v_phot(:,nphot+14) = jion(induv_n2,:,1) 239 v_phot(:,nphot+15) = jion(induv_n2,:,2) 240 v_phot(:,nphot+16) = jion(induv_n2,:,2) 241 !N photoionization 242 v_phot(:,nphot+17) = jion(induv_n,:,1) 243 !H photoionization 244 v_phot(:,nphot+18) = jion(induv_h,:,1) 245 end if 219 246 220 247 else ! night … … 277 304 278 305 hetero_dust = .false. 279 hetero_ice = . true.306 hetero_ice = .false. 280 307 281 308 call reactionrates(nlayer, ionchem, nb_reaction_3_max, nb_reaction_4_max, & … … 758 785 c002(:) = 1.8e-11*exp(180./t(:)) 759 786 787 ! c002(:) = c002(:)*1.12 ! FL li et al. 2007 788 760 789 ! robertson and smith, j. chem. phys. a 110, 6673, 2006 761 790 … … 809 838 810 839 c007(:) = 4.8e-11*exp(250./t(:)) 840 841 ! c007(:) = c007(:)*0.9 ! FL li et al. 2007 811 842 812 843 nb_reaction_4 = nb_reaction_4 + 1 … … 851 882 852 883 do ilev = 1,lswitch-1 884 ! ak0 = 3.1*2.4*4.4e-32*(t(ilev)/300.)**(-1.3) ! FL li et al 2017 853 885 ak0 = 2.4*4.4e-32*(t(ilev)/300.)**(-1.3) 854 886 ak1 = 7.5e-11*(t(ilev)/300.)**(0.2) … … 1058 1090 ! jpl 2015 1059 1091 1092 do ilev = 1,lswitch-1 1093 1094 ! branch 1 : oh + co -> h + co2 1095 1096 rate1 = 1.5e-13*(t(ilev)/300.)**(0.0) 1097 1098 ! branch 2 : oh + co + m -> hoco + m 1099 1100 ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0) 1101 ak1 = 1.1e-12*(t(ilev)/300.)**(1.3) 1102 rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1) 1103 xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2) 1104 1105 e001(ilev) = rate1 + rate2*0.6**xpo 1106 end do 1107 1108 ! joshi et al., 2006 1109 1060 1110 ! do ilev = 1,lswitch-1 1061 1062 ! branch 1 : oh + co -> h + co2 1063 1064 ! rate1 = 1.5e-13*(t(ilev)/300.)**(0.0) 1065 1066 ! branch 2 : oh + co + m -> hoco + m 1067 1068 ! ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0) 1069 ! ak1 = 1.1e-12*(t(ilev)/300.)**(1.3) 1070 ! rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1) 1071 ! xpo = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2) 1072 1073 ! e001(ilev) = rate1 + rate2*0.6**xpo 1111 ! k1a0 = 1.34*2.5*dens(ilev) & 1112 ! *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev))) & 1113 ! + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev)))) ! typo in paper corrected 1114 ! k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev)) & 1115 ! + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev)) 1116 ! k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev)) & 1117 ! + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev)) 1118 ! x = k1a0/(k1ainf - k1b0) 1119 ! y = k1b0/(k1ainf - k1b0) 1120 ! fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev)) & 1121 ! + exp(-t(ilev)/255.) 1122 ! fx = fc**(1./(1. + (alog(x))**2)) ! typo in paper corrected 1123 ! k1a = k1a0*((1. + y)/(1. + x))*fx 1124 ! k1b = k1b0*(1./(1.+x))*fx 1125 ! e001(ilev) = k1a + k1b 1074 1126 ! end do 1075 1076 ! joshi et al., 20061077 1078 do ilev = 1,lswitch-11079 k1a0 = 1.34*2.5*dens(ilev) &1080 *1/(1/(3.62e-26*t(ilev)**(-2.739)*exp(-20./t(ilev))) &1081 + 1/(6.48e-33*t(ilev)**(0.14)*exp(-57./t(ilev)))) ! typo in paper corrected1082 k1b0 = 1.17e-19*t(ilev)**(2.053)*exp(139./t(ilev)) &1083 + 9.56e-12*t(ilev)**(-0.664)*exp(-167./t(ilev))1084 k1ainf = 1.52e-17*t(ilev)**(1.858)*exp(28.8/t(ilev)) &1085 + 4.78e-8*t(ilev)**(-1.851)*exp(-318./t(ilev))1086 x = k1a0/(k1ainf - k1b0)1087 y = k1b0/(k1ainf - k1b0)1088 fc = 0.628*exp(-1223./t(ilev)) + (1. - 0.628)*exp(-39./t(ilev)) &1089 + exp(-t(ilev)/255.)1090 fx = fc**(1./(1. + (alog(x))**2)) ! typo in paper corrected1091 k1a = k1a0*((1. + y)/(1. + x))*fx1092 k1b = k1b0*(1./(1.+x))*fx1093 e001(ilev) = k1a + k1b1094 end do1095 1127 1096 1128 nb_reaction_4 = nb_reaction_4 + 1 … … 2016 2048 2017 2049 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2d, 1.0, i_n) 2050 2051 !=========================================================== 2052 ! HDO + hv -> products 2053 !=========================================================== 2054 2055 !nb_phot = nb_phot + 1 2056 2057 !indice_phot(nb_phot) = z3spec(1.0, i_hdo, 0.0, i_dummy, 0.0, i_dummy) 2018 2058 2019 2059 !Only if ion chemistry included -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_mod.F90
r2170 r2433 6 6 7 7 integer, parameter :: nphot = 13 ! number of photolysis 8 ! integer, parameter :: nphot = 14 ! number of photolysis (with hdo) 8 9 integer, parameter :: nabs = 10 ! number of absorbing gases 9 10 … … 37 38 real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2) 38 39 real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield 40 real, dimension(nw), save :: xshdo ! hdo absorption cross-section (cm2) 39 41 real, dimension(nw), save :: albedo ! surface albedo 40 42 … … 97 99 98 100 call rdxsn2(nw,wl,xsn2,yieldn2) 101 102 ! read and grid hdo cross-sections 103 104 call rdxshdo(nw,wl,xshdo) 99 105 100 106 ! set surface albedo … … 1247 1253 !============================================================================== 1248 1254 1255 subroutine rdxshdo(nw, wl, yg) 1256 1257 !-----------------------------------------------------------------------------* 1258 != PURPOSE: =* 1259 != Read HDO molecular absorption cross section. Re-grid data to match =* 1260 != specified wavelength working grid. =* 1261 !-----------------------------------------------------------------------------* 1262 != PARAMETERS: =* 1263 != NW - INTEGER, number of specified intervals + 1 in working (I)=* 1264 != wavelength grid =* 1265 != WL - REAL, vector of lower limits of wavelength intervals in (I)=* 1266 != working wavelength grid =* 1267 != YG - REAL, molecular absoprtion cross section (cm^2) of HDO at (O)=* 1268 != each specified wavelength =* 1269 !-----------------------------------------------------------------------------* 1270 1271 use datafile_mod, only: datadir 1272 1273 IMPLICIT NONE 1274 1275 ! input 1276 1277 integer :: nw ! number of wavelength grid points 1278 real, dimension(nw) :: wl ! lower and central wavelength for each interval 1279 1280 ! output 1281 1282 real, dimension(nw) :: yg ! hdo cross-sections (cm2) 1283 1284 ! local 1285 1286 integer, parameter :: kdata = 900 1287 real, parameter :: deltax = 1.e-4 1288 REAL x1(kdata) 1289 REAL y1(kdata) 1290 INTEGER ierr 1291 INTEGER i, n 1292 CHARACTER*100 fil 1293 integer :: kin, kout ! input/output logical units 1294 1295 kin = 10 1296 1297 fil = trim(datadir)//'/cross_sections/hdo_composite_295K.txt' 1298 print*, 'section efficace HDO: ', fil 1299 1300 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1301 1302 DO i = 1,17 1303 read(kin,*) 1304 END DO 1305 1306 n = 806 1307 DO i = 1, n 1308 READ(kin,*) x1(i), y1(i) 1309 END DO 1310 CLOSE (kin) 1311 1312 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 1313 CALL addpnt(x1,y1,kdata,n, 0.,0.) 1314 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 1315 CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) 1316 CALL inter2(nw,wl,yg,n,x1,y1,ierr) 1317 IF (ierr .NE. 0) THEN 1318 WRITE(*,*) ierr, fil 1319 STOP 1320 ENDIF 1321 1322 end subroutine rdxshdo 1323 1324 !============================================================================== 1325 1249 1326 subroutine rdxsh2o2(nw, wl, xsh2o2) 1250 1327 -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
r2170 r2433 60 60 61 61 integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o, 62 $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2 62 $ j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2, j_hdo 63 63 64 64 integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, … … 96 96 j_no2 = 12 ! no2 + hv -> no + o 97 97 j_n2 = 13 ! n2 + hv -> n + n 98 j_hdo = 14 ! hdo + hv -> products 98 99 99 100 !==== air column increments and rayleigh optical depth … … 167 168 sj(ilay,iw,j_no) = xsno(iw)*yieldno(iw) ! no 168 169 sj(ilay,iw,j_n2) = xsn2(iw)*yieldn2(iw) ! n2 170 ! sj(ilay,iw,j_hdo) = xshdo(iw) ! hdo 169 171 end do 170 172 end do
Note: See TracChangeset
for help on using the changeset viewer.
