Index: trunk/mars/libf/phymars/aeroptproperties.F
===================================================================
--- trunk/mars/libf/phymars/aeroptproperties.F	(revision 88)
+++ trunk/mars/libf/phymars/aeroptproperties.F	(revision 89)
@@ -37,5 +37,5 @@
 
 c     Min. and max radius of the interpolation grid (in METERS)
-      REAL, PARAMETER :: refftabmin = 5e-8 !2e-8
+      REAL, PARAMETER :: refftabmin = 1e-8 !5e-8
       REAL, PARAMETER :: refftabmax = 35e-6
 c     Log of the min and max variance of the interpolation grid
@@ -60,10 +60,22 @@
       REAL kint
       REAL drad
-      REAL weightgaus(5),radgaus(5)
+      INTEGER, PARAMETER :: ngau = 10
+      REAL weightgaus(ngau),radgaus(ngau)
       SAVE weightgaus,radgaus
-      DATA weightgaus/.2955242247,.2692667193,.2190863625,
-     &                .1494513491,.0666713443/
-      DATA radgaus/.1488743389,.4333953941,.6794095682,
-     &             .8650633666,.9739065285/
+c     DATA weightgaus/.2955242247,.2692667193,.2190863625,
+c    &                .1494513491,.0666713443/
+c     DATA radgaus/.1488743389,.4333953941,.6794095682,
+c    &             .8650633666,.9739065285/
+      DATA    radgaus/0.07652652113350,0.22778585114165,
+     &                0.37370608871528,0.51086700195146,
+     &                0.63605368072468,0.74633190646476,
+     &                0.83911697181213,0.91223442826796,
+     &                0.96397192726078,0.99312859919241/
+
+      DATA weightgaus/0.15275338723120,0.14917298659407,
+     &                0.14209610937519,0.13168863843930,
+     &                0.11819453196154,0.10193011980823,
+     &                0.08327674160932,0.06267204829828,
+     &                0.04060142982019,0.01761400714091/
 
 c     Indices
@@ -106,31 +118,31 @@
 c     Variables used by the Gauss-Legendre integration:
       REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
-      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,5)
-      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,5)
-
-      REAL,SAVE :: radGAUSa(5,naerkind,2)
-      REAL,SAVE :: radGAUSb(5,naerkind,2)
-
-      REAL,SAVE :: qsqrefVISa(nsun,5,naerkind)
-      REAL,SAVE :: qrefVISa(5,naerkind)
-      REAL,SAVE :: qsqrefVISb(nsun,5,naerkind)
-      REAL,SAVE :: qrefVISb(5,naerkind)
-      REAL,SAVE :: omegVISa(nsun,5,naerkind)
-      REAL,SAVE :: omegrefVISa(5,naerkind)
-      REAL,SAVE :: omegVISb(nsun,5,naerkind)
-      REAL,SAVE :: omegrefVISb(5,naerkind)
-      REAL,SAVE :: gVISa(nsun,5,naerkind)
-      REAL,SAVE :: gVISb(nsun,5,naerkind)
-
-      REAL,SAVE :: qsqrefIRa(nir,5,naerkind)
-      REAL,SAVE :: qrefIRa(5,naerkind)
-      REAL,SAVE :: qsqrefIRb(nir,5,naerkind)
-      REAL,SAVE :: qrefIRb(5,naerkind)
-      REAL,SAVE :: omegIRa(nir,5,naerkind)
-      REAL,SAVE :: omegrefIRa(5,naerkind)
-      REAL,SAVE :: omegIRb(nir,5,naerkind)
-      REAL,SAVE :: omegrefIRb(5,naerkind)
-      REAL,SAVE :: gIRa(nir,5,naerkind)
-      REAL,SAVE :: gIRb(nir,5,naerkind)
+      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
+      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
+
+      REAL,SAVE :: radGAUSa(ngau,naerkind,2)
+      REAL,SAVE :: radGAUSb(ngau,naerkind,2)
+
+      REAL,SAVE :: qsqrefVISa(nsun,ngau,naerkind)
+      REAL,SAVE :: qrefVISa(ngau,naerkind)
+      REAL,SAVE :: qsqrefVISb(nsun,ngau,naerkind)
+      REAL,SAVE :: qrefVISb(ngau,naerkind)
+      REAL,SAVE :: omegVISa(nsun,ngau,naerkind)
+      REAL,SAVE :: omegrefVISa(ngau,naerkind)
+      REAL,SAVE :: omegVISb(nsun,ngau,naerkind)
+      REAL,SAVE :: omegrefVISb(ngau,naerkind)
+      REAL,SAVE :: gVISa(nsun,ngau,naerkind)
+      REAL,SAVE :: gVISb(nsun,ngau,naerkind)
+
+      REAL,SAVE :: qsqrefIRa(nir,ngau,naerkind)
+      REAL,SAVE :: qrefIRa(ngau,naerkind)
+      REAL,SAVE :: qsqrefIRb(nir,ngau,naerkind)
+      REAL,SAVE :: qrefIRb(ngau,naerkind)
+      REAL,SAVE :: omegIRa(nir,ngau,naerkind)
+      REAL,SAVE :: omegrefIRa(ngau,naerkind)
+      REAL,SAVE :: omegIRb(nir,ngau,naerkind)
+      REAL,SAVE :: omegrefIRb(ngau,naerkind)
+      REAL,SAVE :: gIRa(nir,ngau,naerkind)
+      REAL,SAVE :: gIRb(nir,ngau,naerkind)
 
       REAL :: radiusm
@@ -162,4 +174,15 @@
       REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
       REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
+
+c     Tests
+c     -----
+
+      LOGICAL,SAVE :: out_qwg = .false.
+      INTEGER, PARAMETER :: out_iaer = 1
+      REAL :: out_qext(ngridmx,nlayermx)
+      REAL :: out_omeg(ngridmx,nlayermx)
+      REAL :: out_g(ngridmx,nlayermx)
+      INTEGER :: out_nchannel
+      CHARACTER*1 :: out_str
 
       DO iaer = 1, naerkind ! Loop on aerosol kind
@@ -225,5 +248,5 @@
 
 c       1.5 Interpolating data at the Gauss quadrature points:
-        DO gausind=1,5
+        DO gausind=1,ngau
           drad=radiusr*radgaus(gausind)
           radGAUSa(gausind,iaer,idomain)=radiusm-drad
@@ -286,5 +309,5 @@
         ENDDO
 
-        DO gausind=1,5
+        DO gausind=1,ngau
           drad=radiusr*radgaus(gausind)
           radGAUSb(gausind,iaer,idomain)=radiusm+drad
@@ -391,5 +414,5 @@
 
                 normd(j,1,iaer,idomain) = 1e-30
-                DO gausind=1,5
+                DO gausind=1,ngau
                   drad=radiusr*radgaus(gausind)
                   dista(j,1,iaer,idomain,gausind) =
@@ -435,5 +458,5 @@
                   omegrefVISgrid(j,1,iaer)=0.
 
-                  DO gausind=1,5
+                  DO gausind=1,ngau
                     DO m=1,nsun
 c                     Convolution:
@@ -551,5 +574,5 @@
                   omegrefIRgrid(j,1,iaer)=0.
 
-                  DO gausind=1,5
+                  DO gausind=1,ngau
                     DO m=1,nir
 c                     Convolution:
@@ -772,5 +795,5 @@
 
                 normd(j,k,iaer,idomain) = 1e-30
-                DO gausind=1,5
+                DO gausind=1,ngau
                   drad=radiusr*radgaus(gausind)
 
@@ -817,5 +840,5 @@
                   omegrefVISgrid(j,k,iaer)=0.
 
-                  DO gausind=1,5
+                  DO gausind=1,ngau
                     DO m=1,nsun
 c                     Convolution:
@@ -932,5 +955,5 @@
                   omegrefIRgrid(j,k,iaer)=0.
 
-                  DO gausind=1,5
+                  DO gausind=1,ngau
                     DO m=1,nir
 c                     Convolution:
@@ -1114,8 +1137,57 @@
       ENDDO ! iaer (loop on aerosol kind)
 
+c=====Radiative properties - TESTS=================================
+      IF (out_qwg) THEN
+      IF (ngrid.NE.1) THEN
+        DO out_nchannel = 1, 2
+c     -------------------------------------------------------------
+          DO lg = 1, nlayer
+            DO ig = 1, ngrid
+              out_qext(ig,lg) = 
+     &           QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
+     &           QREFvis3d(ig,lg,out_iaer)
+              out_omeg(ig,lg) =
+     &           omegaVIS3d(ig,lg,out_nchannel,out_iaer)
+              out_g(ig,lg)    = gVIS3d(ig,lg,out_nchannel,out_iaer)
+            ENDDO ! ig
+          ENDDO ! lg
+          write(out_str(1:1),'(i1.1)') out_nchannel
+          call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
+     &                     3,out_qext)
+          call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
+     &                     3,out_omeg)
+          call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
+     &                     3,out_g)
+c     -------------------------------------------------------------
+        ENDDO ! out_nchannel 
+        DO out_nchannel = 2, 4
+c     -------------------------------------------------------------
+          DO lg = 1, nlayer
+            DO ig = 1, ngrid
+              out_qext(ig,lg) =
+     &          QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
+     &          QREFir3d(ig,lg,out_iaer)
+              out_omeg(ig,lg) =
+     &          omegaIR3d(ig,lg,out_nchannel,out_iaer)
+              out_g(ig,lg)    = gIR3d(ig,lg,out_nchannel,out_iaer)
+            ENDDO ! ig
+          ENDDO ! lg
+          write(out_str(1:1),'(i1.1)') out_nchannel
+          call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
+     &                     3,out_qext)
+          call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
+     &                     3,out_omeg)
+          call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
+     &                     3,out_g)
+c     -------------------------------------------------------------
+        ENDDO ! out_nchannel 
+        call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
+     &                   3,omegaREFvis3d(1,1,out_iaer))
+        call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
+     &                   3,omegaREFir3d(1,1,out_iaer))
+      ENDIF ! ngrid.EQ.1
+      ENDIF ! out_qwg
+c==================================================================
+
       RETURN
       END
-
-
-
-
