source: trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 @ 1834

Last change on this file since 1834 was 1397, checked in by milmd, 10 years ago

In LMDZ.GENERIC replacement of all phystd .h files by module files.
File size: 5.3 KB
Line 
1      subroutine calc_cpp_mugaz
2
3!==================================================================
4!     Purpose
5!     -------
6!     Check to see if the atmospheric specific heat capacity and
7!     mean molar mass for the gas mixture defined in gases.def
8!     corresponds to what we're using. If it doesn't, abort run
9!     unless option 'check_cpp_match' is set to false in
10!     callphys.def.
11!
12!     Authors
13!     -------
14!     Robin Wordsworth (2009)
15!     A. Spiga: make the routine OK with latest changes in rcm1d
16!
17!==================================================================
18
19      use gases_h
20      use comcstfi_mod, only: cpp, mugaz
21      use callkeys_mod, only: check_cpp_match,force_cpp
22      implicit none
23
24      real cpp_c   
25      real mugaz_c
26
27      integer igas
28
29      cpp_c   = 0.0
30      mugaz_c = 0.0
31
32
33! compute mugaz
34      do igas=1,ngasmx
35
36         if(igas.eq.vgas)then
37            ! ignore variable gas in cpp calculation
38         else
39            ! all values at 300 K from Engineering Toolbox
40            if(igas.eq.igas_CO2)then
41               mugaz_c = mugaz_c + 44.01*gfrac(igas)
42            elseif(igas.eq.igas_N2)then
43               mugaz_c = mugaz_c + 28.01*gfrac(igas)
44            elseif(igas.eq.igas_H2)then
45               mugaz_c = mugaz_c + 2.01*gfrac(igas)
46            elseif(igas.eq.igas_He)then
47               mugaz_c = mugaz_c + 4.003*gfrac(igas)
48            elseif(igas.eq.igas_H2O)then
49               mugaz_c = mugaz_c + 18.02*gfrac(igas)
50            elseif(igas.eq.igas_SO2)then
51               mugaz_c = mugaz_c + 64.066*gfrac(igas)
52            elseif(igas.eq.igas_H2S)then
53               mugaz_c = mugaz_c + 34.08*gfrac(igas)
54            elseif(igas.eq.igas_CH4)then
55               mugaz_c = mugaz_c + 16.04*gfrac(igas)
56            elseif(igas.eq.igas_NH3)then
57               mugaz_c = mugaz_c + 17.03*gfrac(igas)
58            elseif(igas.eq.igas_C2H6)then
59               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
60               mugaz_c = mugaz_c + 30.07*gfrac(igas)
61            elseif(igas.eq.igas_C2H2)then
62               ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
63               mugaz_c = mugaz_c + 26.04*gfrac(igas)
64            else
65               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
66               call abort
67            endif
68         endif
69
70      enddo
71
72!compute cpp
73      do igas=1,ngasmx
74
75         if(igas.eq.vgas)then
76            ! ignore variable gas in cpp calculation
77         else
78            ! all values at 300 K from Engineering Toolbox
79            if(igas.eq.igas_CO2)then
80               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
81               !Mars conditions)
82               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
83            elseif(igas.eq.igas_N2)then
84               cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
85            elseif(igas.eq.igas_H2)then
86               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
87            elseif(igas.eq.igas_He)then
88               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
89            elseif(igas.eq.igas_H2O)then
90               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
91            elseif(igas.eq.igas_SO2)then
92               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
93            elseif(igas.eq.igas_H2S)then
94               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
95            elseif(igas.eq.igas_CH4)then
96               cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
97            elseif(igas.eq.igas_NH3)then
98               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
99               print*,'WARNING, cpp for NH3 may be for liquid'
100            elseif(igas.eq.igas_C2H6)then
101               ! C2H6
102               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
103               cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
104            elseif(igas.eq.igas_C2H2)then
105               ! C2H2
106               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
107               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
108            else
109               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
110               call abort
111            endif
112         endif
113
114      enddo
115
116
117
118
119
120      cpp_c = 1000.0*cpp_c
121
122      print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
123      print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
124      print*,'Predefined Cp in physics is ',cpp,'J kg^-1 K^-1'
125      print*,'Predefined Mg in physics is ',mugaz,'amu'
126
127      if (check_cpp_match) then
128         print*,'REQUEST TO CHECK cpp_match :'
129         if((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
130              (abs(1.-mugaz/mugaz_c).gt.1.e-6)) then
131            ! Ehouarn: tolerate a small mismatch between computed/stored values
132            print*,'--> Values do not match!'
133            print*,'    Either adjust cpp / mugaz via newstart to calculated values,'
134            print*,'    or set check_cpp_match to .false. in callphys.def.'
135            stop
136         else
137            print*,'--> OK. Settings match composition.'
138         endif
139      endif
140
141      if (.not.force_cpp) then
142          print*,'*** Setting cpp & mugaz to computations in calc_cpp_mugaz.'
143          mugaz = mugaz_c
144          cpp = cpp_c
145      else
146          print*,'*** Setting cpp & mugaz to predefined values.'
147      endif
148
149
150      return
151    end subroutine calc_cpp_mugaz
Note: See TracBrowser for help on using the repository browser.