Changeset 1430 for trunk/LMDZ.MARS/libf/aeronomars
- Timestamp:
- May 27, 2015, 1:51:37 PM (10 years ago)
- Location:
- trunk/LMDZ.MARS/libf/aeronomars
- Files:
-
- 2 added
- 1 deleted
- 2 edited
-
chimiedata.h (modified) (2 diffs)
-
photochemistry.F (modified) (4 diffs)
-
photolysis.F90 (added)
-
read_phototable.F (deleted)
-
read_phototable.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/chimiedata.h
r655 r1430 1 ! 2 ! --- Data for photolysis lookup table ---3 ! 1 !-------------------------------------------- 2 ! data for photochemistry 3 !-------------------------------------------- 4 4 5 integer nd, nozo, nr, nsza, ntemp, ntau 5 !-------------------------------------------- 6 ! dimensions of photolysis lookup table 7 !-------------------------------------------- 6 8 7 parameter (nd = 29)8 parameter (nozo = 7)9 parameter (nr = nd + 28)10 parameter (nsza = 27)11 parameter (ntemp = 4)12 parameter (ntau = 8)9 integer, parameter :: nd = 13 ! species 10 integer, parameter :: nz = 143 ! altitude 11 integer, parameter :: nozo = 7 ! ozone 12 integer, parameter :: nsza = 27 ! solar zenith angle 13 integer, parameter :: ntemp = 4 ! temperature 14 integer, parameter :: ntau = 8 ! dust 13 15 14 real kb 15 parameter (kb = 1.3806e-23) 16 !-------------------------------------------- 17 18 real, parameter :: kb = 1.3806e-23 16 19 17 20 common/chimiedata/jphot,colairtab,table_ozo 18 21 19 real jphot(ntemp,nsza, 0:200,nozo,ntau,nd)20 real colairtab( 0:200)22 real jphot(ntemp,nsza,nz,nozo,ntau,nd) 23 real colairtab(nz) 21 24 real szatab(nsza) 22 25 real table_ozo(nozo) … … 30 33 31 34 data tautab/0., 0.2, 0.4, 0.6, 0.8, 1., 2., 4./ 35 36 !-------------------------------------------- 37 ! number of reactions in ASIS solver 38 !-------------------------------------------- 39 40 integer, parameter :: nb_phot_max = 18 41 integer, parameter :: nb_reaction_3_max = 6 42 integer, parameter :: nb_reaction_4_max = 30 -
trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F
r1266 r1430 133 133 end do 134 134 135 call phot (nlayer, lswitch, press, temp, sza, tau, dist_sol,136 $ rmco2, rmo3, j)135 call photolysis(nlayer, lswitch, press, temp, sza, tau, dist_sol, 136 $ rmco2, rmo3, j) 137 137 138 138 j_o3_o1d = 5 … … 245 245 integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, 246 246 $ j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2, 247 $ j_ho2, j_no 2, j_ch4_ch3_h, j_ch4_1ch2_h2,247 $ j_ho2, j_no, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2, 248 248 $ j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h, 249 249 $ j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl, … … 295 295 j_o3_o = 6 ! o3 + hv -> o2 + o 296 296 j_h2o = 7 ! h2o + hv -> h + oh 297 j_h do = 8 ! hdo + hv -> d+ oh298 j_h 2o2 = 9 ! h2o2 + hv -> oh + oh299 j_ ho2 = 10 ! ho2 + hv -> oh+ o297 j_h2o2 = 8 ! h2o2 + hv -> oh + oh 298 j_ho2 = 9 ! ho2 + hv -> oh + o 299 j_no = 10 ! no + hv -> n + o 300 300 j_no2 = 11 ! no2 + hv -> no + o 301 j_ch4_ch3_h = 12 ! ch4 + hv -> ch3 + h 302 j_ch4_1ch2_h2 = 13 ! ch4 + hv -> 1ch2 + h2 303 j_ch4_3ch2_h_h = 14 ! ch4 + hv -> 3ch2 + h + h 304 j_ch4_ch_h2_h = 15 ! ch4 + hv -> ch + h2 + h 305 j_ch3o2h = 16 ! ch3o2h + hv -> ch3o + oh 306 j_ch2o_hco = 17 ! ch2o + hv -> h + hco 307 j_ch2o_co = 18 ! ch2o + hv -> h2 + co 308 j_ch3oh = 19 ! ch3oh + hv -> ch3o + h 309 j_c2h6 = 20 ! c2h6 + hv -> products 310 j_hcl = 21 ! hcl + hv -> h + cl 311 j_hocl = 22 ! hocl + hv -> oh + cl 312 j_clo = 23 ! clo + hv -> cl + o 313 j_so2 = 24 ! so2 + hv -> so + o 314 j_so = 25 ! so + hv -> s + o 315 j_h2s = 26 ! h2s + hv -> hs + s 316 j_so3 = 27 ! so2 + hv -> so2 + o 317 j_hno3 = 28 ! hno3 + hv -> oh + no2 318 j_hno4 = 29 ! hno4 + hv -> ho2 + no2 301 j_hno3 = 12 ! hno3 + hv -> oh + no2 302 j_hno4 = 13 ! hno4 + hv -> ho2 + no2 319 303 320 304 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc … … 776 760 c***************************************************************** 777 761 778 subroutine phot(nlayer,779 $ lswitch, press, temp, sza, tauref, dist_sol,780 $ rmco2, rmo3, j)781 782 c*****************************************************************783 784 USE comcstfi_h785 implicit none786 787 #include "chimiedata.h"788 789 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc790 c inputs:791 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc792 793 integer, intent(in) :: nlayer ! number of atmospheric layers794 integer :: lswitch ! interface level between chemistries795 real :: press(nlayer) ! pressure (hPa)796 real :: temp(nlayer) ! temperature (K)797 real :: sza ! solar zenith angle (deg)798 real :: tauref ! optical depth at 7 hpa799 real :: dist_sol ! sun distance (AU)800 real :: rmco2(nlayer) ! co2 mixing ratio801 real :: rmo3(nlayer) ! ozone mixing ratio802 803 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc804 c output:805 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc806 807 real :: j(nlayer,nd) ! interpolated photolysis rates (s-1)808 809 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc810 c local:811 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc812 813 integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,814 $ iozo, isza, itau, it, l815 816 real :: col(nlayer) ! overhead air column (molecule cm-2)817 real :: colo3(nlayer) ! overhead ozone column (molecule cm-2)818 real :: poids(2,2,2,2,2) ! 5D interpolation weights819 real :: tref ! temperature at 1.9 hPa in the gcm (K)820 real :: table_temp(ntemp) ! temperatures at 1.9 hPa in jmars (K)821 real :: cinf, csup, cicol,822 $ ciozo, cisza, citemp, citau823 real :: colo3min, dp, coef824 real :: ratio_o3(nlayer)825 real :: tau826 827 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc828 c day/night criterion829 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc830 831 if (sza .le. 95.) then832 833 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc834 c temperatures at 1.9 hPa in lookup table835 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc836 837 table_temp(1) = 226.2838 table_temp(2) = 206.2839 table_temp(3) = 186.2840 table_temp(4) = 169.8841 842 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc843 c interpolation in solar zenith angle844 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc845 846 indsza = nsza - 1847 do isza = 1,nsza848 if (szatab(isza) .ge. sza) then849 indsza = min(indsza,isza - 1)850 indsza = max(indsza, 1)851 end if852 end do853 cisza = (sza - szatab(indsza))854 $ /(szatab(indsza + 1) - szatab(indsza))855 856 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc857 c interpolation in dust (tau)858 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc859 860 tau = min(tauref, tautab(ntau))861 tau = max(tau, tautab(1))862 863 indtau = ntau - 1864 do itau = 1,ntau865 if (tautab(itau) .ge. tau) then866 indtau = min(indtau,itau - 1)867 indtau = max(indtau, 1)868 end if869 end do870 citau = (tau - tautab(indtau))871 $ /(tautab(indtau + 1) - tautab(indtau))872 873 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc874 c co2 and ozone columns875 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc876 877 c co2 column at model top (molecule.cm-2)878 879 col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.880 $ /(mugaz*g)881 882 c ozone column at model top883 884 colo3(lswitch-1) = 0.885 886 c co2 and ozone columns for other levels (molecule.cm-2)887 888 do l = lswitch-2,1,-1889 dp = (press(l) - press(l+1))*100.890 col(l) = col(l+1)891 $ + (rmco2(l+1) + rmco2(l))*0.5892 $ *6.022e22*dp/(mugaz*g)893 col(l) = min(col(l), colairtab(0))894 colo3(l) = colo3(l+1)895 $ + (rmo3(l+1) + rmo3(l))*0.5896 $ *6.022e22*dp/(mugaz*g)897 end do898 899 c ratio ozone column/minimal theoretical column (0.1 micron-atm)900 c901 c ro3 = 7.171e-10 is the o3 mixing ratio for a uniform902 c profile giving a column 0.1 micron-atmosphere at903 c a surface pressure of 10 hpa.904 905 do l = 1,lswitch-1906 colo3min = col(l)*7.171e-10907 ratio_o3(l) = colo3(l)/colo3min908 ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)909 ratio_o3(l) = max(ratio_o3(l), 1.)910 end do911 912 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc913 c temperature dependence914 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc915 916 c 1) search for temperature at 1.9 hPa (tref): vertical interpolation917 918 tref = temp(1)919 do l = (lswitch-1)-1,1,-1920 if (press(l) .gt. 1.9) then921 cinf = (press(l) - 1.9)922 $ /(press(l) - press(l+1))923 csup = 1. - cinf924 tref = cinf*temp(l+1) + csup*temp(l)925 go to 10926 end if927 end do928 10 continue929 930 c 2) interpolation in temperature931 932 tref = min(tref,table_temp(1))933 tref = max(tref,table_temp(ntemp))934 935 do it = 2, ntemp936 if (table_temp(it) .le. tref) then937 citemp = (log(tref) - log(table_temp(it)))938 $ /(log(table_temp(it-1)) - log(table_temp(it)))939 indtemp = it - 1940 goto 20941 end if942 end do943 20 continue944 945 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc946 c loop over vertical levels947 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc948 949 do l = 1,lswitch-1950 951 c interpolation in air column952 953 do icol = 0,200954 if (colairtab(icol) .lt. col(l)) then955 cicol = (log(col(l)) - log(colairtab(icol)))956 $ /(log(colairtab(icol-1)) - log(colairtab(icol)))957 indcol = icol - 1958 goto 30959 end if960 end do961 30 continue962 963 cc interpolation in ozone column964 965 indozo = nozo - 1966 do iozo = 1,nozo967 if (table_ozo(iozo)*10. .ge. ratio_o3(l)) then968 indozo = min(indozo, iozo - 1)969 indozo = max(indozo, 1)970 end if971 end do972 ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)973 $ /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)974 975 cc 4-dimensional interpolation weights976 977 c poids(temp,sza,co2,o3,tau)978 979 poids(1,1,1,1,1) = citemp980 $ *(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)981 poids(1,1,1,2,1) = citemp982 $ *(1.-cisza)*cicol*ciozo*(1.-citau)983 poids(1,1,2,1,1) = citemp984 $ *(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)985 poids(1,1,2,2,1) = citemp986 $ *(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)987 poids(1,2,1,1,1) = citemp988 $ *cisza*cicol*(1.-ciozo)*(1.-citau)989 poids(1,2,1,2,1) = citemp990 $ *cisza*cicol*ciozo*(1.-citau)991 poids(1,2,2,1,1) = citemp992 $ *cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)993 poids(1,2,2,2,1) = citemp994 $ *cisza*(1.-cicol)*ciozo*(1.-citau)995 poids(2,1,1,1,1) = (1.-citemp)996 $ *(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)997 poids(2,1,1,2,1) = (1.-citemp)998 $ *(1.-cisza)*cicol*ciozo*(1.-citau)999 poids(2,1,2,1,1) = (1.-citemp)1000 $ *(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)1001 poids(2,1,2,2,1) = (1.-citemp)1002 $ *(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)1003 poids(2,2,1,1,1) = (1.-citemp)1004 $ *cisza*cicol*(1.-ciozo)*(1.-citau)1005 poids(2,2,1,2,1) = (1.-citemp)1006 $ *cisza*cicol*ciozo*(1.-citau)1007 poids(2,2,2,1,1) = (1.-citemp)1008 $ *cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)1009 poids(2,2,2,2,1) = (1.-citemp)1010 $ *cisza*(1.-cicol)*ciozo*(1.-citau)1011 1012 poids(1,1,1,1,2) = citemp1013 $ *(1.-cisza)*cicol*(1.-ciozo)*citau1014 poids(1,1,1,2,2) = citemp1015 $ *(1.-cisza)*cicol*ciozo*citau1016 poids(1,1,2,1,2) = citemp1017 $ *(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau1018 poids(1,1,2,2,2) = citemp1019 $ *(1.-cisza)*(1.-cicol)*ciozo*citau1020 poids(1,2,1,1,2) = citemp1021 $ *cisza*cicol*(1.-ciozo)*citau1022 poids(1,2,1,2,2) = citemp1023 $ *cisza*cicol*ciozo*citau1024 poids(1,2,2,1,2) = citemp1025 $ *cisza*(1.-cicol)*(1.-ciozo)*citau1026 poids(1,2,2,2,2) = citemp1027 $ *cisza*(1.-cicol)*ciozo*citau1028 poids(2,1,1,1,2) = (1.-citemp)1029 $ *(1.-cisza)*cicol*(1.-ciozo)*citau1030 poids(2,1,1,2,2) = (1.-citemp)1031 $ *(1.-cisza)*cicol*ciozo*citau1032 poids(2,1,2,1,2) = (1.-citemp)1033 $ *(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau1034 poids(2,1,2,2,2) = (1.-citemp)1035 $ *(1.-cisza)*(1.-cicol)*ciozo*citau1036 poids(2,2,1,1,2) = (1.-citemp)1037 $ *cisza*cicol*(1.-ciozo)*citau1038 poids(2,2,1,2,2) = (1.-citemp)1039 $ *cisza*cicol*ciozo*citau1040 poids(2,2,2,1,2) = (1.-citemp)1041 $ *cisza*(1.-cicol)*(1.-ciozo)*citau1042 poids(2,2,2,2,2) = (1.-citemp)1043 $ *cisza*(1.-cicol)*ciozo*citau1044 1045 cc 4-dimensional interpolation in the lookup table1046 1047 do ij = 1,nd1048 j(l,ij) =1049 $ poids(1,1,1,1,1)1050 $ *jphot(indtemp,indsza,indcol,indozo,indtau,ij)1051 $ + poids(1,1,1,2,1)1052 $ *jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)1053 $ + poids(1,1,2,1,1)1054 $ *jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)1055 $ + poids(1,1,2,2,1)1056 $ *jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)1057 $ + poids(1,2,1,1,1)1058 $ *jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)1059 $ + poids(1,2,1,2,1)1060 $ *jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)1061 $ + poids(1,2,2,1,1)1062 $ *jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)1063 $ + poids(1,2,2,2,1)1064 $ *jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)1065 $ + poids(2,1,1,1,1)1066 $ *jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)1067 $ + poids(2,1,1,2,1)1068 $ *jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)1069 $ + poids(2,1,2,1,1)1070 $ *jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)1071 $ + poids(2,1,2,2,1)1072 $ *jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)1073 $ + poids(2,2,1,1,1)1074 $ *jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)1075 $ + poids(2,2,1,2,1)1076 $ *jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)1077 $ + poids(2,2,2,1,1)1078 $ *jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)1079 $ + poids(2,2,2,2,1)1080 $ *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)1081 1082 $ + poids(1,1,1,1,2)1083 $ *jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)1084 $ + poids(1,1,1,2,2)1085 $ *jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)1086 $ + poids(1,1,2,1,2)1087 $ *jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)1088 $ + poids(1,1,2,2,2)1089 $ *jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)1090 $ + poids(1,2,1,1,2)1091 $ *jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)1092 $ + poids(1,2,1,2,2)1093 $ *jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)1094 $ + poids(1,2,2,1,2)1095 $ *jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)1096 $ + poids(1,2,2,2,2)1097 $ *jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)1098 $ + poids(2,1,1,1,2)1099 $ *jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)1100 $ + poids(2,1,1,2,2)1101 $ *jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)1102 $ + poids(2,1,2,1,2)1103 $ *jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)1104 $ + poids(2,1,2,2,2)1105 $ *jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)1106 $ + poids(2,2,1,1,2)1107 $ *jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)1108 $ + poids(2,2,1,2,2)1109 $ *jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)1110 $ + poids(2,2,2,1,2)1111 $ *jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)1112 $ + poids(2,2,2,2,2)1113 $ *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)1114 end do1115 1116 cc correction for sun distance1117 1118 do ij = 1,nd1119 j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.1120 end do1121 1122 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc1123 c end of loop over vertical levels1124 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc1125 1126 end do1127 1128 else1129 1130 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc1131 c night1132 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc1133 1134 do ij = 1,nd1135 do l = 1,lswitch-11136 j(l,ij) = 0.1137 end do1138 end do1139 1140 end if1141 1142 return1143 end1144 1145 c*****************************************************************1146 1147 762 subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm) 1148 763
Note: See TracChangeset
for help on using the changeset viewer.
