source: LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/llnl_stats.F90 @ 1329

Last change on this file since 1329 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

File size: 5.7 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  IMPLICIT NONE
27
28CONTAINS
29
30!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31!-------------------- FUNCTION COSP_CFAD ------------------------
32!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33FUNCTION COSP_CFAD(Npoints,Ncolumns,Nlevels,Nbins,x,xmin,xmax,bmin,bwidth)
34   ! Input arguments
35   integer,intent(in) :: Npoints,Ncolumns,Nlevels,Nbins
36   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: x
37   real,intent(in) :: xmin,xmax
38   real,intent(in) :: bmin,bwidth
39   
40   real,dimension(Npoints,Nbins,Nlevels) :: cosp_cfad
41   ! Local variables
42   integer :: i, j, k
43   integer :: ibin
44   
45   !--- Input arguments
46   ! Npoints: Number of horizontal points
47   ! Ncolumns: Number of subcolumns
48   ! Nlevels: Number of levels
49   ! Nbins: Number of x axis bins
50   ! x: variable to process (Npoints,Ncolumns,Nlevels)
51   ! xmin: minimum value allowed for x
52   ! xmax: minimum value allowed for x
53   ! bmin: mimumum value of first bin
54   ! bwidth: bin width
55   !
56   ! Output: 2D histogram on each horizontal point (Npoints,Nbins,Nlevels)
57   
58   cosp_cfad = 0.0
59   ! bwidth intervals in the range [bmin,bmax=bmin+Nbins*hwidth]
60   ! Valid x values smaller than bmin and larger than bmax are set
61   ! into the smallest bin and largest bin, respectively.
62   do j = 1, Nlevels, 1
63      do k = 1, Ncolumns, 1
64         do i = 1, Npoints, 1
65            if ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then
66               ibin = ceiling((x(i,k,j) - bmin)/bwidth)
67               if (ibin > Nbins) ibin = Nbins
68               if (ibin < 1)     ibin = 1
69               cosp_cfad(i,ibin,j) = cosp_cfad(i,ibin,j) + 1.0
70            end if
71         enddo  !i
72      enddo  !k
73   enddo  !j
74   cosp_cfad = cosp_cfad / Ncolumns
75END FUNCTION COSP_CFAD
76
77!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78!------------- SUBROUTINE COSP_LIDAR_ONLY_CLOUD -----------------
79!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80SUBROUTINE COSP_LIDAR_ONLY_CLOUD(Npoints,Ncolumns,Nlevels,beta_tot,beta_mol,Ze_tot,lidar_only_freq_cloud,tcc)
81   ! Input arguments
82   integer,intent(in) :: Npoints,Ncolumns,Nlevels
83   real,dimension(Npoints,Nlevels),intent(in) :: beta_mol   ! Molecular backscatter
84   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: beta_tot   ! Total backscattered signal
85   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: Ze_tot     ! Radar reflectivity
86   ! Output arguments
87   real,dimension(Npoints,Nlevels),intent(out) :: lidar_only_freq_cloud
88   real,dimension(Npoints),intent(out) :: tcc
89   
90   ! local variables
91   real :: sc_ratio
92   real :: s_cld, s_att
93!      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
94   parameter (S_cld = 5.0)  ! New (dec 2008) thresold for cloud detection
95   parameter (s_att = 0.01)
96   integer :: flag_sat !first saturated level encountered from top
97   integer :: flag_cld !cloudy column
98   integer :: pr,i,j
99   
100!    lidar_only_freq_cloud = 0.0
101!    tcc = 0.0
102   do pr=1,Npoints
103     do i=1,Ncolumns
104       flag_sat = 0
105       flag_cld = 0
106       do j=Nlevels,1,-1 !top->surf
107        sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
108!         if ((pr == 1).and.(j==8)) print *, pr,i,j,sc_ratio,Ze_tot(pr,i,j)
109        if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
110        if (Ze_tot(pr,i,j) .lt. -30.) then  !radar can't detect cloud
111         if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then  !lidar sense cloud
112!             if ((pr == 1).and.(j==8)) print *, 'L'
113            lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
114            flag_cld=1
115         endif
116        else  !radar sense cloud (z%Ze_tot(pr,i,j) .ge. -30.)
117!            if ((pr == 1).and.(j==8)) print *, 'R'
118           flag_cld=1
119        endif
120       enddo !levels
121       if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1.
122     enddo !columns
123!      if (tcc(pr) > Ncolumns) then
124!      print *, 'tcc(',pr,'): ', tcc(pr)
125!      tcc(pr) = Ncolumns
126!      endif
127   enddo !points
128   lidar_only_freq_cloud=lidar_only_freq_cloud/Ncolumns
129   tcc=tcc/Ncolumns
130
131END SUBROUTINE COSP_LIDAR_ONLY_CLOUD
132END MODULE MOD_LLNL_STATS
Note: See TracBrowser for help on using the repository browser.