Changeset 3019

Timestamp:

Jul 28, 2023, 7:08:13 PM (16 months ago)

Author:

llange

Message:

Mars PEM

• Following -r 2963, adapt the PEM so that the density of water wapor is correctly computed
• I've also modified the way Tsoil is initialized, but it might be changed in the future

LL

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• pem.F90 (modified) (5 diffs)
• pemetat0.F90 (modified) (4 diffs)

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -48 +48 @@
       use dimradmars_mod, only: totcloudfrac, albedo
       use dust_param_mod, only: tauscaling
-      use tracer_mod,     only: noms,igcm_h2o_ice,igcm_co2 ! tracer names
+      use tracer_mod,     only: noms,igcm_h2o_ice,igcm_co2,mmol,igcm_h2o_vap ! tracer names and molar masses
 #else 
 ! 1b: Modules specific from the Generic physiq
@@ -83 +83 @@
                            major_slope,ini_comslope_h
 
+#ifndef CPP_STD
+      USE comcstfi_h, only: r, mugaz
+#else
+      USE comcstfi_mod, only: r, mugaz
+#endif
 
       USE logic_mod,        ONLY: iflag_phys
       USE mod_const_mpi,    ONLY: COMM_LMDZ
       use time_phylmdz_mod, only: daysec,dtphys
-      USE comconst_mod,     ONLY: rad,g,r,cpp,pi
+      USE comconst_mod,     ONLY: rad,g,cpp,pi
       USE infotrac
       USE geometry_mod,     only: latitude_deg 
@@ -517 +522 @@
 #endif
 
-     DO nnq=1,nqtot
+     DO nnq=1,nqtot  ! Why not using ini_tracer ?
        if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
+       if(noms(nnq).eq."h2o_vap") then
+          igcm_h2o_vap = nnq
+          mmol(igcm_h2o_vap)=18.
+       endif
        if(noms(nnq).eq."co2") igcm_co2 = nnq
      ENDDO 
@@ -1058 +1067 @@
        do ig = 1,ngrid
          do isoil = 1,nsoilmx_PEM
-          watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)
+          watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r)
           if(isnan(Tsoil_locslope(ig,isoil))) then
             call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1) 
@@ -1089 +1098 @@
                              tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
                              h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
-     endif
+  
+          
+      totmassco2_adsorbded = 0.
+      totmassh2o_adsorbded = 0.
+      do ig = 1,ngrid
+       do islope =1, nslope
+        do l = 1,nsoilmx_PEM - 1
+           totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
+           cell_area(ig)
+           totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
+           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
+           cell_area(ig)
+        enddo
+       enddo
+      enddo
+      write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
+      write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded                
+      endif
   endif !soil_pem
 

trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

@@ -13 +13 @@
                                    TI_breccia,TI_bedrock
    use soil_thermalproperties_mod, only: update_soil_thermalproperties
+   use tracer_mod, only: mmol,igcm_h2o_vap ! tracer names and molar masses
+
+#ifndef CPP_STD
+      USE comcstfi_h, only: r, mugaz
+#else
+      USE comcstfi_mod, only: r, mugaz
+#endif
+
+
 #ifndef CPP_STD   
    use surfdat_h, only: watercaptag
@@ -248 +257 @@
       do isoil = nsoil_GCM+1,nsoil_PEM
         do ig = 1,ngrid
-        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)
+        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
         enddo
       enddo
@@ -446 +455 @@
     do it = 1,timelen
         do isoil = nsoil_GCM+1,nsoil_PEM
-        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
+      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
         enddo
      enddo
@@ -452 +461 @@
       do isoil = nsoil_GCM+1,nsoil_PEM
         do ig = 1,ngrid
-        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)
+        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
         enddo
       enddo