Changeset 737


Ignore:
Timestamp:
Jul 24, 2012, 6:55:41 PM (12 years ago)
Author:
jbmadeleine
Message:

Corrected aeroptproperties.F to avoid having a floating underflow in
aeroptproperties.F due to a too small minimum radius of the
interpolation grid. This can lead to wrong values of the scattering
parameters or even NaNs?.

refftabmin and refftabmax are now computed automatically.

Previously, the user had to specify them, which is dangerous for the
above mentioned reason and not flexible in the long run.

Added two dimensions (iaer and idomain) to refftab and nuefftab

so that the grid can be changed depending on the range of radiuses
contained in the lookup tables for each aerosol.

Also added two dimensions to logvratgrid and delete vratgrid, which

is only used once.

Added a message at firstcall giving the range of radiuses used for

the interpolation grid or each aerosol.

normd is now checked to make sure that it is not equal to zero (in

which case there is a division by zero at the end of the routine).
This should not occur, but we never know...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/aeroptproperties.F

    r690 r737  
    3737
    3838c     Min. and max radius of the interpolation grid (in METERS)
    39       REAL, PARAMETER :: refftabmin = 1e-8 !5e-8
    40       REAL, PARAMETER :: refftabmax = 35e-6
     39      REAL, SAVE :: refftabmin(naerkind,2)
     40      REAL, SAVE :: refftabmax(naerkind,2)
    4141c     Log of the min and max variance of the interpolation grid
    4242      REAL, PARAMETER :: nuefftabmin = -4.6
    … …  
    8787
    8888c     Radius axis of the interpolation grid
    89       REAL,SAVE :: refftab(refftabsize)
     89      REAL,SAVE :: refftab(refftabsize,naerkind,2)
    9090c     Variance axis of the interpolation grid
    91       REAL,SAVE :: nuefftab(nuefftabsize)
     91      REAL,SAVE :: nuefftab(nuefftabsize,naerkind,2)
    9292c     Volume ratio of the grid
    93       REAL,SAVE :: logvratgrid,vratgrid
     93      REAL,SAVE :: logvratgrid(naerkind,2)
    9494c     Grid used to remember which calculation is done
    9595      LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2)
    … …  
    186186      CHARACTER*1 :: out_str
    187187
     188c     Creating the effective radius and variance grid
     189c     -----------------------------------------------
     190      IF (firstcall) THEN
     191c       0.1 Pi!
     192        pi = 2. * asin(1.e0)
     193
     194        WRITE(*,*) "aeroptproperties: interpolation grid"
     195        DO iaer = 1, naerkind ! Loop on aerosol kind
     196          DO idomain = 1, 2 ! Loop on visible or infrared channel
     197
     198c           0.2 Effective radius
     199            radiusm=
     200     &             0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))+
     201     &                  radiustab(iaer,idomain,1))
     202            radiusr=
     203     &             0.5*(radiustab(iaer,idomain,nsize(iaer,idomain))-
     204     &                  radiustab(iaer,idomain,1))
     205            refftabmin(iaer,idomain) =
     206     &        radiusm-radiusr*radgaus(ngau)
     207            refftabmax(iaer,idomain) =
     208     &        radiustab(iaer,idomain,nsize(iaer,idomain))
     209
     210            WRITE(*,*) "Scatterer: ",iaer
     211            WRITE(*,*) "Domain: ",idomain
     212            WRITE(*,*) "Min radius (m): ", refftabmin(iaer,idomain)
     213            WRITE(*,*) "Max radius (m): ", refftabmax(iaer,idomain)
     214
     215            refftab(1,iaer,idomain) =
     216     &        refftabmin(iaer,idomain)
     217            refftab(refftabsize,iaer,idomain) =
     218     &        refftabmax(iaer,idomain)
     219
     220            logvratgrid(iaer,idomain) =
     221     &                    log(refftabmax(iaer,idomain)/
     222     &                      refftabmin(iaer,idomain)) /
     223     &                    float(refftabsize-1)*3.
     224            do i = 2, refftabsize-1
     225              refftab(i,iaer,idomain) = refftab(i-1,iaer,idomain)*
     226     &                             exp(1./3.*logvratgrid(iaer,idomain))
     227            enddo
     228
     229c           0.3 Effective variance
     230            do i = 0, nuefftabsize-1
     231              nuefftab(i+1,iaer,idomain) = exp( nuefftabmin +
     232     &                 i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
     233            enddo
     234
     235          ENDDO
     236        ENDDO
     237        firstcall = .false.
     238      ENDIF
     239
    188240      DO iaer = 1, naerkind ! Loop on aerosol kind
    189241        IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
    … …  
    212264      DO idomain = 1, 2 ! Loop on visible or infrared channel
    213265c==================================================================
    214 c     1. Creating the effective radius and variance grid
    215 c     --------------------------------------------------
    216       IF (firstcall) THEN
    217 
    218 c       1.1 Pi!
    219         pi = 2. * asin(1.e0)
    220 
    221 c       1.2 Effective radius
    222         refftab(1)    = refftabmin
    223         refftab(refftabsize) = refftabmax
    224 
    225         logvratgrid = log(refftabmax/refftabmin) /
    226      &                float(refftabsize-1)*3.
    227         vratgrid = exp(logvratgrid)
    228 
    229         do i = 2, refftabsize-1
    230           refftab(i) = refftab(i-1)*vratgrid**(1./3.)
    231         enddo
    232 
    233 c       1.3 Effective variance
    234         do i = 0, nuefftabsize-1
    235           nuefftab(i+1) = exp( nuefftabmin +
    236      &             i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
    237         enddo
    238         firstcall = .false.
    239       ENDIF
    240 
    241 c       1.4 Radius middle point and range for Gauss integration
     266
     267c       1.1 Radius middle point and range for Gauss integration
    242268        radiusm=
    243269     &         0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
    … …  
    247273     &              radiustab(iaer,idomain,1))
    248274
    249 c       1.5 Interpolating data at the Gauss quadrature points:
     275c       1.2 Interpolating data at the Gauss quadrature points:
    250276        DO gausind=1,ngau
    251277          drad=radiusr*radgaus(gausind)
    … …  
    379405        DO ig = 1,ngrid
    380406c         2.1 Effective radius index and kx calculation
    381           var_tmp=reffrad(ig,lg,iaer)/refftabmin
     407          var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
    382408          var_tmp=log(var_tmp)*3.
    383           var_tmp=var_tmp/logvratgrid+1.
     409          var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
    384410          grid_i=floor(var_tmp)
    385411          IF (grid_i.GE.refftabsize) THEN
    … …  
    398424            kx = 0.
    399425          ELSE
    400             kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
    401      &           ( refftab(grid_i+1)-refftab(grid_i) )
     426            kx = ( reffrad(ig,lg,iaer)-
     427     &             refftab(grid_i,iaer,idomain) ) /
     428     &           ( refftab(grid_i+1,iaer,idomain)-
     429     &             refftab(grid_i,iaer,idomain) )
    402430          ENDIF
    403431c         2.3 Integration
    … …  
    411439                sizedistk2 = log(1.+nueffrad(1,1,iaer))
    412440                sizedistk1 = exp(2.5*sizedistk2)
    413                 sizedistk1 = refftab(j) / sizedistk1
     441                sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
    414442
    415443                normd(j,1,iaer,idomain) = 1e-30
    … …  
    447475     &              )
    448476                ENDDO
     477                IF (normd(j,1,iaer,idomain).EQ.1e-30) THEN
     478                  WRITE(*,*)"normd:", normd(j,1,iaer,idomain)
     479                  WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)"
     480                  WRITE(*,*)"Check the size of the interpolation grid."
     481                  CALL ABORT
     482                ENDIF
    449483                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
    450484c                 2.3.3.vis Initialization
    … …  
    754788            ky = 0.
    755789          ELSE
    756             ky = ( nueffrad(ig,lg,iaer)-nuefftab(grid_j) ) /
    757      &           ( nuefftab(grid_j+1)-nuefftab(grid_j) )
     790            ky = ( nueffrad(ig,lg,iaer)-
     791     &             nuefftab(grid_j,iaer,idomain) ) /
     792     &           ( nuefftab(grid_j+1,iaer,idomain)-
     793     &             nuefftab(grid_j,iaer,idomain) )
    758794          ENDIF
    759795c         3.2 Effective radius index and kx calculation
    760           var_tmp=reffrad(ig,lg,iaer)/refftabmin
     796          var_tmp=reffrad(ig,lg,iaer)/refftabmin(iaer,idomain)
    761797          var_tmp=log(var_tmp)*3.
    762           var_tmp=var_tmp/logvratgrid+1.
     798          var_tmp=var_tmp/logvratgrid(iaer,idomain)+1.
    763799          grid_i=floor(var_tmp)
    764800          IF (grid_i.GE.refftabsize) THEN
    … …  
    777813            kx = 0.
    778814          ELSE
    779             kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
    780      &           ( refftab(grid_i+1)-refftab(grid_i) )
     815            kx = ( reffrad(ig,lg,iaer)-
     816     &             refftab(grid_i,iaer,idomain) ) /
     817     &           ( refftab(grid_i+1,iaer,idomain)-
     818     &             refftab(grid_i,iaer,idomain) )
    781819          ENDIF
    782820c         3.3 Integration
    … …  
    790828c                 atmospheres", Space Science Reviews 16 527-610.
    791829c                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
    792                 sizedistk2 = log(1.+nuefftab(k))
     830                sizedistk2 = log(1.+nuefftab(k,iaer,idomain))
    793831                sizedistk1 = exp(2.5*sizedistk2)
    794                 sizedistk1 = refftab(j) / sizedistk1
     832                sizedistk1 = refftab(j,iaer,idomain) / sizedistk1
    795833
    796834                normd(j,k,iaer,idomain) = 1e-30
    … …  
    829867     &              )
    830868                ENDDO
     869                IF (normd(j,k,iaer,idomain).EQ.1e-30) THEN
     870                  WRITE(*,*)"normd:", normd(j,k,iaer,idomain)
     871                  WRITE(*,*)"Risk of division by 0 (aeroptproperties.F)"
     872                  WRITE(*,*)"Check the size of the interpolation grid."
     873                  CALL ABORT
     874                ENDIF
    831875                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
    832876c                 2.3.3.vis Initialization
Note: See TracChangeset for help on using the changeset viewer.