- Timestamp:
- Oct 26, 2018, 8:25:32 PM (6 years ago)
- Location:
- trunk/LMDZ.MARS/libf/aeronomars
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
r2027 r2028 93 93 real :: dt_guess ! first-guess timestep (s) 94 94 real :: dt_corrected ! corrected timestep (s) 95 real :: j(nlayer,nd) ! interpolated photolysis rates (s-1)96 95 real :: time ! internal time (between 0 and ptimestep, in s) 97 96 … … 147 146 if (jonline) then 148 147 tau = tau*press(1)/6.1 ! temporary 149 call photolysis_online(nlayer, lswitch, alt, press, temp,&150 zmmean,i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, &151 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, 152 i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, 148 call photolysis_online(nlayer, alt, press, temp, zmmean, & 149 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 150 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 151 i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm, & 153 152 tau, sza, dist_sol, v_phot) 154 153 else … … 1825 1824 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1826 1825 1827 do l = 1, lswitch-11826 do l = 1,nlayer 1828 1827 rm(l,i_co2) = zycol(l, igcm_co2) 1829 1828 rm(l,i_co) = zycol(l, igcm_co) -
trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
r2027 r2028 1 1 !============================================================================== 2 2 3 subroutine photolysis_online(nlayer, lswitch,alt, press, temp,3 subroutine photolysis_online(nlayer, alt, press, temp, 4 4 $ zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, 5 5 $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, … … 14 14 15 15 integer :: nesp ! total number of chemical species 16 integer :: nlayer , lswitch16 integer :: nlayer 17 17 integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, 18 18 $ i_h2, i_oh, i_ho2, i_h2o2, i_h2o, … … 34 34 ! integer, parameter :: nw = 3789 ! number of spectral intervals (high-res) 35 35 integer, parameter :: nw = 152 ! number of spectral intervals (low-res) 36 real, dimension(nw) :: wl, wc, wu ! lower, center, upper wavelength for each interval 37 save wl, wc, wu 36 real, dimension(nw), save :: wl, wc, wu ! lower, center, upper wavelength for each interval 38 37 39 38 ! solar flux 40 39 41 real, dimension(nw) :: f! solar flux (w.m-2.nm-1) at 1 au40 real, dimension(nw), save :: f ! solar flux (w.m-2.nm-1) at 1 au 42 41 real :: factor 43 save f44 42 45 43 ! cross-sections and quantum yields 46 44 47 integer, parameter :: nabs = 10 ! number of absorbing gases 48 real, dimension(nlayer,nw,nd) :: sj ! general cross-section array (cm2) 49 real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2) 50 real, dimension(nw) :: yieldco2 ! co2 photodissociation yield 51 real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 absorption cross-section at 150-200-250-300 k (cm2) 52 real, dimension(nw) :: yieldo2 ! o2 photodissociation yield 53 real, dimension(nw) :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2) 54 real, dimension(nw) :: xsh2o ! h2o absorption cross-section (cm2) 55 real, dimension(nw) :: xsh2o2 ! h2o2 absorption cross-section (cm2) 56 real, dimension(nw) :: xsho2 ! ho2 absorption cross-section (cm2) 57 real, dimension(nw) :: xsh2 ! h2 absorption cross-section (cm2) 58 real, dimension(nw) :: yieldh2 ! h2 photodissociation yield 59 real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2) 60 real, dimension(nw) :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k 61 real, dimension(nw) :: xsno ! no absorption cross-section (cm2) 62 real, dimension(nw) :: yieldno ! no photodissociation yield 63 real, dimension(nw) :: xsn2 ! n2 absorption cross-section (cm2) 64 real, dimension(nw) :: yieldn2 ! n2 photodissociation yield 65 real, dimension(nw) :: albedo ! surface albedo 45 integer, parameter :: nabs = 10 ! number of absorbing gases 46 integer, parameter :: nphot = 13 ! number of photolysis 47 real, dimension(nlayer,nw,nphot) :: sj ! general cross-section array (cm2) 48 real, dimension(nw), save :: xsco2_195, xsco2_295, xsco2_370 ! co2 absorption cross-section at 195-295-370 k (cm2) 49 real, dimension(nw), save :: yieldco2 ! co2 photodissociation yield 50 real, dimension(nw), save :: xso2_150, xso2_200, ! o2 absorption cross-section at 150-200-250-300 k (cm2) 51 $ xso2_250, xso2_300 52 real, dimension(nw), save :: yieldo2 ! o2 photodissociation yield 53 real, dimension(nw), save :: xso3_218, xso3_298 ! o3 absorption cross-section at 218-298 k (cm2) 54 real, dimension(nw), save :: xsh2o ! h2o absorption cross-section (cm2) 55 real, dimension(nw), save :: xsh2o2 ! h2o2 absorption cross-section (cm2) 56 real, dimension(nw), save :: xsho2 ! ho2 absorption cross-section (cm2) 57 real, dimension(nw), save :: xsh2 ! h2 absorption cross-section (cm2) 58 real, dimension(nw), save :: yieldh2 ! h2 photodissociation yield 59 real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294 ! no2 absorption cross-section at 220-294 k (cm2) 60 real, dimension(nw), save :: yldno2_248, yldno2_298 ! no2 quantum yield at 248-298 k 61 real, dimension(nw), save :: xsno ! no absorption cross-section (cm2) 62 real, dimension(nw), save :: yieldno ! no photodissociation yield 63 real, dimension(nw), save :: xsn2 ! n2 absorption cross-section (cm2) 64 real, dimension(nw), save :: yieldn2 ! n2 photodissociation yield 65 real, dimension(nw), save :: albedo ! surface albedo 66 66 67 67 ! atmosphere … … 89 89 integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, 90 90 $ a_no2, a_n2 91 integer :: i, ilay, iw, ialt, mopt 91 integer :: i, ilay, iw, ialt 92 integer, save :: mopt 92 93 real :: deltaj 93 94 94 logical, save :: firstcall 95 data firstcall / .true. / 96 97 ! high-res/low-res switch 98 ! 99 ! mopt = 1 high-resolution 100 ! mopt = 2 low-resolution (recommended for gcm use) 101 102 if (nw >= 3789) then 103 mopt = 1 104 else 105 mopt = 2 106 end if 95 logical, save :: firstcall = .true. 107 96 108 97 ! absorbing gas numbering … … 119 108 a_n2 = 10 ! no2 120 109 121 ! photodissociation rates numbering 110 ! photodissociation rates numbering. 111 ! photodissociations must be ordered the same way in subroutine "indice" 122 112 123 113 j_o2_o = 1 ! o2 + hv -> o + o … … 143 133 if (firstcall) then 144 134 135 ! high-res/low-res switch 136 ! 137 ! mopt = 1 high-resolution 138 ! mopt = 2 low-resolution (recommended for gcm use) 139 140 mopt = 2 141 145 142 ! set wavelength grid 146 143 … … 208 205 ! o2: 209 206 210 call seto2(n d, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,207 call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200, 211 208 $ xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d, 212 209 $ colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj) … … 214 211 ! co2: 215 212 216 call setco2(n d, nlayer, nw, wc, temp, xsco2_195, xsco2_295,213 call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295, 217 214 $ xsco2_370, yieldco2, j_co2_o, j_co2_o1d, 218 215 $ colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj) … … 220 217 ! o3: 221 218 222 call seto3(n d, nlayer, nw, wc, temp, xso3_218, xso3_298,219 call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298, 223 220 $ j_o3_o, j_o3_o1d, 224 221 $ colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj) … … 226 223 ! h2o2: 227 224 228 call seth2o2(n d, nlayer, nw, wc, temp, xsh2o2, j_h2o2,225 call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2, 229 226 $ colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj) 230 227 231 228 ! no2: 232 229 233 call setno2(n d, nlayer, nw, wc, temp, xsno2, xsno2_220,230 call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220, 234 231 $ xsno2_294, yldno2_248, yldno2_298, j_no2, 235 232 $ colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj) … … 301 298 !==== initialise photolysis rates 302 299 303 v_phot(:,1:n d) = 0.300 v_phot(:,1:nphot) = 0. 304 301 305 302 !==== start of wavelength lopp … … 320 317 ! spherical actinic flux 321 318 322 do ilay = 1, lswitch-1319 do ilay = 1,nlayer 323 320 saflux(ilay) = f(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay)) 324 321 end do … … 326 323 ! photolysis rate integration 327 324 328 do i = 1,n d329 do ilay = 1, lswitch-1325 do i = 1,nphot 326 do ilay = 1,nlayer 330 327 deltaj = saflux(ilay)*sj(ilay,iw,i) 331 328 v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw)) … … 335 332 ! eliminate small values 336 333 337 where (v_phot(:,1:n d) < 1.e-30)338 v_phot(:,1:n d) = 0.334 where (v_phot(:,1:nphot) < 1.e-30) 335 v_phot(:,1:nphot) = 0. 339 336 end where 340 337 … … 343 340 else ! night 344 341 345 v_phot(:,1:n d) = 0.342 v_phot(:,1:nphot) = 0. 346 343 347 344 end if ! sza … … 945 942 ! local 946 943 947 integer, parameter :: kdata = 36000944 integer, parameter :: kdata = 42000 948 945 real, parameter :: deltax = 1.e-4 949 946 real, dimension(kdata) :: x1, y1, y2, y3, xion, ion … … 964 961 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 965 962 c 966 c 195K: Chan et al. (1993) + Starck et al. (2006)963 c 195K: huestis and berkowitz (2010) + Starck et al. (2006) 967 964 c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation 968 c 11 lignes d'en-tete et n1 = 34537969 c fichier: co2_euv_uv_195k.txt970 965 c 971 c 295K: Chan et al. (1993) + Starck et al. (2006)966 c 295K: huestis and berkowitz (2010) + Starck et al. (2006) 972 967 c + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation 973 c 11 lignes d'en-tete et n2 = 35360974 c fichier: co2_euv_uv_295k.txt975 968 c 976 c 370K: Chan et al. (1993) + Starck et al. (2006)969 c 370K: huestis and berkowitz (2010) + Starck et al. (2006) 977 970 c + Lewis and Carver (1983) + extrapolation 978 c 11 lignes d'en-tete et n3 = 3884979 c fichier: co2_euv_uv_370k.txt980 971 c 981 972 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 982 c 983 n1 = 34537984 n2 = 35360985 n3 = 3884986 c 987 c195K:988 c 989 fil = trim(datadir)//'/cross_sections/co2_euv_uv_ 195k.txt'973 974 n1 = 40769 975 n2 = 41586 976 n3 = 10110 977 978 ! 195K: 979 980 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_195k.txt' 990 981 print*, 'section efficace CO2 195K: ', fil 991 982 … … 1013 1004 xsco2_195(l) = yg(l) 1014 1005 END DO 1015 c 1016 c295K:1017 c 1018 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2 95k.txt'1006 1007 ! 295K: 1008 1009 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_295k.txt' 1019 1010 print*, 'section efficace CO2 295K: ', fil 1020 c 1011 1021 1012 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1022 1013 DO i = 1,11 1023 1014 read(kin,*) 1024 1015 END DO 1025 c 1016 1026 1017 DO i = 1, n2 1027 1018 READ(kin,*) x1(i), y1(i) … … 1038 1029 STOP 1039 1030 ENDIF 1040 c 1031 1041 1032 DO l = 1, nw-1 1042 1033 xsco2_295(l) = yg(l) 1043 1034 END DO 1044 c 1045 c370K:1046 c 1047 fil = trim(datadir)//'/cross_sections/co2_euv_uv_ 370k.txt'1035 1036 ! 370K: 1037 1038 fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_370k.txt' 1048 1039 print*, 'section efficace CO2 370K: ', fil 1049 c 1040 1050 1041 OPEN(UNIT=kin,FILE=fil,STATUS='old') 1051 1042 DO i = 1,11 1052 1043 read(kin,*) 1053 1044 END DO 1054 c 1045 1055 1046 DO i = 1, n3 1056 1047 READ(kin,*) x1(i), y1(i) 1057 1048 END DO 1058 1049 CLOSE (kin) 1059 c 1050 1060 1051 CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.) 1061 1052 CALL addpnt(x1,y1,kdata,n3, 0.,0.)
Note: See TracChangeset
for help on using the changeset viewer.