Changeset 1396 for LMDZ4/trunk


Ignore:
Timestamp:
May 31, 2010, 11:29:51 AM (14 years ago)
Author:
Laurent Fairhead
Message:

Optimisation of routine for NEC SX8. Routine gives the same results as before.


Optimisation de la routine pour la NEC SX8. La nouvelle routine donne les
mêmes résultats que la précédente

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/trunk/libf/cosp/cosp_stats.F90

    r1279 r1396  
    2929! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
    3030! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
     31! May 2010 - L. Fairhead      - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer
    3132!
    3233!
     
    7172   Ncolumns = gbx%Ncolumns
    7273   Nlr      = vgrid%Nlvgrid
    73  
     74
    7475   if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
    7576
     
    8889        !++++++++++++ Radar CFAD ++++++++++++++++
    8990        if (cfg%Lradar_sim) then
    90             call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
     91           call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
    9192                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
    9293            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
     
    161162
    162163   ! Local variables
    163    integer :: i,j,k
     164   integer :: i,j,k,c
    164165   logical :: lunits
    165    real :: ws
    166    real,dimension(Nlevels) :: xl,xu ! Lower and upper boundaries of model grid
     166   real :: ws(Npoints,M), ws_temp(Npoints,M)
     167   real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid
    167168   real,dimension(M) :: dz          ! Layer depth
    168    real,dimension(Nlevels,M) :: w   ! Weights to do the mean at each point
    169    real,dimension(Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
     169   real,dimension(Npoints,Nlevels,M) :: w   ! Weights to do the mean at each point
     170   real,dimension(Npoints,Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
    170171                                           ! Local copy at a particular point.
    171172                                           ! This allows for change of units.
    172    
     173 
     174
    173175   lunits=.false.
    174176   if (present(log_units)) lunits=log_units
    175    
    176    r = 0.0
    177    do i=1,Npoints
    178      ! Vertical grid at that point
    179      xl = zhalf(i,:)
    180      xu(1:Nlevels-1) = xl(2:Nlevels)
    181      xu(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
    182      dz = zu - zl
    183      yp = y(i,:,:) ! Temporary variable to regrid
    184      ! Find weights
    185      w = 0.0
    186      do k=1,M
    187        do j=1,Nlevels
    188          if ((xl(j) < zl(k)).and.(xu(j) > zl(k)).and.(xu(j) <= zu(k))) then
     177
     178   r(:,:,:) = 0.0
     179   yp(:,:,:) = y(:,:,:)
     180   w(:,:,:) = 0.0
     181   ws(:,:) = 0.0
     182   ws_temp(:,:) = 0.0
     183   dz = zu - zl
     184
     185   ! Check for dBZ and change if necessary
     186   if (lunits) then
     187     where (y /= R_UNDEF)
     188        yp = 10.0**(y/10.0)
     189     end where
     190   endif
     191
     192   ! Vertical grids
     193   xl(:,:) = zhalf(:,:)
     194   xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels)
     195   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
     196  ! Find weights
     197   do k=1, M
     198     do j=1, Nlevels
     199       do i=1, Npoints
     200         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
    189201           !xl(j)-----------------xu(j)
    190202           !      zl(k)------------------------------zu(k)
    191            w(j,k) = xu(j) - zl(k)
    192          else if ((xl(j) >= zl(k)).and.(xu(j) <= zu(k))) then
     203           w(i,j,k) = xu(i,j) - zl(k)
     204         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
    193205           !           xl(j)-----------------xu(j)
    194206           !      zl(k)------------------------------zu(k)
    195            w(j,k) = xu(j) - xl(j)
    196          else if ((xl(j) >= zl(k)).and.(xl(j) < zu(k)).and.(xu(j) >= zu(k))) then
     207           w(i,j,k) =  xl(i,j+1) - xl(i,j)
     208         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
    197209           !                           xl(j)-----------------xu(j)
    198210           !      zl(k)------------------------------zu(k)
    199            w(j,k) = zu(k) - xl(j)
    200          else if ((xl(j) <= zl(k)).and.(xu(j) >= zu(k))) then
     211           w(i,j,k) = zu(k) - xl(i,j)
     212         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
    201213           !  xl(j)---------------------------xu(j)
    202214           !        zl(k)--------------zu(k)
    203            w(j,k) = dz(j)
     215           w(i,j,k) = dz(j)
    204216         endif
    205217       enddo
    206218     enddo
    207      ! Check for dBZ and change if necessary
    208      if (lunits) then
    209         where (yp /= R_UNDEF)
    210           yp = 10.0**(yp/10.0)
    211         elsewhere
    212           yp = 0.0
    213         end where
    214      endif
    215      ! Do the weighted mean
    216      do j=1,Ncolumns
    217        do k=1,M
    218           ws = sum(w(:,k))
    219           if (ws > 0.0) r(i,j,k) = sum(w(:,k)*yp(j,:))/ws
     219   enddo
     220
     221   ! Do the weighted mean
     222   do k=1,M
     223     do i = 1, Npoints
     224        ws(i,k) = sum(w(i,:,k))
     225     enddo
     226   enddo
     227
     228   ws_temp = 1.
     229   where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:)
     230
     231   do c=1,Ncolumns
     232     do k=1,M
     233       do i = 1, Npoints
     234          r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k)
    220235       enddo
    221236     enddo
    222      ! Check for dBZ and change if necessary
    223      if (lunits) then
    224         where (r(i,:,:) <= 0.0)
    225           r(i,:,:) = R_UNDEF
    226         elsewhere
    227           r(i,:,:) = 10.0*log10(r(i,:,:))
    228         end where
    229      endif
     237   enddo
     238
     239   do k=1,M
     240     do i = 1, Npoints
     241        if (ws(i,k) <= 0.0) r(i,:,k)=0.0
     242     enddo
    230243   enddo
    231244 
    232  
     245   ! Check for dBZ and change if necessary
     246   if (lunits) then
     247      where (r(:,:,:) <= 0.0)
     248        r(:,:,:) = R_UNDEF
     249      elsewhere
     250        r(:,:,:) = 10.0*log10(r(:,:,:))
     251      end where
     252   endif
    233253   
    234254END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
Note: See TracChangeset for help on using the changeset viewer.