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

Last change on this file since 3784 was 3641, checked in by mturbet, 5 months ago

generic continuum database
File size: 7.0 KB
RevLine 
[253]1      subroutine calc_cpp_mugaz
2
3!==================================================================
4!     Purpose
5!     -------
[2635]6!     Computes atmospheric specific heat capacity and
[538]7!     mean molar mass for the gas mixture defined in gases.def
[2635]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.
[253]11!
12!     Authors
13!     -------
14!     Robin Wordsworth (2009)
[589]15!     A. Spiga: make the routine OK with latest changes in rcm1d
[2635]16!     Jeremy Leconte (2022)
[253]17!
18!==================================================================
19
[471]20      use gases_h
[1384]21      use comcstfi_mod, only: cpp, mugaz
[2635]22      use callkeys_mod, only: cpp_mugaz_mode
23      use mod_phys_lmdz_para, only : is_master
[253]24      implicit none
25
[2831]26      character(len=20),parameter :: myname="calc_cpp_mugaz"
[589]27      real cpp_c   
[538]28      real mugaz_c
29
[253]30      integer igas
31
[538]32      cpp_c   = 0.0
33      mugaz_c = 0.0
[253]34
[1364]35
36! compute mugaz
[253]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
[869]43            if(igas.eq.igas_CO2)then
[538]44               mugaz_c = mugaz_c + 44.01*gfrac(igas)
[869]45            elseif(igas.eq.igas_N2)then
[538]46               mugaz_c = mugaz_c + 28.01*gfrac(igas)
[3641]47            elseif(igas.eq.igas_O2)then
48               mugaz_c = mugaz_c + 31.999*gfrac(igas)
[869]49            elseif(igas.eq.igas_H2)then
[538]50               mugaz_c = mugaz_c + 2.01*gfrac(igas)
[869]51            elseif(igas.eq.igas_He)then
[716]52               mugaz_c = mugaz_c + 4.003*gfrac(igas)
[869]53            elseif(igas.eq.igas_H2O)then
[538]54               mugaz_c = mugaz_c + 18.02*gfrac(igas)
[869]55            elseif(igas.eq.igas_SO2)then
[716]56               mugaz_c = mugaz_c + 64.066*gfrac(igas)
[869]57            elseif(igas.eq.igas_H2S)then
[716]58               mugaz_c = mugaz_c + 34.08*gfrac(igas)
[869]59            elseif(igas.eq.igas_CH4)then
[538]60               mugaz_c = mugaz_c + 16.04*gfrac(igas)
[869]61            elseif(igas.eq.igas_NH3)then
[538]62               mugaz_c = mugaz_c + 17.03*gfrac(igas)
[869]63            elseif(igas.eq.igas_C2H6)then
[868]64               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
65               mugaz_c = mugaz_c + 30.07*gfrac(igas)
[869]66            elseif(igas.eq.igas_C2H2)then
[868]67               ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
68               mugaz_c = mugaz_c + 26.04*gfrac(igas)
[2831]69            ! GG MODIF JAN2019
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)
[253]79            else
[2831]80              if (is_master) then
[253]81               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
[2831]82              endif
83              call abort_physic(myname,'Gas species not recognised!',1)
[253]84            endif
85         endif
86
87      enddo
88
[1364]89!compute cpp
90      do igas=1,ngasmx
91
92         if(igas.eq.vgas)then
93            ! ignore variable gas in cpp calculation
94         else
95            ! all values at 300 K from Engineering Toolbox
96            if(igas.eq.igas_CO2)then
97               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
98               !Mars conditions)
99               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
100            elseif(igas.eq.igas_N2)then
101               cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
[3641]102            elseif(igas.eq.igas_O2)then
103               cpp_c   = cpp_c   + 0.918*gfrac(igas)*31.999/mugaz_c
[1364]104            elseif(igas.eq.igas_H2)then
105               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
106            elseif(igas.eq.igas_He)then
107               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
108            elseif(igas.eq.igas_H2O)then
109               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
110            elseif(igas.eq.igas_SO2)then
111               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
112            elseif(igas.eq.igas_H2S)then
113               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
114            elseif(igas.eq.igas_CH4)then
115               cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
116            elseif(igas.eq.igas_NH3)then
117               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
118               print*,'WARNING, cpp for NH3 may be for liquid'
119            elseif(igas.eq.igas_C2H6)then
120               ! C2H6
121               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
122               cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
123            elseif(igas.eq.igas_C2H2)then
124               ! C2H2
125               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
126               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
[2831]127            !!!!! MODIF GG JAN 2019  (check source values !!)
128            elseif(igas.eq.igas_CO)then
129               cpp_c   = cpp_c   + 1.0425*gfrac(igas)*28.01/mugaz_c
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
[1364]136            else
[2831]137              if (is_master) then
[1364]138               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
[2831]139              endif
140              call abort_physic(myname,'Gas species not recognised!',1)
[1364]141            endif
142         endif
143
144      enddo
145
146
147
148
149
[538]150      cpp_c = 1000.0*cpp_c
[253]151
[2831]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
[253]156
[2635]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).'
[538]166      endif
[253]167
[2635]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.'
[589]170          mugaz = mugaz_c
171          cpp = cpp_c
172      else
[2635]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)).'
[589]175      endif
176
177
[253]178    end subroutine calc_cpp_mugaz
Note: See TracBrowser for help on using the repository browser.