source: LMDZ5/branches/LMDZ5_AR5/libf/cosp/llnl_stats.F90 @ 3990

Last change on this file since 3990 was 1414, checked in by idelkadi, 14 years ago

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

File size: 5.9 KB
Line 
1! (c) 2008, Lawrence Livermore National Security Limited Liability Corporation.
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without modification, are permitted
5! provided that the following conditions are met:
6!
7!     * Redistributions of source code must retain the above copyright notice, this list
8!       of conditions and the following disclaimer.
9!     * Redistributions in binary form must reproduce the above copyright notice, this list
10!       of conditions and the following disclaimer in the documentation and/or other materials
11!       provided with the distribution.
12!     * Neither the name of the Lawrence Livermore National Security Limited Liability Corporation
13!       nor the names of its contributors may be used to endorse or promote products derived from
14!       this software without specific prior written permission.
15!
16! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
17! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
18! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
19! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
22! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
23! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25MODULE MOD_LLNL_STATS
26  USE MOD_COSP_CONSTANTS
27  IMPLICIT NONE
28
29CONTAINS
30
31!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32!-------------------- FUNCTION COSP_CFAD ------------------------
33!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34FUNCTION COSP_CFAD(Npoints,Ncolumns,Nlevels,Nbins,x,xmin,xmax,bmin,bwidth)
35   ! Input arguments
36   integer,intent(in) :: Npoints,Ncolumns,Nlevels,Nbins
37   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: x
38   real,intent(in) :: xmin,xmax
39   real,intent(in) :: bmin,bwidth
40   
41   real,dimension(Npoints,Nbins,Nlevels) :: cosp_cfad
42   ! Local variables
43   integer :: i, j, k
44   integer :: ibin
45   
46   !--- Input arguments
47   ! Npoints: Number of horizontal points
48   ! Ncolumns: Number of subcolumns
49   ! Nlevels: Number of levels
50   ! Nbins: Number of x axis bins
51   ! x: variable to process (Npoints,Ncolumns,Nlevels)
52   ! xmin: minimum value allowed for x
53   ! xmax: minimum value allowed for x
54   ! bmin: mimumum value of first bin
55   ! bwidth: bin width
56   !
57   ! Output: 2D histogram on each horizontal point (Npoints,Nbins,Nlevels)
58   
59   cosp_cfad = 0.0
60   ! bwidth intervals in the range [bmin,bmax=bmin+Nbins*hwidth]
61   ! Valid x values smaller than bmin and larger than bmax are set
62   ! into the smallest bin and largest bin, respectively.
63   do j = 1, Nlevels, 1
64      do k = 1, Ncolumns, 1
65         do i = 1, Npoints, 1
66            if (x(i,k,j) == R_GROUND) then
67               cosp_cfad(i,:,j) = R_UNDEF
68            elseif ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then
69               ibin = ceiling((x(i,k,j) - bmin)/bwidth)
70               if (ibin > Nbins) ibin = Nbins
71               if (ibin < 1)     ibin = 1
72               cosp_cfad(i,ibin,j) = cosp_cfad(i,ibin,j) + 1.0
73            end if
74         enddo  !i
75      enddo  !k
76   enddo  !j
77   where ((cosp_cfad /= R_UNDEF).and.(cosp_cfad /= 0.0)) cosp_cfad = cosp_cfad / Ncolumns
78END FUNCTION COSP_CFAD
79
80!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81!------------- SUBROUTINE COSP_LIDAR_ONLY_CLOUD -----------------
82!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83SUBROUTINE COSP_LIDAR_ONLY_CLOUD(Npoints,Ncolumns,Nlevels,beta_tot,beta_mol,Ze_tot,lidar_only_freq_cloud,tcc)
84   ! Input arguments
85   integer,intent(in) :: Npoints,Ncolumns,Nlevels
86   real,dimension(Npoints,Nlevels),intent(in) :: beta_mol   ! Molecular backscatter
87   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: beta_tot   ! Total backscattered signal
88   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: Ze_tot     ! Radar reflectivity
89   ! Output arguments
90   real,dimension(Npoints,Nlevels),intent(out) :: lidar_only_freq_cloud
91   real,dimension(Npoints),intent(out) :: tcc
92   
93   ! local variables
94   real :: sc_ratio
95   real :: s_cld, s_att
96!      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
97   parameter (S_cld = 5.0)  ! New (dec 2008) thresold for cloud detection
98   parameter (s_att = 0.01)
99   integer :: flag_sat !first saturated level encountered from top
100   integer :: flag_cld !cloudy column
101   integer :: pr,i,j
102   
103   lidar_only_freq_cloud = 0.0
104   tcc = 0.0
105   do pr=1,Npoints
106     do i=1,Ncolumns
107       flag_sat = 0
108       flag_cld = 0
109       do j=Nlevels,1,-1 !top->surf
110        sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
111!         if ((pr == 1).and.(j==8)) print *, pr,i,j,sc_ratio,Ze_tot(pr,i,j)
112        if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
113        if (Ze_tot(pr,i,j) .lt. -30.) then  !radar can't detect cloud
114         if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then  !lidar sense cloud
115!             if ((pr == 1).and.(j==8)) print *, 'L'
116            lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
117            flag_cld=1
118         endif
119        else  !radar sense cloud (z%Ze_tot(pr,i,j) .ge. -30.)
120!            if ((pr == 1).and.(j==8)) print *, 'R'
121           flag_cld=1
122        endif
123       enddo !levels
124       if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1.
125     enddo !columns
126!      if (tcc(pr) > Ncolumns) then
127!      print *, 'tcc(',pr,'): ', tcc(pr)
128!      tcc(pr) = Ncolumns
129!      endif
130   enddo !points
131   lidar_only_freq_cloud=lidar_only_freq_cloud/Ncolumns
132   tcc=tcc/Ncolumns
133
134END SUBROUTINE COSP_LIDAR_ONLY_CLOUD
135END MODULE MOD_LLNL_STATS
Note: See TracBrowser for help on using the repository browser.