source: LMDZ6/trunk/libf/phylmd/cosp/mod_cosp_stats.F90 @ 5300

Last change on this file since 5300 was 5271, checked in by abarral, 5 weeks ago

Move dimensions.h into a module
Nb: doesn't compile yet

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 14.3 KB
Line 
1! (c) British Crown Copyright 2008, the Met Office.
2! All rights reserved.
3! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
4! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_stats.F90 $
5!
6! Redistribution and use in source and binary forms, with or without modification, are permitted
7! provided that the following conditions are met:
8!
9!     * Redistributions of source code must retain the above copyright notice, this list
10!       of conditions and the following disclaimer.
11!     * 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 materials
13!       provided with the distribution.
14!     * Neither the name of the Met Office nor the names of its contributors may be used
15!       to endorse or promote products derived from this software without specific prior written
16!       permission.
17!
18! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
19! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
20! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
21! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
24! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
25! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26
27!
28! History:
29! Jul 2007 - A. Bodas-Salcedo - Initial version
30! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid
31! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
32! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
33! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
34! Jan 2013 - G. Cesana        - Added betaperp and temperature arguments
35!                             - Added phase 3D/3Dtemperature/Map output variables in diag_lidar
36!
37!
38INCLUDE "cosp_defs.h"
39MODULE MOD_COSP_STATS
40  USE MOD_COSP_CONSTANTS
41  USE MOD_COSP_TYPES
42  USE MOD_LLNL_STATS
43  USE MOD_LMD_IPSL_STATS
44  IMPLICIT NONE
45
46CONTAINS
47
48!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49!------------------- SUBROUTINE COSP_STATS ------------------------
50!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
52
53   ! Input arguments
54   type(cosp_gridbox),intent(in) :: gbx
55   type(cosp_subgrid),intent(in) :: sgx
56   type(cosp_config),intent(in)  :: cfg
57   type(cosp_sgradar),intent(in) :: sgradar
58   type(cosp_sglidar),intent(in) :: sglidar
59   type(cosp_vgrid),intent(in)   :: vgrid
60   ! Output arguments
61   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
62   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
63
64   ! Local variables
65   integer :: Npoints  !# of grid points
66   integer :: Nlevels  !# of levels
67   integer :: Nhydro   !# of hydrometeors
68   integer :: Ncolumns !# of columns
69   integer :: Nlr
70   logical :: ok_lidar_cfad = .false.
71   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
72   real,dimension(:,:),allocatable :: ph_c,betamol_c
73   real,dimension(:,:,:),allocatable ::  betaperptot_out, temp_in, temp_out
74   real,dimension(:,:),allocatable :: temp_c
75
76   Npoints  = gbx%Npoints
77   Nlevels  = gbx%Nlevels
78   Nhydro   = gbx%Nhydro
79   Ncolumns = gbx%Ncolumns
80   Nlr      = vgrid%Nlvgrid
81
82   if (cfg%LcfadLidarsr532) ok_lidar_cfad=.true.
83
84   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
85        allocate(Ze_out(Npoints,Ncolumns,Nlr),betatot_out(Npoints,Ncolumns,Nlr), &
86                 betamol_in(Npoints,1,Nlevels),betamol_out(Npoints,1,Nlr),betamol_c(Npoints,Nlr), &
87                 ph_in(Npoints,1,Nlevels),ph_out(Npoints,1,Nlr),ph_c(Npoints,Nlr))
88        Ze_out = 0.0
89        betatot_out  = 0.0
90        betamol_out= 0.0
91        betamol_c  = 0.0
92        ph_in(:,1,:)  = gbx%ph(:,:)
93        ph_out  = 0.0
94        ph_c    = 0.0
95        allocate(betaperptot_out(Npoints,Ncolumns,Nlr),temp_in(Npoints,1,Nlevels),temp_out(Npoints,1,Nlr), &
96                 temp_c(Npoints,Nlr))
97        betaperptot_out = 0.0
98        temp_in = 0.0
99        temp_out = 0.0
100        temp_c = 0.0
101
102        !++++++++++++ Radar CFAD ++++++++++++++++
103        if (cfg%Lradar_sim) then
104            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
105                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
106            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
107                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
108        endif
109        !++++++++++++ Lidar CFAD ++++++++++++++++
110        if (cfg%Llidar_sim) then
111            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
112            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
113                                           Nlr,vgrid%zl,vgrid%zu,betamol_out)
114            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%beta_tot, &
115                                           Nlr,vgrid%zl,vgrid%zu,betatot_out)
116
117            temp_in(:,1,:) = gbx%T(:,:)
118            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%betaperp_tot, &
119                                           Nlr,vgrid%zl,vgrid%zu,betaperptot_out)
120            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,temp_in, &
121                                           Nlr,vgrid%zl,vgrid%zu,temp_out)
122            temp_c(:,:) = temp_out(:,1,:)
123            stlidar%proftemp = temp_c                                     !TIBO
124            where (stlidar%proftemp  < 150.) stlidar%proftemp   = R_UNDEF !TIBO
125            where (stlidar%proftemp  > 350.) stlidar%proftemp   = R_UNDEF !TIBO
126
127            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,ph_in, &
128                                           Nlr,vgrid%zl,vgrid%zu,ph_out)
129            ph_c(:,:) = ph_out(:,1,:)
130            betamol_c(:,:) = betamol_out(:,1,:)
131            ! Stats from lidar_stat_summary
132            call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
133                            ,temp_c,betatot_out,betaperptot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
134                            ,LIDAR_UNDEF,ok_lidar_cfad &
135                            ,stlidar%cfad_sr,stlidar%srbval &
136                            ,LIDAR_NCAT,LIDAR_NTYPE,stlidar%lidarcld,stlidar%lidarcldtype & !OPAQ
137                            ,stlidar%lidarcldphase,stlidar%cldlayer,stlidar%cldtype &       !OPAQ
138                            ,stlidar%cldlayerphase,stlidar%lidarcldtmp &                    !OPAQ
139                            ,stlidar%parasolrefl,vgrid%z,stlidar%profSR)                    !OPAQ !TIBO
140        endif
141
142        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
143        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
144                                    temp_c,betatot_out,betaperptot_out,betamol_c,Ze_out, &
145                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
146        deallocate(temp_in,temp_out,temp_c,betaperptot_out) !TIBO +temp_in
147
148        ! Deallocate arrays at coarse resolution
149        deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
150   else ! Statistics in model levels
151        !++++++++++++ Radar CFAD ++++++++++++++++
152        if (cfg%Lradar_sim) stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,sgradar%Ze_tot, &
153                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
154        !++++++++++++ Lidar CFAD ++++++++++++++++
155        ! Stats from lidar_stat_summary
156        if (cfg%Llidar_sim) call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
157                        ,sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
158                        ,LIDAR_UNDEF,ok_lidar_cfad &
159                        ,stlidar%cfad_sr,stlidar%srbval &
160                        ,LIDAR_NCAT,LIDAR_NTYPE,stlidar%lidarcld,stlidar%lidarcldtype & !OPAQ
161                        ,stlidar%lidarcldphase,stlidar%cldlayer,stlidar%cldtype &       !OPAQ
162                        ,stlidar%cldlayerphase,stlidar%lidarcldtmp &                    !OPAQ
163                        ,stlidar%parasolrefl,vgrid%z,stlidar%profSR)                    !OPAQ !TIBO
164        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
165        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
166                                    sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sgradar%Ze_tot, &
167                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
168   endif
169   ! Replace undef
170   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
171   where (stlidar%profSR   == LIDAR_UNDEF) stlidar%profSR   = R_UNDEF !TIBO
172   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
173   where (stlidar%lidarcldtype  == LIDAR_UNDEF) stlidar%lidarcldtype  = R_UNDEF !OPAQ
174   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
175   where (stlidar%cldtype  == LIDAR_UNDEF) stlidar%cldtype  = R_UNDEF           !OPAQ
176   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
177   where (stlidar%cldlayerphase  == LIDAR_UNDEF) stlidar%cldlayerphase  = R_UNDEF
178   where (stlidar%lidarcldphase  == LIDAR_UNDEF) stlidar%lidarcldphase  = R_UNDEF
179   where (stlidar%lidarcldtmp  == LIDAR_UNDEF) stlidar%lidarcldtmp  = R_UNDEF
180
181END SUBROUTINE COSP_STATS
182
183!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
185!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
187   implicit none
188   ! Input arguments
189   integer,intent(in) :: Npoints  !# of grid points
190   integer,intent(in) :: Nlevels  !# of levels
191   integer,intent(in) :: Ncolumns !# of columns
192   real,dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
193   real,dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
194   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
195   integer,intent(in) :: Nglevels  !# levels in the new grid
196   real,dimension(Nglevels),intent(in) :: newgrid_bot ! Lower boundary of new levels  [m]
197   real,dimension(Nglevels),intent(in) :: newgrid_top ! Upper boundary of new levels  [m]
198   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
199   ! Output
200   real,dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
201
202   ! Local variables
203   integer :: i,j,k
204   logical :: lunits
205   integer :: l
206   real :: w ! Weight
207   real :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
208   integer :: Nw  ! Number of weights
209   real :: wt  ! Sum of weights
210   real,dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
211   real :: yp ! Local copy of y at a particular point.
212              ! This allows for change of units.
213
214   lunits=.false.
215   if (present(log_units)) lunits=log_units
216
217   r = 0.0
218
219   do i=1,Npoints
220     ! Calculate tops and bottoms of new and old grids
221     oldgrid_bot = zhalf(i,:)
222     oldgrid_top(1:Nlevels-1) = oldgrid_bot(2:Nlevels)
223     oldgrid_top(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
224     l = 0 ! Index of level in the old grid
225     ! Loop over levels in the new grid
226     do k = 1,Nglevels
227       Nw = 0 ! Number of weigths
228       wt = 0.0 ! Sum of weights
229       ! Loop over levels in the old grid and accumulate total for weighted average
230       do
231         l = l + 1
232         w = 0.0 ! Initialise weight to 0
233         ! Distances between edges of both grids
234         dbb = oldgrid_bot(l) - newgrid_bot(k)
235         dtb = oldgrid_top(l) - newgrid_bot(k)
236         dbt = oldgrid_bot(l) - newgrid_top(k)
237         dtt = oldgrid_top(l) - newgrid_top(k)
238         if (dbt >= 0.0) exit ! Do next level in the new grid
239         if (dtb > 0.0) then
240           if (dbb <= 0.0) then
241             if (dtt <= 0) then
242               w = dtb
243             else
244               w = newgrid_top(k) - newgrid_bot(k)
245             endif
246           else
247             if (dtt <= 0) then
248               w = oldgrid_top(l) - oldgrid_bot(l)
249             else
250               w = -dbt
251             endif
252           endif
253           ! If layers overlap (w/=0), then accumulate
254           if (w /= 0.0) then
255             Nw = Nw + 1
256             wt = wt + w
257             do j=1,Ncolumns
258               if (lunits) then
259                 if (y(i,j,l) /= R_UNDEF) then
260                   yp = 10.0**(y(i,j,l)/10.0)
261                 else
262                   yp = 0.0
263                 endif
264               else
265                 yp = y(i,j,l)
266               endif
267               r(i,j,k) = r(i,j,k) + w*yp
268             enddo
269           endif
270         endif
271       enddo
272       l = l - 2
273       if (l < 1) l = 0
274       ! Calculate average in new grid
275       if (Nw > 0) then
276         do j=1,Ncolumns
277           r(i,j,k) = r(i,j,k)/wt
278         enddo
279       endif
280     enddo
281   enddo
282
283   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
284   do k=1,Nglevels
285     do j=1,Ncolumns
286       do i=1,Npoints
287         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
288           if (lunits) then
289             if (r(i,j,k) <= 0.0) then
290               r(i,j,k) = R_UNDEF
291             else
292               r(i,j,k) = 10.0*log10(r(i,j,k))
293             endif
294           endif
295         else ! Level below surface
296           r(i,j,k) = R_GROUND
297         endif
298       enddo
299     enddo
300   enddo
301
302END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
303
304END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.