Changeset 2998

Timestamp:

Jul 18, 2023, 2:23:06 PM (16 months ago)

Author:

llange

Message:

MARS PEM

• Fix some bug in water conservations (transfer from watercap to waterreservoir was wrong)
• Add the H2O glaciers in the main program.
• Add the distinction between co2 frost (qsurf) and perenial co2 ice. Will be developed in the Mars PCM in a future commit
• Thickness thresolhd to pass from frost to glaciers is now not harcoded but a constant fixed in the mars_pem_constants

LL

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• constants_marspem_mod.F90 (modified) (1 diff)
• pem.F90 (modified) (13 diffs)

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

@@ -37 +37 @@
       REAL,PARAMETER ::  Lco2 =  5.71E5                ! Pilorget and Forget 2016
 
+! Conversion H2O/CO2 frost to perenial frost and vice versa
+      REAL,PARAMETER :: threshold_water_frost2perenial = 1000. !~ 1m
+      REAL,PARAMETER :: threshold_co2_frost2perenial = 3200.   !~ 2m
+
 
       

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

@@ -45 +45 @@
                            albedodat, zmea, zstd, zsig, zgam, zthe, &
                            hmons, summit, base,albedo_h2o_frost, &
-                           frost_albedo_threshold,emissiv,watercaptag
+                           frost_albedo_threshold,emissiv,watercaptag,perenial_co2ice
       use dimradmars_mod, only: totcloudfrac, albedo
       use dust_param_mod, only: tauscaling
@@ -93 +93 @@
       use conf_pem_mod,     only: conf_pem
       use pemredem,         only:  pemdem0,pemdem1
-      use glaciers_mod,        only: co2glaciers_evol,h2oglaciersflow,co2glaciersflow,h2oglaciersflow
+      use glaciers_mod,        only: co2glaciers_evol,h2oglaciers_evol,co2glaciersflow,h2oglaciersflow
       use criterion_pem_stop_mod, only: criterion_waterice_stop,criterion_co2_stop
       use constants_marspem_mod,  only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, &
-                                      m_co2,m_noco2
+                                      m_co2,m_noco2,threshold_water_frost2perenial,threshold_co2_frost2perenial
       use evol_co2_ice_s_mod, only: evol_co2_ice_s
       use evol_h2o_ice_s_mod, only: evol_h2o_ice_s
@@ -209 +209 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
-      REAL                               :: watercap_saved          ! Value saved at the previous time step
       REAL, dimension(:,:),  allocatable :: min_co2_ice_1           ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
       REAL, dimension(:,:),  allocatable :: min_co2_ice_2           ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
@@ -221 +220 @@
       REAL, dimension(:,:),  allocatable :: tendencies_co2_ice_ini  ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
       REAL, dimension(:,:),  allocatable :: tendencies_h2o_ice      ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
-      REAL, dimension(:,:),  allocatable :: flag_co2flow(:,:)       !(ngrid,nslope): Flag where there is a glacier flow
-      REAL, dimension(:),    allocatable :: flag_co2flow_mesh(:)    !(ngrid)       : Flag where there is a glacier flow
+      REAL, dimension(:,:),  allocatable :: flag_co2flow(:,:)       !(ngrid,nslope): Flag where there is a CO2 glacier flow
+      REAL, dimension(:),    allocatable :: flag_co2flow_mesh(:)    !(ngrid)       : Flag where there is a CO2 glacier flow
+      REAL, dimension(:,:),  allocatable :: flag_h2oflow(:,:)       !(ngrid,nslope): Flag where there is a H2O glacier flow
+      REAL, dimension(:),    allocatable :: flag_h2oflow_mesh(:)    !(ngrid)       : Flag where there is a H2O glacier flow
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL 
@@ -255 +256 @@
     REAL :: timestep                     ! timestep [s]
     REAL :: watercap_sum                 ! total mass of water cap [kg/m^2] 
+    REAL :: water_sum                    ! total mass of water in the mesh [kg/m^2] 
 
 #ifdef CPP_STD
@@ -444 +446 @@
               tsurf,tsoil,albedo,emis,   &
               q2,qsurf,tauscaling,totcloudfrac,wstar,     &
-              watercap,def_slope,def_slope_mean,subslope_dist)
+              watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
 
 ! Remove unphysical values of surface tracer
@@ -545 +547 @@
      allocate(flag_co2flow(ngrid,nslope))
      allocate(flag_co2flow_mesh(ngrid))
-
-       flag_co2flow(:,:) = 0.     
-       flag_co2flow_mesh(:) = 0.
-
+     allocate(flag_h2oflow(ngrid,nslope))
+     allocate(flag_h2oflow_mesh(ngrid))
+
+     flag_co2flow(:,:) = 0     
+     flag_co2flow_mesh(:) = 0
+     flag_h2oflow(:,:) = 0     
+     flag_h2oflow_mesh(:) = 0
 !---------------------------- READ GCM data --------------------- 
 
@@ -909 +914 @@
 
 ! II.b. Evolution of the ice
-      print *, "Evolution of h2o ice"
-     
+     print *, "Evolution of h2o ice"
      call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
      
+     print *, "Evolution of co2 ice"
+     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
+
      DO islope=1, nslope
        write(str2(1:2),'(i2.2)') islope
        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
        call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
+       call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
+       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
      ENDDO
-
-      print *, "Evolution of co2 ice"
-     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
 
 !------------------------
@@ -927 +933 @@
 !    II_a update pressure, ice and tracers 
 !    II_b Evolution of the ice
-!    II_c CO2 glaciers flows 
+!    II_c CO2  & H2O glaciers flows 
 
 !------------------------
@@ -936 +942 @@
        call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
                          global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
+      endif
+
+      print *, "H2O glacier flows"
+
+      if (h2oglaciersflow) then
+        call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
       endif
 
@@ -1167 +1179 @@
          watercap_sum=0.
          DO islope=1,nslope
-           watercap_saved = watercap(ig,islope)
-           if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then
-             qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
+           if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
+             qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
            else
+!             2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
              watercap(ig,islope)=watercap(ig,islope)+qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
              qsurf(ig,igcm_h2o_ice,islope)=0.
            endif
-           watercap_sum=watercap_sum+(watercap(ig,islope)-watercap_saved)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+           watercap_sum=watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
            watercap(ig,islope)=0.
          enddo
@@ -1182 +1194 @@
 
      DO ig=1,ngrid
+       water_sum = 0.
        DO islope=1,nslope
-         if(qsurf(ig,igcm_h2o_ice,islope).GT.500) then
+         water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+       ENDDO
+       if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming. 
+         if(water_sum.gt.threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
             watercaptag(ig)=.true.
-            qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-250
-            water_reservoir(ig)=water_reservoir(ig)+250
+            water_reservoir(ig)=water_reservoir(ig)+threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
+            DO islope = 1,nslope
+              qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.)
+            ENDDO    
          endif
-       enddo
-     enddo
+       else ! let's check that the infinite source of water disapear
+         if((water_sum + water_reservoir(ig)).lt.threshold_water_frost2perenial) then
+           watercaptag(ig)=.false.
+           DO islope = 1,nslope
+              qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
+           ENDDO
+           water_reservoir(ig) = 0.
+         endif
+       endif
+     enddo
+
+! CO2 ice
 
      DO ig=1,ngrid
-       if(water_reservoir(ig).LT. 10) then
-         watercaptag(ig)=.false.
-         DO islope=1,nslope
-           qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
-         ENDDO
-       endif
-     enddo
+       DO islope=1,nslope
+         if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then
+           perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
+           qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
+         endif
+       ENDDO
+     ENDDO
+       
 
 ! III_a.2  Tsoil update (for startfi)
@@ -1309 +1338 @@
                      tsurf,tsoil,inertiesoil,albedo,             &
                      emis,q2,qsurf,tauscaling,totcloudfrac,wstar,&
-                     watercap)
+                     watercap,perenial_co2ice)
 #else
      call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &