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

Last change on this file since 3566 was 2831, checked in by emillour, 2 years ago

Generic PCM:
Add the possibility to include Venus-like aerosols (triggered by option
aerovenus=.true. in callphys.def); baseline is to use 5 distinct scatterers
but each may be turned on/off (via aerovenus1, aerovenus2, aerovenus2p,
aerovenus3, aerovenusUV flags which may be specified in callphys.def).
GG
File size: 6.8 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)
[869]47            elseif(igas.eq.igas_H2)then
[538]48               mugaz_c = mugaz_c + 2.01*gfrac(igas)
[869]49            elseif(igas.eq.igas_He)then
[716]50               mugaz_c = mugaz_c + 4.003*gfrac(igas)
[869]51            elseif(igas.eq.igas_H2O)then
[538]52               mugaz_c = mugaz_c + 18.02*gfrac(igas)
[869]53            elseif(igas.eq.igas_SO2)then
[716]54               mugaz_c = mugaz_c + 64.066*gfrac(igas)
[869]55            elseif(igas.eq.igas_H2S)then
[716]56               mugaz_c = mugaz_c + 34.08*gfrac(igas)
[869]57            elseif(igas.eq.igas_CH4)then
[538]58               mugaz_c = mugaz_c + 16.04*gfrac(igas)
[869]59            elseif(igas.eq.igas_NH3)then
[538]60               mugaz_c = mugaz_c + 17.03*gfrac(igas)
[869]61            elseif(igas.eq.igas_C2H6)then
[868]62               ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
63               mugaz_c = mugaz_c + 30.07*gfrac(igas)
[869]64            elseif(igas.eq.igas_C2H2)then
[868]65               ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
66               mugaz_c = mugaz_c + 26.04*gfrac(igas)
[2831]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)
[253]77            else
[2831]78              if (is_master) then
[253]79               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
[2831]80              endif
81              call abort_physic(myname,'Gas species not recognised!',1)
[253]82            endif
83         endif
84
85      enddo
86
[1364]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_CO2)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
[2831]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
[1364]132            else
[2831]133              if (is_master) then
[1364]134               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
[2831]135              endif
136              call abort_physic(myname,'Gas species not recognised!',1)
[1364]137            endif
138         endif
139
140      enddo
141
142
143
144
145
[538]146      cpp_c = 1000.0*cpp_c
[253]147
[2831]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
[253]152
[2635]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).'
[538]162      endif
[253]163
[2635]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.'
[589]166          mugaz = mugaz_c
167          cpp = cpp_c
168      else
[2635]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)).'
[589]171      endif
172
173
[253]174    end subroutine calc_cpp_mugaz
Note: See TracBrowser for help on using the repository browser.