source: LMDZ6/trunk/libf/phylmd/Dust/neutral.f90 @ 5296

Last change on this file since 5296 was 5271, checked in by abarral, 5 weeks ago

Move dimensions.h into a module
Nb: doesn't compile yet

File size: 2.8 KB
Line 
1!***********************************************************************
2  subroutine neutral(u10_mps,ustar_mps,obklen_m, &
3          u10n_mps )
4  !-----------------------------------------------------------------------     
5  ! subroutine to compute u10 neutral wind speed
6  ! inputs
7  ! u10_mps - wind speed at 10 m (m/s)
8  ! ustar_mps - friction velocity (m/s)
9  ! obklen_m - monin-obukhov length scale (m)
10  ! outputs
11  ! u10n_mps - wind speed at 10 m under neutral conditions (m/s)
12  ! following code assumes reference height Z is 10m, consistent with use
13  ! of u10 and u10_neutral.  If not, code
14  ! should be changed so that constants of 50. and 160. in equations
15  ! below are changed to -5 * Z and -16 * Z respectively.
16  ! Reference:  G. L. Geernaert.  'Bulk parameterizations for the
17  ! wind stress and heat fluxes,' in Surface Waves and Fluxes, Vol. I,
18  ! Current Theory, Geernaert and W.J. Plant, editors, Kluwer Academic
19  ! Publishers, Boston, MA, 1990.
20  ! subroutine written Feb 2001 by eg chapman
21  ! adapted to LMD-ZT by E. Cosme 310801
22  ! Following Will Shaw (PNL, Seattle) the theory applied for flux
23  ! calculation with the scheme of Nightingale et al. (2000) does not
24  ! hold anymore when -1<obklen<20. In this case, u10n is set to 0,
25  ! so that the transfer velocity  computed in nightingale.F will also
26  ! be 0. The flux is then set to 0.
27  !----------------------------------------------------------------------               
28  !
29  USE dimphy
30  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
31  !
32    real :: u10_mps(klon),ustar_mps(klon),obklen_m(klon)
33    real :: u10n_mps(klon)
34    real :: pi,von_karman
35    ! parameter (pi = 3.141592653589793, von_karman = 0.4)     
36  ! pour etre coherent avec vk de bl_for_dms.F
37    parameter (pi = 3.141592653589793, von_karman = 0.35)
38  !
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.
54  ! 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 )
63  ! 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
71end subroutine neutral
72!***********************************************************************
Note: See TracBrowser for help on using the repository browser.