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