| 1 | subroutine 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 | |
|---|
| 126 | end subroutine gradients_kcm |
|---|
