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

Last change on this file since 2613 was 1714, checked in by mturbet, 8 years ago

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