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