source: LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_utils.F90 @ 5218

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