| 1 | subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,eps23,eps33,eps76) |
|---|
| 2 | |
|---|
| 3 | use comgeomfi_h |
|---|
| 4 | use callkeys_mod, only: strobel |
|---|
| 5 | use tracer_h, only: igcm_n2,mmol |
|---|
| 6 | use comcstfi_mod, only: pi, g, cpp,r |
|---|
| 7 | implicit none |
|---|
| 8 | |
|---|
| 9 | !================================================================== |
|---|
| 10 | ! Purpose |
|---|
| 11 | ! ------- |
|---|
| 12 | ! Calculation of epsilon - non LTE efficiency factor |
|---|
| 13 | ! For CH4 heating and cooling rate calculation |
|---|
| 14 | ! Cf Strobel et al 1996 |
|---|
| 15 | ! |
|---|
| 16 | ! Inputs |
|---|
| 17 | ! ------ |
|---|
| 18 | ! ngrid Number of vertical columns |
|---|
| 19 | ! nlayer Number of layers |
|---|
| 20 | ! pplay(ngrid,nlayer) Pressure layers |
|---|
| 21 | ! pt |
|---|
| 22 | ! |
|---|
| 23 | ! Outputs |
|---|
| 24 | ! ------- |
|---|
| 25 | ! |
|---|
| 26 | ! eps_sw ! efficiency factor for heating rate |
|---|
| 27 | ! eps_lw ! efficiency factor for cooling rate |
|---|
| 28 | ! |
|---|
| 29 | ! Authors |
|---|
| 30 | ! ------- |
|---|
| 31 | ! Tanguy Bertrand (2015,2022) |
|---|
| 32 | ! |
|---|
| 33 | !================================================================== |
|---|
| 34 | |
|---|
| 35 | !----------------------------------------------------------------------- |
|---|
| 36 | |
|---|
| 37 | ! Arguments |
|---|
| 38 | |
|---|
| 39 | INTEGER ngrid, nlayer, nq |
|---|
| 40 | REAL pplay(ngrid,nlayer) |
|---|
| 41 | REAL pplev(ngrid,nlayer) |
|---|
| 42 | REAL pt(ngrid,nlayer) |
|---|
| 43 | REAL vmrch4(ngrid,nlayer) |
|---|
| 44 | |
|---|
| 45 | !----------------------------------------------------------------------- |
|---|
| 46 | ! Local variables |
|---|
| 47 | |
|---|
| 48 | INTEGER l,ig |
|---|
| 49 | |
|---|
| 50 | REAL cbol ! Boltzman constant |
|---|
| 51 | REAL avog ! avogadro constant |
|---|
| 52 | REAL beta23 |
|---|
| 53 | REAL beta33 |
|---|
| 54 | REAL nd(ngrid,nlayer) ! density number cm-3 |
|---|
| 55 | REAL na(ngrid,nlayer) ! column density number cm-2 |
|---|
| 56 | REAL rho(ngrid,nlayer) ! density of atmosphere |
|---|
| 57 | REAL kn3,k3,k2n4,k3n4,k4 |
|---|
| 58 | REAL An3,A3,A2n4,A3n4,A4 |
|---|
| 59 | REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4 |
|---|
| 60 | REAL phin3(ngrid,nlayer) |
|---|
| 61 | REAL phi3(ngrid,nlayer) |
|---|
| 62 | REAL phi2n4(ngrid,nlayer) |
|---|
| 63 | REAL phi3n4(ngrid,nlayer) |
|---|
| 64 | REAL phi4(ngrid,nlayer) |
|---|
| 65 | REAL phi4T(ngrid,nlayer) |
|---|
| 66 | |
|---|
| 67 | ! DATA from Strobel et al 1996 |
|---|
| 68 | REAL, DIMENSION(22) :: gam2 ! band escape fonction 2.3 |
|---|
| 69 | REAL, DIMENSION(22) :: Naval2 ! column density number 2.3 |
|---|
| 70 | REAL, DIMENSION(22) :: gam3 ! band escape fonction 3.3 |
|---|
| 71 | REAL, DIMENSION(22) :: Naval3 ! column density number 3.3 |
|---|
| 72 | REAL gammafi23(nlayer),gammafi33(nlayer) ! Band espace fonction |
|---|
| 73 | REAL gammafi23mod(nlayer),gammafi33mod(nlayer) ! Band espace fonction |
|---|
| 74 | REAL vmr |
|---|
| 75 | REAL k4zT ! function to calculate k4 depending on T and vmr |
|---|
| 76 | REAL log_square ! log function |
|---|
| 77 | REAL log_affine ! log function |
|---|
| 78 | REAL ratio_affine ! linear approx function to calculate ratio for gamma calculation |
|---|
| 79 | REAL ratio_square ! |
|---|
| 80 | REAL ratio_sq,ratio_af |
|---|
| 81 | |
|---|
| 82 | ! Output |
|---|
| 83 | REAL,intent(out) :: eps23(ngrid,nlayer) |
|---|
| 84 | REAL,intent(out) :: eps33(ngrid,nlayer) |
|---|
| 85 | REAL,intent(out) :: eps76(ngrid,nlayer) |
|---|
| 86 | |
|---|
| 87 | ! Option simplified coefficient |
|---|
| 88 | REAL klw, ksw, eps_sw1Pa |
|---|
| 89 | |
|---|
| 90 | !----------------------------------------------------------------------- |
|---|
| 91 | eps23(:,:)=1. |
|---|
| 92 | eps33(:,:)=1. |
|---|
| 93 | eps76(:,:)=1. |
|---|
| 94 | |
|---|
| 95 | IF (strobel) then |
|---|
| 96 | |
|---|
| 97 | ! Constantes |
|---|
| 98 | cbol=1.3806488e-23 !Boltzman constant |
|---|
| 99 | avog=6.022141e23 |
|---|
| 100 | |
|---|
| 101 | ! band 3 |
|---|
| 102 | k3=2.8e-11 !cm-3 |
|---|
| 103 | A3=25.2 |
|---|
| 104 | dphi3=k3/A3 |
|---|
| 105 | ! band n3 |
|---|
| 106 | kn3=2.8e-11 !cm-3 |
|---|
| 107 | An3=27.3 |
|---|
| 108 | dphin3=kn3/An3 |
|---|
| 109 | ! band 3n4 |
|---|
| 110 | k3n4=1.16e-11 !cm-3 |
|---|
| 111 | A3n4=2.12 |
|---|
| 112 | dphi3n4=k3n4/(3.*A3n4) |
|---|
| 113 | ! band 4 |
|---|
| 114 | !! Here we use function k4zT instead of constant k4 |
|---|
| 115 | k4=4.5e-16 !cm-3 |
|---|
| 116 | A4=2.12 |
|---|
| 117 | dphi4=k4/A4 |
|---|
| 118 | ! band 2n4 |
|---|
| 119 | k2n4=1.16e-11 |
|---|
| 120 | A2n4=2.12 |
|---|
| 121 | dphi2n4=k2n4/(2.*A2n4) |
|---|
| 122 | |
|---|
| 123 | ! 2.3 band : |
|---|
| 124 | beta23=0.09 |
|---|
| 125 | ! 3.3 band : |
|---|
| 126 | beta33=0.13 |
|---|
| 127 | |
|---|
| 128 | ! Escape Function Gamma (Data Strobel et al) : |
|---|
| 129 | ! 100 K 2.3 |
|---|
| 130 | Naval2=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, & |
|---|
| 131 | 7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, & |
|---|
| 132 | 9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, & |
|---|
| 133 | 1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/) |
|---|
| 134 | gam2=(/1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,9.90000000e-01, & |
|---|
| 135 | 9.80000000e-01,9.50000000e-0,9.00000000e-01,8.00000000e-01,7.00000000e-01,5.00000000e-01, & |
|---|
| 136 | 3.50000000e-01,2.20000000e-01,1.20000000e-01,6.50000000e-02,3.20000000e-02,1.40000000e-02, & |
|---|
| 137 | 6.20000000e-03,2.50000000e-03,1.05000000e-03,4.83227951e-06/) |
|---|
| 138 | |
|---|
| 139 | ! 100 K 3.3 |
|---|
| 140 | Naval3=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, & |
|---|
| 141 | 7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, & |
|---|
| 142 | 9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, & |
|---|
| 143 | 1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/) |
|---|
| 144 | gam3=(/1.00000000e+00,9.80000000e-01,9.00000000e-01,8.50000000e-01,7.00000000e-01,5.00000000e-01, & |
|---|
| 145 | 3.00000000e-01,1.50000000e-01,1.00000000e-01,6.00000000e-02,4.00000000e-02,2.30000000e-02, & |
|---|
| 146 | 1.45000000e-02,9.00000000e-03,5.00000000e-03,3.00000000e-03,1.50000000e-03,8.00000000e-04, & |
|---|
| 147 | 3.80000000e-04,1.80000000e-04,8.30000000e-05,6.81737544e-07/) |
|---|
| 148 | |
|---|
| 149 | DO ig=1,ngrid |
|---|
| 150 | |
|---|
| 151 | phi3n4(ig,:)=1. |
|---|
| 152 | phi2n4(ig,:)=1. |
|---|
| 153 | phin3(ig,:)=1. |
|---|
| 154 | phi4(ig,:)=1. |
|---|
| 155 | nd(ig,:)=1. |
|---|
| 156 | na(ig,:)=1. |
|---|
| 157 | |
|---|
| 158 | !write(*,*) 'TB22 pplev=',pplev(1,:) |
|---|
| 159 | !write(*,*) 'TB22 tlev=',pt(1,:) |
|---|
| 160 | |
|---|
| 161 | DO l=1,nlayer |
|---|
| 162 | !!calculation of density |
|---|
| 163 | rho(ig,l)=pplay(ig,l)/(r*pt(ig,l)) |
|---|
| 164 | |
|---|
| 165 | !!calculation of column number density for ch4 (cm-2) |
|---|
| 166 | na(ig,l)=vmrch4(ig,l)/100./(mmol(igcm_n2)*1.e-3)* & |
|---|
| 167 | pplev(ig,l)/g*avog*1.e-4 |
|---|
| 168 | |
|---|
| 169 | !!calculation of atmospheric number density (cm-3) for each layer |
|---|
| 170 | nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)*1.e-3)*avog*1.e-6 |
|---|
| 171 | |
|---|
| 172 | !!calculation of the phi : |
|---|
| 173 | phin3(ig,l)=dphin3*nd(ig,l) |
|---|
| 174 | phi3(ig,l)=dphi3*nd(ig,l) |
|---|
| 175 | phi2n4(ig,l)=dphi2n4*nd(ig,l) |
|---|
| 176 | phi3n4(ig,l)=dphi3n4*nd(ig,l) |
|---|
| 177 | phi4(ig,l)=dphi4*nd(ig,l) |
|---|
| 178 | !! Getting phi4 depending on T and vmr |
|---|
| 179 | vmr=vmrch4(ig,l)/100. |
|---|
| 180 | phi4T(ig,l)=k4zT(pt(ig,l),vmr)/A4*nd(ig,l) |
|---|
| 181 | ENDDO |
|---|
| 182 | !write(*,*) 'TB22 nd=',nd(1,:) |
|---|
| 183 | !write(*,*) 'TB22 vmr=',vmrch4 |
|---|
| 184 | !write(*,*) 'TB22 phi4T=',phi4T |
|---|
| 185 | |
|---|
| 186 | !! interpolation of na values in cm-2 from Strobel values for |
|---|
| 187 | !the 100K model: |
|---|
| 188 | CALL interp_line(Naval2,gam2,22,na(ig,:),gammafi23,nlayer) |
|---|
| 189 | CALL interp_line(Naval3,gam3,22,na(ig,:),gammafi33,nlayer) |
|---|
| 190 | |
|---|
| 191 | !write(*,*) 'TB22rad na = ',na(1,:) |
|---|
| 192 | !write(*,*) 'TB22rad gammafi = ',gammafi23 |
|---|
| 193 | |
|---|
| 194 | DO l=1,nlayer |
|---|
| 195 | |
|---|
| 196 | !! Getting ratio depending on temperature and n (cm-2) |
|---|
| 197 | ratio_sq=ratio_square(pt(ig,l),na(ig,l)) |
|---|
| 198 | ratio_af=ratio_affine(pt(ig,l),na(ig,l)) |
|---|
| 199 | !! Modified gamma depending on T |
|---|
| 200 | gammafi23mod(l)=min(ratio_sq*gammafi23(l),1.) |
|---|
| 201 | gammafi33mod(l)=min(ratio_af*gammafi33(l),1.) |
|---|
| 202 | |
|---|
| 203 | !! Heating : cf Strobel 1996 and Zalucha 2011 |
|---|
| 204 | ! eps 2.3 |
|---|
| 205 | eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23mod(l))* & |
|---|
| 206 | (beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)* & |
|---|
| 207 | phi4T(ig,l)/(phi4T(ig,l)+gammafi23mod(l))* & |
|---|
| 208 | (3.*phi2n4(ig,l)+1.)/(1.+phi2n4(ig,l))*(1-beta23)/3.) |
|---|
| 209 | |
|---|
| 210 | ! eps 3.3 |
|---|
| 211 | eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33mod(l))* & |
|---|
| 212 | (beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)* & |
|---|
| 213 | phi4T(ig,l)/(phi4T(ig,l)+gammafi33mod(l))* & |
|---|
| 214 | (1-beta33)) |
|---|
| 215 | |
|---|
| 216 | ! Cooling from Zalucha et al 2013 eq10 |
|---|
| 217 | eps76(ig,l)=1./(1.+gammafi23mod(l)/(2.*phi4T(ig,l))) |
|---|
| 218 | |
|---|
| 219 | ENDDO |
|---|
| 220 | !write(*,*) 'TB22rad eps23 = ',eps23(1,:) |
|---|
| 221 | !write(*,*) 'TB22rad eps33 = ',eps33(1,:) |
|---|
| 222 | !write(*,*) 'TB22rad gamma33 = ',gammafi33mod |
|---|
| 223 | !write(*,*) 'TB22rad eps76 = ',eps76(1,:) |
|---|
| 224 | |
|---|
| 225 | ENDDO ! ngrid |
|---|
| 226 | |
|---|
| 227 | ELSE ! if (strobel) |
|---|
| 228 | !if strobel = false |
|---|
| 229 | |
|---|
| 230 | klw = 3.6E-4 ! fit strobel |
|---|
| 231 | !eps_sw1Pa = 0.6 ! fit strobel |
|---|
| 232 | eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5 |
|---|
| 233 | !eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25 |
|---|
| 234 | |
|---|
| 235 | ksw = 8.5E-5*(0.9/(eps_sw1Pa-0.1)-1.) |
|---|
| 236 | do l=1, nlayer |
|---|
| 237 | do ig=1, ngrid |
|---|
| 238 | rho(ig,l)=pplay(ig,l)/(r*pt(ig,l)) |
|---|
| 239 | eps23(ig,l) = 0.1+ 0.9/(1.+ksw/rho(ig,l)) |
|---|
| 240 | eps76(ig,l) = 1./(1.+klw/rho(ig,l)) |
|---|
| 241 | end do |
|---|
| 242 | end do |
|---|
| 243 | |
|---|
| 244 | END IF |
|---|
| 245 | |
|---|
| 246 | end subroutine nlte_ch4 |
|---|
| 247 | |
|---|
| 248 | real function k4zT(temp,vmr) |
|---|
| 249 | implicit none |
|---|
| 250 | real, intent(in) :: temp,vmr |
|---|
| 251 | k4zT=1.9e-15*exp((temp-240.)/70.)*(1+10.*vmr) |
|---|
| 252 | return |
|---|
| 253 | end function k4zT |
|---|
| 254 | |
|---|
| 255 | !! Log log affine and square function for ratio adjustment and extrapolation |
|---|
| 256 | |
|---|
| 257 | real function log_square(x,y1,y2,x1,x2,coef) |
|---|
| 258 | implicit none |
|---|
| 259 | real, intent(in) :: x,y1,y2,x1,x2,coef |
|---|
| 260 | real aa,bb |
|---|
| 261 | aa=(log(y2)-log(y1))/((log(coef*x2))**4-(log(coef*x1)**4)) |
|---|
| 262 | bb=log(y2)-aa*(log(coef*x2))**4 |
|---|
| 263 | log_square=exp(aa*(log(coef*x))**4+bb) |
|---|
| 264 | return |
|---|
| 265 | end function log_square |
|---|
| 266 | |
|---|
| 267 | real function log_affine(x,y1,y2,x1,x2) |
|---|
| 268 | real, intent(in) :: x,y1,y2,x1,x2 |
|---|
| 269 | real aa,bb |
|---|
| 270 | aa=(log(y2)-log(y1))/(log(x2)-log(x1)) |
|---|
| 271 | bb=log(y2)-aa*log(x2) |
|---|
| 272 | log_affine=exp(aa*log(x)+bb) |
|---|
| 273 | return |
|---|
| 274 | end function log_affine |
|---|
| 275 | |
|---|
| 276 | !! Ratio of gamma depending on T. Linear approximation |
|---|
| 277 | real function ratio_affine(temp,n) |
|---|
| 278 | implicit none |
|---|
| 279 | real, intent(in) :: temp,n |
|---|
| 280 | real aa,bb,ratio0 |
|---|
| 281 | real log_affine |
|---|
| 282 | aa=(1.-0.2)/(100.-40.) |
|---|
| 283 | bb=1.-aa*100. |
|---|
| 284 | ratio0=aa*temp+bb |
|---|
| 285 | ratio_affine=log_affine(n,1.,ratio0,1e15,1e22) |
|---|
| 286 | return |
|---|
| 287 | end function ratio_affine |
|---|
| 288 | |
|---|
| 289 | !! Ratio of gamma depending on T. Linear approximation |
|---|
| 290 | real function ratio_square(temp,n) |
|---|
| 291 | implicit none |
|---|
| 292 | real, intent(in) :: temp,n |
|---|
| 293 | real aa,bb,ratio0,coef |
|---|
| 294 | real log_square |
|---|
| 295 | coef=1.e-12 |
|---|
| 296 | aa=(1.-0.2)/(100.-40.) |
|---|
| 297 | bb=1.-aa*100. |
|---|
| 298 | ratio0=aa*temp+bb |
|---|
| 299 | ratio_square=log_square(n,1.,ratio0,1e15,1e22,coef) |
|---|
| 300 | return |
|---|
| 301 | end function ratio_square |
|---|
| 302 | |
|---|
| 303 | |
|---|
| 304 | |
|---|
