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

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

Several new files added as part of the climate evolution model
(main program kcm.F90). Some general cleanup in physiq.F90 and
callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh
scattering corrected.
  • Property svn:executable set to *
File size: 3.0 KB
RevLine 
[305]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(gnom(2).eq.'H2O')then
45
46        call psat_H2O(T-2d-1,psat_minus)
47        call psat_H2O(T+2d-1,psat_plus)
48        call psat_H2O(T,press)
49
50        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
51        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
52
53        call therm(T-2d-1,rho_minus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
54                   nul,nul,press,s_minus,nul)
55        call therm(T+2d-1,rho_plus*1d-3,nul,nul,nul,nul,nul,nul,nul,&
56                   nul,nul,press,s_plus,nul)
57        s_c = 2.06 * log(T/273.15)
58
59        s_plus  = s_plus  * 1d3
60        s_minus = s_minus * 1d3
61        s_c     = s_c     * 1d3 ! convert to SI
62
63        if(T.lt.280.0)then
64           dpsat = press*1d6 * ( 1730.63*log(10.) / (T-39.714)**2 )
65        else
66           call tdpsdt(T,dpsat)
67           dpsat = dpsat * 1d6 / T
68        endif
69
70     elseif(gnom(2).eq.'NH3')then
71
72        call psat_NH3(T-2d-1,psat_minus)
73        call psat_NH3(T+2d-1,psat_plus)
74        call psat_NH3(T,press)
75
76        rho_minus = m_v*psat_minus*1d6/(Rc*(T-2d-1))
77        rho_plus  = m_v*psat_plus*1d6/(Rc*(T+2d-1))
78
79        call latheat_NH3(T-2d-1,nul,s_minus)
80        call latheat_NH3(T+2d-1,nul,s_plus)
81        call latheat_NH3(T,s_c,nul)
82
83        dpsat = press*1d6 * (-2*1.5609d-4*T + 0.1236)
84
85     endif
86
87     dsv  = (s_plus-s_minus)/4d-1  ! dsv*T = ds / d ln[T]
88     s_v  = (s_plus+s_minus)/2d0
89     dlnr = T/rho_v * (rho_plus-rho_minus)/4d-1  ! d ln[rho_v] / d ln[T]
90
91     if(rho_n/rho_v.lt.1e-5)then
92        dlna = -T*dsv/(s_v-s_c)
93     else
94        a_v  = rho_v/rho_n
95        dlna = (rmn*dlnr - cp_n + rmn - a_v*T*dsv)/(a_v*(s_v-s_c)+rmn) 
96        ! d ln[alpha_v] / d ln[T]
97        ! note cp_n + rmn = cv_n, which is what's required
98     endif
99
100     dTdp  = 1d0 / (dpsat + rho_n*rmn*(1d0 + dlnr - dlna)) ! c.f. Marcq S2.2.2
101     dPvdp = dTdp * dpsat
102     dPndp = 1d0 - dPvdp ! from p = p_v + p_n
103
104  case(0) ! dry
105
106     if(gnom(2).eq.'H2O')then
107        cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v
108     elseif(gnom(2).eq.'NH3')then
109        cp_v = 2.058d3
110     elseif(gnom(2).eq.'CH4')then
111        cp_v = 2.226d3
112     endif
113
114     dTdp  = 1/(rho_n*cp_n+rho_v*cp_v)
115     dPndp = 1/(1d0+m_n/m_v*rho_v/rho_n)
116     dPvdp = 1/(1d0+m_v/m_n*rho_n/rho_v)
117     
118  end select
119
120
121end subroutine gradients_kcm
Note: See TracBrowser for help on using the repository browser.