source: LMDZ5/branches/testing/libf/phylmd/cosp/cosp_utils.F90 @ 2879

Last change on this file since 2879 was 2435, checked in by Laurent Fairhead, 9 years ago

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