source: LMDZ6/trunk/libf/phylmd/Dust/neutral.F @ 5134

Last change on this file since 5134 was 4593, checked in by yann meurdesoif, 18 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

File size: 2.6 KB
Line 
1c***********************************************************************
2        subroutine neutral(u10_mps,ustar_mps,obklen_m,
3     +         u10n_mps )
4c-----------------------------------------------------------------------       
5c subroutine to compute u10 neutral wind speed
6c inputs
7c       u10_mps - wind speed at 10 m (m/s)
8c       ustar_mps - friction velocity (m/s)
9c       obklen_m - monin-obukhov length scale (m)
10c outputs
11c       u10n_mps - wind speed at 10 m under neutral conditions (m/s)
12c following code assumes reference height Z is 10m, consistent with use
13c of u10 and u10_neutral.  If not, code
14c should be changed so that constants of 50. and 160. in equations
15c below are changed to -5 * Z and -16 * Z respectively.
16c Reference:  G. L. Geernaert.  'Bulk parameterizations for the
17c wind stress and heat fluxes,' in Surface Waves and Fluxes, Vol. I,
18c Current Theory, Geernaert and W.J. Plant, editors, Kluwer Academic
19c Publishers, Boston, MA, 1990.
20c subroutine written Feb 2001 by eg chapman
21c adapted to LMD-ZT by E. Cosme 310801
22c Following Will Shaw (PNL, Seattle) the theory applied for flux
23c calculation with the scheme of Nightingale et al. (2000) does not
24c hold anymore when -1<obklen<20. In this case, u10n is set to 0,
25c so that the transfer velocity  computed in nightingale.F will also
26c be 0. The flux is then set to 0.
27c----------------------------------------------------------------------         
28c
29      USE dimphy
30      INCLUDE "dimensions.h"
31c
32        real u10_mps(klon),ustar_mps(klon),obklen_m(klon)
33        real u10n_mps(klon)
34        real pi,von_karman
35c       parameter (pi = 3.141592653589793, von_karman = 0.4)   
36c pour etre coherent avec vk de bl_for_dms.F
37        parameter (pi = 3.141592653589793, von_karman = 0.35)   
38c
39        real phi, phi_inv, phi_inv_sq, f1, f2, f3, dum1, psi
40        integer i
41
42
43        psi = 0.
44        do i=1,klon
45
46        if (u10_mps(i) .lt. 0.) u10_mps(i) = 0.0
47       
48        if  (obklen_m(i) .lt. 0.) then
49                phi = (1. - 160./obklen_m(i))**(-0.25)
50                phi_inv = 1./phi
51                phi_inv_sq = 1./phi * 1./phi
52                f1 = (1. + phi_inv) / 2.
53                f2 = (1. + phi_inv_sq)/2.
54c following to avoid numerical overruns. recall tan(90deg)=infinity
55                dum1 = min (1.e24, phi_inv)
56                f3 = atan(dum1)
57                psi = 2.*log(f1) + log(f2) - 2.*f3 + pi/2.   
58        else if (obklen_m(i) .gt. 0.) then
59                psi = -50. / obklen_m(i)
60        end if
61
62        u10n_mps(i) = u10_mps(i) + (ustar_mps(i) * psi /von_karman )
63c u10n set to 0. if -1 < obklen < 20
64        if ((obklen_m(i).gt.-1.).and.(obklen_m(i).lt.20.)) then
65            u10n_mps(i) = 0.
66        endif
67        if (u10n_mps(i) .lt. 0.) u10n_mps(i) = 0.0
68
69        enddo
70        return
71        end
72c***********************************************************************
Note: See TracBrowser for help on using the repository browser.