Changeset 1364 for trunk/LMDZ.GENERIC

Timestamp:

Nov 6, 2014, 1:02:18 AM (10 years ago)

Author:

bclmd

Message:

Bug fix for calculating cpp

File:

• trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90 (modified) (4 diffs)

trunk/LMDZ.GENERIC/libf/phystd/calc_cpp_mugaz.F90

@@ -20 +20 @@
       implicit none
 
-!#include "dimensions.h"
-!#include "dimphys.h"
+#include "dimensions.h"
+#include "dimphys.h"
 #include "comcstfi.h"
 #include "callkeys.h"
@@ -33 +33 @@
       mugaz_c = 0.0
 
+
+! compute mugaz
       do igas=1,ngasmx
 
@@ -40 +42 @@
             ! all values at 300 K from Engineering Toolbox
             if(igas.eq.igas_CO2)then
-               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for Mars conditions) 
-               cpp_c   = cpp_c   + 0.846*gfrac(igas)
                mugaz_c = mugaz_c + 44.01*gfrac(igas)
             elseif(igas.eq.igas_N2)then
-               cpp_c   = cpp_c   + 1.040*gfrac(igas)
                mugaz_c = mugaz_c + 28.01*gfrac(igas)
             elseif(igas.eq.igas_H2)then
-               cpp_c   = cpp_c   + 14.31*gfrac(igas)
                mugaz_c = mugaz_c + 2.01*gfrac(igas)
             elseif(igas.eq.igas_He)then
-               cpp_c   = cpp_c   + 5.19*gfrac(igas)
                mugaz_c = mugaz_c + 4.003*gfrac(igas)
             elseif(igas.eq.igas_H2O)then
-               cpp_c   = cpp_c   + 1.864*gfrac(igas)
                mugaz_c = mugaz_c + 18.02*gfrac(igas)
             elseif(igas.eq.igas_SO2)then
-               cpp_c   = cpp_c   + 0.64*gfrac(igas)
                mugaz_c = mugaz_c + 64.066*gfrac(igas)
             elseif(igas.eq.igas_H2S)then
-               cpp_c   = cpp_c   + 1.003*gfrac(igas) ! from wikipedia...
                mugaz_c = mugaz_c + 34.08*gfrac(igas)
             elseif(igas.eq.igas_CH4)then
-               cpp_c   = cpp_c   + 2.226*gfrac(igas)
                mugaz_c = mugaz_c + 16.04*gfrac(igas)
             elseif(igas.eq.igas_NH3)then
-               cpp_c   = cpp_c   + 2.175*gfrac(igas)
                mugaz_c = mugaz_c + 17.03*gfrac(igas)
-               print*,'WARNING, cpp for NH3 may be for liquid'
             elseif(igas.eq.igas_C2H6)then 
                ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
-               cpp_c   = cpp_c   + 1.763*gfrac(igas)
                mugaz_c = mugaz_c + 30.07*gfrac(igas)
             elseif(igas.eq.igas_C2H2)then
                ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
-               cpp_c   = cpp_c   + 1.575*gfrac(igas)
                mugaz_c = mugaz_c + 26.04*gfrac(igas)
             else
@@ -83 +72 @@
 
       enddo
+
+!compute cpp
+      do igas=1,ngasmx
+
+         if(igas.eq.vgas)then
+            ! ignore variable gas in cpp calculation
+         else
+            ! all values at 300 K from Engineering Toolbox
+            if(igas.eq.igas_CO2)then
+               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for
+               !Mars conditions) 
+               cpp_c   = cpp_c   + 0.846*gfrac(igas)*44.01/mugaz_c
+            elseif(igas.eq.igas_N2)then
+               cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
+            elseif(igas.eq.igas_H2)then
+               cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
+            elseif(igas.eq.igas_He)then
+               cpp_c   = cpp_c   + 5.19*gfrac(igas)*4.003/mugaz_c
+            elseif(igas.eq.igas_H2O)then
+               cpp_c   = cpp_c   + 1.864*gfrac(igas)*18.02/mugaz_c
+            elseif(igas.eq.igas_SO2)then
+               cpp_c   = cpp_c   + 0.64*gfrac(igas)*64.066/mugaz_c
+            elseif(igas.eq.igas_H2S)then
+               cpp_c   = cpp_c   + 1.003*gfrac(igas)*34.08/mugaz_c ! from wikipedia...
+            elseif(igas.eq.igas_CH4)then
+               cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
+            elseif(igas.eq.igas_NH3)then
+               cpp_c   = cpp_c   + 2.175*gfrac(igas)*17.03/mugaz_c
+               print*,'WARNING, cpp for NH3 may be for liquid'
+            elseif(igas.eq.igas_C2H6)then
+               ! C2H6
+               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
+               cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
+            elseif(igas.eq.igas_C2H2)then
+               ! C2H2
+               ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
+               cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
+            else
+               print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
+               call abort
+            endif
+         endif
+
+      enddo
+
+
+
+
 
       cpp_c = 1000.0*cpp_c