source: trunk/LMDZ.GENERIC/libf/phystd/calcenergy_kcm.F90 @ 603

Last change on this file since 603 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: 1.7 KB
Line 
1subroutine calcenergy_kcm(Tsurf,T,Play,Plev,Qsurf,Q,muvar,Eatmtot)
2
3
4use params_h, only : Rc
5use watercommon_h, only : mH2O
6implicit none
7
8!     ----------------------------------------------------------------
9!     Purpose: Calculate total energy of the (steam) atmosphere
10!     Authour: R. Wordsworth (2011)
11!     ----------------------------------------------------------------
12
13#include "dimensions.h"
14#include "dimphys.h"
15#include "comcstfi.h"
16!#include "callkeys.h"
17
18
19
20  ! inputs
21  real Tsurf,T(1:nlayermx)
22  real Play(1:nlayermx),Plev(1:nlayermx+1)
23  real Qsurf,Q(1:nlayermx)
24  real muvar(ngridmx,nlayermx+1)
25
26  ! internal
27  integer il
28  real m_n,m_v,cp_n,cp_v,cp_a,Eatmlay
29  real s_c,rho_v,L
30  double precision p_v, s_v, nul
31  real VMR(1:nlayermx)
32
33  ! output
34  real Eatmtot
35
36  ! functions
37  double precision cp_neutral
38
39  m_n = mugaz/1000.
40  m_v = mH2O/1000.
41
42
43  do il=1,nlayermx
44     VMR(il)=Q(il)*(muvar(1,il+1)/mH2O)
45  end do
46
47
48  Eatmtot = 0.0
49  do il=1,nlayermx
50
51
52     cp_n = cp_neutral(dble(T(il)))
53     cp_v = (32.24+1.923d-3*T(il) + 1.055d-5*T(il)**2-3.511d-9*T(il)**3)/m_v
54     ! / by m_v makes it per kg not per mol
55
56     call psat_H2O(dble(T(il)),p_v)
57     p_v   = p_v*1d6
58     rho_v = real(p_v)*m_v/(real(Rc)*T(il))
59
60     call therm(dble(T(il)),dble(rho_v)*1d-3,nul,nul,nul,nul,nul,nul,nul,&
61                nul,nul,nul,s_v,nul)
62
63     s_c = 2.06 * log(T(il)/273.15)
64     s_v = s_v * 1d3
65     s_c = s_c * 1d3
66     L   = (real(s_v)-s_c)*T(il)
67
68     cp_a    = (1-VMR(il))*cp_n + VMR(il)*cp_v
69     !cp_a    = (1-Q(il))*cp_n + Q(il)*cp_v
70
71     Eatmlay = ( cp_a*T(il) + L*VMR(il) ) * (Plev(il)-Plev(il+1))/g
72
73     Eatmtot = Eatmtot + Eatmlay
74  enddo
75
76  print*,'WARNING: E-calc currently assumes H2O saturation!'
77
78end subroutine calcenergy_kcm
79
Note: See TracBrowser for help on using the repository browser.