| 1 | SUBROUTINE nirco2abs(nlon,nlev,nplay,dist_sol,nq,pq, |
|---|
| 2 | $ mu0,fract,pdtnirco2) |
|---|
| 3 | |
|---|
| 4 | use dimphy |
|---|
| 5 | use geometry_mod, only: longitude_deg, latitude_deg |
|---|
| 6 | use chemparam_mod, only: i_co2, i_o |
|---|
| 7 | c use compo_hedin83_mod2 |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | c======================================================================= |
|---|
| 12 | c subject: |
|---|
| 13 | c -------- |
|---|
| 14 | c Computing heating rate due to |
|---|
| 15 | c absorption by CO2 in the near-infrared |
|---|
| 16 | c This version includes NLTE effects |
|---|
| 17 | c |
|---|
| 18 | c (Scheme to be described in Forget et al., JGR, 2003) |
|---|
| 19 | c (old Scheme described in Forget et al., JGR, 1999) |
|---|
| 20 | c |
|---|
| 21 | c This version updated with a new functional fit, |
|---|
| 22 | c see NLTE correction-factor of Lopez-Valverde et al (1998) |
|---|
| 23 | c Stephen Lewis 2000 |
|---|
| 24 | c |
|---|
| 25 | c apr 2019 d.quirino Improving NLTE params, SOIR/SPICAV Temp comparison |
|---|
| 26 | c oct 2014 g.gilli Coupling with photochemical model |
|---|
| 27 | C jan 2014 g.gilli Revision (following martian non-lte param) |
|---|
| 28 | C jun 2013 l.salmi First adaptation to Venus and NIR NLTE param |
|---|
| 29 | c jul 2011 malv+fgg New corrections for NLTE implemented |
|---|
| 30 | c 08/2002 : correction for bug when running with diurnal=F |
|---|
| 31 | c |
|---|
| 32 | c author: Frederic Hourdin 1996 |
|---|
| 33 | c ------ |
|---|
| 34 | c Francois Forget 1999 |
|---|
| 35 | c |
|---|
| 36 | c input: |
|---|
| 37 | c ----- |
|---|
| 38 | c nlon number of gridpoint of horizontal grid |
|---|
| 39 | c nlev Number of layer |
|---|
| 40 | c dist_sol sun-Venus distance (AU) |
|---|
| 41 | c mu0(nlon) |
|---|
| 42 | c fract(nlon) day fraction of the time interval |
|---|
| 43 | c declin latitude of subslar point |
|---|
| 44 | c |
|---|
| 45 | c output: |
|---|
| 46 | c ------- |
|---|
| 47 | c |
|---|
| 48 | c pdtnirco2(nlon,nlev) Heating rate (K/sec) |
|---|
| 49 | c |
|---|
| 50 | c |
|---|
| 51 | c======================================================================= |
|---|
| 52 | c |
|---|
| 53 | c 0. Declarations : |
|---|
| 54 | c ------------------ |
|---|
| 55 | c |
|---|
| 56 | |
|---|
| 57 | #include "YOMCST.h" |
|---|
| 58 | #include "clesphys.h" |
|---|
| 59 | c#include "comdiurn.h" |
|---|
| 60 | #include "nirdata.h" |
|---|
| 61 | c#include "tracer.h" |
|---|
| 62 | #include "mmol.h" |
|---|
| 63 | c----------------------------------------------------------------------- |
|---|
| 64 | c Input/Output |
|---|
| 65 | c ------------ |
|---|
| 66 | integer,intent(in) :: nlon ! number of (horizontal) grid points |
|---|
| 67 | integer,intent(in) :: nlev ! number of atmospheric layers |
|---|
| 68 | |
|---|
| 69 | real,intent(in) :: nplay(nlon,nlev) ! Pressure |
|---|
| 70 | real,intent(in) :: dist_sol ! Sun-Venus distance (in AU) |
|---|
| 71 | integer,intent(in) :: nq ! number of tracers |
|---|
| 72 | real,intent(in) :: pq(nlon,nlev,nq) ! mass mixing ratio tracers |
|---|
| 73 | real,intent(in) :: mu0(nlon) ! solar angle |
|---|
| 74 | real,intent(in) :: fract(nlon) ! day fraction of the time interval |
|---|
| 75 | c real,intent(in) :: declin ! latitude of sub-solar point |
|---|
| 76 | real :: co2vmr_gcm(nlon,nlev), o3pvmr_gcm(nlon,nlev) |
|---|
| 77 | |
|---|
| 78 | real,intent(out) :: pdtnirco2(nlon,nlev) ! heating rate (K/sec) |
|---|
| 79 | |
|---|
| 80 | c |
|---|
| 81 | c Local variables : |
|---|
| 82 | c ----------------- |
|---|
| 83 | INTEGER l,ig, n, nstep,i |
|---|
| 84 | REAL co2heat0, zmu(nlon) |
|---|
| 85 | |
|---|
| 86 | c special diurnal=F |
|---|
| 87 | real mu0_int(nlon),fract_int(nlon),zday_int |
|---|
| 88 | real ztim1,ztim2,ztim3,step |
|---|
| 89 | |
|---|
| 90 | logical onepeak |
|---|
| 91 | parameter (onepeak=.false.) |
|---|
| 92 | c parameter (onepeak=.true.) |
|---|
| 93 | c |
|---|
| 94 | c local saved variables |
|---|
| 95 | c --------------------- |
|---|
| 96 | logical,save :: firstcall=.true. |
|---|
| 97 | integer,save :: ico2=0 ! index of "co2" tracer |
|---|
| 98 | integer,save :: io=0 ! index of "o" tracer |
|---|
| 99 | |
|---|
| 100 | ccc================================================= |
|---|
| 101 | cccc parameters for CO2 heating fit |
|---|
| 102 | ccc================================================= |
|---|
| 103 | |
|---|
| 104 | c-------------------------------------------------- |
|---|
| 105 | c One-peak martian-type fit => Gabriella (2014+) |
|---|
| 106 | c-------------------------------------------------- |
|---|
| 107 | c n_a = heating rate for Venusian day at p0, r0, mu =0 [K day-1] |
|---|
| 108 | c Here p0 = p_cloud top [Pa] |
|---|
| 109 | c n_p0 = is a pressure below which non LTE effects are significant [Pa] |
|---|
| 110 | c n_a Solar heating [K/Eday] at the cloud top, taken from Crisps table |
|---|
| 111 | |
|---|
| 112 | real n_a, n_p0, n_b, p_ctop |
|---|
| 113 | |
|---|
| 114 | cc "Nominal" values used in Gilli+2017 |
|---|
| 115 | c parameter (n_a = 18.13/86400.0) !c K/Eday ---> K/sec |
|---|
| 116 | c parameter (p_ctop=13.2e2) |
|---|
| 117 | c parameter (n_p0=0.008) |
|---|
| 118 | |
|---|
| 119 | cc "New" values used to improve SPICAV/SOIR Temperature comparision (D.Quirino) |
|---|
| 120 | cc Gilli+2021 |
|---|
| 121 | parameter (n_a = 15.92/86400.0) !c K/Eday ---> K/sec |
|---|
| 122 | parameter (p_ctop=19.85e2) |
|---|
| 123 | parameter (n_p0=0.1) |
|---|
| 124 | parameter (n_b=1.362) |
|---|
| 125 | |
|---|
| 126 | c -- NLTE Param v2 -- |
|---|
| 127 | C parameter (n_p0=0.01) |
|---|
| 128 | c parameter (n_b = 1.3) |
|---|
| 129 | |
|---|
| 130 | c-------------------------------------------------- |
|---|
| 131 | c Multi-peaks Roldan-type fit => Laura (2013) |
|---|
| 132 | c New paramaters (Param9*0.5) => Enora (2021) |
|---|
| 133 | c-------------------------------------------------- |
|---|
| 134 | c ENORA FINE TUNING used for VCD 1.1 |
|---|
| 135 | c (fit to fig 12 Roldan-2000) |
|---|
| 136 | real n_coFB, n_aFB, n_bFB, n_p0FB, n_eFB |
|---|
| 137 | real n_coISO, n_aISO, n_bISO, n_p0ISO, n_eISO |
|---|
| 138 | real n_coFH, n_aFH, n_bFH, n_p0FH, n_eFH |
|---|
| 139 | real n_co43, n_a43, n_b43, n_p043, n_e43 |
|---|
| 140 | real n_co43b, n_a43b, n_b43b, n_p043b, n_e43b |
|---|
| 141 | real n_conir, n_anir, n_bnir, n_p0nir, n_enir |
|---|
| 142 | |
|---|
| 143 | parameter (n_coFB=119./86400.0) !c K/Eday ---> K/sec |
|---|
| 144 | parameter (n_aFB=0.185) |
|---|
| 145 | parameter (n_bFB=3.7) |
|---|
| 146 | parameter (n_p0FB=2.9e-4) |
|---|
| 147 | parameter (n_eFB=0.76) |
|---|
| 148 | |
|---|
| 149 | parameter (n_coISO=265./86400.0) !c K/Eday ---> K/sec |
|---|
| 150 | parameter (n_aISO=0.313) |
|---|
| 151 | parameter (n_bISO=1.65) |
|---|
| 152 | parameter (n_p0ISO=0.076) |
|---|
| 153 | parameter (n_eISO=0.99) |
|---|
| 154 | |
|---|
| 155 | parameter (n_coFH=2.5/86400.0) !c K/Eday ---> K/sec |
|---|
| 156 | parameter (n_aFH=3.98) |
|---|
| 157 | parameter (n_bFH=2.9) |
|---|
| 158 | parameter (n_p0FH=0.17) |
|---|
| 159 | parameter (n_eFH=2.16) |
|---|
| 160 | |
|---|
| 161 | parameter (n_co43=55./86400.0) !c K/Eday ---> K/sec |
|---|
| 162 | parameter (n_a43=0.625) |
|---|
| 163 | parameter (n_b43=2.6) |
|---|
| 164 | parameter (n_p043=0.043) |
|---|
| 165 | parameter (n_e43=1.654) |
|---|
| 166 | |
|---|
| 167 | ! parameter (n_co43b=100./86400.0) !c K/Eday ---> K/sec |
|---|
| 168 | ! => fine tuning: not affected by the *0.5 below (see ENORA FINE TUNING) |
|---|
| 169 | parameter (n_co43b=200./86400.0) !c K/Eday ---> K/sec |
|---|
| 170 | parameter (n_a43b=5.5) |
|---|
| 171 | parameter (n_b43b=2.3) |
|---|
| 172 | parameter (n_p043b=1.) |
|---|
| 173 | parameter (n_e43b=0.4) |
|---|
| 174 | |
|---|
| 175 | parameter (n_conir=6.5/86400.0) !c K/Eday ---> K/sec |
|---|
| 176 | parameter (n_anir=35.65) |
|---|
| 177 | parameter (n_bnir=2.1) |
|---|
| 178 | parameter (n_p0nir=0.046) |
|---|
| 179 | parameter (n_enir=0.9) |
|---|
| 180 | |
|---|
| 181 | real :: picFB(nlon,nlev), picISO(nlon,nlev), picFH(nlon,nlev) |
|---|
| 182 | real :: pic43(nlon,nlev), pic43b(nlon,nlev), picnir(nlon,nlev) |
|---|
| 183 | |
|---|
| 184 | ccc================================================= |
|---|
| 185 | |
|---|
| 186 | c Variables added to implement NLTE correction factor (feb 2011) |
|---|
| 187 | real pyy(nlev) |
|---|
| 188 | real cor1(nlev),oldoco2(nlev),alfa2(nlev) |
|---|
| 189 | real p2011,cociente1,merge |
|---|
| 190 | real cor0,oco2gcm |
|---|
| 191 | |
|---|
| 192 | c---------------------------------------------------------------------- |
|---|
| 193 | c Initialisation |
|---|
| 194 | c -------------- |
|---|
| 195 | if (firstcall) then |
|---|
| 196 | if (nircorr.eq.1) then |
|---|
| 197 | c ! we will need co2 and o tracers |
|---|
| 198 | ico2= i_co2 |
|---|
| 199 | if (ico2==0) then |
|---|
| 200 | write(*,*) "nirco2abs error: I need a CO2 tracer" |
|---|
| 201 | write(*,*) " when running with nircorr==1" |
|---|
| 202 | stop |
|---|
| 203 | endif |
|---|
| 204 | io=i_o |
|---|
| 205 | if (io==0) then |
|---|
| 206 | write(*,*) "nirco2abs error: I need an O tracer" |
|---|
| 207 | write(*,*) " when running with nircorr==1" |
|---|
| 208 | stop |
|---|
| 209 | endif |
|---|
| 210 | endif |
|---|
| 211 | firstcall=.false. |
|---|
| 212 | endif |
|---|
| 213 | c -------------- |
|---|
| 214 | c co2heat0 is correction for dist_sol (is 1 for Venus) |
|---|
| 215 | co2heat0=(0.7233/dist_sol)**2 |
|---|
| 216 | |
|---|
| 217 | pdtnirco2(:,:)=0. |
|---|
| 218 | c---------------------------------------------------------------------- |
|---|
| 219 | |
|---|
| 220 | c |
|---|
| 221 | c Simple calcul for a given sun incident angle (if cycle_diurne=T) |
|---|
| 222 | c -------------------------------------------- |
|---|
| 223 | |
|---|
| 224 | IF (cycle_diurne) THEN |
|---|
| 225 | |
|---|
| 226 | do ig=1,nlon |
|---|
| 227 | zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35. |
|---|
| 228 | |
|---|
| 229 | c--------------------------- |
|---|
| 230 | if (onepeak) then |
|---|
| 231 | c--------------------------- |
|---|
| 232 | if(nircorr.eq.1) then |
|---|
| 233 | do l=1,nlev |
|---|
| 234 | pyy(l)=nplay(ig,l) |
|---|
| 235 | enddo |
|---|
| 236 | call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres) |
|---|
| 237 | call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres) |
|---|
| 238 | call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres) |
|---|
| 239 | endif |
|---|
| 240 | do l=1,nlev |
|---|
| 241 | c Calculations for the O/CO2 correction |
|---|
| 242 | if(nircorr.eq.1) then |
|---|
| 243 | cor0=1./(1.+n_p0/nplay(ig,l))**n_b |
|---|
| 244 | if(pq(ig,l,ico2) .gt. 1.e-6) then |
|---|
| 245 | oco2gcm=pq(ig,l,io)/pq(ig,l,ico2) |
|---|
| 246 | ! handle the rare cases when pq(ig,l,io)<0 |
|---|
| 247 | if (pq(ig,l,io).lt.0) then |
|---|
| 248 | write(*,*) "nirco2abs: warning ig=",ig," l=",l, |
|---|
| 249 | & " pq(ig,l,io)=",pq(ig,l,io) |
|---|
| 250 | oco2gcm=1.e6 |
|---|
| 251 | endif |
|---|
| 252 | else |
|---|
| 253 | oco2gcm=1.e6 |
|---|
| 254 | endif |
|---|
| 255 | cociente1=oco2gcm/oldoco2(l) |
|---|
| 256 | |
|---|
| 257 | c WRITE(*,*) "nirco2abs line 211", l, cociente1 |
|---|
| 258 | |
|---|
| 259 | merge=alog10(cociente1)*alfa2(l)+alog10(cor0)* |
|---|
| 260 | $ (1.-alfa2(l)) |
|---|
| 261 | merge=10**merge |
|---|
| 262 | p2011=sqrt(merge)*cor0 |
|---|
| 263 | |
|---|
| 264 | else if (nircorr.eq.0) then |
|---|
| 265 | p2011=1. |
|---|
| 266 | cor1(l)=1. |
|---|
| 267 | endif |
|---|
| 268 | |
|---|
| 269 | if(fract(ig).gt.0.) pdtnirco2(ig,l)= |
|---|
| 270 | & co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l)) |
|---|
| 271 | & /(1.+n_p0/nplay(ig,l))**n_b |
|---|
| 272 | c Corrections from tabulation |
|---|
| 273 | $ * cor1(l) * p2011 |
|---|
| 274 | |
|---|
| 275 | enddo !nlev |
|---|
| 276 | c--------------------------- |
|---|
| 277 | else ! multipeak |
|---|
| 278 | c--------------------------- |
|---|
| 279 | do l=1,nlev |
|---|
| 280 | if(fract(ig).gt.0.) then |
|---|
| 281 | picFB(ig,l)=n_coFB |
|---|
| 282 | & *((n_aFB/nplay(ig,l))**n_eFB) |
|---|
| 283 | & *zmu(ig)**0.82 |
|---|
| 284 | & /(1.+n_p0FB/nplay(ig,l))**n_bFB |
|---|
| 285 | |
|---|
| 286 | picISO(ig,l)=n_coISO |
|---|
| 287 | & *((n_aISO/nplay(ig,l))**n_eISO) |
|---|
| 288 | & *zmu(ig)**0.55 |
|---|
| 289 | & /(1.+n_p0ISO/nplay(ig,l))**n_bISO |
|---|
| 290 | |
|---|
| 291 | picFH(ig,l)=n_coFH |
|---|
| 292 | & *((n_aFH/nplay(ig,l))**n_eFH) |
|---|
| 293 | & *zmu(ig)**0.55 |
|---|
| 294 | & /(1.+n_p0FH/nplay(ig,l))**n_bFH |
|---|
| 295 | |
|---|
| 296 | pic43(ig,l)=n_co43 |
|---|
| 297 | & *((n_a43/nplay(ig,l))**n_e43) |
|---|
| 298 | & *zmu(ig)**0.55 |
|---|
| 299 | & /(1.+n_p043/nplay(ig,l))**n_b43 |
|---|
| 300 | |
|---|
| 301 | pic43b(ig,l)=n_co43b |
|---|
| 302 | & *((n_a43b/nplay(ig,l))**n_e43b) |
|---|
| 303 | & *zmu(ig)**0.55 |
|---|
| 304 | & /(1.+n_p043b/nplay(ig,l))**n_b43b |
|---|
| 305 | |
|---|
| 306 | picnir(ig,l)=n_conir |
|---|
| 307 | & *((n_anir/nplay(ig,l))**n_enir) |
|---|
| 308 | & *zmu(ig)**0.55 |
|---|
| 309 | & /(1.+n_p0nir/nplay(ig,l))**n_bnir |
|---|
| 310 | |
|---|
| 311 | pdtnirco2(ig,l)=co2heat0* |
|---|
| 312 | & (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l) |
|---|
| 313 | & +pic43b(ig,l)+picnir(ig,l))*0.5 ! *0.5 = ENORA FINE TUNING |
|---|
| 314 | |
|---|
| 315 | endif |
|---|
| 316 | enddo !nlev |
|---|
| 317 | c--------------------------- |
|---|
| 318 | endif |
|---|
| 319 | c--------------------------- |
|---|
| 320 | enddo !nlon |
|---|
| 321 | |
|---|
| 322 | |
|---|
| 323 | c Averaging over diurnal cycle (if diurnal=F) |
|---|
| 324 | c ------------------------------------------- |
|---|
| 325 | c NIR CO2 abs is slightly non linear. To remove the diurnal |
|---|
| 326 | c cycle, it is better to average the heating rate over 1 day rather |
|---|
| 327 | c than using the mean mu0 computed by mucorr in physiq.F (FF, 1998) |
|---|
| 328 | |
|---|
| 329 | ELSE ! if (.not.diurnal) then |
|---|
| 330 | nstep = 20 ! number of integration step /sol |
|---|
| 331 | do n=1,nstep |
|---|
| 332 | |
|---|
| 333 | zday_int = (n-1)/float(nstep) |
|---|
| 334 | |
|---|
| 335 | CALL zenang(0.,zday_int,RDAY/nstep, |
|---|
| 336 | & latitude_deg,longitude_deg, |
|---|
| 337 | & mu0_int,fract_int) |
|---|
| 338 | |
|---|
| 339 | do ig=1,nlon |
|---|
| 340 | zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35. |
|---|
| 341 | |
|---|
| 342 | c--------------------------- |
|---|
| 343 | if (onepeak) then |
|---|
| 344 | c--------------------------- |
|---|
| 345 | if(nircorr.eq.1) then |
|---|
| 346 | do l=1,nlev |
|---|
| 347 | pyy(l)=nplay(ig,l) |
|---|
| 348 | enddo |
|---|
| 349 | call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres) |
|---|
| 350 | call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres) |
|---|
| 351 | call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres) |
|---|
| 352 | endif |
|---|
| 353 | do l=1,nlev |
|---|
| 354 | c Calculations for the O/CO2 correction |
|---|
| 355 | if(nircorr.eq.1) then |
|---|
| 356 | cor0=1./(1.+n_p0/nplay(ig,l))**n_b |
|---|
| 357 | if(pq(ig,l,ico2) .gt. 1.e-6) then |
|---|
| 358 | oco2gcm=pq(ig,l,io)/pq(ig,l,ico2) |
|---|
| 359 | ! handle the rare cases when pq(ig,l,io)<0 |
|---|
| 360 | if (pq(ig,l,io).lt.0) then |
|---|
| 361 | write(*,*) "nirco2abs: warning ig=",ig," l=",l, |
|---|
| 362 | & " pq(ig,l,io)=",pq(ig,l,io) |
|---|
| 363 | oco2gcm=1.e6 |
|---|
| 364 | endif |
|---|
| 365 | else |
|---|
| 366 | oco2gcm=1.e6 |
|---|
| 367 | endif |
|---|
| 368 | cociente1=oco2gcm/oldoco2(l) |
|---|
| 369 | |
|---|
| 370 | c WRITE(*,*) "nirco2abs line 211", l, cociente1 |
|---|
| 371 | |
|---|
| 372 | merge=alog10(cociente1)*alfa2(l)+alog10(cor0)* |
|---|
| 373 | $ (1.-alfa2(l)) |
|---|
| 374 | merge=10**merge |
|---|
| 375 | p2011=sqrt(merge)*cor0 |
|---|
| 376 | |
|---|
| 377 | else if (nircorr.eq.0) then |
|---|
| 378 | p2011=1. |
|---|
| 379 | cor1(l)=1. |
|---|
| 380 | endif |
|---|
| 381 | |
|---|
| 382 | if(fract(ig).gt.0.) pdtnirco2(ig,l)= |
|---|
| 383 | & pdtnirco2(ig,l) + (1/float(nstep))* |
|---|
| 384 | & co2heat0*n_a*sqrt((p_ctop*zmu(ig))/nplay(ig,l)) |
|---|
| 385 | & /(1.+n_p0/nplay(ig,l))**n_b |
|---|
| 386 | c Corrections from tabulation |
|---|
| 387 | $ * cor1(l) * p2011 |
|---|
| 388 | |
|---|
| 389 | enddo !nlev |
|---|
| 390 | c--------------------------- |
|---|
| 391 | else ! multipeak |
|---|
| 392 | c--------------------------- |
|---|
| 393 | do l=1,nlev |
|---|
| 394 | if(fract(ig).gt.0.) then |
|---|
| 395 | picFB(ig,l)=n_coFB |
|---|
| 396 | & *((n_aFB/nplay(ig,l))**n_eFB) |
|---|
| 397 | & *zmu(ig)**0.82 |
|---|
| 398 | & /(1.+n_p0FB/nplay(ig,l))**n_bFB |
|---|
| 399 | |
|---|
| 400 | picISO(ig,l)=n_coISO |
|---|
| 401 | & *((n_aISO/nplay(ig,l))**n_eISO) |
|---|
| 402 | & *zmu(ig)**0.55 |
|---|
| 403 | & /(1.+n_p0ISO/nplay(ig,l))**n_bISO |
|---|
| 404 | |
|---|
| 405 | picFH(ig,l)=n_coFH |
|---|
| 406 | & *((n_aFH/nplay(ig,l))**n_eFH) |
|---|
| 407 | & *zmu(ig)**0.55 |
|---|
| 408 | & /(1.+n_p0FH/nplay(ig,l))**n_bFH |
|---|
| 409 | |
|---|
| 410 | pic43(ig,l)=n_co43 |
|---|
| 411 | & *((n_a43/nplay(ig,l))**n_e43) |
|---|
| 412 | & *zmu(ig)**0.55 |
|---|
| 413 | & /(1.+n_p043/nplay(ig,l))**n_b43 |
|---|
| 414 | |
|---|
| 415 | pic43b(ig,l)=n_co43b |
|---|
| 416 | & *((n_a43b/nplay(ig,l))**n_e43b) |
|---|
| 417 | & *zmu(ig)**0.55 |
|---|
| 418 | & /(1.+n_p043b/nplay(ig,l))**n_b43b |
|---|
| 419 | |
|---|
| 420 | picnir(ig,l)=n_conir |
|---|
| 421 | & *((n_anir/nplay(ig,l))**n_enir) |
|---|
| 422 | & *zmu(ig)**0.55 |
|---|
| 423 | & /(1.+n_p0nir/nplay(ig,l))**n_bnir |
|---|
| 424 | |
|---|
| 425 | pdtnirco2(ig,l)= |
|---|
| 426 | & pdtnirco2(ig,l)+(1/float(nstep))*co2heat0* |
|---|
| 427 | & (picFB(ig,l)+picISO(ig,l)+picFH(ig,l)+pic43(ig,l) |
|---|
| 428 | & +pic43b(ig,l)+picnir(ig,l))*0.5 ! *0.5 = ENORA FINE TUNING |
|---|
| 429 | |
|---|
| 430 | endif |
|---|
| 431 | enddo !nlev |
|---|
| 432 | c--------------------------- |
|---|
| 433 | endif |
|---|
| 434 | c--------------------------- |
|---|
| 435 | enddo !nlon |
|---|
| 436 | enddo !nstep |
|---|
| 437 | |
|---|
| 438 | END IF ! diurnal cycle |
|---|
| 439 | |
|---|
| 440 | return |
|---|
| 441 | end |
|---|
| 442 | |
|---|
| 443 | |
|---|
| 444 | subroutine interpnir(escout,p,nlev,escin,pin,nl) |
|---|
| 445 | C |
|---|
| 446 | C subroutine to perform linear interpolation in pressure from 1D profile |
|---|
| 447 | C escin(nl) sampled on pressure grid pin(nl) to profile |
|---|
| 448 | C escout(nlev) on pressure grid p(nlev). |
|---|
| 449 | C |
|---|
| 450 | real escout(nlev),p(nlev) |
|---|
| 451 | real escin(nl),pin(nl),wm,wp |
|---|
| 452 | integer nl,nlev,n1,n,nm,np |
|---|
| 453 | do n1=1,nlev |
|---|
| 454 | if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then |
|---|
| 455 | c escout(n1) = 0.0 |
|---|
| 456 | escout(n1) = 1.e-15 |
|---|
| 457 | else |
|---|
| 458 | do n = 1,nl-1 |
|---|
| 459 | if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then |
|---|
| 460 | nm=n |
|---|
| 461 | np=n+1 |
|---|
| 462 | wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np)) |
|---|
| 463 | wp=1.0 - wm |
|---|
| 464 | endif |
|---|
| 465 | enddo |
|---|
| 466 | escout(n1) = escin(nm)*wm + escin(np)*wp |
|---|
| 467 | endif |
|---|
| 468 | enddo |
|---|
| 469 | return |
|---|
| 470 | end |
|---|
