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

Last change on this file since 3918 was 3233, checked in by idelkadi, 7 years ago

Nettoyage du code :

  • changement de nomes des routines pour avoir les memes nome des modules
  • corrections
  • dos2unix appliquee aux fichiers
  • 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 $
[1414]5!
6! Redistribution and use in source and binary forms, with or without modification, are permitted
[1262]7! provided that the following conditions are met:
[1414]8!
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.
[1414]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
[1262]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
[1414]33! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
[2428]34! Jan 2013 - G. Cesana        - Added betaperp and temperature arguments
35!                             - Added phase 3D/3Dtemperature/Map output variables in diag_lidar
[1262]36!
[1414]37!
[2428]38#include "cosp_defs.h"
[1262]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)
[1414]52
[1262]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
[1414]62   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
63
64   ! Local variables
[1262]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
[2428]73   real,dimension(:,:,:),allocatable ::  betaperptot_out, temp_in, temp_out
74   real,dimension(:,:),allocatable :: temp_c
[1414]75
[1262]76   Npoints  = gbx%Npoints
77   Nlevels  = gbx%Nlevels
78   Nhydro   = gbx%Nhydro
79   Ncolumns = gbx%Ncolumns
80   Nlr      = vgrid%Nlvgrid
[1396]81
[2428]82   if (cfg%LcfadLidarsr532) ok_lidar_cfad=.true.
[1262]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
[2428]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
[1262]102        !++++++++++++ Radar CFAD ++++++++++++++++
103        if (cfg%Lradar_sim) then
[1414]104            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
[1262]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
[1414]111            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
[1262]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)
[2428]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,:)
[2955]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
[2428]126
[1262]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 &
[2428]133                            ,temp_c,betatot_out,betaperptot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
[1262]134                            ,LIDAR_UNDEF,ok_lidar_cfad &
135                            ,stlidar%cfad_sr,stlidar%srbval &
[2955]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
[1262]140        endif
[2428]141
[1262]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, &
[2428]144                                    temp_c,betatot_out,betaperptot_out,betamol_c,Ze_out, &
[1414]145                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
[2955]146        deallocate(temp_in,temp_out,temp_c,betaperptot_out) !TIBO +temp_in
[2428]147
[1262]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 &
[2428]157                        ,sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
[1262]158                        ,LIDAR_UNDEF,ok_lidar_cfad &
159                        ,stlidar%cfad_sr,stlidar%srbval &
[2955]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
[1262]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, &
[2428]166                                    sglidar%temp_tot,sglidar%beta_tot,sglidar%betaperp_tot,sglidar%beta_mol,sgradar%Ze_tot, &
[1414]167                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
[1262]168   endif
169   ! Replace undef
[1414]170   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
[2955]171   where (stlidar%profSR   == LIDAR_UNDEF) stlidar%profSR   = R_UNDEF !TIBO
[1414]172   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
[2955]173   where (stlidar%lidarcldtype  == LIDAR_UNDEF) stlidar%lidarcldtype  = R_UNDEF !OPAQ
[1414]174   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
[2955]175   where (stlidar%cldtype  == LIDAR_UNDEF) stlidar%cldtype  = R_UNDEF           !OPAQ
[1414]176   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
[2428]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
[1262]180
181END SUBROUTINE COSP_STATS
182
183!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
185!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]186SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units)
[1262]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
[2428]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]
[1262]198   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
199   ! Output
[2428]200   real,dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid
[1262]201
202   ! Local variables
[1414]203   integer :: i,j,k
[1262]204   logical :: lunits
[1414]205   integer :: l
[2428]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.
[1414]213
[1262]214   lunits=.false.
215   if (present(log_units)) lunits=log_units
[1396]216
[2428]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
[1262]270         endif
271       enddo
[2428]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
[1262]280     enddo
[2428]281   enddo
282
283   ! Set points under surface to R_UNDEF, and change to dBZ if necessary
284   do k=1,Nglevels
[1414]285     do j=1,Ncolumns
286       do i=1,Npoints
[2428]287         if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level
288           if (lunits) then
[1414]289             if (r(i,j,k) <= 0.0) then
[2428]290               r(i,j,k) = R_UNDEF
[1414]291             else
[2428]292               r(i,j,k) = 10.0*log10(r(i,j,k))
[1414]293             endif
294           endif
[2428]295         else ! Level below surface
296           r(i,j,k) = R_GROUND
297         endif
[1262]298       enddo
299     enddo
[2428]300   enddo
[1262]301
[1414]302END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
303
[1262]304END MODULE MOD_COSP_STATS
Note: See TracBrowser for help on using the repository browser.