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