Changeset 3779 for trunk/LMDZ.COMMON/libf

Timestamp:

May 26, 2025, 4:21:44 PM (4 weeks ago)

Author:

jbclement

Message:

PEM:
Information about checking of CO2 mass balance.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (12 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt ¶

@@ -670 +670 @@
 - 'deposits' is renamed as 'layerings_map'
 - Few cleanings
+
+== 26/05/2025 == JBC
+Information about checking of CO2 mass balance.

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

@@ -88 +88 @@
     use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
     use surfini_mod,        only: surfini
-    use comconst_mod,       only: pi, rad, g, cpp, kappa
     use comcstfi_h,         only: mugaz
 #else
@@ -97 +96 @@
     use aerosol_mod,        only: iniaerosol
     use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
-    use comcstfi_mod,       only: pi, rad, g, mugaz, r
+    use comcstfi_mod,       only: pi, rad, g, r, cpp, rcp, mugaz
 #endif
 
 #ifndef CPP_1D
-    use comconst_mod,             only: r
+    use comconst_mod,             only: pi, rad, g, r, cpp, rcp => kappa
     use iniphysiq_mod,            only: iniphysiq
     use control_mod,              only: iphysiq, day_step, nsplit_phys
 #else
-    use comcstfi_h,               only: r
+    use comcstfi_h,               only: pi, rad, g, r, cpp, rcp
     use time_phylmdz_mod,         only: iphysiq, steps_per_sol
     use regular_lonlat_mod,       only: init_regular_lonlat
@@ -194 +193 @@
 real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
 real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
+real(kind = 16)                      :: totmass_co2ice, totmass_atmco2         ! Current total CO2 masses
+real(kind = 16)                      :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses
 
 ! Variables for the evolution of layered layerings_map
@@ -207 +208 @@
 
 ! Variables for surface and soil
-real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Average surface temperature [K]
-real, dimension(:,:),     allocatable :: tsurf_dev                        ! Grid points x Slope field: Surface temperature deviation [K]
-real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
-real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Average Soil Temperature [K]
-real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
-real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries             ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old         ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
-real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Average water surface density [kg/m^3]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
-real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
-real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
-real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
-real                                  :: totmassco2_adsorbed              ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
-real                                  :: totmassh2o_adsorbed              ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
-logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
-real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
-real, dimension(:,:,:),   allocatable :: ice_porefilling_old              ! ngrid x nslope: Ice pore filling at the previous iteration [m]
-real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
-real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
-real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
-real, dimension(:,:),     allocatable :: zshift_surf                      ! Elevation shift for the surface [m]
-real, dimension(:,:),     allocatable :: zlag                             ! Newly built lag thickness [m]
-real, dimension(:,:),     allocatable :: icetable_depth_old               ! Old depth of the ice table
+real, dimension(:,:),     allocatable :: tsurf_avg                          ! Grid points x Slope field: Average surface temperature [K]
+real, dimension(:,:),     allocatable :: tsurf_dev                          ! Grid points x Slope field: Surface temperature deviation [K]
+real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Grid points x Slope field: Average surface temperature of first call of the PCM [K]
+real, dimension(:,:,:),   allocatable :: tsoil_avg                          ! Grid points x Soil x Slope field: Average Soil Temperature [K]
+real, dimension(:,:),     allocatable :: tsoil_avg_old                      ! Grid points x Soil field: Average Soil Temperature at the previous time step [K]
+real, dimension(:,:,:),   allocatable :: tsoil_dev                          ! Grid points x Soil x Slope field: Soil temperature deviation [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                   ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries               ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old           ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
+real, dimension(:,:),     allocatable :: watersurf_density_avg              ! Grid points x Slope: Average water surface density [kg/m^3]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
+real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg          ! Grid points x Soil x Slopes: Average water soil density [kg/m^3]
+real, dimension(:),       allocatable :: delta_co2_adsorbed                 ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
+real, dimension(:),       allocatable :: delta_h2o_adsorbed                 ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
+real(kind = 16)                       :: totmass_adsco2, totmass_adsco2_ini ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
+real                                  :: totmass_adsh2o                     ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
+logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
+real, dimension(:,:),     allocatable :: icetable_thickness_old             ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
+real, dimension(:,:,:),   allocatable :: ice_porefilling_old                ! ngrid x nslope: Ice pore filling at the previous iteration [m]
+real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
+real, dimension(:),       allocatable :: porefill                           ! Pore filling (output) to compute the dynamic ice table
+real                                  :: ssi_depth                          ! Ice table depth (output) to compute the dynamic ice table
+real, dimension(:,:),     allocatable :: zshift_surf                        ! Elevation shift for the surface [m]
+real, dimension(:,:),     allocatable :: zlag                               ! Newly built lag thickness [m]
+real, dimension(:,:),     allocatable :: icetable_depth_old                 ! Old depth of the ice table
 
 ! Some variables for the PEM run
@@ -664 +665 @@
 is_co2ice_ini = .false.
 co2ice_disappeared = .false.
+totmass_co2ice_ini = 0.
+totmass_atmco2_ini = 0.
 if (layering_algo) then
     do ig = 1,ngrid
@@ -673 +676 @@
 endif
 do i = 1,ngrid
+    totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
     do islope = 1,nslope
+        totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)
         if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
         if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
@@ -693 +698 @@
 end where
 
+totmass_adsco2_ini = 0.
+totmass_adsh2o = 0.
 if (adsorption_pem) then
-    totmassco2_adsorbed = 0.
-    totmassh2o_adsorbed = 0.
     do ig = 1,ngrid
         do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
                 if (l == 1) then
-                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
+                   totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
+                   totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                 else
-                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
+                   totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                   totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
+                   totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                 endif
@@ -713 +718 @@
         enddo
     enddo
-    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed
-    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
+    totmass_adsco2 = totmass_adsco2_ini
+    write(*,*) "Tot mass of CO2 in the regolith =", totmass_adsco2
+    write(*,*) "Tot mass of H2O in the regolith =", totmass_adsh2o
 endif ! adsorption
 
@@ -967 +973 @@
 
 ! II_d.6 Update the mass of the regolith adsorbed
+        totmass_adsco2 = 0.
+        totmass_adsh2o = 0.
         if (adsorption_pem) then
             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
                                      tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,       &
                                      h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
-
-            totmassco2_adsorbed = 0.
-            totmassh2o_adsorbed = 0.
             do ig = 1,ngrid
                 do islope = 1,nslope
                     do l = 1,nsoilmx_PEM
                         if (l == 1) then
-                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
+                            totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
+                            totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         else
-                            totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
+                            totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                            totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
+                            totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         endif
@@ -991 +996 @@
                 enddo
             enddo
-            write(*,*) "Total mass of CO2 in the regolith =", totmassco2_adsorbed
-            write(*,*) "Total mass of H2O in the regolith =", totmassh2o_adsorbed
+            write(*,*) "Total mass of CO2 in the regolith =", totmass_adsco2
+            write(*,*) "Total mass of H2O in the regolith =", totmass_adsh2o
         endif
     endif !soil_pem
@@ -1029 +1034 @@
     enddo
     deallocate(flag_co2flow,flag_h2oflow)
+
+    ! Checking mass balance for CO2
+    totmass_co2ice = 0.
+    totmass_atmco2 = 0.
+    do ig = 1,ngrid
+        totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
+        do islope = 1,nslope
+            totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+        enddo 
+    enddo
+    write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini), ' %'
+    if ((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini) > 0.01) then
+        write(*,*) '  /!\ Warning: mass balance is not conseved!'
+        write(*,'(a,f8.3,a)') '       Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
+        write(*,'(a,f8.3,a)') '       CO2 ice mass balance         = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
+        write(*,'(a,f8.3,a)') '       Adsorbed CO2 mass balance    = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
+    endif
 
 !------------------------
@@ -1267 +1289 @@
     call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
     do i = 1,ip1jmp1
-        teta(i,:) = teta(i,:)*pdyn(i)**kappa
+        teta(i,:) = teta(i,:)*pdyn(i)**rcp
     enddo
     ! Correction on atmospheric pressure