source: trunk/libf/cosp/cosp_stats.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 11.3 KB
RevLine 
[1]1! (c) British Crown Copyright 2008, the Met Office.
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 Met Office nor the names of its contributors may be used
13!       to endorse or promote products derived from this software without specific prior written
14!       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
25!
26! History:
27! Jul 2007 - A. Bodas-Salcedo - Initial version
28! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid
29! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
30! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
31! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
32!
33!
34MODULE MOD_COSP_STATS
35  USE MOD_COSP_CONSTANTS
36  USE MOD_COSP_TYPES
37  USE MOD_LLNL_STATS
38  USE MOD_LMD_IPSL_STATS
39  IMPLICIT NONE
40
41CONTAINS
42
43!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44!------------------- SUBROUTINE COSP_STATS ------------------------
45!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
47
48   ! Input arguments
49   type(cosp_gridbox),intent(in) :: gbx
50   type(cosp_subgrid),intent(in) :: sgx
51   type(cosp_config),intent(in)  :: cfg
52   type(cosp_sgradar),intent(in) :: sgradar
53   type(cosp_sglidar),intent(in) :: sglidar
54   type(cosp_vgrid),intent(in)   :: vgrid
55   ! Output arguments
56   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
57   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
58
59   ! Local variables
60   integer :: Npoints  !# of grid points
61   integer :: Nlevels  !# of levels
62   integer :: Nhydro   !# of hydrometeors
63   integer :: Ncolumns !# of columns
64   integer :: Nlr
65   logical :: ok_lidar_cfad = .false.
66   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
67   real,dimension(:,:),allocatable :: ph_c,betamol_c
68
69   Npoints  = gbx%Npoints
70   Nlevels  = gbx%Nlevels
71   Nhydro   = gbx%Nhydro
72   Ncolumns = gbx%Ncolumns
73   Nlr      = vgrid%Nlvgrid
74
75   if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true.
76
77   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
78        allocate(Ze_out(Npoints,Ncolumns,Nlr),betatot_out(Npoints,Ncolumns,Nlr), &
79                 betamol_in(Npoints,1,Nlevels),betamol_out(Npoints,1,Nlr),betamol_c(Npoints,Nlr), &
80                 ph_in(Npoints,1,Nlevels),ph_out(Npoints,1,Nlr),ph_c(Npoints,Nlr))
81        Ze_out = 0.0
82        betatot_out  = 0.0
83        betamol_out= 0.0
84        betamol_c  = 0.0
85        ph_in(:,1,:)  = gbx%ph(:,:)
86        ph_out  = 0.0
87        ph_c    = 0.0
88        !++++++++++++ Radar CFAD ++++++++++++++++
89        if (cfg%Lradar_sim) then
90            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
91                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
92            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
93                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
94        endif
95        !++++++++++++ Lidar CFAD ++++++++++++++++
96        if (cfg%Llidar_sim) then
97            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
98            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
99                                           Nlr,vgrid%zl,vgrid%zu,betamol_out)
100            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%beta_tot, &
101                                           Nlr,vgrid%zl,vgrid%zu,betatot_out)
102            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,ph_in, &
103                                           Nlr,vgrid%zl,vgrid%zu,ph_out)
104            ph_c(:,:) = ph_out(:,1,:)
105            betamol_c(:,:) = betamol_out(:,1,:)
106            ! Stats from lidar_stat_summary
107            call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
108                            ,betatot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
109                            ,LIDAR_UNDEF,ok_lidar_cfad &
110                            ,stlidar%cfad_sr,stlidar%srbval &
111                            ,LIDAR_NCAT,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
112        endif
113        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
114        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
115                                    betatot_out,betamol_c,Ze_out, &
116                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
117        ! Deallocate arrays at coarse resolution
118        deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
119   else ! Statistics in model levels
120        !++++++++++++ Radar CFAD ++++++++++++++++
121        if (cfg%Lradar_sim) stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,sgradar%Ze_tot, &
122                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
123        !++++++++++++ Lidar CFAD ++++++++++++++++
124        ! Stats from lidar_stat_summary
125        if (cfg%Llidar_sim) call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
126                        ,sglidar%beta_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
127                        ,LIDAR_UNDEF,ok_lidar_cfad &
128                        ,stlidar%cfad_sr,stlidar%srbval &
129                        ,LIDAR_NCAT,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
130        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
131        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
132                                    sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
133                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
134   endif
135   ! Replace undef
136   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
137   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
138   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
139   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
140
141END SUBROUTINE COSP_STATS
142
143
144!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
146!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,M,zl,zu,r,log_units)
148   implicit none
149   ! Input arguments
150   integer,intent(in) :: Npoints  !# of grid points
151   integer,intent(in) :: Nlevels  !# of levels
152   integer,intent(in) :: Ncolumns !# of columns
153   real,dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
154   real,dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
155   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
156   integer,intent(in) :: M  !# levels in the new grid
157   real,dimension(M),intent(in) :: zl ! Lower boundary of new levels  [m]
158   real,dimension(M),intent(in) :: zu ! Upper boundary of new levels  [m]
159   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
160   ! Output
161   real,dimension(Npoints,Ncolumns,M),intent(out) :: r ! Variable on new grid
162
163   ! Local variables
164   integer :: i,j,k
165   logical :: lunits
166
167   integer :: l
168   real,dimension(Npoints) :: ws,sumwyp
169   real,dimension(Npoints,Nlevels) :: xl,xu
170   real,dimension(Npoints,Nlevels) :: w
171   real,dimension(Npoints,Ncolumns,Nlevels) :: yp
172
173   lunits=.false.
174   if (present(log_units)) lunits=log_units
175
176   r(:,:,:) = R_GROUND
177   ! Vertical grid at that point
178   xl(:,:) = zhalf(:,:)
179   xu(:,1:Nlevels-1) = xl(:,2:Nlevels)
180   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
181   yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid
182   ! Check for dBZ and change if necessary
183   if (lunits) then
184     where (y /= R_UNDEF)
185       yp = 10.0**(y/10.0)
186     elsewhere
187       yp = 0.0
188     end where
189   endif
190   do k=1,M
191     ! Find weights
192     w(:,:) = 0.0
193     do j=1,Nlevels
194       do i=1,Npoints
195         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
196           !xl(j)-----------------xu(j)
197           !      zl(k)------------------------------zu(k)
198           w(i,j) = xu(i,j) - zl(k)
199         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
200           !           xl(j)-----------------xu(j)
201           !      zl(k)------------------------------zu(k)
202           w(i,j) = xu(i,j) - xl(i,j)
203         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
204           !                           xl(j)-----------------xu(j)
205           !      zl(k)------------------------------zu(k)
206           w(i,j) = zu(k) - xl(i,j)
207         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
208           !  xl(j)---------------------------xu(j)
209           !        zl(k)--------------zu(k)
210           w(i,j) = zu(k) - zl(k)
211         endif
212       enddo
213     enddo
214     ! Do the weighted mean
215     do j=1,Ncolumns
216       ws    (:) = 0.0
217       sumwyp(:) = 0.0
218       do l=1,Nlevels
219         do i=1,Npoints
220           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
221             ws    (i) = ws    (i) + w(i,l)
222             sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l)
223           endif
224         enddo
225       enddo
226       do i=1,Npoints
227         if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
228           if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i)
229         endif
230       enddo
231     enddo
232   enddo
233   ! Check for dBZ and change if necessary
234   if (lunits) then
235     do k=1,M
236       do j=1,Ncolumns
237         do i=1,Npoints
238           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
239             if (r(i,j,k) <= 0.0) then
240                 r(i,j,k) = R_UNDEF
241             else
242                 r(i,j,k) = 10.0*log10(r(i,j,k))
243             endif
244           endif
245         enddo
246       enddo
247     enddo
248   endif
249
250
251
252END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
253
254END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.