Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/cosp_stats.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/cosp_stats.F90	(revision 1394)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/cosp_stats.F90	(revision 1395)
@@ -29,4 +29,5 @@
 ! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
 ! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
+! May 2010 - L. Fairhead      - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer
 !
 ! 
@@ -71,5 +72,5 @@
    Ncolumns = gbx%Ncolumns
    Nlr      = vgrid%Nlvgrid
-  
+
    if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
 
@@ -88,5 +89,5 @@
         !++++++++++++ Radar CFAD ++++++++++++++++
         if (cfg%Lradar_sim) then
-            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
+           call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
                                            Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
             stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
@@ -161,74 +162,93 @@
 
    ! Local variables
-   integer :: i,j,k
+   integer :: i,j,k,c
    logical :: lunits
-   real :: ws
-   real,dimension(Nlevels) :: xl,xu ! Lower and upper boundaries of model grid
+   real :: ws(Npoints,M), ws_temp(Npoints,M)
+   real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid
    real,dimension(M) :: dz          ! Layer depth
-   real,dimension(Nlevels,M) :: w   ! Weights to do the mean at each point
-   real,dimension(Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
+   real,dimension(Npoints,Nlevels,M) :: w   ! Weights to do the mean at each point
+   real,dimension(Npoints,Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
                                            ! Local copy at a particular point.
                                            ! This allows for change of units.
-   
+  
+
    lunits=.false.
    if (present(log_units)) lunits=log_units
-   
-   r = 0.0
-   do i=1,Npoints
-     ! Vertical grid at that point
-     xl = zhalf(i,:)
-     xu(1:Nlevels-1) = xl(2:Nlevels)
-     xu(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
-     dz = zu - zl
-     yp = y(i,:,:) ! Temporary variable to regrid
-     ! Find weights
-     w = 0.0
-     do k=1,M
-       do j=1,Nlevels
-         if ((xl(j) < zl(k)).and.(xu(j) > zl(k)).and.(xu(j) <= zu(k))) then
+
+   r(:,:,:) = 0.0
+   yp(:,:,:) = y(:,:,:)
+   w(:,:,:) = 0.0
+   ws(:,:) = 0.0
+   ws_temp(:,:) = 0.0
+   dz = zu - zl
+
+   ! Check for dBZ and change if necessary
+   if (lunits) then
+     where (y /= R_UNDEF)
+        yp = 10.0**(y/10.0)
+     end where
+   endif
+
+   ! Vertical grids
+   xl(:,:) = zhalf(:,:)
+   xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels)
+   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
+  ! Find weights
+   do k=1, M
+     do j=1, Nlevels
+       do i=1, Npoints
+         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
            !xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(j,k) = xu(j) - zl(k)
-         else if ((xl(j) >= zl(k)).and.(xu(j) <= zu(k))) then
+           w(i,j,k) = xu(i,j) - zl(k)
+         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
            !           xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(j,k) = xu(j) - xl(j)
-         else if ((xl(j) >= zl(k)).and.(xl(j) < zu(k)).and.(xu(j) >= zu(k))) then
+           w(i,j,k) =  xl(i,j+1) - xl(i,j)
+         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
            !                           xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(j,k) = zu(k) - xl(j)
-         else if ((xl(j) <= zl(k)).and.(xu(j) >= zu(k))) then
+           w(i,j,k) = zu(k) - xl(i,j)
+         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
            !  xl(j)---------------------------xu(j)
            !        zl(k)--------------zu(k)
-           w(j,k) = dz(j)
+           w(i,j,k) = dz(j)
          endif
        enddo
      enddo
-     ! Check for dBZ and change if necessary
-     if (lunits) then
-        where (yp /= R_UNDEF)
-          yp = 10.0**(yp/10.0)
-        elsewhere
-          yp = 0.0
-        end where
-     endif
-     ! Do the weighted mean
-     do j=1,Ncolumns
-       do k=1,M
-          ws = sum(w(:,k))
-          if (ws > 0.0) r(i,j,k) = sum(w(:,k)*yp(j,:))/ws
+   enddo
+
+   ! Do the weighted mean
+   do k=1,M
+     do i = 1, Npoints
+        ws(i,k) = sum(w(i,:,k))
+     enddo
+   enddo
+
+   ws_temp = 1.
+   where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:)
+
+   do c=1,Ncolumns
+     do k=1,M
+       do i = 1, Npoints
+          r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k)
        enddo
      enddo
-     ! Check for dBZ and change if necessary
-     if (lunits) then
-        where (r(i,:,:) <= 0.0)
-          r(i,:,:) = R_UNDEF
-        elsewhere
-          r(i,:,:) = 10.0*log10(r(i,:,:))
-        end where
-     endif
+   enddo
+
+   do k=1,M
+     do i = 1, Npoints
+        if (ws(i,k) <= 0.0) r(i,:,k)=0.0
+     enddo
    enddo
  
- 
+   ! Check for dBZ and change if necessary
+   if (lunits) then
+      where (r(:,:,:) <= 0.0)
+        r(:,:,:) = R_UNDEF
+      elsewhere
+        r(:,:,:) = 10.0*log10(r(:,:,:))
+      end where
+   endif
    
 END SUBROUTINE COSP_CHANGE_VERTICAL_GRID