| 1 | SUBROUTINE chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress, |
|---|
| 2 | $ zlocal,zenit,ptimestep,zday) |
|---|
| 3 | |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | c======================================================================= |
|---|
| 6 | c subject: |
|---|
| 7 | c -------- |
|---|
| 8 | c Computing chemical variations in the thermosphere |
|---|
| 9 | c |
|---|
| 10 | c author: MAC July 2003 |
|---|
| 11 | c ------ |
|---|
| 12 | c modifications: |
|---|
| 13 | c ------------- |
|---|
| 14 | c Ehouarn Sept 2008: added handling of tracers by their names |
|---|
| 15 | c======================================================================= |
|---|
| 16 | c |
|---|
| 17 | c 0. Declarations : |
|---|
| 18 | c ------------------ |
|---|
| 19 | c |
|---|
| 20 | #include "dimensions.h" |
|---|
| 21 | #include "dimphys.h" |
|---|
| 22 | #include "comcstfi.h" |
|---|
| 23 | #include "callkeys.h" |
|---|
| 24 | #include "comdiurn.h" |
|---|
| 25 | #include "param.h" |
|---|
| 26 | #include "param_v3.h" |
|---|
| 27 | #include "chimiedata.h" |
|---|
| 28 | #include "conc.h" |
|---|
| 29 | #include"tracer.h" |
|---|
| 30 | c----------------------------------------------------------------------- |
|---|
| 31 | c Input/Output |
|---|
| 32 | c ------------ |
|---|
| 33 | INTEGER lswitch,ig |
|---|
| 34 | |
|---|
| 35 | REAL zday,zycol(nlayermx,nqmx) |
|---|
| 36 | REAL ptimestep |
|---|
| 37 | real zenit |
|---|
| 38 | real ztemp(nlayermx) |
|---|
| 39 | real zdens(nlayermx) |
|---|
| 40 | real zpress(nlayermx) ! in mbar |
|---|
| 41 | c |
|---|
| 42 | c Local variables : |
|---|
| 43 | c ----------------- |
|---|
| 44 | INTEGER nlayer,l,nesptherm,iq |
|---|
| 45 | parameter (nesptherm = 11) |
|---|
| 46 | real tmean |
|---|
| 47 | real aux1(nlayermx),aux2(nlayermx) |
|---|
| 48 | real zlocal(nlayermx) |
|---|
| 49 | real rm(nlayermx,nesptherm) !number density (cm-3) |
|---|
| 50 | |
|---|
| 51 | logical,save :: firstcall=.true. |
|---|
| 52 | |
|---|
| 53 | ! Tracer indexes in the GCM: |
|---|
| 54 | integer,save :: g_co2=0 |
|---|
| 55 | integer,save :: g_co=0 |
|---|
| 56 | integer,save :: g_o=0 |
|---|
| 57 | integer,save :: g_o1d=0 |
|---|
| 58 | integer,save :: g_o2=0 |
|---|
| 59 | ! integer,save :: g_o3=0 |
|---|
| 60 | integer,save :: g_h=0 |
|---|
| 61 | integer,save :: g_h2=0 |
|---|
| 62 | integer,save :: g_oh=0 |
|---|
| 63 | integer,save :: g_ho2=0 |
|---|
| 64 | integer,save :: g_h2o2=0 |
|---|
| 65 | ! integer,save :: g_n2=0 |
|---|
| 66 | integer,save :: g_h2o_vap=0 |
|---|
| 67 | ! Tracer indexes in the thermospheric chemistry: |
|---|
| 68 | integer,parameter :: i_co2=1 |
|---|
| 69 | integer,parameter :: i_o2=2 |
|---|
| 70 | integer,parameter :: i_o=3 |
|---|
| 71 | integer,parameter :: i_co=4 |
|---|
| 72 | integer,parameter :: i_h=5 |
|---|
| 73 | integer,parameter :: i_oh=6 |
|---|
| 74 | integer,parameter :: i_ho2=7 |
|---|
| 75 | integer,parameter :: i_h2=8 |
|---|
| 76 | integer,parameter :: i_h2o=9 |
|---|
| 77 | integer,parameter :: i_h2o2=10 |
|---|
| 78 | integer,parameter :: i_o1d=11 |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | ! Initializations at first call |
|---|
| 82 | if (firstcall) then |
|---|
| 83 | c call param_read |
|---|
| 84 | ! get the indexes of the tracers we'll need |
|---|
| 85 | g_co2=igcm_co2 |
|---|
| 86 | if (g_co2.eq.0) then |
|---|
| 87 | write(*,*) "chemthermos: Error; no CO2 tracer !!!" |
|---|
| 88 | stop |
|---|
| 89 | endif |
|---|
| 90 | g_co=igcm_co |
|---|
| 91 | if (g_co.eq.0) then |
|---|
| 92 | write(*,*) "chemthermos: Error; no CO tracer !!!" |
|---|
| 93 | stop |
|---|
| 94 | endif |
|---|
| 95 | g_o=igcm_o |
|---|
| 96 | if (g_o.eq.0) then |
|---|
| 97 | write(*,*) "chemthermos: Error; no O tracer !!!" |
|---|
| 98 | stop |
|---|
| 99 | endif |
|---|
| 100 | g_o1d=igcm_o1d |
|---|
| 101 | if (g_o1d.eq.0) then |
|---|
| 102 | write(*,*) "chemthermos: Error; no O1D tracer !!!" |
|---|
| 103 | stop |
|---|
| 104 | endif |
|---|
| 105 | g_o2=igcm_o2 |
|---|
| 106 | if (g_o2.eq.0) then |
|---|
| 107 | write(*,*) "chemthermos: Error; no O2 tracer !!!" |
|---|
| 108 | stop |
|---|
| 109 | endif |
|---|
| 110 | g_h=igcm_h |
|---|
| 111 | if (g_h.eq.0) then |
|---|
| 112 | write(*,*) "chemthermos: Error; no H tracer !!!" |
|---|
| 113 | stop |
|---|
| 114 | endif |
|---|
| 115 | g_h2=igcm_h2 |
|---|
| 116 | if (g_h2.eq.0) then |
|---|
| 117 | write(*,*) "chemthermos: Error; no H2 tracer !!!" |
|---|
| 118 | stop |
|---|
| 119 | endif |
|---|
| 120 | g_oh=igcm_oh |
|---|
| 121 | if (g_oh.eq.0) then |
|---|
| 122 | write(*,*) "chemthermos: Error; no OH tracer !!!" |
|---|
| 123 | stop |
|---|
| 124 | endif |
|---|
| 125 | g_ho2=igcm_ho2 |
|---|
| 126 | if (g_ho2.eq.0) then |
|---|
| 127 | write(*,*) "chemthermos: Error; no HO2 tracer !!!" |
|---|
| 128 | stop |
|---|
| 129 | endif |
|---|
| 130 | g_h2o2=igcm_h2o2 |
|---|
| 131 | if (g_h2o2.eq.0) then |
|---|
| 132 | write(*,*) "chemthermos: Error; no H2O2 tracer !!!" |
|---|
| 133 | stop |
|---|
| 134 | endif |
|---|
| 135 | g_h2o_vap=igcm_h2o_vap |
|---|
| 136 | if (g_h2o_vap.eq.0) then |
|---|
| 137 | write(*,*) "chemthermos: Error; no water vapor tracer !!!" |
|---|
| 138 | stop |
|---|
| 139 | endif |
|---|
| 140 | |
|---|
| 141 | firstcall= .false. |
|---|
| 142 | endif ! of if (firstcall) |
|---|
| 143 | |
|---|
| 144 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 145 | |
|---|
| 146 | nlayer=nlayermx |
|---|
| 147 | c zlocal(1)=0.00625 |
|---|
| 148 | c do l=2,nlayer |
|---|
| 149 | c tmean=ztemp(l) |
|---|
| 150 | c if(ztemp(l).ne.ztemp(l-1)) |
|---|
| 151 | c & tmean=(ztemp(l)-ztemp(l-1))/log(ztemp(l)/ztemp(l-1)) |
|---|
| 152 | c zlocal(l)= zlocal(l-1)-log(zpress(l)/zpress(l-1)) |
|---|
| 153 | c & *Rnew(ig,l-1)*tmean/g/1000. |
|---|
| 154 | c enddo |
|---|
| 155 | |
|---|
| 156 | do l=1,nlayer |
|---|
| 157 | aux1(l)=0. |
|---|
| 158 | aux2(l)=0. |
|---|
| 159 | rm(l,i_co2) = zycol(l,g_co2) *zdens(l) |
|---|
| 160 | rm(l,i_co) = zycol(l,g_co) *zdens(l) |
|---|
| 161 | rm(l,i_o) = zycol(l,g_o) *zdens(l) |
|---|
| 162 | rm(l,i_o1d) = zycol(l,g_o1d) *zdens(l) |
|---|
| 163 | rm(l,i_o2) = zycol(l,g_o2) *zdens(l) |
|---|
| 164 | rm(l,i_h) = zycol(l,g_h) *zdens(l) |
|---|
| 165 | rm(l,i_h2) = zycol(l,g_h2) *zdens(l) |
|---|
| 166 | rm(l,i_oh) = zycol(l,g_oh) *zdens(l) |
|---|
| 167 | rm(l,i_ho2) = zycol(l,g_ho2) *zdens(l) |
|---|
| 168 | rm(l,i_h2o2) = zycol(l,g_h2o2) *zdens(l) |
|---|
| 169 | rm(l,i_h2o) = zycol(l,g_h2o_vap)*zdens(l) |
|---|
| 170 | enddo |
|---|
| 171 | |
|---|
| 172 | call flujo(solarcondate+zday/365.) |
|---|
| 173 | |
|---|
| 174 | call jthermcalc |
|---|
| 175 | $ (rm(1,i_co2),rm(1,i_o2),rm(1,i_o),rm(1,i_h2),rm(1,i_h2o), |
|---|
| 176 | & rm(1,i_h2o2),aux1,aux2,ztemp,nlayermx,zlocal, |
|---|
| 177 | & solarcondate+zday/365.,zenit) |
|---|
| 178 | |
|---|
| 179 | call paramfoto(lswitch,zdens,ztemp,ptimestep/3600.,zenit, |
|---|
| 180 | & nlayer,rm(1,i_co2),rm(1,i_o2),rm(1,i_o), |
|---|
| 181 | & rm(1,i_co),rm(1,i_h),rm(1,i_oh),rm(1,i_ho2), |
|---|
| 182 | & rm(1,i_h2),rm(1,i_h2o),rm(1,i_h2o2),rm(1,i_o1d)) |
|---|
| 183 | |
|---|
| 184 | do l=lswitch,nlayer |
|---|
| 185 | zycol(l,g_co2) = max(rm(l,i_co2) / zdens(l) , 1.e-30) |
|---|
| 186 | zycol(l,g_co) = max(rm(l,i_co) / zdens(l) , 1.e-30) |
|---|
| 187 | zycol(l,g_o2) = max(rm(l,i_o2) / zdens(l) , 1.e-30) |
|---|
| 188 | zycol(l,g_h2) = max(rm(l,i_h2) / zdens(l) , 1.e-30) |
|---|
| 189 | zycol(l,g_h) = max(rm(l,i_h) / zdens(l) , 1.e-30) |
|---|
| 190 | zycol(l,g_oh) = max(rm(l,i_oh) / zdens(l) , 1.e-30) |
|---|
| 191 | zycol(l,g_ho2) = max(rm(l,i_ho2) / zdens(l) , 1.e-30) |
|---|
| 192 | zycol(l,g_h2o_vap) = max(rm(l,i_h2o) / zdens(l) , 1.e-30) |
|---|
| 193 | zycol(l,g_h2o2) = max(rm(l,i_h2o2) / zdens(l) , 1.e-30) |
|---|
| 194 | zycol(l,g_o1d) = max(rm(l,i_o1d) / zdens(l) , 1.e-30) |
|---|
| 195 | zycol(l,g_o) = max(rm(l,i_o) / zdens(l) , 1.e-30) |
|---|
| 196 | enddo !nlayer |
|---|
| 197 | |
|---|
| 198 | |
|---|
| 199 | return |
|---|
| 200 | end |
|---|
