source: LMDZ6/trunk/libf/phylmdiso/cosp/mod_cosp_utils.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

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