source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/cospv2/cosp_stats.F90 @ 5134

Last change on this file since 5134 was 3491, checked in by idelkadi, 6 years ago

Integration of version 2 of the COSP simulator in LMDZ
This line, and those below, will be ignored--

M makegcm
M makelmdz
M makelmdz_fcm
M libf/phylmd/physiq_mod.F90
A libf/phylmd/cospv2
A libf/phylmd/cospv2/mo_rng.F90
A libf/phylmd/cospv2/quickbeam_optics.F90
A libf/phylmd/cospv2/cosp_cloudsat_interface.F90
A libf/phylmd/cospv2/cosp_config.F90
A libf/phylmd/cospv2/lidar_simulator.F90
A libf/phylmd/cospv2/prec_scops.F90
A libf/phylmd/cospv2/mrgrnk.F90
A libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
A libf/phylmd/cospv2/cosp_atlid_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
A libf/phylmd/cospv2/cosp_math_constants.F90
A libf/phylmd/cospv2/MISR_simulator.F90
A libf/phylmd/cospv2/modis_simulator.F90
A libf/phylmd/cospv2/math_lib.F90
A libf/phylmd/cospv2/cosp_grLidar532_interface.F90
A libf/phylmd/cospv2/cosp_errorHandling.F90
A libf/phylmd/cospv2/cosp_stats.F90
A libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
A libf/phylmd/cospv2/cosp_utils.F90
A libf/phylmd/cospv2/cosp_optics.F90
A libf/phylmd/cospv2/icarus.F90
A libf/phylmd/cospv2/scops.F90
A libf/phylmd/cospv2/optics_lib.F90
A libf/phylmd/cospv2/cosp_kinds.F90
A libf/phylmd/cospv2/cosp_calipso_interface.F90
A libf/phylmd/cospv2/quickbeam.F90
A libf/phylmd/cospv2/parasol.F90
A libf/phylmd/cospv2/cosp_phys_constants.F90
A libf/phylmd/cospv2/cosp.F90
A libf/phylmd/cospv2/array_lib.F90
A libf/phylmd/cospv2/cosp_isccp_interface.F90
A libf/phylmd/cospv2/cosp_parasol_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_interface.F90
A libf/phylmd/cospv2/cosp_misr_interface.F90
A libf/phylmd/cospv2/cosp_modis_interface.F90

File size: 12.2 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
44  IMPLICIT NONE
45CONTAINS
46
47  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48  !---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
49  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
51   implicit none
52   ! Input arguments
53   integer,intent(in) :: Npoints  !# of grid points
54   integer,intent(in) :: Nlevels  !# of levels
55   integer,intent(in) :: Ncolumns !# of columns
56   real(wp),dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
57   real(wp),dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
58   real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
59   integer,intent(in) :: Nglevels  !# levels in the new grid
60   real(wp),dimension(Nglevels),intent(in) :: newgrid_bot ! Lower boundary of new levels  [m]
61   real(wp),dimension(Nglevels),intent(in) :: newgrid_top ! Upper boundary of new levels  [m]
62   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
63   ! Output
64   real(wp),dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
65
66   ! Local variables
67   integer :: i,j,k
68   logical :: lunits
69   integer :: l
70   real(wp) :: w ! Weight
71   real(wp) :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
72   integer :: Nw  ! Number of weights
73   real(wp) :: wt  ! Sum of weights
74   real(wp),dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
75   real(wp) :: yp ! Local copy of y at a particular point.
76              ! This allows for change of units.
77
78   lunits=.false.
79   if (present(log_units)) lunits=log_units
80
81   r = 0._wp
82
83   do i=1,Npoints
84     ! Calculate tops and bottoms of new and old grids
85     oldgrid_bot = zhalf(i,:)
86     oldgrid_top(1:Nlevels-1) = oldgrid_bot(2:Nlevels)
87     oldgrid_top(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
88     l = 0 ! Index of level in the old grid
89     ! Loop over levels in the new grid
90     do k = 1,Nglevels
91       Nw = 0 ! Number of weigths
92       wt = 0._wp ! Sum of weights
93       ! Loop over levels in the old grid and accumulate total for weighted average
94       do
95         l = l + 1
96         w = 0.0 ! Initialise weight to 0
97         ! Distances between edges of both grids
98         dbb = oldgrid_bot(l) - newgrid_bot(k)
99         dtb = oldgrid_top(l) - newgrid_bot(k)
100         dbt = oldgrid_bot(l) - newgrid_top(k)
101         dtt = oldgrid_top(l) - newgrid_top(k)
102         if (dbt >= 0.0) exit ! Do next level in the new grid
103         if (dtb > 0.0) then
104           if (dbb <= 0.0) then
105             if (dtt <= 0) then
106               w = dtb
107             else
108               w = newgrid_top(k) - newgrid_bot(k)
109             endif
110           else
111             if (dtt <= 0) then
112               w = oldgrid_top(l) - oldgrid_bot(l)
113             else
114               w = -dbt
115             endif
116           endif
117           ! If layers overlap (w/=0), then accumulate
118           if (w /= 0.0) then
119             Nw = Nw + 1
120             wt = wt + w
121             do j=1,Ncolumns
122               if (lunits) then
123                 if (y(i,j,l) /= R_UNDEF) then
124                   yp = 10._wp**(y(i,j,l)/10._wp)
125                 else
126                   yp = 0._wp
127                 endif
128               else
129                 yp = y(i,j,l)
130               endif
131               r(i,j,k) = r(i,j,k) + w*yp
132             enddo
133           endif
134         endif
135       enddo
136       l = l - 2
137       if (l < 1) l = 0
138       ! Calculate average in new grid
139       if (Nw > 0) then
140         do j=1,Ncolumns
141           r(i,j,k) = r(i,j,k)/wt
142         enddo
143       endif
144     enddo
145   enddo
146
147   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
148   do k=1,Nglevels
149     do j=1,Ncolumns
150       do i=1,Npoints
151         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
152           if (lunits) then
153             if (r(i,j,k) <= 0.0) then
154               r(i,j,k) = R_UNDEF
155             else
156               r(i,j,k) = 10._wp*log10(r(i,j,k))
157             endif
158           endif
159         else ! Level below surface
160           r(i,j,k) = R_GROUND
161         endif
162       enddo
163     enddo
164   enddo
165
166END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
167
168  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
169  !------------- SUBROUTINE COSP_LIDAR_ONLY_CLOUD -----------------
170  ! (c) 2008, Lawrence Livermore National Security Limited Liability Corporation.
171  ! All rights reserved.
172  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
173  SUBROUTINE COSP_LIDAR_ONLY_CLOUD(Npoints, Ncolumns, Nlevels, beta_tot, beta_mol,       &
174     Ze_tot, lidar_only_freq_cloud, tcc, radar_tcc, radar_tcc2)
175    ! Inputs
176    integer,intent(in) :: &
177         Npoints,       & ! Number of horizontal gridpoints
178         Ncolumns,      & ! Number of subcolumns
179         Nlevels          ! Number of vertical layers
180    real(wp),dimension(Npoints,Nlevels),intent(in) :: &
181         beta_mol         ! Molecular backscatter
182    real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: &
183         beta_tot,      & ! Total backscattered signal
184         Ze_tot           ! Radar reflectivity
185    ! Outputs
186    real(wp),dimension(Npoints,Nlevels),intent(out) :: &
187         lidar_only_freq_cloud
188    real(wp),dimension(Npoints),intent(out) ::&
189         tcc,       & !
190         radar_tcc, & !
191         radar_tcc2   !
192   
193    ! local variables
194    real(wp) :: sc_ratio
195    real(wp),parameter :: &
196         s_cld=5.0, &
197         s_att=0.01
198    integer :: flag_sat,flag_cld,pr,i,j,flag_radarcld,flag_radarcld_no1km,j_1km
199   
200    lidar_only_freq_cloud = 0._wp
201    tcc = 0._wp
202    radar_tcc = 0._wp
203    radar_tcc2 = 0._wp
204    do pr=1,Npoints
205       do i=1,Ncolumns
206          flag_sat = 0
207          flag_cld = 0
208          flag_radarcld = 0 !+JEK
209          flag_radarcld_no1km=0 !+JEK
210          ! look for j_1km from bottom to top
211          j = 1
212          do while (Ze_tot(pr,i,j) .eq. R_GROUND)
213             j = j+1
214          enddo
215          j_1km = j+1  !this is the vertical index of 1km above surface 
216         
217          do j=1,Nlevels
218             sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j)
219             if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j
220             if (Ze_tot(pr,i,j) .lt. -30.) then  !radar can't detect cloud
221                if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then  !lidar sense cloud
222                   lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf
223                   flag_cld=1
224                endif
225             else  !radar sense cloud (z%Ze_tot(pr,i,j) .ge. -30.)
226                flag_cld=1
227                flag_radarcld=1
228                if (j .gt. j_1km) flag_radarcld_no1km=1             
229             endif
230          enddo !levels
231          if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1._wp
232          if (flag_radarcld .eq. 1) radar_tcc(pr)=radar_tcc(pr)+1.
233          if (flag_radarcld_no1km .eq. 1) radar_tcc2(pr)=radar_tcc2(pr)+1.       
234       enddo !columns
235    enddo !points
236    lidar_only_freq_cloud=lidar_only_freq_cloud/Ncolumns
237    tcc=tcc/Ncolumns
238    radar_tcc=radar_tcc/Ncolumns
239    radar_tcc2=radar_tcc2/Ncolumns
240   
241    ! Unit conversion
242    where(lidar_only_freq_cloud /= R_UNDEF) &
243            lidar_only_freq_cloud = lidar_only_freq_cloud*100._wp
244    where(tcc /= R_UNDEF) tcc = tcc*100._wp
245    where(radar_tcc /= R_UNDEF) radar_tcc = radar_tcc*100._wp
246    where(radar_tcc2 /= R_UNDEF) radar_tcc2 = radar_tcc2*100._wp
247   
248  END SUBROUTINE COSP_LIDAR_ONLY_CLOUD
249 
250  ! ######################################################################################
251  ! FUNCTION hist1D
252  ! ######################################################################################
253  function hist1d(Npoints,var,nbins,bins)
254    ! Inputs
255    integer,intent(in) :: &
256         Npoints, & ! Number of points in input array
257         Nbins      ! Number of bins for sorting
258    real(wp),intent(in),dimension(Npoints) :: &
259         var        ! Input variable to be sorted
260    real(wp),intent(in),dimension(Nbins+1) :: &
261         bins       ! Histogram bins [lowest,binTops] 
262    ! Outputs
263    real(wp),dimension(Nbins) :: &
264         hist1d     ! Output histogram     
265    ! Local variables
266    integer :: ij
267   
268    do ij=2,Nbins+1 
269       hist1D(ij-1) = count(var .ge. bins(ij-1) .and. var .lt. bins(ij))
270       if (count(var .eq. R_GROUND) .ge. 1) hist1D(ij-1)=R_UNDEF
271    enddo
272   
273  end function hist1D
274 
275  ! ######################################################################################
276  ! SUBROUTINE hist2D
277  ! ######################################################################################
278  subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
279    implicit none
280   
281    ! INPUTS
282    integer, intent(in) :: &
283         npts,  & ! Number of data points to be sorted
284         nbin1, & ! Number of bins in histogram direction 1
285         nbin2    ! Number of bins in histogram direction 2
286    real(wp),intent(in),dimension(npts) :: &
287         var1,  & ! Variable 1 to be sorted into bins
288         var2     ! variable 2 to be sorted into bins
289    real(wp),intent(in),dimension(nbin1+1) :: &
290         bin1     ! Histogram bin 1 boundaries
291    real(wp),intent(in),dimension(nbin2+1) :: &
292         bin2     ! Histogram bin 2 boundaries
293    ! OUTPUTS
294    real(wp),intent(out),dimension(nbin1,nbin2) :: &
295         jointHist
296   
297    ! LOCAL VARIABLES
298    integer :: ij,ik
299   
300    do ij=2,nbin1+1
301       do ik=2,nbin2+1
302          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
303               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))       
304       enddo
305    enddo
306  end subroutine hist2D
307END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.