source: LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_stats.F90

Last change on this file was 5185, checked in by abarral, 9 days ago

Replace REPROBUS CPP KEY by logical using handmade wonky wrapper

  • 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! History:
28! Jul 2007 - A. Bodas-Salcedo - Initial version
29! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid
30! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
31! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
32! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
33! Jan 2013 - G. Cesana        - Added betaperp and temperature arguments
34!                             - Added phase 3D/3Dtemperature/Map output variables in diag_lidar
35
36
37#include "cosp_defs.h"
38MODULE MOD_COSP_STATS
39  USE MOD_COSP_CONSTANTS
40  USE MOD_COSP_TYPES
41  USE MOD_LLNL_STATS
42  USE MOD_LMD_IPSL_STATS
43  IMPLICIT NONE
44
45CONTAINS
46
47!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48!------------------- SUBROUTINE COSP_STATS ------------------------
49!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
50SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
51
52   ! Input arguments
53   type(cosp_gridbox),intent(in) :: gbx
54   type(cosp_subgrid),intent(in) :: sgx
55   type(cosp_config),intent(in)  :: cfg
56   type(cosp_sgradar),intent(in) :: sgradar
57   type(cosp_sglidar),intent(in) :: sglidar
58   type(cosp_vgrid),intent(in)   :: vgrid
59   ! Output arguments
60   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
61   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
62
63   ! Local variables
64   integer :: Npoints  !# of grid points
65   integer :: Nlevels  !# of levels
66   integer :: Nhydro   !# of hydrometeors
67   integer :: Ncolumns !# of columns
68   integer :: Nlr
69   logical :: ok_lidar_cfad = .false.
70   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
71   real,dimension(:,:),allocatable :: ph_c,betamol_c
72   real,dimension(:,:,:),allocatable ::  betaperptot_out, temp_in, temp_out
73   real,dimension(:,:),allocatable :: temp_c
74
75   Npoints  = gbx%Npoints
76   Nlevels  = gbx%Nlevels
77   Nhydro   = gbx%Nhydro
78   Ncolumns = gbx%Ncolumns
79   Nlr      = vgrid%Nlvgrid
80
81   if (cfg%LcfadLidarsr532) ok_lidar_cfad=.true.
82
83   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
84        allocate(Ze_out(Npoints,Ncolumns,Nlr),betatot_out(Npoints,Ncolumns,Nlr), &
85                 betamol_in(Npoints,1,Nlevels),betamol_out(Npoints,1,Nlr),betamol_c(Npoints,Nlr), &
86                 ph_in(Npoints,1,Nlevels),ph_out(Npoints,1,Nlr),ph_c(Npoints,Nlr))
87        Ze_out = 0.0
88        betatot_out  = 0.0
89        betamol_out= 0.0
90        betamol_c  = 0.0
91        ph_in(:,1,:)  = gbx%ph(:,:)
92        ph_out  = 0.0
93        ph_c    = 0.0
94        allocate(betaperptot_out(Npoints,Ncolumns,Nlr),temp_in(Npoints,1,Nlevels),temp_out(Npoints,1,Nlr), &
95                 temp_c(Npoints,Nlr))
96        betaperptot_out = 0.0
97        temp_in = 0.0
98        temp_out = 0.0
99        temp_c = 0.0
100
101        !++++++++++++ Radar CFAD ++++++++++++++++
102        if (cfg%Lradar_sim) then
103            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
104                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
105            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
106                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
107        endif
108        !++++++++++++ Lidar CFAD ++++++++++++++++
109        if (cfg%Llidar_sim) then
110            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
111            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
112                                           Nlr,vgrid%zl,vgrid%zu,betamol_out)
113            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%beta_tot, &
114                                           Nlr,vgrid%zl,vgrid%zu,betatot_out)
115
116            temp_in(:,1,:) = gbx%T(:,:)
117            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sglidar%betaperp_tot, &
118                                           Nlr,vgrid%zl,vgrid%zu,betaperptot_out)
119            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,temp_in, &
120                                           Nlr,vgrid%zl,vgrid%zu,temp_out)
121            temp_c(:,:) = temp_out(:,1,:)
122            stlidar%proftemp = temp_c                                     !TIBO
123            where (stlidar%proftemp  < 150.) stlidar%proftemp   = R_UNDEF !TIBO
124            where (stlidar%proftemp  > 350.) stlidar%proftemp   = R_UNDEF !TIBO
125
126            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,ph_in, &
127                                           Nlr,vgrid%zl,vgrid%zu,ph_out)
128            ph_c(:,:) = ph_out(:,1,:)
129            betamol_c(:,:) = betamol_out(:,1,:)
130            ! Stats from lidar_stat_summary
131            call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
132                            ,temp_c,betatot_out,betaperptot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
133                            ,LIDAR_UNDEF,ok_lidar_cfad &
134                            ,stlidar%cfad_sr,stlidar%srbval &
135                            ,LIDAR_NCAT,LIDAR_NTYPE,stlidar%lidarcld,stlidar%lidarcldtype & !OPAQ
136                            ,stlidar%lidarcldphase,stlidar%cldlayer,stlidar%cldtype &       !OPAQ
137                            ,stlidar%cldlayerphase,stlidar%lidarcldtmp &                    !OPAQ
138                            ,stlidar%parasolrefl,vgrid%z,stlidar%profSR)                    !OPAQ !TIBO
139        endif
140
141        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
142        if (cfg%Lradar_sim.AND.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
143                                    temp_c,betatot_out,betaperptot_out,betamol_c,Ze_out, &
144                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
145        deallocate(temp_in,temp_out,temp_c,betaperptot_out) !TIBO +temp_in
146
147        ! Deallocate arrays at coarse resolution
148        deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
149   else ! Statistics in model levels
150        !++++++++++++ Radar CFAD ++++++++++++++++
151        if (cfg%Lradar_sim) stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,sgradar%Ze_tot, &
152                                        DBZE_MIN,DBZE_MAX,CFAD_ZE_MIN,CFAD_ZE_WIDTH)
153        !++++++++++++ Lidar CFAD ++++++++++++++++
154        ! Stats from lidar_stat_summary
155        if (cfg%Llidar_sim) call diag_lidar(Npoints,Ncolumns,Nlr,SR_BINS,PARASOL_NREFL &
156                        ,sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
157                        ,LIDAR_UNDEF,ok_lidar_cfad &
158                        ,stlidar%cfad_sr,stlidar%srbval &
159                        ,LIDAR_NCAT,LIDAR_NTYPE,stlidar%lidarcld,stlidar%lidarcldtype & !OPAQ
160                        ,stlidar%lidarcldphase,stlidar%cldlayer,stlidar%cldtype &       !OPAQ
161                        ,stlidar%cldlayerphase,stlidar%lidarcldtmp &                    !OPAQ
162                        ,stlidar%parasolrefl,vgrid%z,stlidar%profSR)                    !OPAQ !TIBO
163        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
164        if (cfg%Lradar_sim.AND.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
165                                    sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sgradar%Ze_tot, &
166                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
167   endif
168   ! Replace undef
169   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
170   where (stlidar%profSR   == LIDAR_UNDEF) stlidar%profSR   = R_UNDEF !TIBO
171   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
172   where (stlidar%lidarcldtype  == LIDAR_UNDEF) stlidar%lidarcldtype  = R_UNDEF !OPAQ
173   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
174   where (stlidar%cldtype  == LIDAR_UNDEF) stlidar%cldtype  = R_UNDEF           !OPAQ
175   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
176   where (stlidar%cldlayerphase  == LIDAR_UNDEF) stlidar%cldlayerphase  = R_UNDEF
177   where (stlidar%lidarcldphase  == LIDAR_UNDEF) stlidar%lidarcldphase  = R_UNDEF
178   where (stlidar%lidarcldtmp  == LIDAR_UNDEF) stlidar%lidarcldtmp  = R_UNDEF
179
180END SUBROUTINE COSP_STATS
181
182!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
183!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
184!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
186   implicit none
187   ! Input arguments
188   integer,intent(in) :: Npoints  !# of grid points
189   integer,intent(in) :: Nlevels  !# of levels
190   integer,intent(in) :: Ncolumns !# of columns
191   real,dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
192   real,dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
193   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
194   integer,intent(in) :: Nglevels  !# levels in the new grid
195   real,dimension(Nglevels),intent(in) :: newgrid_bot ! Lower boundary of new levels  [m]
196   real,dimension(Nglevels),intent(in) :: newgrid_top ! Upper boundary of new levels  [m]
197   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
198   ! Output
199   real,dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
200
201   ! Local variables
202   integer :: i,j,k
203   logical :: lunits
204   integer :: l
205   REAL :: w ! Weight
206   REAL :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
207   integer :: Nw  ! Number of weights
208   REAL :: wt  ! Sum of weights
209   real,dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
210   REAL :: yp ! Local copy of y at a particular point.
211              ! This allows for change of units.
212
213   lunits=.false.
214   if (present(log_units)) lunits=log_units
215
216   r = 0.0
217
218   DO i=1,Npoints
219     ! Calculate tops and bottoms of new and old grids
220     oldgrid_bot = zhalf(i,:)
221     oldgrid_top(1:Nlevels-1) = oldgrid_bot(2:Nlevels)
222     oldgrid_top(Nlevels) = zfull(i,Nlevels) +  zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric
223     l = 0 ! Index of level in the old grid
224     ! Loop over levels in the new grid
225     DO k = 1,Nglevels
226       Nw = 0 ! Number of weigths
227       wt = 0.0 ! Sum of weights
228       ! Loop over levels in the old grid and accumulate total for weighted average
229       DO
230         l = l + 1
231         w = 0.0 ! Initialise weight to 0
232         ! Distances between edges of both grids
233         dbb = oldgrid_bot(l) - newgrid_bot(k)
234         dtb = oldgrid_top(l) - newgrid_bot(k)
235         dbt = oldgrid_bot(l) - newgrid_top(k)
236         dtt = oldgrid_top(l) - newgrid_top(k)
237         if (dbt >= 0.0) exit ! Do next level in the new grid
238         if (dtb > 0.0) then
239           if (dbb <= 0.0) then
240             if (dtt <= 0) then
241               w = dtb
242             else
243               w = newgrid_top(k) - newgrid_bot(k)
244             endif
245           else
246             if (dtt <= 0) then
247               w = oldgrid_top(l) - oldgrid_bot(l)
248             else
249               w = -dbt
250             endif
251           endif
252           ! If layers overlap (w/=0), then accumulate
253           if (w /= 0.0) then
254             Nw = Nw + 1
255             wt = wt + w
256             DO j=1,Ncolumns
257               if (lunits) then
258                 if (y(i,j,l) /= R_UNDEF) then
259                   yp = 10.0**(y(i,j,l)/10.0)
260                 else
261                   yp = 0.0
262                 endif
263               else
264                 yp = y(i,j,l)
265               endif
266               r(i,j,k) = r(i,j,k) + w*yp
267             enddo
268           endif
269         endif
270       enddo
271       l = l - 2
272       if (l < 1) l = 0
273       ! Calculate average in new grid
274       if (Nw > 0) then
275         DO j=1,Ncolumns
276           r(i,j,k) = r(i,j,k)/wt
277         enddo
278       endif
279     enddo
280   enddo
281
282   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
283   DO k=1,Nglevels
284     DO j=1,Ncolumns
285       DO i=1,Npoints
286         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
287           if (lunits) then
288             if (r(i,j,k) <= 0.0) then
289               r(i,j,k) = R_UNDEF
290             else
291               r(i,j,k) = 10.0*log10(r(i,j,k))
292             endif
293           endif
294         else ! Level below surface
295           r(i,j,k) = R_GROUND
296         endif
297       enddo
298     enddo
299   enddo
300
301END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
302
303END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.