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
RevLine 
[1262]1! (c) British Crown Copyright 2008, the Met Office.
2! All rights reserved.
[2428]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 $
[5099]5
[1414]6! Redistribution and use in source and binary forms, with or without modification, are permitted
[1262]7! provided that the following conditions are met:
[5099]8
[1414]9!     * Redistributions of source code must retain the above copyright notice, this list
[1262]10!       of conditions and the following disclaimer.
11!     * Redistributions in binary form must reproduce the above copyright notice, this list
[1414]12!       of conditions and the following disclaimer in the documentation and/or other materials
[1262]13!       provided with the distribution.
[1414]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
[1262]16!       permission.
[5099]17
[1414]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
[1262]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
[1414]32! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
[2428]33! Jan 2013 - G. Cesana        - Added betaperp and temperature arguments
34!                             - Added phase 3D/3Dtemperature/Map output variables in diag_lidar
[5099]35
36
[4619]37#include "cosp_defs.h"
[1262]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)
[1414]51
[1262]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
[1414]61   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
62
63   ! Local variables
[1262]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
[2428]72   real,dimension(:,:,:),allocatable ::  betaperptot_out, temp_in, temp_out
73   real,dimension(:,:),allocatable :: temp_c
[1414]74
[1262]75   Npoints  = gbx%Npoints
76   Nlevels  = gbx%Nlevels
77   Nhydro   = gbx%Nhydro
78   Ncolumns = gbx%Ncolumns
79   Nlr      = vgrid%Nlvgrid
[1396]80
[2428]81   if (cfg%LcfadLidarsr532) ok_lidar_cfad=.true.
[1262]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
[2428]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
[1262]101        !++++++++++++ Radar CFAD ++++++++++++++++
102        if (cfg%Lradar_sim) then
[1414]103            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
[1262]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
[1414]110            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
[1262]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)
[2428]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,:)
[2955]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
[2428]125
[1262]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 &
[2428]132                            ,temp_c,betatot_out,betaperptot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
[1262]133                            ,LIDAR_UNDEF,ok_lidar_cfad &
134                            ,stlidar%cfad_sr,stlidar%srbval &
[2955]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
[1262]139        endif
[2428]140
[1262]141        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
[5185]142        if (cfg%Lradar_sim.AND.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
[2428]143                                    temp_c,betatot_out,betaperptot_out,betamol_c,Ze_out, &
[1414]144                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
[2955]145        deallocate(temp_in,temp_out,temp_c,betaperptot_out) !TIBO +temp_in
[2428]146
[1262]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 &
[2428]156                        ,sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
[1262]157                        ,LIDAR_UNDEF,ok_lidar_cfad &
158                        ,stlidar%cfad_sr,stlidar%srbval &
[2955]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
[1262]163        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
[5185]164        if (cfg%Lradar_sim.AND.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
[2428]165                                    sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sgradar%Ze_tot, &
[1414]166                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
[1262]167   endif
168   ! Replace undef
[1414]169   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
[2955]170   where (stlidar%profSR   == LIDAR_UNDEF) stlidar%profSR   = R_UNDEF !TIBO
[1414]171   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
[2955]172   where (stlidar%lidarcldtype  == LIDAR_UNDEF) stlidar%lidarcldtype  = R_UNDEF !OPAQ
[1414]173   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
[2955]174   where (stlidar%cldtype  == LIDAR_UNDEF) stlidar%cldtype  = R_UNDEF           !OPAQ
[1414]175   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
[2428]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
[1262]179
180END SUBROUTINE COSP_STATS
181
182!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
183!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
184!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]185SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
[1262]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
[2428]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]
[1262]197   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
198   ! Output
[2428]199   real,dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
[1262]200
201   ! Local variables
[1414]202   integer :: i,j,k
[1262]203   logical :: lunits
[1414]204   integer :: l
[5158]205   REAL :: w ! Weight
206   REAL :: dbb, dtb, dbt, dtt ! Distances between edges of both grids
[2428]207   integer :: Nw  ! Number of weights
[5158]208   REAL :: wt  ! Sum of weights
[2428]209   real,dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid
[5158]210   REAL :: yp ! Local copy of y at a particular point.
[2428]211              ! This allows for change of units.
[1414]212
[1262]213   lunits=.false.
214   if (present(log_units)) lunits=log_units
[1396]215
[2428]216   r = 0.0
217
[5158]218   DO i=1,Npoints
[2428]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
[5158]225     DO k = 1,Nglevels
[2428]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
[5158]229       DO
[2428]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
[5158]256             DO j=1,Ncolumns
[2428]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
[1262]269         endif
270       enddo
[2428]271       l = l - 2
272       if (l < 1) l = 0
273       ! Calculate average in new grid
274       if (Nw > 0) then
[5158]275         DO j=1,Ncolumns
[2428]276           r(i,j,k) = r(i,j,k)/wt
277         enddo
278       endif
[1262]279     enddo
[2428]280   enddo
281
282   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
[5158]283   DO k=1,Nglevels
284     DO j=1,Ncolumns
285       DO i=1,Npoints
[2428]286         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
287           if (lunits) then
[1414]288             if (r(i,j,k) <= 0.0) then
[2428]289               r(i,j,k) = R_UNDEF
[1414]290             else
[2428]291               r(i,j,k) = 10.0*log10(r(i,j,k))
[1414]292             endif
293           endif
[2428]294         else ! Level below surface
295           r(i,j,k) = R_GROUND
296         endif
[1262]297       enddo
298     enddo
[2428]299   enddo
[1262]300
[1414]301END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
302
[1262]303END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.