Changeset 2836 for trunk/LMDZ.VENUS/libf
- Timestamp:
- Nov 30, 2022, 4:46:57 PM (2 years ago)
- Location:
- trunk/LMDZ.VENUS/libf/phyvenus
- Files:
-
- 1 added
- 15 edited
-
chemparam_mod.F90 (modified) (5 diffs)
-
clesphys.h (modified) (2 diffs)
-
concentrations2.F (modified) (3 diffs)
-
conf_phys.F90 (modified) (2 diffs)
-
euvheat.F90 (modified) (6 diffs)
-
hrtherm.F (modified) (2 diffs)
-
iono_h.F90 (added)
-
jthermcalc_e107.F (modified) (20 diffs)
-
moldiff_red.F90 (modified) (1 diff)
-
nonoro_gwd_ran_mod.F90 (modified) (3 diffs)
-
param_read_e107.F (modified) (6 diffs)
-
photochemistry_venus.F90 (modified) (34 diffs)
-
photolysis_mod.F90 (modified) (2 diffs)
-
photolysis_online.F (modified) (2 diffs)
-
physiq_mod.F (modified) (13 diffs)
-
phytrac_chimie.F (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90
r2795 r2836 13 13 i_cocl2, i_s, i_so, i_so2, i_so3, & 14 14 i_s2o2, i_ocs, i_hso3, i_h2so4, i_s2, & 15 i_clso2, i_oscl, i_n2, i_he, i_no2, & 16 i_no, i_n, i_n2d 15 i_clso2, i_oscl, i_n2, i_he, i_n, i_no, & 16 i_no2, i_n2d, & 17 i_co2plus, i_coplus, i_oplus, i_o2plus, & 18 i_n2plus, i_hplus, i_h2oplus, i_nplus, & 19 i_ohplus, i_cplus, i_noplus, i_h3oplus, & 20 i_hcoplus, i_hco2plus, i_elec 17 21 18 22 INTEGER, SAVE :: i_h2oliq, i_h2so4liq … … 26 30 INTEGER, SAVE :: nmicro ! number of microphysical tracers 27 31 28 REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr 32 REAL, DIMENSION(:), SAVE, ALLOCATABLE :: M_tr ! Molecular Mass of tracers 33 REAL, DIMENSION(:), SAVE, ALLOCATABLE :: type_tr ! type of tracer 29 34 30 35 REAL, DIMENSION(:,:),SAVE, ALLOCATABLE :: no_emission … … 697 702 INTEGER :: i 698 703 699 ALLOCATE(M_tr(nqtot)) 700 704 ALLOCATE(M_tr(nqtot)) ! Molecular Mass of Chemistry Tracers 705 ALLOCATE(type_tr(nqtot)) ! Type of chemistry tracers 1: neutral, 2: ion, 3: liquide 706 707 ! Initialisation index chemistry tracers 708 709 ! Neutral gas Tracer 710 i_co2 = 0 711 i_co = 0 712 i_h2 = 0 713 i_h2o = 0 714 i_o1d = 0 715 i_o = 0 716 i_o2 = 0 717 i_o2dg = 0 718 i_o3 = 0 719 i_h = 0 720 i_oh = 0 721 i_ho2 = 0 722 i_h2o2 = 0 723 i_cl = 0 724 i_clo = 0 725 i_cl2 = 0 726 i_hcl = 0 727 i_hocl = 0 728 i_clco = 0 729 i_clco3 = 0 730 i_cocl2 = 0 731 i_s = 0 732 i_so = 0 733 i_so2 = 0 734 i_so3 = 0 735 i_s2o2 = 0 736 i_ocs = 0 737 i_hso3 = 0 738 i_h2so4 = 0 739 i_s2 = 0 740 i_clso2 = 0 741 i_oscl = 0 742 i_n2 = 0 743 i_he = 0 744 i_n2d = 0 745 i_n = 0 746 i_no = 0 747 i_no2 = 0 748 ! ION TRACERS FOR IONCHEM = TRUE 749 i_co2plus = 0 750 i_coplus = 0 751 i_oplus = 0 752 i_o2plus = 0 753 i_n2plus = 0 754 i_hplus = 0 755 i_h2oplus = 0 756 i_nplus = 0 757 i_ohplus = 0 758 i_cplus = 0 759 i_noplus = 0 760 i_h3oplus = 0 761 i_hcoplus = 0 762 i_hco2plus = 0 763 i_elec = 0 764 ! MICROPHYSICAL TRACERS FOR CL_SCHEME=1 765 i_h2oliq = 0 766 i_h2so4liq = 0 767 701 768 DO i=1, nqtot 702 769 … … 705 772 706 773 SELECT CASE(tname(i)) 774 ! NEUTRAL TRACERS 707 775 CASE('co2') 708 776 i_co2=i 709 777 PRINT*,'co2',i_co2 710 M_tr(i_co2) = 44.0095 778 M_tr(i_co2) = 44.0095 779 type_tr(i_co2) = 1 711 780 CASE('co') 712 781 i_co=i 713 782 PRINT*,'co',i_co 714 M_tr(i_co)=28.0101 783 M_tr(i_co)=28.0101 784 type_tr(i_co) = 1 715 785 CASE('h2') 716 786 i_h2=i 717 787 PRINT*,'h2',i_h2 718 788 M_tr(i_h2)= 2.01588 789 type_tr(i_h2) = 1 719 790 CASE('h2o') 720 791 i_h2o=i 721 792 PRINT*,'h2o',i_h2o 722 793 M_tr(i_h2o)=18.0153 794 type_tr(i_h2o) = 1 723 795 CASE('o1d') 724 796 i_o1d=i 725 797 PRINT*,'o1d',i_o1d 726 M_tr(i_o1d)=15.994 798 M_tr(i_o1d)=15.994 799 type_tr(i_o1d) = 1 727 800 CASE('o') 728 801 i_o=i 729 802 PRINT*,'o',i_o 730 M_tr(i_o)=15.994 803 M_tr(i_o)=15.994 804 type_tr(i_o) = 1 731 805 CASE('o2') 732 806 i_o2=i 733 807 PRINT*,'o2',i_o2 734 808 M_tr(i_o2)=31.9988 809 type_tr(i_o2) = 1 735 810 CASE('o2dg') 736 811 i_o2dg=i 737 812 PRINT*,'o2dg',i_o2dg 738 813 M_tr(i_o2dg)=31.9988 814 type_tr(i_o2dg) = 1 739 815 CASE('o3') 740 816 i_o3=i 741 817 PRINT*,'o3',i_o3 742 M_tr(i_o3)= 47.9982 818 M_tr(i_o3)= 47.9982 819 type_tr(i_o3) = 1 743 820 CASE('h') 744 821 i_h=i 745 822 PRINT*,'h',i_h 746 M_tr(i_h)= 1.00794 823 M_tr(i_h)= 1.00794 824 type_tr(i_h) = 1 747 825 CASE('oh') 748 826 i_oh=i 749 827 PRINT*,'oh',i_oh 750 828 M_tr(i_oh)=17.0073 829 type_tr(i_oh) = 1 751 830 CASE('ho2') 752 831 i_ho2=i 753 832 PRINT*,'ho2',i_ho2 754 833 M_tr(i_ho2)=33.0067 834 type_tr(i_ho2) = 1 755 835 CASE('h2o2') 756 836 i_h2o2=i 757 837 PRINT*,'h2o2',i_h2o2 758 M_tr(i_h2o2)=34.0147 838 M_tr(i_h2o2)=34.0147 839 type_tr(i_h2o2) = 1 759 840 CASE('cl') 760 841 i_cl=i 761 842 PRINT*,'cl',i_cl 762 843 M_tr(i_cl)=35.453 844 type_tr(i_cl) = 1 763 845 CASE('clo') 764 846 i_clo=i 765 847 PRINT*,'clo',i_clo 766 848 M_tr(i_clo)=51.452 849 type_tr(i_clo) = 1 767 850 CASE('cl2') 768 851 i_cl2=i 769 852 PRINT*,'cl2',i_cl2 770 853 M_tr(i_cl2)=70.906 854 type_tr(i_cl2) = 1 771 855 CASE('hcl') 772 856 i_hcl=i 773 857 PRINT*,'hcl',i_hcl 774 858 M_tr(i_hcl)=36.461 859 type_tr(i_hcl) = 1 775 860 CASE('hocl') 776 861 i_hocl=i 777 862 PRINT*,'hocl',i_hocl 778 863 M_tr(i_hocl)=52.46 779 CASE('clco') 864 type_tr(i_hocl) = 1 865 CASE('clco') 780 866 i_clco=i 781 867 PRINT*,'clco',i_clco 782 868 M_tr(i_clco)=63.463 783 CASE('clco3') 869 type_tr(i_clco) = 1 870 CASE('clco3') 784 871 i_clco3=i 785 872 PRINT*,'clco3',i_clco3 786 873 M_tr(i_clco3)=95.462 787 CASE('cocl2') 874 type_tr(i_clco3) = 1 875 CASE('cocl2') 788 876 i_cocl2=i 789 877 PRINT*,'cocl2',i_cocl2 790 878 M_tr(i_cocl2)=98.916 791 CASE('s') 879 type_tr(i_cocl2) = 1 880 CASE('s') 792 881 i_s=i 793 882 PRINT*,'s',i_s 794 883 M_tr(i_s)=32.065 795 CASE('so') 884 type_tr(i_s) = 1 885 CASE('so') 796 886 i_so=i 797 887 PRINT*,'so',i_so 798 888 M_tr(i_so)=48.0644 799 CASE('so2') 889 type_tr(i_so) = 1 890 CASE('so2') 800 891 i_so2=i 801 892 PRINT*,'so2',i_so2 802 893 M_tr(i_so2)=64.064 803 CASE('so3') 894 type_tr(i_so2) = 1 895 CASE('so3') 804 896 i_so3=i 805 897 PRINT*,'so3',i_so3 806 898 M_tr(i_so3)=80.063 807 CASE('s2o2') 899 type_tr(i_so3) = 1 900 CASE('s2o2') 808 901 i_s2o2=i 809 902 PRINT*,'s2o2',i_s2o2 810 903 M_tr(i_s2o2)= 96.1288 811 CASE('ocs') 904 type_tr(i_s2o2) = 1 905 CASE('ocs') 812 906 i_ocs=i 813 907 PRINT*,'ocs',i_ocs 814 908 M_tr(i_ocs)=60.0751 815 CASE('hso3') 909 type_tr(i_ocs) = 1 910 CASE('hso3') 816 911 i_hso3=i 817 912 PRINT*,'hso3',i_hso3 818 913 M_tr(i_hso3)=81.071 819 CASE('h2so4') 914 type_tr(i_hso3) = 1 915 CASE('h2so4') 820 916 i_h2so4=i 821 917 PRINT*,'h2so4',i_h2so4 822 918 M_tr(i_h2so4)=98.078 823 CASE('s2') 919 type_tr(i_h2so4) = 1 920 CASE('s2') 824 921 i_s2=i 825 922 PRINT*,'s2',i_s2 826 923 M_tr(i_s2)=64.13 827 CASE('clso2') 924 type_tr(i_s2) = 1 925 CASE('clso2') 828 926 i_clso2=i 829 927 PRINT*,'clso2',i_clso2 830 928 M_tr(i_clso2)=99.517 831 CASE('oscl') 929 type_tr(i_clso2) = 1 930 CASE('oscl') 832 931 i_oscl=i 833 932 PRINT*,'oscl',i_oscl 834 933 M_tr(i_oscl)=83.517 835 CASE('n2') 836 i_n2=i 837 M_tr(i_n2)=28.013 838 CASE('he') 839 i_he=i 840 M_tr(i_he)=4.0026 841 CASE('no2') 842 i_no2=i 843 M_tr(i_no2)=46.0055 844 CASE('no') 845 i_no = i 846 M_tr(i_no)=30.0061 847 CASE('n') 848 i_n = i 849 M_tr(i_n)=14.0067 850 CASE('n2d') 851 i_n2d = i 852 M_tr(i_n2d)=14.0067 934 type_tr(i_oscl) = 1 935 CASE('n2') 936 i_n2=i 937 PRINT*,'n2',i_n2 938 M_tr(i_n2)=28.013 939 type_tr(i_n2) = 1 940 CASE('he') 941 i_he=i 942 PRINT*,'he',i_he 943 M_tr(i_he)=4.0026 944 type_tr(i_he) = 1 945 CASE('n2d') 946 i_n2d=i 947 PRINT*,'n2d',i_n2d 948 M_tr(i_n2d)=14.0067 949 type_tr(i_n2d) = 1 950 CASE('n') 951 i_n=i 952 PRINT*,'n',i_n 953 M_tr(i_n)=14.0067 954 type_tr(i_n) = 1 955 CASE('no') 956 i_no=i 957 PRINT*,'no',i_no 958 M_tr(i_no)=30.0061 959 type_tr(i_no) = 1 960 CASE('no2') 961 i_no2=i 962 PRINT*,'no2',i_no2 963 M_tr(i_no2)=46.0055 964 type_tr(i_no2) = 1 965 ! ION TRACERS FOR OK_IONCHEM = TRUE 966 CASE('co2plus') 967 i_co2plus=i 968 PRINT*,'co2plus',i_co2plus 969 M_tr(i_co2plus)=44.0095 970 type_tr(i_co2plus) = 2 971 CASE('coplus') 972 i_coplus=i 973 PRINT*,'coplus',i_coplus 974 M_tr(i_coplus)=28.0101 975 type_tr(i_coplus) = 2 976 CASE('oplus') 977 i_oplus=i 978 PRINT*,'oplus',i_oplus 979 M_tr(i_oplus)=15.994 980 type_tr(i_oplus) = 2 981 CASE('o2plus') 982 i_o2plus=i 983 PRINT*,'o2plus',i_o2plus 984 M_tr(i_o2plus)=31.9988 985 type_tr(i_o2plus) = 2 986 CASE('n2plus') 987 i_n2plus=i 988 PRINT*,'n2plus',i_n2plus 989 M_tr(i_n2plus)=28.013 990 type_tr(i_n2plus) = 2 991 CASE('hplus') 992 i_hplus=i 993 PRINT*,'hplus',i_hplus 994 M_tr(i_hplus)=1.00794 995 type_tr(i_hplus) = 2 996 CASE('h2oplus') 997 i_h2oplus=i 998 PRINT*,'h2oplus',i_h2oplus 999 M_tr(i_h2oplus)=18.0153 1000 type_tr(i_h2oplus) = 2 1001 CASE('nplus') 1002 i_nplus=i 1003 PRINT*,'nplus',i_nplus 1004 M_tr(i_nplus)=14.0067 1005 type_tr(i_nplus) = 2 1006 CASE('ohplus') 1007 i_ohplus=i 1008 PRINT*,'ohplus',i_ohplus 1009 M_tr(i_ohplus)=17.0073 1010 type_tr(i_ohplus) = 2 1011 CASE('cplus') 1012 i_cplus=i 1013 PRINT*,'cplus',i_cplus 1014 M_tr(i_cplus)=12.011 1015 type_tr(i_cplus) = 2 1016 CASE('noplus') 1017 i_noplus=i 1018 PRINT*,'noplus',i_noplus 1019 M_tr(i_noplus)=30.0061 1020 type_tr(i_noplus) = 2 1021 CASE('h3oplus') 1022 i_h3oplus=i 1023 PRINT*,'h3oplus',i_h3oplus 1024 M_tr(i_h3oplus)=19.0232 1025 type_tr(i_h3oplus) = 2 1026 CASE('hcoplus') 1027 i_hcoplus=i 1028 PRINT*,'hcoplus',i_hcoplus 1029 M_tr(i_hcoplus)=29.0180 1030 type_tr(i_hcoplus) = 2 1031 CASE('hco2plus') 1032 i_hco2plus=i 1033 PRINT*,'hco2plus',i_hco2plus 1034 M_tr(i_hco2plus)=45. 1035 type_tr(i_hco2plus) = 2 1036 CASE('elec') 1037 i_elec=i 1038 PRINT*,'elec',i_elec 1039 M_tr(i_elec)=1./1822.89 1040 type_tr(i_elec) = 2 853 1041 ! MICROPHYSICAL TRACERS FOR CL_SCHEME=1 854 1042 CASE('h2oliq') … … 856 1044 PRINT*,'h2oliq',i_h2oliq 857 1045 M_tr(i_h2oliq)=18.0153 1046 type_tr(i_h2oliq) = 3 858 1047 CASE('h2so4liq') 859 1048 i_h2so4liq=i 860 1049 PRINT*,'h2so4liq',i_h2so4liq 861 1050 M_tr(i_h2so4liq)=98.078 1051 type_tr(i_h2so4liq) = 3 862 1052 ! MICROPHYSICAL TRACERS FOR CL_SCHEME=2 863 1053 CASE('M0_aer') 864 1054 i_m0_aer=i 1055 type_tr(i_m0_aer) = 10 865 1056 PRINT*,'M0_aer',i_m0_aer 866 1057 CASE('M3_aer') 867 1058 i_m3_aer=i 1059 type_tr(i_m3_aer) = 10 868 1060 PRINT*,'M3_aer',i_m3_aer 869 1061 CASE('M0_m1drop') 870 1062 i_m0_mode1drop=i 1063 type_tr(i_m0_mode1drop) = 10 871 1064 PRINT*,'M0_m1drop',i_m0_mode1drop 872 1065 CASE('M0_m1ccn') 873 1066 i_m0_mode1ccn=i 1067 type_tr(i_m0_mode1ccn) = 10 874 1068 PRINT*,'M0_m1ccn',i_m0_mode1ccn 875 1069 CASE('M3_m1sa') 876 1070 i_m3_mode1sa=i 1071 type_tr(i_m3_mode1sa) = 10 877 1072 PRINT*,'M3_m1sa',i_m3_mode1sa 878 1073 CASE('M3_m1w') 879 1074 i_m3_mode1w=i 1075 type_tr(i_m3_mode1w) = 10 880 1076 PRINT*,'M3_m1w',i_m3_mode1w 881 1077 CASE('M3_m1ccn') 882 1078 i_m3_mode1ccn=i 1079 type_tr(i_m3_mode1ccn) = 10 883 1080 PRINT*,'M3_m1ccn',i_m3_mode1ccn 884 1081 CASE('M0_m2drop') 885 1082 i_m0_mode2drop=i 1083 type_tr(i_m0_mode2drop) = 10 886 1084 PRINT*,'M0_m2drop',i_m0_mode2drop 887 1085 CASE('M0_m2ccn') 888 1086 i_m0_mode2ccn=i 1087 type_tr(i_m0_mode2ccn) = 10 889 1088 PRINT*,'M0_m2ccn',i_m0_mode2ccn 890 1089 CASE('M3_m2sa') 891 1090 i_m3_mode2sa=i 1091 type_tr(i_m3_mode2sa) = 10 892 1092 PRINT*,'M3_m2sa',i_m3_mode2sa 893 1093 CASE('M3_m2w') 894 1094 i_m3_mode2w=i 1095 type_tr(i_m3_mode2w) = 10 895 1096 PRINT*,'M3_m2w',i_m3_mode2w 896 1097 CASE('M3_m2ccn') 897 1098 i_m3_mode2ccn=i 1099 type_tr(i_m3_mode2ccn) = 10 898 1100 PRINT*,'M3_m2ccn',i_m3_mode2ccn 899 1101 END SELECT -
trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h
r2580 r2836 10 10 LOGICAL cycle_diurne,soil_model 11 11 LOGICAL ok_orodr,ok_orolf,ok_gw_nonoro 12 LOGICAL ok_kzmin,tuneupperatm 12 LOGICAL ok_kzmin,tuneupperatm, ok_ionchem, ok_jonline 13 13 LOGICAL callnlte,callnirco2,callthermos 14 14 LOGICAL ok_cloud, ok_chem, reinit_trac, ok_sedim … … 29 29 & tuneupperatm,callnlte,callnirco2,callthermos, & 30 30 & ok_cloud, ok_chem, reinit_trac, ok_sedim, & 31 & ok_clmain, physideal, startphy_file 31 & ok_clmain, physideal, startphy_file, ok_ionchem, ok_jonline 32 32 33 33 COMMON/clesphys_i/ nbapp_rad, nbapp_chim, & -
trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F
r2686 r2836 65 65 nbq = 0 ! to count number of tracers used in this subroutine 66 66 67 ! aki values comes from Aeronomy part B by P.M. BANKS / G. KOCKARTS – page 12-20-24) 68 69 ! Ions are not here because the sum of all ions abundace is not 70 ! above 10^-4 until 250 km and we don't have their cpi and aki. 71 72 ! Heat capacity for H, He, N, N2, O, O2, CO2, CO: 73 ! Circular of the bureau of Standards no. 564: tables of thermal properties of gases comprising [...] 74 ! Tables of Thermodynamic and Transport Properties of Air, Argon, Carbon Dioxide, Carbon Monoxide, 75 ! Hydrogen (molecular and atomic), Nitrogen (molecular and atomic), Oxygen (molecular and atomic), and Steam 76 ! https://www.govinfo.gov/content/pkg/GOVPUB-C13-89baf9f9b4a43e5f25820bd51b0f3f11/pdf/GOVPUB-C13-89baf9f9b4a43e5f25820bd51b0f3f11.pdf 77 78 67 79 if (i_co2 /= 0) then 68 80 nbq = nbq + 1 … … 144 156 cpi(nbq) = 1.870e3 145 157 end if 146 c if (i_n /= 0) then 147 c nbq = nbq + 1 148 c niq(nbq) = i_n 149 c aki(nbq) = 0.0 150 c cpi(nbq) = 0.0 151 c endif 152 c if(i_no /= 0) then 153 c nbq = nbq + 1 154 c niq(nbq) = i_no 155 c aki(nbq) = 0.0 156 c cpi(nbq) = 0.0 157 c endif 158 c if(i_no2 /= 0) then 159 c nbq = nbq + 1 160 c niq(nbq) = i_no2 161 c aki(nbq) = 0.0 162 c cpi(nbq) = 0.0 163 c endif 164 c if(i_n2d /= 0) then 165 c nbq = nbq + 1 166 c niq(nbq) = i_n2d 167 c aki(nbq) = 0.0 168 c cpi(nbq) = 0.0 169 c endif 158 if (i_he /= 0) then 159 nbq = nbq + 1 160 niq(nbq) = i_he 161 aki(nbq) = 29.9e-4 162 cpi(nbq) = 5.2e3 163 end if 164 if (i_n /= 0) then 165 nbq = nbq + 1 166 niq(nbq) = i_n 167 aki(nbq) = 0.0 168 cpi(nbq) = 1.4844e3 169 endif 170 if(i_no /= 0) then 171 nbq = nbq + 1 172 niq(nbq) = i_no 173 aki(nbq) = 0.0 174 cpi(nbq) = 0.0 175 endif 176 if(i_no2 /= 0) then 177 nbq = nbq + 1 178 niq(nbq) = i_no2 179 aki(nbq) = 0.0 180 cpi(nbq) = 0.0 181 endif 182 if(i_n2d /= 0) then 183 nbq = nbq + 1 184 niq(nbq) = i_n2d 185 aki(nbq) = 0.0 186 cpi(nbq) = 1.4844e3 !? 187 endif 170 188 if (i_he /= 0) then 171 189 nbq = nbq + 1 … … 255 273 end do 256 274 257 258 275 cpnew(ig,l) = cpnew(ig,l)/ntot 259 276 akknew(ig,l)= akknew(ig,l)/ntot -
trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90
r2580 r2836 485 485 tuneupperatm = .false. 486 486 call getin('tuneupperatm',tuneupperatm) 487 488 ! 489 !Config Key = ok_jonline 490 !Config Desc = 491 !Config Def = .false. 492 !Config Help = 493 ! 494 ok_jonline = .false. 495 call getin('ok_jonline',ok_jonline) 496 497 ! 498 !Config Key = ok_ionchem 499 !Config Desc = 500 !Config Def = .false. 501 !Config Help = 502 ! 503 ok_ionchem = .false. 504 call getin('ok_ionchem',ok_ionchem) 505 506 if ((ok_jonline.eq..false.).and.(.true..eq.ok_ionchem)) then 507 write(*,*) "Attention, incoherence :" 508 write(*,*) "ok_jonline=",ok_jonline," / ok_ionchem=",ok_ionchem 509 write(*,*) "Si vous souhaitez les ions, ok_jonline==.true." 510 write(*,*) "Si vous ne voulez pas des ions, ok_ionchem==.false." 511 write(*,*) "Verifiez votre physiq.def" 512 stop 513 endif 487 514 488 515 ! … … 547 574 write(lunout,*)' euveff = ',euveff 548 575 write(lunout,*)' tuneupperatm = ',tuneupperatm 549 576 write(lunout,*)' ok_jonline = ',ok_jonline 577 write(lunout,*)' ok_ionchem = ',ok_ionchem 578 550 579 end subroutine conf_phys 551 580 -
trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90
r2464 r2836 76 76 !!! If the values are changed there, the same has to be done here !!! 77 77 78 integer,parameter :: ix_co2=1 79 integer,parameter :: ix_o=3 80 integer,parameter :: ix_co=4 81 integer,parameter :: ix_n2=13 82 78 ! integer,parameter :: ix_co2=1 79 ! integer,parameter :: ix_o=3 80 ! integer,parameter :: ix_co=4 81 ! integer,parameter :: ix_n2=13 82 83 integer,parameter :: ix_co2 = 1 84 integer,parameter :: ix_co = 2 85 integer,parameter :: ix_o = 3 86 integer,parameter :: ix_o1d = 4 87 integer,parameter :: ix_o2 = 5 88 integer,parameter :: ix_o3 = 6 89 integer,parameter :: ix_h = 7 90 integer,parameter :: ix_h2 = 8 91 integer,parameter :: ix_oh = 9 92 integer,parameter :: ix_ho2 = 10 93 integer,parameter :: ix_h2o2 = 11 94 integer,parameter :: ix_h2o = 12 95 integer,parameter :: ix_n = 13 96 integer,parameter :: ix_n2d = 14 97 integer,parameter :: ix_no = 15 98 integer,parameter :: ix_no2 = 16 99 integer,parameter :: ix_n2 = 17 83 100 84 101 ! Tracer indexes in the GCM: 85 102 integer,save :: g_co2=0 86 103 integer,save :: g_o=0 104 integer,save :: g_o2=0 105 integer,save :: g_h2=0 106 integer,save :: g_h2o2=0 107 integer,save :: g_h2o=0 108 integer,save :: g_o3=0 109 integer,save :: g_n2=0 110 integer,save :: g_n=0 111 integer,save :: g_no=0 87 112 integer,save :: g_co=0 88 integer,save :: g_n2=0 113 integer,save :: g_h=0 114 integer,save :: g_no2=0 115 integer,save :: g_oh=0 116 integer,save :: g_ho2=0 117 integer,save :: g_o1d=0 118 integer,save :: g_n2d=0 89 119 90 120 logical,save :: firstcall=.true. … … 110 140 stop 111 141 else 112 nspeuv_vgcm=nspeuv_vgcm+1142 nspeuv_vgcm=nspeuv_vgcm+1 113 143 endif 114 144 g_co=i_co … … 122 152 ! n2 123 153 g_n2=i_n2 124 if (g_n2.eq.0) then 125 write(*,*) "euvheat: Error; no N2 tracer !!!" 126 ! write(*,*) "N2 needed if NO is in traceur.def" 127 stop 128 else 129 nspeuv_vgcm=nspeuv_vgcm+1 130 endif 131 132 ! g_o2=igcm_o2 133 ! if (g_o2.eq.0) then 134 ! write(*,*) "euvheat: Error; no O2 tracer !!!" 154 if (g_n2.eq.0) then 155 write(*,*) "euvheat: Error; no N2 tracer !!!" 156 ! write(*,*) "N2 needed if NO is in traceur.def" 157 stop 158 else 159 nspeuv_vgcm=nspeuv_vgcm+1 160 endif 161 g_o2=i_o2 162 if (g_o2.eq.0) then 163 write(*,*) "euvheat: Error; no O2 tracer !!!" 135 164 ! write(*,*) "O2 is always needed if calleuv=.true." 136 !stop137 !else138 ! nespeuv=nespeuv+1139 !endif140 ! g_h2=igcm_h2141 !if (g_h2.eq.0) then142 !write(*,*) "euvheat: Error; no H2 tracer !!!"165 stop 166 else 167 nspeuv_vgcm=nspeuv_vgcm+1 168 endif 169 g_h2=i_h2 170 if (g_h2.eq.0) then 171 write(*,*) "euvheat: Error; no H2 tracer !!!" 143 172 ! write(*,*) "H2 is always needed if calleuv=.true." 144 !stop145 !else146 ! nespeuv=nespeuv+1147 !endif148 ! g_oh=igcm_oh149 !if (g_oh.eq.0) then150 !write(*,*) "euvheat: Error; no OH tracer !!!"173 stop 174 else 175 nspeuv_vgcm=nspeuv_vgcm+1 176 endif 177 g_oh=i_oh 178 if (g_oh.eq.0) then 179 write(*,*) "euvheat: Error; no OH tracer !!!" 151 180 ! write(*,*) "OH must always be present if thermochem=T" 152 !stop153 !else154 ! nespeuv=nespeuv+1155 !endif156 ! g_ho2=igcm_ho2157 !if (g_ho2.eq.0) then158 !write(*,*) "euvheat: Error; no HO2 tracer !!!"181 stop 182 else 183 nspeuv_vgcm=nspeuv_vgcm+1 184 endif 185 g_ho2=i_ho2 186 if (g_ho2.eq.0) then 187 write(*,*) "euvheat: Error; no HO2 tracer !!!" 159 188 ! write(*,*) "HO2 must always be present if thermochem=T" 160 !stop161 !else162 ! nespeuv=nespeuv+1163 !endif164 ! g_h2o2=igcm_h2o2165 !if (g_h2o2.eq.0) then166 !write(*,*) "euvheat: Error; no H2O2 tracer !!!"189 stop 190 else 191 nspeuv_vgcm=nspeuv_vgcm+1 192 endif 193 g_h2o2=i_h2o2 194 if (g_h2o2.eq.0) then 195 write(*,*) "euvheat: Error; no H2O2 tracer !!!" 167 196 ! write(*,*) "H2O2 is always needed if calleuv=.true." 168 !stop169 !else170 ! nespeuv=nespeuv+1171 !endif172 ! g_h2o=igcm_h2o_vap 173 !if (g_h2o.eq.0) then174 !write(*,*) "euvheat: Error; no water vapor tracer !!!"197 stop 198 else 199 nspeuv_vgcm=nspeuv_vgcm+1 200 endif 201 g_h2o=i_h2o 202 if (g_h2o.eq.0) then 203 write(*,*) "euvheat: Error; no water vapor tracer !!!" 175 204 ! write(*,*) "H2O is always needed if calleuv=.true." 176 !stop177 !else178 ! nespeuv=nespeuv+1179 !endif180 ! g_o1d=igcm_o1d181 !if (g_o1d.eq.0) then182 !write(*,*) "euvheat: Error; no O1D tracer !!!"205 stop 206 else 207 nspeuv_vgcm=nspeuv_vgcm+1 208 endif 209 g_o1d=i_o1d 210 if (g_o1d.eq.0) then 211 write(*,*) "euvheat: Error; no O1D tracer !!!" 183 212 ! write(*,*) "O1D must always be present if thermochem=T" 184 !stop185 !else186 ! nespeuv=nespeuv+1187 !endif188 ! g_h=igcm_h189 !if (g_h.eq.0) then190 !write(*,*) "euvheat: Error; no H tracer !!!"213 stop 214 else 215 nspeuv_vgcm=nspeuv_vgcm+1 216 endif 217 g_h=i_h 218 if (g_h.eq.0) then 219 write(*,*) "euvheat: Error; no H tracer !!!" 191 220 ! write(*,*) "H is always needed if calleuv=.true." 192 !stop193 !else194 ! nespeuv=nespeuv+1195 !endif221 stop 222 else 223 nspeuv_vgcm=nspeuv_vgcm+1 224 endif 196 225 197 ! euvmod = 3!Default: C/O/H chemistry226 ! euvmod = 1 !Default: C/O/H chemistry 198 227 ! !Check if O3 is present 199 ! g_o3=igcm_o3200 !if (g_o3.eq.0) then201 !write(*,*) "euvheat: Error; no O3 tracer !!!"202 !write(*,*) "O3 must be present if calleuv=.true."203 !stop204 !else205 ! nespeuv=nespeuv+1206 !euvmod=1207 !endif228 g_o3=i_o3 229 if (g_o3.eq.0) then 230 write(*,*) "euvheat: Error; no O3 tracer !!!" 231 write(*,*) "O3 must be present if calleuv=.true." 232 stop 233 else 234 nspeuv_vgcm=nspeuv_vgcm+1 235 euvmod=1 236 endif 208 237 209 238 !Nitrogen species 210 239 !NO is used to determine if N chemistry is wanted 211 240 !euvmod=2 -> N chemistry 212 ! g_no=igcm_no213 !if (g_no.eq.0) then214 !write(*,*) "euvheat: no NO tracer"215 !write(*,*) "No N species in UV heating"216 !else if(g_no.ne.0) then217 ! nespeuv=nespeuv+1218 !euvmod=2219 !endif241 g_no=i_no 242 if (g_no.eq.0) then 243 write(*,*) "euvheat: no NO tracer" 244 write(*,*) "No N species in UV heating" 245 else if(g_no.ne.0) then 246 nspeuv_vgcm=nspeuv_vgcm+1 247 euvmod=2 248 endif 220 249 ! N 221 ! g_n=igcm_n222 !if(euvmod == 2) then223 !if (g_n.eq.0) then224 !write(*,*) "euvheat: Error; no N tracer !!!"225 !write(*,*) "N needed if NO is in traceur.def"226 !stop227 !else if(g_n.ne.0) then228 ! nespeuv=nespeuv+1229 !endif230 !else231 !if(g_n /= 0) then232 !write(*,*) "euvheat: Error: N present, but NO not!!!"233 !write(*,*) "Both must be in traceur.def"234 !stop235 !endif236 !endif !Of if(euvmod==2)237 !!NO2238 ! g_no2=igcm_no2239 !if(euvmod == 2) then240 !if (g_no2.eq.0) then241 !write(*,*) "euvheat: Error; no NO2 tracer !!!"242 !write(*,*) "NO2 needed if NO is in traceur.def"243 !stop244 !else if(g_no2.ne.0) then245 ! nespeuv=nespeuv+1246 !endif247 !else248 !if(g_no2 /= 0) then249 !write(*,*) "euvheat: Error: NO2 present, but NO not!!!"250 !write(*,*) "Both must be in traceur.def"251 !stop252 !endif253 !endif !Of if(euvmod==2)250 g_n=i_n 251 if(euvmod == 2) then 252 if (g_n.eq.0) then 253 write(*,*) "euvheat: Error; no N tracer !!!" 254 write(*,*) "N needed if NO is in traceur.def" 255 stop 256 else if(g_n.ne.0) then 257 nspeuv_vgcm=nspeuv_vgcm+1 258 endif 259 else 260 if(g_n /= 0) then 261 write(*,*) "euvheat: Error: N present, but NO not!!!" 262 write(*,*) "Both must be in traceur.def" 263 stop 264 endif 265 endif !Of if(euvmod==2) 266 !NO2 267 g_no2=i_no2 268 if(euvmod == 2) then 269 if (g_no2.eq.0) then 270 write(*,*) "euvheat: Error; no NO2 tracer !!!" 271 write(*,*) "NO2 needed if NO is in traceur.def" 272 stop 273 else if(g_no2.ne.0) then 274 nspeuv_vgcm=nspeuv_vgcm+1 275 endif 276 else 277 if(g_no2 /= 0) then 278 write(*,*) "euvheat: Error: NO2 present, but NO not!!!" 279 write(*,*) "Both must be in traceur.def" 280 stop 281 endif 282 endif !Of if(euvmod==2) 254 283 !N2D 255 ! g_n2d=igcm_n2d256 !if(euvmod == 2) then257 !if (g_n2d.eq.0) then258 !write(*,*) "euvheat: Error; no N2D tracer !!!"259 !write(*,*) "N2D needed if NO is in traceur.def"260 !stop261 !else262 ! nespeuv=nespeuv+1263 !endif264 !else265 !if(g_n2d /= 0) then266 !write(*,*) "euvheat: Error: N2D present, but NO not!!!"267 !write(*,*) "Both must be in traceur.def"268 !stop269 !endif270 !endif !Of if(euvmod==2)284 g_n2d=i_n2d 285 if(euvmod == 2) then 286 if (g_n2d.eq.0) then 287 write(*,*) "euvheat: Error; no N2D tracer !!!" 288 write(*,*) "N2D needed if NO is in traceur.def" 289 stop 290 else 291 nspeuv_vgcm=nspeuv_vgcm+1 292 endif 293 else 294 if(g_n2d /= 0) then 295 write(*,*) "euvheat: Error: N2D present, but NO not!!!" 296 write(*,*) "Both must be in traceur.def" 297 stop 298 endif 299 endif !Of if(euvmod==2) 271 300 272 301 !Check if nespeuv is appropriate for the value of euvmod … … 299 328 300 329 301 !Allocate density vector302 allocate(rm(nlev,nespeuv))330 !Allocate density vector 331 allocate(rm(nlev,nespeuv)) 303 332 304 333 firstcall= .false. … … 310 339 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 311 340 341 ! euvmod selection security 342 if(euvmod.gt.2.or.euvmod.lt.0) then 343 write(*,*)'euvheat: bad value for euvmod. Stop' 344 stop 345 endif 346 312 347 ! build local updated values of tracers (if any) and temperature 313 314 348 do l=1,nlev 315 349 do ig=1,nlon 316 317 ! chemical species 318 zq(ig,l,g_co2)=pq(ig,l,g_co2) 319 zq(ig,l,g_co)=pq(ig,l,g_co) 320 zq(ig,l,g_o)=pq(ig,l,g_o) 321 zq(ig,l,g_n2)=pq(ig,l,g_n2) 322 323 ! atmospheric temperature 324 zt(ig,l)=pt(ig,l) 325 350 ! chemical species 351 zq(ig,l,g_co2)=pq(ig,l,g_co2) ! CO2 352 zq(ig,l,g_o2)=pq(ig,l,g_o2) ! O2 353 zq(ig,l,g_o)=pq(ig,l,g_o) ! O(3P) 354 zq(ig,l,g_h2)=pq(ig,l,g_h2) ! H2 355 zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2) ! H2O2 356 zq(ig,l,g_h2o)=pq(ig,l,g_h2o) ! H2O 357 zq(ig,l,g_n2)=pq(ig,l,g_n2) ! N2 358 zq(ig,l,g_co)=pq(ig,l,g_co) ! CO 359 zq(ig,l,g_h)=pq(ig,l,g_h) ! H 360 361 !Only if O3, N or ion chemistry requested 362 if(euvmod.ge.1) then 363 zq(ig,l,g_o3)=pq(ig,l,g_o3) ! 03 364 endif 365 !Only if N or ion chemistry requested 366 if(euvmod.ge.2) then 367 zq(ig,l,g_n)=pq(ig,l,g_n) ! N 368 zq(ig,l,g_no)=pq(ig,l,g_no) ! NO 369 zq(ig,l,g_no2)=pq(ig,l,g_no2) ! NO2 370 endif 371 ! atmospheric temperature 372 zt(ig,l)=pt(ig,l) 326 373 ! write(*,*), "CHECK update densities L332 euv", zq(ig,l,g_co2) 327 328 329 374 enddo 330 375 enddo 331 376 332 !Solar flux calculation 333 377 !Solar flux calculation 334 378 do ig=1,nlon 335 379 zenit=acos(mu0(ig))*180./acos(-1.) !convers from rad to deg 336 380 337 do l=1,nlev 338 339 ! Conversion to number density 340 341 !!! use R specific = R/MolarMass 342 381 do l=1,nlev 382 !Conversion to number density 383 !!! use R specific = R/MolarMass 343 384 dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21 ! [g mol-1] [cm-3] 344 385 345 rm(l,ix_co2) = zq(ig,l,g_co2) * dens / M_tr(g_co2) ! [cm-3] 386 rm(l,ix_co2) = zq(ig,l,g_co2) * dens / M_tr(g_co2) ! [cm-3] 387 rm(l,ix_o2) = zq(ig,l,g_o2) * dens / M_tr(g_o2) 346 388 rm(l,ix_o) = zq(ig,l,g_o) * dens / M_tr(g_o) 347 rm(l,ix_co) = zq(ig,l,g_co) * dens / M_tr(g_co) 348 rm(l,ix_n2) = zq(ig,l,g_n2) * dens / M_tr(g_n2) 389 rm(l,ix_h2) = zq(ig,l,g_h2) * dens / M_tr(g_h2) 390 rm(l,ix_h2o) = zq(ig,l,g_h2o) * dens / M_tr(g_h2o) 391 rm(l,ix_h2o2) = zq(ig,l,g_h2o2)* dens / M_tr(g_h2o2) 392 rm(l,ix_co) = zq(ig,l,g_co) * dens / M_tr(g_co) 393 rm(l,ix_n2) = zq(ig,l,g_n2) * dens / M_tr(g_n2) 394 rm(l,ix_h) = zq(ig,l,g_h) * dens / M_tr(g_h) 395 396 !Only if O3, N or ion chemistry requested 397 if(euvmod.ge.1) then 398 rm(l,ix_o3) = zq(ig,l,g_o3) * dens / M_tr(g_o3) 399 endif 400 !Only if N or ion chemistry requested 401 if(euvmod.ge.2) then 402 rm(l,ix_n) = zq(ig,l,g_n) * dens / M_tr(g_n) 403 rm(l,ix_no) = zq(ig,l,g_no) * dens / M_tr(g_no) 404 rm(l,ix_no2) = zq(ig,l,g_no2) * dens / M_tr(g_no2) 405 endif 349 406 350 407 ! write(*,*), "CHECK n density", l, rm(l,ix_co2) 351 352 353 408 enddo 354 409 … … 368 423 369 424 425 !euveff: UV heating efficiency. Following Fox et al. ASR 1996 426 !should vary between 19% and 23%. Lower values 427 !(i.e. 16%) can be used to compensate 428 !underestimation 429 !of 15-um cooling (see Forget et al. JGR 2009 and 430 !Gonzalez-Galindo et al. JGR 2009) for details 431 370 432 !Calculates the UV heating from the total photoabsorption coefficient 371 433 do l=1,nlev -
trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F
r2686 r2836 42 42 !!! If the values are changed there, the same has to be done here !!! 43 43 44 integer,parameter :: i_co2=145 integer,parameter :: i_n2=1346 integer,parameter :: i_n=1447 integer,parameter :: i_o=348 integer,parameter :: i_co=444 ! integer,parameter :: i_co2=1 45 ! integer,parameter :: i_n2=13 46 ! integer,parameter :: i_n=14 47 ! integer,parameter :: i_o=3 48 ! integer,parameter :: i_co=4 49 49 50 integer,parameter :: ix_co2 = 1 51 integer,parameter :: ix_co = 2 52 integer,parameter :: ix_o = 3 53 integer,parameter :: ix_o1d = 4 54 integer,parameter :: ix_o2 = 5 55 integer,parameter :: ix_o3 = 6 56 integer,parameter :: ix_h = 7 57 integer,parameter :: ix_h2 = 8 58 integer,parameter :: ix_oh = 9 59 integer,parameter :: ix_ho2 = 10 60 integer,parameter :: ix_h2o2 = 11 61 integer,parameter :: ix_h2o = 12 62 integer,parameter :: ix_n = 13 63 integer,parameter :: ix_n2d = 14 64 integer,parameter :: ix_no = 15 65 integer,parameter :: ix_no2 = 16 66 integer,parameter :: ix_n2 = 17 50 67 51 68 c*************************PROGRAM STARTS******************************* … … 73 90 ! 74 91 do i=1,klev 75 xabsi(1,i) = rm(i,i _co2)76 xabsi( 3,i) = rm(i,i_o)77 xabsi( 8,i) = rm(i,i_n2)78 xabsi( 11,i) = rm(i,i_co)79 80 c xabsi(6,i) = rm(i,i_h2o2) 92 xabsi(1,i) = rm(i,ix_co2) ! CO2 93 xabsi(2,i) = rm(i,ix_o2) ! O2 94 xabsi(3,i) = rm(i,ix_o) ! O(3P) 95 xabsi(4,i) = rm(i,ix_h2o) ! H2O 96 xabsi(5,i) = rm(i,ix_h2) ! H2 97 xabsi(6,i) = rm(i,ix_h2o2) ! H2O2 81 98 !Only if O3, N or ion chemistry requested 82 c if(euvmod.ge.1) then 83 c xabsi(7,i) = rm(i,i_o3) 84 c endif 99 if(euvmod.ge.1) then 100 xabsi(7,i) = rm(i,ix_o3) ! O3 101 endif 102 xabsi(8,i) = rm(i,ix_n2) ! N2 85 103 !Only if N or ion chemistry requested 86 c if(euvmod.ge.2) then 87 c xabsi(8,i) = rm(i,i_n2) 88 c xabsi(9,i) = rm(i,i_n) 89 c xabsi(10,i) = rm(i,i_no) 90 c xabsi(13,i) = rm(i,i_no2) 91 c endif 104 if(euvmod.ge.2) then 105 xabsi(9,i) = rm(i,ix_n) ! N 106 xabsi(10,i) = rm(i,ix_no) ! NO 107 xabsi(13,i) = rm(i,ix_no2) ! NO2 108 endif 109 xabsi(11,i) = rm(i,ix_co) ! CO 110 xabsi(12,i) = rm(i,ix_h) ! H 92 111 end do 93 112 -
trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F
r2809 r2836 90 90 real*8 limup ! "" 91 91 92 !!!ATTENTION. Here i _co2 has to have the same value than in euvheat.F9092 !!!ATTENTION. Here ix_co2 has to have the same value than in euvheat.F90 93 93 !!!If the value is changed there, if has to be changed also here !!!! 94 integer,parameter :: i _co2=194 integer,parameter :: ix_co2=1 95 95 96 96 character*20 modname … … 123 123 !Calculation of column amounts 124 124 call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit, 125 $ co2colx,o3pcolx,n2colx,cocolx) 125 $ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx, 126 $ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx) 126 127 127 128 !Auxiliar column to include the temperature dependence … … 130 131 do i=nlayer-1,1,-1 131 132 coltemp(i)=!coltemp(i+1)+ PQ SE ELIMINA? REVISAR 132 $ ( rm(i,i _co2) + rm(i+1,i_co2) ) * 0.5133 $ ( rm(i,ix_co2) + rm(i+1,ix_co2) ) * 0.5 133 134 $ * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i)) 134 135 end do … … 1115 1116 1116 1117 subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit, 1117 $ co2colx,o3pcolx, n2colx, cocolx) 1118 1118 $ co2colx,o2colx,o3pcolx,h2colx,h2ocolx,h2o2colx,o3colx, 1119 $ n2colx,ncolx,nocolx,cocolx,hcolx,no2colx) 1120 1121 c Jun 2022 AM readapted to Venus GCM 1119 1122 c mar 2014 gg adapted to Venus GCM 1120 1123 c nov 2002 fgg first version … … 1135 1138 integer ig 1136 1139 integer chemthermod 1137 integer nesptherm !# of species undergoing 1138 chemistry, input 1140 integer nesptherm !# of species undergoing chemistry, input 1139 1141 real rm(klev,nesptherm) !densities (cm-3), input 1140 1142 real tx(klev) !temperature profile, input … … 1142 1144 real zenit !SZA, input 1143 1145 real co2colx(klev) !column density of CO2 (cm^-2), output 1146 real o2colx(klev) !column density of O2(cm^-2), output 1144 1147 real o3pcolx(klev) !column density of O(3P)(cm^-2), output 1148 real h2colx(klev) !H2 column density (cm-2), output 1149 real h2ocolx(klev) !H2O column density (cm-2), output 1150 real h2o2colx(klev) !column density of H2O2(cm^-2), output 1151 real o3colx(klev) !O3 column density (cm-2), output 1145 1152 real n2colx(klev) !N2 column density (cm-2), output 1153 real ncolx(klev) !N column density (cm-2), output 1154 real nocolx(klev) !NO column density (cm-2), output 1146 1155 real cocolx(klev) !CO column density (cm-2), output 1156 real hcolx(klev) !H column density (cm-2), output 1157 real no2colx(klev) !NO2 column density (cm-2), output 1158 1147 1159 1148 1160 c local variables … … 1153 1165 real Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2 1154 1166 1167 ! density 1155 1168 real co2x(klev) 1169 real o2x(klev) 1156 1170 real o3px(klev) 1157 1171 real cox(klev) 1172 real hx(klev) 1173 real h2x(klev) 1174 real h2ox(klev) 1175 real h2o2x(klev) 1176 real o3x(klev) 1158 1177 real n2x(klev) 1159 1178 real nx(klev) 1179 real nox(klev) 1180 real no2x(klev) 1160 1181 1161 1182 integer i,j,k,icol,indexint !indexes … … 1179 1200 !!! If the values are changed there, the same has to be done here !!! 1180 1201 1181 integer,parameter :: ix_co2=1 1182 integer,parameter :: ix_n2=13 1183 integer,parameter :: ix_o=3 1184 integer,parameter :: ix_co=4 1202 ! integer,parameter :: ix_co2=1 1203 ! integer,parameter :: ix_n2=13 1204 ! integer,parameter :: ix_o=3 1205 ! integer,parameter :: ix_co=4 1206 1207 integer,parameter :: ix_co2 = 1 1208 integer,parameter :: ix_co = 2 1209 integer,parameter :: ix_o = 3 1210 integer,parameter :: ix_o1d = 4 1211 integer,parameter :: ix_o2 = 5 1212 integer,parameter :: ix_o3 = 6 1213 integer,parameter :: ix_h = 7 1214 integer,parameter :: ix_h2 = 8 1215 integer,parameter :: ix_oh = 9 1216 integer,parameter :: ix_ho2 = 10 1217 integer,parameter :: ix_h2o2 = 11 1218 integer,parameter :: ix_h2o = 12 1219 integer,parameter :: ix_n = 13 1220 integer,parameter :: ix_n2d = 14 1221 integer,parameter :: ix_no = 15 1222 integer,parameter :: ix_no2 = 16 1223 integer,parameter :: ix_n2 = 17 1185 1224 1186 1225 c*************************PROGRAM STARTS******************************* … … 1195 1234 xx = kboltzman * tx(klev) * n_avog / grav(klev) ! g cm mol-1 1196 1235 1197 Ho3p = xx / mmolo1198 1236 Hco2 = xx / mmolco2 1237 Ho2 = xx / mmolo2 1238 Ho3p = xx / mmolo ! Oxygen 3P 1239 Hh2 = xx / mmolh2 1240 Hh2o2 = xx / mmolh2o2 1241 Hh2o = xx / mmolh2o 1242 1243 !Only if O3 chem. required 1244 if(chemthermod.ge.1) 1245 $ Ho3 = xx / mmolo3 1246 Hn2 = xx / mmoln2 1247 !Only if N or ion chem. 1248 if(chemthermod.ge.2) then 1249 Hn = xx / mmoln 1250 Hno = xx / mmolno 1251 Hno2 = xx / mmolno2 1252 endif 1199 1253 Hco = xx / mmolco 1200 Hn2 = xx / mmoln2 1201 Hn = xx / mmoln 1202 1203 !Only if O3 chem. required 1204 c if(chemthermod.ge.1) 1205 ! $ Ho3 = xx / mmol(igcm_o3) 1206 c $ Ho3 = xx / mmolo3 1207 c !Only if N or ion chem. 1208 c if(chemthermod.ge.2) then 1209 c Hn2 = xx / mmoln2 1210 c Hn = xx / mmoln 1211 c Hno = xx / mmolno 1212 c Hno2 = xx / mmolno2 1213 c endif 1254 Hh = xx / mmolh 1214 1255 ! first loop in altitude : initialisation 1215 1256 do i=klev,1,-1 1216 1257 !Column initialisation 1217 1258 co2colx(i) = 0. 1259 o2colx(i) = 0. 1218 1260 o3pcolx(i) = 0. 1261 h2colx(i) = 0. 1262 h2ocolx(i) = 0. 1263 h2o2colx(i) = 0. 1264 o3colx(i) = 0. 1219 1265 n2colx(i) = 0. 1266 ncolx(i) = 0. 1267 nocolx(i) = 0. 1220 1268 cocolx(i) = 0. 1221 1269 hcolx(i) = 0. 1270 no2colx(i) = 0. 1222 1271 !--Densities [cm-3] 1223 1272 co2x(i) = rm(i,ix_co2) 1273 o2x(i) = rm(i,ix_o2) 1224 1274 o3px(i) = rm(i,ix_o) 1275 h2x(i) = rm(i,ix_h2) 1276 h2ox(i) = rm(i,ix_h2o) 1277 h2o2x(i) = rm(i,ix_h2o2) 1225 1278 cox(i) = rm(i,ix_co) 1279 hx(i) = rm(i,ix_h) ! write(*,*), '--jthermcalc--', co2x(i) 1280 1281 !Only if O3 chem. required 1282 if(chemthermod.ge.1) 1283 $ o3x(i) = rm(i,ix_o3) 1226 1284 n2x(i) = rm(i,ix_n2) 1227 1228 c write(*,*), '--jthermcalc--', co2x(i) 1229 1230 !Only if O3 chem. required 1231 c if(chemthermod.ge.1) 1232 c $ o3x(i) = rm(i,i_o3) 1233 c !Only if Nitrogen of ion chemistry requested 1234 c if(chemthermod.ge.2) then 1235 c n2x(i) = rm(i,i_n2) 1236 c nx(i) = rm(i,i_n) 1237 c nox(i) = rm(i,i_no) 1238 c no2x(i) = rm(i,i_no2) 1239 c endif 1285 !Only if Nitrogen of ion chemistry requested 1286 if(chemthermod.ge.2) then 1287 nx(i) = rm(i,ix_n) 1288 nox(i) = rm(i,ix_no) 1289 no2x(i) = rm(i,ix_no2) 1290 endif 1240 1291 enddo ! end first loop 1241 1292 … … 1247 1298 if(ilayesp(nlayesp).eq.-1) then 1248 1299 co2colx(i)=1.e25 1300 o2colx(i)=1.e25 1249 1301 o3pcolx(i)=1.e25 1302 h2colx(i)=1.e25 1303 h2ocolx(i)=1.e25 1304 h2o2colx(i)=1.e25 1305 o3colx(i)=1.e25 1250 1306 n2colx(i)=1.e25 1251 1307 cocolx(i)=1.e25 1252 1253 c o2colx(i)=1.e25 1254 c o3pcolx(i)=1.e25 1255 c h2colx(i)=1.e25 1256 c h2ocolx(i)=1.e25 1257 c h2o2colx(i)=1.e25 1258 c o3colx(i)=1.e25 1259 c ncolx(i)=1.e25 1260 c nocolx(i)=1.e25 1261 1262 c cocolx(i)=1.e25 1263 c hcolx(i)=1.e25 1264 c no2colx(i)=1.e25 1308 hcolx(i)=1.e25 1309 ncolx(i)=1.e25 1310 nocolx(i)=1.e25 1311 no2colx(i)=1.e25 1265 1312 else 1266 1313 rcmnz = ( radio + iz(klev) ) * 1.e5 ! km --> cm 1267 1314 rcmmini = ( radio + zmini ) * 1.e5 1268 !Column calculation taking into account the geometrical 1269 !depth 1315 !Column calculation taking into account the geometrical depth 1270 1316 !calculated before 1271 1317 do j=1,nlayesp … … 1278 1324 co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j) 1279 1325 $ *1.e-5 1326 h2o2colx(i)=h2o2colx(i)+ 1327 $ h2o2x(klev)*Hh2o2*esp(j)*1.e-5 1328 o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j) 1329 $ *1.e-5 1330 h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j) 1331 $ *1.e-5 1332 h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j) 1333 $ *1.e-5 1334 n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j) 1335 $ *1.e-5 1280 1336 cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j) 1281 1337 $ *1.e-5 1282 n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)1338 hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j) 1283 1339 $ *1.e-5 1284 1340 1285 c h2o2colx(i)=h2o2colx(i)+1286 c $ h2o2x(klev)*Hh2o2*esp(j)*1.e-51287 c o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)1288 c $ *1.e-51289 c h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)1290 c $ *1.e-51291 c h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)1292 c $ *1.e-51293 c cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)1294 c $ *1.e-51295 c hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)1296 c $ *1.e-51297 1341 !Only if O3 chemistry required 1298 cif(chemthermod.ge.1) o3colx(i)=1299 c$ o3colx(i)+o3x(klev)*Ho3*esp(j)1300 c$ *1.e-51342 if(chemthermod.ge.1) o3colx(i)= 1343 $ o3colx(i)+o3x(klev)*Ho3*esp(j) 1344 $ *1.e-5 1301 1345 !Only if N or ion chemistry requested 1302 c if(chemthermod.ge.2) then 1303 c n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j) 1304 c $ *1.e-5 1305 1306 c endif 1346 if(chemthermod.ge.2) then 1347 ncolx(i)=ncolx(i)+nx(klev)*Hn*esp(j) 1348 $ *1.e-5 1349 nocolx(i)=nocolx(i)+nox(klev)*Hno*esp(j) 1350 $ *1.e-5 1351 no2colx(i)=no2colx(i)+no2x(klev)*Hno2*esp(j) 1352 $ *1.e-5 1353 endif 1307 1354 else if(zenit.gt.60.) then 1308 1355 espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j) 1356 espo2 = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j) 1309 1357 espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j) 1358 esph2 = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j) 1359 esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j) 1360 esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j) 1310 1361 espco = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j) 1311 1362 espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j) 1312 espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) 1313 1314 c espo2 = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j) 1315 c esph2 = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j) 1316 c esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j) 1317 c esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j) 1318 c esph = sqrt((rcmnz+Hh)**2 -rcmmini**2) - esp(j) 1363 esph = sqrt((rcmnz+Hh)**2 -rcmmini**2) - esp(j) 1364 1319 1365 !Only if O3 chemistry required 1320 c if(chemthermod.ge.1) 1321 c $ espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j) 1322 c !Only if N or ion chemistry requested 1323 c if(chemthermod.ge.2) then 1324 c espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j) 1325 c espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) 1326 c espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j) 1327 c espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j) 1328 c endif 1366 if(chemthermod.ge.1) 1367 $ espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j) 1368 !Only if N or ion chemistry requested 1369 if(chemthermod.ge.2) then 1370 espn =sqrt((rcmnz+Hn)**2-rcmmini**2) - esp(j) 1371 espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j) 1372 espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j) 1373 endif 1329 1374 1330 1375 co2colx(i) = co2colx(i) + espco2*co2x(klev) 1376 o2colx(i) = o2colx(i) + espo2*o2x(klev) 1331 1377 o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev) 1378 n2colx(i) = n2colx(i) + espn2*n2x(klev) 1379 h2colx(i) = h2colx(i) + esph2*h2x(klev) 1380 h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev) 1381 h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev) 1332 1382 cocolx(i) = cocolx(i) + espco*cox(klev) 1333 n2colx(i) = n2colx(i) + espn2*n2x(klev) 1334 1335 c o2colx(i) = o2colx(i) + espo2*o2x(klev) 1336 c h2colx(i) = h2colx(i) + esph2*h2x(klev) 1337 c h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev) 1338 c h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev) 1339 c cocolx(i) = cocolx(i) + espco*cox(klev) 1340 c hcolx(i) = hcolx(i) + esph*hx(klev) 1383 hcolx(i) = hcolx(i) + esph*hx(klev) 1384 1341 1385 !Only if O3 chemistry required 1342 c if(chemthermod.ge.1) 1343 c $ o3colx(i) = o3colx(i) + espo3*o3x(klev) 1344 c !Only if N or ion chemistry requested 1345 c if(chemthermod.ge.2) then 1346 c n2colx(i) = n2colx(i) + espn2*n2x(klev) 1347 c ncolx(i) = ncolx(i) + espn*nx(klev) 1348 c nocolx(i) = nocolx(i) + espno*nox(klev) 1349 c no2colx(i) = no2colx(i) + espno2*no2x(klev) 1350 c endif 1386 if(chemthermod.ge.1) 1387 $ o3colx(i) = o3colx(i) + espo3*o3x(klev) 1388 !Only if N or ion chemistry requested 1389 if(chemthermod.ge.2) then 1390 ncolx(i) = ncolx(i) + espn*nx(klev) 1391 nocolx(i) = nocolx(i) + espno*nox(klev) 1392 no2colx(i) = no2colx(i) + espno2*no2x(klev) 1393 endif 1351 1394 endif !Of if zenit.lt.60 1352 1395 !Other layers 1353 1396 else 1354 co2colx(i) = co2colx(i) + 1397 co2colx(i) = co2colx(i) + 1355 1398 $ esp(j) * (co2x(jj)+co2x(jj+1)) / 2. 1356 o3pcolx(i) = o3pcolx(i) + 1399 o2colx(i) = o2colx(i) + 1400 $ esp(j) * (o2x(jj)+o2x(jj+1)) / 2. 1401 o3pcolx(i) = o3pcolx(i) + 1357 1402 $ esp(j) * (o3px(jj)+o3px(jj+1)) / 2. 1358 cocolx(i) = cocolx(i) + 1359 $ esp(j) * (cox(jj)+cox(jj+1)) / 2. 1403 h2colx(i) = h2colx(i) + 1404 $ esp(j) * (h2x(jj)+h2x(jj+1)) / 2. 1405 h2ocolx(i) = h2ocolx(i) + 1406 $ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2. 1407 h2o2colx(i) = h2o2colx(i) + 1408 $ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2. 1360 1409 n2colx(i) = n2colx(i) + 1361 1410 $ esp(j) * (n2x(jj)+n2x(jj+1)) / 2. 1362 1363 c 1364 c o2colx(i) = o2colx(i) + 1365 c $ esp(j) * (o2x(jj)+o2x(jj+1)) / 2. 1366 c h2colx(i) = h2colx(i) + 1367 c $ esp(j) * (h2x(jj)+h2x(jj+1)) / 2. 1368 c h2ocolx(i) = h2ocolx(i) + 1369 c $ esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2. 1370 c h2o2colx(i) = h2o2colx(i) + 1371 c $ esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2. 1372 c hcolx(i) = hcolx(i) + 1373 c $ esp(j) * (hx(jj)+hx(jj+1)) / 2. 1411 cocolx(i) = cocolx(i) + 1412 $ esp(j) * (cox(jj)+cox(jj+1)) / 2. 1413 hcolx(i) = hcolx(i) + 1414 $ esp(j) * (hx(jj)+hx(jj+1)) / 2. 1415 1374 1416 !Only if O3 chemistry required 1375 c if(chemthermod.ge.1) 1376 c $ o3colx(i) = o3colx(i) + 1377 c $ esp(j) * (o3x(jj)+o3x(jj+1)) / 2. 1378 c !Only if N or ion chemistry requested 1379 c if(chemthermod.ge.2) then 1380 c n2colx(i) = n2colx(i) + 1381 c $ esp(j) * (n2x(jj)+n2x(jj+1)) / 2. 1382 c ncolx(i) = ncolx(i) + 1383 c $ esp(j) * (nx(jj)+nx(jj+1)) / 2. 1384 c nocolx(i) = nocolx(i) + 1385 c $ esp(j) * (nox(jj)+nox(jj+1)) / 2. 1386 c no2colx(i) = no2colx(i) + 1387 c $ esp(j) * (no2x(jj)+no2x(jj+1)) / 2. 1388 c endif 1389 1417 if(chemthermod.ge.1) 1418 $ o3colx(i) = o3colx(i) + 1419 $ esp(j) * (o3x(jj)+o3x(jj+1)) / 2. 1420 !Only if N or ion chemistry requested 1421 if(chemthermod.ge.2) then 1422 ncolx(i) = ncolx(i) + 1423 $ esp(j) * (nx(jj)+nx(jj+1)) / 2. 1424 nocolx(i) = nocolx(i) + 1425 $ esp(j) * (nox(jj)+nox(jj+1)) / 2. 1426 no2colx(i) = no2colx(i) + 1427 $ esp(j) * (no2x(jj)+no2x(jj+1)) / 2. 1428 endif 1390 1429 endif !Of if jj.eq.klev 1391 1430 end do !Of do j=1,nlayesp … … 1406 1445 C escout(nlayer) on pressure grid p(nlayer). 1407 1446 C 1408 real*8 wm(nlayer),wp(nlayer),p(nlayer) 1409 integer nm(nlayer) 1410 real*8 pin(nl) 1411 real*8 limup,limdown 1412 integer nl,nlayer,n1,n,np,nini 1447 real*8,intent(out) :: wm(nlayer),wp(nlayer) ! interpolation weights 1448 integer,intent(out) :: nm(nlayer) ! index of nearest point 1449 real*8,intent(in) :: pin(nl),p(nlayer) 1450 real*8,intent(in) :: limup,limdown 1451 integer,intent(in) :: nl,nlayer 1452 integer :: n1,n,np,nini 1413 1453 nini=1 1414 1454 do n1=1,nlayer … … 1451 1491 real szadeg ! I. SZA [rad] 1452 1492 real z ! I. altitude of interest [km] 1453 integer nz3,ig ! I. dimension of esp, ylayesp, etc... 1454 ! (=2*klev= max# of layers in 1455 ! ray path) 1456 real iz(klev+1) ! I. Altitude of each layer 1493 integer nz3,ig ! I. dimension of esp, ylayesp, etc... 1494 ! (=2*klev= max# of layers in ray path) 1495 real iz(klev+1) ! I. Altitude of each layer 1457 1496 real*8 esp(nz3) ! O. layer widths after geometrically 1458 ! amplified; in [cm] except 1459 ! at TOA 1460 ! where an auxiliary value is 1461 ! used 1497 ! amplified; in [cm] except at TOA 1498 ! where an auxiliary value is used 1462 1499 real*8 ilayesp(nz3) ! O. Indexes of layers along ray path 1463 1500 real*8 szalayesp(nz3) ! O. SZA [deg] " " " 1464 1501 integer nlayesp 1465 ! real*8 nlayesp ! O. # layers along ray path at 1466 ! this z 1502 ! real*8 nlayesp ! O. # layers along ray path at this z 1467 1503 real*8 zmini ! O. Minimum altitud of ray path [km] 1468 1504 … … 1501 1537 ! The optical thickness will be given by 1/cos(sza) 1502 1538 ! We deal with 2 different regions: 1503 ! 1: First, all layers between z and ztop ("upper part of 1504 ! ray") 1539 ! 1: First, all layers between z and ztop ("upper part of ray") 1505 1540 ! 2: Second, the layer at ztop 1506 1541 if(szadeg.lt.60.d0) then … … 1529 1564 ! The optical thickness is evaluated. 1530 1565 ! (the magnitude of the effect of not using cos(sza) is 3.e-5 1531 ! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, 1532 ! approximately) 1566 ! for z=60km & sza=30, and 5e-4 for z=60km & sza=60, approximately) 1533 1567 ! We deal with 2 different regions: 1534 ! 1: First, all layers between z and ztop ("upper part of 1535 ! ray") 1568 ! 1: First, all layers between z and ztop ("upper part of ray") 1536 1569 ! 2: Second, the layer at ztop ("uppermost layer") 1537 1570 else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then … … 1553 1586 rkmj = radio+diz(j) 1554 1587 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1555 szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592 1556 ! [deg] 1588 szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg] 1557 1589 end do 1558 1590 1470 continue … … 1603 1635 rkmj = radio+diz(j) 1604 1636 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1605 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 1606 ! [deg] 1637 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg] 1607 1638 end do 1608 1639 … … 1629 1660 rkmj = radio+diz(j+1) 1630 1661 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1631 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 1632 ! [deg] 1662 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg] 1633 1663 endif 1634 1664 … … 1644 1674 rkmj = radio+diz(j) 1645 1675 szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad] 1646 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 1647 ! [deg] 1676 szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592! [deg] 1648 1677 end do 1649 1678 … … 1655 1684 esp(nlayesp) = 1656 1685 $ sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) 1657 $ - sqrt( (radio+diz(j))**2 - rkmmini**2 ) 1658 ! [km] 1659 esp(nlayesp) = esp(nlayesp) * 1.d5 1660 ! [cm] 1686 $ - sqrt( (radio+diz(j))**2 - rkmmini**2 )! [km] 1687 esp(nlayesp) = esp(nlayesp) * 1.d5! [cm] 1661 1688 rkmj = radio+diz(j) 1662 szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj ) 1663 ! [rad] 1664 szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592 1665 ! [deg] 1689 szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )! [rad] 1690 szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592! [deg] 1666 1691 end do 1667 1692 -
trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90
r2795 r2836 421 421 ! DIffusion coefficient 422 422 CALL DCOEFF(nn,dij,Praf,Traf,Nraf,Nrafk,Draf,nlraf,ncompdiff) 423 !Draf(:) = max(100.,Draf(:)) 424 !Draf(:) = 0.01*Draf(:) 423 425 Drafmol(:,nn)=Draf(:) 424 426 ! Scale height of the density of the specie -
trunk/LMDZ.VENUS/libf/phyvenus/nonoro_gwd_ran_mod.F90
r2686 r2836 42 42 43 43 USE ioipsl_getin_p_mod, ONLY : getin_p 44 use geometry_mod, only: cell_area44 USE geometry_mod, ONLY: cell_area, latitude_deg 45 45 !#ifdef CPP_XIOS 46 46 ! use xios_output_mod, only: send_xios_field … … 107 107 REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity 108 108 109 !!! Nominal as Gilli2017 110 REAL, parameter:: epflux_max = 0.005 ! Max EP flux value at launching altitude (previous name: RUWMAX) 109 !---------------------------------------- 110 ! Nominal as Gilli2017 111 ! REAL, parameter:: epflux_max = 0.005 ! Max EP flux value at launching altitude (previous name: RUWMAX) 112 !---------------------------------------- 113 ! VCD 2.1 tuning 114 REAL, parameter:: epflux_max = 0.001 ! Max EP flux value at launching altitude (previous name: RUWMAX) 115 !---------------------------------------- 111 116 112 117 REAL cpha ! absolute phase velocity frequency … … 251 256 DO jw = 1, nw 252 257 DO ii = 1, ngrid 253 258 254 259 ran_num_1 = mod(tt(ii, jw) * 10., 1.) 255 260 ran_num_2 = mod(tt(ii, jw) * 100., 1.) 256 261 257 262 !! OPTIONS GENERIC DIFF VENUS !! 258 ! angle (random) 259 zp(jw, ii) = (sign(1., 0.5 - ran_num_1) &263 ! angle (random) - reference 264 zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & 260 265 + 1.) * RPI / 2. 261 ! angle (0 or RPI) 262 !zp(jw, ii) = RPI*mod(jw,2) 266 ! --------- TRY 00 ----------------------- 267 !IF((abs(latitude_deg(ii)) .le. 55.)) THEN 268 ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & 269 ! * (10.) ! +/- 10 deg par rapport equateur 270 !ELSE IF ((abs(latitude_deg(ii)) .le. 75.)) THEN 271 ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & 272 ! * (10. +(90.-10.) & 273 ! * (latitude_deg(ii)-55.)/20.) 274 !ELSE 275 ! zp(jw, ii) = ran_num_1 * RPI 276 !ENDIF 277 !zp(jw, ii) = mod(zp(jw, ii),2.*RPI) 278 ! ------ TRY 01------------------------------- 279 !IF((abs(latitude_deg(ii)) .le. 55)) THEN 280 ! zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & 281 ! + 1.) * RPI / 2. 282 !ELSE 283 ! zp(jw, ii) = ran_num_1 * RPI 284 !ENDIF 285 ! ---- angle (0 or RPI) ----- 286 !!zp(jw, ii) = RPI*mod(jw,2) 263 287 264 288 ! horizontal wavenumber amplitude -
trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F
r2464 r2836 47 47 write(*,*)'It should be in :', trim(datafile),'/' 48 48 write(*,*)'1) You can change this directory address in ' 49 write(*,*)' callphys.def with data dir=/path/to/dir'49 write(*,*)' callphys.def with datafile=/path/to/dir' 50 50 write(*,*)'2) If necessary, EUVDAT (and other datafiles)' 51 51 write(*,*)' can be obtained online on:' … … 279 279 c dissociation and ionization efficiencies 280 280 281 ! do inter=1,ninter 282 ! efdisco2(inter)=0. 283 ! efdiso2(inter)=0. 284 ! efdish2(inter)=0. 285 ! efdish2o(inter)=0. 286 ! efdish2o2(inter)=0. 287 ! efdiso3(inter)=0. 288 ! efdisco(inter)=0. 289 ! efdisn2(inter)=0. 290 ! efdisno(inter)=0. 291 ! efdisno2(inter)=0. 292 ! efionco2(inter,1)=0. 293 ! efionco2(inter,2)=0. 294 ! efionco2(inter,3)=0. 295 ! efionco2(inter,4)=0. 296 ! efiono3p(inter)=0. 297 ! efionn2(inter,1)=0. 298 ! efionn2(inter,2)=0. 299 ! efionco(inter,1)=0. 300 ! efionco(inter,2)=0. 301 ! efionco(inter,3)=0. 302 ! efionn(inter)=0. 303 ! efionh(inter)=0. 304 ! efionno(inter)=0. 305 ! enddo 281 do inter=1,ninter 282 efdisco2(inter)=0. 283 efdiso2(inter)=0. 284 efdish2(inter)=0. 285 efdish2o(inter)=0. 286 efdish2o2(inter)=0. 287 efdiso3(inter)=0. 288 efdisco(inter)=0. 289 efdisn2(inter)=0. 290 efdisno(inter)=0. 291 efdisno2(inter)=0. 292 efionco2(inter,1)=0. 293 efionco2(inter,2)=0. 294 efionco2(inter,3)=0. 295 efionco2(inter,4)=0. 296 efiono2(inter,1)=0. 297 efiono2(inter,2)=0. 298 efiono3p(inter)=0. 299 efionn2(inter,1)=0. 300 efionn2(inter,2)=0. 301 efionco(inter,1)=0. 302 efionco(inter,2)=0. 303 efionco(inter,3)=0. 304 efionn(inter)=0. 305 efionh(inter)=0. 306 efionno(inter)=0. 307 enddo 306 308 307 309 … … 320 322 ! efdisno(inter)=1. 321 323 ! enddo 324 ! close(120) 322 325 323 326 324 327 c N2 325 328 326 !efdisn2(15)=0.1327 !do inter=16,ninter328 !efdisn2(inter)=1.329 !enddo329 efdisn2(15)=0.1 330 do inter=16,ninter 331 efdisn2(inter)=1. 332 enddo 330 333 331 334 332 335 c CO 333 336 334 !efdisco(16)=0.5335 !do inter=17,ninter336 !efdisco(inter)=1.337 !enddo337 efdisco(16)=0.5 338 do inter=17,ninter 339 efdisco(inter)=1. 340 enddo 338 341 339 342 340 343 c O, N, H 341 344 342 !do inter=1,ninter343 !efdiso(inter)=0.344 !efdisn(inter)=0.345 !efdish(inter)=0.346 !enddo345 do inter=1,ninter 346 efdiso(inter)=0. 347 efdisn(inter)=0. 348 efdish(inter)=0. 349 enddo 347 350 348 351 349 352 c H2O, H2O2, O3, NO2 350 353 351 !do inter=25,31352 !efdish2o(inter)=1.353 !enddo354 !do inter=25,35355 !efdish2o2(inter)=1.356 !enddo357 !do inter=34,36358 !efdiso3(inter)=1.359 !enddo360 !do inter=27,36361 !efdisno2(inter)=1.362 !enddo363 !do inter=1,15364 !efdish2(inter)=1.365 !enddo354 do inter=25,31 355 efdish2o(inter)=1. 356 enddo 357 do inter=25,35 358 efdish2o2(inter)=1. 359 enddo 360 do inter=34,36 361 efdiso3(inter)=1. 362 enddo 363 do inter=27,36 364 efdisno2(inter)=1. 365 enddo 366 do inter=1,15 367 efdish2(inter)=1. 368 enddo 366 369 367 370 !4 possible channels for CO2 ionization 371 open(130,file=trim(datafile)//'/EUVDAT'// 372 $ '/co2ion_branchingratio_schunkandnagy2000_param.dat') 373 do inter=1,16 374 read(130,*)i,nada,efionco2(inter,1),efionco2(inter,2), 375 $ efionco2(inter,3),efionco2(inter,4) 376 !Multiply the relative efficiency of each channel by the total ionization efficiency (second column) 377 efdisco2(inter)=1.-nada 378 efionco2(inter,1)=(1.-efdisco2(inter))*efionco2(inter,1) 379 efionco2(inter,2)=(1.-efdisco2(inter))*efionco2(inter,2) 380 efionco2(inter,3)=(1.-efdisco2(inter))*efionco2(inter,3) 381 efionco2(inter,4)=(1.-efdisco2(inter))*efionco2(inter,4) 382 enddo 383 close(130) 384 385 do inter=17,36 386 efdisco2(inter)=1. 387 enddo 388 368 389 ! do inter=14,16 369 390 ! efionco2(inter,1)=1.-efdisco2(inter) … … 380 401 ! enddo 381 402 403 !2 possible channels for O2 ionization 404 open(131,file=trim(datafile)//'/EUVDAT'// 405 $ '/o2ion_branchingratio_schunkandnagy2000_param.dat') 406 do inter=1,23 407 read(131,*)i,nada,efiono2(inter,1),efiono2(inter,2) 408 !Multiply the relative efficiency of each channel by the total ionization efficiency (second column) 409 efdiso2(inter)=1.-nada 410 efiono2(inter,1)=(1.-efdiso2(inter))*efiono2(inter,1) 411 efiono2(inter,2)=(1.-efdiso2(inter))*efiono2(inter,2) 412 enddo 413 close(131) 414 415 do inter=24,36 416 efdiso2(inter)=1. 417 enddo 418 419 382 420 !For O(3p) total ionization under 91.1 nm 383 !do inter=1,16384 !efiono3p(inter)=1.d0385 !enddo421 do inter=1,16 422 efiono3p(inter)=1.d0 423 enddo 386 424 387 425 !2 channels for N2 ionization 426 open(132,file=trim(datafile)//'/EUVDAT'// 427 $ '/n2ion_branchingratio_schunkandnagy2000_param.dat') 428 do inter=1,15 429 read(132,*)i,nada,efionn2(inter,1),efionn2(inter,2) 430 !Multiply the relative efficiency of each channel by the total ionization efficiency (second column) 431 efdisn2(inter)=1.-nada 432 efionn2(inter,1)=(1.-efdisn2(inter))*efionn2(inter,1) 433 efionn2(inter,2)=(1.-efdisn2(inter))*efionn2(inter,2) 434 enddo 435 close(132) 436 437 do inter=16,36 438 efdisn2(inter)=1. 439 enddo 440 388 441 ! do inter=9,15 389 442 ! efionn2(inter,1)=1.-efdisn2(inter) … … 394 447 395 448 !3 channels for CO ionization 449 open(133,file=trim(datafile)//'/EUVDAT'// 450 $ '/coion_branchingratio_schunkandnagy2000_param.dat') 451 do inter=1,16 452 read(133,*)i,nada,efionco(inter,1),efionco(inter,2), 453 $ efionco(inter,3) 454 !Multiply the relative efficiency of each channel by the total ionization efficiency (second column) 455 efdisco(inter)=1.-nada 456 efionco(inter,1)=(1.-efdisco(inter))*efionco(inter,1) 457 efionco(inter,2)=(1.-efdisco(inter))*efionco(inter,2) 458 efionco(inter,3)=(1.-efdisco(inter))*efionco(inter,3) 459 enddo 460 close(133) 461 462 463 464 do inter=17,36 465 efdisco(inter)=1. 466 enddo 467 396 468 ! do inter=11,16 397 469 ! efionco(inter,1)=1.-efdisco(inter) … … 408 480 ! enddo 409 481 482 410 483 !Total ionization under 85 nm for N 411 !do inter=1,16412 !efionn(inter)=1.413 !enddo484 do inter=1,16 485 efionn(inter)=1. 486 enddo 414 487 415 488 !NO 416 !do inter=2,28417 !efionno(inter)=1.-efdisno(inter)418 !enddo489 do inter=2,28 490 efionno(inter)=1.-efdisno(inter) 491 enddo 419 492 420 493 !Total ionization under 90 nm for H 421 !do inter=3,16422 !efionh(inter)=1.423 !enddo494 do inter=3,16 495 efionh(inter)=1. 496 enddo 424 497 425 498 -
trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90
r2795 r2836 1 subroutine photochemistry_venus(nz, n_lon, zlocal, ptimestep, p, t, tr, mumean, sza_input, lon, lat, nesp, iter, prod_tr, loss_tr, em_no, em_o2) 1 !**************************************************************** 2 ! 3 ! Photochemical routine 4 ! 5 ! Authors: Franck Lefevre, Francisco Gonzalez-Galindo 6 ! ------- 7 ! 8 ! Version: 14/11/2020 9 ! 10 ! ASIS scheme : for details on the method see 11 ! Cariolle et al., Geosci. Model Dev., 10, 1467-1485, 2017. 12 ! 13 ! ------- 14 ! 15 ! 2022/09/15: adding ion chemistry by Antoine Martinez 16 ! 17 !***************************************************************** 18 19 subroutine photochemistry_venus(nz, n_lon, zlocal, ptimestep, & 20 ok_jonline,ok_ionchem,tuneupperatm, & 21 nb_reaction_3_max,nb_reaction_4_max, & 22 nb_phot_max,nphotion,ig, & 23 p, t, t_elect, tr, vmr_dens_euv, mumean, & 24 sza_input, lon, lat, nesp, nespeuv, iter, & 25 prod_tr, loss_tr, em_no, em_o2) 2 26 3 27 use chemparam_mod 4 28 use photolysis_mod 29 use param_v4_h, only: jion 30 use iono_h, only: phdisrate 31 !nphot ! number of photodissociations, include in "use photolysis_mod" 5 32 6 33 implicit none … … 10 37 !=================================================================== 11 38 12 integer, intent(in) :: nz ! number of atmospheric layers 13 integer, intent(in) :: nesp ! number of tracers in traceur.def 14 15 real, dimension(nz) :: p ! pressure (hpa) 16 real, dimension(nz) :: t ! temperature (k) 17 real, dimension(nz) :: zlocal ! altitude (km) 18 real, dimension(nz) :: mumean ! mean molecular mass (g/mol) 19 real :: ptimestep ! physics timestep (s) 20 real :: sza_input ! solar zenith angle (degrees) 21 real :: lon, lat ! longitude and latitude (degrees) 22 23 integer :: n_lon ! for 1D test 39 integer, intent(in) :: nz ! number of atmospheric layers 40 integer, intent(in) :: nesp ! number of tracers in traceur.def 41 integer, intent(in) :: nespeuv ! number of tracers for jthermcal_e107.F 42 43 integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions 44 integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions 45 integer, intent(in) :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions 46 integer, intent(in) :: nphotion ! number of photoionizations 47 48 logical, intent(in) :: ok_ionchem ! switch for ion chemistry 49 logical, intent(in) :: ok_jonline ! switch for on-line calculation of photolysis rates 50 logical, intent(in) :: tuneupperatm! upper atmosphere chemistry tuning 51 52 real, dimension(nz) :: p ! pressure (hpa) 53 real, dimension(nz) :: t ! temperature (k) 54 real, dimension(nz) :: t_elect ! electronic temperature (k) 55 real, dimension(nz) :: zlocal ! altitude (km) 56 real, dimension(nz) :: mumean ! mean molecular mass (g/mol) 57 real, dimension(nz,nespeuv) :: vmr_dens_euv ! tracer mixing ratio for jthermcal_e107 routine 58 real :: ptimestep ! physics timestep (s) 59 real :: sza_input ! solar zenith angle (degrees) 60 real :: lon, lat ! longitude and latitude (degrees) 61 62 integer :: ig ! grid point index 63 64 integer :: n_lon ! for 1D test 24 65 25 66 !=================================================================== … … 30 71 real, dimension(nz,nesp) :: prod_tr ! production (cm-3.s-1) 31 72 real, dimension(nz,nesp) :: loss_tr ! loss (cm-3.s-1) 32 real, dimension (nz) :: em_no! volume emission rate of no33 real, dimension (nz) :: em_o2! volume emission rate of o2(deltag)73 real, dimension (nz) :: em_no ! volume emission rate of no 74 real, dimension (nz) :: em_o2 ! volume emission rate of o2(deltag) 34 75 35 76 !=================================================================== … … 43 84 !=================================================================== 44 85 45 ! jonline86 ! ok_ jonline: see physiq.def 46 87 ! true : on-line calculation of photodissociation rates ! false : lookup table 47 48 logical, save :: jonline = .true.49 88 50 89 logical, save :: firstcall = .true. … … 60 99 61 100 integer, parameter :: nj = 22, nztable = 281, nsza = 27, nso2 = 13 62 real, dimension(nso2,nsza,nztable,nj), save :: jphot 101 !integer, parameter :: nztable = 281, nsza = 27, nso2 = 13 102 real, dimension(nso2,nsza,nztable,nj), save :: jphot ! nj must be equal to nphot 63 103 real, dimension(nztable), save :: table_colair 64 104 real, dimension(nso2,nztable), save :: table_colso2 … … 68 108 ! number densities 69 109 70 real, dimension(nesp) :: cold ! number densities at previous timestep (molecule.cm-3) 71 real, dimension(nz,nesp) :: c ! number densities at current timestep (molecule.cm-3) 72 real, dimension(nesp) :: cnew ! number densities at next timestep (molecule.cm-3) 73 110 real (kind = 8), dimension(nesp) :: cold ! number densities at previous timestep (molecule.cm-3) 111 real (kind = 8), dimension(nz,nesp) :: c ! number densities at current timestep (molecule.cm-3) 112 real (kind = 8), dimension(nz,nespeuv) :: c_euv! number densities for jthermcal_e107 at current timestep (molecule.cm-3) 113 real (kind = 8), dimension(nesp) :: cnew ! number densities at next timestep (molecule.cm-3) 114 74 115 ! timesteps 75 116 … … 80 121 integer :: phychemrat 81 122 123 !Tracer indexes for photionization coeffs 124 125 integer,parameter :: induv_co2 = 1 126 integer,parameter :: induv_o2 = 2 127 integer,parameter :: induv_o = 3 128 integer,parameter :: induv_h2o = 4 129 integer,parameter :: induv_h2 = 5 130 integer,parameter :: induv_h2o2= 6 131 integer,parameter :: induv_o3 = 7 132 integer,parameter :: induv_n2 = 8 133 integer,parameter :: induv_n = 9 134 integer,parameter :: induv_no = 10 135 integer,parameter :: induv_co = 11 136 integer,parameter :: induv_h = 12 137 integer,parameter :: induv_no2 = 13 138 82 139 ! reaction rates 83 140 84 integer, parameter :: nb_phot_max = 3585 integer, parameter :: nb_phot = 2286 integer, parameter :: nb_reaction_3_max = 1287 integer, parameter :: nb_reaction_4_max = 9888 89 real , dimension(nz,nb_phot_max):: v_phot90 real , dimension(nz,nb_reaction_3_max) :: v_391 real , dimension(nz,nb_reaction_4_max) :: v_4141 !integer, parameter :: nb_phot_max = 35 142 !integer, parameter :: nb_phot = 22 143 !integer, parameter :: nb_reaction_3_max = 12 144 !integer, parameter :: nb_reaction_4_max = 98 145 146 real (kind = 8), dimension(nz, nb_phot_max) :: v_phot 147 real (kind = 8), dimension(nz,nb_reaction_3_max) :: v_3 148 real (kind = 8), dimension(nz,nb_reaction_4_max) :: v_4 92 149 93 150 logical,parameter :: hetero_ice = .false. … … 96 153 ! matrix 97 154 98 real , dimension(nesp,nesp) :: mat, mat199 integer, dimension(nesp) :: indx100 integer :: code155 real (kind = 8), dimension(nesp,nesp) :: mat, mat1 156 integer, dimension(nesp) :: indx 157 integer :: code 101 158 102 159 ! production and loss terms (for first-guess solution only) 103 160 104 real , dimension(nesp) :: prod, loss, lossconc161 real (kind = 8), dimension(nesp) :: prod, loss, lossconc 105 162 106 163 ! indexes … … 113 170 !=================================================================== 114 171 115 if (jonline) then 172 ! ok jonline 173 ! true : on-line calculation of photodissociation rates ! false : lookup table 174 if (ok_jonline) then 116 175 print*, 'photochemistry: Read UV absorption cross-sections:' 117 176 call init_photolysis 118 177 else 119 178 print*, 'photochemistry: Read photolysis lookup table:' 120 call init_chimie(n j, nztable, nsza, nso2, jphot, table_colair, table_colso2, table_sza)179 call init_chimie(nphot, nztable, nsza, nso2, jphot, table_colair, table_colso2, table_sza) 121 180 end if 122 181 … … 125 184 !=================================================================== 126 185 127 call indice( nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)186 call indice(ok_ionchem, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max) 128 187 129 188 firstcall = .false. … … 140 199 141 200 do iz = 1,nz 142 conc(iz) = p(iz)/(1.38E-19*t(iz))143 c(iz,:) = tr(iz,:)*conc(iz)201 conc(iz) = p(iz)/(1.38E-19*t(iz)) 202 c(iz,:) = tr(iz,:)*conc(iz) 144 203 end do 145 204 … … 152 211 dist_sol = 0.72333 153 212 154 if ( jonline) then155 if (sza_input <= 95.) then ! day at 30 km213 if (ok_jonline) then 214 if (sza_input <= 95.) then ! day at 300 km 156 215 call photolysis_online(nz, nb_phot_max, zlocal, p, & 157 216 t, mumean, i_co2,i_co, i_o, i_o1d, & … … 161 220 i_no2, i_no, i_n2, i_n2d, & 162 221 nesp, tr, sza_input, dist_sol, v_phot) 222 223 !Calculation of photoionization rates, if needed 224 if (ok_ionchem) then 225 do iz = 1,nz 226 c_euv(iz,:) = vmr_dens_euv(iz,:)*conc(iz) 227 end do 228 !! FAIRE ULTRA ATTENTION LE c_euv et vmr_dens_euv NE SONT PAS DANS LE MEME REGIME 229 !! DE i_ESPECE QUE C, CNEW, COLD et TR 230 call jthermcalc_e107(ig,nz,2,c_euv,nespeuv,t,zlocal,sza_input) 231 do iz=1,nz 232 call phdisrate(ig,nz,2,sza_input,iz) 233 end do 234 !CO2 photoionization 235 v_phot(:,nphot+ 1) = jion(induv_co2,:,1) 236 v_phot(:,nphot+ 2) = jion(induv_co2,:,2) 237 v_phot(:,nphot+ 3) = jion(induv_co2,:,2) 238 v_phot(:,nphot+ 4) = jion(induv_co2,:,3) 239 v_phot(:,nphot+ 5) = jion(induv_co2,:,3) 240 v_phot(:,nphot+ 6) = jion(induv_co2,:,4) 241 v_phot(:,nphot+ 7) = jion(induv_co2,:,4) 242 !O2 photoionization 243 v_phot(:,nphot+ 8) = jion(induv_o2,:,1) 244 !O photoionization 245 v_phot(:,nphot+ 9) = jion(induv_o,:,1) 246 !NO photoionization 247 v_phot(:,nphot+10) = jion(induv_no,:,1) 248 !CO photoionization 249 v_phot(:,nphot+11) = jion(induv_co,:,1) 250 v_phot(:,nphot+12) = jion(induv_co,:,2) 251 v_phot(:,nphot+13) = jion(induv_co,:,2) 252 !N2 photoionization 253 v_phot(:,nphot+14) = jion(induv_n2,:,1) 254 v_phot(:,nphot+15) = jion(induv_n2,:,2) 255 v_phot(:,nphot+16) = jion(induv_n2,:,2) 256 !N photoionization 257 v_phot(:,nphot+17) = jion(induv_n,:,1) 258 !H photoionization 259 v_phot(:,nphot+18) = jion(induv_h,:,1) 260 end if 163 261 else ! night 164 262 v_phot(:,:) = 0. 165 263 end if 166 264 else 167 call phot( nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2), &265 call phot( nj, nztable, nsza, nso2, sza_input, dist_sol, mumean, tr(:,i_co2), tr(:,i_so2), & 168 266 jphot, table_colair, table_colso2, table_sza, nz, nb_phot_max, t, p, v_phot) 169 267 end if … … 173 271 !=================================================================== 174 272 175 call krates(hetero_ice, hetero_dust, nz, nesp, nb_phot, c, conc,t, p, nb_phot_max, nb_reaction_3_max, &176 nb_reaction_4_max, v_3, v_4, v_phot, sza_input, ind_norec, ind_orec)273 call krates(hetero_ice,hetero_dust,ok_ionchem, nphotion, nz, nesp, c, conc, t, t_elect, p, nb_phot_max, nb_reaction_3_max, & 274 nb_reaction_4_max, tuneupperatm, v_3, v_4, v_phot, sza_input, ind_norec, ind_orec) 177 275 178 276 !=================================================================== … … 248 346 c(iz,:) = cnew(:) 249 347 cnew(:) = 0. 348 349 ! force charge neutrality (mod fgg, july 2019) 350 351 if (ok_ionchem) then 352 if(c(iz,i_elec).ne.c(iz,i_co2plus)+c(iz,i_oplus)+c(iz,i_o2plus)+ & 353 c(iz,i_noplus)+c(iz,i_coplus)+c(iz,i_cplus)+c(iz,i_n2plus)+ & 354 c(iz,i_nplus)+c(iz,i_hplus)+c(iz,i_hco2plus)+ & 355 c(iz,i_hcoplus)+c(iz,i_h2oplus)+c(iz,i_h3oplus)+ & 356 c(iz,i_ohplus)) then 357 c(iz,i_elec) = c(iz,i_co2plus)+c(iz,i_oplus)+c(iz,i_o2plus)+ & 358 c(iz,i_noplus)+c(iz,i_coplus)+c(iz,i_cplus)+ & 359 c(iz,i_n2plus)+c(iz,i_nplus)+c(iz,i_hplus)+ & 360 c(iz,i_hco2plus)+c(iz,i_hcoplus)+c(iz,i_h2oplus)+ & 361 c(iz,i_h3oplus)+c(iz,i_ohplus) 362 ! write(*,*)'photochemistry/359' 363 ! write(*,*)'Forcing charge neutrality at ilev,',ilev,' ig=',ig 364 end if 365 end if 250 366 251 367 ! increment internal time … … 323 439 !====================================================================== 324 440 325 subroutine indice( nb_phot_max, nb_reaction_3_max, nb_reaction_4_max)441 subroutine indice(ok_ionchem, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max) 326 442 327 443 !================================================================ … … 344 460 ! input 345 461 346 integer :: nb_phot_max, nb_reaction_3_max, nb_reaction_4_max 462 integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions 463 integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions 464 integer, intent(in) :: nb_phot_max ! number of processes treated numerically as photodissociations 465 logical, intent(in) :: ok_ionchem ! True: Ion reaction 347 466 348 467 ! local … … 537 656 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n, 1.0, i_n2d) 538 657 539 !=========================================================== 540 ! a001 : O + O2 + CO2 -> O3 + CO2 658 !Only if ion chemistry included 659 if (ok_ionchem) then 660 661 !=========================================================== 662 ! CO2 + hv -> CO2+ + e- 663 !=========================================================== 664 665 nb_phot = nb_phot + 1 666 667 indice_phot(nb_phot) = z3spec(1.0, i_co2, 1.0, i_co2plus, 1.0, i_elec) 668 669 !=========================================================== 670 ! CO2 + hv -> O+ + CO + e- 671 !=========================================================== 672 !We divide this reaction in two 673 674 !0.5 CO2 + hv -> CO 675 nb_phot = nb_phot + 1 676 677 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_co, 0.0, i_dummy) 678 679 !0.5 CO2 + hv -> O+ + e- 680 nb_phot = nb_phot + 1 681 682 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_oplus, 1.0, i_elec) 683 684 !=========================================================== 685 ! CO2 + hv -> CO+ + O + e- 686 !=========================================================== 687 !We divide this reaction in two 688 689 !0.5 CO2 + hv -> O 690 nb_phot = nb_phot + 1 691 692 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o, 0.0, i_dummy) 693 694 !0.5 CO2 + hv -> CO+ + e- 695 nb_phot = nb_phot + 1 696 697 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_coplus, 1.0, i_elec) 698 699 !=========================================================== 700 ! CO2 + hv -> C+ + O2 + e- 701 !=========================================================== 702 !We divide this reaction in two 703 704 !0.5 CO2 + hv -> O2 705 nb_phot = nb_phot + 1 706 707 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_o2, 0.0, i_dummy) 708 709 !0.5 CO2 + hv -> C+ + e- 710 nb_phot = nb_phot + 1 711 712 indice_phot(nb_phot) = z3spec(0.5, i_co2, 1.0, i_cplus, 1.0, i_elec) 713 714 !=========================================================== 715 ! O2 + hv -> O2+ + e- 716 !=========================================================== 717 718 nb_phot = nb_phot + 1 719 720 indice_phot(nb_phot) = z3spec(1.0, i_o2, 1.0, i_o2plus, 1.0, i_elec) 721 722 !=========================================================== 723 ! O + hv -> O+ + e- 724 !=========================================================== 725 726 nb_phot = nb_phot + 1 727 728 indice_phot(nb_phot) = z3spec(1.0, i_o, 1.0, i_oplus, 1.0, i_elec) 729 730 !=========================================================== 731 ! NO + hv -> NO+ + e- 732 !=========================================================== 733 734 nb_phot = nb_phot + 1 735 736 indice_phot(nb_phot) = z3spec(1.0, i_no, 1.0, i_noplus, 1.0, i_elec) 737 738 !=========================================================== 739 ! CO + hv -> CO+ + e- 740 !=========================================================== 741 742 nb_phot = nb_phot + 1 743 744 indice_phot(nb_phot) = z3spec(1.0, i_co, 1.0, i_coplus, 1.0, i_elec) 745 746 !=========================================================== 747 ! CO + hv -> C+ + O + e- 748 !=========================================================== 749 !We divide this reaction in two 750 751 !0.5 CO + hv -> O 752 nb_phot = nb_phot + 1 753 754 indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_o, 0.0, i_dummy) 755 756 !0.5 CO + hv -> C+ + e- 757 nb_phot = nb_phot + 1 758 759 indice_phot(nb_phot) = z3spec(0.5, i_co, 1.0, i_cplus, 1.0, i_elec) 760 761 !=========================================================== 762 ! N2 + hv -> N2+ + e- 763 !=========================================================== 764 765 nb_phot = nb_phot + 1 766 767 indice_phot(nb_phot) = z3spec(1.0, i_n2, 1.0, i_n2plus, 1.0, i_elec) 768 769 !=========================================================== 770 ! N2 + hv -> N+ + N + e- 771 !=========================================================== 772 !We divide this reaction in two 773 774 !0.5 N2 + hv -> N 775 nb_phot = nb_phot + 1 776 777 indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_n, 0.0, i_dummy) 778 779 !0.5 N2 + hv -> N+ + e- 780 nb_phot = nb_phot + 1 781 782 indice_phot(nb_phot) = z3spec(0.5, i_n2, 1.0, i_nplus, 1.0, i_elec) 783 784 !=========================================================== 785 ! N + hv -> N+ + e- 786 !=========================================================== 787 788 nb_phot = nb_phot + 1 789 790 indice_phot(nb_phot) = z3spec(1.0, i_n, 1.0, i_nplus, 1.0, i_elec) 791 792 !=========================================================== 793 ! H + hv -> H+ + e- 794 !=========================================================== 795 796 nb_phot = nb_phot + 1 797 798 indice_phot(nb_phot) = z3spec(1.0, i_h, 1.0, i_hplus, 1.0, i_elec) 799 800 end if !ok_ionchem 801 802 !=========================================================== 803 ! a001 : O + O2 + (CO2 or M) -> O3 + (CO2 or M) 541 804 !=========================================================== 542 805 … … 546 809 547 810 !=========================================================== 548 ! a002 : O + O + CO2 -> O2(Dg) + CO2811 ! a002 : O + O + (CO2 or M) -> O2(Dg) + (CO2 or M) 549 812 !=========================================================== 550 813 … … 1502 1765 indice_phot(nb_phot) = z3spec(1.0, i_h2o2, 1.0, i_h2o, 0.5, i_o2) 1503 1766 1504 !=========================================================== 1505 ! i001: O2(Dg) + CO2 -> O2 + CO2 1767 !Only if ion chemistry 1768 if (ok_ionchem) then 1769 1770 !=========================================================== 1771 ! i001 : CO2+ + O2 -> O2+ + CO2 1772 !=========================================================== 1773 1774 nb_reaction_4 = nb_reaction_4 + 1 1775 1776 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_co2) 1777 1778 !=========================================================== 1779 ! i002 : CO2+ + O -> O+ + CO2 1780 !=========================================================== 1781 1782 nb_reaction_4 = nb_reaction_4 + 1 1783 1784 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co2) 1785 1786 !=========================================================== 1787 ! i003 : CO2+ + O -> O2+ + CO 1788 !=========================================================== 1789 1790 nb_reaction_4 = nb_reaction_4 + 1 1791 1792 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_co) 1793 1794 !=========================================================== 1795 ! i004 : O2+ + e- -> O + O 1796 !=========================================================== 1797 1798 nb_reaction_4 = nb_reaction_4 + 1 1799 1800 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_elec, 2.0, i_o, 0.0, i_dummy) 1801 1802 !=========================================================== 1803 ! i005 : O+ + CO2 -> O2+ + CO 1804 !=========================================================== 1805 1806 nb_reaction_4 = nb_reaction_4 + 1 1807 1808 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_co2, 1.0, i_o2plus, 1.0, i_co) 1809 1810 !=========================================================== 1811 ! i006 : CO2+ + e -> CO + O 1812 !=========================================================== 1813 1814 nb_reaction_4 = nb_reaction_4 + 1 1815 1816 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_elec, 1.0, i_co, 1.0, i_o) 1817 1818 !=========================================================== 1819 ! i007 : CO2+ + NO -> NO+ + CO2 1820 !=========================================================== 1821 1822 nb_reaction_4 = nb_reaction_4 + 1 1823 1824 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_co2) 1825 1826 !=========================================================== 1827 ! i008 : O2+ + NO -> NO+ + O2 1828 !=========================================================== 1829 1830 nb_reaction_4 = nb_reaction_4 + 1 1831 1832 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_no, 1.0, i_noplus, 1.0, i_o2) 1833 1834 !=========================================================== 1835 ! i009 : O2+ + N2 -> NO+ + NO 1836 !=========================================================== 1837 1838 nb_reaction_4 = nb_reaction_4 + 1 1839 1840 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_no) 1841 1842 !=========================================================== 1843 ! i010 : O2+ + N -> NO+ + O 1844 !=========================================================== 1845 1846 nb_reaction_4 = nb_reaction_4 + 1 1847 1848 indice_4(nb_reaction_4) = z4spec(1.0, i_o2plus, 1.0, i_n, 1.0, i_noplus, 1.0, i_o) 1849 1850 !=========================================================== 1851 ! i011 : O+ + N2 -> NO+ + N 1852 !=========================================================== 1853 1854 nb_reaction_4 = nb_reaction_4 + 1 1855 1856 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_n2, 1.0, i_noplus, 1.0, i_n) 1857 1858 !=========================================================== 1859 ! i012 : NO+ + e -> N + O 1860 !=========================================================== 1861 1862 nb_reaction_4 = nb_reaction_4 + 1 1863 1864 indice_4(nb_reaction_4) = z4spec(1.0, i_noplus, 1.0, i_elec, 1.0, i_n, 1.0, i_o) 1865 1866 !=========================================================== 1867 ! i013 : CO+ + CO2 -> CO2+ + CO 1868 !=========================================================== 1869 1870 nb_reaction_4 = nb_reaction_4 + 1 1871 1872 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_co) 1873 1874 !=========================================================== 1875 ! i014 : CO+ + O -> O+ + CO 1876 !=========================================================== 1877 1878 nb_reaction_4 = nb_reaction_4 + 1 1879 1880 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_co) 1881 1882 !=========================================================== 1883 ! i015 : C+ + CO2 -> CO+ + CO 1884 !=========================================================== 1885 1886 nb_reaction_4 = nb_reaction_4 + 1 1887 1888 indice_4(nb_reaction_4) = z4spec(1.0, i_cplus, 1.0, i_co2, 1.0, i_coplus, 1.0, i_co) 1889 1890 !=========================================================== 1891 ! i016 : N2+ + CO2 -> CO2+ + N2 1892 !=========================================================== 1893 1894 nb_reaction_4 = nb_reaction_4 + 1 1895 1896 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n2) 1897 1898 !=========================================================== 1899 ! i017 : N2+ + O -> NO+ + N 1900 !=========================================================== 1901 1902 nb_reaction_4 = nb_reaction_4 + 1 1903 1904 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_noplus, 1.0, i_n) 1905 1906 !=========================================================== 1907 ! i018 : N2+ + CO -> CO+ + N2 1908 !=========================================================== 1909 1910 nb_reaction_4 = nb_reaction_4 + 1 1911 1912 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_co, 1.0, i_coplus, 1.0, i_n2) 1913 1914 !=========================================================== 1915 ! i019 : N2+ + e -> N + N 1916 !=========================================================== 1917 1918 nb_reaction_4 = nb_reaction_4 + 1 1919 1920 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_elec, 2.0, i_n, 0.0, i_dummy) 1921 1922 !=========================================================== 1923 ! i020 : N2+ + O -> O+ + N2 1924 !=========================================================== 1925 1926 nb_reaction_4 = nb_reaction_4 + 1 1927 1928 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_o, 1.0, i_oplus, 1.0, i_n2) 1929 1930 !=========================================================== 1931 ! i021 : N+ + CO2 -> CO2+ + N 1932 !=========================================================== 1933 1934 nb_reaction_4 = nb_reaction_4 + 1 1935 1936 indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_co2, 1.0, i_co2plus, 1.0, i_n) 1937 1938 !=========================================================== 1939 ! i022 : CO+ + H -> H+ + CO 1940 !=========================================================== 1941 1942 nb_reaction_4 = nb_reaction_4 + 1 1943 1944 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_co) 1945 1946 !=========================================================== 1947 ! i023 : O+ + H -> H+ + O 1948 !=========================================================== 1949 1950 nb_reaction_4 = nb_reaction_4 + 1 1951 1952 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h, 1.0, i_hplus, 1.0, i_o) 1953 1954 !=========================================================== 1955 ! i024 : H+ + O -> O+ + H 1956 !=========================================================== 1957 1958 nb_reaction_4 = nb_reaction_4 + 1 1959 1960 indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_o, 1.0, i_oplus, 1.0, i_h) 1961 1962 !=========================================================== 1963 ! i025 : CO2+ + H2 -> HCO2+ + H 1964 !=========================================================== 1965 1966 nb_reaction_4 = nb_reaction_4 + 1 1967 1968 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2, 1.0, i_hco2plus, 1.0, i_h) 1969 1970 !=========================================================== 1971 ! i026 : HCO2+ + e -> H + CO2 1972 !=========================================================== 1973 1974 nb_reaction_4 = nb_reaction_4 + 1 1975 1976 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2) 1977 1978 !=========================================================== 1979 ! i027 : HCO2+ + e -> H + O + CO 1980 !=========================================================== 1981 !We divide this reaction in two 1982 1983 !0.5HCO2+ + 0.5e -> H 1984 1985 nb_reaction_4 = nb_reaction_4 + 1 1986 1987 indice_4(nb_reaction_4) = z4spec(.5, i_hco2plus, 0.5, i_elec, 1.0, i_h, 0.0, i_dummy) 1988 1989 !0.5 HCO2+ + 0.5 e -> O + CO 1990 1991 nb_reaction_4 = nb_reaction_4 + 1 1992 1993 indice_4(nb_reaction_4) = z4spec(0.5, i_hco2plus, 0.5, i_elec, 1.0, i_o, 1.0, i_co) 1994 1995 !=========================================================== 1996 ! i029 : HCO2+ + e -> OH + CO 1997 !=========================================================== 1998 1999 nb_reaction_4 = nb_reaction_4 + 1 2000 2001 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_oh, 1.0, i_co) 2002 2003 2004 !=========================================================== 2005 ! i030 : HCO2+ + e -> H + CO2 2006 !=========================================================== 2007 2008 nb_reaction_4 = nb_reaction_4 + 1 2009 2010 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_elec, 1.0, i_h, 1.0, i_co2) 2011 2012 2013 !=========================================================== 2014 ! i031 : HCO2+ + O -> HCO+ + O2 2015 !=========================================================== 2016 2017 nb_reaction_4 = nb_reaction_4 + 1 2018 2019 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_o, 1.0, i_hcoplus, 1.0, i_o2) 2020 2021 2022 !=========================================================== 2023 ! i032 : HCO2+ + CO -> HCO+ + CO2 2024 !=========================================================== 2025 2026 nb_reaction_4 = nb_reaction_4 + 1 2027 indice_4(nb_reaction_4) = z4spec(1.0, i_hco2plus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_co2) 2028 2029 2030 !=========================================================== 2031 ! i033 : H+ + CO2 -> HCO+ + O 2032 !=========================================================== 2033 2034 nb_reaction_4 = nb_reaction_4 + 1 2035 indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_co2, 1.0, i_hcoplus, 1.0, i_o) 2036 2037 2038 !=========================================================== 2039 ! i034 : CO2+ + H -> HCO+ + O 2040 !=========================================================== 2041 2042 nb_reaction_4 = nb_reaction_4 + 1 2043 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h, 1.0, i_hcoplus, 1.0, i_o) 2044 2045 2046 !=========================================================== 2047 ! i035 : CO+ + H2 -> HCO+ + H 2048 !=========================================================== 2049 2050 nb_reaction_4 = nb_reaction_4 + 1 2051 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2, 1.0, i_hcoplus, 1.0, i_h) 2052 2053 2054 !=========================================================== 2055 ! i036 : HCO+ + e- -> CO + H 2056 !=========================================================== 2057 2058 nb_reaction_4 = nb_reaction_4 + 1 2059 indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_elec, 1.0, i_co, 1.0, i_h) 2060 2061 !=========================================================== 2062 ! i037 : CO2+ + H2O -> H2O+ + CO2 2063 !=========================================================== 2064 2065 nb_reaction_4 = nb_reaction_4 + 1 2066 indice_4(nb_reaction_4) = z4spec(1.0, i_co2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co2) 2067 2068 !=========================================================== 2069 ! i038 : CO+ + H2O -> H2O+ + CO 2070 !=========================================================== 2071 2072 nb_reaction_4 = nb_reaction_4 + 1 2073 indice_4(nb_reaction_4) = z4spec(1.0, i_coplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_co) 2074 2075 !=========================================================== 2076 ! i039 : O+ + H2O -> H2O+ + O 2077 !=========================================================== 2078 2079 nb_reaction_4 = nb_reaction_4 + 1 2080 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_o) 2081 2082 !=========================================================== 2083 ! i040 : N2+ + H2O -> H2O+ + N2 2084 !=========================================================== 2085 2086 nb_reaction_4 = nb_reaction_4 + 1 2087 indice_4(nb_reaction_4) = z4spec(1.0, i_n2plus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n2) 2088 2089 !=========================================================== 2090 ! i041 : N+ + H2O -> H2O+ + N 2091 !=========================================================== 2092 2093 nb_reaction_4 = nb_reaction_4 + 1 2094 indice_4(nb_reaction_4) = z4spec(1.0, i_nplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_n) 2095 2096 !=========================================================== 2097 ! i042 : H+ + H2O -> H2O+ + H 2098 !=========================================================== 2099 2100 nb_reaction_4 = nb_reaction_4 + 1 2101 indice_4(nb_reaction_4) = z4spec(1.0, i_hplus, 1.0, i_h2o, 1.0, i_h2oplus, 1.0, i_h) 2102 2103 !=========================================================== 2104 ! i043 : H2O+ + O2 -> O2+ + H2O 2105 !=========================================================== 2106 2107 nb_reaction_4 = nb_reaction_4 + 1 2108 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_h2o) 2109 2110 !=========================================================== 2111 ! i044 : H2O+ + CO -> HCO+ + OH 2112 !=========================================================== 2113 2114 nb_reaction_4 = nb_reaction_4 + 1 2115 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_oh) 2116 2117 !=========================================================== 2118 ! i045 : H2O+ + O -> O2+ + H2 2119 !=========================================================== 2120 2121 nb_reaction_4 = nb_reaction_4 + 1 2122 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h2) 2123 2124 !=========================================================== 2125 ! i046 : H2O+ + NO -> NO+ + H2O 2126 !=========================================================== 2127 2128 nb_reaction_4 = nb_reaction_4 + 1 2129 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_h2o) 2130 2131 !=========================================================== 2132 ! i047 : H2O+ + e- -> H + H + O 2133 !=========================================================== 2134 2135 nb_reaction_4 = nb_reaction_4 + 1 2136 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 2.0, i_h, 1.0, i_o) 2137 2138 !=========================================================== 2139 ! i048 : H2O+ + e- -> H + OH 2140 !=========================================================== 2141 2142 nb_reaction_4 = nb_reaction_4 + 1 2143 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h, 1.0, i_oh) 2144 2145 !=========================================================== 2146 ! i049 : H2O+ + e- -> H2 + O 2147 !=========================================================== 2148 2149 nb_reaction_4 = nb_reaction_4 + 1 2150 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_elec, 1.0, i_h2, 1.0, i_o) 2151 2152 !=========================================================== 2153 ! i050 : H2O+ + H2O -> H3O+ + OH 2154 !=========================================================== 2155 2156 nb_reaction_4 = nb_reaction_4 + 1 2157 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_oh) 2158 2159 !=========================================================== 2160 ! i051 : H2O+ + H2 -> H3O+ + H 2161 !=========================================================== 2162 2163 nb_reaction_4 = nb_reaction_4 + 1 2164 indice_4(nb_reaction_4) = z4spec(1.0, i_h2oplus, 1.0, i_h2, 1.0, i_h3oplus, 1.0, i_h) 2165 2166 !=========================================================== 2167 ! i052 : HCO+ + H2O -> H3O+ + CO 2168 !=========================================================== 2169 2170 nb_reaction_4 = nb_reaction_4 + 1 2171 indice_4(nb_reaction_4) = z4spec(1.0, i_hcoplus, 1.0, i_h2o, 1.0, i_h3oplus, 1.0, i_co) 2172 2173 !=========================================================== 2174 ! i053: H3O+ + e -> OH + H + H 2175 !=========================================================== 2176 2177 nb_reaction_4 = nb_reaction_4 + 1 2178 indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 2.0, i_h) 2179 2180 !=========================================================== 2181 ! i054: H3O+ + e -> H2O + H 2182 !=========================================================== 2183 2184 nb_reaction_4 = nb_reaction_4 + 1 2185 indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_h2o, 1.0, i_h) 2186 2187 !=========================================================== 2188 ! i055: H3O+ + e -> HO + H2 2189 !=========================================================== 2190 2191 nb_reaction_4 = nb_reaction_4 + 1 2192 indice_4(nb_reaction_4) = z4spec(1.0, i_h3oplus, 1.0, i_elec, 1.0, i_oh, 1.0, i_h2) 2193 2194 !=========================================================== 2195 ! i056: H3O+ + e -> O + H2 + H 2196 !=========================================================== 2197 !We divide this reaction in two 2198 2199 !0.5H3O+ + 0.5e -> O 2200 2201 nb_reaction_4 = nb_reaction_4 + 1 2202 indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_o, 0.0, i_dummy) 2203 2204 !0.5H3O+ + 0.5e -> H2 + H 2205 2206 nb_reaction_4 = nb_reaction_4 + 1 2207 indice_4(nb_reaction_4) = z4spec(0.5, i_h3oplus, 0.5, i_elec, 1.0, i_h2, 1.0, i_h) 2208 2209 !=========================================================== 2210 ! i057: O+ + H2 -> OH+ + H 2211 !=========================================================== 2212 2213 nb_reaction_4 = nb_reaction_4 + 1 2214 indice_4(nb_reaction_4) = z4spec(1.0, i_oplus, 1.0, i_h2, 1.0, i_ohplus, 1.0, i_h) 2215 2216 !=========================================================== 2217 ! i058: OH+ + O -> O2+ + H 2218 !=========================================================== 2219 2220 nb_reaction_4 = nb_reaction_4 + 1 2221 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o, 1.0, i_o2plus, 1.0, i_h) 2222 2223 !=========================================================== 2224 ! i059: OH+ + CO2 -> HCO2+ + O 2225 !=========================================================== 2226 2227 nb_reaction_4 = nb_reaction_4 + 1 2228 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co2, 1.0, i_hco2plus, 1.0, i_o) 2229 2230 !=========================================================== 2231 ! i060: OH+ + CO -> HCO+ + O 2232 !=========================================================== 2233 2234 nb_reaction_4 = nb_reaction_4 + 1 2235 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_co, 1.0, i_hcoplus, 1.0, i_o) 2236 2237 !=========================================================== 2238 ! i061: OH+ + NO -> NO+ + OH 2239 !=========================================================== 2240 2241 nb_reaction_4 = nb_reaction_4 + 1 2242 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_no, 1.0, i_noplus, 1.0, i_oh) 2243 2244 !=========================================================== 2245 ! i062: OH+ + H2 -> H2O+ + H 2246 !=========================================================== 2247 2248 nb_reaction_4 = nb_reaction_4 + 1 2249 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_h2, 1.0, i_h2oplus, 1.0, i_h) 2250 2251 !=========================================================== 2252 ! i063: OH+ + O2 -> O2+ + OH 2253 !=========================================================== 2254 2255 nb_reaction_4 = nb_reaction_4 + 1 2256 indice_4(nb_reaction_4) = z4spec(1.0, i_ohplus, 1.0, i_o2, 1.0, i_o2plus, 1.0, i_oh) 2257 2258 end if !ok_ionchem 2259 2260 !=========================================================== 2261 ! j001: O2(Dg) + (CO2 or O) -> O2 + (CO2 or O) 1506 2262 !=========================================================== 1507 2263 … … 1511 2267 1512 2268 !=========================================================== 1513 ! i002: O2(Dg) -> O2 + hv2269 ! j002: O2(Dg) -> O2 + hv 1514 2270 !=========================================================== 1515 2271 … … 1732 2488 !-- TuneE 1733 2489 ! VCD 2.0 tuning 1734 v_phot(65:90,4) = v_phot(65:90,4)*10. 2490 v_phot(65:90,4) = v_phot(65:90,4)*10. ! CO2 + hv ==> O(1D) + CO 1735 2491 !-- 1736 2492 ! v_phot(:,4) = v_phot(:,4)*10. … … 1753 2509 !====================================================================== 1754 2510 1755 subroutine krates(hetero_ice,hetero_dust, nz, nesp, nj, c, conc, t, p, nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, v_3, v_4, v_phot,sza_input, ind_norec, ind_orec) 2511 subroutine krates(hetero_ice,hetero_dust,ok_ionchem, nphotion, & 2512 nz, nesp, c, conc, t, t_elect, p, & 2513 nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, & 2514 tuneupperatm, v_3, v_4, v_phot,sza_input, & 2515 ind_norec, ind_orec) 1756 2516 1757 2517 !================================================================ … … 1768 2528 1769 2529 USE chemparam_mod 2530 USE photolysis_mod, only : nphot 1770 2531 implicit none 1771 2532 1772 real, INTENT(IN), dimension(nz) :: t, p, conc 1773 integer, INTENT(IN) :: nesp, nj, nz 1774 real, INTENT(IN) :: sza_input 2533 !---------------------------------------------------------------------- 2534 ! input 2535 !---------------------------------------------------------------------- 2536 2537 integer, INTENT(IN) :: nesp, nz 2538 real, INTENT(IN) :: sza_input ! [degree] 2539 real, INTENT(IN), dimension(nz) :: t_elect, t, p, conc 2540 logical, INTENT(IN) :: hetero_ice, hetero_dust, tuneupperatm 2541 2542 real, INTENT(IN), dimension(nz,nesp) :: c 2543 2544 integer, intent(in) :: nb_reaction_3_max ! number of quadratic reactions 2545 integer, intent(in) :: nb_reaction_4_max ! number of bimolecular reactions 2546 integer, intent(in) :: nb_phot_max ! number of reactions treated numerically as photodissociations 2547 integer, intent(in) :: nphotion ! number of photoionizations 2548 logical, intent(in) :: ok_ionchem ! if .true. then ionchem reaction used 2549 2550 !---------------------------------------------------------------------- 2551 ! output 2552 !---------------------------------------------------------------------- 2553 2554 real (kind = 8), dimension(nz, nb_phot_max) :: v_phot 2555 real (kind = 8), dimension(nz, nb_reaction_3_max) :: v_3 2556 real (kind = 8), dimension(nz, nb_reaction_4_max) :: v_4 2557 integer :: ind_norec 2558 integer :: ind_orec 2559 2560 !---------------------------------------------------------------------- 2561 ! local 2562 !---------------------------------------------------------------------- 2563 1775 2564 integer :: iz 1776 logical, INTENT(IN) :: hetero_ice, hetero_dust1777 integer :: ind_norec, ind_orec1778 2565 real :: ak0, ak1, xpo, rate, rate1, rate2, pi, gam 1779 2566 real :: k1a0, k1b0, k1ainf, k1a, k1b, fc, fx, x, y 2567 integer :: nb_phot, nb_reaction_3, nb_reaction_4 1780 2568 real, dimension(nz) :: surfice1d, surfdust1d 1781 2569 real, dimension(nz) :: deq … … 1807 2595 g029, g030, g031, g032, g033, & 1808 2596 h001, h002, h003, & 1809 i001, i002 1810 1811 real, INTENT(IN), dimension(nz,nesp) :: c 1812 1813 integer :: nb_phot_max, nb_reaction_3_max, nb_reaction_4_max, nb_reaction_3, nb_reaction_4, nb_phot 1814 1815 real, dimension(nz,nb_phot_max) :: v_phot 1816 real, dimension(nz,nb_reaction_3_max) :: v_3 1817 real, dimension(nz,nb_reaction_4_max) :: v_4 1818 1819 pi = acos(-1.) 1820 1821 nb_phot = nj 1822 nb_reaction_3 = 0 1823 nb_reaction_4 = 0 2597 i001, i002, i003, i004, i005, i006, & 2598 i007, i008, i009, i010, i011, i012, & 2599 i013, i014, i015, i016, i017, i018, i019, & 2600 i020, i021, i022, i023, i024, i025, i026, & 2601 i027, i028, i029, i030, i031, i032, i033, & 2602 i034, i035, i036, i037, i038, i039, i040, & 2603 i041, i042, i043, i044, i045, i046, i047, & 2604 i048, i049, i050, i051, i052, i053, i054, & 2605 i055, i056, i057, i058, i059, i060, i061, & 2606 i062, i063, & 2607 j001, j002 2608 2609 !---------------------------------------------------------------------- 2610 ! initialisation 2611 !---------------------------------------------------------------------- 2612 2613 pi = acos(-1.) 2614 2615 nb_phot = nphot + nphotion ! initialised to the number of photolysis + number of photoionization rates 2616 nb_reaction_3 = 0 2617 nb_reaction_4 = 0 1824 2618 1825 2619 !---------------------------------------------------------------------- … … 1827 2621 !---------------------------------------------------------------------- 1828 2622 1829 !--- a001: o + o2 + co2 -> o3 + co22623 !--- a001: o + o2 + (co2 or M) -> o3 + (co2 or M) 1830 2624 1831 2625 ! jpl 2003 1832 1833 a001(:) = 2.5*6.0E-34*(t(:)/300.)**(-2.4)*conc(:) 2626 ! ------ BEFORE VCD 2.1 ------- 2627 ! a001(:) = 2.5*6.0E-34*(t(:)/300.)**(-2.4)*conc(:) 2628 ! ----------------------------- 2629 2630 ! ! nist expression + take account of the N2 and CO2 factor 2631 ! ------ VCD 2.1 update ------- 2632 a001(:) = 6.0E-34*(t(:)/300.)**(-2.4) * & 2633 ( 2.5 * c(:,i_co2) + & 2634 1.0 * (conc(:)-c(:,i_co2)) ) 2635 ! ----------------------------- 1834 2636 1835 2637 nb_reaction_4 = nb_reaction_4 + 1 1836 2638 v_4(:,nb_reaction_4) = a001(:) 1837 2639 1838 !--- a002: o + o + co2 -> o2 + co22640 !--- a002: o + o + (co2 or M) -> o2(Dg) + (co2 or M) 1839 2641 1840 2642 ! Tsang and Hampson, J. Chem. Phys. Ref. Data, 15, 1087, 1986 … … 1843 2645 1844 2646 ! Campbell and Gray, Chem. Phys. Lett., 18, 607, 1973 1845 1846 a002(:) = 2.5*9.46E-34*exp(485./t(:))*conc(:) ! nist expression 1847 2647 ! ------ BEFORE VCD 2.1 ------- 2648 ! a002(:) = 2.5*9.46E-34*exp(485./t(:))*conc(:) ! nist expression 2649 ! ----------------------------- 2650 2651 ! ! nist expression + take account of the N2 and CO2 factor 2652 ! a002(:) = 9.46E-34*exp(485./t(:)) * & 2653 ! ( 2.5 * c(:,i_co2) + & 2654 ! 1.0 * (conc(:)-c(:,i_co2)) ) 2655 2656 ! Baulch et al., 1976 (confirmed by smith and robertson, 2008) 2657 ! ------ VCD 2.1 update ------- 2658 2659 a002(:) = 2.76E-34*exp(720./t(:)) * & 2660 ( 2.5 * c(:,i_co2) + & 2661 1.0 * (conc(:)-c(:,i_co2)) ) 2662 ! ----------------------------- 2663 1848 2664 nb_reaction_3 = nb_reaction_3 + 1 1849 2665 v_3(:,nb_reaction_3) = a002(:) … … 3126 3942 3127 3943 !---------------------------------------------------------------------- 3944 ! ionospheric reactions 3945 ! only if ok_ionchem=true 3946 !---------------------------------------------------------------------- 3947 3948 if (ok_ionchem) then 3949 3950 !--- i001: co2+ + o2 -> o2+ + co2 3951 3952 ! aninich, j. phys. chem. ref. data 1993 3953 3954 i001(:) = 5.5e-11*(300./t_elect(:))**0.5 3955 3956 nb_reaction_4 = nb_reaction_4 + 1 3957 v_4(:,nb_reaction_4) = i001(:) 3958 3959 !--- i002: co2+ + o -> o+ + co2 3960 3961 ! UMIST database 3962 3963 i002(:) = 9.6e-11 3964 3965 nb_reaction_4 = nb_reaction_4 + 1 3966 v_4(:,nb_reaction_4) = i002(:) 3967 3968 !--- i003: co2+ + o -> o2+ + co 3969 3970 ! UMIST database 3971 3972 i003(:) = 1.64e-10 3973 3974 nb_reaction_4 = nb_reaction_4 + 1 3975 v_4(:,nb_reaction_4) = i003(:) 3976 3977 !--- i004: o2+ + e- -> o + o 3978 3979 ! Alge et al., J. Phys. B, At. Mol. Phys. 1983 3980 !i004(:) = 2.0e-7*(300./t_elect(:))**0.7 3981 3982 do iz = 1,nz 3983 if (t_elect(iz)<1200.) then 3984 i004(iz) = 2.0e-7*(300./t_elect(iz))**0.7 3985 else 3986 i004(iz) = 7.4e-8*(1200./t_elect(iz))**0.56 3987 end if 3988 end do 3989 3990 nb_reaction_4 = nb_reaction_4 + 1 3991 v_4(:,nb_reaction_4) = i004(:) 3992 3993 !--- i005: o+ + co2 -> o2+ + co 3994 3995 ! UMIST database 3996 3997 i005(:) = 9.4e-10 3998 3999 nb_reaction_4 = nb_reaction_4 + 1 4000 v_4(:,nb_reaction_4) = i005(:) 4001 4002 4003 !--- i006: co2+ + e- -> co + o 4004 4005 ! UMIST database 4006 4007 i006(:) = 3.8e-7*(300./t_elect(:))**0.5 4008 4009 nb_reaction_4 = nb_reaction_4 + 1 4010 v_4(:,nb_reaction_4) = i006(:) 4011 4012 4013 !--- i007: co2+ + no -> no+ + co2 4014 4015 ! UMIST database 4016 4017 i007(:) = 1.2e-10 4018 4019 nb_reaction_4 = nb_reaction_4 + 1 4020 v_4(:,nb_reaction_4) = i007(:) 4021 4022 !--- i008: o2+ + no -> no+ + o2 4023 4024 ! UMIST database 4025 4026 i008(:) = 4.6e-10 4027 4028 nb_reaction_4 = nb_reaction_4 + 1 4029 v_4(:,nb_reaction_4) = i008(:) 4030 4031 !--- i009: o2+ + n2 -> no+ + no 4032 4033 ! Fox & Sung 2001 4034 4035 i009(:) = 1.0e-15 4036 4037 nb_reaction_4 = nb_reaction_4 + 1 4038 v_4(:,nb_reaction_4) = i009(:) 4039 4040 !--- i010: o2+ + n -> no+ + o 4041 4042 ! Fox & Sung 2001 4043 4044 i010(:) = 1.0e-10 4045 4046 nb_reaction_4 = nb_reaction_4 + 1 4047 v_4(:,nb_reaction_4) = i010(:) 4048 4049 !--- i011: o+ + n2 -> no+ + n 4050 4051 ! Fox & Sung 2001 4052 4053 i011(:) = 1.2e-12 * (300./t_elect(:))**0.45 4054 4055 nb_reaction_4 = nb_reaction_4 + 1 4056 v_4(:,nb_reaction_4) = i011(:) 4057 4058 !--- i012: no+ + e -> n + o 4059 4060 ! UMIST database 4061 4062 i012(:) = 4.3e-7*(300./t_elect(:))**0.37 4063 4064 nb_reaction_4 = nb_reaction_4 + 1 4065 v_4(:,nb_reaction_4) = i012(:) 4066 4067 4068 !--- i013: co+ + co2 -> co2+ + co 4069 4070 ! UMIST database 4071 4072 i013(:) = 1.0e-9 4073 4074 nb_reaction_4 = nb_reaction_4 + 1 4075 v_4(:,nb_reaction_4) = i013(:) 4076 4077 4078 !--- i014: co+ + o -> o+ + co 4079 4080 ! UMIST database 4081 4082 i014(:) = 1.4e-10 4083 4084 nb_reaction_4 = nb_reaction_4 + 1 4085 v_4(:,nb_reaction_4) = i014(:) 4086 4087 !--- i015: c+ + co2 -> co+ + co 4088 4089 ! UMIST database 4090 4091 i015(:) = 1.1e-9 4092 4093 nb_reaction_4 = nb_reaction_4 + 1 4094 v_4(:,nb_reaction_4) = i015(:) 4095 4096 4097 !--- i016: N2+ + co2 -> co2+ + N2 4098 4099 ! Fox & Song 2001 4100 4101 i016(:) = 9.0e-10*(300./t_elect(:))**0.23 4102 4103 nb_reaction_4 = nb_reaction_4 + 1 4104 v_4(:,nb_reaction_4) = i016(:) 4105 4106 4107 !--- i017: N2+ + o -> no+ + N 4108 4109 ! Fox & Song 2001 4110 4111 i017(:) = 1.33e-10*(300./t_elect(:))**0.44 4112 4113 nb_reaction_4 = nb_reaction_4 + 1 4114 v_4(:,nb_reaction_4) = i017(:) 4115 4116 !--- i018: N2+ + co -> co+ + N2 4117 4118 ! UMIST 4119 4120 i018(:) = 7.4e-11 4121 4122 nb_reaction_4 = nb_reaction_4 + 1 4123 v_4(:,nb_reaction_4) = i018(:) 4124 4125 !--- i019: N2+ + e -> N + N 4126 4127 ! UMIST 4128 4129 i019(:) = 7.7e-7*(300./t_elect(:))**0.3 4130 4131 nb_reaction_4 = nb_reaction_4 + 1 4132 v_4(:,nb_reaction_4) = i016(:) 4133 4134 !--- i020: N2+ + o -> o+ + N2 4135 4136 ! Fox & Song 2001 4137 4138 i020(:) = 7.0e-12*(300./t_elect(:))**0.23 4139 4140 nb_reaction_4 = nb_reaction_4 + 1 4141 v_4(:,nb_reaction_4) = i020(:) 4142 4143 !--- i021: N+ + co2 -> co2+ + N 4144 4145 ! UMIST 4146 4147 i021(:) = 7.5e-10 4148 4149 nb_reaction_4 = nb_reaction_4 + 1 4150 v_4(:,nb_reaction_4) = i021(:) 4151 4152 !--- i022: CO+ + H -> H+ + CO 4153 4154 ! Fox & Sung 2001 4155 4156 i022(:) = 4.0e-10 4157 4158 nb_reaction_4 = nb_reaction_4 + 1 4159 v_4(:,nb_reaction_4) = i022(:) 4160 4161 !--- i023: O+ + H -> H+ + O 4162 4163 ! UMIST 4164 4165 i023(:) = 5.66e-10*((t_elect(:)/300.)**0.36)*exp(8.6/t_elect(:)) 4166 4167 nb_reaction_4 = nb_reaction_4 + 1 4168 v_4(:,nb_reaction_4) = i023(:) 4169 4170 !--- i024: H+ + O -> O+ + H 4171 4172 ! UMIST 4173 4174 i024(:) = 6.86e-10*((t_elect(:)/300.)**0.26)*exp(-224.3/t_elect(:)) 4175 4176 nb_reaction_4 = nb_reaction_4 + 1 4177 v_4(:,nb_reaction_4) = i024(:) 4178 4179 !--- i025: CO+ + H2 -> HCO2+ + H 4180 4181 ! UMIST 4182 4183 i025(:) = 9.5e-10 4184 4185 nb_reaction_4 = nb_reaction_4 + 1 4186 v_4(:,nb_reaction_4) = i025(:) 4187 4188 !--- i026: HCO2+ + e -> H + CO2 4189 4190 ! UMIST 4191 4192 i026(:) = 1.75e-8*((300./t_elect(:))**0.5) 4193 4194 nb_reaction_4 = nb_reaction_4 + 1 4195 v_4(:,nb_reaction_4) = i026(:) 4196 4197 !--- i027+i028: HCO2+ + e -> H + O + CO 4198 4199 ! UMIST 4200 !Reaction splitted in 2: i027: 0.5 (HCO2+ + e-) -> H 4201 !i028: 0.5 (HCO2+ + e-) -> O + CO 4202 4203 i027(:) = 8.1e-7*((300./t_elect(:))**0.64) 4204 4205 nb_reaction_4 = nb_reaction_4 + 1 4206 v_4(:,nb_reaction_4) = i027(:) 4207 4208 nb_reaction_4 = nb_reaction_4 + 1 4209 v_4(:,nb_reaction_4) = i027(:) 4210 4211 !--- i029: HCO2+ + e -> OH + CO 4212 4213 ! UMIST 4214 4215 i029(:) = 3.2e-7*((300./t_elect(:))**0.64) 4216 4217 nb_reaction_4 = nb_reaction_4 + 1 4218 v_4(:,nb_reaction_4) = i029(:) 4219 4220 !--- i030: HCO2+ + e -> H + CO2 4221 4222 i030(:) = 6.0e-8*((300./t_elect(:))**0.64) 4223 nb_reaction_4 = nb_reaction_4 + 1 4224 v_4(:,nb_reaction_4) = i030(:) 4225 4226 !--- i031: HCO2+ + O -> HCO+ + O2 4227 4228 ! UMIST 4229 4230 i031(:) = 1.e-9 4231 nb_reaction_4 = nb_reaction_4 + 1 4232 v_4(:,nb_reaction_4) = i031(:) 4233 4234 !--- i032: HCO2+ + CO -> HCO+ + CO2 4235 4236 ! UMIST, from Prassad & Huntress 1980 4237 4238 i032(:) = 7.8e-10 4239 nb_reaction_4 = nb_reaction_4 + 1 4240 v_4(:,nb_reaction_4) = i032(:) 4241 4242 !--- i033: H+ + CO2 -> HCO+ + O 4243 4244 ! UMIST, from Smith et al., Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992) 4245 4246 i033(:) = 3.5e-9 4247 nb_reaction_4 = nb_reaction_4 + 1 4248 v_4(:,nb_reaction_4) = i033(:) 4249 4250 4251 !--- i034: CO2+ + H -> HCO+ + O 4252 4253 ! Seen in Fox 2015, from Borodi et al., Int. J. Mass Spectrom. 280, 218-225, 2009 4254 4255 i034(:) = 4.5e-10 4256 nb_reaction_4 = nb_reaction_4 + 1 4257 v_4(:,nb_reaction_4) = i034(:) 4258 4259 !--- i035: CO+ + H2 -> HCO+ + H 4260 4261 !UMIST, from Scott et al., J. Chem. Phys., 106, 3982-3987(1997) 4262 4263 i035(:) = 7.5e-10 4264 nb_reaction_4 = nb_reaction_4 + 1 4265 v_4(:,nb_reaction_4) = i035(:) 4266 4267 !--- i036: HCO+ + e- -> CO + H 4268 4269 !UMIST, from Mitchell, Phys. Rep., 186, 215 (1990) 4270 4271 i036(:) = 2.4e-7 *((300./t_elect(:))**0.69) 4272 nb_reaction_4 = nb_reaction_4 + 1 4273 v_4(:,nb_reaction_4) = i036(:) 4274 4275 !--- i037: CO2+ + H2O -> H2O+ + CO2 4276 4277 !UMIST, from Karpas, Z., Anicich, V.G., and Huntress, W.T., Chem. Phys. Lett., 59, 84 (1978) 4278 4279 i037(:) = 2.04e-9 *((300./t_elect(:))**0.5) 4280 nb_reaction_4 = nb_reaction_4 + 1 4281 v_4(:,nb_reaction_4) = i037(:) 4282 4283 !--- i038: CO+ + H2O -> H2O+ + CO 4284 4285 !UMIST, from Huntress, W.T., McEwan, M.J., Karpas, Z., and Anicich, V.G., Astrophys. J. Supp. Series, 44, 481 (1980) 4286 4287 i038(:) = 1.72e-9*((300./t_elect(:))**0.5) 4288 nb_reaction_4 = nb_reaction_4 + 1 4289 v_4(:,nb_reaction_4) = i038(:) 4290 4291 !--- i039: O+ + H2O -> H2O+ + O 4292 4293 !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978) 4294 4295 i039(:) = 3.2e-9*((300./t_elect(:))**0.5) 4296 nb_reaction_4 = nb_reaction_4 + 1 4297 v_4(:,nb_reaction_4) = i039(:) 4298 4299 !--- i040: N2+ + H2O -> H2O+ + N2 4300 4301 !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978) 4302 4303 i040(:) = 2.3e-9*((300./t_elect(:))**0.5) 4304 nb_reaction_4 = nb_reaction_4 + 1 4305 v_4(:,nb_reaction_4) = i040(:) 4306 4307 !--- i041: N+ + H2O -> H2O+ + N 4308 4309 !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978) 4310 4311 i041(:) = 2.8e-9*((300./t_elect(:))**0.5) 4312 nb_reaction_4 = nb_reaction_4 + 1 4313 v_4(:,nb_reaction_4) = i041(:) 4314 4315 4316 !--- i042: H+ + H2O -> H2O+ + H 4317 4318 !UMIST, from D. Smith, P. Spanel and C. A. Mayhew, Int. J. Mass Spectrom. Ion Proc., 117, 457-473(1992) 4319 4320 i042(:) = 6.9e-9*((300./t_elect(:))**0.5) 4321 nb_reaction_4 = nb_reaction_4 + 1 4322 v_4(:,nb_reaction_4) = i042(:) 4323 4324 !--- i043: H2O+ + O2 -> O2+ + H2O 4325 4326 !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980) 4327 4328 i043(:) = 4.6e-10 4329 nb_reaction_4 = nb_reaction_4 + 1 4330 v_4(:,nb_reaction_4) = i043(:) 4331 4332 !--- i044: H2O+ + CO -> HCO+ + OH 4333 4334 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4335 4336 i044(:) = 5.0e-10 4337 nb_reaction_4 = nb_reaction_4 + 1 4338 v_4(:,nb_reaction_4) = i044(:) 4339 4340 !--- i045: H2O+ + O -> O2+ + H2 4341 4342 !UMIST, from Viggiano, A.A, Howarka, F., Albritton, D.L., Fehsenfeld, F.C., Adams, N.G., and Smith, D., Astrophys. J., 236, 492 (1980) 4343 4344 i045(:) = 4.0e-11 4345 nb_reaction_4 = nb_reaction_4 + 1 4346 v_4(:,nb_reaction_4) = i045(:) 4347 4348 !--- i046: H2O+ + NO -> NO+ + H2O 4349 4350 !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980) 4351 4352 i046(:) = 2.7e-10 4353 nb_reaction_4 = nb_reaction_4 + 1 4354 v_4(:,nb_reaction_4) = i046(:) 4355 4356 !--- i047: H2O+ + e- -> H + H + O 4357 4358 !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295 4359 4360 i047(:) = 3.05e-7*((300./t_elect(:))**0.5) 4361 nb_reaction_4 = nb_reaction_4 + 1 4362 v_4(:,nb_reaction_4) = i047(:) 4363 4364 !--- i048: H2O+ + e- -> H + OH 4365 4366 !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295 4367 4368 i048(:) = 8.6e-8*((300./t_elect(:))**0.5) 4369 nb_reaction_4 = nb_reaction_4 + 1 4370 v_4(:,nb_reaction_4) = i048(:) 4371 4372 !--- i049: H2O+ + e- -> O + H2 4373 4374 !UMIST, from Rosen, S., Derkatch, A., Semaniak, J., et al., 2000, Far. Disc., 115, 295 4375 4376 i049(:) = 3.9e-8*((300./t_elect(:))**0.5) 4377 nb_reaction_4 = nb_reaction_4 + 1 4378 v_4(:,nb_reaction_4) = i049(:) 4379 4380 !--- i050: H2O+ + H2O -> H3O+ + OH 4381 4382 !UMIST, from Huntress, W.T. and Pinizzotto, R.F., J. Chem. Phys., 59, 4742 (1973) 4383 4384 i050(:) = 2.1e-9*((300./t_elect(:))**0.5) 4385 nb_reaction_4 = nb_reaction_4 + 1 4386 v_4(:,nb_reaction_4) = i050(:) 4387 4388 4389 !--- i051: H2O+ + H2 -> H3O+ + H 4390 4391 !UMIST, from A. B. Raksit and P. Warneck, J. Chem. Soc. Faraday Trans., 76, 1084-1092(1980) 4392 4393 i051(:) = 6.4e-10 4394 nb_reaction_4 = nb_reaction_4 + 1 4395 v_4(:,nb_reaction_4) = i051(:) 4396 4397 !--- i052: HCO+ + H2O -> H3O+ + CO 4398 4399 !UMIST, from Adams, N.G., Smith, D., and Grief, D., Int. J. Mass Spectrom. Ion Phys., 26, 405 (1978) 4400 4401 i052(:) = 2.5e-9*((300./t_elect(:))**0.5) 4402 nb_reaction_4 = nb_reaction_4 + 1 4403 v_4(:,nb_reaction_4) = i052(:) 4404 4405 !--- i053: H3O+ + e -> OH + H + H 4406 4407 !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874 4408 4409 i053(:) = 3.05e-7*((300./t_elect(:))**0.5) 4410 nb_reaction_4 = nb_reaction_4 + 1 4411 v_4(:,nb_reaction_4) = i053(:) 4412 4413 !--- i054: H3O+ + e -> H2O + H 4414 4415 !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874 4416 4417 i054(:) = 7.09e-8*((300./t_elect(:))**0.5) 4418 nb_reaction_4 = nb_reaction_4 + 1 4419 v_4(:,nb_reaction_4) = i054(:) 4420 4421 !--- i055: H3O+ + e -> OH + H2 4422 4423 !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874 4424 4425 i055(:) = 5.37e-8*((300./t_elect(:))**0.5) 4426 nb_reaction_4 = nb_reaction_4 + 1 4427 v_4(:,nb_reaction_4) = i055(:) 4428 4429 !--- i056: H3O+ + e -> O + H2 + H 4430 4431 !UMIST, from Novotny, O. et al., J. Phys. Chem. A 2010, 114, 14, 4870-4874 4432 4433 i056(:) = 5.6e-9*((300./t_elect(:))**0.5) 4434 nb_reaction_4 = nb_reaction_4 + 1 4435 v_4(:,nb_reaction_4) = i056(:) 4436 4437 nb_reaction_4 = nb_reaction_4 + 1 4438 v_4(:,nb_reaction_4) = i056(:) 4439 4440 !--- i057: O+ + H2 -> OH+ + H 4441 4442 !UMIST, from Adams, N.G., Smith, D., and Paulson, J.F., J. Chem. Phys., 72, 288 (1980); Smith, D., Adams, N.G., and Miller, T.M., J. Chem. Phys.., 69, 308 (1978) 4443 4444 i057(:) = 1.7e-9 4445 nb_reaction_4 = nb_reaction_4 + 1 4446 v_4(:,nb_reaction_4) = i057(:) 4447 4448 !--- i058: OH+ + O -> O2+ + H 4449 4450 !UMIST, from Prasad & Huntress, 1980, ApJS, 43, 1 4451 4452 i058(:) = 7.1e-10 4453 nb_reaction_4 = nb_reaction_4 + 1 4454 v_4(:,nb_reaction_4) = i058(:) 4455 4456 !--- i059: OH+ + CO2 -> HCO2+ + O 4457 4458 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4459 4460 i059(:) = 1.44e-9 4461 nb_reaction_4 = nb_reaction_4 + 1 4462 v_4(:,nb_reaction_4) = i059(:) 4463 4464 !--- i060: OH+ + CO -> HCO+ + O 4465 4466 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4467 4468 i060(:) = 1.05e-9 4469 nb_reaction_4 = nb_reaction_4 + 1 4470 v_4(:,nb_reaction_4) = i060(:) 4471 4472 !--- i061: OH+ + NO -> NO+ + OH (tasa de reacción UMIST 3.59e-10) 4473 4474 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4475 4476 i061(:) = 3.59e-10 4477 nb_reaction_4 = nb_reaction_4 + 1 4478 v_4(:,nb_reaction_4) = i061(:) 4479 4480 !--- i062: OH+ + H2 -> H2O+ + H (tasa de reacción UMIST 1.01e-9, 4481 4482 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4483 4484 i062(:) = 1.01e-9 4485 nb_reaction_4 = nb_reaction_4 + 1 4486 v_4(:,nb_reaction_4) = i062(:) 4487 4488 !--- i063: OH+ + O2 -> O2+ + OH (tasa de reacción UMIST 5.9e-10 4489 4490 !UMIST, from Jones, J.D.C., Birkinshaw, K., and Twiddy, N.D., Chem. Phys. Lett., 77, 484 (1981) 4491 4492 i063(:) = 5.9e-10 4493 nb_reaction_4 = nb_reaction_4 + 1 4494 v_4(:,nb_reaction_4) = i063(:) 4495 4496 end if ! ok_ionchem 4497 4498 !---------------------------------------------------------------------- 3128 4499 ! reactions avec 02(Dg) 3129 4500 !---------------------------------------------------------------------- 3130 4501 3131 !--- i001: O2(Dg) + CO2 -> O2 + CO2+ hv3132 3133 ! Krasnopolsky (2010a) 3134 3135 i001(:) = 1.E-204502 !--- j001: O2(Dg) + (CO2 and O) -> O2 + (CO2 and O) + hv 4503 4504 ! Krasnopolsky (2010a) for CO2 & Clark and Wayne, 1969 for O (JPL) 4505 4506 j001(:) = 1.E-20 * c(:,i_co2) !+ 2.E-16 * c(:,i_o) 3136 4507 3137 4508 nb_phot = nb_phot + 1 3138 v_phot(:,nb_phot) = i001(:)*c(:,i_co2)3139 3140 !--- i002: O2(Dg) -> O2 + hv4509 v_phot(:,nb_phot) = j001(:) 4510 4511 !--- j002: O2(Dg) -> O2 + hv 3141 4512 3142 4513 ! Lafferty et al; (1998) 3143 4514 3144 i002(:) = 2.2E-44515 j002(:) = 2.2E-4 3145 4516 3146 4517 nb_phot = nb_phot + 1 3147 v_phot(:,nb_phot) = i002(:) 4518 v_phot(:,nb_phot) = j002(:) 4519 4520 !!! TEST: artificial increase of CO2 photodissociation 4521 if (tuneupperatm) then 4522 !-- TuneA 4523 ! v_phot(65:78,4) = v_phot(65:78,4)*10. 4524 ! v_phot(60:64,4) = v_phot(60:64,4)*3. 4525 ! v_phot(55:59,3) = v_phot(55:59,3)*2. 4526 !-- 4527 !-- TuneB 4528 ! v_phot(65:78,4) = v_phot(65:78,4)*10. 4529 ! v_phot(55:59,3) = v_phot(55:59,3)*5. 4530 !-- 4531 !-- TuneC 4532 ! VCD 1.1 tuning 4533 ! v_phot(65:78,4) = v_phot(65:78,4)*10. 4534 ! v_phot(52:59,3) = v_phot(52:59,3)*5. 4535 !-- 4536 !-- TuneE 4537 ! VCD 2.0 tuning 4538 ! v_phot(65:90,4) = v_phot(65:90,4)*10. ! CO2 + hv ==> O(1D) + CO 4539 !-- 4540 ! VCD 2.1 tuning 4541 v_phot(65:90,4) = v_phot(65:90,4)*8. ! CO2 + hv ==> O(1D) + CO 4542 !-- 4543 ! v_phot(:,4) = v_phot(:,4)*10. 4544 !do ij=3,4 4545 ! v_phot(:,ij) = v_phot(:,ij)*10. 4546 !end do 4547 endif 4548 !!!!!!!!!!!!!!! 3148 4549 3149 4550 return … … 3176 4577 ! number of processes treated numerically as photodissociations 3177 4578 3178 real , dimension(nlayer,nesp) :: c ! number densities3179 real , dimension(nlayer, nb_phot_max) :: v_phot3180 real , dimension(nlayer,nb_reaction_3_max) :: v_33181 real , dimension(nlayer,nb_reaction_4_max) :: v_44579 real (kind = 8), dimension(nlayer,nesp) :: c ! number densities 4580 real (kind = 8), dimension(nlayer, nb_phot_max) :: v_phot 4581 real (kind = 8), dimension(nlayer,nb_reaction_3_max) :: v_3 4582 real (kind = 8), dimension(nlayer,nb_reaction_4_max) :: v_4 3182 4583 3183 4584 ! output 3184 4585 3185 real , dimension(nesp,nesp), intent(out) :: mat ! matrix3186 real , dimension(nesp), intent(out) :: prod, loss, lossconc4586 real (kind = 8), dimension(nesp,nesp), intent(out) :: mat ! matrix 4587 real (kind = 8), dimension(nesp), intent(out) :: prod, loss, lossconc 3187 4588 3188 4589 ! local … … 3194 4595 integer :: iphot,i3,i4 3195 4596 3196 real :: eps, eps_4 ! implicit/explicit coefficient4597 real(kind = 8) :: eps, eps_4 ! implicit/explicit coefficient 3197 4598 3198 4599 ! initialisations … … 3299 4700 3300 4701 real :: dtold, ctimestep 3301 real , dimension(nesp) :: cold, ccur3302 real , dimension(nesp,nesp) :: mat13303 real , dimension(nesp) :: prod, loss4702 real (kind = 8), dimension(nesp) :: cold, ccur 4703 real (kind = 8), dimension(nesp,nesp) :: mat1 4704 real (kind = 8), dimension(nesp) :: prod, loss 3304 4705 real :: dens 3305 4706 real :: lon, lat … … 3311 4712 ! local 3312 4713 3313 real , dimension(nesp) :: cnew3314 real , dimension(nesp,nesp) :: mat3315 real :: atol, ratio, e, es, coef4714 real (kind = 8), dimension(nesp) :: cnew 4715 real (kind = 8), dimension(nesp,nesp) :: mat 4716 real (kind = 8) :: atol, ratio, e, es, coef 3316 4717 3317 4718 integer :: code, iesp, iter … … 3323 4724 ! parameters 3324 4725 3325 real, parameter :: dtmin = 10. ! minimum time step (s) 3326 real, parameter :: vmrtol = 1.e-11 ! absolute tolerance on vmr 3327 real, parameter :: rtol = 0.05 ! rtol recommended value : 0.1-0.02 3328 integer, parameter :: niter = 3 ! number of iterations 3329 real, parameter :: coefmax = 2. 3330 real, parameter :: coefmin = 0.1 3331 logical :: fast_guess = .true. 4726 real (kind = 8), parameter :: dtmin = 10. ! minimum time step (s) 4727 real (kind = 8), parameter :: vmrtol = 1.e-11 ! absolute tolerance on vmr 4728 real (kind = 8), parameter :: rtol = 0.05 ! rtol recommended value : 0.1-0.02 4729 integer, parameter :: niter = 3 ! number of iterations 4730 real (kind = 8), parameter :: coefmax = 2. 4731 real (kind = 8), parameter :: coefmin = 0.1 4732 logical :: fast_guess = .true. 4733 3332 4734 3333 4735 dttest = dtold ! dttest = dtold = dt_guess … … 3403 4805 3404 4806 !====================================================================== 4807 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4808 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4809 ! CE CODE EST OBSOLETE !! A NE PAS UTILISER !!!!!!!!!!!!!!!!!!!!!!!!!!! 4810 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4811 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3405 4812 3406 4813 SUBROUTINE rate_save( & -
trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90
r2799 r2836 569 569 implicit none 570 570 571 #include "clesphys.h" ! fixed_euv_value 572 571 573 ! input 572 574 … … 582 584 integer, parameter :: kdata = 20000 ! max dimension of input solar flux 583 585 integer :: msun ! choice of solar flux 584 integer :: iw, nhead, ihead, n, i, ierr, kin 586 integer :: iw, nhead, ihead, n, i, ierr, kin1, kin2 585 587 586 588 real, parameter :: deltax = 1.e-4 587 real, dimension(kdata) :: x1, y1 ! input solar flux 588 real, dimension(nw) :: yg1 ! gridded solar flux 589 590 character(len=100) :: fil 591 592 kin = 10 ! input logical unit 593 589 ! atlas1_thuillier_tuv.txt 590 real, dimension(kdata) :: x1, y1 ! input solar flux - HIGH SOLAR ACTIVITY 591 real :: E107_max = 196 592 ! atlas3_thuillier_tuv.txt 593 real, dimension(kdata) :: x2, y2 ! input solar flux - "LOW" SOLAR ACTIVITY 594 real :: E107_min = 97 595 ! be careful, x2 need to be equal to x1 596 real, dimension(kdata) :: x0, y0 ! input solar flux - interpolated solar activity 597 real, dimension(nw) :: yg0 ! gridded solar flux 598 real :: factor_interp 599 600 character(len=100) :: fil1, fil2 601 602 kin1 = 10 ! input logical unit 603 kin2 = 11 ! input logical unit 604 594 605 ! select desired extra-terrestrial solar irradiance, using msun: 595 606 596 ! 18 = atlas1_thuillier_tuv.txt 0-900 nm March 1992 607 ! 18 = atlas1_thuillier_tuv.txt 0-900 nm March 29, 1992 608 ! Article: F10.7 = 192 s.f.u | <F10.7>81d = 171 s.f.u | <Rz> = 121 (sunsport number) 597 609 ! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004 598 599 msun = 18 600 601 if (msun == 18) THEN 602 603 fil = 'solar_fluxes/atlas1_thuillier_tuv.txt' 604 print*, 'solar flux : ', fil 605 606 if(is_master) then 607 608 open(kin, file=fil, status='old', iostat=ierr) 609 610 if (ierr /= 0) THEN 611 write(*,*)'cant find solar flux : ', fil 612 write(*,*)'It should be in : INPUT/solar_fluxes' 613 write(*,*)'1) You can change this directory address in ' 614 write(*,*)' callphys.def with datadir=/path/to/dir' 615 write(*,*)'2) If necessary, /solar fluxes (and other datafiles)' 616 write(*,*)' can be obtained online on:' 617 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 618 stop 619 end if 620 621 nhead = 9 622 n = 19193 623 DO ihead = 1, nhead 624 READ(kin,*) 625 ENDDO 626 DO i = 1, n 627 READ(kin,*) x1(i), y1(i) 628 y1(i) = y1(i)*1.e-3 ! mw -> w 629 ENDDO 630 CLOSE (kin) 631 CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.) 632 CALL addpnt(x1,y1,kdata,n, 0.,0.) 633 CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.) 634 CALL addpnt(x1,y1,kdata,n, 1.e+38,0.) 635 CALL inter2(nw,wl,yg1,n,x1,y1,ierr) 636 637 IF (ierr .NE. 0) THEN 638 WRITE(*,*) ierr, fil 639 STOP 640 ENDIF 610 ! Model SOLAR2000 version 2015/12: E10.7 = 195.8 s.f.u | <E10.7>81d = 196.9 s.f.u 611 ! Choix de la limite haute: E10.7 = 196 s.f.u 612 613 ! 20 = atlas3_thuillier_tuv.txt 0-900 nm November 11, 1994 614 ! Article: F10.7 = 77.5 s.f.u | <F10.7>81d = 83.5 s.f.u | <Rz> = 20 (sunsport number) 615 ! Thuillier et al., Adv. Space. Res., 34, 256-261, 2004 616 ! Model SOLAR2000 version 2015/12: E10.7 = 96.9 s.f.u | <E10.7>81d = 100.0 s.f.u 617 ! Choix de la limite basse: E10.7 = 97 s.f.u 618 619 if (fixed_euv_value .ge. 196) THEN 620 msun = 18 621 print*, 'Atlas1 spectrum Thuiller chosen' 622 else 623 msun = 19 624 print*, 'interpolated Spectrum case' 625 end if 626 627 IF (is_master) THEN 628 629 fil1 = 'solar_fluxes/atlas1_thuillier_tuv.txt' 630 print*, 'solar flux : ', fil1 631 632 open(kin1, file=fil1, status='old', iostat=ierr) 633 634 if (ierr /= 0) THEN 635 write(*,*)'cant find solar flux : ', fil1 636 write(*,*)'It should be in : INPUT/solar_fluxes' 637 write(*,*)'1) You can change this directory address in ' 638 write(*,*)' callphys.def with datadir=/path/to/dir' 639 write(*,*)'2) If necessary, /solar fluxes (and other datafiles)' 640 write(*,*)' can be obtained online on:' 641 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 642 stop 643 end if 644 645 nhead = 9 646 n = 19193 647 DO ihead = 1, nhead 648 READ(kin1,*) 649 ENDDO 650 DO i = 1, n 651 READ(kin1,*) x1(i), y1(i) 652 y1(i) = y1(i)*1.e-3 ! mw -> w 653 ENDDO 654 CLOSE (kin1) 655 656 fil2 = 'solar_fluxes/atlas3_thuillier_tuv.txt' 657 print*, 'solar flux : ', fil2 658 659 open(kin2, file=fil2, status='old', iostat=ierr) 660 661 if (ierr /= 0) THEN 662 write(*,*)'cant find solar flux : ', fil2 663 write(*,*)'It should be in : INPUT/solar_fluxes' 664 write(*,*)'1) You can change this directory address in ' 665 write(*,*)' callphys.def with datadir=/path/to/dir' 666 write(*,*)'2) If necessary, /solar fluxes (and other datafiles)' 667 write(*,*)' can be obtained online on:' 668 write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir' 669 stop 670 end if 671 672 nhead = 9 673 n = 19193 674 DO ihead = 1, nhead 675 READ(kin2,*) 676 ENDDO 677 DO i = 1, n 678 READ(kin2,*) x2(i), y2(i) 679 y2(i) = y2(i)*1.e-3 ! mw -> w 680 ENDDO 681 CLOSE (kin2) 682 683 IF (msun .eq. 18) THEN 684 DO i = 1, n 685 x0(i) = x1(i) 686 y0(i) = y1(i) 687 ENDDO 688 ELSE 689 ! INTERPOLATION BETWEEN E107_min to E107_max and extrapolation below E107_min 690 factor_interp=(fixed_euv_value-E107_min)/(E107_max-E107_min) 691 DO i = 1, n 692 x0(i) = x1(i) 693 y0(i) = y2(i) + (y1(i) - y2(i))* factor_interp 694 ENDDO 695 ENDIF ! msun .ne. 18 696 697 CALL addpnt(x0,y0,kdata,n,x0(1)*(1.-deltax),0.) 698 CALL addpnt(x0,y0,kdata,n, 0.,0.) 699 CALL addpnt(x0,y0,kdata,n,x0(n)*(1.+deltax),0.) 700 CALL addpnt(x0,y0,kdata,n, 1.e+38,0.) 701 CALL inter2(nw,wl,yg0,n,x0,y0,ierr) 702 703 IF (ierr .NE. 0) THEN 704 WRITE(*,*) ierr, fil1, fil2 705 STOP 706 ENDIF 641 707 642 708 ! convert to photon.s-1.nm-1.cm-2 643 709 ! 5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8) 644 710 645 646 f(iw) = yg1(iw)*wc(iw)*5.039e11647 711 DO iw = 1, nw-1 712 f(iw) = yg0(iw)*wc(iw)*5.039e11 713 ENDDO 648 714 649 715 endif !is_master 650 716 651 717 call bcast(f) 652 653 end if654 718 655 719 end subroutine rdsolarflux -
trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F
r2795 r2836 13 13 14 14 implicit none 15 15 16 16 ! input 17 17 … … 433 433 434 434 ! production of o(1d) for wavelengths shorter than 167 nm 435 436 yieldco2(:) = 1437 435 438 436 if (wc(l) >= 167.) then -
trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
r2795 r2836 80 80 use nonoro_gwd_ran_mod, only: nonoro_gwd_ran 81 81 use sed_and_prod_mad, only: aer_sedimentation, drop_sedimentation 82 use iono_h, only: temp_elect, temp_ion 82 83 #ifdef CPP_XIOS 83 84 use xios_output_mod, only: initialize_xios_output, … … 237 238 c Declaration des procedures appelees 238 239 c 239 EXTERNAL ajsec ! ajustement sec240 EXTERNAL clmain ! couche limite241 EXTERNAL clmain_ideal 242 EXTERNAL hgardfou ! verifier les temperatures243 c EXTERNAL orbite ! calculer l'orbite244 EXTERNAL phyetat0 ! lire l'etat initial de la physique245 EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique246 EXTERNAL radlwsw ! rayonnements solaire et infrarouge247 ! EXTERNAL suphec ! initialiser certaines constantes248 EXTERNAL transp ! transport total de l'eau et de l'energie240 EXTERNAL ajsec ! ajustement sec 241 EXTERNAL clmain ! couche limite 242 EXTERNAL clmain_ideal ! couche limite simple 243 EXTERNAL hgardfou ! verifier les temperatures 244 c EXTERNAL orbite ! calculer l'orbite 245 EXTERNAL phyetat0 ! lire l'etat initial de la physique 246 EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique 247 EXTERNAL radlwsw ! rayonnements solaire et infrarouge 248 ! EXTERNAL suphec ! initialiser certaines constantes 249 EXTERNAL transp ! transport total de l'eau et de l'energie 249 250 EXTERNAL printflag 250 251 EXTERNAL zenang … … 270 271 c 271 272 CXXX PB 272 REAL fluxt(klon,klev) ! flux turbulent de chaleur273 REAL fluxu(klon,klev) ! flux turbulent de vitesse u274 REAL fluxv(klon,klev) ! flux turbulent de vitesse v273 REAL fluxt(klon,klev) ! flux turbulent de chaleur 274 REAL fluxu(klon,klev) ! flux turbulent de vitesse u 275 REAL fluxv(klon,klev) ! flux turbulent de vitesse v 275 276 c 276 277 REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique … … 278 279 REAL flux_ec(klon,klev) ! flux de chaleur Ec 279 280 c 280 REAL tmpout(klon,klev) ! K s-1281 REAL tmpout(klon,klev) ! [K/s] 281 282 c 282 283 REAL dist, rmu0(klon), fract(klon) … … 334 335 c Tendencies due to molecular viscosity and conduction 335 336 real d_t_conduc(klon,klev) ! [K/s] 336 real d_u_molvis(klon,klev) ! (m/s)/s337 real d_v_molvis(klon,klev) ! (m/s)/s337 real d_u_molvis(klon,klev) ! [m/s] /s 338 real d_v_molvis(klon,klev) ! [m/s] /s 338 339 339 340 c Tendencies due to molecular diffusion … … 820 821 IF (ancien_ok) THEN 821 822 DO k = 1, klev 822 DO i = 1, klon823 DO i = 1, klon 823 824 d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime 824 825 d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime 825 ENDDO826 ENDDO 826 827 ENDDO 827 828 … … 860 861 c 861 862 DO k = 1, klev 862 DO i = 1, klon863 DO i = 1, klon 863 864 zphi(i,k) = pphi(i,k) + pphis(i) 864 ENDDO865 ENDDO 865 866 ENDDO 866 867 … … 996 997 DO ig=1,klon 997 998 t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K] 998 u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s999 v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s999 u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! [m/s] 1000 v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! [m/s] 1000 1001 ENDDO 1001 1002 ENDDO … … 1609 1610 DO ig=1,klon 1610 1611 d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k) 1611 1612 1612 ENDDO 1613 1613 ENDDO … … 2072 2072 CALL send_xios_field("dugwo",d_u_oro) 2073 2073 CALL send_xios_field("dugwno",d_u_hin) 2074 CALL send_xios_field("dvgwno",-1.*d_v_hin) 2074 2075 CALL send_xios_field("dumolvis",d_u_molvis) 2075 2076 c VENUS: regardee a l envers!!!!!!!!!!!!!!! … … 2102 2103 ! production and destruction rate, cm-3.s-1 2103 2104 ! Beware of the context*.xml file !! 2104 ! if ((tr_scheme == 3) .and. (ok_chem)) THEN 2105 ! do iq = 1,nqmax - nmicro 2106 ! !if ((iq == i_co) .or. (iq == i_o)) then 2107 ! ! call send_xios_field('prod_'+tname(iq), prod_tr(:,:,iq)) 2108 ! ! call send_xios_field('loss_'+tname(iq), loss_tr(:,:,iq)) 2109 ! !end if 2110 ! if (iq == i_o) then 2111 ! call send_xios_field('prod_o', prod_tr(:,:,iq)) 2112 ! call send_xios_field('loss_o', loss_tr(:,:,iq)) 2113 ! end if 2114 ! if (iq == i_co) then 2115 ! call send_xios_field('prod_co', prod_tr(:,:,iq)) 2116 ! call send_xios_field('loss_co', loss_tr(:,:,iq)) 2117 ! end if 2118 ! end do 2119 ! end if 2105 if ((tr_scheme == 3) .and. (ok_chem)) THEN 2106 do iq = 1,nqmax - nmicro 2107 if ((iq.eq.i_o).or.(iq.eq.i_co)) THEN 2108 call send_xios_field("prod_"//tname(iq), 2109 $ prod_tr(:,:,iq)) 2110 call send_xios_field("loss_"//tname(iq), 2111 $ loss_tr(:,:,iq)) 2112 end if 2113 !if (iq.eq.i_o) then 2114 ! call send_xios_field('prod_o', prod_tr(:,:,iq)) 2115 ! call send_xios_field('loss_o', loss_tr(:,:,iq)) 2116 !end if 2117 !if (iq.eq.i_co) then 2118 ! call send_xios_field('prod_co', prod_tr(:,:,iq)) 2119 ! call send_xios_field('loss_co', loss_tr(:,:,iq)) 2120 !end if 2121 end do 2122 end if 2120 2123 2121 2124 ! tracers in gas phase, volume mixing ratio … … 2130 2133 do iq = 1,nqmax - nmicro 2131 2134 col_dens_tr(:,iq)=0. 2132 do k = 1, klev 2133 col_dens_tr(:,iq) = col_dens_tr(:,iq) + 2134 $ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG 2135 end do 2136 call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq)) 2135 if (type_tr(iq).eq.1) THEN 2136 do k = 1, klev 2137 col_dens_tr(:,iq) = col_dens_tr(:,iq) + 2138 $ tr_seri(:,k,iq)* (paprs(:,k)-paprs(:,k+1)) / RG 2139 end do 2140 call send_xios_field("col_"//tname(iq),col_dens_tr(:,iq)) 2141 end if 2137 2142 end do 2138 2143 … … 2151 2156 ! aeronomical emissions 2152 2157 2153 ! call send_xios_field("no_emission",no_emission)2154 ! call send_xios_field("o2_emission",o2_emission)2158 call send_xios_field("no_emis",no_emission) 2159 call send_xios_field("o2_emis",o2_emission) 2155 2160 2156 2161 ! chemical iterations -
trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
r2795 r2836 21 21 use chemparam_mod 22 22 use conc, only: mmean,rnew 23 use photolysis_mod, only: init_photolysis, nphot 24 use iono_h, only: temp_elect 23 25 #ifdef CPP_XIOS 24 26 use xios_output_mod, only: send_xios_field … … 63 65 real :: lon_sun 64 66 real :: zlocal(nlev) ! altitude for photochem (km) 65 67 real :: t_elec(nlev) ! electron temperature [K] 68 69 integer, parameter :: t_elec_origin=2 70 !Electronic temperature. Index 1 -> Theis et al. 1980 - model data ; Index 2-> Theis et al. 1984 - model data 71 66 72 integer :: i, iq 67 73 integer :: ilon, ilev … … 76 82 real, dimension(nlev) :: no_em 77 83 real, dimension(nlev) :: o2_em 78 84 85 integer, save :: nb_reaction_3_max ! number of quadratic reactions 86 integer, save :: nb_reaction_4_max ! number of bimolecular reactions 87 integer, save :: nquench ! number of quenching + heterogeneous reactions 88 integer, save :: nphotion ! number of photoionizations 89 integer, save :: nb_reaction_4_ion ! quadratic reactions for ionosphere 90 ! integer, save :: nb_reaction_4_deut ! quadratic reactions for deuterium chem 91 integer, save :: nb_phot_max ! total number of photolysis+photoionizations+quenching reactions 92 93 ! tracer indexes for the EUV heating: 94 !!! ATTENTION. These values have to be identical to those in euvheat.F90 95 !!! If the values are changed there, the same has to be done here !!! 96 97 ! integer,parameter :: i_co2=1 98 ! integer,parameter :: i_n2=13 99 ! integer,parameter :: i_n=14 100 ! integer,parameter :: i_o=3 101 ! integer,parameter :: i_co=4 102 103 integer,parameter :: ix_co2 = 1 104 integer,parameter :: ix_co = 2 105 integer,parameter :: ix_o = 3 106 integer,parameter :: ix_o1d = 4 107 integer,parameter :: ix_o2 = 5 108 integer,parameter :: ix_o3 = 6 109 integer,parameter :: ix_h = 7 110 integer,parameter :: ix_h2 = 8 111 integer,parameter :: ix_oh = 9 112 integer,parameter :: ix_ho2 = 10 113 integer,parameter :: ix_h2o2 = 11 114 integer,parameter :: ix_h2o = 12 115 integer,parameter :: ix_n = 13 116 integer,parameter :: ix_n2d = 14 117 integer,parameter :: ix_no = 15 118 integer,parameter :: ix_no2 = 16 119 integer,parameter :: ix_n2 = 17 120 121 ! NEED TO BE THE SAME THAN IN EUVHEAT.F90 122 integer,parameter :: nespeuv=17 ! Number of species considered (11, 12 or 17 (with nitrogen)) 123 124 real :: vmr_dens_euv(nlon,nlev,nespeuv) ! local species density for EUV heating 125 79 126 !=================================================================== 80 127 ! first call : initialisations … … 90 137 91 138 !------------------------------------------------------------------- 139 ! Determination of chemistry reaction number 140 !------------------------------------------------------------------- 141 142 if (ok_chem) then 143 ! set number of reactions, depending on ion chemistry or not 144 nb_reaction_4_ion = 64 145 !nb_reaction_4_deut = 35 146 147 !Default numbers if no ion and no deuterium chemistry included 148 149 nb_reaction_4_max = 98 ! set number of bimolecular reactions 150 nb_reaction_3_max = 12 ! set number of quadratic reactions 151 nquench = 13 ! set number of quenching + heterogeneous 152 nphotion = 0 ! set number of photoionizations 153 154 if (ok_ionchem) then 155 nb_reaction_4_max = nb_reaction_4_max+nb_reaction_4_ion 156 nphotion = 18 ! set number of photoionizations 157 endif 158 !if(deutchem) then 159 ! nb_reaction_4_max = nb_reaction_4_max + nb_reaction_4_deut 160 !end if 161 162 !nb_phot_max is the total number of processes that are treated 163 !numerically as a photolysis: 164 165 nb_phot_max = nphot + nphotion + nquench 166 167 !------------------------------------------------------------------- 92 168 ! case of tracers re-initialisation with chemistry 93 169 !------------------------------------------------------------------- 94 95 if (reinit_trac .and. ok_chem) then 96 97 !!! in this reinitialisation, trac is VOLUME mixing ratio 98 ! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 99 ! convert mass to volume mixing ratio 100 do iq = 1,nqmax - nmicro 101 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq) 102 end do 103 print*, "SO2 is re-initialised" 104 if (i_so2 /= 0) then 105 trac(:,1:22,i_so2) = 10.e-6 106 107 ! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!! 108 ! print*, "Tracers are re-initialised" 109 ! trac(:,:,:) = 1.0e-30 110 111 ! if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 112 ! $ .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0) 113 ! $ .and. (i_co2 /= 0)) then 114 115 ! trac(:,1:22,i_ocs) = 3.e-6 116 ! trac(:,1:22,i_co) = 25.e-6 117 ! trac(:,:,i_hcl) = 0.4e-6 118 ! trac(:,1:22,i_so2) = 7.e-6 119 ! trac(:,1:22,i_h2o) = 30.e-6 120 ! trac(:,:,i_n2) = 0.35e-1 121 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 122 123 ! ensure that sum of mixing ratios = 1 124 125 trac_sum(:,:) = 0. 126 127 do iq = 1,nqmax - nmicro 128 if (iq /= i_co2) then 129 trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq) 130 end if 131 end do 132 133 ! initialise co2 134 135 trac(:,:,i_co2) = 1. - trac_sum(:,:) 136 137 else 138 write(*,*) "at least one tracer is missing : stop" 139 stop 140 end if 141 142 ! update mmean 143 144 mmean(:,:) = 0. 145 do iq = 1,nqmax - nmicro 146 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq) 147 enddo 148 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 149 150 ! convert volume to mass mixing ratio 151 152 do iq = 1,nqmax - nmicro 153 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:) 154 end do 155 156 end if ! reinit_trac 170 171 if (reinit_trac .and. ok_chem) then 172 173 !!! in this reinitialisation, trac is VOLUME mixing ratio 174 ! ONLY SO2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 175 ! convert mass to volume mixing ratio 176 do iq = 1,nqmax - nmicro 177 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq) 178 end do 179 print*, "SO2 is re-initialised" 180 if (i_so2 /= 0) then 181 trac(:,1:22,i_so2) = 10.e-6 182 183 ! AL TRACERS!!!!!!!!!!!!!!!!!!!!!!!!!!!! 184 ! print*, "Tracers are re-initialised" 185 ! trac(:,:,:) = 1.0e-30 186 187 ! if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0) 188 ! $ .and. (i_so2 /= 0) .and. (i_h2o /= 0) 189 ! $ .and. (i_n2 /= 0) .and. (i_co2 /= 0)) then 190 191 ! trac(:,1:22,i_ocs) = 3.e-6 192 ! trac(:,1:22,i_co) = 25.e-6 193 ! trac(:,:,i_hcl) = 0.4e-6 194 ! trac(:,1:22,i_so2) = 7.e-6 195 ! trac(:,1:22,i_h2o) = 30.e-6 196 ! trac(:,:,i_n2) = 0.35e-1 197 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 198 199 ! ensure that sum of mixing ratios = 1 200 201 trac_sum(:,:) = 0. 202 203 do iq = 1,nqmax - nmicro 204 if (iq /= i_co2) then 205 trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq) 206 end if 207 end do 208 209 ! initialise co2 210 211 trac(:,:,i_co2) = 1. - trac_sum(:,:) 212 213 else 214 write(*,*) "at least one tracer is missing : stop" 215 stop 216 end if 217 218 ! update mmean 219 220 mmean(:,:) = 0. 221 do iq = 1,nqmax - nmicro 222 mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq) 223 enddo 224 rnew(:,:) = 8.314/mmean(:,:)*1.e3 ! J/kg K 225 226 ! convert volume to mass mixing ratio 227 228 do iq = 1,nqmax - nmicro 229 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:) 230 end do 231 232 end if ! reinit_trac 233 234 end if ! ok_chem 157 235 158 236 !------------------------------------------------------------------- … … 286 364 lon_local(:) = lon(:)*rpi/180. 287 365 lat_local(:) = lat(:)*rpi/180. 288 366 367 if (ok_ionchem) then 368 369 vmr_dens_euv(:,:,ix_co2) = ztrac(:,:,i_co2) ! CO2 370 vmr_dens_euv(:,:,ix_co) = ztrac(:,:,i_co) ! CO 371 vmr_dens_euv(:,:,ix_o) = ztrac(:,:,i_o) ! O 372 vmr_dens_euv(:,:,ix_o1d) = ztrac(:,:,i_o1d) ! O(1D) 373 vmr_dens_euv(:,:,ix_o2) = ztrac(:,:,i_o2) ! O2 374 vmr_dens_euv(:,:,ix_o3) = ztrac(:,:,i_o3) ! O3 375 vmr_dens_euv(:,:,ix_h) = ztrac(:,:,i_h) ! H 376 vmr_dens_euv(:,:,ix_h2) = ztrac(:,:,i_h2) ! H2 377 vmr_dens_euv(:,:,ix_oh) = ztrac(:,:,i_oh) ! OH 378 vmr_dens_euv(:,:,ix_ho2) = ztrac(:,:,i_ho2) ! HO2 379 vmr_dens_euv(:,:,ix_h2o2)= ztrac(:,:,i_h2o2) ! H2O2 380 vmr_dens_euv(:,:,ix_h2o) = ztrac(:,:,i_h2o) ! H2O 381 vmr_dens_euv(:,:,ix_n) = ztrac(:,:,i_n) ! N 382 vmr_dens_euv(:,:,ix_n2d) = ztrac(:,:,i_n2d) ! N(2D) 383 vmr_dens_euv(:,:,ix_no) = ztrac(:,:,i_no) ! NO 384 vmr_dens_euv(:,:,ix_no2) = ztrac(:,:,i_no2) ! NO2 385 vmr_dens_euv(:,:,ix_n2) = ztrac(:,:,i_n2) ! N2 386 387 end if 388 289 389 do ilon = 1,nlon 290 390 zlocal(:)=zzlay(ilon,:)/1000. 291 391 292 392 ! solar zenith angle 293 393 … … 296 396 $ *sin(lon_local(ilon))*sin(lon_sun))*180./rpi 297 397 398 ! electron temperature 399 do ilev = 1, nlev 400 t_elec(ilev) = temp_elect(zlocal(ilev),temp(ilon,ilev), 401 $ sza_local,t_elec_origin) 402 end do 403 298 404 call photochemistry_venus(nlev, nlon, zlocal, pdtphys, 299 $ pplay(ilon,:)/100., 300 $ temp(ilon,:), 301 $ ztrac(ilon,:,:), 302 $ mmean(ilon,:), 303 $ sza_local, 304 $ lon(ilon), lat(ilon), 305 $ nqmax, iter(ilon,:), 306 $ prod_tr(ilon,:,:), 307 $ loss_tr(ilon,:,:), 308 $ no_em, o2_em) 405 $ ok_jonline,ok_ionchem,tuneupperatm, 406 $ nb_reaction_3_max,nb_reaction_4_max, 407 $ nb_phot_max,nphotion,ilon, 408 $ pplay(ilon,:)/100., 409 $ temp(ilon,:), t_elec(:), 410 $ ztrac(ilon,:,:),vmr_dens_euv(ilon,:,:), 411 $ mmean(ilon,:), 412 $ sza_local, 413 $ lon(ilon), lat(ilon), 414 $ nqmax, nespeuv, iter(ilon,:), 415 $ prod_tr(ilon,:,:), 416 $ loss_tr(ilon,:,:), 417 $ no_em, o2_em) 309 418 310 419 no_emi(ilon,:) = no_em(:)
Note: See TracChangeset
for help on using the changeset viewer.
