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

Last change on this file since 3193 was 3184, checked in by afalco, 10 months ago

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