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

Last change on this file since 3686 was 3669, checked in by afalco, 13 months ago

Pluto: add C2H4 as gas.
AF
File size: 7.1 KB
RevLine 
[3184]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
[3353]13!     -------
[3184]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"
[3353]27      real cpp_c
[3184]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            ! all values at 300 K from Engineering Toolbox
[3353]43            if(igas.eq.igas_CO2)then
[3184]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)
[3353]61            elseif(igas.eq.igas_C2H6)then
[3184]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
[3669]65               ! C2H2 https://encyclopedia.airliquide.com/fr/acetylene
[3184]66               mugaz_c = mugaz_c + 26.04*gfrac(igas)
67            ! GG MODIF JAN2019
[3669]68            elseif(igas.eq.igas_C2H4)then
69               mugaz_c = mugaz_c + 28.054*gfrac(igas)
[3184]70            elseif(igas.eq.igas_CO)then
71               mugaz_c = mugaz_c + 28.01*gfrac(igas)
72            elseif(igas.eq.igas_OCS)then
73               ! https://pubchem.ncbi.nlm.nih.gov/compound/carbonyl_sulfide
74               mugaz_c = mugaz_c + 60.0751*gfrac(igas)
75            elseif(igas.eq.igas_HCl)then
76               mugaz_c = mugaz_c + 36.46*gfrac(igas)
77            elseif(igas.eq.igas_HF)then
78               mugaz_c = mugaz_c + 20.01*gfrac(igas)
79            else
80              if (is_master) then
81               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
82              endif
83              call abort_physic(myname,'Gas species not recognised!',1)
84            endif
85         endif
86
87      enddo
88
89!compute cpp
90      do igas=1,ngasmx
91
92         if(igas.eq.vgas)then
93            ! ignore variable gas in cpp calculation
94         else
[3669]95            ! values at minimum temperature from Engineering Toolbox
96            ! (better for Pluto)
[3353]97            if(igas.eq.igas_CO2)then
[3669]98            ! https://www.engineeringtoolbox.com/carbon-dioxide-d_974.html
99               cpp_c   = cpp_c   + 0.709*gfrac(igas)*44.01/mugaz_c
[3184]100            elseif(igas.eq.igas_N2)then
[3669]101               cpp_c   = cpp_c   + 1.039*gfrac(igas)*28.01/mugaz_c
[3184]102            elseif(igas.eq.igas_H2)then
103               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
104            elseif(igas.eq.igas_He)then
105               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
106            elseif(igas.eq.igas_H2O)then
107               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
108            elseif(igas.eq.igas_SO2)then
109               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
110            elseif(igas.eq.igas_H2S)then
111               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
112            elseif(igas.eq.igas_CH4)then
[3669]113               cpp_c   = cpp_c   + 2.087*gfrac(igas)*16.04/mugaz_c
[3184]114            elseif(igas.eq.igas_NH3)then
115               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
116               print*,'WARNING, cpp for NH3 may be for liquid'
117            elseif(igas.eq.igas_C2H6)then
118               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
[3669]119               cpp_c   = cpp_c   + 1.535*gfrac(igas)*30.07/mugaz_c
[3184]120            elseif(igas.eq.igas_C2H2)then
[3669]121               ! https://encyclopedia.airliquide.com/fr/acetylene
[3184]122               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
123            !!!!! MODIF GG JAN 2019  (check source values !!)
[3669]124            elseif(igas.eq.igas_C2H4)then
125               ! https://www.engineeringtoolbox.com/ethylene-ethene-C2H4-properties-d_2104.html
126               cpp_c   = cpp_c   + 1.53*gfrac(igas)*28.054/mugaz_c
127            !!!!! MODIF GG JAN 2019  (check source values !!)
[3184]128            elseif(igas.eq.igas_CO)then
[3669]129               cpp_c   = cpp_c   + 1.039*gfrac(igas)*28.01/mugaz_c
[3184]130            elseif(igas.eq.igas_OCS)then
131               cpp_c   = cpp_c   + 0.6909*gfrac(igas)*60.07/mugaz_c
132            elseif(igas.eq.igas_HCl)then
133               cpp_c   = cpp_c   + 1.7087*gfrac(igas)*36.46/mugaz_c
134            elseif(igas.eq.igas_HF)then
135               cpp_c   = cpp_c   + 1.4553*gfrac(igas)*20.01/mugaz_c
136            else
137              if (is_master) then
138               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
139              endif
140              call abort_physic(myname,'Gas species not recognised!',1)
141            endif
142         endif
143
144      enddo
145
146
147
148
149
150      cpp_c = 1000.0*cpp_c
151
152      if (is_master) then
153        print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
154        print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
155      endif
156
157      if(((abs(1.-cpp/cpp_c).gt.1.e-6) .or.  &
158           (abs(1.-mugaz/mugaz_c).gt.1.e-6)).and. is_master ) then
159         ! Ehouarn: tolerate a small mismatch between computed/stored values
160         print*,'--> Values do not match with the predefined one ! (',cpp,',',mugaz,')'
161         print*,'    Because cp varies with temperature and that some gases may not appear in gases.def,'
162         print*,'    a small discrepancy might be completely normal.'
163         print*,'    But you might want to check that!'
164         print*,'    If you want to use the values calculated here, adjust cpp / mugaz in the dynamics via newstart (3d)'
165         print*,'    or use cpp_mugaz_mode=2 (if you are in 1d).'
166      endif
167
168      if (cpp_mugaz_mode==2) then
169          if (is_master) print*,'*** cpp_mugaz_mode==2, so setting cpp & mugaz to computations in calc_cpp_mugaz.'
170          mugaz = mugaz_c
171          cpp = cpp_c
172      else
173          if (is_master) print*,'*** Leaving cpp & mugaz equal to predefined values'
174          if (is_master) print*,'(either from dynamics (cpp_mugaz_mode=0) or callfis (cpp_mugaz_mode=1)).'
175      endif
176
177
178    end subroutine calc_cpp_mugaz
Note: See TracBrowser for help on using the repository browser.