source: trunk/LMDZ.PLUTO/libf/phypluto/calc_cpp_mugaz.F90 @ 3947

Last change on this file since 3947 was 3691, checked in by afalco, 9 months ago

Pluto: added HCN;
some comments about what variables actually do in optcv.
AF
File size: 7.3 KB
Line 
1      subroutine calc_cpp_mugaz
2
3!==================================================================
4!     Purpose
5!     -------
6!     Computes atmospheric specific heat capacity and
7!     mean molar mass for the gas mixture defined in gases.def
8!     with standard values.
9!     In 1d, one can force the model to adopt these standard values
10!     by setting cpp_mugaz_mode=2 in callfis.def.
11!
12!     Authors
13!     -------
14!     Robin Wordsworth (2009)
15!     A. Spiga: make the routine OK with latest changes in rcm1d
16!     Jeremy Leconte (2022)
17!
18!==================================================================
19
20      use gases_h
21      use comcstfi_mod, only: cpp, mugaz
22      use callkeys_mod, only: cpp_mugaz_mode
23      use mod_phys_lmdz_para, only : is_master
24      implicit none
25
26      character(len=20),parameter :: myname="calc_cpp_mugaz"
27      real cpp_c
28      real mugaz_c
29
30      integer igas
31
32      cpp_c   = 0.0
33      mugaz_c = 0.0
34
35
36! compute mugaz
37      do igas=1,ngasmx
38
39         if(igas.eq.vgas)then
40            ! ignore variable gas in cpp calculation
41         else
42            !  Fill in specific heat cp (J/mol/K) for each gas
43            if(igas.eq.igas_CO2)then
44               mugaz_c = mugaz_c + 44.01*gfrac(igas)
45            elseif(igas.eq.igas_N2)then
46               mugaz_c = mugaz_c + 28.01*gfrac(igas)
47            elseif(igas.eq.igas_H2)then
48               mugaz_c = mugaz_c + 2.01*gfrac(igas)
49            elseif(igas.eq.igas_He)then
50               mugaz_c = mugaz_c + 4.003*gfrac(igas)
51            elseif(igas.eq.igas_H2O)then
52               mugaz_c = mugaz_c + 18.02*gfrac(igas)
53            elseif(igas.eq.igas_SO2)then
54               mugaz_c = mugaz_c + 64.066*gfrac(igas)
55            elseif(igas.eq.igas_H2S)then
56               mugaz_c = mugaz_c + 34.08*gfrac(igas)
57            elseif(igas.eq.igas_CH4)then
58               mugaz_c = mugaz_c + 16.04*gfrac(igas)
59            elseif(igas.eq.igas_NH3)then
60               mugaz_c = mugaz_c + 17.03*gfrac(igas)
61            elseif(igas.eq.igas_C2H6)then
62               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
63               mugaz_c = mugaz_c + 30.07*gfrac(igas)
64            elseif(igas.eq.igas_C2H2)then
65               ! C2H2 https://encyclopedia.airliquide.com/fr/acetylene
66               mugaz_c = mugaz_c + 26.04*gfrac(igas)
67            elseif(igas.eq.igas_C2H4)then
68               mugaz_c = mugaz_c + 28.054*gfrac(igas)
69            elseif(igas.eq.igas_CO)then
70               mugaz_c = mugaz_c + 28.01*gfrac(igas)
71            elseif(igas.eq.igas_OCS)then
72               ! https://pubchem.ncbi.nlm.nih.gov/compound/carbonyl_sulfide
73               mugaz_c = mugaz_c + 60.0751*gfrac(igas)
74            elseif(igas.eq.igas_HCl)then
75               mugaz_c = mugaz_c + 36.46*gfrac(igas)
76            elseif(igas.eq.igas_HCN)then
77               mugaz_c = mugaz_c + 35.85*gfrac(igas)
78            elseif(igas.eq.igas_HF)then
79               mugaz_c = mugaz_c + 20.01*gfrac(igas)
80            else
81              if (is_master) then
82               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
83              endif
84              call abort_physic(myname,'Gas species not recognised!',1)
85            endif
86         endif
87
88      enddo
89
90!compute cpp
91      do igas=1,ngasmx
92
93         if(igas.eq.vgas)then
94            ! ignore variable gas in cpp calculation
95         else
96            ! values at minimum temperature from Engineering Toolbox
97            ! (better for Pluto)
98            if(igas.eq.igas_CO2)then
99            ! https://www.engineeringtoolbox.com/carbon-dioxide-d_974.html
100               cpp_c   = cpp_c   + 0.709*gfrac(igas)*44.01/mugaz_c
101            elseif(igas.eq.igas_N2)then
102               cpp_c   = cpp_c   + 1.039*gfrac(igas)*28.01/mugaz_c
103            elseif(igas.eq.igas_H2)then
104               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
105            elseif(igas.eq.igas_He)then
106               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
107            elseif(igas.eq.igas_H2O)then
108               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
109            elseif(igas.eq.igas_SO2)then
110               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
111            elseif(igas.eq.igas_H2S)then
112               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
113            elseif(igas.eq.igas_CH4)then
114               cpp_c   = cpp_c   + 2.087*gfrac(igas)*16.04/mugaz_c
115            elseif(igas.eq.igas_NH3)then
116               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
117               print*,'WARNING, cpp for NH3 may be for liquid'
118            elseif(igas.eq.igas_C2H6)then
119               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
120               cpp_c   = cpp_c   + 1.535*gfrac(igas)*30.07/mugaz_c
121            elseif(igas.eq.igas_C2H2)then
122               ! https://encyclopedia.airliquide.com/fr/acetylene
123               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
124            !!!!! MODIF GG JAN 2019  (check source values !!)
125            elseif(igas.eq.igas_C2H4)then
126               ! https://www.engineeringtoolbox.com/ethylene-ethene-C2H4-properties-d_2104.html
127               cpp_c   = cpp_c   + 1.295*gfrac(igas)*28.054/mugaz_c
128            !!!!! MODIF GG JAN 2019  (check source values !!)
129            elseif(igas.eq.igas_CO)then
130               cpp_c   = cpp_c   + 1.039*gfrac(igas)*28.01/mugaz_c
131            elseif(igas.eq.igas_OCS)then
132               cpp_c   = cpp_c   + 0.6909*gfrac(igas)*60.07/mugaz_c
133            elseif(igas.eq.igas_HCl)then
134               cpp_c   = cpp_c   + 1.7087*gfrac(igas)*36.46/mugaz_c
135            elseif(igas.eq.igas_HCN)then
136               cpp_c   = cpp_c   + 1.7087*gfrac(igas)*36.46/mugaz_c
137            elseif(igas.eq.igas_HF)then
138               cpp_c   = cpp_c   + 1.4553*gfrac(igas)*20.01/mugaz_c
139            else
140              if (is_master) then
141               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
142              endif
143              call abort_physic(myname,'Gas species not recognised!',1)
144            endif
145         endif
146
147      enddo
148
149
150
151
152
153      cpp_c = 1000.0*cpp_c
154
155      if (is_master) then
156        print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
157        print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
158      endif
159
160      if(((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
161           (abs(1.-mugaz/mugaz_c).gt.1.e-6)).and. is_master ) then
162         ! Ehouarn: tolerate a small mismatch between computed/stored values
163         print*,'--> Values do not match with the predefined one ! (',cpp,',',mugaz,')'
164         print*,'    Because cp varies with temperature and that some gases may not appear in gases.def,'
165         print*,'    a small discrepancy might be completely normal.'
166         print*,'    But you might want to check that!'
167         print*,'    If you want to use the values calculated here, adjust cpp / mugaz in the dynamics via newstart (3d)'
168         print*,'    or use cpp_mugaz_mode=2 (if you are in 1d).'
169      endif
170
171      if (cpp_mugaz_mode==2) then
172          if (is_master) print*,'*** cpp_mugaz_mode==2, so setting cpp & mugaz to computations in calc_cpp_mugaz.'
173          mugaz = mugaz_c
174          cpp = cpp_c
175      else
176          if (is_master) print*,'*** Leaving cpp & mugaz equal to predefined values'
177          if (is_master) print*,'(either from dynamics (cpp_mugaz_mode=0) or callfis (cpp_mugaz_mode=1)).'
178      endif
179
180
181    end subroutine calc_cpp_mugaz
Note: See TracBrowser for help on using the repository browser.