source: LMDZ6/trunk/libf/phylmd/cosp2/cosp_stats.F90 @ 5069

Last change on this file since 5069 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"

File size: 11.3 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2015, Regents of the University of Colorado
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History:
30! Jul 2007 - A. Bodas-Salcedo - Initial version
31! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid
32! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns
33!                               in COSP_STATS
34! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
35! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
36! Jan 2013 - G. Cesana        - Added betaperp and temperature arguments
37!                             - Added phase 3D/3Dtemperature/Map output variables in diag_lidar
38! May 2015 - D. Swales        - Modified for cosp2.0
39! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
40MODULE MOD_COSP_STATS
41  USE COSP_KINDS, ONLY: wp
42  USE MOD_COSP_CONFIG, ONLY: R_UNDEF,R_GROUND
43  IMPLICIT NONE
44CONTAINS
45
46  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47  !---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
48  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
50   implicit none
51   ! Input arguments
52   integer,intent(in) :: Npoints  !# of grid points
53   integer,intent(in) :: Nlevels  !# of levels
54   integer,intent(in) :: Ncolumns !# of columns
55   real(wp),dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
56   real(wp),dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
57   real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
58   integer,intent(in) :: Nglevels  !# levels in the new grid
59   real(wp),dimension(Nglevels),intent(in) :: newgrid_bot ! Lower boundary of new levels  [m]
60   real(wp),dimension(Nglevels),intent(in) :: newgrid_top ! Upper boundary of new levels  [m]
61   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
62   ! Output
63   real(wp),dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
64
65   ! Local variables
66   integer :: i,j,k
67   logical :: lunits
68   integer :: l
69   real(wp) :: w ! Weight
70   real(wp) :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
71   integer :: Nw  ! Number of weights
72   real(wp) :: wt  ! Sum of weights
73   real(wp),dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
74   real(wp) :: yp ! Local copy of y at a particular point.
75              ! This allows for change of units.
76
77   lunits=.false.
78   if (present(log_units)) lunits=log_units
79
80   r = 0._wp
81
82   do i=1,Npoints
83     ! Calculate tops and bottoms of new and old grids
84     oldgrid_bot = zhalf(i,:)
85     oldgrid_top(1:Nlevels-1) = oldgrid_bot(2:Nlevels)
86     oldgrid_top(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
87     l = 0 ! Index of level in the old grid
88     ! Loop over levels in the new grid
89     do k = 1,Nglevels
90       Nw = 0 ! Number of weigths
91       wt = 0._wp ! Sum of weights
92       ! Loop over levels in the old grid and accumulate total for weighted average
93       do
94         l = l + 1
95         w = 0.0 ! Initialise weight to 0
96         ! Distances between edges of both grids
97         dbb = oldgrid_bot(l) - newgrid_bot(k)
98         dtb = oldgrid_top(l) - newgrid_bot(k)
99         dbt = oldgrid_bot(l) - newgrid_top(k)
100         dtt = oldgrid_top(l) - newgrid_top(k)
101         if (dbt >= 0.0) exit ! Do next level in the new grid
102         if (dtb > 0.0) then
103           if (dbb <= 0.0) then
104             if (dtt <= 0) then
105               w = dtb
106             else
107               w = newgrid_top(k) - newgrid_bot(k)
108             endif
109           else
110             if (dtt <= 0) then
111               w = oldgrid_top(l) - oldgrid_bot(l)
112             else
113               w = -dbt
114             endif
115           endif
116           ! If layers overlap (w/=0), then accumulate
117           if (w /= 0.0) then
118             Nw = Nw + 1
119             wt = wt + w
120             do j=1,Ncolumns
121               if (lunits) then
122                 if (y(i,j,l) /= R_UNDEF) then
123                   yp = 10._wp**(y(i,j,l)/10._wp)
124                 else
125                   yp = 0._wp
126                 endif
127               else
128                 yp = y(i,j,l)
129               endif
130               r(i,j,k) = r(i,j,k) + w*yp
131             enddo
132           endif
133         endif
134       enddo
135       l = l - 2
136       if (l < 1) l = 0
137       ! Calculate average in new grid
138       if (Nw > 0) then
139         do j=1,Ncolumns
140           r(i,j,k) = r(i,j,k)/wt
141         enddo
142       endif
143     enddo
144   enddo
145
146   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
147   do k=1,Nglevels
148     do j=1,Ncolumns
149       do i=1,Npoints
150         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
151           if (lunits) then
152             if (r(i,j,k) <= 0.0) then
153               r(i,j,k) = R_UNDEF
154             else
155               r(i,j,k) = 10._wp*log10(r(i,j,k))
156             endif
157           endif
158         else ! Level below surface
159           r(i,j,k) = R_GROUND
160         endif
161       enddo
162     enddo
163   enddo
164
165END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
166
167  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
168  !------------- SUBROUTINE COSP_LIDAR_ONLY_CLOUD -----------------
169  ! (c) 2008, Lawrence Livermore National Security Limited Liability Corporation.
170  ! All rights reserved.
171  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
172  SUBROUTINE COSP_LIDAR_ONLY_CLOUD(Npoints,Ncolumns,Nlevels,beta_tot, &
173                                   beta_mol,Ze_tot,lidar_only_freq_cloud,tcc)
174    ! Inputs
175    integer,intent(in) :: &
176         Npoints,       & ! Number of horizontal gridpoints
177         Ncolumns,      & ! Number of subcolumns
178         Nlevels          ! Number of vertical layers
179    real(wp),dimension(Npoints,Nlevels),intent(in) :: &
180         beta_mol         ! Molecular backscatter
181    real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: &
182         beta_tot,      & ! Total backscattered signal
183         Ze_tot           ! Radar reflectivity
184    ! Outputs
185    real(wp),dimension(Npoints,Nlevels),intent(out) :: &
186         lidar_only_freq_cloud
187    real(wp),dimension(Npoints),intent(out) ::&
188         tcc
189   
190    ! local variables
191    real(wp) :: sc_ratio
192    real(wp),parameter :: &
193         s_cld=5.0, &
194         s_att=0.01
195    integer :: flag_sat,flag_cld,pr,i,j
196   
197    lidar_only_freq_cloud = 0._wp
198    tcc = 0._wp
199    do pr=1,Npoints
200       do i=1,Ncolumns
201          flag_sat = 0
202          flag_cld = 0
203          do j=1,Nlevels
204             sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
205             if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
206             if (Ze_tot(pr,i,j) .lt. -30.) then  !radar can't detect cloud
207                if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then  !lidar sense cloud
208                   lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
209                   flag_cld=1
210                endif
211             else  !radar sense cloud (z%Ze_tot(pr,i,j) .ge. -30.)
212                flag_cld=1
213             endif
214          enddo !levels
215          if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1._wp
216       enddo !columns
217    enddo !points
218    lidar_only_freq_cloud=lidar_only_freq_cloud/Ncolumns
219    tcc=tcc/Ncolumns
220   
221    ! Unit conversion
222    where(lidar_only_freq_cloud /= R_UNDEF) &
223            lidar_only_freq_cloud = lidar_only_freq_cloud*100._wp
224    where(tcc /= R_UNDEF) tcc = tcc*100._wp
225   
226  END SUBROUTINE COSP_LIDAR_ONLY_CLOUD
227 
228  ! ######################################################################################
229  ! FUNCTION hist1D
230  ! ######################################################################################
231  function hist1d(Npoints,var,nbins,bins)
232    ! Inputs
233    integer,intent(in) :: &
234         Npoints, & ! Number of points in input array
235         Nbins      ! Number of bins for sorting
236    real(wp),intent(in),dimension(Npoints) :: &
237         var        ! Input variable to be sorted
238    real(wp),intent(in),dimension(Nbins+1) :: &
239         bins       ! Histogram bins [lowest,binTops] 
240    ! Outputs
241    real(wp),dimension(Nbins) :: &
242         hist1d     ! Output histogram     
243    ! Local variables
244    integer :: ij
245   
246    do ij=2,Nbins+1 
247       hist1D(ij-1) = count(var .ge. bins(ij-1) .and. var .lt. bins(ij))
248       if (count(var .eq. R_GROUND) .ge. 1) hist1D(ij-1)=R_UNDEF
249    enddo
250   
251  end function hist1D
252 
253  ! ######################################################################################
254  ! SUBROUTINE hist2D
255  ! ######################################################################################
256  subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
257    implicit none
258   
259    ! INPUTS
260    integer, intent(in) :: &
261         npts,  & ! Number of data points to be sorted
262         nbin1, & ! Number of bins in histogram direction 1
263         nbin2    ! Number of bins in histogram direction 2
264    real(wp),intent(in),dimension(npts) :: &
265         var1,  & ! Variable 1 to be sorted into bins
266         var2     ! variable 2 to be sorted into bins
267    real(wp),intent(in),dimension(nbin1+1) :: &
268         bin1     ! Histogram bin 1 boundaries
269    real(wp),intent(in),dimension(nbin2+1) :: &
270         bin2     ! Histogram bin 2 boundaries
271    ! OUTPUTS
272    real(wp),intent(out),dimension(nbin1,nbin2) :: &
273         jointHist
274   
275    ! LOCAL VARIABLES
276    integer :: ij,ik
277   
278    do ij=2,nbin1+1
279       do ik=2,nbin2+1
280          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
281               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))       
282       enddo
283    enddo
284  end subroutine hist2D
285END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.