[253] | 1 | Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, & |
---|
| 2 | aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky) |
---|
[135] | 3 | |
---|
[726] | 4 | use radinc_h, only : L_TAUMAX,naerkind |
---|
| 5 | use aerosol_mod |
---|
[787] | 6 | USE comgeomfi_h |
---|
[858] | 7 | USE tracer_h, only: noms,rho_co2,rho_ice |
---|
[135] | 8 | |
---|
| 9 | implicit none |
---|
| 10 | |
---|
| 11 | !================================================================== |
---|
| 12 | ! |
---|
| 13 | ! Purpose |
---|
| 14 | ! ------- |
---|
| 15 | ! Compute aerosol optical depth in each gridbox. |
---|
| 16 | ! |
---|
| 17 | ! Authors |
---|
| 18 | ! ------- |
---|
| 19 | ! F. Forget |
---|
| 20 | ! F. Montmessin (water ice scheme) |
---|
| 21 | ! update J.-B. Madeleine (2008) |
---|
| 22 | ! dust removal, simplification by Robin Wordsworth (2009) |
---|
| 23 | ! |
---|
| 24 | ! Input |
---|
| 25 | ! ----- |
---|
| 26 | ! ngrid Number of horizontal gridpoints |
---|
| 27 | ! nlayer Number of layers |
---|
| 28 | ! nq Number of tracers |
---|
| 29 | ! pplev Pressure (Pa) at each layer boundary |
---|
| 30 | ! pq Aerosol mixing ratio |
---|
| 31 | ! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius |
---|
[787] | 32 | ! QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients |
---|
| 33 | ! QREFir3d(ngrid,nlayermx,naerkind) / at reference wavelengths |
---|
[135] | 34 | ! |
---|
| 35 | ! Output |
---|
| 36 | ! ------ |
---|
| 37 | ! aerosol Aerosol optical depth in layer l, grid point ig |
---|
| 38 | ! tau_col Total column optical depth at grid point ig |
---|
| 39 | ! |
---|
| 40 | !======================================================================= |
---|
| 41 | |
---|
| 42 | #include "dimensions.h" |
---|
| 43 | #include "dimphys.h" |
---|
| 44 | #include "callkeys.h" |
---|
| 45 | #include "comcstfi.h" |
---|
[726] | 46 | #include "comvert.h" |
---|
[135] | 47 | |
---|
[858] | 48 | INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns |
---|
| 49 | INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers |
---|
| 50 | INTEGER,INTENT(IN) :: nq ! number of tracers |
---|
| 51 | REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) |
---|
| 52 | REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) |
---|
| 53 | REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air) |
---|
| 54 | REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth |
---|
| 55 | REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius |
---|
| 56 | REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayermx,naerkind) ! extinction coefficient in the visible |
---|
| 57 | REAL,INTENT(IN) :: QREFir3d(ngrid,nlayermx,naerkind) |
---|
| 58 | REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth |
---|
| 59 | ! BENJAMIN MODIFS |
---|
| 60 | real,intent(in) :: cloudfrac(ngrid,nlayermx) ! cloud fraction |
---|
| 61 | real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction |
---|
| 62 | logical,intent(in) :: clearsky |
---|
[135] | 63 | |
---|
[253] | 64 | real aerosol0 |
---|
| 65 | |
---|
[135] | 66 | INTEGER l,ig,iq,iaer |
---|
| 67 | |
---|
[858] | 68 | LOGICAL,SAVE :: firstcall=.true. |
---|
[135] | 69 | |
---|
| 70 | REAL CBRT |
---|
| 71 | EXTERNAL CBRT |
---|
| 72 | |
---|
| 73 | INTEGER,SAVE :: i_co2ice=0 ! co2 ice |
---|
| 74 | INTEGER,SAVE :: i_h2oice=0 ! water ice |
---|
| 75 | CHARACTER(LEN=20) :: tracername ! to temporarily store text |
---|
| 76 | |
---|
[253] | 77 | ! for fixed dust profiles |
---|
| 78 | real topdust, expfactor, zp |
---|
[787] | 79 | REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling |
---|
| 80 | REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling |
---|
[858] | 81 | |
---|
[728] | 82 | real CLFtot |
---|
[253] | 83 | |
---|
| 84 | ! identify tracers |
---|
[135] | 85 | IF (firstcall) THEN |
---|
| 86 | |
---|
[726] | 87 | write(*,*) "Tracers found in aeropacity:" |
---|
[787] | 88 | do iq=1,nq |
---|
[135] | 89 | tracername=noms(iq) |
---|
| 90 | if (tracername.eq."co2_ice") then |
---|
| 91 | i_co2ice=iq |
---|
[726] | 92 | write(*,*) "i_co2ice=",i_co2ice |
---|
| 93 | |
---|
[135] | 94 | endif |
---|
| 95 | if (tracername.eq."h2o_ice") then |
---|
| 96 | i_h2oice=iq |
---|
[726] | 97 | write(*,*) "i_h2oice=",i_h2oice |
---|
[135] | 98 | endif |
---|
| 99 | enddo |
---|
| 100 | |
---|
[741] | 101 | if (noaero) then |
---|
| 102 | print*, "No active aerosols found in aeropacity" |
---|
| 103 | else |
---|
[726] | 104 | print*, "If you would like to use aerosols, make sure any old" |
---|
| 105 | print*, "start files are updated in newstart using the option" |
---|
| 106 | print*, "q=0" |
---|
[741] | 107 | write(*,*) "Active aerosols found in aeropacity:" |
---|
[726] | 108 | endif |
---|
[253] | 109 | |
---|
[741] | 110 | if ((iaero_co2.ne.0).and.(.not.noaero)) then |
---|
[726] | 111 | print*, 'iaero_co2= ',iaero_co2 |
---|
| 112 | endif |
---|
| 113 | if (iaero_h2o.ne.0) then |
---|
| 114 | print*,'iaero_h2o= ',iaero_h2o |
---|
| 115 | endif |
---|
| 116 | if (iaero_dust.ne.0) then |
---|
| 117 | print*,'iaero_dust= ',iaero_dust |
---|
| 118 | endif |
---|
| 119 | if (iaero_h2so4.ne.0) then |
---|
| 120 | print*,'iaero_h2so4= ',iaero_h2so4 |
---|
| 121 | endif |
---|
| 122 | |
---|
[135] | 123 | firstcall=.false. |
---|
| 124 | ENDIF ! of IF (firstcall) |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | ! --------------------------------------------------------- |
---|
| 128 | !================================================================== |
---|
[726] | 129 | ! CO2 ice aerosols |
---|
[135] | 130 | !================================================================== |
---|
| 131 | |
---|
[728] | 132 | if (iaero_co2.ne.0) then |
---|
[726] | 133 | iaer=iaero_co2 |
---|
[135] | 134 | ! 1. Initialization |
---|
[787] | 135 | aerosol(1:ngrid,1:nlayermx,iaer)=0.0 |
---|
[135] | 136 | ! 2. Opacity calculation |
---|
[804] | 137 | if (noaero) then ! aerosol set to zero |
---|
| 138 | aerosol(1:ngrid,1:nlayermx,iaer)=0.0 |
---|
[741] | 139 | elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed |
---|
[787] | 140 | aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9 |
---|
| 141 | !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option |
---|
[253] | 142 | else |
---|
| 143 | DO ig=1, ngrid |
---|
| 144 | DO l=1,nlayer-1 ! to stop the rad tran bug |
---|
[135] | 145 | |
---|
[253] | 146 | aerosol0 = & |
---|
| 147 | ( 0.75 * QREFvis3d(ig,l,iaer) / & |
---|
| 148 | ( rho_co2 * reffrad(ig,l,iaer) ) ) * & |
---|
| 149 | ( pq(ig,l,i_co2ice) + 1.E-9 ) * & |
---|
| 150 | ( pplev(ig,l) - pplev(ig,l+1) ) / g |
---|
| 151 | aerosol0 = max(aerosol0,1.e-9) |
---|
| 152 | aerosol0 = min(aerosol0,L_TAUMAX) |
---|
| 153 | aerosol(ig,l,iaer) = aerosol0 |
---|
| 154 | ! aerosol(ig,l,iaer) = 0.0 |
---|
[726] | 155 | ! print*, aerosol(ig,l,iaer) |
---|
[253] | 156 | ! using cloud fraction |
---|
| 157 | ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) |
---|
| 158 | ! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX) |
---|
| 159 | |
---|
| 160 | |
---|
| 161 | ENDDO |
---|
| 162 | ENDDO |
---|
[726] | 163 | end if ! if fixed or varying |
---|
[728] | 164 | end if ! if CO2 aerosols |
---|
[135] | 165 | !================================================================== |
---|
[726] | 166 | ! Water ice / liquid |
---|
[135] | 167 | !================================================================== |
---|
| 168 | |
---|
[728] | 169 | if (iaero_h2o.ne.0) then |
---|
[726] | 170 | iaer=iaero_h2o |
---|
[135] | 171 | ! 1. Initialization |
---|
[787] | 172 | aerosol(1:ngrid,1:nlayermx,iaer)=0.0 |
---|
[135] | 173 | ! 2. Opacity calculation |
---|
[726] | 174 | if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then |
---|
[787] | 175 | aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9 |
---|
[305] | 176 | |
---|
| 177 | ! put cloud at cloudlvl |
---|
| 178 | if(kastprof.and.(cloudlvl.ne.0.0))then |
---|
| 179 | ig=1 |
---|
| 180 | do l=1,nlayer |
---|
| 181 | if(int(cloudlvl).eq.l)then |
---|
| 182 | !if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then |
---|
| 183 | print*,'Inserting cloud at level ',l |
---|
| 184 | !aerosol(ig,l,iaer)=10.0 |
---|
| 185 | |
---|
| 186 | rho_ice=920.0 |
---|
| 187 | |
---|
| 188 | ! the Kasting approximation |
---|
| 189 | aerosol(ig,l,iaer) = & |
---|
| 190 | ( 0.75 * QREFvis3d(ig,l,iaer) / & |
---|
| 191 | ( rho_ice * reffrad(ig,l,iaer) ) ) * & |
---|
| 192 | !( pq(ig,l,i_h2oice) + 1.E-9 ) * & |
---|
| 193 | ( 4.0e-4 + 1.E-9 ) * & |
---|
| 194 | ( pplev(ig,l) - pplev(ig,l+1) ) / g |
---|
| 195 | |
---|
| 196 | |
---|
| 197 | open(115,file='clouds.out',form='formatted') |
---|
| 198 | write(115,*) l,aerosol(ig,l,iaer) |
---|
| 199 | close(115) |
---|
| 200 | |
---|
| 201 | return |
---|
| 202 | endif |
---|
| 203 | end do |
---|
| 204 | |
---|
| 205 | call abort |
---|
| 206 | endif |
---|
| 207 | |
---|
[253] | 208 | else |
---|
[728] | 209 | |
---|
[253] | 210 | do ig=1, ngrid |
---|
| 211 | do l=1,nlayer-1 ! to stop the rad tran bug |
---|
[135] | 212 | |
---|
[728] | 213 | aerosol(ig,l,iaer) = & !modification by BC |
---|
[253] | 214 | ( 0.75 * QREFvis3d(ig,l,iaer) / & |
---|
| 215 | ( rho_ice * reffrad(ig,l,iaer) ) ) * & |
---|
[716] | 216 | ! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same |
---|
[728] | 217 | !( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the |
---|
| 218 | ! clear=false/pq=0 case |
---|
[716] | 219 | ( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW) |
---|
[728] | 220 | ( pplev(ig,l) - pplev(ig,l+1) ) / g |
---|
[135] | 221 | |
---|
[253] | 222 | enddo |
---|
| 223 | enddo |
---|
[135] | 224 | |
---|
[728] | 225 | if(CLFvarying)then |
---|
[787] | 226 | call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer)) |
---|
[728] | 227 | do ig=1, ngrid |
---|
| 228 | do l=1,nlayer-1 ! to stop the rad tran bug |
---|
| 229 | CLFtot = max(totcloudfrac(ig),0.01) |
---|
| 230 | aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot |
---|
| 231 | aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) |
---|
| 232 | enddo |
---|
| 233 | enddo |
---|
| 234 | else |
---|
| 235 | do ig=1, ngrid |
---|
| 236 | do l=1,nlayer-1 ! to stop the rad tran bug |
---|
| 237 | CLFtot = CLFfixval |
---|
| 238 | aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot |
---|
| 239 | aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9) |
---|
| 240 | enddo |
---|
| 241 | enddo |
---|
| 242 | end if!(CLFvarying) |
---|
| 243 | endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky) |
---|
| 244 | |
---|
| 245 | end if ! End if h2o aerosol |
---|
| 246 | |
---|
[726] | 247 | !================================================================== |
---|
| 248 | ! Dust |
---|
| 249 | !================================================================== |
---|
[728] | 250 | if (iaero_dust.ne.0) then |
---|
[726] | 251 | iaer=iaero_dust |
---|
| 252 | ! 1. Initialization |
---|
[787] | 253 | aerosol(1:ngrid,1:nlayermx,iaer)=0.0 |
---|
[726] | 254 | |
---|
[728] | 255 | topdust=30.0 ! km (used to be 10.0 km) LK |
---|
[726] | 256 | |
---|
| 257 | ! 2. Opacity calculation |
---|
| 258 | |
---|
| 259 | ! expfactor=0. |
---|
| 260 | DO l=1,nlayer-1 |
---|
| 261 | DO ig=1,ngrid |
---|
| 262 | ! Typical mixing ratio profile |
---|
| 263 | |
---|
| 264 | zp=(preff/pplay(ig,l))**(70./topdust) |
---|
| 265 | expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) |
---|
| 266 | |
---|
| 267 | ! Vertical scaling function |
---|
| 268 | aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) & |
---|
| 269 | *expfactor |
---|
| 270 | |
---|
| 271 | |
---|
| 272 | ENDDO |
---|
| 273 | ENDDO |
---|
| 274 | |
---|
| 275 | ! Rescaling each layer to reproduce the choosen (or assimilated) |
---|
| 276 | ! dust extinction opacity at visible reference wavelength, which |
---|
| 277 | ! is scaled to the "preff" reference surface pressure available in comvert.h |
---|
| 278 | ! and stored in startfi.nc |
---|
| 279 | |
---|
| 280 | taudusttmp(1:ngrid)=0. |
---|
| 281 | DO l=1,nlayer |
---|
| 282 | DO ig=1,ngrid |
---|
| 283 | taudusttmp(ig) = taudusttmp(ig) & |
---|
| 284 | + aerosol(ig,l,iaer) |
---|
| 285 | ENDDO |
---|
| 286 | ENDDO |
---|
| 287 | DO l=1,nlayer-1 |
---|
| 288 | DO ig=1,ngrid |
---|
| 289 | aerosol(ig,l,iaer) = max(1E-20, & |
---|
| 290 | dusttau & |
---|
| 291 | * pplev(ig,1) / preff & |
---|
| 292 | * aerosol(ig,l,iaer) & |
---|
| 293 | / taudusttmp(ig)) |
---|
| 294 | |
---|
| 295 | ENDDO |
---|
| 296 | ENDDO |
---|
[728] | 297 | end if ! If dust aerosol |
---|
[726] | 298 | |
---|
[135] | 299 | !================================================================== |
---|
[726] | 300 | ! H2SO4 |
---|
[253] | 301 | !================================================================== |
---|
[726] | 302 | ! added by LK |
---|
| 303 | if (iaero_h2so4.ne.0) then |
---|
[728] | 304 | iaer=iaero_h2so4 |
---|
[135] | 305 | |
---|
[253] | 306 | ! 1. Initialization |
---|
[787] | 307 | aerosol(1:ngrid,1:nlayermx,iaer)=0.0 |
---|
[253] | 308 | |
---|
| 309 | |
---|
| 310 | ! 2. Opacity calculation |
---|
| 311 | |
---|
[726] | 312 | ! expfactor=0. |
---|
[728] | 313 | DO l=1,nlayer-1 |
---|
| 314 | DO ig=1,ngrid |
---|
| 315 | ! Typical mixing ratio profile |
---|
[253] | 316 | |
---|
[728] | 317 | zp=(preff/pplay(ig,l))**(70./30) !emulating topdust |
---|
| 318 | expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) |
---|
[135] | 319 | |
---|
[726] | 320 | ! Vertical scaling function |
---|
[728] | 321 | aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor |
---|
[135] | 322 | |
---|
[728] | 323 | ENDDO |
---|
| 324 | ENDDO |
---|
| 325 | tauh2so4tmp(1:ngrid)=0. |
---|
| 326 | DO l=1,nlayer |
---|
| 327 | DO ig=1,ngrid |
---|
| 328 | tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer) |
---|
| 329 | ENDDO |
---|
| 330 | ENDDO |
---|
| 331 | DO l=1,nlayer-1 |
---|
| 332 | DO ig=1,ngrid |
---|
| 333 | aerosol(ig,l,iaer) = max(1E-20, & |
---|
[726] | 334 | 1 & |
---|
| 335 | * pplev(ig,1) / preff & |
---|
| 336 | * aerosol(ig,l,iaer) & |
---|
| 337 | / tauh2so4tmp(ig)) |
---|
| 338 | |
---|
| 339 | ENDDO |
---|
[728] | 340 | ENDDO |
---|
[726] | 341 | |
---|
| 342 | ! 1/700. is assuming a "sulfurtau" of 1 |
---|
| 343 | ! Sulfur aerosol routine to be improved. |
---|
| 344 | ! aerosol0 = & |
---|
| 345 | ! ( 0.75 * QREFvis3d(ig,l,iaer) / & |
---|
| 346 | ! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * & |
---|
| 347 | ! ( pq(ig,l,i_h2so4) + 1.E-9 ) * & |
---|
| 348 | ! ( pplev(ig,l) - pplev(ig,l+1) ) / g |
---|
| 349 | ! aerosol0 = max(aerosol0,1.e-9) |
---|
| 350 | ! aerosol0 = min(aerosol0,L_TAUMAX) |
---|
| 351 | ! aerosol(ig,l,iaer) = aerosol0 |
---|
| 352 | |
---|
| 353 | ! ENDDO |
---|
| 354 | ! ENDDO |
---|
| 355 | end if |
---|
| 356 | |
---|
| 357 | |
---|
| 358 | |
---|
[135] | 359 | ! -------------------------------------------------------------------------- |
---|
| 360 | ! Column integrated visible optical depth in each point (used for diagnostic) |
---|
| 361 | |
---|
[253] | 362 | tau_col(:)=0.0 |
---|
[135] | 363 | do iaer = 1, naerkind |
---|
| 364 | do l=1,nlayer |
---|
| 365 | do ig=1,ngrid |
---|
| 366 | tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) |
---|
| 367 | end do |
---|
| 368 | end do |
---|
| 369 | end do |
---|
| 370 | |
---|
[787] | 371 | do ig=1,ngrid |
---|
[253] | 372 | do l=1,nlayer |
---|
| 373 | do iaer = 1, naerkind |
---|
| 374 | if(aerosol(ig,l,iaer).gt.1.e3)then |
---|
| 375 | print*,'WARNING: aerosol=',aerosol(ig,l,iaer) |
---|
| 376 | print*,'at ig=',ig,', l=',l,', iaer=',iaer |
---|
| 377 | print*,'QREFvis3d=',QREFvis3d(ig,l,iaer) |
---|
| 378 | print*,'reffrad=',reffrad(ig,l,iaer) |
---|
| 379 | endif |
---|
| 380 | end do |
---|
| 381 | end do |
---|
| 382 | end do |
---|
| 383 | |
---|
[787] | 384 | do ig=1,ngrid |
---|
[253] | 385 | if(tau_col(ig).gt.1.e3)then |
---|
| 386 | print*,'WARNING: tau_col=',tau_col(ig) |
---|
| 387 | print*,'at ig=',ig |
---|
| 388 | print*,'aerosol=',aerosol(ig,:,:) |
---|
| 389 | print*,'QREFvis3d=',QREFvis3d(ig,:,:) |
---|
| 390 | print*,'reffrad=',reffrad(ig,:,:) |
---|
| 391 | endif |
---|
| 392 | end do |
---|
[135] | 393 | return |
---|
| 394 | end subroutine aeropacity |
---|
| 395 | |
---|