| 1 | SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad, |
|---|
| 2 | & QVISsQREF3d,omegaVIS3d,gVIS3d, |
|---|
| 3 | & QIRsQREF3d,omegaIR3d,gIR3d, |
|---|
| 4 | & QREFvis3d,QREFir3d, |
|---|
| 5 | & omegaREFvis3d,omegaREFir3d) |
|---|
| 6 | use dimradmars_mod, only: nir, nsun, naerkind, |
|---|
| 7 | & radiustab, nsize, QVISsQREF, omegavis, gvis, |
|---|
| 8 | & QIRsQREF, omegaIR, gIR, QREFvis, QREFir, |
|---|
| 9 | & omegaREFvis, omegaREFir |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | c ============================================================= |
|---|
| 12 | c Aerosol Optical Properties |
|---|
| 13 | c |
|---|
| 14 | c Description: |
|---|
| 15 | c Compute the scattering parameters in each grid |
|---|
| 16 | c box, depending on aerosol grain sizes. Log-normal size |
|---|
| 17 | c distribution and Gauss-Legendre integration are used. |
|---|
| 18 | |
|---|
| 19 | c Parameters: |
|---|
| 20 | c Don't forget to set the value of varyingnueff below; If |
|---|
| 21 | c the effective variance of the distribution for the given |
|---|
| 22 | c aerosol is considered homogeneous in the atmosphere, please |
|---|
| 23 | c set varyingnueff(iaer) to .false. Resulting computational |
|---|
| 24 | c time will be much better. |
|---|
| 25 | |
|---|
| 26 | c Authors: J.-B. Madeleine, F. Forget, F. Montmessin |
|---|
| 27 | c ============================================================= |
|---|
| 28 | |
|---|
| 29 | include "callkeys.h" |
|---|
| 30 | |
|---|
| 31 | c Local variables |
|---|
| 32 | c --------------- |
|---|
| 33 | |
|---|
| 34 | c ============================================================= |
|---|
| 35 | LOGICAL,SAVE,ALLOCATABLE :: varyingnueff(:) |
|---|
| 36 | c ============================================================= |
|---|
| 37 | |
|---|
| 38 | c Min. and max radius of the interpolation grid (in METERS) |
|---|
| 39 | REAL, SAVE, ALLOCATABLE :: refftabmin(:,:) |
|---|
| 40 | REAL, SAVE, ALLOCATABLE :: refftabmax(:,:) |
|---|
| 41 | |
|---|
| 42 | !$OMP THREADPRIVATE(varyingnueff,refftabmin,refftabmax) |
|---|
| 43 | |
|---|
| 44 | c Log of the min and max variance of the interpolation grid |
|---|
| 45 | REAL, PARAMETER :: nuefftabmin = -4.6 |
|---|
| 46 | REAL, PARAMETER :: nuefftabmax = 0. |
|---|
| 47 | c Number of effective radius of the interpolation grid |
|---|
| 48 | INTEGER, PARAMETER :: refftabsize = 100 |
|---|
| 49 | c Number of effective variances of the interpolation grid |
|---|
| 50 | INTEGER, PARAMETER :: nuefftabsize = 100 |
|---|
| 51 | c Interpolation grid indices (reff,nueff) |
|---|
| 52 | INTEGER :: grid_i,grid_j |
|---|
| 53 | c Intermediate variable |
|---|
| 54 | REAL :: var_tmp,var3d_tmp(ngrid,nlayer) |
|---|
| 55 | c Bilinear interpolation factors |
|---|
| 56 | REAL :: kx,ky,k1,k2,k3,k4 |
|---|
| 57 | c Size distribution parameters |
|---|
| 58 | REAL :: sizedistk1,sizedistk2 |
|---|
| 59 | c Pi! |
|---|
| 60 | REAL,SAVE :: pi |
|---|
| 61 | |
|---|
| 62 | !$OMP THREADPRIVATE(pi) |
|---|
| 63 | |
|---|
| 64 | c Variables used by the Gauss-Legendre integration: |
|---|
| 65 | INTEGER radius_id,gausind |
|---|
| 66 | REAL kint |
|---|
| 67 | REAL drad |
|---|
| 68 | INTEGER, PARAMETER :: ngau = 10 |
|---|
| 69 | REAL weightgaus(ngau),radgaus(ngau) |
|---|
| 70 | SAVE weightgaus,radgaus |
|---|
| 71 | |
|---|
| 72 | !$OMP THREADPRIVATE(weightgaus,radgaus) |
|---|
| 73 | |
|---|
| 74 | c DATA weightgaus/.2955242247,.2692667193,.2190863625, |
|---|
| 75 | c & .1494513491,.0666713443/ |
|---|
| 76 | c DATA radgaus/.1488743389,.4333953941,.6794095682, |
|---|
| 77 | c & .8650633666,.9739065285/ |
|---|
| 78 | DATA radgaus/0.07652652113350,0.22778585114165, |
|---|
| 79 | & 0.37370608871528,0.51086700195146, |
|---|
| 80 | & 0.63605368072468,0.74633190646476, |
|---|
| 81 | & 0.83911697181213,0.91223442826796, |
|---|
| 82 | & 0.96397192726078,0.99312859919241/ |
|---|
| 83 | |
|---|
| 84 | DATA weightgaus/0.15275338723120,0.14917298659407, |
|---|
| 85 | & 0.14209610937519,0.13168863843930, |
|---|
| 86 | & 0.11819453196154,0.10193011980823, |
|---|
| 87 | & 0.08327674160932,0.06267204829828, |
|---|
| 88 | & 0.04060142982019,0.01761400714091/ |
|---|
| 89 | |
|---|
| 90 | c Indices |
|---|
| 91 | INTEGER :: i,j,k,l,m,iaer,idomain |
|---|
| 92 | INTEGER :: ig,lg,chg |
|---|
| 93 | |
|---|
| 94 | c Local saved variables |
|---|
| 95 | c --------------------- |
|---|
| 96 | |
|---|
| 97 | c Radius axis of the interpolation grid |
|---|
| 98 | REAL,SAVE,ALLOCATABLE :: refftab(:,:,:) |
|---|
| 99 | c Variance axis of the interpolation grid |
|---|
| 100 | REAL,SAVE,ALLOCATABLE :: nuefftab(:,:,:) |
|---|
| 101 | c Volume ratio of the grid |
|---|
| 102 | REAL,SAVE,ALLOCATABLE :: logvratgrid(:,:) |
|---|
| 103 | |
|---|
| 104 | !$OMP THREADPRIVATE(refftab,nuefftab,logvratgrid) |
|---|
| 105 | |
|---|
| 106 | c Grid used to remember which calculation is done |
|---|
| 107 | LOGICAL,SAVE,ALLOCATABLE :: checkgrid(:,:,:,:) |
|---|
| 108 | c Optical properties of the grid (VISIBLE) |
|---|
| 109 | REAL,SAVE,ALLOCATABLE :: qsqrefVISgrid(:,:,:,:) |
|---|
| 110 | REAL,SAVE,ALLOCATABLE :: qextVISgrid(:,:,:,:) |
|---|
| 111 | REAL,SAVE,ALLOCATABLE :: qscatVISgrid(:,:,:,:) |
|---|
| 112 | REAL,SAVE,ALLOCATABLE :: omegVISgrid(:,:,:,:) |
|---|
| 113 | REAL,SAVE,ALLOCATABLE :: gVISgrid(:,:,:,:) |
|---|
| 114 | |
|---|
| 115 | !$OMP THREADPRIVATE(checkgrid,qsqrefVISgrid,qextVISgrid, |
|---|
| 116 | !$OMP& qscatVISgrid,omegVISgrid,gVISgrid) |
|---|
| 117 | |
|---|
| 118 | |
|---|
| 119 | c Optical properties of the grid (INFRARED) |
|---|
| 120 | REAL,SAVE,ALLOCATABLE :: qsqrefIRgrid(:,:,:,:) |
|---|
| 121 | REAL,SAVE,ALLOCATABLE :: qextIRgrid(:,:,:,:) |
|---|
| 122 | REAL,SAVE,ALLOCATABLE :: qscatIRgrid(:,:,:,:) |
|---|
| 123 | REAL,SAVE,ALLOCATABLE :: omegIRgrid(:,:,:,:) |
|---|
| 124 | REAL,SAVE,ALLOCATABLE :: gIRgrid(:,:,:,:) |
|---|
| 125 | |
|---|
| 126 | !$OMP THREADPRIVATE(qsqrefIRgrid,qextIRgrid,qscatIRgrid, |
|---|
| 127 | !$OMP& omegIRgrid,gIRgrid) |
|---|
| 128 | |
|---|
| 129 | c Optical properties of the grid (REFERENCE WAVELENGTHS) |
|---|
| 130 | REAL,SAVE,ALLOCATABLE :: qrefVISgrid(:,:,:) |
|---|
| 131 | REAL,SAVE,ALLOCATABLE :: qscatrefVISgrid(:,:,:) |
|---|
| 132 | REAL,SAVE,ALLOCATABLE :: qrefIRgrid(:,:,:) |
|---|
| 133 | REAL,SAVE,ALLOCATABLE :: qscatrefIRgrid(:,:,:) |
|---|
| 134 | REAL,SAVE,ALLOCATABLE :: omegrefVISgrid(:,:,:) |
|---|
| 135 | REAL,SAVE,ALLOCATABLE :: omegrefIRgrid(:,:,:) |
|---|
| 136 | c Firstcall |
|---|
| 137 | LOGICAL,SAVE :: firstcall = .true. |
|---|
| 138 | |
|---|
| 139 | !$OMP THREADPRIVATE(qrefVISgrid,qscatrefVISgrid,qrefIRgrid , |
|---|
| 140 | !$OMP& qscatrefIRgrid,omegrefVISgrid,omegrefIRgrid, |
|---|
| 141 | !$OMP& firstcall) |
|---|
| 142 | |
|---|
| 143 | c Variables used by the Gauss-Legendre integration: |
|---|
| 144 | REAL,SAVE,ALLOCATABLE :: normd(:,:,:,:) |
|---|
| 145 | REAL,SAVE,ALLOCATABLE :: dista(:,:,:,:,:) |
|---|
| 146 | REAL,SAVE,ALLOCATABLE :: distb(:,:,:,:,:) |
|---|
| 147 | REAL,SAVE,ALLOCATABLE :: radGAUSa(:,:,:) |
|---|
| 148 | REAL,SAVE,ALLOCATABLE :: radGAUSb(:,:,:) |
|---|
| 149 | |
|---|
| 150 | !$OMP THREADPRIVATE(normd,dista,distb,radGAUSa,radGAUSb) |
|---|
| 151 | |
|---|
| 152 | REAL,SAVE,ALLOCATABLE :: qsqrefVISa(:,:,:) |
|---|
| 153 | REAL,SAVE,ALLOCATABLE :: qrefVISa(:,:) |
|---|
| 154 | REAL,SAVE,ALLOCATABLE :: qsqrefVISb(:,:,:) |
|---|
| 155 | REAL,SAVE,ALLOCATABLE :: qrefVISb(:,:) |
|---|
| 156 | REAL,SAVE,ALLOCATABLE :: omegVISa(:,:,:) |
|---|
| 157 | REAL,SAVE,ALLOCATABLE :: omegrefVISa(:,:) |
|---|
| 158 | REAL,SAVE,ALLOCATABLE :: omegVISb(:,:,:) |
|---|
| 159 | REAL,SAVE,ALLOCATABLE :: omegrefVISb(:,:) |
|---|
| 160 | REAL,SAVE,ALLOCATABLE :: gVISa(:,:,:) |
|---|
| 161 | REAL,SAVE,ALLOCATABLE :: gVISb(:,:,:) |
|---|
| 162 | |
|---|
| 163 | !$OMP THREADPRIVATE(qsqrefVISa,qrefVISa,qsqrefVISb,qrefVISb, |
|---|
| 164 | !$OMP& omegVISa,omegrefVISa,omegVISb,omegrefVISb, |
|---|
| 165 | !$OMP& gVISa,gVISb) |
|---|
| 166 | |
|---|
| 167 | REAL,SAVE,ALLOCATABLE :: qsqrefIRa(:,:,:) |
|---|
| 168 | REAL,SAVE,ALLOCATABLE :: qrefIRa(:,:) |
|---|
| 169 | REAL,SAVE,ALLOCATABLE :: qsqrefIRb(:,:,:) |
|---|
| 170 | REAL,SAVE,ALLOCATABLE :: qrefIRb(:,:) |
|---|
| 171 | REAL,SAVE,ALLOCATABLE :: omegIRa(:,:,:) |
|---|
| 172 | REAL,SAVE,ALLOCATABLE :: omegrefIRa(:,:) |
|---|
| 173 | REAL,SAVE,ALLOCATABLE :: omegIRb(:,:,:) |
|---|
| 174 | REAL,SAVE,ALLOCATABLE :: omegrefIRb(:,:) |
|---|
| 175 | REAL,SAVE,ALLOCATABLE :: gIRa(:,:,:) |
|---|
| 176 | REAL,SAVE,ALLOCATABLE :: gIRb(:,:,:) |
|---|
| 177 | |
|---|
| 178 | !$OMP THREADPRIVATE(qsqrefIRa,qrefIRa,qsqrefIRb,qrefIRb, |
|---|
| 179 | !$OMP& omegIRa,omegrefIRa,omegIRb,omegrefIRb,gIRa,gIRb) |
|---|
| 180 | |
|---|
| 181 | REAL :: radiusm |
|---|
| 182 | REAL :: radiusr |
|---|
| 183 | |
|---|
| 184 | c Inputs |
|---|
| 185 | c ------ |
|---|
| 186 | |
|---|
| 187 | INTEGER,INTENT(IN) :: ngrid,nlayer |
|---|
| 188 | c Aerosol effective radius used for radiative transfer (meter) |
|---|
| 189 | REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) |
|---|
| 190 | c Aerosol effective variance used for radiative transfer (n.u.) |
|---|
| 191 | REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) |
|---|
| 192 | |
|---|
| 193 | c Outputs |
|---|
| 194 | c ------- |
|---|
| 195 | |
|---|
| 196 | REAL,INTENT(OUT) :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind) |
|---|
| 197 | REAL,INTENT(OUT) :: omegaVIS3d(ngrid,nlayer,nsun,naerkind) |
|---|
| 198 | REAL,INTENT(OUT) :: gVIS3d(ngrid,nlayer,nsun,naerkind) |
|---|
| 199 | |
|---|
| 200 | REAL,INTENT(OUT) :: QIRsQREF3d(ngrid,nlayer,nir,naerkind) |
|---|
| 201 | REAL,INTENT(OUT) :: omegaIR3d(ngrid,nlayer,nir,naerkind) |
|---|
| 202 | REAL,INTENT(OUT) :: gIR3d(ngrid,nlayer,nir,naerkind) |
|---|
| 203 | |
|---|
| 204 | REAL,INTENT(OUT) :: QREFvis3d(ngrid,nlayer,naerkind) |
|---|
| 205 | REAL,INTENT(OUT) :: QREFir3d(ngrid,nlayer,naerkind) |
|---|
| 206 | |
|---|
| 207 | REAL,INTENT(OUT) :: omegaREFvis3d(ngrid,nlayer,naerkind) |
|---|
| 208 | REAL,INTENT(OUT) :: omegaREFir3d(ngrid,nlayer,naerkind) |
|---|
| 209 | |
|---|
| 210 | c Tests |
|---|
| 211 | c ----- |
|---|
| 212 | |
|---|
| 213 | LOGICAL,SAVE :: out_qwg = .false. |
|---|
| 214 | |
|---|
| 215 | !$OMP THREADPRIVATE(out_qwg) |
|---|
| 216 | |
|---|
| 217 | INTEGER, PARAMETER :: out_iaer = 2 |
|---|
| 218 | INTEGER :: out_ndim |
|---|
| 219 | REAL :: out_qext(ngrid,nlayer) |
|---|
| 220 | REAL :: out_omeg(ngrid,nlayer) |
|---|
| 221 | REAL :: out_g(ngrid,nlayer) |
|---|
| 222 | INTEGER :: out_nchannel |
|---|
| 223 | CHARACTER*1 :: out_str |
|---|
| 224 | |
|---|
| 225 | c Creating the effective radius and variance grid |
|---|
| 226 | c ----------------------------------------------- |
|---|
| 227 | ! AS: OK firstcall absolute |
|---|
| 228 | IF (firstcall) THEN |
|---|
| 229 | |
|---|
| 230 | c 0.0 Allocate all local saved arrays: |
|---|
| 231 | allocate(refftab(refftabsize,naerkind,2)) |
|---|
| 232 | allocate(nuefftab(nuefftabsize,naerkind,2)) |
|---|
| 233 | ! Optical properties of the grid (VISIBLE) |
|---|
| 234 | allocate(qsqrefVISgrid(refftabsize,nuefftabsize,nsun,naerkind)) |
|---|
| 235 | allocate(qextVISgrid(refftabsize,nuefftabsize,nsun,naerkind)) |
|---|
| 236 | allocate(qscatVISgrid(refftabsize,nuefftabsize,nsun,naerkind)) |
|---|
| 237 | allocate(omegVISgrid(refftabsize,nuefftabsize,nsun,naerkind)) |
|---|
| 238 | allocate(gVISgrid(refftabsize,nuefftabsize,nsun,naerkind)) |
|---|
| 239 | ! Optical properties of the grid (INFRARED) |
|---|
| 240 | allocate(qsqrefIRgrid(refftabsize,nuefftabsize,nir,naerkind)) |
|---|
| 241 | allocate(qextIRgrid(refftabsize,nuefftabsize,nir,naerkind)) |
|---|
| 242 | allocate(qscatIRgrid(refftabsize,nuefftabsize,nir,naerkind)) |
|---|
| 243 | allocate(omegIRgrid(refftabsize,nuefftabsize,nir,naerkind)) |
|---|
| 244 | allocate(gIRgrid(refftabsize,nuefftabsize,nir,naerkind)) |
|---|
| 245 | |
|---|
| 246 | allocate(qsqrefVISa(nsun,ngau,naerkind)) |
|---|
| 247 | allocate(qrefVISa(ngau,naerkind)) |
|---|
| 248 | allocate(qsqrefVISb(nsun,ngau,naerkind)) |
|---|
| 249 | allocate(qrefVISb(ngau,naerkind)) |
|---|
| 250 | allocate(omegVISa(nsun,ngau,naerkind)) |
|---|
| 251 | allocate(omegrefVISa(ngau,naerkind)) |
|---|
| 252 | allocate(omegVISb(nsun,ngau,naerkind)) |
|---|
| 253 | allocate(omegrefVISb(ngau,naerkind)) |
|---|
| 254 | allocate(gVISa(nsun,ngau,naerkind)) |
|---|
| 255 | allocate(gVISb(nsun,ngau,naerkind)) |
|---|
| 256 | |
|---|
| 257 | allocate(qsqrefIRa(nir,ngau,naerkind)) |
|---|
| 258 | allocate(qrefIRa(ngau,naerkind)) |
|---|
| 259 | allocate(qsqrefIRb(nir,ngau,naerkind)) |
|---|
| 260 | allocate(qrefIRb(ngau,naerkind)) |
|---|
| 261 | allocate(omegIRa(nir,ngau,naerkind)) |
|---|
| 262 | allocate(omegrefIRa(ngau,naerkind)) |
|---|
| 263 | allocate(omegIRb(nir,ngau,naerkind)) |
|---|
| 264 | allocate(omegrefIRb(ngau,naerkind)) |
|---|
| 265 | allocate(gIRa(nir,ngau,naerkind)) |
|---|
| 266 | allocate(gIRb(nir,ngau,naerkind)) |
|---|
| 267 | |
|---|
| 268 | allocate(qrefVISgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 269 | allocate(qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 270 | allocate(qrefIRgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 271 | allocate(qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 272 | allocate(omegrefVISgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 273 | allocate(omegrefIRgrid(refftabsize,nuefftabsize,naerkind)) |
|---|
| 274 | |
|---|
| 275 | allocate(normd(refftabsize,nuefftabsize,naerkind,2)) |
|---|
| 276 | allocate(dista(refftabsize,nuefftabsize,naerkind,2,ngau)) |
|---|
| 277 | allocate(distb(refftabsize,nuefftabsize,naerkind,2,ngau)) |
|---|
| 278 | allocate(radGAUSa(ngau,naerkind,2)) |
|---|
| 279 | allocate(radGAUSb(ngau,naerkind,2)) |
|---|
| 280 | |
|---|
| 281 | allocate(checkgrid(refftabsize,nuefftabsize,naerkind,2)) |
|---|
| 282 | |
|---|
| 283 | allocate(logvratgrid(naerkind,2)) |
|---|
| 284 | |
|---|
| 285 | allocate(refftabmin(naerkind,2)) |
|---|
| 286 | allocate(refftabmax(naerkind,2)) |
|---|
| 287 | |
|---|
| 288 | allocate(varyingnueff(naerkind)) |
|---|
| 289 | |
|---|
| 290 | checkgrid(1:refftabsize,1:nuefftabsize,1:naerkind,1:2) = .false. |
|---|
| 291 | varyingnueff(1:naerkind) = .false. |
|---|
| 292 | |
|---|
| 293 | c 0.1 Pi! |
|---|
| 294 | pi = 2. * asin(1.e0) |
|---|
| 295 | |
|---|
| 296 | WRITE(*,*) "aeroptproperties: interpolation grid" |
|---|
| 297 | DO iaer = 1, naerkind ! Loop on aerosol kind |
|---|
| 298 | DO idomain = 1, 2 ! Loop on visible or infrared channel |
|---|
| 299 | |
|---|
| 300 | c 0.2 Effective radius |
|---|
| 301 | radiusm= |
|---|
| 302 | & 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))+ |
|---|
| 303 | & radiustab(iaer,idomain,1)) |
|---|
| 304 | radiusr= |
|---|
| 305 | & 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))- |
|---|
| 306 | & radiustab(iaer,idomain,1)) |
|---|
| 307 | refftabmin(iaer,idomain) = |
|---|
| 308 | & radiusm-radiusr*radgaus(ngau) |
|---|
| 309 | refftabmax(iaer,idomain) = |
|---|
| 310 | & radiustab(iaer,idomain,nsize(iaer,idomain)) |
|---|
| 311 | |
|---|
| 312 | WRITE(*,*) "Scatterer: ",iaer |
|---|
| 313 | WRITE(*,*) "Domain: ",idomain |
|---|
| 314 | WRITE(*,*) "Min radius (m): ", refftabmin(iaer,idomain) |
|---|
| 315 | WRITE(*,*) "Max radius (m): ", refftabmax(iaer,idomain) |
|---|
| 316 | |
|---|
| 317 | refftab(1,iaer,idomain) = |
|---|
| 318 | & refftabmin(iaer,idomain) |
|---|
| 319 | refftab(refftabsize,iaer,idomain) = |
|---|
| 320 | & refftabmax(iaer,idomain) |
|---|
| 321 | |
|---|
| 322 | logvratgrid(iaer,idomain) = |
|---|
| 323 | & log(refftabmax(iaer,idomain)/ |
|---|
| 324 | & refftabmin(iaer,idomain)) / |
|---|
| 325 | & float(refftabsize-1)*3. |
|---|
| 326 | do i = 2, refftabsize-1 |
|---|
| 327 | refftab(i,iaer,idomain) = refftab(i-1,iaer,idomain)* |
|---|
| 328 | & exp(1./3.*logvratgrid(iaer,idomain)) |
|---|
| 329 | enddo |
|---|
| 330 | |
|---|
| 331 | c 0.3 Effective variance |
|---|
| 332 | do i = 0, nuefftabsize-1 |
|---|
| 333 | nuefftab(i+1,iaer,idomain) = exp( nuefftabmin + |
|---|
| 334 | & i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) ) |
|---|
| 335 | enddo |
|---|
| 336 | |
|---|
| 337 | ENDDO |
|---|
| 338 | ENDDO |
|---|
| 339 | firstcall = .false. |
|---|
| 340 | ENDIF |
|---|
| 341 | |
|---|
| 342 | DO iaer = 1, naerkind ! Loop on aerosol kind |
|---|
| 343 | |
|---|
| 344 | IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN |
|---|
| 345 | c================================================================== |
|---|
| 346 | c If there is one single particle size, optical |
|---|
| 347 | c properties of the considered aerosol are homogeneous |
|---|
| 348 | DO lg = 1, nlayer |
|---|
| 349 | DO ig = 1, ngrid |
|---|
| 350 | DO chg = 1, nsun |
|---|
| 351 | QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1) |
|---|
| 352 | omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1) |
|---|
| 353 | gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1) |
|---|
| 354 | ENDDO |
|---|
| 355 | DO chg = 1, nir |
|---|
| 356 | QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1) |
|---|
| 357 | omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1) |
|---|
| 358 | gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1) |
|---|
| 359 | ENDDO |
|---|
| 360 | QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1) |
|---|
| 361 | QREFir3d(ig,lg,iaer)=QREFir(iaer,1) |
|---|
| 362 | omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1) |
|---|
| 363 | omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1) |
|---|
| 364 | ENDDO |
|---|
| 365 | ENDDO |
|---|
| 366 | ELSE ! Varying effective radius and variance |
|---|
| 367 | DO idomain = 1, 2 ! Loop on visible or infrared channel |
|---|
| 368 | c================================================================== |
|---|
| 369 | |
|---|
| 370 | c 1.1 Radius middle point and range for Gauss integration |
|---|
| 371 | radiusm= |
|---|
| 372 | & 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + |
|---|
| 373 | & radiustab(iaer,idomain,1)) |
|---|
| 374 | radiusr= |
|---|
| 375 | & 0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - |
|---|
| 376 | & radiustab(iaer,idomain,1)) |
|---|
| 377 | |
|---|
| 378 | c 1.2 Interpolating data at the Gauss quadrature points: |
|---|
| 379 | DO gausind=1,ngau |
|---|
| 380 | drad=radiusr*radgaus(gausind) |
|---|
| 381 | radGAUSa(gausind,iaer,idomain)=radiusm-drad |
|---|
| 382 | |
|---|
| 383 | radius_id=minloc(abs(radiustab(iaer,idomain,:) - |
|---|
| 384 | & (radiusm-drad)),DIM=1) |
|---|
| 385 | IF ((radiustab(iaer,idomain,radius_id) - |
|---|
| 386 | & (radiusm-drad)).GT.0) THEN |
|---|
| 387 | radius_id=radius_id-1 |
|---|
| 388 | ENDIF |
|---|
| 389 | IF (radius_id.GE.nsize(iaer,idomain)) THEN |
|---|
| 390 | radius_id=nsize(iaer,idomain)-1 |
|---|
| 391 | kint = 1. |
|---|
| 392 | ELSEIF (radius_id.LT.1) THEN |
|---|
| 393 | radius_id=1 |
|---|
| 394 | kint = 0. |
|---|
| 395 | ELSE |
|---|
| 396 | kint = ( (radiusm-drad) - |
|---|
| 397 | & radiustab(iaer,idomain,radius_id) ) / |
|---|
| 398 | & ( radiustab(iaer,idomain,radius_id+1) - |
|---|
| 399 | & radiustab(iaer,idomain,radius_id) ) |
|---|
| 400 | ENDIF |
|---|
| 401 | IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- |
|---|
| 402 | DO m=1,nsun |
|---|
| 403 | qsqrefVISa(m,gausind,iaer)= |
|---|
| 404 | & (1-kint)*QVISsQREF(m,iaer,radius_id) + |
|---|
| 405 | & kint*QVISsQREF(m,iaer,radius_id+1) |
|---|
| 406 | omegVISa(m,gausind,iaer)= |
|---|
| 407 | & (1-kint)*omegaVIS(m,iaer,radius_id) + |
|---|
| 408 | & kint*omegaVIS(m,iaer,radius_id+1) |
|---|
| 409 | gVISa(m,gausind,iaer)= |
|---|
| 410 | & (1-kint)*gVIS(m,iaer,radius_id) + |
|---|
| 411 | & kint*gVIS(m,iaer,radius_id+1) |
|---|
| 412 | ENDDO |
|---|
| 413 | qrefVISa(gausind,iaer)= |
|---|
| 414 | & (1-kint)*QREFvis(iaer,radius_id) + |
|---|
| 415 | & kint*QREFvis(iaer,radius_id+1) |
|---|
| 416 | omegrefVISa(gausind,iaer)= |
|---|
| 417 | & (1-kint)*omegaREFvis(iaer,radius_id) + |
|---|
| 418 | & kint*omegaREFvis(iaer,radius_id+1) |
|---|
| 419 | ELSE ! INFRARED DOMAIN ---------------------------------- |
|---|
| 420 | DO m=1,nir |
|---|
| 421 | qsqrefIRa(m,gausind,iaer)= |
|---|
| 422 | & (1-kint)*QIRsQREF(m,iaer,radius_id) + |
|---|
| 423 | & kint*QIRsQREF(m,iaer,radius_id+1) |
|---|
| 424 | omegIRa(m,gausind,iaer)= |
|---|
| 425 | & (1-kint)*omegaIR(m,iaer,radius_id) + |
|---|
| 426 | & kint*omegaIR(m,iaer,radius_id+1) |
|---|
| 427 | gIRa(m,gausind,iaer)= |
|---|
| 428 | & (1-kint)*gIR(m,iaer,radius_id) + |
|---|
| 429 | & kint*gIR(m,iaer,radius_id+1) |
|---|
| 430 | ENDDO |
|---|
| 431 | qrefIRa(gausind,iaer)= |
|---|
| 432 | & (1-kint)*QREFir(iaer,radius_id) + |
|---|
| 433 | & kint*QREFir(iaer,radius_id+1) |
|---|
| 434 | omegrefIRa(gausind,iaer)= |
|---|
| 435 | & (1-kint)*omegaREFir(iaer,radius_id) + |
|---|
| 436 | & kint*omegaREFir(iaer,radius_id+1) |
|---|
| 437 | ENDIF |
|---|
| 438 | ENDDO |
|---|
| 439 | |
|---|
| 440 | DO gausind=1,ngau |
|---|
| 441 | drad=radiusr*radgaus(gausind) |
|---|
| 442 | radGAUSb(gausind,iaer,idomain)=radiusm+drad |
|---|
| 443 | |
|---|
| 444 | radius_id=minloc(abs(radiustab(iaer,idomain,:) - |
|---|
| 445 | & (radiusm+drad)),DIM=1) |
|---|
| 446 | IF ((radiustab(iaer,idomain,radius_id) - |
|---|
| 447 | & (radiusm+drad)).GT.0) THEN |
|---|
| 448 | radius_id=radius_id-1 |
|---|
| 449 | ENDIF |
|---|
| 450 | IF (radius_id.GE.nsize(iaer,idomain)) THEN |
|---|
| 451 | radius_id=nsize(iaer,idomain)-1 |
|---|
| 452 | kint = 1. |
|---|
| 453 | ELSEIF (radius_id.LT.1) THEN |
|---|
| 454 | radius_id=1 |
|---|
| 455 | kint = 0. |
|---|
| 456 | ELSE |
|---|
| 457 | kint = ( (radiusm+drad) - |
|---|
| 458 | & radiustab(iaer,idomain,radius_id) ) / |
|---|
| 459 | & ( radiustab(iaer,idomain,radius_id+1) - |
|---|
| 460 | & radiustab(iaer,idomain,radius_id) ) |
|---|
| 461 | ENDIF |
|---|
| 462 | IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------------- |
|---|
| 463 | DO m=1,nsun |
|---|
| 464 | qsqrefVISb(m,gausind,iaer)= |
|---|
| 465 | & (1-kint)*QVISsQREF(m,iaer,radius_id) + |
|---|
| 466 | & kint*QVISsQREF(m,iaer,radius_id+1) |
|---|
| 467 | omegVISb(m,gausind,iaer)= |
|---|
| 468 | & (1-kint)*omegaVIS(m,iaer,radius_id) + |
|---|
| 469 | & kint*omegaVIS(m,iaer,radius_id+1) |
|---|
| 470 | gVISb(m,gausind,iaer)= |
|---|
| 471 | & (1-kint)*gVIS(m,iaer,radius_id) + |
|---|
| 472 | & kint*gVIS(m,iaer,radius_id+1) |
|---|
| 473 | ENDDO |
|---|
| 474 | qrefVISb(gausind,iaer)= |
|---|
| 475 | & (1-kint)*QREFvis(iaer,radius_id) + |
|---|
| 476 | & kint*QREFvis(iaer,radius_id+1) |
|---|
| 477 | omegrefVISb(gausind,iaer)= |
|---|
| 478 | & (1-kint)*omegaREFvis(iaer,radius_id) + |
|---|
| 479 | & kint*omegaREFvis(iaer,radius_id+1) |
|---|
| 480 | ELSE ! INFRARED DOMAIN ---------------------------------- |
|---|
| 481 | DO m=1,nir |
|---|
| 482 | qsqrefIRb(m,gausind,iaer)= |
|---|
| 483 | & (1-kint)*QIRsQREF(m,iaer,radius_id) + |
|---|
| 484 | & kint*QIRsQREF(m,iaer,radius_id+1) |
|---|
| 485 | omegIRb(m,gausind,iaer)= |
|---|
| 486 | & (1-kint)*omegaIR(m,iaer,radius_id) + |
|---|
| 487 | & kint*omegaIR(m,iaer,radius_id+1) |
|---|
| 488 | gIRb(m,gausind,iaer)= |
|---|
| 489 | & (1-kint)*gIR(m,iaer,radius_id) + |
|---|
| 490 | & kint*gIR(m,iaer,radius_id+1) |
|---|
| 491 | ENDDO |
|---|
| 492 | qrefIRb(gausind,iaer)= |
|---|
| 493 | & (1-kint)*QREFir(iaer,radius_id) + |
|---|
| 494 | & kint*QREFir(iaer,radius_id+1) |
|---|
| 495 | omegrefIRb(gausind,iaer)= |
|---|
| 496 | & (1-kint)*omegaREFir(iaer,radius_id) + |
|---|
| 497 | & kint*omegaREFir(iaer,radius_id+1) |
|---|
| 498 | ENDIF |
|---|
| 499 | ENDDO |
|---|
| 500 | |
|---|
| 501 | c================================================================== |
|---|
| 502 | IF ( .NOT.varyingnueff(iaer) ) THEN ! CONSTANT NUEFF |
|---|
| 503 | c================================================================== |
|---|
| 504 | c 2. Compute the scattering parameters using linear |
|---|
| 505 | c interpolation over grain sizes and constant nueff |
|---|
| 506 | c --------------------------------------------------- |
|---|
| 507 | |
|---|
| 508 | DO lg = 1,nlayer |
|---|
| 509 | DO ig = 1,ngrid |
|---|
| 510 | c 2.1 Effective radius index and kx calculation |
|---|
| 511 | var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain) |
|---|
| 512 | var_tmp=log(var_tmp)*3. |
|---|
| 513 | var_tmp=var_tmp/logvratgrid(iaer,idomain)+1. |
|---|
| 514 | grid_i=floor(var_tmp) |
|---|
| 515 | IF (grid_i.GE.refftabsize) THEN |
|---|
| 516 | c WRITE(*,*) 'Warning: particle size in grid box #' |
|---|
| 517 | c WRITE(*,*) ig,' is too large to be used by the ' |
|---|
| 518 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 519 | c WRITE(*,*) 'interpolation grid to larger grain sizes.' |
|---|
| 520 | grid_i=refftabsize-1 |
|---|
| 521 | kx = 1. |
|---|
| 522 | ELSEIF (grid_i.LT.1) THEN |
|---|
| 523 | c WRITE(*,*) 'Warning: particle size in grid box #' |
|---|
| 524 | c WRITE(*,*) ig,' is too small to be used by the ' |
|---|
| 525 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 526 | c WRITE(*,*) 'interpolation grid to smaller grain sizes.' |
|---|
| 527 | grid_i=1 |
|---|
| 528 | kx = 0. |
|---|
| 529 | ELSE |
|---|
| 530 | kx = ( reffrad(ig,lg,iaer)- |
|---|
| 531 | & refftab(grid_i,iaer,idomain) ) / |
|---|
| 532 | & ( refftab(grid_i+1,iaer,idomain)- |
|---|
| 533 | & refftab(grid_i,iaer,idomain) ) |
|---|
| 534 | ENDIF |
|---|
| 535 | c 2.3 Integration |
|---|
| 536 | DO j=grid_i,grid_i+1 |
|---|
| 537 | c 2.3.1 Check if the calculation has been done |
|---|
| 538 | IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN |
|---|
| 539 | c 2.3.2 Log-normal dist., r_g and sigma_g are defined |
|---|
| 540 | c in [hansen_1974], "Light scattering in planetary |
|---|
| 541 | c atmospheres", Space Science Reviews 16 527-610. |
|---|
| 542 | c Here, sizedistk1=r_g and sizedistk2=sigma_g^2 |
|---|
| 543 | sizedistk2 = log(1.+nueffrad(1,1,iaer)) |
|---|
| 544 | sizedistk1 = exp(2.5*sizedistk2) |
|---|
| 545 | sizedistk1 = refftab(j,iaer,idomain) / sizedistk1 |
|---|
| 546 | |
|---|
| 547 | normd(j,1,iaer,idomain) = 1e-30 |
|---|
| 548 | DO gausind=1,ngau |
|---|
| 549 | drad=radiusr*radgaus(gausind) |
|---|
| 550 | dista(j,1,iaer,idomain,gausind) = |
|---|
| 551 | & LOG((radiusm-drad)/sizedistk1) |
|---|
| 552 | dista(j,1,iaer,idomain,gausind) = |
|---|
| 553 | & EXP(-dista(j,1,iaer,idomain,gausind) * |
|---|
| 554 | & dista(j,1,iaer,idomain,gausind) * |
|---|
| 555 | & 0.5e0/sizedistk2)/(radiusm-drad) |
|---|
| 556 | dista(j,1,iaer,idomain,gausind) = |
|---|
| 557 | & dista(j,1,iaer,idomain,gausind) / |
|---|
| 558 | & (sqrt(2e0*pi*sizedistk2)) |
|---|
| 559 | |
|---|
| 560 | distb(j,1,iaer,idomain,gausind) = |
|---|
| 561 | & LOG((radiusm+drad)/sizedistk1) |
|---|
| 562 | distb(j,1,iaer,idomain,gausind) = |
|---|
| 563 | & EXP(-distb(j,1,iaer,idomain,gausind) * |
|---|
| 564 | & distb(j,1,iaer,idomain,gausind) * |
|---|
| 565 | & 0.5e0/sizedistk2)/(radiusm+drad) |
|---|
| 566 | distb(j,1,iaer,idomain,gausind) = |
|---|
| 567 | & distb(j,1,iaer,idomain,gausind) / |
|---|
| 568 | & (sqrt(2e0*pi*sizedistk2)) |
|---|
| 569 | |
|---|
| 570 | normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) + |
|---|
| 571 | & weightgaus(gausind) * |
|---|
| 572 | & ( |
|---|
| 573 | & distb(j,1,iaer,idomain,gausind) * pi * |
|---|
| 574 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 575 | & radGAUSb(gausind,iaer,idomain) + |
|---|
| 576 | & dista(j,1,iaer,idomain,gausind) * pi * |
|---|
| 577 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 578 | & radGAUSa(gausind,iaer,idomain) |
|---|
| 579 | & ) |
|---|
| 580 | ENDDO |
|---|
| 581 | |
|---|
| 582 | IF (normd(j,1,iaer,idomain).EQ.1e-30) THEN |
|---|
| 583 | WRITE(*,*)"normd:", normd(j,1,iaer,idomain) |
|---|
| 584 | WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)" |
|---|
| 585 | WRITE(*,*)"Check the size of the interpolation grid." |
|---|
| 586 | CALL abort_physic("aeroptproperties", |
|---|
| 587 | & "Risk of division by 0",1) |
|---|
| 588 | ENDIF |
|---|
| 589 | IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------- |
|---|
| 590 | c 2.3.3.vis Initialization |
|---|
| 591 | qsqrefVISgrid(j,1,:,iaer)=0. |
|---|
| 592 | qextVISgrid(j,1,:,iaer)=0. |
|---|
| 593 | qscatVISgrid(j,1,:,iaer)=0. |
|---|
| 594 | omegVISgrid(j,1,:,iaer)=0. |
|---|
| 595 | gVISgrid(j,1,:,iaer)=0. |
|---|
| 596 | qrefVISgrid(j,1,iaer)=0. |
|---|
| 597 | qscatrefVISgrid(j,1,iaer)=0. |
|---|
| 598 | omegrefVISgrid(j,1,iaer)=0. |
|---|
| 599 | |
|---|
| 600 | DO gausind=1,ngau |
|---|
| 601 | DO m=1,nsun |
|---|
| 602 | c Convolution: |
|---|
| 603 | qextVISgrid(j,1,m,iaer) = |
|---|
| 604 | & qextVISgrid(j,1,m,iaer) + |
|---|
| 605 | & weightgaus(gausind) * |
|---|
| 606 | & ( |
|---|
| 607 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 608 | & qrefVISb(gausind,iaer) * |
|---|
| 609 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 610 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 611 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 612 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 613 | & qrefVISa(gausind,iaer) * |
|---|
| 614 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 615 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 616 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 617 | & ) |
|---|
| 618 | qscatVISgrid(j,1,m,iaer) = |
|---|
| 619 | & qscatVISgrid(j,1,m,iaer) + |
|---|
| 620 | & weightgaus(gausind) * |
|---|
| 621 | & ( |
|---|
| 622 | & omegVISb(m,gausind,iaer) * |
|---|
| 623 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 624 | & qrefVISb(gausind,iaer) * |
|---|
| 625 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 626 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 627 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 628 | & omegVISa(m,gausind,iaer) * |
|---|
| 629 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 630 | & qrefVISa(gausind,iaer) * |
|---|
| 631 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 632 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 633 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 634 | & ) |
|---|
| 635 | gVISgrid(j,1,m,iaer) = |
|---|
| 636 | & gVISgrid(j,1,m,iaer) + |
|---|
| 637 | & weightgaus(gausind) * |
|---|
| 638 | & ( |
|---|
| 639 | & omegVISb(m,gausind,iaer) * |
|---|
| 640 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 641 | & qrefVISb(gausind,iaer) * |
|---|
| 642 | & gVISb(m,gausind,iaer) * |
|---|
| 643 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 644 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 645 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 646 | & omegVISa(m,gausind,iaer) * |
|---|
| 647 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 648 | & qrefVISa(gausind,iaer) * |
|---|
| 649 | & gVISa(m,gausind,iaer) * |
|---|
| 650 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 651 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 652 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 653 | & ) |
|---|
| 654 | ENDDO |
|---|
| 655 | qrefVISgrid(j,1,iaer) = |
|---|
| 656 | & qrefVISgrid(j,1,iaer) + |
|---|
| 657 | & weightgaus(gausind) * |
|---|
| 658 | & ( |
|---|
| 659 | & qrefVISb(gausind,iaer) * |
|---|
| 660 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 661 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 662 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 663 | & qrefVISa(gausind,iaer) * |
|---|
| 664 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 665 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 666 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 667 | & ) |
|---|
| 668 | qscatrefVISgrid(j,1,iaer) = |
|---|
| 669 | & qscatrefVISgrid(j,1,iaer) + |
|---|
| 670 | & weightgaus(gausind) * |
|---|
| 671 | & ( |
|---|
| 672 | & omegrefVISb(gausind,iaer) * |
|---|
| 673 | & qrefVISb(gausind,iaer) * |
|---|
| 674 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 675 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 676 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 677 | & omegrefVISa(gausind,iaer) * |
|---|
| 678 | & qrefVISa(gausind,iaer) * |
|---|
| 679 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 680 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 681 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 682 | & ) |
|---|
| 683 | ENDDO |
|---|
| 684 | |
|---|
| 685 | qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) / |
|---|
| 686 | & normd(j,1,iaer,idomain) |
|---|
| 687 | qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / |
|---|
| 688 | & normd(j,1,iaer,idomain) |
|---|
| 689 | omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) / |
|---|
| 690 | & qrefVISgrid(j,1,iaer) |
|---|
| 691 | DO m=1,nsun |
|---|
| 692 | qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / |
|---|
| 693 | & normd(j,1,iaer,idomain) |
|---|
| 694 | qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / |
|---|
| 695 | & normd(j,1,iaer,idomain) |
|---|
| 696 | gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) / |
|---|
| 697 | & qscatVISgrid(j,1,m,iaer) / |
|---|
| 698 | & normd(j,1,iaer,idomain) |
|---|
| 699 | |
|---|
| 700 | qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) / |
|---|
| 701 | & qrefVISgrid(j,1,iaer) |
|---|
| 702 | omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) / |
|---|
| 703 | & qextVISgrid(j,1,m,iaer) |
|---|
| 704 | ENDDO |
|---|
| 705 | ELSE ! INFRARED DOMAIN ---------- |
|---|
| 706 | c 2.3.3.ir Initialization |
|---|
| 707 | qsqrefIRgrid(j,1,:,iaer)=0. |
|---|
| 708 | qextIRgrid(j,1,:,iaer)=0. |
|---|
| 709 | qscatIRgrid(j,1,:,iaer)=0. |
|---|
| 710 | omegIRgrid(j,1,:,iaer)=0. |
|---|
| 711 | gIRgrid(j,1,:,iaer)=0. |
|---|
| 712 | qrefIRgrid(j,1,iaer)=0. |
|---|
| 713 | qscatrefIRgrid(j,1,iaer)=0. |
|---|
| 714 | omegrefIRgrid(j,1,iaer)=0. |
|---|
| 715 | |
|---|
| 716 | DO gausind=1,ngau |
|---|
| 717 | DO m=1,nir |
|---|
| 718 | c Convolution: |
|---|
| 719 | qextIRgrid(j,1,m,iaer) = |
|---|
| 720 | & qextIRgrid(j,1,m,iaer) + |
|---|
| 721 | & weightgaus(gausind) * |
|---|
| 722 | & ( |
|---|
| 723 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 724 | & qrefVISb(gausind,iaer) * |
|---|
| 725 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 726 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 727 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 728 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 729 | & qrefVISa(gausind,iaer) * |
|---|
| 730 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 731 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 732 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 733 | & ) |
|---|
| 734 | qscatIRgrid(j,1,m,iaer) = |
|---|
| 735 | & qscatIRgrid(j,1,m,iaer) + |
|---|
| 736 | & weightgaus(gausind) * |
|---|
| 737 | & ( |
|---|
| 738 | & omegIRb(m,gausind,iaer) * |
|---|
| 739 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 740 | & qrefVISb(gausind,iaer) * |
|---|
| 741 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 742 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 743 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 744 | & omegIRa(m,gausind,iaer) * |
|---|
| 745 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 746 | & qrefVISa(gausind,iaer) * |
|---|
| 747 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 748 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 749 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 750 | & ) |
|---|
| 751 | gIRgrid(j,1,m,iaer) = |
|---|
| 752 | & gIRgrid(j,1,m,iaer) + |
|---|
| 753 | & weightgaus(gausind) * |
|---|
| 754 | & ( |
|---|
| 755 | & omegIRb(m,gausind,iaer) * |
|---|
| 756 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 757 | & qrefVISb(gausind,iaer) * |
|---|
| 758 | & gIRb(m,gausind,iaer) * |
|---|
| 759 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 760 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 761 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 762 | & omegIRa(m,gausind,iaer) * |
|---|
| 763 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 764 | & qrefVISa(gausind,iaer) * |
|---|
| 765 | & gIRa(m,gausind,iaer) * |
|---|
| 766 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 767 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 768 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 769 | & ) |
|---|
| 770 | ENDDO |
|---|
| 771 | qrefIRgrid(j,1,iaer) = |
|---|
| 772 | & qrefIRgrid(j,1,iaer) + |
|---|
| 773 | & weightgaus(gausind) * |
|---|
| 774 | & ( |
|---|
| 775 | & qrefIRb(gausind,iaer) * |
|---|
| 776 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 777 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 778 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 779 | & qrefIRa(gausind,iaer) * |
|---|
| 780 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 781 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 782 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 783 | & ) |
|---|
| 784 | qscatrefIRgrid(j,1,iaer) = |
|---|
| 785 | & qscatrefIRgrid(j,1,iaer) + |
|---|
| 786 | & weightgaus(gausind) * |
|---|
| 787 | & ( |
|---|
| 788 | & omegrefIRb(gausind,iaer) * |
|---|
| 789 | & qrefIRb(gausind,iaer) * |
|---|
| 790 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 791 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 792 | & distb(j,1,iaer,idomain,gausind) + |
|---|
| 793 | & omegrefIRa(gausind,iaer) * |
|---|
| 794 | & qrefIRa(gausind,iaer) * |
|---|
| 795 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 796 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 797 | & dista(j,1,iaer,idomain,gausind) |
|---|
| 798 | & ) |
|---|
| 799 | ENDDO |
|---|
| 800 | |
|---|
| 801 | qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) / |
|---|
| 802 | & normd(j,1,iaer,idomain) |
|---|
| 803 | qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / |
|---|
| 804 | & normd(j,1,iaer,idomain) |
|---|
| 805 | omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) / |
|---|
| 806 | & qrefIRgrid(j,1,iaer) |
|---|
| 807 | DO m=1,nir |
|---|
| 808 | qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / |
|---|
| 809 | & normd(j,1,iaer,idomain) |
|---|
| 810 | qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / |
|---|
| 811 | & normd(j,1,iaer,idomain) |
|---|
| 812 | gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) / |
|---|
| 813 | & qscatIRgrid(j,1,m,iaer) / |
|---|
| 814 | & normd(j,1,iaer,idomain) |
|---|
| 815 | |
|---|
| 816 | qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) / |
|---|
| 817 | & qrefVISgrid(j,1,iaer) |
|---|
| 818 | omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) / |
|---|
| 819 | & qextIRgrid(j,1,m,iaer) |
|---|
| 820 | ENDDO |
|---|
| 821 | ENDIF ! -------------------------- |
|---|
| 822 | checkgrid(j,1,iaer,idomain) = .true. |
|---|
| 823 | ENDIF !checkgrid |
|---|
| 824 | ENDDO !grid_i |
|---|
| 825 | c 2.4 Linear interpolation |
|---|
| 826 | k1 = (1-kx) |
|---|
| 827 | k2 = kx |
|---|
| 828 | IF (idomain.EQ.1) THEN ! VISIBLE ------------------------ |
|---|
| 829 | DO m=1,nsun |
|---|
| 830 | QVISsQREF3d(ig,lg,m,iaer) = |
|---|
| 831 | & k1*qsqrefVISgrid(grid_i,1,m,iaer) + |
|---|
| 832 | & k2*qsqrefVISgrid(grid_i+1,1,m,iaer) |
|---|
| 833 | omegaVIS3d(ig,lg,m,iaer) = |
|---|
| 834 | & k1*omegVISgrid(grid_i,1,m,iaer) + |
|---|
| 835 | & k2*omegVISgrid(grid_i+1,1,m,iaer) |
|---|
| 836 | gVIS3d(ig,lg,m,iaer) = |
|---|
| 837 | & k1*gVISgrid(grid_i,1,m,iaer) + |
|---|
| 838 | & k2*gVISgrid(grid_i+1,1,m,iaer) |
|---|
| 839 | ENDDO !nsun |
|---|
| 840 | QREFvis3d(ig,lg,iaer) = |
|---|
| 841 | & k1*qrefVISgrid(grid_i,1,iaer) + |
|---|
| 842 | & k2*qrefVISgrid(grid_i+1,1,iaer) |
|---|
| 843 | omegaREFvis3d(ig,lg,iaer) = |
|---|
| 844 | & k1*omegrefVISgrid(grid_i,1,iaer) + |
|---|
| 845 | & k2*omegrefVISgrid(grid_i+1,1,iaer) |
|---|
| 846 | ELSE ! INFRARED ----------------------- |
|---|
| 847 | DO m=1,nir |
|---|
| 848 | QIRsQREF3d(ig,lg,m,iaer) = |
|---|
| 849 | & k1*qsqrefIRgrid(grid_i,1,m,iaer) + |
|---|
| 850 | & k2*qsqrefIRgrid(grid_i+1,1,m,iaer) |
|---|
| 851 | omegaIR3d(ig,lg,m,iaer) = |
|---|
| 852 | & k1*omegIRgrid(grid_i,1,m,iaer) + |
|---|
| 853 | & k2*omegIRgrid(grid_i+1,1,m,iaer) |
|---|
| 854 | gIR3d(ig,lg,m,iaer) = |
|---|
| 855 | & k1*gIRgrid(grid_i,1,m,iaer) + |
|---|
| 856 | & k2*gIRgrid(grid_i+1,1,m,iaer) |
|---|
| 857 | ENDDO !nir |
|---|
| 858 | QREFir3d(ig,lg,iaer) = |
|---|
| 859 | & k1*qrefIRgrid(grid_i,1,iaer) + |
|---|
| 860 | & k2*qrefIRgrid(grid_i+1,1,iaer) |
|---|
| 861 | omegaREFir3d(ig,lg,iaer) = |
|---|
| 862 | & k1*omegrefIRgrid(grid_i,1,iaer) + |
|---|
| 863 | & k2*omegrefIRgrid(grid_i+1,1,iaer) |
|---|
| 864 | ENDIF ! -------------------------------- |
|---|
| 865 | ENDDO !nlayer |
|---|
| 866 | ENDDO !ngrid |
|---|
| 867 | c================================================================== |
|---|
| 868 | ELSE ! VARYING NUEFF |
|---|
| 869 | c================================================================== |
|---|
| 870 | c 3. Compute the scattering parameters in each grid box |
|---|
| 871 | c using bilinear interpolation over the grain sizes |
|---|
| 872 | c and the effective variances; |
|---|
| 873 | c ----------------------------------------------------- |
|---|
| 874 | |
|---|
| 875 | DO lg = 1,nlayer |
|---|
| 876 | DO ig = 1,ngrid |
|---|
| 877 | c 3.1 Effective variance index and ky calculation |
|---|
| 878 | var_tmp=log(nueffrad(ig,lg,iaer)) |
|---|
| 879 | grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)* |
|---|
| 880 | & (var_tmp-nuefftabmin)+1. ) |
|---|
| 881 | IF (grid_j.GE.nuefftabsize) THEN |
|---|
| 882 | c WRITE(*,*) 'Warning: effective variance ' |
|---|
| 883 | c WRITE(*,*) 'is too large to be used by the ' |
|---|
| 884 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 885 | c WRITE(*,*) 'interpolation grid to larger values.' |
|---|
| 886 | grid_j=nuefftabsize-1 |
|---|
| 887 | ky = 1. |
|---|
| 888 | ELSEIF (grid_j.LT.1) THEN |
|---|
| 889 | c WRITE(*,*) 'Warning: effective variance ' |
|---|
| 890 | c WRITE(*,*) 'is too small to be used by the ' |
|---|
| 891 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 892 | c WRITE(*,*) 'interpolation grid to smaller values.' |
|---|
| 893 | grid_j=1 |
|---|
| 894 | ky = 0. |
|---|
| 895 | ELSE |
|---|
| 896 | ky = ( nueffrad(ig,lg,iaer)- |
|---|
| 897 | & nuefftab(grid_j,iaer,idomain) ) / |
|---|
| 898 | & ( nuefftab(grid_j+1,iaer,idomain)- |
|---|
| 899 | & nuefftab(grid_j,iaer,idomain) ) |
|---|
| 900 | ENDIF |
|---|
| 901 | c 3.2 Effective radius index and kx calculation |
|---|
| 902 | var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain) |
|---|
| 903 | var_tmp=log(var_tmp)*3. |
|---|
| 904 | var_tmp=var_tmp/logvratgrid(iaer,idomain)+1. |
|---|
| 905 | grid_i=floor(var_tmp) |
|---|
| 906 | IF (grid_i.GE.refftabsize) THEN |
|---|
| 907 | c WRITE(*,*) 'Warning: particle size in grid box #' |
|---|
| 908 | c WRITE(*,*) ig,' is too large to be used by the ' |
|---|
| 909 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 910 | c WRITE(*,*) 'interpolation grid to larger grain sizes.' |
|---|
| 911 | grid_i=refftabsize-1 |
|---|
| 912 | kx = 1. |
|---|
| 913 | ELSEIF (grid_i.LT.1) THEN |
|---|
| 914 | c WRITE(*,*) 'Warning: particle size in grid box #' |
|---|
| 915 | c WRITE(*,*) ig,' is too small to be used by the ' |
|---|
| 916 | c WRITE(*,*) 'radiative transfer; please extend the ' |
|---|
| 917 | c WRITE(*,*) 'interpolation grid to smaller grain sizes.' |
|---|
| 918 | grid_i=1 |
|---|
| 919 | kx = 0. |
|---|
| 920 | ELSE |
|---|
| 921 | kx = ( reffrad(ig,lg,iaer)- |
|---|
| 922 | & refftab(grid_i,iaer,idomain) ) / |
|---|
| 923 | & ( refftab(grid_i+1,iaer,idomain)- |
|---|
| 924 | & refftab(grid_i,iaer,idomain) ) |
|---|
| 925 | ENDIF |
|---|
| 926 | c 3.3 Integration |
|---|
| 927 | DO j=grid_i,grid_i+1 |
|---|
| 928 | DO k=grid_j,grid_j+1 |
|---|
| 929 | c 3.3.1 Check if the calculation has been done |
|---|
| 930 | IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN |
|---|
| 931 | |
|---|
| 932 | c 3.3.2 Log-normal dist., r_g and sigma_g are defined |
|---|
| 933 | c in [hansen_1974], "Light scattering in planetary |
|---|
| 934 | c atmospheres", Space Science Reviews 16 527-610. |
|---|
| 935 | c Here, sizedistk1=r_g and sizedistk2=sigma_g^2 |
|---|
| 936 | sizedistk2 = log(1.+nuefftab(k,iaer,idomain)) |
|---|
| 937 | sizedistk1 = exp(2.5*sizedistk2) |
|---|
| 938 | sizedistk1 = refftab(j,iaer,idomain) / sizedistk1 |
|---|
| 939 | |
|---|
| 940 | normd(j,k,iaer,idomain) = 1e-30 |
|---|
| 941 | DO gausind=1,ngau |
|---|
| 942 | drad=radiusr*radgaus(gausind) |
|---|
| 943 | |
|---|
| 944 | dista(j,k,iaer,idomain,gausind) = |
|---|
| 945 | & LOG((radiusm-drad)/sizedistk1) |
|---|
| 946 | dista(j,k,iaer,idomain,gausind) = |
|---|
| 947 | & EXP(-dista(j,k,iaer,idomain,gausind) * |
|---|
| 948 | & dista(j,k,iaer,idomain,gausind) * |
|---|
| 949 | & 0.5e0/sizedistk2)/(radiusm-drad) |
|---|
| 950 | dista(j,k,iaer,idomain,gausind) = |
|---|
| 951 | & dista(j,k,iaer,idomain,gausind) / |
|---|
| 952 | & (sqrt(2e0*pi*sizedistk2)) |
|---|
| 953 | |
|---|
| 954 | distb(j,k,iaer,idomain,gausind) = |
|---|
| 955 | & LOG((radiusm+drad)/sizedistk1) |
|---|
| 956 | distb(j,k,iaer,idomain,gausind) = |
|---|
| 957 | & EXP(-distb(j,k,iaer,idomain,gausind) * |
|---|
| 958 | & distb(j,k,iaer,idomain,gausind) * |
|---|
| 959 | & 0.5e0/sizedistk2)/(radiusm+drad) |
|---|
| 960 | distb(j,k,iaer,idomain,gausind) = |
|---|
| 961 | & distb(j,k,iaer,idomain,gausind) / |
|---|
| 962 | & (sqrt(2e0*pi*sizedistk2)) |
|---|
| 963 | |
|---|
| 964 | normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) + |
|---|
| 965 | & weightgaus(gausind) * |
|---|
| 966 | & ( |
|---|
| 967 | & distb(j,k,iaer,idomain,gausind) * pi * |
|---|
| 968 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 969 | & radGAUSb(gausind,iaer,idomain) + |
|---|
| 970 | & dista(j,k,iaer,idomain,gausind) * pi * |
|---|
| 971 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 972 | & radGAUSa(gausind,iaer,idomain) |
|---|
| 973 | & ) |
|---|
| 974 | ENDDO |
|---|
| 975 | IF (normd(j,k,iaer,idomain).EQ.1e-30) THEN |
|---|
| 976 | WRITE(*,*)"normd:", normd(j,k,iaer,idomain) |
|---|
| 977 | WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)" |
|---|
| 978 | WRITE(*,*)"Check the size of the interpolation grid." |
|---|
| 979 | CALL abort_physic("aeroptproperties", |
|---|
| 980 | & "Risk of division by 0",1) |
|---|
| 981 | ENDIF |
|---|
| 982 | IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN ----------- |
|---|
| 983 | c 2.3.3.vis Initialization |
|---|
| 984 | qsqrefVISgrid(j,k,:,iaer)=0. |
|---|
| 985 | qextVISgrid(j,k,:,iaer)=0. |
|---|
| 986 | qscatVISgrid(j,k,:,iaer)=0. |
|---|
| 987 | omegVISgrid(j,k,:,iaer)=0. |
|---|
| 988 | gVISgrid(j,k,:,iaer)=0. |
|---|
| 989 | qrefVISgrid(j,k,iaer)=0. |
|---|
| 990 | qscatrefVISgrid(j,k,iaer)=0. |
|---|
| 991 | omegrefVISgrid(j,k,iaer)=0. |
|---|
| 992 | |
|---|
| 993 | DO gausind=1,ngau |
|---|
| 994 | DO m=1,nsun |
|---|
| 995 | c Convolution: |
|---|
| 996 | qextVISgrid(j,k,m,iaer) = |
|---|
| 997 | & qextVISgrid(j,k,m,iaer) + |
|---|
| 998 | & weightgaus(gausind) * |
|---|
| 999 | & ( |
|---|
| 1000 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 1001 | & qrefVISb(gausind,iaer) * |
|---|
| 1002 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1003 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1004 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1005 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 1006 | & qrefVISa(gausind,iaer) * |
|---|
| 1007 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1008 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1009 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1010 | & ) |
|---|
| 1011 | qscatVISgrid(j,k,m,iaer) = |
|---|
| 1012 | & qscatVISgrid(j,k,m,iaer) + |
|---|
| 1013 | & weightgaus(gausind) * |
|---|
| 1014 | & ( |
|---|
| 1015 | & omegVISb(m,gausind,iaer) * |
|---|
| 1016 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 1017 | & qrefVISb(gausind,iaer) * |
|---|
| 1018 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1019 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1020 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1021 | & omegVISa(m,gausind,iaer) * |
|---|
| 1022 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 1023 | & qrefVISa(gausind,iaer) * |
|---|
| 1024 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1025 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1026 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1027 | & ) |
|---|
| 1028 | gVISgrid(j,k,m,iaer) = |
|---|
| 1029 | & gVISgrid(j,k,m,iaer) + |
|---|
| 1030 | & weightgaus(gausind) * |
|---|
| 1031 | & ( |
|---|
| 1032 | & omegVISb(m,gausind,iaer) * |
|---|
| 1033 | & qsqrefVISb(m,gausind,iaer) * |
|---|
| 1034 | & qrefVISb(gausind,iaer) * |
|---|
| 1035 | & gVISb(m,gausind,iaer) * |
|---|
| 1036 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1037 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1038 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1039 | & omegVISa(m,gausind,iaer) * |
|---|
| 1040 | & qsqrefVISa(m,gausind,iaer) * |
|---|
| 1041 | & qrefVISa(gausind,iaer) * |
|---|
| 1042 | & gVISa(m,gausind,iaer) * |
|---|
| 1043 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1044 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1045 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1046 | & ) |
|---|
| 1047 | ENDDO |
|---|
| 1048 | qrefVISgrid(j,k,iaer) = |
|---|
| 1049 | & qrefVISgrid(j,k,iaer) + |
|---|
| 1050 | & weightgaus(gausind) * |
|---|
| 1051 | & ( |
|---|
| 1052 | & qrefVISb(gausind,iaer) * |
|---|
| 1053 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1054 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1055 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1056 | & qrefVISa(gausind,iaer) * |
|---|
| 1057 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1058 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1059 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1060 | & ) |
|---|
| 1061 | qscatrefVISgrid(j,k,iaer) = |
|---|
| 1062 | & qscatrefVISgrid(j,k,iaer) + |
|---|
| 1063 | & weightgaus(gausind) * |
|---|
| 1064 | & ( |
|---|
| 1065 | & omegrefVISb(gausind,iaer) * |
|---|
| 1066 | & qrefVISb(gausind,iaer) * |
|---|
| 1067 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1068 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1069 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1070 | & omegrefVISa(gausind,iaer) * |
|---|
| 1071 | & qrefVISa(gausind,iaer) * |
|---|
| 1072 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1073 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1074 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1075 | & ) |
|---|
| 1076 | ENDDO |
|---|
| 1077 | qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) / |
|---|
| 1078 | & normd(j,k,iaer,idomain) |
|---|
| 1079 | qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) / |
|---|
| 1080 | & normd(j,k,iaer,idomain) |
|---|
| 1081 | omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) / |
|---|
| 1082 | & qrefVISgrid(j,k,iaer) |
|---|
| 1083 | DO m=1,nsun |
|---|
| 1084 | qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) / |
|---|
| 1085 | & normd(j,k,iaer,idomain) |
|---|
| 1086 | qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) / |
|---|
| 1087 | & normd(j,k,iaer,idomain) |
|---|
| 1088 | gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) / |
|---|
| 1089 | & qscatVISgrid(j,k,m,iaer) / |
|---|
| 1090 | & normd(j,k,iaer,idomain) |
|---|
| 1091 | |
|---|
| 1092 | qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) / |
|---|
| 1093 | & qrefVISgrid(j,k,iaer) |
|---|
| 1094 | omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) / |
|---|
| 1095 | & qextVISgrid(j,k,m,iaer) |
|---|
| 1096 | ENDDO |
|---|
| 1097 | ELSE ! INFRARED DOMAIN ---------- |
|---|
| 1098 | c 2.3.3.ir Initialization |
|---|
| 1099 | qsqrefIRgrid(j,k,:,iaer)=0. |
|---|
| 1100 | qextIRgrid(j,k,:,iaer)=0. |
|---|
| 1101 | qscatIRgrid(j,k,:,iaer)=0. |
|---|
| 1102 | omegIRgrid(j,k,:,iaer)=0. |
|---|
| 1103 | gIRgrid(j,k,:,iaer)=0. |
|---|
| 1104 | qrefIRgrid(j,k,iaer)=0. |
|---|
| 1105 | qscatrefIRgrid(j,k,iaer)=0. |
|---|
| 1106 | omegrefIRgrid(j,k,iaer)=0. |
|---|
| 1107 | |
|---|
| 1108 | DO gausind=1,ngau |
|---|
| 1109 | DO m=1,nir |
|---|
| 1110 | c Convolution: |
|---|
| 1111 | qextIRgrid(j,k,m,iaer) = |
|---|
| 1112 | & qextIRgrid(j,k,m,iaer) + |
|---|
| 1113 | & weightgaus(gausind) * |
|---|
| 1114 | & ( |
|---|
| 1115 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 1116 | & qrefVISb(gausind,iaer) * |
|---|
| 1117 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1118 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1119 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1120 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 1121 | & qrefVISa(gausind,iaer) * |
|---|
| 1122 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1123 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1124 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1125 | & ) |
|---|
| 1126 | qscatIRgrid(j,k,m,iaer) = |
|---|
| 1127 | & qscatIRgrid(j,k,m,iaer) + |
|---|
| 1128 | & weightgaus(gausind) * |
|---|
| 1129 | & ( |
|---|
| 1130 | & omegIRb(m,gausind,iaer) * |
|---|
| 1131 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 1132 | & qrefVISb(gausind,iaer) * |
|---|
| 1133 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1134 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1135 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1136 | & omegIRa(m,gausind,iaer) * |
|---|
| 1137 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 1138 | & qrefVISa(gausind,iaer) * |
|---|
| 1139 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1140 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1141 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1142 | & ) |
|---|
| 1143 | gIRgrid(j,k,m,iaer) = |
|---|
| 1144 | & gIRgrid(j,k,m,iaer) + |
|---|
| 1145 | & weightgaus(gausind) * |
|---|
| 1146 | & ( |
|---|
| 1147 | & omegIRb(m,gausind,iaer) * |
|---|
| 1148 | & qsqrefIRb(m,gausind,iaer) * |
|---|
| 1149 | & qrefVISb(gausind,iaer) * |
|---|
| 1150 | & gIRb(m,gausind,iaer) * |
|---|
| 1151 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1152 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1153 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1154 | & omegIRa(m,gausind,iaer) * |
|---|
| 1155 | & qsqrefIRa(m,gausind,iaer) * |
|---|
| 1156 | & qrefVISa(gausind,iaer) * |
|---|
| 1157 | & gIRa(m,gausind,iaer) * |
|---|
| 1158 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1159 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1160 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1161 | & ) |
|---|
| 1162 | ENDDO |
|---|
| 1163 | qrefIRgrid(j,k,iaer) = |
|---|
| 1164 | & qrefIRgrid(j,k,iaer) + |
|---|
| 1165 | & weightgaus(gausind) * |
|---|
| 1166 | & ( |
|---|
| 1167 | & qrefIRb(gausind,iaer) * |
|---|
| 1168 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1169 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1170 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1171 | & qrefIRa(gausind,iaer) * |
|---|
| 1172 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1173 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1174 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1175 | & ) |
|---|
| 1176 | qscatrefIRgrid(j,k,iaer) = |
|---|
| 1177 | & qscatrefIRgrid(j,k,iaer) + |
|---|
| 1178 | & weightgaus(gausind) * |
|---|
| 1179 | & ( |
|---|
| 1180 | & omegrefIRb(gausind,iaer) * |
|---|
| 1181 | & qrefIRb(gausind,iaer) * |
|---|
| 1182 | & pi*radGAUSb(gausind,iaer,idomain) * |
|---|
| 1183 | & radGAUSb(gausind,iaer,idomain) * |
|---|
| 1184 | & distb(j,k,iaer,idomain,gausind) + |
|---|
| 1185 | & omegrefIRa(gausind,iaer) * |
|---|
| 1186 | & qrefIRa(gausind,iaer) * |
|---|
| 1187 | & pi*radGAUSa(gausind,iaer,idomain) * |
|---|
| 1188 | & radGAUSa(gausind,iaer,idomain) * |
|---|
| 1189 | & dista(j,k,iaer,idomain,gausind) |
|---|
| 1190 | & ) |
|---|
| 1191 | ENDDO |
|---|
| 1192 | qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) / |
|---|
| 1193 | & normd(j,k,iaer,idomain) |
|---|
| 1194 | qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) / |
|---|
| 1195 | & normd(j,k,iaer,idomain) |
|---|
| 1196 | omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) / |
|---|
| 1197 | & qrefIRgrid(j,k,iaer) |
|---|
| 1198 | DO m=1,nir |
|---|
| 1199 | qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) / |
|---|
| 1200 | & normd(j,k,iaer,idomain) |
|---|
| 1201 | qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) / |
|---|
| 1202 | & normd(j,k,iaer,idomain) |
|---|
| 1203 | gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) / |
|---|
| 1204 | & qscatIRgrid(j,k,m,iaer) / |
|---|
| 1205 | & normd(j,k,iaer,idomain) |
|---|
| 1206 | |
|---|
| 1207 | qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) / |
|---|
| 1208 | & qrefVISgrid(j,k,iaer) |
|---|
| 1209 | omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) / |
|---|
| 1210 | & qextIRgrid(j,k,m,iaer) |
|---|
| 1211 | ENDDO |
|---|
| 1212 | ENDIF ! -------------------------- |
|---|
| 1213 | checkgrid(j,k,iaer,idomain) = .true. |
|---|
| 1214 | ENDIF !checkgrid |
|---|
| 1215 | ENDDO !grid_j |
|---|
| 1216 | ENDDO !grid_i |
|---|
| 1217 | c 3.4 Bilinear interpolation |
|---|
| 1218 | k1 = (1-kx)*(1-ky) |
|---|
| 1219 | k2 = kx*(1-ky) |
|---|
| 1220 | k3 = kx*ky |
|---|
| 1221 | k4 = (1-kx)*ky |
|---|
| 1222 | IF (idomain.EQ.1) THEN ! VISIBLE ------------------------ |
|---|
| 1223 | DO m=1,nsun |
|---|
| 1224 | QVISsQREF3d(ig,lg,m,iaer) = |
|---|
| 1225 | & k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1226 | & k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1227 | & k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1228 | & k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1229 | omegaVIS3d(ig,lg,m,iaer) = |
|---|
| 1230 | & k1*omegVISgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1231 | & k2*omegVISgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1232 | & k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1233 | & k4*omegVISgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1234 | gVIS3d(ig,lg,m,iaer) = |
|---|
| 1235 | & k1*gVISgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1236 | & k2*gVISgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1237 | & k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1238 | & k4*gVISgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1239 | ENDDO !nsun |
|---|
| 1240 | QREFvis3d(ig,lg,iaer) = |
|---|
| 1241 | & k1*qrefVISgrid(grid_i,grid_j,iaer) + |
|---|
| 1242 | & k2*qrefVISgrid(grid_i+1,grid_j,iaer) + |
|---|
| 1243 | & k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) + |
|---|
| 1244 | & k4*qrefVISgrid(grid_i,grid_j+1,iaer) |
|---|
| 1245 | omegaREFvis3d(ig,lg,iaer) = |
|---|
| 1246 | & k1*omegrefVISgrid(grid_i,grid_j,iaer) + |
|---|
| 1247 | & k2*omegrefVISgrid(grid_i+1,grid_j,iaer) + |
|---|
| 1248 | & k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) + |
|---|
| 1249 | & k4*omegrefVISgrid(grid_i,grid_j+1,iaer) |
|---|
| 1250 | ELSE ! INFRARED ----------------------- |
|---|
| 1251 | DO m=1,nir |
|---|
| 1252 | QIRsQREF3d(ig,lg,m,iaer) = |
|---|
| 1253 | & k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1254 | & k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1255 | & k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1256 | & k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1257 | omegaIR3d(ig,lg,m,iaer) = |
|---|
| 1258 | & k1*omegIRgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1259 | & k2*omegIRgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1260 | & k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1261 | & k4*omegIRgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1262 | gIR3d(ig,lg,m,iaer) = |
|---|
| 1263 | & k1*gIRgrid(grid_i,grid_j,m,iaer) + |
|---|
| 1264 | & k2*gIRgrid(grid_i+1,grid_j,m,iaer) + |
|---|
| 1265 | & k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) + |
|---|
| 1266 | & k4*gIRgrid(grid_i,grid_j+1,m,iaer) |
|---|
| 1267 | ENDDO !nir |
|---|
| 1268 | QREFir3d(ig,lg,iaer) = |
|---|
| 1269 | & k1*qrefIRgrid(grid_i,grid_j,iaer) + |
|---|
| 1270 | & k2*qrefIRgrid(grid_i+1,grid_j,iaer) + |
|---|
| 1271 | & k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) + |
|---|
| 1272 | & k4*qrefIRgrid(grid_i,grid_j+1,iaer) |
|---|
| 1273 | omegaREFir3d(ig,lg,iaer) = |
|---|
| 1274 | & k1*omegrefIRgrid(grid_i,grid_j,iaer) + |
|---|
| 1275 | & k2*omegrefIRgrid(grid_i+1,grid_j,iaer) + |
|---|
| 1276 | & k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) + |
|---|
| 1277 | & k4*omegrefIRgrid(grid_i,grid_j+1,iaer) |
|---|
| 1278 | ENDIF ! -------------------------------- |
|---|
| 1279 | ENDDO !nlayer |
|---|
| 1280 | ENDDO !ngrid |
|---|
| 1281 | |
|---|
| 1282 | ENDIF ! varyingnueff |
|---|
| 1283 | c================================================================== |
|---|
| 1284 | ENDDO ! idomain |
|---|
| 1285 | |
|---|
| 1286 | ENDIF ! nsize = 1 |
|---|
| 1287 | |
|---|
| 1288 | ENDDO ! iaer (loop on aerosol kind) |
|---|
| 1289 | |
|---|
| 1290 | c=====Radiative properties - TESTS================================= |
|---|
| 1291 | IF (out_qwg) THEN |
|---|
| 1292 | c ------------------------------------------------------------- |
|---|
| 1293 | IF (ngrid.NE.1) THEN |
|---|
| 1294 | out_ndim = 3 |
|---|
| 1295 | ELSE |
|---|
| 1296 | out_ndim = 1 |
|---|
| 1297 | ENDIF |
|---|
| 1298 | c ------------------------------------------------------------- |
|---|
| 1299 | DO out_nchannel = 1, 2 |
|---|
| 1300 | c ------------------------------------------------------------- |
|---|
| 1301 | DO lg = 1, nlayer |
|---|
| 1302 | DO ig = 1, ngrid |
|---|
| 1303 | out_qext(ig,lg) = |
|---|
| 1304 | & QVISsQREF3d(ig,lg,out_nchannel,out_iaer)* |
|---|
| 1305 | & QREFvis3d(ig,lg,out_iaer) |
|---|
| 1306 | out_omeg(ig,lg) = |
|---|
| 1307 | & omegaVIS3d(ig,lg,out_nchannel,out_iaer) |
|---|
| 1308 | out_g(ig,lg) = gVIS3d(ig,lg,out_nchannel,out_iaer) |
|---|
| 1309 | ENDDO ! ig |
|---|
| 1310 | ENDDO ! lg |
|---|
| 1311 | #ifndef MESOSCALE |
|---|
| 1312 | write(out_str(1:1),'(i1.1)') out_nchannel |
|---|
| 1313 | call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","", |
|---|
| 1314 | & out_ndim,out_qext) |
|---|
| 1315 | call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","", |
|---|
| 1316 | & out_ndim,out_omeg) |
|---|
| 1317 | call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","", |
|---|
| 1318 | & out_ndim,out_g) |
|---|
| 1319 | #endif |
|---|
| 1320 | c ------------------------------------------------------------- |
|---|
| 1321 | ENDDO ! out_nchannel |
|---|
| 1322 | DO out_nchannel = 2, 4 |
|---|
| 1323 | c ------------------------------------------------------------- |
|---|
| 1324 | DO lg = 1, nlayer |
|---|
| 1325 | DO ig = 1, ngrid |
|---|
| 1326 | out_qext(ig,lg) = |
|---|
| 1327 | & QIRsQREF3d(ig,lg,out_nchannel,out_iaer)* |
|---|
| 1328 | & QREFir3d(ig,lg,out_iaer) |
|---|
| 1329 | out_omeg(ig,lg) = |
|---|
| 1330 | & omegaIR3d(ig,lg,out_nchannel,out_iaer) |
|---|
| 1331 | out_g(ig,lg) = gIR3d(ig,lg,out_nchannel,out_iaer) |
|---|
| 1332 | ENDDO ! ig |
|---|
| 1333 | ENDDO ! lg |
|---|
| 1334 | #ifndef MESOSCALE |
|---|
| 1335 | write(out_str(1:1),'(i1.1)') out_nchannel |
|---|
| 1336 | call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","", |
|---|
| 1337 | & out_ndim,out_qext) |
|---|
| 1338 | call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","", |
|---|
| 1339 | & out_ndim,out_omeg) |
|---|
| 1340 | call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","", |
|---|
| 1341 | & out_ndim,out_g) |
|---|
| 1342 | #endif |
|---|
| 1343 | c ------------------------------------------------------------- |
|---|
| 1344 | ENDDO ! out_nchannel |
|---|
| 1345 | #ifndef MESOSCALE |
|---|
| 1346 | call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","", |
|---|
| 1347 | & out_ndim,omegaREFvis3d(1,1,out_iaer)) |
|---|
| 1348 | call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","", |
|---|
| 1349 | & out_ndim,omegaREFir3d(1,1,out_iaer)) |
|---|
| 1350 | #endif |
|---|
| 1351 | ENDIF ! out_qwg |
|---|
| 1352 | c================================================================== |
|---|
| 1353 | |
|---|
| 1354 | RETURN |
|---|
| 1355 | END |
|---|
