source: trunk/LMDZ.GENERIC/libf/phystd/gradients_kcm.F90 @ 374

Last change on this file since 374 was 366, checked in by rwordsworth, 13 years ago

OSR output bugs fixed.
Improvements to kcm for pure H2 atmospheres.
  • Property svn:executable set to *
File size: 3.2 KB
Line 
1subroutine gradients_kcm(profil_flag,rho_v,rho_n,T,dTdp,dPvdp,dPndp)
2
3  use params_h
4  implicit none
5
6#include "gases.h"
7
8  ! inputs
9  integer profil_flag ! 0 = dry, 1 = moist, 2 = isothermal
10  double precision rho_v,rho_n,T
11
12  ! outputs
13  double precision dTdp,dPndp,dPvdp
14  double precision a_v
15
16  ! internal
17  double precision cp_n,cp_v
18
19  double precision press, rho_plus, rho_minus, dVdT, rho_c
20  double precision dlnr,dlna,dpsat,dsv
21  double precision s_minus,s_plus
22  double precision s_v,s_c,L
23  double precision psat_plus,psat_minus,Pn
24  double precision nul
25
26  ! functions
27  double precision cp_neutral
28
29  cp_n = cp_neutral(T)
30
31
32  select case(profil_flag)
33  case(2) ! isothermal
34
35     dTdp  = 0.
36     a_v   = rho_v/rho_n ! constant here
37     dPndp = 1/(1d0+m_n/m_v*a_v)
38     dPvdp = 1 - dPndp
39
40  case(1) ! moist
41
42     Pn = rho_n*T*rmn
43
44     if(ngasmx.eq.1)then
45        print*,'Cannot have moist adiabat with one gas...'
46        stop
47     endif
48
49     if(gnom(ngasmx).eq.'H2O')then
50
51        call psat_H2O(T-2d-1,psat_minus)
52        call psat_H2O(T+2d-1,psat_plus)
53        call psat_H2O(T,press)
54
55        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
56        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
57
58        call therm(T-2d-1,rho_minus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
59                   nul,nul,press,s_minus,nul)
60        call therm(T+2d-1,rho_plus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
61                   nul,nul,press,s_plus,nul)
62        s_c = 2.06 * log(T/273.15)
63
64        s_plus  = s_plus  * 1d3
65        s_minus = s_minus * 1d3
66        s_c     = s_c     * 1d3 ! convert to SI
67
68        if(T.lt.280.0)then
69           dpsat = press*1d6 * ( 1730.63*log(10.) / (T-39.714)**2 )
70        else
71           call tdpsdt(T,dpsat)
72           dpsat = dpsat * 1d6 / T
73        endif
74
75     elseif(gnom(ngasmx).eq.'NH3')then
76
77        call psat_NH3(T-2d-1,psat_minus)
78        call psat_NH3(T+2d-1,psat_plus)
79        call psat_NH3(T,press)
80
81        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
82        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
83
84        call latheat_NH3(T-2d-1,nul,s_minus)
85        call latheat_NH3(T+2d-1,nul,s_plus)
86        call latheat_NH3(T,s_c,nul)
87
88        dpsat = press*1d6 * (-2*1.5609d-4*T + 0.1236)
89
90     endif
91
92     dsv  = (s_plus-s_minus)/4d-1  ! dsv*T = ds / d ln[T]
93     s_v  = (s_plus+s_minus)/2d0
94     dlnr = T/rho_v * (rho_plus-rho_minus)/4d-1  ! d ln[rho_v] / d ln[T]
95
96     if(rho_n/rho_v.lt.1e-5)then
97        dlna = -T*dsv/(s_v-s_c)
98     else
99        a_v  = rho_v/rho_n
100        dlna = (rmn*dlnr - cp_n + rmn - a_v*T*dsv)/(a_v*(s_v-s_c)+rmn) 
101        ! d ln[alpha_v] / d ln[T]
102        ! note cp_n + rmn = cv_n, which is what's required
103     endif
104
105     dTdp  = 1d0 / (dpsat + rho_n*rmn*(1d0 + dlnr - dlna)) ! c.f. Marcq S2.2.2
106     dPvdp = dTdp * dpsat
107     dPndp = 1d0 - dPvdp ! from p = p_v + p_n
108
109  case(0) ! dry
110
111     cp_v=0.0
112     if(gnom(ngasmx).eq.'H2O')then
113        cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v
114     elseif(gnom(ngasmx).eq.'NH3')then
115        cp_v = 2.058d3
116     elseif(gnom(ngasmx).eq.'CH4')then
117        cp_v = 2.226d3
118     endif
119
120     dTdp  = 1/(rho_n*cp_n+rho_v*cp_v)
121     dPndp = 1/(1d0+m_n/m_v*rho_v/rho_n)
122     dPvdp = 1/(1d0+m_v/m_n*rho_n/rho_v)
123     
124  end select
125
126
127end subroutine gradients_kcm
Note: See TracBrowser for help on using the repository browser.