source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/neutral.F @ 2303

Last change on this file since 2303 was 2175, checked in by jescribano, 10 years ago

SPLA code included for first time

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 #include "../phylmd/dimphy.h"
32c
33        real u10_mps(klon),ustar_mps(klon),obklen_m(klon)
34        real u10n_mps(klon)
35        real pi,von_karman
36c       parameter (pi = 3.141592653589793, von_karman = 0.4)   
37c pour etre coherent avec vk de bl_for_dms.F
38        parameter (pi = 3.141592653589793, von_karman = 0.35)   
39c
40        real phi, phi_inv, phi_inv_sq, f1, f2, f3, dum1, psi
41        integer i
42
43
44        psi = 0.
45        do i=1,klon
46
47        if (u10_mps(i) .lt. 0.) u10_mps(i) = 0.0
48       
49        if  (obklen_m(i) .lt. 0.) then
50                phi = (1. - 160./obklen_m(i))**(-0.25)
51                phi_inv = 1./phi
52                phi_inv_sq = 1./phi * 1./phi
53                f1 = (1. + phi_inv) / 2.
54                f2 = (1. + phi_inv_sq)/2.
55c following to avoid numerical overruns. recall tan(90deg)=infinity
56                dum1 = min (1.e24, phi_inv)
57                f3 = atan(dum1)
58                psi = 2.*log(f1) + log(f2) - 2.*f3 + pi/2.   
59        else if (obklen_m(i) .gt. 0.) then
60                psi = -50. / obklen_m(i)
61        end if
62
63        u10n_mps(i) = u10_mps(i) + (ustar_mps(i) * psi /von_karman )
64c u10n set to 0. if -1 < obklen < 20
65        if ((obklen_m(i).gt.-1.).and.(obklen_m(i).lt.20.)) then
66            u10n_mps(i) = 0.
67        endif
68        if (u10n_mps(i) .lt. 0.) u10n_mps(i) = 0.0
69
70        enddo
71        return
72        end
73c***********************************************************************
Note: See TracBrowser for help on using the repository browser.