[305] | 1 | subroutine gradients_kcm(profil_flag,rho_v,rho_n,T,dTdp,dPvdp,dPndp) |
---|
| 2 | |
---|
| 3 | use params_h |
---|
[471] | 4 | use gases_h |
---|
[305] | 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 |
---|
[716] | 16 | !double precision cp_n,cp_v |
---|
| 17 | double precision cp_v |
---|
[305] | 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 | |
---|
[716] | 29 | !cp_n = cp_neutral(T) |
---|
[305] | 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 | |
---|
[366] | 43 | if(ngasmx.eq.1)then |
---|
| 44 | print*,'Cannot have moist adiabat with one gas...' |
---|
| 45 | stop |
---|
| 46 | endif |
---|
[305] | 47 | |
---|
[869] | 48 | if(trim(gnom(ngasmx)).eq.'H2O')then |
---|
[366] | 49 | |
---|
[305] | 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 | |
---|
[869] | 74 | elseif(trim(gnom(ngasmx)).eq.'NH3')then |
---|
[305] | 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 | |
---|
[366] | 110 | cp_v=0.0 |
---|
[869] | 111 | if(trim(gnom(ngasmx)).eq.'H2O')then |
---|
[305] | 112 | cp_v = (32.24+1.923d-3*T+1.055d-5*T**2-3.511d-9*T**3)/m_v |
---|
[869] | 113 | elseif(trim(gnom(ngasmx)).eq.'NH3')then |
---|
[305] | 114 | cp_v = 2.058d3 |
---|
[869] | 115 | elseif(trim(gnom(ngasmx)).eq.'CH4')then |
---|
[305] | 116 | cp_v = 2.226d3 |
---|
| 117 | endif |
---|
[1714] | 118 | |
---|
[305] | 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 | |
---|
| 126 | end subroutine gradients_kcm |
---|