Changeset 1414 for LMDZ4/trunk/libf/cosp/cosp_stats.F90

Timestamp:

Jul 15, 2010, 5:21:22 PM (14 years ago)

Author:

idelkadi

Message:

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

File:

• LMDZ4/trunk/libf/cosp/cosp_stats.F90 (modified) (12 diffs)

LMDZ4/trunk/libf/cosp/cosp_stats.F90

@@ -1 +1 @@
 ! (c) British Crown Copyright 2008, the Met Office.
 ! All rights reserved.
-! 
-! Redistribution and use in source and binary forms, with or without modification, are permitted 
+!
+! Redistribution and use in source and binary forms, with or without modification, are permitted
 ! provided that the following conditions are met:
-! 
-!     * Redistributions of source code must retain the above copyright notice, this list 
+!
+!     * Redistributions of source code must retain the above copyright notice, this list
 !       of conditions and the following disclaimer.
 !     * Redistributions in binary form must reproduce the above copyright notice, this list
-!       of conditions and the following disclaimer in the documentation and/or other materials 
+!       of conditions and the following disclaimer in the documentation and/or other materials
 !       provided with the distribution.
-!     * Neither the name of the Met Office nor the names of its contributors may be used 
-!       to endorse or promote products derived from this software without specific prior written 
+!     * Neither the name of the Met Office nor the names of its contributors may be used
+!       to endorse or promote products derived from this software without specific prior written
 !       permission.
-! 
-! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
-! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
-! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
-! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
-! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
-! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
-! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
+! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
+! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
@@ -29 +29 @@
 ! 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
-!
-! 
+! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
+!
+!
 MODULE MOD_COSP_STATS
   USE MOD_COSP_CONSTANTS
@@ -45 +45 @@
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
-  
+
    ! Input arguments
    type(cosp_gridbox),intent(in) :: gbx
@@ -55 +55 @@
    ! Output arguments
    type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
-   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar 
-   
-   ! Local variables 
+   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
+
+   ! Local variables
    integer :: Npoints  !# of grid points
    integer :: Nlevels  !# of levels
@@ -66 +66 @@
    real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
    real,dimension(:,:),allocatable :: ph_c,betamol_c
- 
+
    Npoints  = gbx%Npoints
    Nlevels  = gbx%Nlevels
@@ -73 +73 @@
    Nlr      = vgrid%Nlvgrid
 
-   if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
+   if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true.
 
    if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
@@ -81 +81 @@
         Ze_out = 0.0
         betatot_out  = 0.0
-        betamol_in(:,1,:) = sglidar%beta_mol(:,:)
         betamol_out= 0.0
         betamol_c  = 0.0
@@ -89 +88 @@
         !++++++++++++ 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, &
@@ -96 +95 @@
         !++++++++++++ Lidar CFAD ++++++++++++++++
         if (cfg%Llidar_sim) then
+            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
             call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
                                            Nlr,vgrid%zl,vgrid%zu,betamol_out)
@@ -114 +114 @@
         if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
                                     betatot_out,betamol_c,Ze_out, &
-                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
+                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
         ! Deallocate arrays at coarse resolution
         deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
@@ -131 +131 @@
         if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
                                     sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
-                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
+                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
    endif
    ! Replace undef
-   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF 
-   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF 
-   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF 
-   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF 
+   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
+   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
+   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
+   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
 
 END SUBROUTINE COSP_STATS
@@ -162 +162 @@
 
    ! Local variables
-   integer :: i,j,k,c
+   integer :: i,j,k
    logical :: lunits
-   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(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.
-  
+
+   integer :: l
+   real,dimension(Npoints) :: ws,sumwyp
+   real,dimension(Npoints,Nlevels) :: xl,xu
+   real,dimension(Npoints,Nlevels) :: w
+   real,dimension(Npoints,Ncolumns,Nlevels) :: yp
 
    lunits=.false.
    if (present(log_units)) lunits=log_units
 
-   r(:,:,:) = 0.0
-   yp(:,:,:) = y(:,:,:)
-   w(:,:,:) = 0.0
-   ws(:,:) = 0.0
-   ws_temp(:,:) = 0.0
-   dz = zu - zl
-
+   r(:,:,:) = R_GROUND
+   ! Vertical grid at that point
+   xl(:,:) = zhalf(:,:)
+   xu(:,1:Nlevels-1) = xl(:,2:Nlevels)
+   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
+   yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid
    ! Check for dBZ and change if necessary
    if (lunits) then
      where (y /= R_UNDEF)
-        yp = 10.0**(y/10.0)
+       yp = 10.0**(y/10.0)
+     elsewhere
+       yp = 0.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
+   do k=1,M
+     ! Find weights
+     w(:,:) = 0.0
+     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(i,j,k) = xu(i,j) - zl(k)
+           w(i,j) = 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(i,j,k) =  xl(i,j+1) - xl(i,j)
+           w(i,j) = xu(i,j) - 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(i,j,k) = zu(k) - xl(i,j)
+           w(i,j) = 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(i,j,k) = dz(j)
+           w(i,j) = zu(k) - zl(k)
+         endif
+       enddo
+     enddo
+     ! Do the weighted mean
+     do j=1,Ncolumns
+       ws    (:) = 0.0
+       sumwyp(:) = 0.0
+       do l=1,Nlevels
+         do i=1,Npoints
+           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+             ws    (i) = ws    (i) + w(i,l)
+             sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l)
+           endif
+         enddo
+       enddo
+       do i=1,Npoints
+         if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+           if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i)
          endif
        enddo
      enddo
    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
-   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
+     do k=1,M
+       do j=1,Ncolumns
+         do i=1,Npoints
+           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+             if (r(i,j,k) <= 0.0) then
+                 r(i,j,k) = R_UNDEF
+             else
+                 r(i,j,k) = 10.0*log10(r(i,j,k))
+             endif
+           endif
+         enddo
+       enddo
+     enddo
    endif
-   
-END SUBROUTINE COSP_CHANGE_VERTICAL_GRID 
+
+
+
+END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
 
 END MODULE MOD_COSP_STATS