source: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90 @ 1798

Last change on this file since 1798 was 1415, checked in by idelkadi, 14 years ago

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

File size: 12.4 KB
RevLine 
[1262]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!
29
30MODULE MOD_COSP_UTILS
31  USE MOD_COSP_CONSTANTS
32  IMPLICIT NONE
33
34  INTERFACE Z_TO_DBZ
35    MODULE PROCEDURE Z_TO_DBZ_2D,Z_TO_DBZ_3D,Z_TO_DBZ_4D
36  END INTERFACE
37
38  INTERFACE COSP_CHECK_INPUT
39    MODULE PROCEDURE COSP_CHECK_INPUT_1D,COSP_CHECK_INPUT_2D,COSP_CHECK_INPUT_3D
40  END INTERFACE
41CONTAINS
42
[1415]43!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44!------------------- SUBROUTINE COSP_PRECIP_MXRATIO --------------
45!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46SUBROUTINE COSP_PRECIP_MXRATIO(Npoints,Nlevels,Ncolumns,p,T,prec_frac,prec_type, &
47                          n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2, &
48                          flux,mxratio)
[1262]49
[1415]50    ! Input arguments, (IN)
51    integer,intent(in) :: Npoints,Nlevels,Ncolumns
52    real,intent(in),dimension(Npoints,Nlevels) :: p,T,flux
53    real,intent(in),dimension(Npoints,Ncolumns,Nlevels) :: prec_frac
54    real,intent(in) :: n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2,prec_type
55    ! Input arguments, (OUT)
56    real,intent(out),dimension(Npoints,Ncolumns,Nlevels) :: mxratio
57    ! Local variables
58    integer :: i,j,k
59    real :: sigma,one_over_xip1,xi,rho0,rho
60   
61    mxratio = 0.0
62
63    if (n_ax >= 0.0) then ! N_ax is used to control which hydrometeors need to be computed
64        !gamma1  = gamma(alpha_x + b_x + d_x + 1.0)
65        !gamma2  = gamma(alpha_x + b_x + 1.0)
66        xi      = d_x/(alpha_x + b_x - n_bx + 1.0)
67        rho0    = 1.29
68        sigma   = (gamma2/(gamma1*c_x))*(n_ax*a_x*gamma2)**xi
69        one_over_xip1 = 1.0/(xi + 1.0)
70       
71        do k=1,Nlevels
72            do j=1,Ncolumns
73                do i=1,Npoints
74                    if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then
75                        rho = p(i,k)/(287.05*T(i,k))
76                        mxratio(i,j,k)=(flux(i,k)*((rho/rho0)**g_x)*sigma)**one_over_xip1
77                        mxratio(i,j,k)=mxratio(i,j,k)/rho
78                    endif
79                enddo
80            enddo
81        enddo
82    endif
83END SUBROUTINE COSP_PRECIP_MXRATIO
84
85
[1262]86!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87!------------------- SUBROUTINE ZERO_INT -------------------------
88!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89ELEMENTAL SUBROUTINE ZERO_INT(x,y01,y02,y03,y04,y05,y06,y07,y08,y09,y10, &
90                                 y11,y12,y13,y14,y15,y16,y17,y18,y19,y20, &
91                                 y21,y22,y23,y24,y25,y26,y27,y28,y29,y30)
92
93  integer,intent(inout) :: x
94  integer,intent(inout),optional :: y01,y02,y03,y04,y05,y06,y07,y08,y09,y10, &
95                                    y11,y12,y13,y14,y15,y16,y17,y18,y19,y20, &
96                                    y21,y22,y23,y24,y25,y26,y27,y28,y29,y30
97  x = 0
98  if (present(y01)) y01 = 0
99  if (present(y02)) y02 = 0
100  if (present(y03)) y03 = 0
101  if (present(y04)) y04 = 0
102  if (present(y05)) y05 = 0
103  if (present(y06)) y06 = 0
104  if (present(y07)) y07 = 0
105  if (present(y08)) y08 = 0
106  if (present(y09)) y09 = 0
107  if (present(y10)) y10 = 0
108  if (present(y11)) y11 = 0
109  if (present(y12)) y12 = 0
110  if (present(y13)) y13 = 0
111  if (present(y14)) y14 = 0
112  if (present(y15)) y15 = 0
113  if (present(y16)) y16 = 0
114  if (present(y17)) y17 = 0
115  if (present(y18)) y18 = 0
116  if (present(y19)) y19 = 0
117  if (present(y20)) y20 = 0
118  if (present(y21)) y21 = 0
119  if (present(y22)) y22 = 0
120  if (present(y23)) y23 = 0
121  if (present(y24)) y24 = 0
122  if (present(y25)) y25 = 0
123  if (present(y26)) y26 = 0
124  if (present(y27)) y27 = 0
125  if (present(y28)) y28 = 0
126  if (present(y29)) y29 = 0
127  if (present(y30)) y30 = 0
128END SUBROUTINE  ZERO_INT
129
130!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131!------------------- SUBROUTINE ZERO_REAL ------------------------
132!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133ELEMENTAL SUBROUTINE ZERO_REAL(x,y01,y02,y03,y04,y05,y06,y07,y08,y09,y10, &
134                                 y11,y12,y13,y14,y15,y16,y17,y18,y19,y20, &
135                                 y21,y22,y23,y24,y25,y26,y27,y28,y29,y30)
136
137  real,intent(inout) :: x
138  real,intent(inout),optional :: y01,y02,y03,y04,y05,y06,y07,y08,y09,y10, &
139                                 y11,y12,y13,y14,y15,y16,y17,y18,y19,y20, &
140                                 y21,y22,y23,y24,y25,y26,y27,y28,y29,y30
141  x = 0.0
142  if (present(y01)) y01 = 0.0
143  if (present(y02)) y02 = 0.0
144  if (present(y03)) y03 = 0.0
145  if (present(y04)) y04 = 0.0
146  if (present(y05)) y05 = 0.0
147  if (present(y06)) y06 = 0.0
148  if (present(y07)) y07 = 0.0
149  if (present(y08)) y08 = 0.0
150  if (present(y09)) y09 = 0.0
151  if (present(y10)) y10 = 0.0
152  if (present(y11)) y11 = 0.0
153  if (present(y12)) y12 = 0.0
154  if (present(y13)) y13 = 0.0
155  if (present(y14)) y14 = 0.0
156  if (present(y15)) y15 = 0.0
157  if (present(y16)) y16 = 0.0
158  if (present(y17)) y17 = 0.0
159  if (present(y18)) y18 = 0.0
160  if (present(y19)) y19 = 0.0
161  if (present(y20)) y20 = 0.0
162  if (present(y21)) y21 = 0.0
163  if (present(y22)) y22 = 0.0
164  if (present(y23)) y23 = 0.0
165  if (present(y24)) y24 = 0.0
166  if (present(y25)) y25 = 0.0
167  if (present(y26)) y26 = 0.0
168  if (present(y27)) y27 = 0.0
169  if (present(y28)) y28 = 0.0
170  if (present(y29)) y29 = 0.0
171  if (present(y30)) y30 = 0.0
172END SUBROUTINE  ZERO_REAL
173
174!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
175!--------------------- SUBROUTINE Z_TO_DBZ_2D --------------------
176!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
177  SUBROUTINE Z_TO_DBZ_2D(mdi,z)
178    real,intent(in) :: mdi
179    real,dimension(:,:),intent(inout) :: z
180    ! Reflectivity Z:
181    ! Input in [m3]
182    ! Output in dBZ, with Z in [mm6 m-3]
183   
184    ! 1.e18 to convert from [m3] to [mm6 m-3]
185    z = 1.e18*z
186    where (z > 1.0e-6) ! Limit to -60 dBZ
187      z = 10.0*log10(z)
188    elsewhere
189      z = mdi
190    end where 
191  END SUBROUTINE Z_TO_DBZ_2D
192!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
193!--------------------- SUBROUTINE Z_TO_DBZ_3D --------------------
194!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
195  SUBROUTINE Z_TO_DBZ_3D(mdi,z)
196    real,intent(in) :: mdi
197    real,dimension(:,:,:),intent(inout) :: z
198    ! Reflectivity Z:
199    ! Input in [m3]
200    ! Output in dBZ, with Z in [mm6 m-3]
201   
202    ! 1.e18 to convert from [m3] to [mm6 m-3]
203    z = 1.e18*z
204    where (z > 1.0e-6) ! Limit to -60 dBZ
205      z = 10.0*log10(z)
206    elsewhere
207      z = mdi
208    end where 
209  END SUBROUTINE Z_TO_DBZ_3D
210!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211!--------------------- SUBROUTINE Z_TO_DBZ_4D --------------------
212!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213  SUBROUTINE Z_TO_DBZ_4D(mdi,z)
214    real,intent(in) :: mdi
215    real,dimension(:,:,:,:),intent(inout) :: z
216    ! Reflectivity Z:
217    ! Input in [m3]
218    ! Output in dBZ, with Z in [mm6 m-3]
219   
220    ! 1.e18 to convert from [m3] to [mm6 m-3]
221    z = 1.e18*z
222    where (z > 1.0e-6) ! Limit to -60 dBZ
223      z = 10.0*log10(z)
224    elsewhere
225      z = mdi
226    end where 
227  END SUBROUTINE Z_TO_DBZ_4D
228
229!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
230!----------------- SUBROUTINES COSP_CHECK_INPUT_1D ---------------
231!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
232  SUBROUTINE COSP_CHECK_INPUT_1D(vname,x,min_val,max_val)
233    character(len=*) :: vname
234    real,intent(inout) :: x(:)
235    real,intent(in),optional :: min_val,max_val
236    logical :: l_min,l_max
237    character(len=128) :: pro_name='COSP_CHECK_INPUT_1D'
238   
239    l_min=.false.
240    l_max=.false.
241   
242    if (present(min_val)) then
243!       if (x < min_val) x = min_val
244      if (any(x < min_val)) then
245      l_min = .true.
246        where (x < min_val)
247          x = min_val
248        end where
249      endif
250    endif   
251    if (present(max_val)) then
252!       if (x > max_val) x = max_val
253      if (any(x > max_val)) then
254        l_max = .true.
255        where (x > max_val)
256          x = max_val
257        end where 
258      endif   
259    endif   
260   
261    if (l_min) print *,'----- WARNING: '//trim(pro_name)//': minimum value of '//trim(vname)//' set to: ',min_val
262    if (l_max) print *,'----- WARNING: '//trim(pro_name)//': maximum value of '//trim(vname)//' set to: ',max_val
263  END SUBROUTINE COSP_CHECK_INPUT_1D
264!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265!----------------- SUBROUTINES COSP_CHECK_INPUT_2D ---------------
266!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267  SUBROUTINE COSP_CHECK_INPUT_2D(vname,x,min_val,max_val)
268    character(len=*) :: vname
269    real,intent(inout) :: x(:,:)
270    real,intent(in),optional :: min_val,max_val
271    logical :: l_min,l_max
272    character(len=128) :: pro_name='COSP_CHECK_INPUT_2D'
273   
274    l_min=.false.
275    l_max=.false.
276   
277    if (present(min_val)) then
278!       if (x < min_val) x = min_val
279      if (any(x < min_val)) then
280      l_min = .true.
281        where (x < min_val)
282          x = min_val
283        end where
284      endif
285    endif   
286    if (present(max_val)) then
287!       if (x > max_val) x = max_val
288      if (any(x > max_val)) then
289        l_max = .true.
290        where (x > max_val)
291          x = max_val
292        end where 
293      endif   
294    endif   
295   
296    if (l_min) print *,'----- WARNING: '//trim(pro_name)//': minimum value of '//trim(vname)//' set to: ',min_val
297    if (l_max) print *,'----- WARNING: '//trim(pro_name)//': maximum value of '//trim(vname)//' set to: ',max_val
298  END SUBROUTINE COSP_CHECK_INPUT_2D
299!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
300!----------------- SUBROUTINES COSP_CHECK_INPUT_3D ---------------
301!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
302  SUBROUTINE COSP_CHECK_INPUT_3D(vname,x,min_val,max_val)
303    character(len=*) :: vname
304    real,intent(inout) :: x(:,:,:)
305    real,intent(in),optional :: min_val,max_val
306    logical :: l_min,l_max
307    character(len=128) :: pro_name='COSP_CHECK_INPUT_3D'
308   
309    l_min=.false.
310    l_max=.false.
311   
312    if (present(min_val)) then
313!       if (x < min_val) x = min_val
314      if (any(x < min_val)) then
315      l_min = .true.
316        where (x < min_val)
317          x = min_val
318        end where
319      endif
320    endif   
321    if (present(max_val)) then
322!       if (x > max_val) x = max_val
323      if (any(x > max_val)) then
324        l_max = .true.
325        where (x > max_val)
326          x = max_val
327        end where 
328      endif   
329    endif   
330   
331    if (l_min) print *,'----- WARNING: '//trim(pro_name)//': minimum value of '//trim(vname)//' set to: ',min_val
332    if (l_max) print *,'----- WARNING: '//trim(pro_name)//': maximum value of '//trim(vname)//' set to: ',max_val
333  END SUBROUTINE COSP_CHECK_INPUT_3D
334
335
336END MODULE MOD_COSP_UTILS
Note: See TracBrowser for help on using the repository browser.