source: LMDZ4/branches/LMDZ4-dev-20091210/libf/cosp/cosp_stats.F90 @ 3655

Last change on this file since 3655 was 1262, checked in by idelkadi, 15 years ago

-Rajout des fichiers cosp_input_nl.txt et cosp_output_nl.txt pour controler les parametres entrees sorties COSP
-Rajout du repertoire cosp contenant les sources de COSP

File size: 10.9 KB
Line 
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!
32!
33MODULE MOD_COSP_STATS
34  USE MOD_COSP_CONSTANTS
35  USE MOD_COSP_TYPES
36  USE MOD_LLNL_STATS
37  USE MOD_LMD_IPSL_STATS
38  IMPLICIT NONE
39
40CONTAINS
41
42!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43!------------------- SUBROUTINE COSP_STATS ------------------------
44!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
46 
47   ! Input arguments
48   type(cosp_gridbox),intent(in) :: gbx
49   type(cosp_subgrid),intent(in) :: sgx
50   type(cosp_config),intent(in)  :: cfg
51   type(cosp_sgradar),intent(in) :: sgradar
52   type(cosp_sglidar),intent(in) :: sglidar
53   type(cosp_vgrid),intent(in)   :: vgrid
54   ! Output arguments
55   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
56   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
57   
58   ! Local variables
59   integer :: Npoints  !# of grid points
60   integer :: Nlevels  !# of levels
61   integer :: Nhydro   !# of hydrometeors
62   integer :: Ncolumns !# of columns
63   integer :: Nlr
64   logical :: ok_lidar_cfad = .false.
65   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
66   real,dimension(:,:),allocatable :: ph_c,betamol_c
67 
68   Npoints  = gbx%Npoints
69   Nlevels  = gbx%Nlevels
70   Nhydro   = gbx%Nhydro
71   Ncolumns = gbx%Ncolumns
72   Nlr      = vgrid%Nlvgrid
73 
74   if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
75
76   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
77        allocate(Ze_out(Npoints,Ncolumns,Nlr),betatot_out(Npoints,Ncolumns,Nlr), &
78                 betamol_in(Npoints,1,Nlevels),betamol_out(Npoints,1,Nlr),betamol_c(Npoints,Nlr), &
79                 ph_in(Npoints,1,Nlevels),ph_out(Npoints,1,Nlr),ph_c(Npoints,Nlr))
80        Ze_out = 0.0
81        betatot_out  = 0.0
82        betamol_in(:,1,:) = sglidar%beta_mol(:,:)
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            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
98                                           Nlr,vgrid%zl,vgrid%zu,betamol_out)
99            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%beta_tot, &
100                                           Nlr,vgrid%zl,vgrid%zu,betatot_out)
101            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,ph_in, &
102                                           Nlr,vgrid%zl,vgrid%zu,ph_out)
103            ph_c(:,:) = ph_out(:,1,:)
104            betamol_c(:,:) = betamol_out(:,1,:)
105            ! Stats from lidar_stat_summary
106            call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
107                            ,betatot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
108                            ,LIDAR_UNDEF,ok_lidar_cfad &
109                            ,stlidar%cfad_sr,stlidar%srbval &
110                            ,LIDAR_NCAT,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
111        endif
112        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
113        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
114                                    betatot_out,betamol_c,Ze_out, &
115                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
116        ! Deallocate arrays at coarse resolution
117        deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
118   else ! Statistics in model levels
119        !++++++++++++ Radar CFAD ++++++++++++++++
120        if (cfg%Lradar_sim) stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,sgradar%Ze_tot, &
121                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
122        !++++++++++++ Lidar CFAD ++++++++++++++++
123        ! Stats from lidar_stat_summary
124        if (cfg%Llidar_sim) call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
125                        ,sglidar%beta_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
126                        ,LIDAR_UNDEF,ok_lidar_cfad &
127                        ,stlidar%cfad_sr,stlidar%srbval &
128                        ,LIDAR_NCAT,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
129        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
130        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
131                                    sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
132                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
133   endif
134   ! Replace undef
135   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
136   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
137   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
138   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
139
140END SUBROUTINE COSP_STATS
141
142
143!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
145!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
146SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,M,zl,zu,r,log_units)
147   implicit none
148   ! Input arguments
149   integer,intent(in) :: Npoints  !# of grid points
150   integer,intent(in) :: Nlevels  !# of levels
151   integer,intent(in) :: Ncolumns !# of columns
152   real,dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
153   real,dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
154   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
155   integer,intent(in) :: M  !# levels in the new grid
156   real,dimension(M),intent(in) :: zl ! Lower boundary of new levels  [m]
157   real,dimension(M),intent(in) :: zu ! Upper boundary of new levels  [m]
158   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
159   ! Output
160   real,dimension(Npoints,Ncolumns,M),intent(out) :: r ! Variable on new grid
161
162   ! Local variables
163   integer :: i,j,k
164   logical :: lunits
165   real :: ws
166   real,dimension(Nlevels) :: xl,xu ! Lower and upper boundaries of model grid
167   real,dimension(M) :: dz          ! Layer depth
168   real,dimension(Nlevels,M) :: w   ! Weights to do the mean at each point
169   real,dimension(Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
170                                           ! Local copy at a particular point.
171                                           ! This allows for change of units.
172   
173   lunits=.false.
174   if (present(log_units)) lunits=log_units
175   
176   r = 0.0
177   do i=1,Npoints
178     ! Vertical grid at that point
179     xl = zhalf(i,:)
180     xu(1:Nlevels-1) = xl(2:Nlevels)
181     xu(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
182     dz = zu - zl
183     yp = y(i,:,:) ! Temporary variable to regrid
184     ! Find weights
185     w = 0.0
186     do k=1,M
187       do j=1,Nlevels
188         if ((xl(j) < zl(k)).and.(xu(j) > zl(k)).and.(xu(j) <= zu(k))) then
189           !xl(j)-----------------xu(j)
190           !      zl(k)------------------------------zu(k)
191           w(j,k) = xu(j) - zl(k)
192         else if ((xl(j) >= zl(k)).and.(xu(j) <= zu(k))) then
193           !           xl(j)-----------------xu(j)
194           !      zl(k)------------------------------zu(k)
195           w(j,k) = xu(j) - xl(j)
196         else if ((xl(j) >= zl(k)).and.(xl(j) < zu(k)).and.(xu(j) >= zu(k))) then
197           !                           xl(j)-----------------xu(j)
198           !      zl(k)------------------------------zu(k)
199           w(j,k) = zu(k) - xl(j)
200         else if ((xl(j) <= zl(k)).and.(xu(j) >= zu(k))) then
201           !  xl(j)---------------------------xu(j)
202           !        zl(k)--------------zu(k)
203           w(j,k) = dz(j)
204         endif
205       enddo
206     enddo
207     ! Check for dBZ and change if necessary
208     if (lunits) then
209        where (yp /= R_UNDEF)
210          yp = 10.0**(yp/10.0)
211        elsewhere
212          yp = 0.0
213        end where
214     endif
215     ! Do the weighted mean
216     do j=1,Ncolumns
217       do k=1,M
218          ws = sum(w(:,k))
219          if (ws > 0.0) r(i,j,k) = sum(w(:,k)*yp(j,:))/ws
220       enddo
221     enddo
222     ! Check for dBZ and change if necessary
223     if (lunits) then
224        where (r(i,:,:) <= 0.0)
225          r(i,:,:) = R_UNDEF
226        elsewhere
227          r(i,:,:) = 10.0*log10(r(i,:,:))
228        end where
229     endif
230   enddo
231 
232 
233   
234END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
235
236END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.