source: LMDZ4/branches/LMDZ4V5.0-dev/libf/cosp/cosp_stats.F90 @ 1401

Last change on this file since 1401 was 1395, checked in by Laurent Fairhead, 14 years ago

Optimisation of routine for NEC SX8. Routine gives the same results as before.


Optimisation de la routine pour la NEC SX8. La nouvelle routine donne les
mêmes résultats que la précédente

File size: 11.3 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! May 2010 - L. Fairhead      - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer
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_in(:,1,:) = sglidar%beta_mol(:,:)
84        betamol_out= 0.0
85        betamol_c  = 0.0
86        ph_in(:,1,:)  = gbx%ph(:,:)
87        ph_out  = 0.0
88        ph_c    = 0.0
89        !++++++++++++ Radar CFAD ++++++++++++++++
90        if (cfg%Lradar_sim) then
91           call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
92                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
93            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
94                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
95        endif
96        !++++++++++++ Lidar CFAD ++++++++++++++++
97        if (cfg%Llidar_sim) then
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,c
165   logical :: lunits
166   real :: ws(Npoints,M), ws_temp(Npoints,M)
167   real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid
168   real,dimension(M) :: dz          ! Layer depth
169   real,dimension(Npoints,Nlevels,M) :: w   ! Weights to do the mean at each point
170   real,dimension(Npoints,Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
171                                           ! Local copy at a particular point.
172                                           ! This allows for change of units.
173 
174
175   lunits=.false.
176   if (present(log_units)) lunits=log_units
177
178   r(:,:,:) = 0.0
179   yp(:,:,:) = y(:,:,:)
180   w(:,:,:) = 0.0
181   ws(:,:) = 0.0
182   ws_temp(:,:) = 0.0
183   dz = zu - zl
184
185   ! Check for dBZ and change if necessary
186   if (lunits) then
187     where (y /= R_UNDEF)
188        yp = 10.0**(y/10.0)
189     end where
190   endif
191
192   ! Vertical grids
193   xl(:,:) = zhalf(:,:)
194   xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels)
195   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
196  ! Find weights
197   do k=1, M
198     do j=1, Nlevels
199       do i=1, Npoints
200         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
201           !xl(j)-----------------xu(j)
202           !      zl(k)------------------------------zu(k)
203           w(i,j,k) = xu(i,j) - zl(k)
204         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
205           !           xl(j)-----------------xu(j)
206           !      zl(k)------------------------------zu(k)
207           w(i,j,k) =  xl(i,j+1) - xl(i,j)
208         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
209           !                           xl(j)-----------------xu(j)
210           !      zl(k)------------------------------zu(k)
211           w(i,j,k) = zu(k) - xl(i,j)
212         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
213           !  xl(j)---------------------------xu(j)
214           !        zl(k)--------------zu(k)
215           w(i,j,k) = dz(j)
216         endif
217       enddo
218     enddo
219   enddo
220
221   ! Do the weighted mean
222   do k=1,M
223     do i = 1, Npoints
224        ws(i,k) = sum(w(i,:,k))
225     enddo
226   enddo
227
228   ws_temp = 1.
229   where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:)
230
231   do c=1,Ncolumns
232     do k=1,M
233       do i = 1, Npoints
234          r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k)
235       enddo
236     enddo
237   enddo
238
239   do k=1,M
240     do i = 1, Npoints
241        if (ws(i,k) <= 0.0) r(i,:,k)=0.0
242     enddo
243   enddo
244 
245   ! Check for dBZ and change if necessary
246   if (lunits) then
247      where (r(:,:,:) <= 0.0)
248        r(:,:,:) = R_UNDEF
249      elsewhere
250        r(:,:,:) = 10.0*log10(r(:,:,:))
251      end where
252   endif
253   
254END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
255
256END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.